l1/2 regularization for wavelet frames based few-view CT reconstruction

Reducing the radiation exposure in computed tomography (CT) is always a significant research topic in radiology. Image reconstruction from few-view projection is a reasonable and effective way to decrease the number of rays to lower the radiation exposure. But how to maintain high image reconstruction quality while reducing radiation exposure is a major challenge. To solve this problem, several researchers are absorbed in l0 or l1 regularization based optimization models to deal with it. However, the solution of l1 regularization based optimization model is not sparser than that of l1/2 or l0 regularization, and solving the l0 regularization is more difficult than solving the l1/2 regularization. In this paper, we develop l1/2 regularization for wavelet frames based image reconstruction model to research the few-view problem. First, the existence of the solution of the corresponding model is demonstrated. Second, an alternate direction method (ADM) is utilized to separate the original problem into two subproblems, where the former subproblem about the image is solved using the idea of the proximal mapping, the simultaneous algebraic reconstruction technique (SART) and the projection and contraction (PC) algorithm, and the later subproblem about the wavelet coefficients is solved using the half thresholding (HT) algorithm. Furthermore, the convergence analysis of our method is given by the simulated implementions. Simulated and real experiments confirm the effectiveness of our method.


Introduction
Nowadays, X-ray computed tomography (CT) has been extensively applied in medical diagnosis, biomedical research, non-destructive testing and so forth [1]. However, the inherent radiation dose of CT may induce cancer and other diseases in medical diagnosis [2,3]. The reduction of X-ray radiation is a more and more urgent issue. Generally, there are two strategies to reduce radiation dose: one is to lower the tube current or voltage values in data acquisition protocols, and the other is to decrease the number of the X-ray attenuation measurements through an object to be reconstructed. The former situation often introduces excessive noise into the projection data [4,5]. The latter situation necessarily results in insufficient projection data, which leads to few-view reconstruction [8][9][10][11], limited-angle reconstruction [8,13,14], etc.. In this paper, we focus on the few-view reconstruction within short-scan. Iconically, The sketch map of the scanning geometry for the few-view reconstruction problem within short-scan can be presented as Figure 1, where the red point denotes the sampled X-ray source. Restricted by the noise and the number of the projection data, the image reconstruction problem is usually an ill-posed inverse problem from a mathematical point of view [1] formulated as: . b e Ax   (1) where x denotes the discrete attenuation coefficients of the object to be reconstructed (that is the image); A is the system matrix; e represents the noise with level (that is ) and b is the measurement of the calibrated and log-transformed projection data. Generally, the solution of (1) can be found by minimizing the following optimization problem: . min When the projection data is complete, the commercial Filtered Back-Projection (FBP) [6] method is a common choice for two-dimensional reconstruction. However, if the projection data is incomplete or contains high noise, the reconstructed images using commercial FBP method will suffer from artifacts and noise. Fortunately, the iterative algorithms (such as the simultaneous algebraic reconstruction technique (SART) [6,7]) may perform better than FBP method. SART can suppress the noise well when the projection data with the noise is complete. However, this technique will lead to obvious artifacts in reconstructed images when projection data is incomplete. Several attempts have been carried out to settle such problem and overcome its shortcomings [4,5,8,11,15], which are often considered as optimization algorithms.
However, the problem (2) is an ill-posed problem and its solution will be instable when the projection data is incomplete. In order to surmount the instability of an ill-posed problem, numerous regularization methods have been researched. The thought of compressed sensing (CS) has been widely applied in CT image processing [16][17][18][19]. Typically, motivated by the advantage of the total variation (TV) in image denoising [20], TV based image reconstruction model as shown in (3) was proposed to solve few-view and limited-angle problems [8].
x , is the intensity value at the position (i, j).
Subsequently, adaptive steepest descent projection onto convex sets (ASD-POCS) was proposed to obtain relatively good image with suppressing streak artifacts [9], but it did not work for slope artifacts well. Derivative algorithms of TV [12][13][14] were put forward to further improve the slope artifacts for limited-angle reconstruction, however, the edge of an object was more or less contortive. Recently, regularization based image reconstruction optimization models are springing up [19,21,22], for example, based TV image reconstruction model was developed [19] as follows: There were several methods for solving lp based optimization problems such as the reweighted l1 norm method [23] and the thresholding method [24][25][26][27]. Although these TV based algorithms are successful in a large number of situations, the power of them are still limited for preserving gradually varied edges.
In order to preserve the detailed information of an image, some other forms of sparsifying transforms have been developed, such as wavelet frames [28][29][30][31], Haar transform [32], S-transform [33], etc.. The core of wavelet frames is on account of the sparsity of some features in an object to be reconstructed. Sparse transform takes some prior knowledge of the object to be reconstructed into account. To measure the transformed result is the second important aspect. As it is known that l0 based prior can obtain a sparser representation than l1 based prior [20]. However, l0 based image reconstruction is often an NP-hard (nondeterministic polynomial-time hard) problem and the objective function will be non-convex and noncontinuous [34]. In the reference [24], the authors proposed an l1/2 regularizer which has many promising properties such as unbiasedness, sparsity and oracle properties and it can be considered as a representative of lp (0 < p < 1) regularizers. Nevertheless, l1/2 based optimization problem is a non-convex problem, and most algorithms for solving that can only provide an approximate local minimizer [35]. Some authors investigated the existence of non-smooth and nonconvex optimization problems [36], and the other authors introduced a proximal alternating linearized minimization (PALM) algorithm to solve such problems [37], where the sequence generated by PALM can be convergent to a stationary point under some conditions.
In this study, we mainly concentrate on l1/2 regularization for wavelet frames based optimization problem to settle the few-view reconstruction within short-scan. To settle this problem, we utilize an alternate direction method (ADM) [37,39,42] to separate the original problem into two subproblems. In the first problem, considering the storage of the system matrix A, we utilize the SART [6] to obtain a proximal point and then use the projection and contraction (PC) [43,44] algorithm to settle the proximal problem. Ultimately, we use the half thresholding (HT) [27] algorithm to settle the second subproblem.
The main contents in this paper are shown below: 1. We use the l1/2 quasi-norm of the wavelet coefficients for the image to develop the image optimization model. After that we give the existence of a solution for the presented model.

