Junfeng Zhang

Yining Hu

Jian Yang

Yang Chen

Jean-Louis Coatrieux

Limin Luo

Sparse-view X-ray CT Reconstruction with Gamma Regularization

Keywords: computer tomography, Gamma regularization, TV regularization, l 2 -norm

By providing fast scanning with low radiation doses, sparse-view (or sparse-projection) reconstruction has attracted much research attention in X-ray computerized tomography (CT) imaging. Recent contributions have demonstrated that the total variation (TV) constraint can lead to improved solution by regularizing the underdetermined ill-posed problem of sparse-view reconstruction. However, when the projection views are reduced below certain numbers, the performance of TV regularization tends to deteriorate with severe artifacts. In this paper, we explore the applicability of Gamma regularization for the sparse-view CT reconstruction. Experiments on simulated data and clinical data demonstrate that the Gamma regularization can lead to good performance in sparse-view reconstruction.

Introduction.

By non-invasively providing high-contrast human anatomic information, X-ray CT has become an important tool in radiological examination since its introduction in 1970s. However, accumulated radiation doses due to repeated X-ray CT examinations can induce cancer and other diseases in human bodies [START_REF] Brenner | Estimated risks of radiation-induced fatal cancer from pediatric CT[END_REF][START_REF] Brenner | Computed tomography-an increasing source of radiation exposure[END_REF][START_REF] De González | Risk of cancer from diagnostic X-rays: estimates for the UK and 14 other countries[END_REF]. The quality of CT images being dependent of the radiation dose [START_REF] Hsieh | Computed Tomography Principles, Design, Artifacts, and Recent Advances[END_REF], how to improve the image quality under low dose CT scanning has become a major research topic. Generally, there are two main approaches to reduce the X-ray radiation dose, which are lowering the tube current (or voltage) values and sparse-view scanning. Lowering tube current (or voltage) unavoidably leads to less effective photons and to increased noise in projection data (or sinogram) and consequently degrades the quality of reconstructed images. Much work has been conducted to obtain clinically acceptable CT images from low dose projection data, and the main solutions include sinogram data restoration [START_REF] Hsieh | Adaptive streak artifact reduction in computed tomography resulting from excessive X-ray photon noise[END_REF][START_REF] Kachelriess | Generalized multidimensional adaptive filtering for conventional and spiral single-slice multi-slice, and cone-beam CT[END_REF], statistical reconstruction [START_REF] Zeng | A filtered back projection MAP algorithm with non-uniform sampling and noise modeling[END_REF][START_REF] Wang | Penalized weighted least-squares approach to sinogram noise reduction and image reconstruction for low-dose x-ray CT[END_REF][START_REF] Chen | Nonlocal prior Bayesian tomographic reconstruction[END_REF], and image based post-processing [START_REF] Chen | Thoracic low-dose CT image processing using an artifact suppressed large-scale nonlocal means[END_REF][START_REF] Chen | Improving abdomen tumor low-dose CT images using a fast dictionary learning based processing[END_REF][START_REF] Chen | Artifact suppressed dictionary learning for low dose CT image processing[END_REF].

Reducing the number of projections is another solution. Sparse-view reconstruction leads to low radiation doses as well as faster reconstruction. But, by doing this, the original inverse problem becomes more and more ill-posed. Severe streak artifacts in CT images can be observed when using for instance the Filtered Back Projection (FBP) reconstruction algorithm. Compared with FBP algorithm, iterative methods can improve the image quality by modeling the measurement statistics and the imaging geometry. However, due to the underdetermined linear system with sparse projections, the reconstruction obtained this way is still far from being satisfactory [START_REF] Gordon | Algebraic reconstruction techniques for three-dimensional electron microscopy and photography[END_REF]. Recently, a special attention has been devoted to compressed sensing (CS) theory [START_REF] Candes | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF][START_REF] Dong | Exponential Wavelet Iterative Shrinkage Thresholding Algorithm for Compressed Sensing Magnetic Resonance Imaging[END_REF][START_REF] Donoho | Compressed sensing[END_REF][START_REF] Zhang | Exponential Wavelet Iterative Shrinkage Thresholding Algorithm with Random Shift for Compressed Sensing Magnetic Resonance Imaging[END_REF]. Studies on CS theory unveil the possibility of recovering from sparse signals when the requirement specified by the Nyquist sampling theorem cannot be satisfied [START_REF] Candes | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF]. CS has been extensive applied in pattern recognition and computer vision, such as machine learning [START_REF] Liu | A neural network for ℓ1-ℓ2 minimization based on scaled gradient projection: Application to compressed sensing[END_REF][START_REF] Wang | Magnetic Resonance Fingerprinting with compressed sensing and distance metric learning[END_REF], sparse coding and regression [START_REF] Lee | Efficient sparse coding algorithms[END_REF][START_REF] Deng | Sparse Coding Inspired Optimal Trading System for HFT Industry[END_REF][START_REF] Zhu | Convolutional sparse coding for trajectory reconstruction[END_REF][START_REF] Tropp | Algorithms for simultaneous sparse approximation, Part I: Greedy pursuit[END_REF][START_REF] Lee | Efficient L1 regularized logistic regression[END_REF], image segmentation [START_REF] Maiora | Random forest active learning for AAA thrombus segmentation in computed tomography angiography images[END_REF][START_REF] Kong | Discriminative Clustering and Feature Selection for Brain MRI Segmentation[END_REF], and so on [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF][START_REF] Shen | Image reconstruction algorithm from compressed sensing measurements by dictionary learning[END_REF]. In the field of medical imaging, CS theory has been successfully considered in Magnetic Resonance Imaging (MRI) with sparse sampling [START_REF] Lustig | Sparse MRI: The application of compressed sensing for rapid MR imaging[END_REF][START_REF] Trzasko | Highly undersampled magnetic resonance image reconstruction via homotopic ell-0minimization[END_REF], spectrally resolved bioluminescence tomography [START_REF] Lu | Source Reconstruction for Spectrally-resolved Bioluminescence Tomography with Sparse A priori Information[END_REF], optical coherence tomography [START_REF] Fang | Fast Acquisition and Reconstruction of Optical Coherence Tomography Images via Sparse Representation[END_REF].

