An improved statistical iterative algorithm for sparse-view and limited-angle CT image reconstruction

Because radiation is harmful to patients, it is important to reduce X-ray exposure in the clinic. For CT, reconstructions from sparse views or limited angle tomography are being used more frequently for low dose imaging. However, insufficient sampling data causes severe streak artifacts in images reconstructed using conventional methods. To solve this issue, various methods have recently been developed. In this paper, we improve a statistical iterative algorithm based on the minimization of the image total variation (TV) for sparse or limited projection views during CT image reconstruction. Considering the statistical nature of the projection data, the TV is performed under a penalized weighted least-squares (PWLS-TV) criterion. During implementation of the proposed method, the image reconstructed using the filtered back-projection (FBP) method is used as the initial value of the first iteration. Next, the feature refinement (FR) step is performed after each PWLS-TV iteration to extract the fine features lost in the TV minimization, which we refer to as ‘PWLS-TV-FR’.

image reconstruction, and the regularization term shows the prior information for the desired image intended to regularize the solution 14 . The traditional regularization term usually produces unfavorable smoothing effects on edge regions. To address this issue, several edge-preserving regularization terms have been presented, such as total variation (TV) regularization with a piecewise constant assumption (PCA) [15][16][17] . Numerous investigations have shown that high-quality CT images can be reconstructed via TV minimization from sparse-view sampling data. However, PCA often produces noticeable patchy artifacts in reconstructed images 5 .
Recently, we developed a statistical interior tomography approach for CT reconstruction 18 . Its potential usefulness have been demonstrated for feature preservation of interior tomography with truncated projection measurements. However, this algorithm is a little bit complex, involving FBP with projection extension, PWLS-TV and feature refinement and the last component is not very well characterized with the CT interior tomography application. Moreover, the algorithm was applied to a scanning configuration which is not practically used (namely the interior problem).
In this paper, we extend this method to a more practical protocol, sparse-view and limited-angle image reconstruction. Furthermore, the FBP algorithm is adopted instead of the FBP-EP algorithm which is time-consuming. To remove the patchy artifacts from the TV-based methods, we improved a statistical iterative image reconstruction algorithm based on minimizing the image TV that is specifically performed using PWLS criteria (PWLS-TV), which models the statistical properties of the projection data. Furthermore, after each PWLS-TV iteration, a feature refinement step (FR) is performed to restore the desired structure details that are lost in the TV minimization. For simplicity, the presented method is termed "PWLS-TV-FR". The proposed PWLS-TV-FR method tends to relieve the PCA requirement by statistical modeling and iterative feature refinement.
The remainder of the paper is organized as follows. Section II briefly reviews the CT imaging model and describes the TV-regularized PWLS criterion in the image domain and the FR step. The experimental setup and evaluation metrics are also presented in this section. In section III and IV, the experimental results on a digital phantom and using micro-CT are reported, followed by a conclusion section.

