Smoothed L0-Constraint Dictionary Learning for Low-Dose X-Ray CT Reconstruction

The iterative algorithms of computed tomography (CT) reconstruction derived from the dictionary learning (DL) regularization have been developed to make high quality recovery from the under-sampled data acquired by a low dose protocol. However, when they are applied to noisy data with low sampling rate, streaking artifacts and bias tends to appear in early iteration results. Since the dictionary is over-complete, the artifacts and bias can also be represented well by the dictionary, resulting in the reservation of these unexpected structures in the final image. We proposed a smoothed L0 norm-constraint dictionary learning (SL0-DL) algorithm to deal with these unexpected structures. For the proposed algorithm, we introduce smoothed-L0 norm regularization to the objective function. In each iteration process, the intermediate image generated by the DL representation will be smoothed using SL0 norm, and then the smoothed image is used to update the output of this iteration. The raw data from both the numerical simulation and actual CT acquisition are used to test the performances of the proposed SL0-DL method. Experimental results demonstrate that the proposed method performs better than other competing algorithms with better noise and artifacts suppression performance while preserving image texture details. And the results show a significant improvement in the quality of the reconstructed image, which demonstrates that the proposed algorithm is really effective.


I. INTRODUCTION
Nowadays, x-ray computed tomography (CT) is still one of the most widely used imaging modalities to provide patients' anatomical structure because of its high spatial resolution and quality of its image. However, it is well known that an overdose of radiation may increase the excess lifetime cancer risk or any other diseases associated with genetic modification [1]. Simply, there are two direct ways to reduce x-ray dose. First approach is to lower mAs in the CT data acquisition protocols and this will increase the quantum noise level and reduce the signal-to-noise ratio (SNR) of the out-The associate editor coordinating the review of this manuscript and approving it for publication was Wen Chen . put signal. The quality of the reconstructed image will be degraded by the noise-contaminated sinogram data. To reduce the sinogram noise, the penalized weighted least-squares (PWLS) approach solves this problem in two dimensions, where the WLS considers the first and second order noise moments and the penalty models the signal spatial correlations [2]. Another way is to reduce the number of projection angles and the image reconstructed by the traditional analytic reconstruction algorithm (FBP etc.) will have a lot of streak artifacts [3]. The analytic-based algorithms are derived from a continuous imaging model so that it really needs dense sampled projections following the Nyquist/Shannon sampling theorem [4]. For the reconstruction problem of incomplete projection data, algebraic algorithms like the simultaneous VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ algebraic reconstruction technique (SART) [5] performs better compared to analytical algorithms. In recent years, compressed sensing theory [6], [7] provides new ideas and theoretical basis for CT image reconstruction, making an improvement on this algebraic algorithm. Compressed sensing (CS) theory shows that various signals can be recovered accurately by far fewer samples than the signal dimension. And Wang and Shim [8] believe that the lossless compressed signal can be completed when the condition of Restricted Isometry Property (RIP) is satisfied. Many compressed sensing-based algorithms have been published and applied to few-view CT reconstruction problem. Most of these methods applied the total variation (TV) [9]- [11] as the regularization constraint. The most commonly used sparse transform based on the TV method is the DGT (Discrete Gradient Transform) of the image. Algorithms like adaptive steepest descent projection onto convex sets (ASD-POCS) [9] and gradient projection Barzilai Borwein (GPBB) [10] are typical TV-based reconstruction algorithms. Another group of these methods generate the sparse constraint based on dictionary learning (DL), which extract many small overlapping patches from the image, and these image patches can be sparsely represented by an overcomplete dictionary. The typical examples are global dictionary based statistical iterative reconstruction (GDSIR) and adaptive dictionary based statistical iterative reconstruction (ADSIR) [11]. Soltani et al. [12] present a tensor dictionary learning method that focuses on solving the tomographic image reconstruction problem. When it comes to the spectral CT reconstruction problem, Zhang et al. [13] propose a tensor-based dictionary learning (TDL) method to sparsely represent the tensor image patches during the iterative process. Since the spectral CT image can be sparsely represented in each of its multiple energy channels, and the channels are highly correlated, the TDL method can produce superior image quality and leads to more accurate material decomposition. Chen et al. [14] develop a platform-independent discriminative dictionary representation method to mitigate and reduce CT truncation artifacts directly in image domain. This discriminative dictionary is composed of an artifact sub dictionary and a non-artifact sub dictionary, which were constructed to achieve selective and sparse representation of artifact and non-artifact components of a CT image respectively. In order to improve the performance of the previous dictionary learning, Lin et al. [15] propose a new robust, differentiated and comprehensive dictionary learning model (RDCDL) in order to fully represent the versatility, particularity and interference of data. Recently, with the progress of deep learning, many scholars have applied deep learning based methods into the field of medical imaging [16]- [21]. Compared with dictionary learning, the image quality by deep learning has been improved to some extent. But it relies heavily on the quality of the training set, and most of the training data currently used in deep learning are simulated via noise insertion, which may not reflect the realistic noise distribution of low dose CT scan. In contrast, dictionary learning does not require a large amount of training data, and can be easier to be integrated into the traditional iterative reconstruction framework.
However, there is one inevitable challenge of a regular DL method especially when the projection data are under-sampled and corrupted with noise. During the iteration process of a regular DL method, since the noise affects the consistency of the projection data and the sampling rate is very low, therefore irregular streaking artifacts and bias will definitely appear in the output results of the early iterations. Depending on the theory of DL method, the result is regularized by the sparse representation according to the over-complete dictionary. However, both the original structures of the image and the unexpected artifact bias can be well sparsely represented by the same dictionary, resulting in the decline in image quality caused by the presence of unexpected structures in the final output image. To solve the sparse-view projection data reconstruction problem, our previous work [22] and Wu et al. [23] introduce L1-norm DL regularization term and L1-norm sparsity constraint respectively into the statistic iterative reconstruction (SIR) framework for the reconstruction problem of sparse-view projection measurements. Mohimani et al. [24] propose a smoothed L0 norm algorithm that approximates the discrete L0 norm with a suitable smoothing function, and conduct a sparsity smoothing experiment. However, all these mentioned methods cannot suppress the artifact associated with regular DL method effectively. Inspired by the fact that the smoothed L 0 norm-constraint can smooth the unexpected structures, we therefore want to find a balanced solution that utilizes the smoothed L 0 norm (SL0) constraint for suppressing the streaking artifacts and preserving the low-contrast edge information at the same time. Therefore, we bring in the SL0 constraint DL regularization to the reconstruction model, in which the intermediate image generated after the dictionary-based sparse representation is smoothed by SL0 constraint before the output image updating step. By adding this smoothing step, we anticipate that it could help the DL method to suppress the streaking artifacts. In order to be sure that the method will have mild impact on the real image edge structures, we perform the smoothing step directly on the intermediate image result generated by DL-based representations rather the reconstructed image.
Therefore, the main contribution of this paper is to propose a smoothed L0 norm-constraint dictionary learning (SL0-DL) algorithm to deal with these unexpected artifacts generated by the DL representation during iteration process. The rest of the article is arranged as follows. In Section II, the paper introduces the reconstruction model and the traditional iteration process description of the SIR and DL methods, and then provides the new version of iteration process and the proposed SL0-DL model, followed by the parameters selection strategy and finally the information of the simulation phantom studies and real animal study are introduced, and the quantitative metric is provided. In Section III, the experimental results are displayed to illustrate the proposed algorithm's efficacy.
Section IV gives the corresponding discussions and further analysis and the final conclusion is presented at the end of this paper.