CS theory assumes that sparse-distributed gradient information can be used to regularize tomographic reconstruction by minimizing a constrained TV function [START_REF] Candes | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF][START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. Although the restricted isometry property (RIP) condition is not often satisfied in practice, CS theory based reconstruction can yield more satisfactory results than the traditional FBP algorithm in sparse-view CT reconstruction [START_REF] Yu | A soft-threshold filtering approach for reconstruction from a limited number of projections[END_REF]. In [START_REF] Sidky | Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT[END_REF], Sidky et al. developed a TV minimization algorithm for 3-D divergent beam CT reconstruction with limited projections. Chen' group suggested a prior image constrained compressed sensing (PICCS) to reconstruct dynamic cardiac CT images from highly under sampled projections, where the gradient difference between the target image and a high-quality prior image is incorporated [START_REF] Chen | Prior image constrained compressed sensing (PICCS): A method to accurately reconstruction dynamic CT image from highly undersampled projection data sets[END_REF][START_REF] Lauzier | Prior image constrained compressed sensing: implementation and performance evaluation[END_REF]. Lu et al. proposed a dual dictionary algorithm for CT image reconstruction from limited projections [START_REF] Lu | Few-view image reconstruction with dual dictionary[END_REF]. It was also pointed out that, as to a piecewise constant object, a region of interest (ROI) could be exactly and stably reconstructed via the TV minimization [START_REF] Yu | Compressed sensing based interior tomography[END_REF]. Recently, Liu et al. proposed a total TV-stokes based method to improve sparse-view Xray CT reconstruction [START_REF] Liu | Total Variation-Stokes Sparse-View X-ray CT Image Reconstruction[END_REF]. Liu et al. proposed an algorithm, which combines a first order primal-dual and a projection onto convex sets (FOPD-POCS) to solve the TV minimization problem [START_REF] Liu | X-ray computed tomography using sparsity based regularization[END_REF].

Theoretically, the l 0 -norm can provide a sparser representation than the TV. But the reconstruction using l 0 -norm regularization is a NP hard problem due to its nonconvex function in discontinuous form [START_REF] Candes | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF]. The application of l 0 -norm regularization has already been explored in sparse-view CT reconstruction. Guo et al. proposed an EdgeCS (edge guided CS reconstruction) where a truncated TV regularization is applied to preserve those edges with gradient values larger than some specified threshold values [START_REF] Guo | EdgeCS: Edge Guided Compressive Sensing Reconstruction[END_REF][START_REF] Guo | Edge guided reconstruction for compressive imaging[END_REF]. Lee et al. formulated the imaging problem as a joint sparse recovery problem under CS framework and proposed an exact inversion algorithm that achieved the l 0 optimality as the measurement rank increases [START_REF] Lee | Compressive Diffuse Optical Tomography: Non-Iterative Exact Reconstruction using Joint Sparsity[END_REF]. Hu et al. developed an approximate l 0 -norm constrained sparse reconstruction algorithm for 3-D CT reconstruction [START_REF] Hu | L0 constrained sparse reconstruction for multi-slice helical CT reconstruction[END_REF]. Zhang et al. proposed a Gamma regularized CT reconstruction [START_REF] Zhang | Gamma regularization based reconstruction for low dose CT[END_REF].

The contributions of the present work can be featured as follows. Firstly, the Gamma regularization is extended to the sparse-view CT reconstruction to get an improved reconstruction performance over the general l 1 -norm based TV reconstruction. Secondly, the convexity property of the regularization term is thoroughly analyzed to show that the cost function can be monotonically decreased with ensured local minimum solution in the proposed algorithm. Thirdly, the proposed reconstruction is shown not sensitive to initial conditions. Experiments point out that the Gamma regularization leads to better reconstruction results than the cases using l 2 -norm and TV regularizations.

In section 2, the sparse-view CT reconstruction model and Gamma regularization are reviewed. The optimization of Gamma regularization based sparse-view reconstruction model and the related parameters setting are given in section 3. The performance of our model is assessed on simulated and clinical data in section 4. A discussion is provided in section 5 and concluding remarks and perspectives are sketched in section 6.

Sparse-view reconstruction with Gamma regularization.

Similar to the normal-view CT reconstruction, sparse-view reconstruction is also the process of finding the attenuation coefficients of the scanned object using available projection data. Intuitively, the mathematical formula of sparse-view CT reconstruction problem can be expressed as follows:

Gf y 

