A parametric level-set method for partially discrete tomography

TL;DR

This paper presents a parametric level-set approach for partially discrete tomography, effectively reconstructing anomalies with known grey-value within a smooth background, even from limited data, outperforming existing methods.

Contribution

The paper introduces a novel parametric level-set method using radial basis functions for anomaly geometry representation in tomography, with a bi-level optimization framework.

Findings

Successfully reconstructs anomaly geometry from limited data

Outperforms Total Variation, DART, and P-DART methods

Effective in numerical phantom tests

Abstract

This paper introduces a parametric level-set method for tomographic reconstruction of partially discrete images. Such images consist of a continuously varying background and an anomaly with a constant (known) grey-value. We represent the geometry of the anomaly using a level-set function, which we represent using radial basis functions. We pose the reconstruction problem as a bi-level optimization problem in terms of the background and coefficients for the level-set function. To constrain the background reconstruction we impose smoothness through Tikhonov regularization. The bi-level optimization problem is solved in an alternating fashion; in each iteration we first reconstruct the background and consequently update the level-set function. We test our method on numerical phantoms and show that we can successfully reconstruct the geometry of the anomaly, even from limited data. On these phantoms, our method outperforms Total Variation reconstruction, DART and P-DART.

Equations

$$p_i={\int }_{\mathcal{D}}u(\mathbf{x})\delta (s_i-\mathbf{n}({\theta }_i)\cdot \mathbf{x})\mathrm{d}\mathbf{x},$$

$$u(\mathbf{x})=\sum _{j=1}^n{u}_{j}b(\mathbf{x}-{\mathbf{x}}_j),$$

$$\mathbf{p}=W\mathbf{u},$$

$$\underset{\mathbf{u}}{\min }\frac{1}{2}\Vert W\mathbf{u}-\mathbf{p}{\Vert }_2^2+\frac{\lambda }{2}\Vert R\mathbf{u}{\Vert }_2^2,$$