II. METHODS AND MATERIALS
A. STATISTICAL ITERATIVE RECONSTRUCTION WITH DICTIONARY LEARNING (SIR-DL) MODEL Generally, we can build a series of linear equations to describe the relationship between CT images and projection data where A = {a ij } ∈ R I ×N 2 represents the geometry of the data acquisition system; the unknown image µ is the linear attenuation coefficient distribution of the X-ray corresponding to the CT image; and g = (g 1 , g 2 , . . . , g I ) ∈ R I ×1 is the logtransformed projection data. Since this is an underdetermined equation with undersampled projection data, the algorithm based on DL reconstructs the CT image by minimizing the energy function with a DL regularization term as where the above optimization problem is proposed by combining statistical iterative reconstruction (SIR) with DL term according to previous work [11]. λ is the regularization parameter; ω i = (y i − γ i ) 2 /y i is the statistical weight for each X-ray projection, y = (y 1 , y 2 , . . . , y I ) T ∈ R I ×1 are the projection data detected by the detectors, γ i represents the scattering and background noise;ĝ = (ĝ 1 ,ĝ 2 , . . . ,ĝ I ) T ∈ R I ×1 is an estimate of g; E s = {e s nj } ∈ R N 2 o ×N 2 is an operator to extract patches with N 0 × N 0 pixels from the image; D = (d 1 , d 2 , . . . , d K ) ∈ R N 2 0 ×K is the over-complete dictionary, the dictionary is over-complete (N 2 0 K ); α = {α s ∈ R K ×1 } are the sparse representations of the corresponding patches under the dictionary basis D; ν s is a regularization parameter relative to the dictionary training process.