2.
We utilize the ADM to separate the l1/2 quasinorm regularization term and l2 norm fidelity term into two subproblems. Then through the idea of the proximal mapping, we utilize SART and PC algorithm to solve the first subproblem, and then use HT algorithm to solve the second subproblem. In addition, we analyze the convergence of special situation of our method for the presented model. 3. To assess the effectiveness of our method, simulated and real experiments are implemented.
This paper is organized as follows. In section 2, we present the corresponding definitions and the fundamental model and then give the existence of the model's solution, and finally we give solving process of our method. In section 3, we represent the convergence analysis of our method, and the simulated and real experiments are given to demonstrate the effectiveness and accuracy of our method. Finally, the conclusions and perspectives are presented in section 4.

Existence of a Solution
Whether an optimization model has a global solution or not is a key issue. In this section, we first give some definitions and then develop an l1/2 regularization for wavelet frames based image reconstruction model for the few-view CT reconstruction problem within shortscan. In addition, we give the existence of a solution for the presented model.

Definitions
Some definitions [40] are prepared for the existence of the model's solution in this paper.
The function , and lower semi-continuous on R m if this holds for every x∈R m .
The level set of G from R m → R can be defined as: The asymptotic function G∞ for any proper function The kernel of  G can be given as: be an lsc and proper function. Then G is said asymptotically level stable (als) if for each ρ > 0, each bounded sequence of reals {λn}, and each sequence {xn}∈R n satisfying , ker , ), , ( there exists n0 such that xn-ρx * ∈ lev(G, λn), ∀n ≥ n0.

Model
Before presenting our model in this paper, we consider the following general problem: ).
( min x G X x (8) where G denotes a real and extended functional on a real space X. It is unnecessary that G is convex or smooth. The problem (8) has several contributions in preferences [36,38]. The authors considered that the existence of a solution to problem (8) can be obtained by utilizing a sequence of the problems in which β > 0 and it tends to 0: The purpose of the term 2 2 2 x  is to guarantee the existence of a solution to the problem (8) whether G is coercive or not.
In this paper, we mainly concentrate on the image optimization problem. Especially, we put emphasis on the l1/2 quasi-norm of the wavelet frame coefficients since the piecewise smooth function (such as an image) can be sparsely represented using tight wavelet frame transform. That is and X = Ω = {x|x = (x1, x2, . . . , xm) T ∈ R m , xi ≥ 0, i = 1, 2, . . . , m} which can be counted as m R  . Then, the presented model can be expressed as where x is the image vector comprised of pixel coefficients; is a projection vector; D denotes the positive definite diagonal matrix; W : R m → R q is a multilevel wavelet tight framelets transform operator; λ is a parameter balancing the fidelity term and regularization term; The piece constant linear B-spline framelets transform W which can be reconstructed by reference [31] is utilized in this paper.

Existence
A family of problems Gβ (β → 0 + ) to approximate the problem (10) are adopted to accomplish the existence of a solution. First, we need a corollary to accomplish the existence of a solution for the presented model. To prove the existence of (10), we should validate the conditions of Lemma 2.1 and consider the problem as follows: For simplicity, we let , 2

Theorem 2.2. (The existence of the presented model):
Let G and Gβ be defined as (12) and (13), respectively.
According to lemma 2.1, the problem (10) admits a global solution.

Solving Process
According to remark 1, we can solve the problem (11) to obtain the solution of the problem (10). In this section, numerical implementation of the problem (11) is presented. The alternating direction method (ADM) is utilized to deal with this problem. Considering the storage of the system matrix A, we will incorporate the SART with the ADM, and then use projection contraction (PC) algorithm to solve the constraint problem. For the l1/2 quasi-norm, we will use the half thresholding (HT) algorithm. we let (xi,j) denote a discretized image on a rectangular grid m1×m2, i = 1, 2, . . . , m1, j = 1, 2, . . . , m2, and (bi,j) denote a discretized image space of the function A on a rectangular grid l1 × l2, i = 1, 2, . . . , l1, j = 1, 2, . . . , l2. Where l1 is the number of the detector bins and l2 is the number of the projection angle. The projection function A can be defined by the basis functions which are given by the seventh formula in reference [7].
The key to dealing with the presented model is to de-couple the D-weighted l2 and l1/2 portions of the model. It is difficult to deal with the problem directly. Fortunately, ADM [37] is widely adopted to separate the original problem (11) into two sub-problems. The problem (11) can be converted into: Using the augmented lagrangian method, it implies that }.
The ADM minimizes the problem (15) by iteratively minimizing x and α alternately. Its (n + 1)th process can be expressed as follows: Although the problem (16) has a closed form solution, the system matrix A is difficult to store in memory considering the practical CT projection. To solve the step 1, in this paper, we utilize SART to get a proximal point by minimizing the objective function where τ is a balancing parameter. Note that there are three parameters to be tuned in (22). In fact, (22) expresses the proximal forward-backward scheme [37]. Using the proximal mapping, it follows that )}. , , ( Afterwards, the PC algorithm [43,44] is utilized to solve (24). The solution process of (16) includes two steps as follows: 1. Obtain a proximal point n x using the SART; 2. Optimize (24) using PC algorithm to obtain That can be expressed as The implementation process of PC algorithm [43,44] is as following: Step 3. To solve the problem (17) The solution of the problem (25) has a closed form and it can be effectively settled using the half thresholding (HT) algorithm [27] as follows.
The step 3 is to update the lagrangian multiplier. The above process including three steps can be regarded as SART-PC-HT method (our method). We denote Nmax the maximum number of the iteration, and is the relative error between 1  n x and n x . The pseudo-code of SART-PC-HT method can be described as shown in Table 2.

Simulated and Real Experiments
In this section, simulated and real experiments are implemented to verify the effectiveness of our method. We utilize the simulated experiments with ten different levels of Guassian noise to do the quantitative statistics interpretations which include RMSE and PSNR [41] as follows: (27) . log 10 2 2 10 RMSE x PSNR r   (28) where x denotes the reconstructed image, xr is the reference image and m is the total number of the image pixels. Then, we use the real experiments to conform the effectiveness of our method, compared with the commercial FBP method [6]. In the experiments, the parameters of our method are chose by trial and error for the best reconstructed image quality. All the experiments are implemented on 2.90 GHz intel(R) Pentium(R) G2020 CPU processor with 4G memory and coded.

Simulated and Real Experiments
We test and verify the effectiveness of our method for few-view reconstruction within short-scan using a digital NURBS based cardiac torso (NCAT) phantom [45]. The simulated projection data are generated by projecting a 256×256 discretized NCAT phantom with adding ten different levels of Gaussian noise, where all the mean value of Gaussian noise are zero, and their standard deviations are 0.1%, 0.2%,..., 1% of the current projection values (cpv), respectively. The geometry scanning parameters for CT imaging system are listed in Table 3. The number of projection views, the parameters of our method and the corresponding statistical analysis RMSE and PSNR are listed in Table  2. The stopping criterion is reaching the maximum iteration number 100 max  N for ten different levels of Gaussian noise. In PC algorithm, maximum iteration number is 2000 and stopping criterion is As shown in Table 4, with the decrease of the standard deviation and the increase of the projection views, RMSE is decreasing and PSNR is increasing. Meanwhile, the main parameters of our method are chose differently with diverse standard deviation. Figure 2 demonstrates some of the reconstructed results from ten different levels of Gaussian noise. As shown in Figure1, (a) is an original image for NCAT, and the rest of them are reconstructed using our method from different projection views by adding Gaussian noise with standard deviation 0.3%cpv, 0.5%cpv, 0.8%cpv and 1%cpv, respectively. From the enlarged images in Figure 2, with the decrease of the standard deviation and the increase of the projection views, the noises are restrained effectively and details are maintained well.
With the analysis of Figure 3 and Figure 4, the RMSE of the reconstructed NCAT phantoms is gradually decreasing and their PSNR are increasing with the standard deviation decreasing and the projection views increasing. That is, the smaller standard deviation and more complete projection views will lead the sequence {xn} to converge to a stationary point which is very close to a global point. Meanwhile, the bigger standard deviation and fewer projection views will make the sequence {xn} converge to a stationary point early, but the reconstructed result is not good for the details. Table 3. Geometry scanning parameters of NCAT for simulated CT imaging system The distance from X-ray source to rotation center 900.0mm The distance from detector to rotation center 400.0mm The angle between two adjacent projection views(interval angle) 0.6667 0 The angle between two adjacent rays 0.0011 0 The number of detector units 256

Reconstruction from Real Walnut Data
To further confirm the effectiveness and accuracy of our method for few-view reconstruction within shortscan, we use the real Walnut projection data which can be well known in reference [46]. Some details of the Walnut are difficult to reconstruct since its projection data contains some noises. In order to highlight our method's effectiveness to few-view reconstruction, we adopt the different projection views in [0,180 o ], compared with commercial FBP method. The size of Walnut image to be reconstructed is 656×656. The real CT imaging system can refer to [46]. In the experiments, the stopping criterion is reaching Nmax = 100, and the main parameters of our method in this subsection are λ = 0.00003 and γ = 0.003. In SART, the parameter ω = 1. In PC algorithm, τ = 1, maximum iteration number is 2000 and stopping criterion is ε = 1 × 10 -5 , and the other stopping criterion is ε0 = 1 × 10 -6 to our method. Figure5 and Figure 6 demonstrates the reconstructed Walnut images using commercial FBP method and our method. As shown in Figure 5, the first row are the images using commercial FBP method from 600, 300 and 200 projection views in angular scope [0, 180 o ], respectively; the third row are the results using our method from 600, 300 and 200 projection views in angular scope [0, 180 o ], respectively; the second row and third row are the enlarged images of the red rectangle box of the first row and third row, respectively; the lower left corner images of the first and third row are the enlarged images of the red rectangle box. As shown in Figure 6, the reconstructed Walnut images using commercial FBP method and our method are reconstructed from 150, 100 and 76 projection views in angular [0, 180 o ], respectively. The layout of the rest in Figure 6 is as shown in Figure 5.  According to Figure 5 and Figure 6, it is observed that the Walnut images using FBP method contain the larger noise and some details of them are distorted along with the decrease of the projection views. Although the number of the projection views are decreasing, the results using our method have relatively high quality and contrast, except the rupture of the filament at the top of the Walnut image. It can be seen our method outperforms the commercial FBP method in suppressing noise and preserving details, and our method has the ability for reducing the projection views to lower the radiation exposure while maintaining relatively high image quality.

Conclusions and Perspectives
To deal with the few-view reconstruction problem within shortscan, we adopt l1/2 regularization for wavelet frames based image reconstruction model, and adopt the SART-PC-HT method to solve it. First, the existence of the presented model is analyzed through the existing corollary. Second, ADM is utilized to separate the original problem into two subproblems, which can be settled using SART, PC algorithm and HT algorithm, respectively. Third, simulated NCAT experiment and real Walnut experiment confirm that our method can obtain the reconstructed images with relatively good quality and high contrast while decreasing the projection views. That means our method has the ability to settle the few-view reconstruction problem within short-scan whose projection data contains noise. However, when the projection data are very incomplete and has high-level noise, the details of the reconstructed images using our method will not be preserved very well.
In our method, the parameters are chose by trial and error, which are not adaptive for any given projection data. In the future, we will analyze the adaption of the parameters in our method. For high level noise and few projection views, our method will be not effective, which will be improved with taking the type of the noise of the projection data into consideration.