$$u(\mathbf{x})=\{\begin{array}{cc}{u}_0(\mathbf{x})& \text{if}\mathbf{x}\in \Omega ,\\ u_1& \text{otherwise}.\end{array}$$

$$\Omega =\{\mathbf{x}|\varphi (\mathbf{x})>0\}.$$

$$u(\mathbf{x})=(1-h(\varphi (\mathbf{x})))u_0(\mathbf{x})+h(\varphi (\mathbf{x}))u_1,$$

$$\frac{\partial \varphi }{\partial t}+v|\nabla \varphi |=0,$$

$$\varphi (\mathbf{x})=\sum _{j=1}^{n^{\prime }}{\alpha }_j\Psi ({\beta }_j\Vert \mathbf{x}-{\mathbf{\chi }}_j{\Vert }_2),$$

$$a_{ij}=\Psi ({\beta }_j\Vert {\mathbf{x}}_i-{\mathbf{\chi }}_j{\Vert }_2),$$

$$\mathbf{u}=(1-h(A(\mathbf{\chi },\mathbf{\beta })\mathbf{\alpha }))\odot {\mathbf{u}}_0+h(A(\mathbf{\chi },\mathbf{\beta })\mathbf{\alpha })u_1,$$

$$\underset{\mathbf{\alpha }}{\min }\left\{f(\mathbf{\alpha })=\Vert W[(u_1-{\mathbf{u}}_0)\odot h_{ϵ}(A\mathbf{\alpha })]-(\mathbf{p}-W{\mathbf{u}}_0){\Vert }_2^2\right\},$$

$$\begin{array}{rl}\nabla f(\mathbf{\alpha })& =A^T{D}_{\mathbf{\alpha }}^T{W}^T\mathbf{r}(\mathbf{\alpha }),\\ H_{GN}(f(\mathbf{\alpha }))& =A^T{D}_{\mathbf{\alpha }}^T{W}^{T}W{D}_{\mathbf{\alpha }}A.\end{array}$$

$$D_{\mathbf{\alpha }}=\mathrm{diag}((u_1-{\mathbf{u}}_0)\odot h_{ϵ}^{\prime }(A\mathbf{\alpha })),$$

$${\mathbf{\alpha }}^{(k+1)}={\mathbf{\alpha }}^{(k)}-{\mu }^{(k)}{\left(H_{GN}(f({\mathbf{\alpha }}^{(k)}))\right)}^{-1}\nabla f({\mathbf{\alpha }}^{(k)}),$$

$$\delta (\mathbf{x})=\{\begin{array}{ll}0& \mathbf{x}\le -ϵ\\ \frac{1}{4(1-\mu )ϵ}\left(1+\frac{\mathbf{x}+(1-\mu )ϵ}{\mu ϵ}+\frac{1}{\pi }\sin (\pi \frac{\mathbf{x}+(1-\mu )ϵ}{\mu ϵ})\right)& -ϵ<\mathbf{x}\le -\mu ϵ\\ \frac{1}{2(1-\mu )ϵ}& -\mu ϵ<\mathbf{x}\le \mu ϵ\\ \frac{1}{4(1-\mu )ϵ}\left(1-\frac{\mathbf{x}-(1-\mu )ϵ}{\mu ϵ}-\frac{1}{\pi }\sin (\pi \frac{\mathbf{x}-(1-\mu )ϵ}{\mu ϵ})\right)& \mu ϵ<\mathbf{x}\le ϵ\\ 0& \mathbf{x}\ge ϵ\end{array}$$

$$\varphi (\mathbf{x})=\varphi ({\mathbf{x}}_0)+(\mathbf{x}-{\mathbf{x}}_0{)}^T\nabla \varphi ({\mathbf{x}}_0)+\mathcal{O}(\Vert \mathbf{x}-{\mathbf{x}}_0{\Vert }^2).$$

$$ϵ=\kappa \left(\frac{\max (\varphi (\mathbf{x}))-\min (\varphi (\mathbf{x}))}{\Delta x}\right)=\kappa \left(\frac{\max (A\mathbf{\alpha })-\min (A\mathbf{\alpha })}{\Delta x}\right),$$

$$\underset{{\mathbf{u}}_0,\mathbf{\alpha }}{\min }\left\{f(\mathbf{\alpha },{\mathbf{u}}_0):=\frac{1}{2}\Vert W[(1-h(A\mathbf{\alpha }){\mathbf{u}}_0+h(A\mathbf{\alpha })u_1]-\mathbf{p}{\Vert }_2^2+\frac{\lambda }{2}\Vert L{\mathbf{u}}_0{\Vert }_2^2\right\},$$

$$\overline{f}(\mathbf{\alpha })=\underset{{\mathbf{u}}_0}{\min }f(\mathbf{\alpha },{\mathbf{u}}_0).$$

$$\nabla \overline{f}(\mathbf{\alpha })$$

$${\nabla }^2\overline{f}(\mathbf{\alpha })$$

$${\mathbf{u}}_0^{(k+1)}$$

$${\mathbf{\alpha }}^{(k+1)}$$

$$\Psi (r)=(1-r{)}_{+}^8(32r^3+25r^2+8r+1).$$

Full text

11footnotetext: Paper submitted to $20^{th}$ International Conference on Discrete Geometry for Computer Imagery

1 Introduction

The need to reconstruct (quantitative) images of an object from tomographic measurements appears in many applications. At the heart of many of these applications is a projection model based on the Radon transform. Characterizing the object under investigation by a function $u(\mathbf{x})$ with $\mathbf{x}\in\mathcal{D}=[0,1]^{2}$, tomographic measurements are modeled as

[TABLE]

where $s_{i}\in[0,1]$ denotes the shift, $\theta_{i}\in[0,2\pi)$ denotes the angle and $\mathbf{n}(\theta)=(\cos\theta,\sin\theta)$. The goal is to retrieve $u$ from a number, $m$, of such measurements for various shifts and directions.

If the shifts and angles are regularly sampled, the transform can be inverted directly by Filtered back-projection or Fourier reconstruction [9]. A common approach for dealing with non-regularly sampled or missing data, is to express $u$ in terms of a basis

[TABLE]

where $b$ are piece-wise polynomial basis functions and $\{\mathbf{x}_{j}\}_{j=1}^{n}$ is a regular (pixel) grid. This leads to a set of $m$ linear equations in $n$ unknowns

[TABLE]

with $w_{ij}=\int_{\mathcal{D}}\,b(\mathbf{x}-\mathbf{x}_{j})\delta(s_{i}-\mathbf{n}(\theta_{i})\cdot\mathbf{x})\,\mathrm{d}\mathbf{x}$. Due to noise in the data or errors in the projection model the system of equations is inconsistent, so a solution may not exist. Furthermore, there may be many solutions that fit the observations equally well because the system is underdetermined. A standard approach to mitigate these issues is to formulate a regularized least-squares problem

[TABLE]

where $R$ is the regularization operator. Such a formulation is popular mainly because very efficient algorithms exist for solving it. Depending on the choice of $R$, however, this formulation forces the solution to have certain properties which may not reflect the truth. For example, setting $R$ to be the discrete Laplace operator will produce a smooth reconstruction, whereas setting $R$ to be the identity matrix forces the individual coefficients $u_{i}$ to be small. In many applications such quadratic regularization terms do not reflect the characteristics of the object we are reconstructing. For example, if we expect $u$ to be piecewise constant, we could use a Total Variation regularization term $\|R\mathbf{u}\|_{1}$ where $R$ is a discrete gradient operator [14]. Recently, a lot of progress has been made in developing efficient algorithms for solving such non-smooth optimization problems [6]. If the object under investigation is known to consist of only two distinct materials, the regularization can be formulated in terms of a non-convex constraint $\mathbf{u}\in\{u_{0},u_{1}\}^{n}$. The latter leads to a combinatorial optimization problem, solutions to which can be approximated using heuristic algorithms [3].

In this paper, we consider tomographic reconstruction of partially discrete objects that consist of a region of constant density embedded in a continuously varying background. In this case, neither the quadratic, Total Variation nor non-convex constraints by themselves are suitable. We therefore propose the following parametrization

[TABLE]

The inverse problem now consists of finding $u_{0}(\mathbf{x})$, $u_{1}$ and the set $\Omega$. We can subsequently apply suitable regularization to $u_{0}$ separately. To formulate a tractable optimization algorithm, we represent the set $\Omega$ using a level-set function $\phi(\mathbf{x})$ such that

[TABLE]

In the following sections, we discuss how to formulate a variational problem to reconstruct $\Omega$ and $u_{0}$ based on a parametric level-set representation of $\Omega$ and assuming we know $u_{1}$.

The outline of the paper is as follows. In section 2 we discuss the parametric level-set method and propose some practical heuristics for choosing various paramaters that occur in the formulation. A joint background-anomaly reconstruction algorithm for partially discrete tomography is discussed in section 3. The results on few moderately complicated numerical phantoms are presented in Section 4. We provide some concluding remarks in Section 5.

2 Level-set methods

In terms of the level-set function, we can express $u$ as

[TABLE]

where $h$ is the Heaviside function and the latter term represents the anomaly.

Level-set methods have received much attention in geometric inverse problems, interface tracking, segmentation and shape optimization. The reason being their ability to handle topological changes. The classical level-set method, introduced by Sethian and Osher [12], solves the Hamiltonian-Jacobi equation, also known as level-set equation.

[TABLE]

where $\phi:\mathbb{R}^{2}\times\mathbb{R}^{+}\rightarrow\mathbb{R}$ denotes the level-set function as a time-dependent quantity for representing the shape and $v$ denotes the normal velocity. In the inverse-problems setting, the velocity $v$ is often derived from the gradient of the cost function with respect to the model parameter [5], [7]. There are various numerical issues associated with the numerical solution of level-set equation, e.g. reinitialization of the level-set. We refer the interested reader to a seminal paper in level-set method [11] and its application to computational tomography [10].

Instead of taking this classical level-set approach, we employ a parametric level-set approach, first introduced by Aghasi et al [1]. In this method, the level-set function is parametrized using radial basis functions:

[TABLE]

where $\Psi(.)$ is a radial basis function, $\{\alpha_{j}\}_{j=1}^{n^{\prime}}$ and $\{\chi_{j}\}_{j=1}^{n^{\prime}}$ are the amplitudes and nodes respectively, and the parameters $\{\beta_{j}\}_{j=1}^{n^{\prime}}$ control the widths. Introducing the kernel matrix $A(\boldsymbol{\chi},\boldsymbol{\beta})$ with elements

[TABLE]

we can now express $\mathbf{u}$ as

[TABLE]

where $h$ is applied element-wise to the vector $A(\boldsymbol{\chi},\boldsymbol{\beta})\boldsymbol{\alpha}$ and $\odot$ denotes the element-wise (Hadamard) product. By choosing the parameters $(\boldsymbol{\chi},\boldsymbol{\beta},\boldsymbol{\alpha})$ appropriately we can represent any (smooth) shape. To simplify matters and make the resulting optimization problem more tractable, we consider a fixed regular grid $\{\boldsymbol{\chi}_{j}\}_{j=1}^{n^{\prime}}$ and a fixed width $\beta_{j}\equiv\beta$. In the following we choose $\beta$ in accordance with the gridspacing $\Delta\chi$ as $\beta=1/(\eta\Delta\chi)$, where $\eta$ determines the influence of RBF on its neighbors.

2.0.1 Example

To show that the reconstruction of level-set with a finitely many radial basis functions, we consider the level-set shown in Figure 1 (a). With ${n^{\prime}}=196$ RBFs, it is possible to reconstruct a smooth shape discretized on a grid with $n=256\times 256$ pixels.

Finally, the discretized reconstruction problem for determining the shape is now formulated as

[TABLE]

where $h_{\epsilon}$ is a smooth approximation of the Heaviside function. The gradient and Gauss-Newton Hessian of ${f(\boldsymbol{\alpha})}$ are given by

[TABLE]

where the diagonal matrix and residual vectors are given by

[TABLE]

Using a Gauss-Newton method, the level-set parameters are updated as

[TABLE]

where $\mu_{k}$ is a suitable stepsize and $\boldsymbol{\alpha}^{(0)}$ is a given initial estimate of the shape.

From equation 4, it can be observed that the ability to update the level-set parameters depends on two main factors: 1) The difference between $\mathbf{u}_{0}$ and $u_{1}$, and 2) the derivative of the Heaviside function. Hence, the support and smoothness of $h^{\prime}_{\epsilon}$ plays a crucial role in the sensitivity. More details on the choice of $h_{\epsilon}$ are discussed in section 2.1.

2.0.2 Example

We demonstrate the parametric level-set method on a (binary) discrete tomography problem. We consider the model described in Figure 2(a). For a full-angle case ($0\leq\theta\leq\pi$) with a large number of samples, Figure 2(c) shows that it is possible to accurately reconstruct a complex shape.

2.1 Approximation to Heaviside function

The update of the level-set function primarily depends on the Heaviside function. Various approximations have been mentioned earlier [1]. These approximations suffer from the variational region of Dirac-Delta function near its peak ($\delta|_{x=0}$) which amplifies the gradient disproportionally. This sometimes results in poor updates for the level-set parameter $\boldsymbol{\alpha}$, and hence ruining the reconstructions. To solve this issue, we propose a new formulation of the Heaviside function. We construct the piecewise Dirac-Delta function shown in equation (5):

[TABLE]

This new approximation has been plotted in Figure 3. The above formulation provides mainly 3 benefits: 1) constant sensitivity in the boundary region controlled by parameter $\mu$, 2) a smooth transition part and 3) the compact support.

