A Smoothed 𝑙 0 -Norm and 𝑙 1 -Norm Regularization Algorithm for Computed Tomography

. The nonmonotone alternating direction algorithm (NADA) was recently proposed for effectively solving a class of equality-constrained nonsmooth optimization problems and applied to the total variation minimization in image reconstruction, but the reconstructed images suffer from the artifacts. Though by the 𝑙 0 -norm regularization the edge can be effectively retained, the problem is NP hard. The smoothed 𝑙 0 -norm approximatesthe 𝑙 0 -norm as a limit of smooth convex functions and provides a smooth measure of sparsity in applications. The smoothed 𝑙 0 -norm regularization has been an attractive research topic in sparse image and signal recovery. In this paper, we present a combined smoothed 𝑙 0 -norm and 𝑙 1 -norm regularization algorithm using the NADA for image reconstruction in computed tomography. We resolve the computation challenge resulting from the smoothed 𝑙 0 -norm minimization. The numerical experiments demonstrate that the proposed algorithm improves the quality of the reconstructed images with the same cost of CPU time and reduces the computation time significantly while maintaining the same image quality compared with the 𝑙 1 -norm regularization in absence of the smoothed 𝑙 0 -norm.


Introduction
Statistical and iterative reconstruction algorithms in computed tomography (CT) are widely applied since they yield more accurate results than analytic approaches for low-dose and limited-view reconstruction.These algorithms involve solving a linear system: where Φ is an  ×  2 projection matrix,  ∈   2 represents a 2D  ×  image to be reconstructed,  is an additive noise with ‖‖ 2 ≤  for some known  > 0, and  ∈   is the noisy projection data.For limited-view reconstruction, the underdetermined system ( ≪  2 ) has infinitely many solutions.An optimal solution representing the original image as well as possible is sought by iteration methods.
The theory of compressed sensing [1][2][3] has recently shown that signals and images that have sparse representations in some orthonormal bases can be reconstructed at high quality from much less data than what the Nyquist sampling theory [4] requires.In many cases in tomography, images can be approximately modeled to be essentially piecewise constant so the gradients are sparse.With the gradient operator as a sparse transform, compressed sensing provides a novel framework for CT image reconstruction [1,3,5].The  0 -norm of a vector, the number of its nonzero elements, is an appropriate measurement for the sparsity of a vector.So the  0 -norm regularization might be an approach for the reconstruction of an image with sparse gradient.However, the  0 -norm is a nonconvex function and the  0 -norm regularization problem is NP hard [6].A pseudoinverse transform of the discrete gradient and the iterative hard thresholding algorithm are used to address the issues of the  0 -norm [7].The smoothed  0 -norm ( 0norm) provides a smooth measure of sparsity and is applied in compressed sensing MRI imaging [8].The  0 -norm is used to find the jointly sparse representation via the lowresolution image [9]. 0 -norm regularization model is proposed for sparse-view X-ray CT reconstruction [10].Two new smoothed functions approximating the  0 -norm are proposed in the mechanism of reweighted regularization [11,12].An edge-preserving image reconstruction method based on  0 -regularized gradient prior for limited-angle CT is investigated [13].The local and global minimizers of  0 gradient regularized model with box constraints for image restoration are discussed [14].A new  0 regularization and wavelet tight framelets to suppress the slope artifacts in the limited-angle X-ray CT reconstruction are addressed [15].An image model based on  0 and  2 regularizations for the limited-angle CT is proposed, and the existence of a solution and a local convergence analysis under certain conditions are proved [16].
The  1 -norm regularization, known as total variation (TV) minimization [17], has been widely used in CT reconstruction.Under certain conditions, the  1 -norm minimization has a unique sparsest solution which is a good approximation to the original image [18].However, the  1 -norm regularization provides less sparsity representation than the  0 -norm regularization and may lose some detailed features and contrast.In addition, the discrete gradient transformation ∇ does not satisfy the restricted isometry property required by the compressed sensing theory.
Researchers in this area have been seeking for other efficient and stable algorithms inspired by the compressed sensing theory, for example, solving an  1 -norm minimization with an  2 -norm constraint using the alternating direction method (ADM) [19,20] and the nonmonotone alternating direction method (NADA) [21,22].Regularization including multiple norms is developed to address the drawbacks of the  1 -norm and  0 -norm.A combined  1 -norm and  0 -norm regularization minimization with an  2 -norm constraint using SART algorithm and the gradient decent method is proposed for sparse-view CT image reconstruction in [23].It is observed in the paper that the convergence is slow and the computation is time consuming because of the alternative minimization of the  1 -norm and  0 -norm.
The purposes of this paper are multifold: (i) Presenting a combined  0 -norm and  1 -norm regularization algorithm for image reconstruction in CT: The proposed algorithm is unique in the following aspects: A new parameter is introduced to balance the  0 -norm and  1 -norm terms; two-norm regularization is combined in one Lagrangian object function to be minimized (ii) Adopting a newly developed alternating direction method NADA to efficiently solve minimization (iii) Resolving the computation challenge problem caused from the  0 -norm minimization The numerical experiments demonstrate that the proposed algorithm improves the quality of the recovered images for the same cost of CPU time and reduces the computation time significantly while maintaining the same image quality compared with the  1 -norm regularization in absence of the  0 -norm.
The rest of the paper is organized as follows.Section 2 includes the background and notations.In Section 3, a combined  0 -norm and  1 -norm regularization algorithm with the NADA is developed.Numerical experiments with the Shepp-Logan phantom and a cardiac image are presented in Section 4. Finally, Section 5 concludes the paper by discussions and conclusions.

