A weighted difference of L1 and L2 on the gradient minimization based on alternating direction method for circular computed tomography

Abstract

Iterative reconstruction algorithms for computed tomography (CT) through total variation (TV) regularization can provide accurate and stable reconstruction results. TV minimization is the L₁-norm of gradient-magnitude images and can be regarded as a convex relaxation method to replace the L₀ norm. In this study, a fast and efficient algorithm, which is named a weighted difference of L₁ and L₂ (L₁ - αL₂) on the gradient minimization, was proposed and investigated. The new algorithm provides a better description of sparsity for the optimization-based algorithms than TV minimization algorithms. The alternating direction method is an efficient method to solve the proposed model, which is utilized in this study. Both simulations and real CT projections were tested to verify the performances of the proposed algorithm. In the simulation experiments, the reconstructions from the proposed method provided better image quality than TV minimization algorithms with only 7 views in 180 degrees, which is also computationally faster. Meanwhile, the new algorithm enabled to achieve the final solution with less iteration numbers.

Keywords

Computed tomography image reconstruction difference of L1 and L2 minimization alternating direction method

1 Introduction

The circular scanning configuration is widely used in many computed tomography (CT) applications [1]. The high radiation dose can increase the health risk, which has become a major clinical concern. Two classes of strategies to reduce radiation dose exist: 1) lower the measured X-ray tube current or voltage; and 2) lower the required number of projection views. The first strategy usually leads to noisy projection data at each view with lower signal to noise ratio (SNR), which causes a low reconstruction image quality compared to the data from higher current-measured or voltage-measured scan. The second strategy focuses on the reconstruction from incomplete-view data. To the analytic reconstruction algorithm, the reconstructions with noticeable artifacts will be obtained from the incomplete-view projection data because it does not satisfy the Tuy–Smith condition [2, 3]. On the other hand, the traditional iterative reconstruction algorithms, such as simultaneous algebraic reconstruction technique(SART) [4], provide noisy reconstructions with noticeable artifacts from the incomplete-view projection data as well.

Iterative reconstruction algorithm based on compressive sensing theory [5, 6] can reconstruct images accurately utilizing the sparsity of signals [7, 8]. The sparse of images is generally not widely applicable in industrial and security scanning applications, but it is often the case that images have sparse gradient-magnitude images (GMI) [8]. Method based on the total variation (TV) minimization, which is the L1-norm of GMI, has been widely studied in the reconstruction from few-view and limited-angle data [7 –13]. These methods can be roughly categorized into different classes of optimization-based algorithms, including the adaptive-steepest-descent-projection-onto-convex-sets (ASD-POCS) framework [8], the Chambolle–Pock (CP) framework [14], the split-Bregman (SB) iteration [15], the alternating direction minimization (ADM) method [16, 17], etc.

It should be noted that TV is not the most direct method to describe the GMI’s sparsity. For a digital signal x , the sparsity is directly measured by totaling the non-zero entities in it, i.e., ∥ x ∥ ₀ (the L₀ quasi-norm of x ). The best way to survey the sparsity of GMI is using the GMI’s L₀ quasi-norm. Unfortunately, the optimization for this problem is NP-hard [18]. TV is actually a convex relaxation of the L₀ quasi-norm of GMI replaced by L₁. Meanwhile, several non-convex relaxation methods also exist to replace the L₀ norm. One method is the L_p norm ${∥ x ∥}_{p}^{p} ≜ \sum_{i} {| x_{i} |}^{p}$ (0 < p < 1). Several practical solvers for L_p minimization have been developed [16 , 20]. CT imaging can obtain greater capabilities than those designed with the conventional TV model with reduced datasets [21 –23]. Another method is the difference of L₁ and L₂ norm minimization [24, 25], which gives better sparse recovery results than L₁ and L_p when unsatisfying restricted isometry property (RIP) regime. Additionally, the associated transfer matrix of sparse signals in the transfer domain is less likely to meet the RIP condition for CT image reconstruction [7, 8]. On the basis of this, a weighted difference of anisotropic and isotropic total variation model, i.e., L₁ - αL₂ is developed for image processing [26], where α is chosen according to gradient distributions within [0, 1]. Its convergence is proved and it gives more accurate results than that of L₁.

In this study, we propose a weighted difference of L₁ and L₂ on the gradient minimization for circular CT. It is closer to the L₀ norm of gradient than TV minimization model. The alternating direction method is utilized to iteratively reconstruct a CT image from sparse-view projection measurements by splitting the augmented Lagrangian function of the optimization into two sub-problems. One sub-problem, containing the image constraints term, can be solved by the shrinkage operator after slacking the L₂ of image gradient. The other sub-problem is solved by a local linearization and proximity technique that avoids pseudo-inverse computation of big matrices, which is a quadratic function. The algorithm is tested by simulation and real CT data and shows outstanding capabilities in sampling reduction. The proposed method can provide better image quality when α increases from 0 to 1. Meanwhile, the convergence speed is also faster with the increase of α. The best reconstruction can be provided for the proposed method, with α = 1 than that of α = 0, which equals TV minimization.

The remainder of this paper is organized as follows. Section 2 presents the proposed reconstruction model and the design of the proposed algorithm. In Section 3, numerical simulations and real data experiments are used to validate the performance and capability of the proposed algorithm. In Section 4, a brief discussion is given of the results obtained from Section 3, and our conclusions are given in the end.

2 Methods

2.1 Proposed reconstruction model

The imaging model can be approximated using the following discrete linear system: $p = W f$ (1) where $p \in ℝ^{l}$ is the data vector, $f \in ℝ^{s}$ is the object vector, and $W \in ℝ^{l \times s}$ is the system matrix. The element in system matrix W represents an intersection of an X-ray detected with a voxel in the image array. With the linear equation expressed above, image reconstruction is performed to solve equation (1). However, even under ideal and consistent condition, finding the solution is not an easy task. This is mainly for two reasons, one of which is that the projection matrix W is tremendously huge thus computing its (pseudo) inverse is computationally intractable. Another reason is that W is often severely ill-conditioned having huge condition number resulting in numerical computing instability. For sparse-view projection data, equation (1) is highly ill-posed, and conventionally an optimization program based on CS is established to solve this problem. x is generally not sparse in CT, but its sparse GMI is often available to be set into optimization program. Here, we take the gradient operator as ∇, so the direct approach is L₀ optimization: $f^{*} = arg min {∥ \nabla f ∥}_{0} s . t . p = W f, f \geq 0 .$ (2)