Definition 2.1**.**

In accordance with the compact approximation of the Heaviside function with width $\epsilon$, a level-set boundary, denoted by $\partial\Omega$, is defined as the set of all the points $\mathbf{x}\in\mathbb{R}^{2}$ satisfying the condition $h_{\epsilon}^{\prime}(\phi(\mathbf{x}))>0$.

Figure 3(c) shows a graphical representation of level-set boundary.

Lemma 2.1**.**

For any smooth and compact approximation of the Heaviside function with finite width $\epsilon$, there exists a relation between level-set boundary and gradient of level-set function, given by $|\delta_{\mathbf{x}}^{T}\nabla\phi(\mathbf{x})|\leq\epsilon$, where, $\delta_{\mathbf{x}}=\max_{x\in\partial\Omega}|\mathbf{x}-\mathbf{x}_{0}|$ and $\mathbf{x}_{0}$ is the point on the level-set.

Proof.

From Taylor series expansion for $\phi(\mathbf{x})$ near the level-set point $\mathbf{x}_{0}$, we get

[TABLE]

$h_{\epsilon}^{\prime}(\phi(\mathbf{x}))>0$ if and only if $|\phi(\mathbf{x})|<\epsilon$. Neglecting higher-order terms, we get $|(\mathbf{x}-\mathbf{x}_{0})^{T}\nabla\phi(\mathbf{x}_{0})|\leq\epsilon$. This implies the above relation.