Background and Notations
The TV minimization of an image  for solving system (1) is mathematically represented as [17] min where the total variation ‖‖  is the  1 -norm of the magnitude of the discrete gradients: Then a 2D image  with sparse gradients and noisy measurements in CT can be reconstructed by solving the following  1 -norm minimization with an  2 -norm constraint: which can be implemented [24] by Here  > 0 is a regularization parameter which controls the trade-off between sparsity of gradients and accuracy of approximation.
For convenience, we denote    as the forward difference of  at a pixel  in both horizontal and vertical directions; i.e., for  = ( − 1)  + .
Then ‖∇‖ 1 = ∑  2 =1 ‖  ‖ 2 and minimization ( 5) is rewritten as With the Lagrangian function [19,20]  (, V, ) minimization ( 7) is converted to which can be solved by the ADM [19][20][21][22].The ADM is a variant of the classic augmented Lagrangian method for optimization.It has been applied to solve different types of  1 -minimization problems for sparse solution recovery.The th step of the ADM for solving (9) involves the procedures Minimization (10) and minimization ( 11) can be successfully solved using the framework of the NADA with a nonmonotone line search scheme recently developed in [21,22].The idea is outlined as follows.
(ii) Select a step size   uniformly bounded above such that ( Combined with the solution of ( 11), the sequence {(  , V  ,   )} is bounded above by a monotonically nonincreasing sequence {  } though {(  , V  ,   )} itself is not decreasing.The advantage of the NADA lies in the fact that (, V, ) is not required to be differentiable while reducing the computation complexity and improving the efficiency.The convergence analysis of the NADA for solving a general model including minimization (11) can be found in [21].
The  0 -norm approximates the  0 -norm by a smooth function [8].The  0 -norm of a vector  ∈   2  , denoted by ‖‖ 0 , is defined as It is easy to see that ‖‖ 0 = ‖‖ 0 .A new approach for solving a more general minimization than the regularization minimization in [23] is proposed in the next section.

𝑠𝑙 0 -Norm and 𝑙 1 -Norm Regularization Algorithm
In this section we present  0 -norm and  1 -norm regularization algorithm for CT reconstruction using the NADA.Consider the following constrained minimization problem for image reconstruction in CT: where the parameter  is used to balance the two terms of the minimization.While the  1 -norm is used for the approximation accuracy, the  0 -norm is introduced in consideration of the sparsity of gradients.It is remarked that minimization ( 14) is an extension of minimization (4) ( = 0) and the minimization ( = 1) in [23].

New Approach.
The term ‖∇‖ 0 in ( 14) can be approximated by a smooth function for a small positive parameter .So minimization (14) Using a positive parameter  and Lagrange vectors   ∈  2 , 1 ≤  ≤  2 , we define a Lagrangian function Finally, minimization ( 18) is converted to the following minimization: It is obvious that   (, V, ) is an extension of (, V, ) in (8).

Calculation of V 𝑘+1
.In this subsection we address how to calculate V +1  in ( 22) for 1 ≤  ≤  2 .Actually, for each , we need to solve for V  from the following equation: It follows from (24) that If ‖V  ‖ 2 is calculated then It is obvious that there is a positive zero  of ℎ() between 0 and  and  →  as  → 0. Equation ( 28) can be rewritten as By substituting  = √, the above equation becomes Thus, we have shown the following lemma.
Before we develop a procedure to calculate a positive root of () = 0, we will list some properties of () below.
The second derivative indicates the following features of ().
Proof.We identify an interval containing a positive root of () = 0 on which   () and   () do not change signs for each subcase.Then the convergence of Newton's method with an initial guess above is guaranteed.
Case 1.For subcase (i), there is a positive solution of () = 0 on ( 2 2 ,  2 ) on which both   () and   () are negative.For subcase (ii), the inflection point  of the graph of () is on the upper half plane.So a positive zero is on ((4/9) 2 ,  2 ) on which both   () and   () are negative.So one should choose  0 <  2 and near  2 with ( 0 ) < 0 for the convergence of Newton's method.

Algorithm.
Based on the above discussion, the calculation process of the proposed  0 +  1 regularization algorithm is summarized in Algorithm 1.

Numerical Experiments and Results
In this section, the proposed combined  0 -norm and  1 -norm regularization is compared with the regularization without  0 -norm for its performance.Both regularization algorithms are implemented using the NADA.The MATLAB code is developed based on the software package TVAL3 [22], and the numerical experiments are conducted on an Intel Core i7 3.40 GHz PC.The 2D Shepp-Logan phantom and a cardiac image [25] of size 128 × 128 are tested.In each test, a random   2, indicates that the  0 +  1norm regularization requires less iterations to achieve the same accuracy.Both types of the numerical experiments demonstrate that the proposed  0 +  1 -norm regularization is superior to the regular  1 -norm regularization minimization.

Discussions and Conclusions
The application of the  0 -norm in CT reconstruction has become one of active research topics recently.A combined  0 -norm and  1 -norm regularization algorithm is proposed in this paper.The  0 -norm of a vector  is the limit of a smooth convex function   ().In iterations of the proposed algorithm, the parameter  is getting smaller and closer to zero.The smooth function   (∇) in the objective function   (, V, ) approaches ‖∇‖ 0 , and the sparsity of the recovered image is fulfilled quickly.So the proposed algorithm improves the quality of the reconstructed images and the efficiency in terms of the iteration number and CPU time even with the extra computation from the added  0 -norm term.The numerical results demonstrate that the proposed algorithm improves the  1 -norm regularization with the NADA in reducing the computation time significantly.The effects of the weight parameter  of the  0 -norm and the  0norm parameter  are also tested.From our experiments, the parameter  from 0.5 to 2 and the initial value  0 from 0.05 to 0.9 almost do not affect the efficiency of the algorithm.The decreasing ratio of parameter  between 0.8 and 0.95 is a good choice, but a ratio smaller than 0.7 is not recommended.
There are some aspects to be further investigated.It is a challenging problem to select good values of other parameters such as  and  in the Lagrangian function   (, V, ).The impact of these parameters on the performance of the algorithm will be further investigated.

Data Availability
The MATLAB numerical data used to support the findings of this study are available from the corresponding author upon request.

Figure 2 :
Figure 2: Relative error vs. iteration number in reconstruction by two algorithms to achieve the same tolerance of relative error 0.05.(a) Shepp-Logan phantom; (b) cardiac image.

Table 2 :
Comparison of two algorithms with the same tolerance of relative error 0.05.The average iteration numbers and CPU time from 100 tests of the  0 +  1 -norm regularization and  1 -norm only regularization are listed in Table2.The comparison indicates that the CPU time is significantly reduced in the proposed algorithm while producing the same image quality.The graph of relative errors vs. the number of iterations for the two algorithms, shown in Figure