where f represents the image to reconstruct corresponding to the attenuation coefficients, and y refers to the projection measurements (or sonogram data). The system matrix G specifies the particular CT imaging geometry [START_REF] Peters | Algorithms for fast back-and re-projection in computed tomography[END_REF][START_REF] Zhuang | Numerical evaluation of methods for computing tomographic projections[END_REF][START_REF] Joseph | An improved algorithm for reprojecting rays through pixel images[END_REF][START_REF] Siddon | Fast calculation of the exact radiological path of a three-dimensional CT array[END_REF][START_REF] Man | Distance-driven projection and backprojection in three dimensions[END_REF].

Iterative reconstruction algorithms have received more attention over the last years. Among all the iterative reconstruction models, Least Mean Square (LMS) model is the simplest one, and can be formulated as follows:

2 2 1 arg min 2 () f f Gf y f         (2)
where () f  is the cost function. In standard CT reconstruction, due to the highly underdetermined property of G , directly solving (2) to get f is already an ill-posed problem. In sparse-view reconstruction, the reduction of projection views makes (2) much more underdetermined and leads to unstable and non-unique solutions.

CS theory can be used to overcome the ill-posed nature of (2) by adding a regularization term based on certain assumptions. The TV norm regularization with a sparse prior on image gradients has been reported in the 3-D sparse-view CT reconstruction [START_REF] Sidky | Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT[END_REF]. However, the performance of TV norm regularization significantly decreases as the number of projections decreases [START_REF] Candes | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF]. Recently, by utilizing the l 0 -norm approximation inherent to cumulative distribution function (CDF) of the Gamma distribution, Zhang et al. proposed a Gamma based regularization model for low dose CT reconstruction. It has been shown that this model yields better reconstruction performance in low dose CT imaging than the l 1 -norm based TV regularization [START_REF] Zhang | Gamma regularization based reconstruction for low dose CT[END_REF][START_REF] Hogg | Introduction to Mathematical Statistics[END_REF]. The corresponding Gamma potential function can be written:

    1 0 ;, t x te x dt              (3)
where     denotes the Gamma function. α and β are called respectively the shape and the scale parameters. By taking the gradient term in (3), the Gamma regularization function can be expressed by:

      , , 1
11 0

1 11 0 v ij h ij t f mn ij t f mn ij te f dt te dt                          (4)
where m and n are the two dimensional sizes of the 2D image to reconstruct.

,

v ij f  and , h ij f
 denote the vertical and horizontal direction gradients, respectively. By combining ( 2) and ( 4), the cost function for the sparse-view reconstruction with Gamma regularization is given by:

  2 2 1 arg min () 2 f f G y f f f             ( 5
)
where λ is the positive balance parameter.

Optimization.

Optimization.

The conjugate gradient method (CG) [START_REF] Boyd | Convex Optimization[END_REF][START_REF] Fletcher | Function minimization by conjugate gradients[END_REF] is used in the optimization of sparse-view reconstruction model [START_REF] Hsieh | Adaptive streak artifact reduction in computed tomography resulting from excessive X-ray photon noise[END_REF]. The partial derivative of ( 5) with respect to image element , ij f contains two terms. The derivative of the data fidelity term can be written:

    2 2 , , 1 2 = T ij ij Gf y G Gf y f     ( ) (6)
while the derivative of the Gamma regularization term is given by [START_REF] Zeng | A filtered back projection MAP algorithm with non-uniform sampling and noise modeling[END_REF]:

    , ij f f                     1 2 2 1, , 1 2 2 1, , 1 2 2 1, , 1, , i j i j ff i j i j i j i j i j i j f f e f f f f                                                1 2 2 , 1 , 1 2 2 , 1 , 1 2 2 , 1 , , 1 , i j i j ff i j i j i j i j i j i j f f e f f f f                                                1 2 2 1, , 1 2 2 1, , 1 2 2 , 1, , 1, i j i j ff i j i j i j i j i j i j f f e f f f f                                                1 2 2 , 1 , 1 2 2 , 1 , 1 2 2 , , 1 , , 1 i j i j ff i j i j i j i j i j i j f f e f f f f                                (7)
where ε is a positive constant. As a consequence, the derivative with respect to the proposed model ( 5) can be expressed by:

        , , , ()
T ij i j i j f G Gf y f f f            ( ) (8)
The overall solution is outlined in Algorithm 1 where the stopping criterion is jointly determined by the minimum norm of gradient 0  and the preset maximum iteration number 0 N . The Armijo-Goldstein condition based step setting parameters , ,  are specified in Table 1 [START_REF] Frank | An algorithm for quadratic programming[END_REF][START_REF] Armijo | Minimization of functions having Lipschitz continuous first partial derivatives[END_REF].

Algorithm 1: Sparse-view reconstruction with Gamma Regularization Input : projection data y

Output: reconstructed image f Set  ,  ,  ,  ,  , , 0 N , 1 k  Initialize 0 f , 0 g , 0 d while 0 k g   or 0 kN  while ( ) ( ) ( ) T kk f d f f d f                 end 1 k k k f f d    11 () kk gf   2 1 2 2 2 k k g g    11 k k k d g d      1 kk  end 3.2. Gamma parameters setting.
The Gamma regularization term contains two parameters, the shape parameter α and the scale parameter β. Their setting strategies have been discussed in [START_REF] Zhang | Gamma regularization based reconstruction for low dose CT[END_REF]. An interesting finding is that the two parameters affect the approximation to l 0 -norm in opposite ways. Based on that, the parameters α and β can be set by: (i) applying the FBP method to get a FBP reconstructed image for parameter estimation; (ii) calculating the 25% quantile of the gradient amplitude distribution of the FBP reconstructed image; (iii) setting α to 1.2 and 5E (E is the expectation the corresponding Gamma distribution with α and β, i.e. E=α/β) to the acquired quantile value to solve β. More details on the choice of these two parameters can be found in [START_REF] Zhang | Gamma regularization based reconstruction for low dose CT[END_REF]. The experiments conducted in this study validate that this parameter setting strategy is effective in providing good reconstruction results.