B. ITERATION PROCESS DESCRIPTION
In this work, we only consider the situation that the K-means Singular Value Decomposition (K-SVD) [25] algorithm is utilized to train the dictionary D. Since µ and α s are both unknown in the optimization problem (2), the traditional plan of converging to an optimal result is an alternating minimization scheme. The iterative process is divided into two alternate steps: update of sparse coefficients and CT images. The final result is obtained when the iterative process reaches a stopping criterion.
When updating the patch sparse representations α, supposing the image µ is fixed, the dictionary learning problem is presented as (3).
The above formula is a sparse expression problem. The orthogonal matching pursuit (OMP) algorithm can be utilized to solve the sparse representation coefficients of image patches [26]. The OMP algorithm first selects the atom with the largest product of the iteration residuals (observed residuals) from the compression matrix (D) according to the principle of selecting atoms in the matching tracking (MP) algorithm, and then the selected atoms are subjected to Gram-Schmidt Orthogonalization to obtain the space formed by these orthogonal atoms. Next, the sampling vector E s µ, (s = 1, . . . , S) is projected onto this space, which can be used to obtain the components and iteration margin of the sampling vector in the orthogonal space. Finally, use the same method to decompose the residual, use the least square method to find the best matching atom, and loop until it approaches the original sampling vector. When updating the reconstructed image, the sparse representation coefficients remain constant; the optimization problem is expressed as follows: The data fidelity term can be replaced by a separable paraboloid surrogate [27], and then the optimization can be iteratively resolved by (5), as shown at the bottom of the next page.