Several viable alternatives are available because minimizing L₀ norm is NP-hard [18]. Greedy methods (matching pursuit [27], orthogonal matching pursuits (OMP) [28], and regularized OMP (ROMP) [29]) work well if the dimension N is not too large. Another strategy is relaxing the L₀ norm by a continuous sparsity promoting penalty functions. Convex relaxation uniquely selects the L₁ norm, which is the TV minimization:

$\begin{matrix} f^{*} = arg min {∥ f ∥}_{tv} s . t . p = W f, f \geq 0, \\ or & f^{*} = arg min {∥ \nabla f ∥}_{1} s . t . p = W f, f \geq 0 . \end{matrix}$ (3)

Here, we select a novel non-convex relaxation approach, called weighted difference of L₁ and L₂: $f^{*} = arg min {∥ \nabla f ∥}_{1} - α {∥ \nabla f ∥}_{2} s . t . p = W f, f \geq 0,$ (4) where α ∈ [0, 1]. The weighted difference of L₁ and L₂ is a closer approximation to the L₀ than the L₁. Meanwhile, it changes gradually from L₁ to L₀ with the increase of α from 0 to 1. It can be observed visually in Fig. 1.

Fig.1

Contours of three sparsity metrics. The level curves of approach the x and y axes as the values decrease, hence promoting sparsity. Panel (a) shows the contour of L₁, panel (b) shows the contour of L₁-(1/2)L₂, and panel (c) shows the contour of L₁-L₂.

The model can also be written as $f^{*} = arg min {∥ z ∥}_{1} - α {∥ z ∥}_{2} s . t . p = W f, \nabla f = z, f \geq 0,$ (5)

2.2 ADMM reconstruction for proposed model

To solve the proposed reconstruction model (5) efficiently, the alternating minimization method (ADM) algorithm is utilized in this paper. For constrained optimization, an influential class of method class seeks for the minimizer or maximizer by approaching the original constrained problem by a sequence of unconstrained sub-problems. We apply the augmented Lagrange (AL) function, which could be regarded as a combination of the Lagrangian and quadratic penalty functions, to convert formula (5) to an unconstrained form: $\begin{matrix} min L_{A} (f, z; ρ, λ) & = & ρ^{T} (\nabla f - z) + \frac{β_{1}}{2} {∥ \nabla f - z ∥}_{2}^{2} + λ^{T} (W f - p) \\ + \frac{β_{2}}{2} {∥ W f - p ∥}_{2}^{2} + {∥ z ∥}_{1} - α {∥ z ∥}_{2}, \end{matrix}$ (6a) $or min L_{A} (f, z; ρ, λ) = \frac{β_{1}}{2} {∥ \nabla f - z + \frac{ρ}{β_{1}} ∥}_{2}^{2} + \frac{β_{2}}{2} {∥ W f - p + \frac{λ}{β_{2}} ∥}_{2}^{2} + {∥ z ∥}_{1} - α {∥ z ∥}_{2},$ (6b)

where β₁ and β₂ are the penalty coefficients, and ρ and λ are the Lagrange multipliers. The non-negative constraint in formula (5) can be maintained by setting all the negative entities to 0 during algorithm iterations. The ADM is used here to solve the minimization of formula (6) efficiently, which splits the optimization into two sub-problems, L _z (f) and L _f (z), alternatingly. So the solution of the proposed reconstruction model can be solved as following iterations: $f^{k + 1} \in \underset{f}{arg min} L_{A} (f, z^{k}; ρ^{k}, λ^{k}),$ (7a) $z^{k + 1} \in \underset{z}{arg min} L_{A} (f^{k + 1}, z; ρ^{k}, λ^{k}),$ (7b) $ρ^{k + 1} = ρ^{k} + β_{1} (\nabla f^{k + 1} - z^{k + 1}),$ (7c) $λ^{k + 1} = λ^{k} + β_{2} (W f^{k + 1} - p) .$ (7d)

Fixing the variable z, the sub-problem L _z ( f ) can be written as $min L_{z} (f) = \frac{β_{1}}{2} {∥ \nabla f - z + \frac{ρ}{β_{1}} ∥}_{2}^{2} + \frac{β_{2}}{2} {∥ W f - p + \frac{λ}{β_{2}} ∥}_{2}^{2} .$ (8)

L _z (f) is a quadratic function, and its gradient is $d_{k} (f) = β_{1} \nabla^{T} (\nabla f - z + \frac{ρ}{β_{1}}) + β_{2} W^{T} (W f - p + \frac{λ}{β_{2}}) .$ (9)

Let d _k ( f ) =0, the exact minimizer of L _z ( f ) is $f^{*} = {(β_{1} \nabla^{T} \nabla + β_{2} W^{T} W)}^{+} (β_{1} \nabla^{T} z - \nabla^{T} ρ + β_{2} W^{T} p - W^{T} λ)$ (10) where M+ denotes the Moore–Penrose pseudo inverse of matrix M. However, the calculation of pseudo inverse for the large matrix β₁ ∇ ^T ∇ + β₂W^TW is too extensive. Here, we adopt an inexact strategy to tackle the sub-problem of f minimization. The fidelity term ∥Wf - p + λ / - β₂ ∥ ² in Equation (8) can be linearized at the current point f ^k. The proximal form is as follows ${∥ W f - p + \frac{λ}{β_{2}} ∥}^{2} \approx {∥ W f^{k} - p + \frac{λ}{β_{2}} ∥}^{2} + 2 g_{k}^{T} (f - f^{k}) + \frac{1}{τ} {∥ f - f^{k} ∥}^{2},$ (11) where g _k = W^T (Wf^k - p + λ / - β₂) is the gradient of ∥Wf - p + λ / - β₂ ∥ ² at the current point of f ^k, τ > 0. Thus, the sub-problem L _z (f) can be converted into the following form