Experiment setting and results.

We conducted experiments on two numeric phantom data and one clinical pelvis projection data. The implementation runs under Matlab 2012a, 64Bit Microsoft Windows 7, Intel Core™ 2.40GHz (Quad-Core), 8Gb RAM. The proposed Gamma regularization is compared to the l 2 -norm and TV regularizations in sparse-view reconstruction. The initialization is performed with a zero-valued image and the parameter setting indicated in Table 1. The values of α and β were selected by applying the procedure described in section 3.2. λ was adjusted to get the best results in terms of PSNR (peak signal to noise ratio):

    2 max 2 , , 1 1
PSNR , 10log10 1

m n i j i j i j P PI PI mn                ( 9
)
where I is the reconstructed image and P the reference image. max P denotes the maximum intensity value in the reference image P . The reference image P is set as the phantom image for the simulated data, and the image reconstructed with TV regularization from the complete projection views for the clinical data.

Simulation Settings.

The numeric phantoms include a Modified Shepp-Logan (MSL) phantom (Fig. 1 (a)) and a simulated nonuniform rational B-splines (NURBS) based cardiac-torso (NACT) phantom (Fig. 1 (b)) [START_REF] Segars | Study of the efficacy of respiratory gating in myocardial SPECT using the new 4-D NCAT phantom[END_REF]. Both phantoms are composed of 256×256 pixels with 1mm×1mm pixel size. A monochromatic single slice CT parallel scanner with a total of 180 scanning views and 367 radial bins per view is modeled in the simulation. As illustrated in Fig. 2, we simulated uniformly-spaced 10, 12, 15, 18 projections for the MSL phantom, and uniformly-spaced 8, 9, 10, 11 projections for the NACT phantom. The total iteration number was set to 5000 to ensure a stable convergence in reconstruction. 2 and3. As it can be seen from ( a1), (b1), (c1), (d1) and the zoomed ROI (az1), (bz1), (cz1), (dz1) in Fig. 3 and Fig. 4, the l 2norm regularization fails to provide satisfying results in sparse-view reconstruction (some structures are blurred in the reconstruction). The reconstructions obtained with the Gamma and the TV regularizations lead to images with better visual quality than those resulting from the reconstruction based on the l 2 -norm regularization. When compared with the TV regularization, the proposed Gamma regularization provides a better contrast preservation (see the arrows in (bz3), (cz3) in Fig. 3 and the arrows in (az3)-(dz3) in Fig. 4). The performances of all the methods decrease as the number of projections decreases. The results of the TV regularization suffer from structure blurring and increased artifacts when the projection number is less than 12 for the MSL phantom and 10 for the NACT phantom (see the (az2), (bz2), (cz2), (dz2) in Fig. 3 and Fig. 4). However, even in these situations, the robustness of the proposed Gamma regularization appears much better as it is shown in the images (cz3) of Fig. 3 and Fig. 4).

l 2 -norm TV Gamma (a1) (a3) (a2) (b3) (b2) (b1) (b1) (a2) (a1) (a3) (b3) (b2) (a4) (b4)
l 2 -norm TV Gamma (c3) (c2) (d3) (d2) (d1) (az1) (az3) (az2) (bz3) (bz2) (bz1) (cz1) (cz2) (cz3) (dz1) (dz2) (dz3) (a1) (a3) (a2) (b3) (b2) (b1) (c1)

Quantitative analysis of the simulation data.

The PSNR values together with the corresponding regularization parameters of the reconstructed images calculated via [START_REF] Chen | Nonlocal prior Bayesian tomographic reconstruction[END_REF] are reported in Tables 2 and3. These results confirm the observations made through visual image inspections. The PSNR values decrease

(c3) (c2) (c1) (d3) (d2) (d1) (az1) (az3) (az2) (bz1) (bz2) (bz3) (cz3) (cz2) (cz1) (dz3) (dz2) (dz1)
8 views as the view number decreases and higher PSNR values are obtained for the Gamma regularization. Note that, for the two phantoms, α was set to 1.2 for the two phantoms and β according to the strategy described in section 3.2. [START_REF] Zhang | Gamma regularization based reconstruction for low dose CT[END_REF].

In Fig. 6 and Fig. 7, we plot the horizontal intensity profiles across the central lines of different reconstructed images (Fig. 5). Although the TV regularization performs already well, the Gamma reconstructions appear slightly better when matched to the ground truth phantom profiles.

Experiments on clinical pelvis projection data

