Modified Level Set Method

Summary: Tomography means reconstructing the internal structure of a physical object using the X-ray images of the object taken from different directions. Mathematically, the problem is to recover a non-negative function f(x) from a collection of the line integrals of f. Here we show how to solve a dynamic tomography problem using a modified level set method, introduced by Niemi, Lassas and Siltanen in 2015. This method enforces regularity in the reconstruction not only in spatial, but also in temporal direction and produces acceptable results with very sparse collections of projection directions. The X-ray attenuation coefficient is modelled by the level set function inside the domain (defined by the level set) and by zero outside.

Figures 1-4: 1) Dimensionally perfect can was found for building the phantom. 2) Candles melting in a water bath. 3) Graphite stick and aluminum stick crossed and embedded with the melted candle wax. 4) The candle wax during the solidification phase.

Author: Salla-Maaria Latva-Äijö

Software: Matlab

Data: The open data set 3D cross

Download Package: LevelSetMethodBlogShare

Literature: [1] Dynamic multi-source X-ray tomography using a spacetime level set method, Niemi E, Lassas M, Kallonen A, Harhanen L, Hämäläinen K and Siltanen S, 2015, Journal of Computational Physics 291, pp. 218-237. [2] Dynamic x-ray tomography with multiple sources, E.Niemi, M. Lassas, S. Siltanen, in: 8th International Symposium on Image and Signal Processing and Analysis (ISPA), Sept. 2013, 2013, pp. 618–621. Documentation of the experimental data can be found on Arxive.

Experimental measurements

A simple target which is changing in time was prepared setting aluminum and graphite sticks in different angles in relation to each other. The target was stabilized with candle wax. When the target rotates and we observe its cross sections along the z-axis, the cross sections of the sticks are moving frame by frame. The dynamic phantom was measured in the X-ray laboratory of the University of Helsinki, which is build by Alexander Meaney. The tube voltage was set to 50 kV and the object was rotated one degree after every taken projection image. A 3D-volume video was reconstructed with the FDK algorithm of Matlab’s ASTRA toolbox. It serves us as a ground truth reconstruction. 

In real life applications, the aim is to minimize the X-ray dose to the target. One way of doing this is to collect a sparse data set, where X-ray images are taken from fewer directions. The level set method can tolerate this kind of highly incomplete data, so we down sampled the data, picking only every sixth column from the original sinogram. Final data could have been produced measuring only 60 angles.

Figures 5-9: 5) The X-ray setup in the industrial mathematics laboratory in the University of Helsinki. Photographed by Markus Juvonen. 6) The final phantom on an imaging platform. 7) One example of a total of 360 projection images taken from the target. 8-9)  An example of one 360 degree sinogram (in color + black and white) from one slice of the produced data.

Principle of the spacetime level set method

Level set methods consider that the attenuation function f(u) is given by the level set function itself inside the level set and the attenuation is zero outside the level set. We model our target as two-dimensional by a non-negative X-ray attenuation function f(u) and add the time variable t as a third dimension u = u(x, y, t) to enforce regularity also in the temporal direction.

We find the level set function u as a minimizer of the functional

where A is an operator modeling 2D Radon transforms measured at several times, β = (β₁, β₂, β₃) is a multi-index with |β| = β₁ + β₂ + β₃, and α>0 is a regularization parameter. For the special case n = 1, F_n is equivalent to non-negativity constrained Tikhonov regularization.

In the case n = 2, we would like to minimize the functional F₂, but drop for simplicity the mixed derivatives from the functional. The resulting discretized functional to be minimized is

Because this functional is non-differentiable due to the singularity of f at zero, we smooth-out the singularity using the approximation

before applying the Barzilai-Borwein gradient algorithm for minimization. We use δ = 10^-2 in numerical computations. We approximate the second derivatives of u with respect to x; y and t using the central difference approximations

Where the spacing h_t adjusts the amount of regularization in temporal direction. We start the Barzilai-Borwein algorithm from the initial guess u₀=0, choose the first step size to be λ₀ = 0.0001 and iterate as

where the step size λ_l is given by

The gradient   is of the form

where   is the diagonal Jacobian matrix of    with diagonal elements ,  and

consists of terms of the form

.

Here we apply the negative boundary condition u = −1 on the boundary of Ω rather than the zero condition, since we would like to have the level set function u to be negative outside the level set.

Figures 10-13: Examples of the reconstructions made with the Level Set method. On the left with colors and on the right with gray scale. By changing the choise of the level, we can highlight different features in the reconstructions. In 10 and 12, you can see the candle wax also, but in 11 and in 13, only the shapes of the moving sticks are visible.