# An Algorithm ofl1-Norm and l0-Norm Regularization Algorithm for CT Image Reconstruction from Limited Projection

Authors: Xiezhang Li; Guocan Feng; Jiehua Zhu
Journal: International Journal of Biomedical Imaging
 (2020)

---

## Abstract

Thel1-norm regularization has attracted attention for image reconstruction in computed tomography. The l0-norm of the gradients of an image provides a measure of the sparsity of gradients of the image. In this paper, we present a new combined l1-norm and l0-norm regularization model for image reconstruction from limited projection data in computed tomography. We also propose an algorithm in the algebraic framework to solve the optimization effectively using the nonmonotone alternating direction algorithm with hard thresholding method. Numerical experiments indicate that this new algorithm makes much improvement by involving l0-norm regularization.


---

## Body

## 1. Introduction

Computed tomography (CT) has been widely applied in medical imaging and industry for over decades. In some applications, the reconstructions do not require a high-quality image, or complete projection data for exact reconstructions is not available because of the scanning conditions or radiation issue. CT reconstruction from limited projection data is of particular importance, and statistical and iterative reconstruction algorithms outperform the analytic approaches in these applications.The theory of compressed sensing [1, 2] has provided a novel framework for CT image reconstruction from limited projection data with the prior knowledge of the sparse gradients for CT images. The total variation (TV) based regularization [3] for CT image reconstruction is connected to compressed sensing under the sparse transform of the gradient operator [4, 5]. The TV of an image is the l1-norm of the gradients of the image, so the TV minimization is also known as the l1-minimization.Thel1-norm-based regularization can greatly reduce the streak artifacts arising from few-view CT [6]. However, the results are over smoothing and lose some detailed features including contrast that causes edge blurred [7]. In addition, the l1-norm regularization does not provide a representation sparse enough. To overcome the drawback, many improved TV-based reconstruction algorithms have been developed. One of them is the reweighted l1 minimization [8]. Algorithms for regularization involving multiple norms are also developed. For example, the alternating direction method (ADM) [9, 10] and the nonmonotone alternating direction algorithm (NADA) [11, 12] are proposed for solving an l1-norm minimization with an l2-norm constraint.Thel0-norm of an image vector, defined as the number of its nonzero elements, measures the sparsity of a vector appropriately. So the l0-norm-based regularization for CT represents the sparsity of the gradient image better than the l1-norm-based regularization and preserves the edge while suppressing the streak artifacts [13]. However, the l0-norm is a nonconvex function and the l0-norm regularization problem is NP hard [14]. So many variants of the regularization algorithms are developed. An edge-preserving image reconstruction method for limited-angle CT is investigated based on l0-norm regularized gradient prior [15]. An imaging reconstruction model for few-view CT based on l0 sparse regularization is proposed [16]. A new l0 regularization with wavelet tight framelets is addressed to suppress the slope artifacts in the limited-angle X-ray CT reconstruction [17]. Regularization involving multiple norms is also developed to address the individual drawbacks of the l1-norm and l0-norm. An image reconstruction model based on l0-norm and l2-norm regularization for the limited-angle CT is proposed [18]. Numerical experiments indicate that the algorithm has the advantage in suppressing slope artifacts. A combined smoothed l0-norm and l1-norm regularization algorithm using the NADA method for CT image reconstruction is proposed and demonstrated to have better performance than the l1-norm regularization [19].In this paper, we will develop a new combinedl1-norm and l0-norm regularization model for CT image reconstruction from limited projection data. The main contributions of this paper are as follows:
(1)
We generalize thel1-norm regularization to an l1-norm and l0-norm regularization for CT image reconstruction from limited projection data. The proposed model can achieve superior performance compared with the existing related methods for multiple reasons. Firstly, due to the l1-norm regularization term, the proposed model can reduce streak artifacts in limited-data CT. Secondly, due to the l0-norm regularization term, the oversmoothing from the l1-norm regularization term is improved. The proposed model is less smoothing, thus better edge-preserving, and provides a better sparsity representation. Thirdly, since the l0-norm is adopted without using the smoothed l0-norm approximation, the proposed model has the advantage in suppressing slope artifacts(2)
We propose an algorithm in the algebraic framework to solve the optimization effectively using the nonmonotone alternating direction algorithm with hard thresholding (NADA-HT) methodThe rest of the paper is organized as follows. In Section2, we present a combined l1-norm and l0-norm-based regularization model and propose an algorithm for CT image reconstruction from limited projection data. Numerical experiments for a comparison of regularization models with/without the l0-norm and the discussion are given in Section 3.