The proposed approach was also tested on a clinical pelvis projection data of a 56-year-old male patient acquired on a monoenergetic Siemens Somatom CT system. The tube voltage and the current are 120 KVp and 100 mA, respectively. The image to reconstruct is composed of 256×256 pixels with 1mm×1mm pixel size. The imaging geometry is set as follows: the distance of source to center of rotation is 408 mm, and the distance of source to detector is 816 mm. A monochromatic single slice CT fan-beam scanner is imaged with a total of 360 scanning views and 512 radial bins per view. The acquired complete projection data is illustrated in Fig. 8 (a). For quantitative analysis purpose, a reference image (shown in Fig. 8 (b)) is reconstructed by means of a TV regularization using the complete projection data. Two ROIs (ROI1 and ROI2) marked in red and green rectangles are specified in Fig. 8 (b).

Fig. 9 to Fig. 12 illustrate the reconstruction results for l 2 -norm, TV, and Gamma regularizations with 1000 iterations for the uniformly-sampled 180, 120, 90 and 60 views, respectively. The regularization parameters λ, α and β are given in the corresponding figure captions below. We also display the two zoomed ROIs for comparison and we can observe much less strip artifacts in the reconstruction using Gamma regularization than other methods. Fig. 13 plots the PSNR values (with respect to the reference image in Fig. 8 (b)) for the whole image and two ROIs. We can see that the proposed method always leads to the highest PSNR among all the methods when the total sampling view is lower than 180. Some ringing artifacts can be observed in Fig. 4 (dz3) and Fig. 8-12 around the edge structures in the reconstructed images. As a matter of fact, the generated ringing artifacts are composed of the streak artifacts' residues. The main reason is due to the regularization effects that are related to the combination of streak artifacts and the inherent image edges,which can be analyzed as follows. On the one hand, when compared to the l 2 -norm or TV regularization based algorithms, the proposed Gamma regularization significantly reduces streak artifacts in the homogeneous areas of the reconstructed images, whether it is the bright bone areas or the dark areas. On the other hand, in the high contrast connected areas, it cannot well distinguish the bright streak artifacts from bright bone areas as the l 2 -norm or TV regularization based

Discussion.

The RIP condition for the TV regularization

The RIP condition defined in [START_REF] Candes | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF][START_REF] Candes | The restricted isometry property and its implications for compressed sensing[END_REF], is a key element in the CS theory. It is an important property of the sample matrix that evaluates the deviation between the reconstruction results and the optimal sparse approximation. This deviation for l 1 -norm regularization has been analyzed in [START_REF] Candes | The restricted isometry property and its implications for compressed sensing[END_REF].

Different from the canonical regularization term within the CS reconstruction model, TV regularization is used to constrain the total gradient fluctuations of the image to reconstruct. TV regularization has been widely employed in CS for image reconstruction [START_REF] Candes | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF]. Furthermore, RIP condition of the sample matrix for the TV regularization has also been analyzed in [START_REF] Needell | Stable image reconstruction using total variation minimization[END_REF] together with the sample condition in order to obtain satisfying reconstructions. More precisely, the relation between the RIP of the sample matrix and the reconstruction error has been established in [START_REF] Needell | Stable image reconstruction using total variation minimization[END_REF]. When compared to the l 1norm regularization, the total reconstruction error of TV regularization is not only up to the error between the gradients of the reconstruction results and the optimal sparse approximation, but also to the logarithmic value related to the pixel number in the target image [START_REF] Needell | Stable image reconstruction using total variation minimization[END_REF].

Why the Gamma regularization works better than l 1 -norm?

The reason why the l 0 -norm regularization is better than the l 1 -norm in the case of highly limited measurements has been already theoretically proven in CS [START_REF] Candes | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF]. Therefore, because the Gamma regularization approaches closer to the l 0 -norm regularization than the l 1 -norm one in those points with larger gradients, it works better than l 1 -norm regularization in sparse-view reconstruction. This mechanism in modulating regularization effect is closely consistent with the intuitive observation that large gradients are often related to image edges while small gradients often come from noise and artifacts. A better preservation of edge structures with improved reconstruction quality can thus be achieved with the Gamma regularization.

Convergence analysis.

The first fidelity term in ( 5) is a function with strict convexity. However, a global minimum of the objection function cannot be guaranteed with Algorithm 1, as explained in the convexity analysis of the regularization term in the Appendix below.

Nevertheless, according to the Armijo-Goldstein based backtracking line search step length setting strategy applied in our CG algorithm [START_REF] Zhang | Gamma regularization based reconstruction for low dose CT[END_REF][START_REF] Frank | An algorithm for quadratic programming[END_REF][START_REF] Armijo | Minimization of functions having Lipschitz continuous first partial derivatives[END_REF], the following relation can be easily established: [START_REF] Chen | Thoracic low-dose CT image processing using an artifact suppressed large-scale nonlocal means[END_REF] where k f denotes the th k iterated image in reconstruction, and