∎

From the lemma 2.1, it is important to choose the Heaviside width in such a way that the level-set boundary exists on model grid. For simplicity, we crudely approximate the gradient of level-set function using upper and lower bounds [8]. Hence, the heaviside width is represented by

[TABLE]

where $\kappa$ controls the number of gridpoints a level-set boundary can have. This formulation of $\epsilon$ solves the re-initialization issue associated with the level-set method. The steepness ($|\nabla\phi(\mathbf{x})|\gg 1$) of the level-set function in the level-set boundary can be handled by this formulation as well, as it adapts the level-set boundary to global change in level-set function.

3 Joint reconstruction algorithm

Reconstructing both the shape and the background parameter can be cast as a bi-level optimization problem

[TABLE]

where $L$ is of form $[L_{x}^{T}\quad L_{y}^{T}]^{T}$. $L_{x}$ and $L_{y}$ are the second-order finite-difference operators in $x$ and $y$ directions respectively. This optimization problem is separable; it is quadratic in $\mathbf{u}_{0}$ and non-linear in $\boldsymbol{\alpha}$. In order to exploit the fact that the problem has a closed-form solution in $\mathbf{u}_{0}$ for each $\boldsymbol{\alpha}$, we introduce a reduced objective

[TABLE]

The gradient and Hessian of this reduced objective are given by