## 2. Materials and Methods

The projection data in CT can be modelled as a linear system,(1)Φf+e=u,where Φ is an m×n2 projection matrix, f∈Rn2 represents a 2D n×n image to be reconstructed, e an additive noise with e2≤ε for some known ε>0, and u∈Rm the noisy projection data. For limited data reconstruction, the underdetermined system (m≪n2) has infinitely many solutions. An optimal solution representing the original image as good as possible can be sought by minimizing the energy function with a TV regularization term [8]:
(2)minffTVsubjecttoΦf+e=u,where the total variation fTV is the l1-norm of the magnitude of the discrete gradients,
(3)fTV=∇f1=∑i,jfi+1,j−fi,j2+fi,j+1−fi,j2.Then, a 2D imagef with sparse gradients and noisy measurements in CT can be reconstructed by solving the following l1-norm minimization with an l2-norm constraint,
(4)minf∇f1subjecttoΦf−u22≤ε.However, the image reconstructed by solving (4) has slope artifacts near the edge of the object for limited-angle projection data [20].We consider the following minimization model(5)minf∇f1+α∇f0+γ2f22subjecttoΦf−u22≤ε,where ∇f0 stands for the number of nonzero gradients ∇fp, 1≤p≤n2. The term γ/2f22, γ⟶0+, makes the energy of reconstructed image to be minimized and reduces the ill-posedness of CT reconstruction [18].For convenience, we denoteDpf as the forward difference of f at a pixel p in both horizontal and vertical directions, i.e.,
(6)Dpf=Di,jf=fi+1,j−fi,j,fi,j+1−fi,jT∈R2,forp=i−1n+j.Then, ∇fp=Dpf for 1≤p≤n2 and ∇f1=∑p=1n2Dpf2.Minimization (5) is equivalent to
(7)minfβ2Φf−u22+∑p=1n2Dpf2+αDpf0+γ2f22,where the notation Dpf0 is set to be 1 for nonzero Dpf, otherwise 0. The above minimization can be converted, with a parameter β>0 and vectors vp′s, to
(8)minfβ2Φf−u22+∑p=1n2vp2+αvp0+γ2f22subjecttoDpf=vp,1≤p≤n2.Adapting Lagrangian vectorsλp∈R2, 1≤p≤n2, we convert Minimization (8) into the following problem using the augmented Lagrangian method
(9)minfβ2Φf−u22+∑p=1n2vp2+αvp0+μ2Dpf−vp+λp22+γ2f22.For simplicity, we writeDf, v, and λ as vectors whose p-th components are Dpf, vp, and λp, 1≤p≤n2, respectively. Then, Df−v+λ22=∑p=1n2Dpf−vp+λp22. Thus, Minimization (9) becomes
(10)minfβ2Φf−u22+∑p=1n2vp2+αvp0+μ2Df−v+λ22+γ2f22.Next, we will develop an algorithm for solving Minimization (10). By convention, Minimization (10) can be solved by ADM iteratively. The k-th step of the ADM involves three procedures,
(11)fk+1∈argminfβ2Φf−u22+μ2Df−vk+λk22+γ2f22,(12)vk+1∈argminv∑p=1n2vp2+αvp0+μ2Dfk+1−v+λk22,(13)λk+1=λk−vk+1−Dfk+1.First, Minimization (11) can be solved by ADM with the gradient decent method involving nonmonotone line search, for example, the NADA.Next, we seek for a solution of (12). Minimization (12) is equivalent to finding
(14)vpk+1∈argminvpvp2+αvp0+μ2Dpfk+1−vp+λpk22,for1≤p≤n2.We introduce a hard thresholds (HT) operatorHκw on R2 with the threshold κ defined as
(15)Hκw=0ifw2≤κwotherwise.A solutionvpk+1 is provided in the following theorem.Theorem 1.
Letw=Dpfk+1+λpk and κ=1/μ1+2μα. If 2μα>1, then a minimizer vpk+1 for Minimization (14) is given by
(16)vpk+1=max1−1μw2,0Hκw.Proof.
For simplicity, denotez=vp. The objective function in Minimization (14) is rewritten as
(17)gz=z2+αz0+μ2z−w22.Without loss of generality, we assume thatw≠0 since zopt=0 in the case of w=0.Forz=0, gz=μ/2w22.Forz≠0, gz=z2+α+μ/2z−w22. The minimizer is derived in the following two cases.
(i)
Case 1.μw2>1We claim thatgmin=w2+α−1/2μ at z=1−1/μw2w. In fact, it follows from
(18)dgdz=1z2z+μz−w=0that a minimizer z is a scalar multiple of w. Let z=tw where t>0 due to the minimization problem. Then, the function
(19)ht≔gz=μ2tw−w22+tw2+α=μ2w22t2−μw22−w2t+μ2w22+αreaches its minimum value
(20)hmin=α+μ2w22−μw22−w222μw22=α+w2−12μat t=1−1/μw2. The claim is proved.
(ii)
Case 2.μw2≤1In this case,h′t=μw22t−1+1/μw2>0 for all t>0 and limt→0+ht=μ/2w22+α>μ/2w22. Thus, the minimizer is zopt=0.Combining the results forz=0 and z≠0, we conclude that if μ/2w22≤α+w2−1/2μ then, zopt=0; otherwise, zopt=max1−1/μw2,0w.Next, we look for conditions under which the inequalityμ/2w22≤α+w2−1/2μ holds. Substituting r=w2 into the inequality, we consider the equation
(21)qr≔μ2r2−r−α−12μ=0in an unknown r. The zeros of the equation are given by r1,2=1/μ1∓2μα. So qr≤0 for r∈r1,r2, where r1≤0 for 2μα>1, and qr>0 for r>r2. Therefore, for 2μα>1, if w2≤1/μ1+2μα, we have μ/2w22≤α+w2−1/2μ, and consequently, zopt=0; otherwise, μ/2w22>α+w2−1/2μ, and consequently, zopt=max1−1/μw2,0w.With a thresholdκ=1/μ1+2μα, if w2≤κ, then zopt=0; otherwise, zopt=max1−1/μw2,0w. Using the HT operator in (14), we have zopt=max1−1/μw2,0Hκw. We complete the proof of the theorem.Remark. In practice, the condition2μα>1 for nonzero α in Theorem 1 is fulfilled. Now we summarize our algorithm as follows.Algorithm 1.
(l1-norm and l0-norm regularization algorithm using NADA-HT method).1. inputΦ,u,α,ε2. initializen,β,γ,μ,k,f0,v0,λ0,maxit,ratio,tol3. whilek<maxit3.1   updatefk+1 by NADA3.2.  for each integerp, 1≤p≤n2, update vpk+1 by HT methodsetw=Dpfk+1+λpk, κ=1/μ1+2μαupdatevpk+1=max1−1/μw2,0Hκw3.3  updateλk+1=λk+Dfk+1−vk+13.4  if error<tol then output fk+1; stop3.5γ=γ×ratio3.6k=k+1end.