$\begin{matrix} min L_{z} (f) & = & \frac{β_{1}}{2} {∥ \nabla f - z + \frac{ρ}{β_{1}} ∥}_{2}^{2} + \frac{β_{2}}{2} {∥ W f^{k} - p + \frac{λ}{β_{2}} ∥}_{2}^{2} \\ + β_{2} g_{k}^{T} (f - f^{k}) + \frac{β_{2}}{2 τ} {∥ f - f^{k} ∥}^{2} . \end{matrix}$ (12)

Setting its gradient to 0 results in $(β_{1} \nabla^{T} \nabla + \frac{β_{2}}{τ}) f = β_{1} \nabla^{T} z - \nabla^{T} ρ + \frac{β_{2}}{τ} f^{k} - β_{2} g_{k} .$ (13)

On the periodic definition condition, the matrix form of ∇^T∇ is block circulant. Therefore, ∇^T∇ can be diagonalized using FFT transform as follows: $Λ = β_{1} 𝔽 \nabla^{T} \nabla 𝔽^{H} + \frac{β_{2}}{τ} I,$ (14) where $𝔽^{H}$ denotes the Hermite form of $𝔽$ . The matrix Λ is diagonal and all the entities in the principal diagonal are all positive which makes it easy to compute its inverse. Thus, the solution can be easily computed as follows: $f^{*} = 𝔽^{- 1} {Λ^{- 1} 𝔽 (β_{1} \nabla^{T} z - \nabla^{T} ρ + \frac{β_{2}}{τ} f^{k} - β_{2} g_{k})} .$ (15)

Fixing the variable f , the sub-problem L _f (z) can be written as follows: $min L_{f} (z) = {∥ z ∥}_{1} - α {∥ z ∥}_{2} + \frac{β_{1}}{2} {∥ \nabla f - z + \frac{ρ}{β_{1}} ∥}_{2}^{2} .$ (16)