[TABLE]

where $\overline{\mathbf{u}}_{0}=\operatornamewithlimits{argmin}_{\mathbf{u}_{0}}f(\boldsymbol{\alpha},\mathbf{u}_{0})$ [2].

Using a modified Gauss-Newton algorithm to find a minimizer of $\overline{f}$, leads to the following alternating algorithm

[TABLE]

where the expressions for the gradient and Gauss-Newton Hessian are given by (4). Convergence of this alternating approach to a local minimum of (7) is guaranteed as long as the step-length satisfies the strong Wolfe conditions [16].

The reconstruction algorithm based on this iterative scheme is presented in Algorithm 1.

We use the LSQR method in step 3, with pre-defined maximum iterations and a tolerance value. A trust-region method is applied to compute $\boldsymbol{\alpha}^{(k+1)}$ in step 4 restricting the conjugate gradient to only 10 iterations.

4 Numerical Experiments

The numerical experiments are performed on 4 phantoms shown in figure 4. Each phantom has a constant gray value of parameter 1. For the first two phantoms, the background varies from 0 to 0.5, while for the next two, it varies from 0 to 0.8. In order to avoid inverse crime, the data is generated using a line Kernel, and the forward model uses a Joseph kernel. We use ASTRA toolbox to compute the forward and backward projections [4]. First, we show the results on the noiseless full-view data and later we compare various methods to proposed method in limited-data case with additive gaussian noise of 10 dB SNR.

For the parametric level-set method, we use compactly supported radial basis functions. The basis functions has the form given below:

[TABLE]