## 3. Results and Discussion

In this section, the performance of the proposed algorithm is compared with that of thel1-norm regularization algorithm. Both algorithms are implemented in MATLAB using the NADA and tested with the 2D Shepp-Logan phantom and a cardiac image [21] of size 128×128 on an Intel Core i7 3.40 GHz PC. In each test, the same system u=Φf+e is created, where Φ∈Rm×n2 (m≈0.3n2) is a random matrix and the noise e=0.02∗meanΦf∗randnm. The weight parameter α is chosen as 0 and 1 for the l0-norm regularization and the proposed algorithm, respectively. Other two important parameters in the objective function are taken as β=28, μ=24 for the Shepp-Logan phantom and β=27, μ=22 for the cardiac image, respectively.Experiments are conducted to compare the reconstruction by the two algorithms after the same number of iterations. The original and reconstructed images for two images after the same numbers of iterations are shown in Figure1. The images demonstrate that the proposed algorithm yields much better reconstruction. The quality of reconstructed images is compared using several commonly used criteria: the relative error, the root-mean-square error (RMSE), the normalized root mean square deviation (NRMSD), the normalized mean absolute deviation (NMAD), and the structural similarity index (SSIM). The experimental results from 100 tests are summarized in Table 1. It shows that the proposed algorithm produces significantly improved evaluation measurement: 88%–167% for the Shepp-Logan phantom and 67%–114% for the cardiac image, respectively, while taking less CPU time. Overall, the new algorithm provides better accuracy after the same iteration number.Figure 1
Original and reconstructed images. First row: Shepp-Logan phantom after 200 iterations. Second row: cardiac image after 500 iterations.Table 1
Comparison of two algorithms after same number of iterations.
(a)Shepp-Logan phantom, after 200 iterationsTime(s)ErrorRMSENRMSDNMADSSIMl1 regularization530.3410.0830.0860.4220.356l0+l1 regularization520.0350.0090.0090.0330.949(b)Cardiac image, after 500 iterationsTime(s)ErrorRMSENRMSDNMADSSIMl1 regularization880.3620.1110.0760.3520.442l0+l1 regularization830.1190.0140.0100.0410.936Figure2 shows the graph of relative error v.s. the number of iterations for the two algorithms. The curves indicate that the proposed algorithm requires much less iterations and time to achieve the same accuracy and yields much faster convergence. After a certain number of iteration numbers, the relative error of the proposed algorithm drops sharply, while the relative error of the l1-norm regularization is slowly decreasing and does not improve much over the iterations.Figure 2
Relative error v.s. iteration number in reconstruction by two algorithms. (a) Shepp-Logan phantom. (b) Cardiac image.
(a)
Shepp-Logan phantom(b)
Cardiac imageThe proposed regularization model introduces anl0-norm term in addition to an l1-norm term. The l0-norm of a vector is the number of nonzero elements in the vector; thus, the l0-norm is appropriate for representing sparsity. Meanwhile, the l0-norm term makes the model less smoothing. Numerical experiments demonstrate that, in the proposed model, the edge is better preserved, and a better sparsity representation is provided. So, the proposed algorithm improves the quality and the efficiency of the reconstruction even with the extra computation from the added l0-norm term.The effect of the weight parameterα of the l0-norm is tested. From our experiments, the parameter α from 0.5 to 1 almost does not not affect the efficiency of the algorithm. Selecting good values of other parameters β, μ, and γ in Minimization (10) is a challenging problem and will be further investigated in the future.High-quality CT image reconstruction from limited projection data is challenging. Learning based algorithms such as structure-aware sparse Bayesian learning [22] could yield improved performance in reconstructing tomographic images from limited data, since structural prior knowledge is exploited. Enhancing the proposed algorithm in a learning-based framework using prior knowledge will be the topic of future research.



---

*Source: 8873865-2020-08-28.xml*