(c) (b)                     2 1 1 11 2 1 2 () 1 2 () k k vh kk k k vh kk Gf y f ff Gf y f ff                                      
1 k f  denotes 1 th k  ( )
. It can be seen that a monotonic decrement in iteration can be obtained for Algorithm 1, which means that local minimum can be reached. The convergence behavior of our algorithm is similar to the majorization-minimization (MM) algorithm investigated in [START_REF] Hunter | A tutorial on MM algorithms[END_REF]. As a consequence, with the condition of accepted error tolerance, a local minimization can be carried out in practice. We analyzed the monotonic property of the overall cost function [START_REF] Hsieh | Adaptive streak artifact reduction in computed tomography resulting from excessive X-ray photon noise[END_REF] in Algorithm 1 via the reconstructions of the MSL phantom with 15 views and the NACT phantom with 10 views. All the reconstruction parameters were set according to Table 1. The evolutions of the cost function values and the PSNR values over iterations are displayed in Fig. 14 (a) and Fig. 14 (b). We can observe monotonic decreases of the cost function values and monotonic increases of the PSNR over the whole 5000 iterations, and this shows that the proposed algorithm provides a stable solution for sparse-view reconstructions. To explore the impact of the iterations, we also conducted experiments with the MSL phantom with 15 views and the NACT phantom with 10 views. The iteration numbers are set to 5000 and 10000, respectively. The reconstructed images are shown in Fig. 15. It can be seen that the reconstructed images are not impacted by the iterations.

Initial condition analysis

Although the Gamma norm is not convex, our algorithm always converges to an ensured local minimization. To examine the influence of the initialization (i.e. zero-valued initial image or FBP initial reconstruction), we analyzed the reconstructions of the MSL phantom with 12 views and the NACT phantom with 10 views. The reconstruction results are depicted in Fig. 16. The balance parameter λ was here selected to get the reconstruction with the PSNR value. We can see that nearly the same reconstructed images with same PSNR values can be obtained for these two different initial conditions. Therefore, the proposed reconstruction method does not appear sensitive to initial conditions.

Computation complexity analysis.

The proposed method needs about 9 minutes for the simulated data experiments (5000 iterations) and the clinical data experiments (1000 iterations) which is more time-consuming than the reconstruction using TV (about 2 minutes) and l 2 -norm (less than 1 minute) based regularizations. This is mainly due to the CDF calculation of Gamma regularization in the step size searching. Approximations of the CDF will be studied in our future work to speed up the computation efficiency.

Comparison with Hu' l 0 -norm regularization.

In [START_REF] Hu | L0 constrained sparse reconstruction for multi-slice helical CT reconstruction[END_REF], Hu et al. proposed an l 0 -norm based helical CT reconstruction method. The so-called l 0 -norm is in fact an approximate mathematical l 0 -norm obtained by introducing a potential function:

 , log (
1) t t    ( 11
)
where ρ is a positive parameter. We compared the performance of the proposed Gamma regularization with Hu's l 0 -norm regularization in our sparse-view CT reconstruction method. The MSL phantom with 12 views and the NACT phantom with 10 views were used for validation. All the involved parameters in Hu' l 0 -norm regularization were set to get the best results in terms of PSNR. The picture in Fig. 17 shows that our algorithm leads to reconstructed images with better edge preservation than Hu's method for sparseview CT reconstruction (see the arrows in the zoomed images in Fig. 17).

Fig. 17. The reconstruction images with Hu's method and the proposed method. (a1) Hu' method, λ=0.5, PSNR=31.37; and (b1), Hu' method, λ=0.3, PSNR=30.98, (a2) the proposed method, λ=0.1, PSNR=33.03, (b2), the proposed method, λ=0.01, PSNR=32.58.

Conclusion.

In this paper, Gamma regularization proposed in [START_REF] Zhang | Gamma regularization based reconstruction for low dose CT[END_REF] was employed for sparse-view CT reconstruction. This Gamma regularization provides a tunable strategy of fractional order norm to improve the reconstruction. The performance of the proposed approach was assessed on simulated phantoms and real clinical data. In the simulations on the MSL phantom, we conducted the reconstructions using 10, 12, 15 and 18 views. For the NACT phantom, the sampled 8, 9, 10 and 11 sinogram views were used. In the clinical data experiments, 180, 120, 90, 60 views were used. Comparisons with the TV and l 2 -norm regularizations were also conducted. Although more clinical applications are still required, the present visual and quantitative evaluations confirmed the benefits brought by our approach and that its advantage becomes more prominent as the projection number decreases. However, some ring artifact residuals remain in the reconstructed images. They are due to regularization effects related to the combination of streak artifacts and the inherent image edges. Operations like edge detection, segmentation [START_REF] Wang | Computer-aided Diagnosis of Abnormal Breasts in Mammogram Images by Weighted-Type Fractional Fourier Transform[END_REF] or classification [START_REF] Zhang | Statistical iterative reconstruction using adaptive fractional order regularization[END_REF] should be incorporated to our model capable to deal with these problems.

Currently, all the parameters involved in our approach were chosen according to the strategy reported in [START_REF] Zhang | Gamma regularization based reconstruction for low dose CT[END_REF] where the shape parameter α was fixed to 1.2, and the scale parameter β set according to the image distribution histogram issued from a preliminary FBP reconstruction. Our future work will develop adaptive strategies for setting these parameters in order to better use the intrinsic image data features, and doing so, we will explore the connection between Gamma regularization and factor regularization [START_REF] Wang | Computer-aided Diagnosis of Abnormal Breasts in Mammogram Images by Weighted-Type Fractional Fourier Transform[END_REF][START_REF] Zhang | Statistical iterative reconstruction using adaptive fractional order regularization[END_REF][START_REF] Zhang | Few-view image reconstruction with fractional-order total variation[END_REF]. Furthermore, this proposed method will be tested in other image processing direction, like segmentation [START_REF] Zheng | Image segmentation by generalized hierarchical fuzzy C-means algorithm[END_REF], classification [START_REF] Wen | A rapid learning algorithm for vehicle classification[END_REF], and advanced medical imaging modalities, such as MRI [START_REF] Qu | Magnetic resonance image reconstruction from undersampled measurements using a patch-based nonlocal operator[END_REF][START_REF] Dong | Detection of subjects and brain regions related to Alzheimer's disease using 3D MRI scans based on eigenbrain and machine learning[END_REF][START_REF] Qu | Undersampled MRI reconstruction with patch-based directional wavelets[END_REF], et al.