C. NEW VERSION OF ITERATION PROCESS
In general, the training dictionary is over-complete, which is conducive to improve the sparse representation of image patches. However, low-dose CT reconstructions will have a large number of artifacts and unknown structures in early iterations. This property of dictionary learning results in that these artifacts and bias can also be well represented, which brings into the final image and affects the image quality. Therefore, we propose a low-dose CT image reconstruction algorithm (SL0-DL) for suppressing artifacts of intermediate image generated by the DL representation.
As shown in Fig.1, in order to obtain the intermediate image using the dictionary-based sparse representation different from formula (5) that adopts the overall, optimization objective function (2), we optimize the data fidelity term and the regularization term separately.
Therefore, different from the traditional alternating minimization scheme, every iterative process of image reconstruction is divided into three steps.
Firstly, the image is updated by SIR process, and the optimization problem is shown as VOLUME 8, 2020 By using the separable paraboloid surrogate method [26], (6) can be iteratively resolved as where tp = 1, . . . , 30 is the SIR iterative index.
Secondly, based on the pre-defined dictionary, sparse representations α of the image patches are calculated by the orthogonal matching pursuit (OMP) algorithm as formula (3) does and then the sparse constraint regularization term becomes where µ old = µ tp .
Similarly, according to the work [11], the optimization of formula (8) can be iteratively resolved by Then according to the above formulas, an intermediate image is generated by using sparse representation α and the predefined dictionary D. Expanding the second item to get the following results and recall , where E T sj represents the transpose operator of E sj , and E sj is the vector of the j-th column of E s ; I N 2 0 is a vector of size N 2 0 with all its elements being 1. After that we have where x j represents the j-th pixel of the intermediate image x. Therefore x can be given by where./ is the operator that divides the corresponding element between the two vectors; E T s represents the transpose operator of E s .
Because both the original structures of the image and the unexpected artifact bias can be well sparsely represented by the same dictionary D, so x needs to be processed for noise reduction and artifact suppression. We use the smoothed L 0 norm (SL0) constraint to process x before the output image updating step. In compressed sensing (CS), the L0 norm is the most idealized regularity, but it is difficult to solve. Therefore, in practice, the L1 norm is usually used to replace the L0 norm. The smoothed L0 norm uses a continuous function to approximate the L0 norm to solve its discontinuous problem, which is theoretically beneficial to reconstruct high-quality CT images. It can be considered as the following minimum optimization problem [28]: . . , W, H and W are respectively the width and height of the image x; σ is a modulation parameter, which determines approximation degree.
Finally, the output of this iteration will be obtained by the formula (14) where x ts j is the j-th pixel of the smoothed intermediate image x ts , ts is the current iteration index of SL0 loop.

D. ALGORITHM IMPLEMENTATION
We added SL0 constraint to the intermediate image x, aiming to reduce the streaking artifacts and bias formed by the noisy under-sampled data, so that the objective function transforms to min µ,α,D Compared to the iteration process of the regular DL method in Section II.B, the first two steps in each iteration remain the same when solving the optimization problem (15) and before the third step which updates the output image, the process that smooths the intermediate image by adding SL0 constraint to the iteration. The pseudo-code of the SL0-DL model is shown as follows, initializes a random image µ 0 , as the predefined dictionary D, β 1 and β 2 are empirically selected values.
Moreover, parameters tp and ts are obtained through testing and experience selection. When the number of tp is small, the proposed method converges slowly. While tp > 30, not only the convergence speed has not increased significantly, but it will greatly increase the reconstruction time. As to ts, a smaller value will result in a lower denoising and artifact reduction effect, and a large number will easily produce some blurring effects. In both phantom and real data experiments, we use tp = 30 and ts = 5.

III. EXPERIMENTAL DESIGN RESULTS
In order to illustrate the effectiveness of the proposed algorithm (SL0-DL) in the reduction of artifacts and bias generated by the under-sampled noisy data. We executed SL0-DL on two different phantom simulation and real animal studies. Then it was compared with the gradient projection Barzilai Borwein (GPBB) [10], the total variation (TV) constraint SIR algorithm (SIR-TV) [29] and the traditional DL [11].
A. DATA SOURCES AND PARAMETER SELECTION 1) SIMULATED DATA In this paper, we utilized two different phantom simulation experiments to demonstrate the effectiveness of the proposed algorithm (SL0-DL). It consists of human head and abdomen slices shown in Fig. 2(a) and Fig. 3(a). The projection data are simulated using fan-beam geometry with the scanning circle spanning through 360 • range around the phantom and the number of angles is 240 with the equal sampling rate as 1.5 • . The head slice phantom covers a region of 20 × 20 cm 2 and the distance from the X-ray source to the system origin is 40cm. For each projection view, 732 detector elements span a FOV of 25cm in diameter. On the other hand, the abdomen slice phantom covers a region of 40×40 cm 2 and the distance from the X-ray source to the system origin is 80cm. For each projection view, 732 detector elements span a FOV of 50cm in diameter. Moreover, all the phantoms are simulated by  photon number counting method that the transmission data can be calculated with y i = b i e −ĝ i = e −ĝ i by setting the entrance x-ray intensity b i as 6.0 × 10 5 photons. We test the regularization parameters ranging from 0 to 1, with an interval of 0.1 to determine final parameters. For these two phantoms, the parameters of GPBB were chosen as 0.6 and 0.6; and the parameters of SIR-TV were chosen as 0.5 and 0.5; and the parameter β for the regular DL method are determined as 0.4 and 0.4; and the parameter β 2 in the SL0-DL method are chosen as 0.4 and 0.4.