Unfortunately, L _f (z) is a non-convex problem, and it cannot be solved directly. Note that ∥z ∥ ₂ is differentiable with gradient $\begin{matrix} z / - {∥ z ∥}_{2} & \forall \end{matrix} z \neq 0$ and that 0 ∈ ∂ ∥ z ∥ ₂ for z = 0. Here, we solve it by an approximate solution, as follows: $min L_{f} (z) = {\begin{matrix} \begin{matrix} \frac{β_{1}}{2} {∥ \nabla f - z + \frac{ρ}{β_{1}} ∥}_{2}^{2} + {∥ z ∥}_{1} & if z^{k} = 0, \end{matrix} \\ \begin{matrix} \frac{β_{1}}{2} {∥ \nabla f - z + \frac{ρ}{β_{1}} ∥}_{2}^{2} + {∥ z ∥}_{1} - 〈 z, α \frac{z^{k}}{{∥ z^{k} ∥}_{2}} 〉 & otherwise, \end{matrix} \end{matrix}$ (17) where k represents the current number of iterations. The above equation can be efficiently solved by the shrinkage operator, and the exact solution is as follows: $z^{*} = {\begin{matrix} \begin{matrix} max {| \nabla f + \frac{ρ}{β_{1}} | - \frac{1}{β_{1}}, 0} • sgn (\nabla f + \frac{ρ}{β_{1}}) & if z^{k} = 0, \end{matrix} \\ \begin{matrix} max {| \nabla f + \frac{α z^{k} / - {∥ z^{k} ∥}_{2} + ρ}{β_{1}} | - \frac{1}{β_{1}}, 0} • sgn (\nabla f + \frac{α z^{k} / - {∥ z^{k} ∥}_{2} + ρ}{β_{1}}) & otherwise . \end{matrix} \end{matrix}$ (18)

In conclusion, the weighted difference of L₁ and L₂ on the gradient minimization with ADM is proposed, and the pseudocode of the proposed algorithm is as follows:

Algorithm 1: Pseudocode for k steps ADM iteration scheme for proposed algorithm
1: β₁ ← (β₁) ₀ > 0, β₂ ← (β₂) ₀ > 0 ; k ← 0 ; α ← α₀ ∈ [0, 1]
2: f ← f ₀ ; z ← z₀ ; ρ ← 0, λ ← 0
3: While “not converged” Do
4: $f^{k + 1} \leftarrow 𝔽^{- 1} {Λ^{- 1} 𝔽 (β_{1} \nabla^{T} z^{k} - \nabla^{T} ρ^{k} + \frac{β_{2}}{τ} f^{k} - β_{2} g_{k})}, g_{k} = W^{T} (W f^{k} - p + λ^{k} / - β_{2}),$
5: ifz^k = 0
6: $z^{k + 1} \leftarrow max {\| \nabla f^{k + 1} + \frac{ρ^{k}}{β_{1}} \| - \frac{1}{β_{1}}, 0} • sgn (\nabla f^{k + 1} + \frac{ρ^{k}}{β_{1}}),$
7: else
8: $z^{k + 1} \leftarrow max {\| \nabla f^{k + 1} + \frac{α z^{k} / - {∥ z^{k} ∥}_{2} + ρ^{k}}{β_{1}} \| - \frac{1}{β_{1}}, 0} • sgn (\nabla f^{k + 1} + \frac{α z^{k} / - {∥ z^{k} ∥}_{2} + ρ^{k}}{β_{1}}),$
9: end if
10: ρ ^k+1 ← ρ ^k + β₁ (∇ f ^k+1 - z^k+1),
11: λ ^k+1 ← λ ^k + β₂ (Wf^k+1 - p),
12: k ← k + 1,
13: end Do
14: return_bfitf.

While under the influences of noise, the observation equation is no longer p = Wf. In this paper, a simple and direct method is introduced to allow for the presence of noise. Assume the noise is denoted by e , then the observation can be modified by adding e in (5), then getting $\begin{matrix} f^{*} = arg min {∥ z ∥}_{1} - α {∥ z ∥}_{2} \\ s . t . {\begin{matrix} p - W f = e (∥ e ∥ \leq δ), \\ \nabla f = z, \\ f \geq 0, \end{matrix} \end{matrix}$ (19) where the positive parameter δ is an estimate of the noise level in the data. So the AL function of the above problem is $min L_{A} (f, z; ρ, λ) = \frac{β_{1}}{2} {∥ \nabla f - z + \frac{ρ}{β_{1}} ∥}_{2}^{2} + \frac{β_{2}}{2} {∥ W f - p + e + \frac{λ}{β_{2}} ∥}_{2}^{2} + {∥ z ∥}_{1} - α {∥ z ∥}_{2} .$ (20)

The implication of parameters in (20) is the same to that in (6). Compare the above function to that in (6), it can be easily found that the solvations of f , z in (20) are almost the same to f , z except for some notations. Here we omit the derivation of these two sub-problems, and the final expressions will be shown in Algorithm 2.

Excepting updating variables f , z and Lagrange multipliers ρ , λ , the remaining vector r should be updated as well as $e^{k + 1} = p - W f^{k + 1} .$ (21)

However, e ^k+1 may fail to satisfy the constraint ∥ e ∥ ≤ δ. To cope with this defect, we should force it to go back to the feasible set by a projection $e^{k + 1} \leftarrow P_{U} [p - W f^{k + 1}],$ (22) where $U = {x \in ℝ^{m} | | | x | |_{2} \leq δ}$ and P_U [·] denotes the projection operator onto the ball U under Euclidean norm.

Again, thus, a new working algorithm for allowing for the noise can be generated as follow

Algorithm 2: Pseudocode for k steps ADM iteration scheme for proposed algorithm with the noise data
1: β₁ ← (β₁) ₀ > 0, β₂ ← (β₂) ₀ > 0 ; k ← 0 ; α ← α₀ ∈ [0, 1]
2: f ← f ₀ ; z ← z₀ ; ρ ← 0, λ ← 0, e ← 0
3: While “not converged” Do
4: $f^{k + 1} \leftarrow 𝔽^{- 1} {Λ^{- 1} 𝔽 (β_{1} \nabla^{T} z^{k} - \nabla^{T} ρ^{k} + \frac{β_{2}}{τ} f^{k} - β_{2} g_{k})}, g_{k} = W^{T} (W f^{k} - p + λ^{k} / - β_{2} + e^{k}),$
5: ifz^k = 0
6: $z^{k + 1} \leftarrow max {\| \nabla f^{k + 1} + \frac{ρ^{k}}{β_{1}} \| - \frac{1}{β_{1}}, 0} • sgn (\nabla f^{k + 1} + \frac{ρ^{k}}{β_{1}}),$
7: else
8: $z^{k + 1} \leftarrow max {\| \nabla f^{k + 1} + \frac{α z^{k} / - {∥ z^{k} ∥}_{2} + ρ^{k}}{β_{1}} \| - \frac{1}{β_{1}}, 0} • sgn (\nabla f^{k + 1} + \frac{α z^{k} / - {∥ z^{k} ∥}_{2} + ρ^{k}}{β_{1}}),$
9: end if
10: e ^k+1 ← P_U [p - Wf^k+1],
11: ρ ^k+1 ← ρ ^k + β₁ (∇ f ^k+1 - z^k+1),
12: λ ^k+1 ← λ ^k + β₂ (Wf^k+1 - p + e^k+1),
13: k ← k + 1,
14: end Do
15: return_bfitf.

2.3 Stopping criterion

It is important to identify an effective stopping criterion for the iteration sequence (7) for algorithm 1. And the Karush-Kuhn-Tucker (KKT) conditions [30] are satisfied by optimal solutions of constrained optimization problems. Here, we take ( f ^*, z^*, ρ ^*, λ ^*) to be the optimal solution of the original minimization problem (5). Thus, we can get the KKT conditions as follow $\nabla^{T} ρ^{*} + W^{T} λ^{*} = 0,$ (23a) $\partial {∥ z ∥}_{1} | z = z^{*} - α \partial {∥ z ∥}_{2} | z = z^{*} - ρ^{*} = 0,$ (23b) $z^{*} - \nabla f^{*} = 0,$ (23c) $p - W f^{*} = 0,$ (23d)

where ∂ H denotes the subdifferential of a convex function H . Based on (23), the proximity of the iteration sequence ( f ^k, z^k, ρ ^k, λ ^k) generated by (7) to an optimal solution of the original minimization problem (5) can be verified. Due to the approximate treatment of sub-problem L _z ( f ) and L _f (z), the optimality condition of (12) and (17) are considered rather than (8) and (16). Combined with (7c) and (7d), we can get $\nabla^{T} ρ^{k + 1} + W^{T} λ^{k + 1} = β_{1} \nabla^{T} (z^{k} - z^{k + 1}) + ((β_{2} / - τ) I - β_{2} W^{T} W) (f^{k} - f^{k + 1}),$ (24a) $\partial {∥ z ∥}_{1} | z = z^{k + 1} - α \partial {∥ z ∥}_{2} | z = z^{k + 1} - ρ^{k + 1} = α (z^{k} / - {∥ z^{k} ∥}_{2} - \partial {∥ z ∥}_{2} | z = z^{k + 1}),$ (24b)

where $\partial {∥ z ∥}_{2} | z = z^{k + 1} = {\begin{matrix} 0 & if z^{k + 1} = 0, \\ z^{k + 1} / - {∥ z^{k + 1} ∥}_{2} & otherwise . \end{matrix}$ (25)

And the detailed derivation of (24) is shown in Appendix 1. Thus, the residuals for the feasibility conditions can be viewed as $r_{f}^{k + 1} = β_{1} \nabla^{T} (z^{k} - z^{k + 1}) + ((β_{2} / - τ) I - β_{2} W^{T} W) (f^{k} - f^{k + 1}),$ (26a) $r_{z}^{k + 1} = α (z^{k} / - {∥ z^{k} ∥}_{2} - \partial {∥ z ∥}_{2} | z = z^{k + 1}) .$ (26b)

Under ideal conditions, the mathematical convergence conditions of Algorithm 1 should be ${\bar{r}}_{f}^{n} \to 0$ and ${\bar{r}}_{z}^{n} \to 0$ as n → ∞, where ${\bar{r}}_{f}^{n} = | r_{f}^{n} / - r_{f}^{1} |$ and ${\bar{r}}_{z}^{n} = | r_{z}^{n} / - r_{z}^{1} |$ [31]. However, they are unachievable in practical, due to limited computer precision and computation time. So we can use $max {{\bar{r}}_{f}^{k + 1}, {\bar{r}}_{z}^{k + 1}}$ to measure the k+1 iteration result of the minimization problem (7). The iteration sequence in algorithm 1 can be terminated when the residuals for the feasibility conditions are small enough in a prescribed sense.

For algorithm 2, the KKT conditions has a little difference because of the addition of noise e . The formula (23a-c) have no change, and the formula (23d) is changed into $p - W f^{*} = e, (∥ e ∥ \leq δ),$ (27)

While it has no effect on the residuals for f and z, leading to the same ${\bar{r}}_{f}^{k + 1}$ and ${\bar{r}}_{z}^{k + 1}$ as formula (26). Thus, the stopping criterion of Algorithm 2 is the same as Algorithm 1.

2.4 Parameter selections

For the proposed method, the selections of parameters have huge impact on the convergence. Parameter β₁ and β₂ are used to balance the quadratic term and penalty function. As a matter of experience, they are better to be set more than 0 and less than 10⁴. Between them, they are also used to balance the weighting of influence by two constraint terms. Generally, they are in the same order of magnitudes. Parameter τ is a positive number, and it is usually set in [1, 2]. In this paper, it is set 1 for all of the experiments. For noisy data, parameter δ in algorithm 2 should be considered. If the noise e ₀ is knowable, we set δ = ∥ e₀ ∥ ₂ or less than that. If the noise e ₀ is not knowable, it is always better to set δ with a small value close to 0.

3 Experimental results

The algorithm was performed with the use of both numerical simulation and real data. All experiments were implemented with MATLAB programming language on an HP workstation with Intel^® Xeon^® CPU E5-2620 v3 @2.40 GHz×2 and 64 GB RAM on a 64-bit Windows 7 system.

3.1 Simulation experiments

In this experiment, a fan beam geometry is set to simulate the data acquisition of CT imaging. The typical Shepp–Logan phantom with a discretization of 256×256 pixels is utilized as the standard image. The distance from the X-ray source to the iso-center is 164.86 mm and the distance from the X-ray source to the center of the line detector is 1052.00 mm. The number of detector bins is 512 with a size of 0.074 mm for each.

In the first experiment, seven projections are evenly acquired within 180° along the circular trajectory. The parameter α is discussed first. We set α from 0 to 1 with interval 0.1. The model equals TV minimization when α = 0, and it satisfies the model of L₁- L₂ in [25] when α = 1. Here we set β₁ = 66 and β₂ = 62 for all of the tests in experiment 1. Through the stopping criterion strategy in 2.3 chapter, we consider that the reconstruction result can be the solution when $max {{\bar{r}}_{f}^{k + 1}, {\bar{r}}_{z}^{k + 1}} < 10^{- 5}$ . Figure 2 shows the curves of the value $max {r_{f}^{k + 1}, r_{z}^{k + 1}}$ for the proposed algorithm with each value of parameter α. We can find that the stopping iteration number is decreasing with the increasing of α. To guarantee the final solutions are convergent, we set the stopping iteration number is larger than the largest one above, which is 2500 here. Figure 3 presents the model of the Shepp–Logan phantom and the reconstruction results obtained via the proposed method with α = 0, 0.5, 1 from the top to the bottom with a grayscale window [0.195, 0.205]. Moreover, from columns 2 to 4 are reconstructions at iteration numbers 1200, 1600, 2000, respectively. Moreover, the root-mean-square-error (RMSE) curves for the proposed algorithm with each value of parameter α are plotted in Fig. 4.

Fig.2

Curves of value $max {{\bar{r}}_{f}^{k + 1}, {\bar{r}}_{z}^{k + 1}}$ for the proposed algorithm with each value of parameter α based on the original Shepp–Logan phantom.

Fig.3

Original phantom and seven-view reconstructions by the proposed algorithm with α = 0, 0.5, 1. Panel (a) shows the original phantom, panel (b) shows results at 1200 iterations, panel (c) shows results at 1600 iterations, and panel (d) shows results at 2000 iterations. Its grayscale window is [0.195, 0.205].

Fig.4

The comparisons of RMSE curves for the proposed algorithm with different α values for seven views of projections based on the original Shepp–Logan phantom.

Figure 3 shows some results in special parameter α that it can provide better reconstructions using less iterations for the larger α. With the iteration increasing, the final results can have a high image quality for all the α in [0, 1]. For a more detailed result, the curves of RMSE for each value of α in Fig. 4 indicate that reconstructions by the proposed method have faster error reduction speed with increasing α. Table 1 shows that the value of RMSE is smaller when the value of α increases at each iteration number. Although some fluctuations exist in the back iterations, the final results improve with increasing α. In other words, the final reconstruction is the best one when α = 1 with the smallest value of RMSE.

Table 1

The numerical comparisons at different iterations of each α for the proposed algorithm with seven views

RMSEs	100 iter.	400 iter.	800 iter.	1200 iter.	1600 iter.	2000 iter.	2500 iter
α = 0	0.08475	0.05462	0.02891	0.01116	0.00292	2.12E–4	5.25E-5
α = 0.2	0.08471	0.05394	0.02674	0.00789	9.08E–4	4.91E–5	4.06E-5
α = 0.4	0.08467	0.05328	0.02449	0.00489	4.29E–4	4.75E–5	2.96E-5
α = 0.6	0.08463	0.05259	0.02213	0.00290	2.76E–4	3.48E–5	3.81E-5
α = 0.8	0.08460	0.05191	0.01973	0.00176	1.17E–4	6.07E–5	4.28E-5
α = 1.0	0.08456	0.05126	0.01736	0.00113	9.76E–5	1.17E–5	2.74E-5

In further studies, the parameter α is discussed with two special values, 0 and 1, based on the results in the first experiment. Other values of α from 0 to 1 cause reconstruction performance between that of α = 0 and 1. In the second experiment, we use the low-contrast Shepp–Logan phantom generated by MATLAB code phantom(’Shepp-Logan’,256), and the system parameters is invariant to that in the first experiment. We also use $max {{\bar{r}}_{f}^{k + 1}, {\bar{r}}_{z}^{k + 1}} < 10^{- 5}$ to be the convergence condition in this experiment. Figure 5 shows the curves of the value $max {{\bar{r}}_{f}^{k + 1}, {\bar{r}}_{z}^{k + 1}}$ for the proposed algorithm with α = 0 and 1 in this experiment. After calculating the stopping criterion, the final stopping iteration number is set 3500. The results for each parameter α at the iteration number of 3500 are shown in Fig. 6 with a grayscale window [0.0195, 0.0205] and the curves of RMSE for each α are plotted in Fig. 7. The results are similar to that of the original Shepp–Logan phantom in the first experiment. The proposed algorithm with α = 1 can provide a better reconstruction with the RMSE of 7.02E–6 than TV minimization when α = 0 with RMSE of 9.09E–6. Meanwhile, it can provide the available reconstruction using less iteration numbers.

Fig.5

Curves of value $max {{\bar{r}}_{f}^{k + 1}, {\bar{r}}_{z}^{k + 1}}$ for the proposed algorithm based on the low-contrast Shepp–Logan phantom. Parameter α = 0 in the left and α = 1 in the right.

Fig.6

Original phantom and seven-view reconstructions with different α values at an iteration number of 3500 for low-contrast Shepp–Logan phantom. Panel (a) shows the original phantom, panel (b) shows the result with α = 0, and its RMSE is 9.09E–6, and panel (c) shows the result with α = 1, and its RMSE is 7.02E–6. Its grayscale window is [0.0195, 0.0205].

Fig.7

The comparisons of RMSE curves for the proposed algorithm with different α values for seven views of projections for low-contrast Shepp–Logan phantom.

3.2 Simulation experiments with Poisson random noise

Reconstructions with noisy projections are conducted and demonstrated in this subsection. The Poisson noise is added into the simulated projections. Although the noise level in the experiment is not very high compared with the original projections, reconstruction with very few views is still very difficult. In this experiment, the reconstructions under seven views with α = 0 and 1 are carried out for the proposed algorithms. Here we set β₁ = 128 and β₂ = 100 for all of the tests in this experiment. In this experiment, we use algorithm 2 to test the performance of the proposed method with α = 0 and 1. In this test, we consider that the reconstruction result can be the solution when $max {{\bar{r}}_{f}^{k + 1}, {\bar{r}}_{z}^{k + 1}} < 10^{- 4}$ . The curves of the value $max {{\bar{r}}_{f}^{k + 1}, {\bar{r}}_{z}^{k + 1}}$ for the proposed algorithm with α = 0 and 1 in this experiment is shown in Fig. 8. Figure 9 shows the original phantom and reconstruction results at an iteration number of 2100 with α = 0 and 1. The result at α = 0 is not convergent, but the result at α = 1 can be the final solution. Furthermore, Fig. 10 shows the original phantom and reconstruction results at an iteration number of 3000 with α = 0 and 1. The results for α = 1 is better than that of α = 0 as well. Meanwhile, the reconstruction at 2100 iterations is better than that of 3000 iterations when α = 1, because the influence of noise will be amplified with increasing iterations after converging. In conclusion, the proposed algorithm can provide better reconstruction with α = 1 than α = 0, and its error reduction speed is much faster, as shown in Fig. 11.

Fig.8

Curves of value $max {{\bar{r}}_{f}^{k + 1}, {\bar{r}}_{z}^{k + 1}}$ for the proposed algorithm based on the original Shepp–Logan phantom adding Poisson random noise. Parameter α = 0 in the left and α = 1 in the right.

Fig.9

Original phantom and seven-view reconstructions with different α values at an iteration number of 2100 for original Shepp–Logan phantom adding Poisson random noise. Panel (a) shows the original phantom, panel (b) shows the result with α = 0, and its RMSE is 2.4E–3, and panel (c) shows the result with α = 1, and its RMSE is 1.08E–4. Its grayscale window is [0.195, 0.205].

Fig.10

Original phantom and seven-view reconstructions with different α values at an iteration number of 3000 for original Shepp–Logan phantom adding Poisson random noise. Panel (a) shows the original phantom, panel (b) shows the result with α = 0, and its RMSE is 2.54E–4, and panel (c) shows the result with α = 1, and its RMSE is 1.38E–4. Its grayscale window is [0.195, 0.205].

Fig.11

The comparisons of RMSE curves for the proposed algorithm with different α values for seven views of projections for the original Shepp–Logan phantom.

3.3 Experiments with real CT projections

In this Section, real data verifications are conducted to test the proposed Algorithm 2. Projections are acquired from an adult male rabbit in vivo on a CT imaging system. The distance from the source to the rotation center is 502.8 mm, and the distance from the source to the center of the detector is 1434.7 mm. The flat panel detector has 1024 bins with a pixel size of 0.381 mm. The size of the reconstructed image is 800×800 with 0.1335 mm×0.1335 mm for each voxel. A total of 360 projections are acquired evenly within the range of 360°.

We first perform the reconstruction with the real CT projections using the standard FBP algorithm. Reconstruction with full views and sparse views (60 projections evenly chosen with angular step of 3° within 180°) by FBP are conducted and the typical results are shown in the two left images in Fig. 13. In this test, we set β₁ = 10 and β₂ = 10 for all of the tests in this experiment. Algorithm 2 is used to test the performance of the proposed method with α = 0 and 1 as well. For iterative algorithms under real CT projections, the iteration should be stopped at an appropriate number, thereby avoiding the influence of noise. Thus, we consider that the reconstruction result can be the solution when $max {{\bar{r}}_{f}^{k + 1}, {\bar{r}}_{z}^{k + 1}} < 10^{- 3}$ . The curves of the value $max {{\bar{r}}_{f}^{k + 1}, {\bar{r}}_{z}^{k + 1}}$ for the proposed algorithm with α = 0 and 1 in this experiment is shown in Fig. 12. Here, the reconstructions by the proposed algorithm with α = 0 and 1 with 60 views evenly chosen within 180° at iteration number 130 are shown in the two right images in Fig. 13. We set the gray scale window to [0.008, 0.14].

Fig.12

Curves of value $max {{\bar{r}}_{f}^{k + 1}, {\bar{r}}_{z}^{k + 1}}$ for the proposed algorithm based on the real data. Parameter α = 0 in the left and α = 1 in the right.

Fig.13

Reconstructions for the rabbit by FBP and the proposed algorithm at iteration number of 130. Panel (a) shows the result of FBP with full views, panel (b) shows the result of FBP with 60 views evenly chosen within 180°, panel (c) shows the result of the proposed algorithm at α = 0 with 60 views evenly chosen within 180°, and panel (d) shows the result of the proposed algorithm at α = 1 with 60 views evenly chosen within 180°. Its grayscale window is [0.008, 0.14].

Images reconstructed from sparse views by FBP suffer severely from the streak artifacts and several tiny structures are degraded, which prove that the analytical reconstruction algorithm is not applied to sparse view projections. No exact image for the real CT exists, even when using the FBP under full views. The result by the full view FBP can be a reference to the result by iterative algorithm. The reconstructions by the proposed algorithm, as shown in the two right images in Fig. 13, are close to the result by full view FBP. For detailed results, the ROI of these four reconstructions are shown in Fig. 14. The reconstruction by the proposed method with α = 1 is closer to the full view FBP than the reconstruction by the proposed method with α = 0 for the tiny structures of bones.

Fig.14

ROI of reconstructions for the rabbit by FBP and the proposed algorithm at iteration number of 130. Panel (a) shows the ROI result of FBP with full views, panel (b) shows the ROI result of FBP with 60 views evenly chosen within 180°, panel (c) shows the ROI result of the proposed algorithm at α = 0 with 60 views evenly chosen within 180°, and panel (d) shows the ROI result of the proposed algorithm at α = 1 with 60 views evenly chosen within 180°. Its grayscale window is [0.008, 0.14].

4 Discussions and conclusions

TV-based methods have been widely used in CT image reconstruction under sparse-view projection measurements. While it is not the best way to describe the sparsity, but the L₀ norm. This paper proposed an efficient reconstruction algorithm named a weighted difference of L₁ and L₂ on the gradient minimization for circular CT. It is closer to the L₀ norm than the TV minimization method, that is the better description of image sparsity. This new method is derived by ADM and can achieve promising results using less iteration numbers under very few views of projections compared with conventional TV minimization method. Considering the noise data, we improve the reconstruction model, and propose a new reconstruction algorithm. Moreover, we provide the stopping criterion through the KKT conditions to guarantee that the reconstructed image can be considered as the final solution. The proposed method can be regarded as the L₁ - αL₂ norm of gradient, and the TV minimization method is the L₁ norm of gradient. Though the proposed method appends only one term compared with TV minimization method, they are different in essence. TV minimization can be regarded as a convex relaxation of L₀ norm of gradient, whereas the proposed method is a non-convex relaxation of L₀ norm of gradient. The proposed method has a better performance in reconstruction under very few views of projections because it is has better description of image sparsity than the TV minimization method.

In the simulation experiments, the stopping iteration number and RMSE curves show that the proposed method can provide better reconstruction results than traditional TV minimization algorithms whether noise-free data or noise data, and requiring less iteration number with increasing α from 0 to 1. The best performance of the proposed method is obtained at α = 1. In the real data reconstruction, we test the proposed method by sparse-view projection data from the skull of adult male rabbit on a CT imaging system. The reconstructions indicate that the proposed method can preserve many fine details that many tissues of structural information are visible, which can be as good as the full-scan FDK reconstruction. While, some of them are lost in the TV reconstruction.

Footnotes

Appendix 1

The optimality condition of (12) and (17) for Algorithm 1.

For the k iteration of (7), we can get following expressions from (12) and (17) (28) $\nabla^{T} (ρ^{k} + β_{1} (\nabla f^{k + 1} - z^{k})) + W^{T} (λ^{k} + β_{2} (W f^{k} - p)) + \frac{β_{2}}{τ} (f^{k + 1} - f^{k}) = 0$ (29) $- ρ^{k} + β_{1} (z^{k + 1} - \nabla f^{k + 1}) + \partial {∥ z ∥}_{1} | z = z^{k + 1} - α \frac{z^{k}}{{∥ z^{k} ∥}_{2}} = 0 .$

Then, combined with (7c) and (7d), they are converted into (30a) $\nabla^{T} ρ^{k + 1} + W^{T} λ^{k + 1} + β_{1} \nabla^{T} (z^{k + 1} - z^{k}) + β_{2} W^{T} W (f^{k} - f^{k + 1}) + \frac{β_{2}}{τ} (f^{k + 1} - f^{k}) = 0,$ (30b) $\partial {∥ z ∥}_{1} | z = z^{k + 1} - α \partial {∥ z ∥}_{2} | z = z^{k + 1} - ρ^{k + 1} + α \partial {∥ z^{k + 1} ∥}_{2} - α z^{k} / - {∥ z^{k} ∥}_{2} = 0 .$

So we can get following expressions (31a) $\nabla^{T} ρ^{k + 1} + W^{T} λ^{k + 1} = β_{1} \nabla^{T} (z^{k} - z^{k + 1}) + ((β_{2} / - τ) I - β_{2} W^{T} W) (f^{k} - f^{k + 1}),$ (31b) $\partial {∥ z ∥}_{1} | z = z^{k + 1} - α \partial {∥ z ∥}_{2} | z = z^{k + 1} - ρ^{k + 1} = α (z^{k} / - {∥ z^{k} ∥}_{2} - \partial {∥ z ∥}_{2} | z = z^{k + 1}),$

where (32) $\partial {∥ z ∥}_{2} | z = z^{k + 1} = {\begin{matrix} 0 & if z^{k + 1} = 0 \\ z^{k + 1} / - {∥ z^{k + 1} ∥}_{2} & otherwise \end{matrix}$

Acknowledgments

This work is supported by the National Natural Science Foundation of China (Nos. 61372172, 61601518). We also want to thank the anonymous referees for giving us useful comments and suggestions to improve this paper.

References

Buzug

T.M.

Computed Tomography: From Photon Statistics to Modern Cone-Beam CT, (Springer-Verlag, Berlin Heidelberg, 2008).

Tuy

H.K.

, Inversion formula for cone-beam reconstruction, Siam J Appl Math 42 (1983), 546–552.

Smith

B.D.

, Image reconstruction from cone-beam projections: Necessary and sufficient conditions and reconstruction methods, IEEE Trans Med Imag 4 (1985), 14–25.

Andersen

A.H.

and Kak

A.C.

, simultaneous algebraic reconstruction technique (SART): A superior implementation of the ART algorithm, Ultrasonic Imaging 6 (1984), 81–94.

Candes

, Romberg

and Tao

, Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information, IEEE Trans Inf Theory 52 (2006), 489–509.

Candes

, Romberg

and Tao

, Stable signal recovery from incomplete and inaccurate measurements, Commun Pure Appl Math 59 (2006), 1207–1223.

Sidky

E.Y.

, Kao

K.M.

and Pan

, Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT, J X-ray Sci Technol 14 (2006), 119–139.

Sidky

E.Y.

and Pan

, Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization, Phys Med Biol 53 (2008), 4777–4807.

Chen

G.H.

, Tang

and Leng

, Prior Image Constrained Compressed Sensing (PICCS), Med Phys 35 (2008), 660.

10.

Sidky

E.Y.

, Pan

, Reiser

I.S.

, Nishikawa

R.M.

, Moore

R.H.

and Kopans

D.B.

, Enhanced imaging of microcalcifications in digital breast tomosynthesis through improved image-reconstruction algorithms, Med Phys 36 (2009), 4920–4932.

11.

, Sun

and Yuan

, Optimization for limited angle tomography in medical image processing, Pattern Recogn 44 (2011), 2427–2435.

12.

Cho

and Lee

, Feasibility study on many-view under-samplingtechnique for low-dose computed tomography, Opt Eng 51 (2012), 527–529.

13.

, Liu

, Zhang

, Peng

, Wu

, Zheng

and Liang

, Image reconstruction from few-view CT data by gradient-domain dictionary learning, J Xray Sci Technol 24 (2016), 627–638.

14.

Sidky

E.Y.

, Jørgensen

H.J.

and Pan

, Convex optimization problem prototyping for image reconstruction incomputed tomography with the Chambolle-Pock algorithm, Phys Med Biol 57 (2012), 3065–3091.

15.

Goldstein

and Osher

, The split Bregman method for l1regularized problems, SIAM J Imag Sci 2 (2009), 323–343.

16.

Chartrand

, exact reconstruction of sparse signals via nonconvex minimization, IEEE Signal Process Lett 14 (2007), 707–710.

17.

Wang

Y.L.

, Yang

J.F.

, Yin

W.T.

and Zhang

, A new alternating minimization algorithm for total variation image reconstruction, SIAM J Imaging Sci 1 (2008), 248–272.

18.

Natarajan

, Sparse approximate solutions to linear systems, SIAM Journal on Computing 24 (1995), 227–234.

19.

Chartrand

, Shrinkage mappings and their induced penalty functions, IEEE International Conference on Acoustics, Speech and Signal Processing (2014), 1026–1029.

20.

Woodworth

and Chartran

, Compressed sensing recovery via nonconvex shrinkage penalties, CAM Rep 2014, 14–78.

21.

Chartrand

, Fast algorithms for nonconvex compressive sensing: MRI reconstruction from very few data, Proceedings IEEE International Symposium on Biomedical Imaging: From nano to macro 70 (2009), 262–265.

22.

Sidky

E.Y.

, Chartrand

, Boone

J.M.

and Pan

, Constrained TpV minimization for enhanced exploitation ofgradientsparsity: Application to CTimage reconstruction, IEEE J-Transl Eng Health Med 2 (2014), 1–18.

23.

Cai

, Wang

, Yan

, Li

, Zhang

and Hu

, Efficient TpV minimization for circular, cone-beam computed tomography reconstruction via non-convex optimization, Computerized Medical Imaging & Graphics 45 (2015), 1–10.

24.

Lou

, Yin

, He

and Xin

, Computing sparse representation in a highly coherent dictionary based on difference of L1 and L2, J Sci Comput 64 (2015), 178–196.

25.

Yin

, Lou

, He

and Xin

, Minimization of L1 - L2 for compressed sensing, SIAM J Sci Comput 37 (2015), A536–A563.

26.

Lou

, Zeng

, Osher

and Xin

, A weighted difference of anisotropic and isotropic total variation model for image processing, SIAM J Imaging Sci 8 (2015), 1798–1823.

27.

Mallat

and Zhang

, Matching pursuits with time-frequency dictionaries, IEEE T Signal Proces 41 (1993), 3397–3415.

28.

Tropp

and Gilbert

, Signal recovery from partial information via orthogonal matching pursuit, IEEE Trans Inf Theory 53 (2007), 4655–4666.

29.

Needell

and Vershynin

, Signal recovery from incomplete and inaccurate measurements via regularized orthogonal matching pursuit, IEEE J Sel Top Signal Proces 4 (2010), 310–316.

30.

Hiriart-Urruty

J.-B.

and Lemarechal

, Convex Analysis and Minimization Algorithms I. Springer-Verlag, Berlin Heidelberg, 1993.

31.

Zhang

, Ye

, Chen

, Perkins

A.E.

, Rose

, Sidky

E.Y.

, Kao

, Xia

, Tung

and Pan

, Investigation of optimization-based reconstruction with an image-total-variation constraint in PET, Phys Med Biol 61 (2016), 6055–6084.