PWLS-TV minimization.
Based on the statistical properties of the projection data, the iterative PWLS approach including a data-fidelity term and a penalty term has been applied to the image domain for dual-energy X-ray CT reconstruction [19][20][21] . But the PWLS method for low-dose X-ray CT reconstruction are flawed in the suppression of noise and artifacts. According to our previous study, the TV-regularized PWLS algorithm we referred to as "PWLS-TV" can address this issue, which can be solved using the steepest descent algorithm. Based on the image gradient magnitude sparseness 22 , TV of the image is one of the most commonly used techniques: in which s and t are the index of the location of the attenuation coefficients of the desired image and α is a small constant used for keeping the equation differentiable with respect to image intensity, which is assigned a value of − 10 8 in this work. Thus, given the current estimated image μ n , a new image μ n+1 can be calculated by minimizing the cost function of PWLS-TV, which can be described as follows: n n n T n n 1 1 in which ∇TV(μ n ) denotes the gradient of TV(μ n ), τ n represents the gradient step size and β is a smoothing parameter balancing the degree of agreement between the estimated and the measured data. This can be computed as follows: n T T T n 1 1 Feature refinement (FR). High-diagnostic CT images reconstructed from limited-angle or sparse-view acquisitions are very important in clinical practice. However, due to the piecewise constant assumption, some details in the reconstruction using the TV minimization are likely to be lost; in particular, the residual image between two successive iterations in PWLS-TV minimization might contain useful structural information, especially when the dose is very low or the projection data are insufficient. Motivated by this discovery, an FR was performed after each PWLS-TV minimization iteration to extract the main features of the loss from the residual image. This is described in the following equation: in which the symbol ⊗ denotes point-wise multiplication, P is a non-negative control notation described as 0 . + f t n 1 denotes a feature descriptor which is defined as: in which the constants C is introduced for numerical stability. Local statistics σ p , σ q and σ pq at pixel i are defined as follows: , p i and q i denote two local image patches with a size of × N N centered at pixel i and extracted from µ + n 1 and a degraded image µ d obtained by a linear Gaussian filter from µ + n 1 , respectively. The nature of the presented model structure descriptor is a structure-related map and includes a contrast variation component and a structure correlation component. The former calculates the reduction of contrast variation due to the degraded operation, and the latter quantifies the structural correction between the original and degraded images 23 . The value of each element in the feature descriptor image f t is in the interval [0, 1], and the larger the value, the more likely it belongs to part of the structure. In this paper, the proposed statistical iterative CT image reconstruction algorithm using the FR approach is referred to as "PWLS-TV-FR" and includes the three main steps shown in Fig. 1. Data availability. The datasets generated or analyzed during the current study are available from the corresponding author on reasonable request.

Experiments
To validate the performance of PWLS-TV-FR, we performed a digital simulation experiment and a physical experiment. We used digital phantom simulation data that spanned 120 views to verify the effect of PWLS-TV-FR on sparse-view data. Additionally, sparse-view anesthetized mouse data and limited-angle mouse bone sample data were acquired from a commercial CT scanner to demonstrate the feasibility of the proposed method in practical applications for sparse-view or limited-angle sampling data. For both phantom simulation and Micro-CT data, Table 1 shows the idea projection number needed for FBP and the projection number used for proposed method.

Digital phantom simulations.
A head model of the digital XCAT phantom was chosen as experimental simulation data (as shown in Fig. 2) to evaluate the performance of the FR step after each PWLS-TV iteration. To simulate sparse sampling, only 120 views were spanned on a circular orbit of 360°, and the number of bins per view was 672. The distance from the X-ray source to the detector was 1040 mm, and the distance from the rotation center to the curved detector was 570 mm.
Micro-CT data. In this experiment, two scans were performed using a Micro-CT Skyscan1076 scanner with a normal-dose (70 kVp, 225 mAs) protocol. One scan was performed on an anesthetized mouse, in which 450 projections were uniformly collected over a 360° range, and another scan was performed on a mouse bone sample, in which 291 projections were uniformly collected over a 197.88° range, which can be considered limited-angle sampling. In this work, 450 projections of the anesthetized mouse were encoded from 0 to 449, and only even-numbered projections were selected for use as under-sampled data for a final total of 225 projections. Only the central slice projection was extracted, which is typical for fan-beam geometry. For each projection, 4000  detector elements were equally angularly distributed with a total length of 50.12 mm. The distance from the object to the curved detector was 165 mm, and the distance from the object to the X-ray source was 121 mm.
Performance evaluation. (1) In this work, the root mean squared error (RMSE) was used to evaluate reconstruction accuracy and is described as: m Q true 1 2 in which µ denotes the reconstructed image, µ true denotes the true image, m indexes the pixels of the image and Q is the total number of pixels in the image.
(2) To evaluate the noise reduction capability of different reconstruction algorithms, the peak signal to noise ratio (PSNR) metric was adopted in this paper and is as follows: in which µ denotes the reconstructed image, µ true denotes the true image, and µ MAX( ) true represents the maximum intensity value of µ true .