2) REAL DATA
The proposed SL0-DL algorithm and other competing algorithms also applied real rat datasets from our laboratory. The X-ray tube voltage and tube current are set to 50 kV and 1 mA, respectively. The projection data are obtained by fan-beam scanning. The distance between the detector and the system origin is 160.5mm, while the distance from the X-ray source to the system origin is set to 241.8mm. We test two sampling rate situations with a total of 360 and 180 projections acquired over 360 • . The number of radial bins per view is 512, and the size of each bin is 0.15 × 0.15 mm 2 . The reconstructed image size is 512 × 512 with the real size of 44 × 44 mm 2 . For the two sampling rate situations, the regularization parameters of GPBB are chosen as 0.4 and 0.6; and the regularization parameters of SIR-TV were chosen as 0.5 and 0.5; and the parameter β for the regular DL method are determined as 0.4 and 0.6; and the parameter β 2 in the SL0-DL method are chosen as 0.4 and 0.6. Furthermore, we also verified the proposed SL0-DL algorithm via real sheep head datasets from our partner hospital. The datasets were acquired using Varian IMG-CBCT device under the full fan mode. The X-ray tube voltage and tube current are set to 100 kV and 20 mA, respectively. The distance from the X-ray source to rotation center is 1000mm, while the distance between the X-ray source and detector is 1500mm. The flat detector resolution is 768 × 1024 with each bin size 0.388 × 0.388 mm 2 . The reconstruction matrix size is 512 × 512 with the real size 250 × 250 mm 2 . We test two sampling scenarios using a total of 328 and 164 projections acquired over 360 • . For both sampling scenarios, the regularization parameters of GPBB are chosen as 0.4 and 0.4; and the regularization parameters of SIR-TV were chosen as 0.5 and 0.5; and the parameter β for the regular DL method are determined as 0.5 and 0.5; and the parameter β 2 in the SL0-DL method are chosen as 0.5 and 0.5.

3) PARAMETER SELECTION IN THE SL0-DL ALGORITHM
Some details about the parameter selection of the DL method are explained as follows. The image patches are overlapping with a sliding distance of one pixel, and the patches are with N 0 × N 0 pixels (N 0 = 7). The reconstructed image is of N × N pixels (N = 512) so the number of the patches is S = (N − N 0 + 1) 2 = 256036. The dictionary, which is pre-trained by K-SVD algorithm, consist of K = 196 atoms and for each study the regular DL and SL0-DL method use the same dictionary. The fixed sparsity of over-complete dictionary is set as L S 0 = 4. The initial value of the relaxation parameter β 1 can be empirically determined as 1. However, the initial values of the other two parameters (β in the regular DL method and β 2 in the SL0-DL method) should be chosen properly by tests to get a good result.
All the algorithms of the above experiments are coded in MATLAB and run on a dual-core PC with 3.10 GHz Intel Core i5-2400.