RBF nodes are placed on a rectangular grid with the gridspacing 5 times the computational (model) gridspacing. The grid extends to two points outside the model grid to compensate for the background effects. The heaviside width parameter $\kappa$ is set to be 0.01 and the its inclination parameter $\mu$ is set to be 0.1.

The level-set parameter $\boldsymbol{\alpha}$ is optimized using the fminunc package (trust-region algorithm) in MATLAB. A total of 50 iterations are performed for predicting the $\boldsymbol{\alpha}$, while 200 iterations are performed for predicting $\mathbf{u}_{0}(x)$ using LSQR at each step.

4.1 Regularization parameter selection

The reconstruction with the proposed algorithm is influenced by the parameter for Tikhonov regularization. In general, there are various strategies to choose this parameter, e.g., [15]. As our problem formulation deals with the non-linearity in the level-set parameter, application of these kinds of strategies is not clear. Instead we analyze the various residuals, introduced below, with respect to the regularization parameter.

We define three measures (all in the least-squares sense) to quantify the residuals: 1) data residual (DR), determines the data fit between the true data and reconstructed data, 2) model residual (MR), determines the fit between reconstructed model and true model, 3) shape residual (SR), determines the fit between the reconstructed and true anomaly shape. In practice, one can only have a data residual measure to figure out the regularization parameter $\lambda$. From Figure 5, it is evident that there exists a sufficient region of $\lambda$ for which the reconstructions almost stays constant. This region is easily identifiable from the data residual plot for various $\lambda$.

4.2 Benchmark test

For the full-view (benchmark) case, the projection data is generated on $256\times 256$ grid with 256 detectors and 180 projections with $0\leq\theta\leq\pi$. The noise is assumed to be zero in this case. The results on the phantoms with the full-view data are shown in Figure 6. Anomaly geometries in all of these models are reconstructed almost perfectly with the proposed method, although the background has been smoothened out with the tikhonov regularization.

4.3 Limited-angle test

In this case, we use only 5 projections with $\theta$ restricted from [math] to $2\pi/3$. The data is now reduced to almost $3\%$ compared to the benchmark test. We also add Gaussian noise of 10 dB SNR to this synthetic data. To check the performance of the proposed method, we compare it to Total-variation method [4], DART [3] and its modified version for partially discrete tomography, P-DART [13]. A total of 200 iterations were performed with regularization parameter determined from shape residual curve. In DART, the background part was modeled using 20 discrete gray-values between its bounds for model A and B, while 30 discrete gray-values for model C and D. 40 DART iterations were perfomed in each case. For P-DART, a total of 150 iterations were performed.

The results on noisy limited-angle with limited data are presented in Figure 7. The proposed method is able to capture most of the fine details (evident from the shape residual) in the phantoms even with the very limited data with moderate noise. The P-DART method achieves the least amount of data residual in all the cases, but fails to capture the complete geometry of the anomaly.

5 Conclusions and Discussion

We discussed a parametric level-set method for partially discrete tomography. We model such objects as a constant-valued shape embedded in a continuously varying background. The shape is represented using a level-set function, which in turn is represented using radial basis functions. The reconstruction problem is posed as a bi-level optimization problem for the background and level-set parameters. This reconstruction problem can be efficiently solved using a variable projection approach, where the shape is iteratively updated. Each iteration requires a full reconstruction of the background. The algorithm includes some practical heuristics for choosing various parameters that are introduced as part of the parametric level-set method. Numerical experiments on a few numerical phantoms show that the proposed approach can outperform other popular methods for (partially) discrete tomography in terms of reconstruction error. As the proposed algorithm requires repeated full reconstructions, future research is directed at making the method more efficient.

Acknowledgments. This work is part of the Industrial Partnership Programme (IPP) ‘Computational sciences for energy research’ of the Foundation for Fundamental Research on Matter (FOM), which is part of the Netherlands Organisation for Scientific Research (NWO). This research programme is co-financed by Shell Global Solutions International B.V. The second and third authors are financially supported by the Netherlands Organisation for Scientific Research (NWO) as part of research programmes 613.009.032 and 639.073.506 respectively.