Results
Digital phantom result. In this experiment, with sparse sampling data, images were reconstructed using three methods, FBP, PWLS-TV and PWLS-TV-FR, which are shown with partial magnification in Fig. 3. From Fig. 3, it is clear that there are severe streak artifacts in the image reconstructed using FBP, whereas the PWLS-TV results exhibit fewer artifacts but degraded details. Moreover, the image reconstructed using PWLS-TV-FR showed noticeable improvements over that reconstructed using PWLS-TV by preserving small-scale and low-contrast structures, as indicated by the arrows in the amplified areas. Figure 4 shows residual reconstructed images (shown in Fig. 3) resulting from the true image subtracted from the reconstructed image. The smaller the features of the residual image, the more accurate the reconstructed image. As is clear from Fig. 4, the residual reconstructed image using PWLS-TV-FR has the smallest features, which indicates that the PWLS-TV-FR method for sparse-view data retains more details than either the PWLS-TV or FBP methods.
In the process of updating images using the PWLS-TV or PWLS-TV-FR functions, the trend is that the image gradually becomes clearer and eventually reaches an extreme value as the number of iterations increases. However, as shown Fig. 5, whether it is evaluated by RMSE or by PSNR, the PWLS-TV-FR method always performs better than PWLS-TV. Figure 6 depicts the horizontal profiles of the reconstructions in Fig. 3 and the corresponding horizontal residual value. The result reconstructed using PWLS-TV-FR was closer to the reference profiles, and the residual values reconstructed using PWLS-TV-FR are much smaller than those using PWLS-TV;    it is obvious that PWLS-TV-FR achieves better performance. Figure 7 shows a histogram of the RMSE and PSNR evaluation results for the reconstructions in Fig. 3. The red column represents the RMSE value, and the blue column represents the PSNR value. It is important that the RMSE value of PWLS-TV-FR is the minimum and the PSNR value of PWLS-TV-FR is the maximum, which verifies the validity of the PWLS-TV-FR algorithm for sparse-view digital phantom data. Figures 8 and 9 show images and zoomed images reconstructed using FBP, PWLS-TV and PWLS-TV-FR. Figure 8 is a reconstruction using sparse-view sampling of the anesthetized mouse projection data with 225 views, which shows that the PWLS-TV-FR method can suppress severe artifacts and preserve fine tissue structure from sparse-view sampling data in the clinic. Figure 9 is the reconstruction using limited-angle sampling data of mouse bone sample projection data using a scanning angle of 197.88°. Image reconstructions using FBP showed severe streak artifacts resulting in low detail visibility. The images reconstructed using the other two algorithms show a lower level of streak artifacts. However, the features reconstructed using PWLS-TV-FR are clearer than those reconstructed by PWLS-TV. In conclusion, the experimental micro-CT data demonstrate that the presented method can recover structural information lost using the PWLS-TV method for sparse-view or limited-angle sampling data.

Conclusion and Discussion
To address the global use of high-diagnostic CT images obtained from sparse-view or limited-angle projection data, we improved on the statistical iterative algorithm by including an FR step. Two key components were included in this algorithm. One is a TV minimization model using PWLS criteria that account for the statistical nature of projections. The other is an FR step that is performed after each iteration of PWLS-TV to recover any  missing fine features. The improved method is designed to relieve the piecewise constant assumption of TV-based methods for sparse-view or limited-angle CT image reconstruction. The reconstruction results from a digital phantom, anesthetized mouse and a mouse bone sample dataset show that the proposed PWLS-TV-FR approach is a significant improvement over the PWLS-TV method in terms of the preservation of structural details and the suppression of undesired patchy artifacts. This result demonstrates its potential to reduce CT radiation doses by decreasing the number of projection views or the rotation angle of the scan.
An iterative reconstruction algorithm usually involves numerous parameters, which can significantly impact reconstruction properties. In this work, we try our best to select the "optimal" parameters for different algorithms in order to make a fair comparison on their performance. For example, the smoothing parameter β balancing the degree of agreement between the estimated and the measured data in function (2) should be a small value to avoid data distortion, however, after testing, we found that β = − e 2 is better. Furthermore, the image patch size in the feature descriptor designed to distinguish structure from noise and artifacts has to be carefully selected because the large one leads to better detail-detection ability, while the small one can benefit the computational load. According our previous study, the patch size of × 3 3 is a good choice for balancing the structure detection ability and computational efficiency. Additionally, the parameter C is introduced for numerical stability to avoid the condition that σ σ + p q 2 2 is very close to zero. Therefore, it is assigned . × − 1 25 10 5 as a small constant in our work.