B. RESULTS OF PHANTOM STUDIES
Under the premise that a reference image is used as the gold standard, the quality of image reconstruction is quantitatively evaluated by the peak-signal-to-noise ratio (PSNR), the normalized mean absolute deviation (NMAD) and the structure similarity index (SSIM), PSNR and NMAD are  defined by (16) and (17), SSIM can be found in literature [30].
where µ truth is the gray value of reference image, L is the maximum gray level of the image. The reference images are shown in Fig. 2(a) and Fig. 3(a). Fig. 2 (b)-(f) and Fig. 3 (b)-(f) show the reconstruction results of FBP, GPBB, SIR-TV, DL and SL0-DL with the two phantom studies. The PSNR, NMAD and SSIM criterions to evaluate the image quality are presented in Table 1. In addition, to better verify the artifact suppressing effect of the proposed algorithm, the enlarged regions of the human head are shown in Fig. 2 and Fig. 4 respectively. The image reconstructed by FBP method shown in Fig. 2 is ruined by the low sampling rate and noise. The results in Table 1 indicate that SL0-DL shows some improvements over regular DL, and second DL based methods are better than GPBB and SIR-TV by comparing PSNR, NMAD and SSIM.
Looking critically at the enlarged images in Fig. 2 (a)-(f), the resolution of SL0-DL and DL is significantly higher than that of GPBB and SIR-TV. But the images reconstructed from the traditional DL method contain some irregular artifacts in the soft tissue region while SL0-DL can be able to suppress the artifacts effectively. Moreover, the DL algorithm still has obvious noise and unexpected structures, as shown by the white arrows in Fig. 2(e), but the SL0-DL algorithm has more real structures and almost none noise. When it comes to GPBB and SIR-TV methods, they are all based on L1 constrained reconstruction algorithms. Therefore, as shown in the reconstruction results and zoomed region of interest (ROIs) of Fig. 2 (c) -(d), although there are no artifact in the reconstructed result, the reconstruction image edge regions are also over-smoothed so that the PSNR is lower than the regular DL and the proposed methods. Likewise, the results with the abdomen slice experiment are shown in Fig. 3, especially the part indicated by the red arrow, the regular DL has some residual artifacts, GPBB, SIR-TV and the proposed method show some substantial improvements in artifacts reduction. However the regular DL and the SL0-DL algorithm are better at protecting the texture structure of images, while GPBB and SIR-TV destroy many CT image details.
In order to better demonstrate the advantages of the proposed method over the traditional DL algorithm in the entire reconstruction process, we compare the results of the traditional DL algorithm and the proposed method after 10, 30, 50, and 80 iterations, respectively. As shown in Fig. 4, the traditional DL method has obvious streaking artifacts and bias in early iteration results, as indicated in Fig. 4(a). As mentioned above, the streaking artifacts and bias can also be represented well by the over-complete dictionary, resulting to the reservation of these unexpected structures in the output image, as seen in Fig. 4(b)-(d). So the traditional DL method still has some residual artifacts around the soft tissue region in the final image, as shown in Fig. 2(e). In comparison, in Fig. 4(e), the proposed method was able to suppress the artifact effectively in the early iterations. In Fig. 4(e)-(h), we can observe the proposed method shows a progressive improvement of image reconstruction quality after 10, 30, 50, and 80 iterations. This illustrates that the smoothing of the intermediate image by the smoothed L 0 norm-constraint can effectively suppress artifacts. Furthermore, we also make a comparison of image details by calculating the differences between the reconstructed images and the reference image. As mentioned above, there are irregular artifacts in the reconstruction result of the traditional DL method (as shown in Fig. 2(e)), so there will be obvious irregular traces in the difference image as shown in Fig. 5 (a1) and (a2). As shown in Fig. 5 (b1) and (b2), the proposed method shows less artifact traces.
On the other hand, this fact can also be established in Fig. 6, which shows the horizontal intensity profiles going through the center of the reconstructed images. As it can be seen, in Fig. 6(a) and (b), the closer the pixel values of the reconstructed image are to the reference image, the more accurate the details of the reconstruction result. The result of SL0-DL is closest to the ideal phantom, and there is less noise in the reconstruction result. The traditional DL method is shown as a turquoise line in Fig. 6, it can be seen that its gray value will jump in stages relative to the ground truth (GT), and this can also reflect the irregular difference image traces caused by the irregular artifacts in the traditional DL method  in Fig. 5 (a1) and (b1). Thus, in conclusion, considering artifacts reduction and image texture protection, the proposed method performs better.
As the iterative process of the proposed algorithm only adds an extra step compared to the traditional DL algorithm, which is the smoothing intermediate image by SL0 constraint, so the reconstruction time of the proposed algorithm is in the same level. The average reconstruction time of one iteration: the traditional DL method is 8.96 seconds and the proposed method is 9.85 seconds. Furthermore, comparative experiments on the statistical performances between the proposed method and traditional DL are also performed, Fig.7 shows the NMAD curves of the reconstruction results of these two algorithms with different noise levels. By the comparison of the NMAD curves of head and abdomen, it can be found that the proposed method performs better than the traditional DL method.  In addition, Fig.8 shows the NMAD curves of the human head slice phantom with different noise levels. It can be found that results produced by 1 × 10 5 photons and 2 × 10 5 photons are poor, 4 × 10 5 photons is better, and 6 × 10 5 photons is enough to produce very good results, then adding more photons may not show obvious improvement of reconstruction quality.