Its applicability in clinical setting requires however a significant acceleration which will be achieved by means of parallel computation based on Graphics Processing Unit (GPU).

Appendix.

We analyze the concavity of the Gamma regularization function ( 5) by evaluating its second order derivative. We can see in [START_REF] Zeng | A filtered back projection MAP algorithm with non-uniform sampling and noise modeling[END_REF] that the first derivative of Gamma regularization function includes four terms, all with the same structure. The second order derivative can be derived by calculating the derivative of the first term in [START_REF] Zeng | A filtered back projection MAP algorithm with non-uniform sampling and noise modeling[END_REF]:

            1 2 2 1, , 1 2 2 1, , , 1,
,

i j i j ff i j i j i j i j ij f f e ff f                                                             1 2 2 1, , 1 2 2 1, , 1 2 2 1, , 2 22 2 1, , , 1, 3 22 2 1, , , 1, 1 2 2 1,
,

1 i j i j i j i j i j i j ff i j i j i j i j ff i j i j i j i j ff i j i j f f e f f e f f f f f f e                                                                                        1 2 2 1, , 2 2 2 1, , 1 2 2 2 2 , 1, 1, , , 1, 2 1, , 1 i
j i j ff i j i j i j i j i j i j i j i j i j i j f f e f f f f f f ff                                                               1 2 2 1, , 2 2 2 1, , 1 22 2 , 1, 1, , 11 i
j i j ff i j i j i j i j i j i j f f e f f f f                                                         1 2 2 1, , 2 2 2 1, , 1 22 2 , 1, 1,
,

i j i j ff i j i j i j i j i j i j f f e f f f f                                       (A1)
Similarly, we can obtain (A2), (A3) and (A4) for the remaining three terms and thus, the overall second order derivative of the regularization term is: From (A5), we can see that if From above we can see that (A8) cannot well hold to ensure the overall convexity of the cost function (5) because the s Max value is determined by the unpredictable intensity variation in the reconstructed images. So, the global minimum cannot be obtained through Algorithm 1, and a local minimum solution can only be attained.

            1 2 2 , 1 , 1 2 2 , 1 , , , 1
, i j i j ff i j i j i j i j ij f f e ff f                                           1 2 2 , 1 , 1 2 2 2 2 22 , 1 , , , 1 , 1 ,
i j i j ff i j i j i j i j i j i j f f e f f f f                                        (A2)             1 2 2 1, , 1 2 2 1, , , 1, , i
j i j ff i j i j i j i j ij f f e ff f                                           1 2 2 1, , 1 2 2 2 2 22 1, , , 1, 1,
                    1 2 2 1, , 2 , 1 2 2 2 2 22 1, , , 1, 1,
                                                       1 2 2 , 1 , 1 2 2 2 2 22 , 1 , , , 1 , 1 ,
1, , 1, , , 1 , , 1 , 0 i

Table 1 :Fig. 1 :Fig. 2 :

 112 Fig. 1: The numeric phantom data used in this paper. (a) MSL phantom, (b) NACT phantom.

Fig. 3 :

 3 Fig. 3: Reconstructed images of the MLS phantom. (a1)-(a3), (b1)-(b3), (c1)-(c3) and (d1)-(d3) are the reconstructed images using the l2-norm regularization, TV regularization and the Gamma regularization with 18 views, 15 views, 12 views and 10 views, respectively. (az1)-(az3), (bz1)-(bz3), (cz1)-(cz3) and (dz1)-(dz3) show the zoomed marked ROI in (a1)-(a3), (b1)-(b3), (c1)-(c3), and (d1)-(d3), respectively.

Fig. 4 :

 4 Fig. 4: Reconstructed images of the NACT phantom. (a1)-(a3), (b1)-(b3), (c1)-(c3) and (d1)-(d3) are the reconstructed images using the l2-norm regularization, TV regularization and the Gamma regularization with 11 views, 10 views, 9 views and 8 views, respectively. (az1)-(az3), (bz1)-(bz3), (cz1)-(cz3) and (dz1)-(dz3) show the zoomed marked ROI in (a1)-(a3), (b1)-(b3), (c1)-(c3), and (d1)-(d3), respectively.

4. 2 .

 2 Visual evaluation of the simulation data.

Fig. 3

 3 Fig. 3 and Fig. 4 display the images reconstructed for performance assessment and comparison. For each phantom, a zoomed region of interest (ROI) is chosen (delineated by red line) to highlight the local image structures. The PSNR values of the reconstructed images are summarized in Tables2 and 3. As it can be seen from (a1), (b1), (c1), (d1) and the zoomed ROI (az1), (bz1), (cz1), (dz1) in Fig.3and Fig.4, the l 2norm regularization fails to provide satisfying results in sparse-view reconstruction (some structures are blurred in the reconstruction). The reconstructions obtained with the Gamma and the TV regularizations lead to images with better visual quality than those resulting from the reconstruction based on the l 2 -norm regularization. When compared with the TV regularization, the proposed Gamma regularization provides a better contrast preservation (see the arrows in (bz3), (cz3) in Fig.3and the arrows in (az3)-(dz3) in Fig.4). The performances of all the methods decrease as the number of projections decreases. The results of the TV regularization suffer from structure blurring and increased artifacts when the projection number is less than 12 for the MSL phantom and 10 for the NACT phantom (see the (az2), (bz2), (cz2), (dz2) in Fig.3and Fig.4). However, even in these situations, the robustness of the proposed Gamma regularization appears much better as it is shown in the images (cz3) of Fig.3and Fig.4).