C. RESULTS OF REAL DATA STUDY 1) MOUSE DATA
The reconstruction results of the real mouse abdomen slice are shown in Fig. 9. Fig. 9 (a1)-(a2) as the same reference images are reconstructed by FBP using 360 projections data. Fig. 9 (b1)-(e1) show the reconstruction results obtained using the GPBB, SIR-TV, DL, and SL0-DL algorithms under 360 projection angle data. Fig. 9 (b2)-(e2) show the reconstruction results obtained using the GPBB, SIR-TV, DL, and SL0-DL algorithms under 180 projection angle data. As it can be seen, similar results compared with the simulation experiments also emerge in the real animal study. As is shown Fig. 9(d1) and (d2), the results of traditional DL method are ruined by the streaking artifacts and bias. In Fig. 9(b1, b2) and Fig.9 (c1, c2), GPBB and SIR-TV have over-smoothing effect on the edge of image structure.
In order to see the reconstruction details more obviously, Fig. 10 shows the zoomed region of interest (ROI) which is marked with yellow square frame in Fig. 9(a1). In Fig. 10, GPBB and SIR-TV are able to suppress the artifact effectively, but the real image edge regions are also over-smoothed so much so that there are stair-step artifacts around the soft tissue edges, which are shown by the red arrows in Fig. 10 (b1, b2) and Fig. 10(c1, c2). In comparison, dictionary learning based methods preserve the image texture very well, which means the tissue edge regions reconstructed by the dictionary learning based methods are closer to the actual tissue distribution and look more realistic, as seen in Fig. 10 (d1), (d2), (e1) and (e2). However, as mentioned VOLUME 8, 2020 above, since the dictionary is over-completed, the artifacts and bias generated by DL in early iteration results can also be represented well by the dictionary, resulting in the reservation of these unexpected structures in the final image. As shown in Fig. 9 and Fig.10, the images reconstructed from the traditional DL method contain obvious artifacts and bias in the soft tissue region, especially as indicated the domain by the red arrows in Fig. 9 (d1) and (d2). In Fig. 9 (d1) and (d2), looking at the results of the traditional DL method under different sampling rates, we observed that the streaking artifacts become thicker and more obvious when the sampling rate declines. The proposed method shows some substantial improvements in artifacts reduction, and the results are more natural and realistic, especially in the soft tissue region. By observing the whole image reconstruction quality, considering the improvements of artifacts suppression and image texture protection of the proposed method when compared with GPBB, SIR-TV and DL in his paper, the proposed method performs significantly better than the traditional DL and GPBB methods, although the proposed algorithm still has little residual artifacts, which are shown by the white arrows in Fig. 9 (e2).
2) SHEEP HEAD DATA Fig. 11 (a1) and (a2) show the same reference images reconstructed by FBP under 656 projections data. Fig. 11 (b1)-(e1) and (b2)-(e2) are the reconstructed images with 328 projections data and 164 projections data, respectively. The lower right corner of each image is a partially enlarged image marked with a red frame. First, from the whole reconstruction results in Fig. 11, GPBB and SIR-TV have more blurred or stair-stepping artifacts in the tissue edge area, it can be observed more clearly in its corresponding ROI area. In Fig. 11 (d1) and (d2), the results of traditional DL algorithm are densely distributed with irregular stripe artifacts, which are very obvious under 164 projections data. However, as shown in Fig. 11 (e1) and (e2), both blur artifacts and irregular artifacts are well suppressed in the results obtained by the proposed algorithm, and the reconstruction results look relatively clear and natural. Then, we can further observe the ROIs of the reconstructed images. There are different levels of irregular artifacts in the soft tissue regions of traditional DL images (as shown by the red arrows in Fig. 11 (d2)), and the proposed algorithm has suppressed the artifacts well. GPBB and SIR-TV show blocky artifacts in some regions near bone tissues, however, the proposed SL0-DL perform better in the preservation of the bone tissue structures, as shown by the yellow arrows in Fig. 11 (b1)-(e1). So no matter the tissue or structure, the sharpness of the proposed method is significantly better. Moreover, as shown in Fig. 11, we can also see that the artifacts of GPBB, SIR-TV and DL reconstruction results will become very serious with the decrease of sampling rate, but the SL0-DL algorithm performs better than other methods in this respect.

IV. DISCUSSION
Two possible methods to reduce the x-ray dose are to lower mAs in CT data acquisition protocols and to decrease the number of projections. The iterative algorithm with DL method demonstrated its potency in CT image reconstruction with under-sampled projection data. However, when it comes across noisy data with low sampling rate, streaking artifacts and bias will appear in early iteration results, which hampers recovery of high-quality image. In this work, we propose the SL0-DL algorithm to solve the problem by introducing a smoothed L0-constraint DL regularization to the objective function. During each iteration process, the intermediate image generated by the DL representation is smoothed by the smoothed L0-constraint before using it to update the output of this iteration. The result indicated that the image reconstructed by SL0-DL is superior to those reconstructed by other competitive algorithms, and the smoothed L0-constraint DL regularization is useful in reducing the artifacts and bias generated by the under-sampled noisy data.
By further analyzing the experimental results, we got more information about these algorithms. It is obvious that FBP performs the worst under the situation when the data is undersampled; this fact explains that the iterative algorithms play a vital role in the under-sampled reconstruction problem. In addition, since both models of regular DL and SL0-DL have similar configuration and the only difference is the added intermediate image smoothing step, the better results achieved by SL0-DL demonstrated that some changes in the regularization term and iterative process definitely improve the image quality.
However, the convergence rate of the proposed algorithm is slightly slower than traditional DL, due to the fact that the added intermediate image smoothing step also performs some smoothing effect on the edge structures of the image. Our future work will focus on speeding up the converging rate without decreasing the image quality.

V. CONCLUSIONS
In this work, an under-sampled CT reconstruction method based on a smoothed L0-constraint DL regularization has been studied. It was observed that the result of reconstructed image by the traditional DL method contains some streaking artifacts and bias when the data are noisy with low sampling rate. Since the dictionary is over-complete, the artifacts and bias can also be represented well by the dictionary. To break the reservation of these unexpected structures, the developed algorithm suppresses the artifact by smoothing the intermediate image generated from the DL representation with the smoothed L0-constraint. By involving the smoothed L0-constraint DL regularization, the proposed algorithm performs better than the regular DL algorithm. Experimental results show that the images reconstructed by SL0-DL are better than DL, GPBB and SIR-TV. The results proved that the proposed algorithm reconstructs the image closer to the ground truth compared with other algorithms.