Fig. 5 :Fig. 6 :Fig. 7 :

 567 Fig. 5: The profiles (the 128 th row from the top) of the two phantoms.

  , these streak artifacts' residues are leaved around the bright bone areas and thus, ringing artifacts are formed. How to overcome this drawback, by means of postprocessing, such as segmentation, edge detection, etc., will be studied in our next work.

Fig. 8 :

 8 Fig. 8: (a) the complete pelvis projection data. (b) image reconstructed with TV regularization. Note that the two ROIs (ROI1 and ROI2) for more thorough analysis are delineated in red and green rectangles.

Fig. 9 :Fig. 10 :Fig. 11 :Fig. 12 :Fig. 13 :

 910111213 Fig. 9: Reconstructed pelvis CT images with 180 projections. (a)-(c): the images reconstructed by l2-norm regularization (λ=0.1), TV regularization (λ=0.1), Gamma regularization (λ=1.8, α=1.2, β=2.4), respectively. (az1)-(cz1): the corresponding zoomed ROI1 in the (a)-(c). (az2)-(cz2): the corresponding zoomed ROI2 in the (a)-(c). (dz1) and (dz2) are the corresponding zoomed ROIs from the reference image in Fig. 8 (b).

Fig. 14 :

 14 Fig. 14: The PSNR (a) and logarithmic cost function values (b) over iterations for MSL phantom with 15 views and NACT phantom with 10 views.

Fig. 15 .

 15 Fig. 15. The reconstruction images with different iterations in the experiments using the MSL phantom with 15 views and the NACT phantom with 10 views. (a1) reconstructed image with 5000 iterations for MSL phantom (λ=0.1), PSNR=45.11. (a2) reconstructed

Fig. 16 .

 16 Fig.16. The reconstruction images with different initial conditions for the experiment using the MSL phantom with 12 views and the NACT phantom with 10 views. (a1) reconstructed image using FBP initial condition for MSL phantom (λ=0.08), PSNR=33.03. (a2) reconstructed image using all zero initial condition for MSL phantom (λ=0.1), PSNR=33.03. (b1) reconstructed image using FBP initial condition for NACT phantom (λ=0.01), PSNR=32.58. (b2) reconstructed image using all zero initial condition for NACT phantom (λ=0.01), PSNR=32.58.

s

  Min denotes the minimum among the four values in S . In the same way, we can define s Max as the maximum among the four values in S and get (A8):

Table 2 :

 2 The PSNR (dB) of the reconstructed CT images with different view numbers for MSL phantom

	Phantom	Methods	10 views	12 views	15 views	18 views
		l2-norm	22.67 (λ=0.1)	24.10 (λ=0.1)	25.56 (λ=0.1)	26.07 (λ=0.1)
		TV	21.77 (λ=0.03)	25.28 (λ=0.03)	32.34 (λ=0.03)	38.67 (λ=0.05)
	MSL		(λ=0.1, α=1.2,	(λ=0.1, α=1.2,	(λ=0.1, α=1.2,	(λ=0.1, α=1.2,
		Gamma	β=8.55)	β=7.59)	β=6.67)	β=5.83)
			27.90	33.03	45.11	49.61

Table 3 :

 3 The PSNR (dB) of the reconstructed CT images with different view numbers for NACT phantom

	Phantom	Methods	8 views	9 views	10 views	11 views
		l2-norm	18.39 (λ=0.1)	19.05 (λ=0.1)	21.14 (λ=0.1)	21.95 (λ=0.1)
		TV	23.73 (λ=0.03)	28.68 (λ=0.03)	30.76 (λ=0.03)	34.16 (λ=0.03)
	NACT		(λ=0.01, α=1.2,	(λ=0.1, α=1.2,	(λ=0.01, α=1.2,	(λ=0.01, α=1.2,
		Gamma	β=9.07)	β=10.85)	β=9.65)	β=8.78)
			24.69	29.62	32.58	42.08

  ,

		 	 	f	i		j		f	i j						e      	i j ff  i j   	i j   i j i j i j f f f f                  	(A3)
	             1 2 2 , 1 , 1 2 2 , 1 , , , 1 i j i j ff i j i j i j i j f f e ff                    , ij     f  
		       	f	, 1 i j 		f	, i j		2				 22          1 2 2 , 1 , 1 2 2 2 , , 1 , 1 , i j i j ff i j i j i j i j e f f f f                              	(A4)

-norm regularization (a) (d) (c) (b) (a)(b)

-norm regularization (a) (c) (d) (b)

Acknowledgements

This research was supported by National Basic Research Program of China under grant (2010CB732503), Natural Science Foundation under grants (81370040, 31100713), and the Open Fund of China-USA Computer Science Research Center (KJR16026) and by the Qing Lan Project of Jiangsu Province.