A feature refinement approach for statistical interior CT reconstruction

Interior tomography is clinically desired to reduce the radiation dose rendered to patients. In this work, a new statistical interior tomography approach for computed tomography is proposed. The developed design focuses on taking into account the statistical nature of local projection data and recovering fine structures which are lost in the conventional total-variation (TV)—minimization reconstruction. The proposed method falls within the compressed sensing framework of TV minimization, which only assumes that the interior ROI is piecewise constant or polynomial and does not need any additional prior knowledge. To integrate the statistical distribution property of projection data, the objective function is built under the criteria of penalized weighed least-square (PWLS-TV). In the implementation of the proposed method, the interior projection extrapolation based FBP reconstruction is first used as the initial guess to mitigate truncation artifacts and also provide an extended field-of-view. Moreover, an interior feature refinement step, as an important processing operation is performed after each iteration of PWLS-TV to recover the desired structure information which is lost during the TV minimization. Here, a feature descriptor is specifically designed and employed to distinguish structure from noise and noise-like artifacts. A modified steepest descent algorithm is adopted to minimize the associated objective function. The proposed method is applied to both digital phantom and in vivo Micro-CT datasets, and compared to FBP, ART-TV and PWLS-TV. The reconstruction results demonstrate that the proposed method performs better than other conventional methods in suppressing noise, reducing truncated and streak artifacts, and preserving features. The proposed approach demonstrates its potential usefulness for feature preservation of interior tomography under truncated projection measurements.

the desired structure information which is lost during the TV minimization. Here, a feature descriptor is specifically designed and employed to distinguish structure from noise and noise-like artifacts. A modified steepest descent algorithm is adopted to minimize the associated objective function. The proposed method is applied to both digital phantom and in vivo Micro-CT datasets, and compared to FBP, ART-TV and PWLS-TV. The reconstruction results demonstrate that the proposed method performs better than other conventional methods in suppressing noise, reducing truncated and streak artifacts, and preserving features. The proposed approach demonstrates its potential usefulness for feature preservation of interior tomography under truncated projection measurements.
Keywords: compressed sensing, computed tomography (CT), interior tomography, penalized weighed least-square (PWLS), feature refinement (Some figures may appear in colour only in the online journal)

Introduction
Excessive radiation exposure in computed tomography (CT) examinations can induce cancerous, genetic and other diseases, which has caused significant concerns to patients (Cohnen et al 2006, Brenner and Hall 2007, McCollough et al 2009, Hu et al 2011, Zou et al 2012, Long and Fessler 2014. In many clinical practices, physicians only pay attention to relatively small interior organs, such as heart, lung and bones. The radiation dose would be reduced significantly if only the interior region of interest (ROI) are scanned using a narrower beam for which x-rays transverse through the ROI (Wang and Yu 2013). Motivated by this need, exact reconstruction of the interior ROI from truncated local projections has been widely studied, and is commonly referred to as the 'interior problem' (Sidky et al 2006, Snyder et al 2006, Sidky and Pan 2008, Han et al 2009, Zeng and Gullberg 2009a, 2009b, Bian et al 2010, Yang et al 2010, Lauzier et al 2012, Yang et al 2012a, Wang et al 2014, Rose et al 2015, Han et al 2015. Generally, due to the inconsistent data measurements for CT reconstruction, the interior problem is ill-posed and does not have a unique solution. Therefore, applying a traditional CT algorithm on interior problem will lead to degraded reconstruction accuracy (Wang and Yu 2013). To address this issue, different additional information is introduced to achieve an exact interior reconstruction. At present, there are three types of such additional information: (1) additional projection data beyond the local truncated parts (Louis and Rieder 1989, Wiegert et al 2005, Snyder et al 2006, Lauzier et al 2012; (2) complete knowledge on a small sub-region inside the ROI is available , Zeng and Gullberg 2009a, 2009b; and (3) the interior ROI is piecewise constant or polynomial (Han et al 2009, Yang et al 2010. Some researchers have used previous high-quality result to induce the current reconstruction with the low-quality projection (Lauzier et al 2012). Such methods were mostly applied to dynamic cone-beam CT (CBCT) imaging with a small sized detector. For example, in image-guided radiation therapy, the high-quality planning CT measurements were used as the prior to reconstruct the target image in CBCT imaging wherein the projection data is incomplete and/or truncated. The performance of this kind of approach heavily depends upon the accuracy of image registration operations between the target-and prior-projection measurements. Meanwhile, designing a dedicated registration algorithm is a very challenging task in practice (Lauzier et al 2012). On the other hand, projection extrapolation method as an alternative method has been presented to fit the missing projection segments (Louis and Rieder 1989, Wiegert et al 2005, Snyder et al 2006. Specifically, the support of the object is first contained in a defined ellipsoid region, and then the limited-size detector is extended to a 'virtual'one. As the reconstruction volume is forward projected, some calculated raysums will be mapped to the extended detector region. Finally, these missing rays on the detector are then given the values of the calculated raysums at the current iteration. However, this extrapolation method cannot yield exact results with theoretical basis and often adds new artifacts in the boundaries. Different from the above mentioned methods, iterative image reconstruction algorithms are capable of handling projection truncation naturally and producing fewer truncation artifacts. And, the interior problem can be exactly and robustly solved based on the assumption that a sub-region in the ROI is known in advance , Zeng and Gullberg 2009a, 2009b. For example, a general interior reconstruction approach was presented by using a truly truncated Hilbert transform (THT) on a line-segment inside a compactly supported object, wherein partial knowledge on one or both neighboring intervals of that segment is used (Noo et al 2004, Defrise et al 2006, Ye et al 2007. To solve the inversion of the THT, the projection onto convex set (POCS) method , Zeng and Gullberg 2009a, 2009b and the singular value decomposition (SVD) method  were adapted respectively. Compared with the SVD method, the POCS method is more computationally expensive, but easier to employ additional constraints. Based on this foundation, an interior reconstruction approach with a truncated limited-angle scan was further extended . However, in many practical situations, especially in the scenario of contrast-enhanced CT studies, prior knowledge on a sub region inside the interior ROI is unavailable in advance, which limits its applications in clinical practice.
Recently, the image reconstruction community experienced a dramatic revitalization when the compressed sensing (CS) theory was developed, which addressed the sampling condition for accurate reconstruction of an image object under proper assumption (Candes et al 2006, Donoho 2006, Chen et al 2008, Lubner et al 2011, 2016, Niu et al 2014b, 2014c, Xi et al 2015. Inspired by this powerful theory, it was mathematically proven that an interior ROI can be exactly and robustly reconstructed via a total variation (TV) minimization (Han et al 2009 or high-order TV (HOT) minimization (Yang et al 2010(Yang et al , 2012a(Yang et al , 2012b, if the ROI is piecewise constant or piecewise polynomial. The promising part is that this approach doesn't need any previous knowledge about projection data or small sub-region inside the ROI, in addition to the piecewise constant or polynomial assumption. However, this assumption may not be valid globally within the interior ROI, especially in some regions where subtle image intensity changes. Therefore, TV or HOT minimization may produce images with low-frequency patch structures, and lose fine structure information (feature) which conventionally is the focus in interior tomography. This phenomenon becomes even more severe especially when the dose is very low (Tang et al 2009, Lauzier and Chen 2013, Niu et al 2014a. Moreover, these methods did not take into account the statistical nature of the local projection data and would not work well in the case of low-dose data Tsui 2011, Xu et al 2011). Motivated by these observations, restoring the fine structures lost in TV minimization and utilizing statistical model are highly desirable in order to achieve exact interior reconstruction.
In this paper, we propose a statistical interior tomography approach to improve interior tomography performance for radiation dose reduction while maintaining a high-quality image reconstruction. Specifically, a TV based interior tomography is performed under the criteria of penalized weighed least-square (PWLS-TV) which takes into account the statistical nature of the local projection data. More importantly, a feature refinement step is conducted after each iteration of PWLS-TV to recover the desired structure information which is lost in the TV minimization. The proposed approach aims to relieve the requirement of piecewise constant or polynomial by statistical modeling and iterative feature refinement.
The remaining sections of the paper are organized as follows. Section 2 describes the CT imaging model, TV regularized PWLS criterion for CT image reconstruction, wherein an interior projection extrapolation based FBP reconstruction and the interior feature refinement step are introduced in detail. Moreover, the experimental setup and evaluation metrics are also outlined in this section. In section 3, the experimental results on digital phantom and in vivo Micro-CT are reported, followed by sections discussion and conclusion, respectively.

CT imaging model
Without loss of generality, under the assumption of a mono-energetic beam, the x-ray CT measurement can be approximately expressed as a discrete linear system: where y represents the obtained sinogram data (projections after system calibration and logarithm transformation), i.e. = y y y y , ,..., M T 1 2 ( ), µ is the vector of attenuation coefficients to be estimated, i.e. , where 'T' denotes the matrix transpose. The operator H represents the system or projection matrix with the size of × M N. The element of H denotes the length of intersection of projection ray i with pixel j. In our implementation, the associated element was pre-calculated by a fast ray-tracing technique stored as a file. The goal for CT image reconstruction is to estimate the attenuation coefficients µ from the measurement y with H.
In this study, the proposed CS based statistical interior tomography approach contains three steps. First, an interior projection extrapolation based FBP reconstruction is conducted. Second, a TV based interior tomography is performed under the PWLS criteria with the FBP reconstruction as the initial guess. Third, a feature refinement step is added to recover the losing structure information and preserve details. The above three steps will be introduced in the following sub-sections, respectively.

PWLS criterion for CT image reconstruction
Statistical iterative reconstruction (SIR) for x-ray CT, utilizing the statistical distribution resulting from the x-ray interaction process, has been extensively explored for radiation dose reduction in CT field (De Man et al 2001, Elbakri and Fessler 2002, Tang et al 2009, Ma et al 2012b, Huang et al 2013, Lauzier and Chen 2013, Niu et al 2014a. Among them, PWLS is one of state-of-the-art approaches. Mathematically, the PWLS criterion for CT image reconstruction can be rewritten as follows (Sukovic and Clinthorne 2000): where H represents the system or projection matrix, y represents the obtained sinogram data.
is the penalty term. β is a hyper-parameter to balance the fidelity term and penalty term. Σ is a diagonal matrix with the ith element of σ i 2 which is the variance of sinogram data y. In the implementation, the variance of σ i 2 was determined by the following mean-variance relationship (Wang et al 2008a, Ma et al 2012a: where I 0 denotes the incident x-ray intensity, p i is the mean of the sinogram data at bin i and σ e 2 is the background electronic noise variance. According to the CS-based interior tomography theory, the interior ROI can be exactly reconstructed via minimizing an objective function which is associated with an appropriate sparsifying transform (Han et al 2009, Yang et al 2010. Among all the existing sparsifying transforms, total variance (TV) is one of most commonly used: here s and t are the indices of the location of the attenuation coefficients of the desired-image. δ is a small constant used for keeping differentiable with respect to the image intensity. In this paper, the TV regularized PWLS statistical formulation is referred to as 'PWLS-TV', which can be solved using the steepest descent algorithm. Two key ingredients for improving the performance of PWLS-TV for interior tomography: projection extrapolation and interior feature refinement are integrated into the PWLS-TV framework.

Projection extrapolation
In the original PWLS method, the initialization is usually with a smooth FBP image, whose low spatial-frequency components are nearly correct (Wang et al 2006). However, the FBP reconstruction from truncated projection has serious artifacts due to truncation. In this regard, the pre-processing on the truncated projection is needed to obtain a good initial estimation. In this paper, the projection data extrapolation step is performed to mitigate the truncation artifacts and to moderately extend the field of view. Specifically, this step approximates the support of the object by an ellipsoid, which is large enough to contain all object pixels, so that they can be possibly assigned to their correct locations (Wiegert et al 2005, Xu and Tsui 2011, Lauzier et al 2012. The typical procedure is as following: the projection data of an ellipsoid is simulated, and the extrapolation is performed between the boundary of the truncated projection and the projection of the simulated ellipsoid. The FBP algorithm is then used to reconstruct the extrapolated projection. This procedure ensures continuity at the boundary between the measured and extrapolated regions. In this paper, the FBP reconstruction with projection extrapolation is referred to as 'FBP-EP', and when we refer to 'PWLS-TV' in this paper, which always includes FBP-EP initialization.

Interior feature refinement (IFR)
In clinical practice, the image object is complicated. Small-scale and low-contrast structures are often very important for diagnosis. However, fine structures, such as small vessels and tissues are likely to be lost in the TV minimization reconstruction, especially when dose is very low. This is mainly due to the piecewise smooth assumption (Sidky et al 2006, Sidky and Pan 2008, Han et al 2015. It means the residual image between two successive iterations in PWLS-TV minimization may contain useful structure information, which has been discarded in the conventional PWLS-TV method. Inspired by this observation, an interior feature refinement (IFR) step is proposed to improve the performance of the TV minimization on detail preservation. This step aims to extract the structure information from the residual image and add back to the TV minimized image, which can be described in following equation: where µ IFR is the feature refined image, µ is the TV minimized image before feature refinement, ν is the residual image between two successive iterations in PWLS-TV minimization. f t denotes a feature descriptor, ν ⊗ f t is the detected feature image and the symbol ⊗ represents the point-wise multiplication. It is clear that the performance of the IFR step heavily depends on the construction of the feature descriptor. This descriptor can be seen as a filter or mask for extracting the meaningful structure information and discarding the noise and noiselike artifacts. motivated by the idea from the computer vision community , a feature descriptor f t is designed as: where local statistics σ p , σ q and σ pq at pixel i are defined as, , p i and q i denote two local image patches with size of × N N centered at pixel i in corresponding interior ROI, extracted from µ and a degraded images µ d obtained by a linear Gaussian filter from µ, respectively. The constants C 1 and / = C C 2 2 1 are introduced for numerical stability. The nature of the proposed model structure descriptor is a structure-related map, and it involves a contrast variation component and a structure correlation component. The former calculates the reduction of contrast variation due to the degraded operation, and the later quantifies the structural correction between the original and degraded images (please refer to  for details of the descriptor). The value of each element in the feature descriptor image f t is in the interval [0, 1], and the larger its value is, the more likely it belongs to part of the structure. In this paper, the proposed CS based statistical interior tomography with feature refinement approach is referred to as 'PWLS-TV-IFR'.

Implementation
Presented in figure 1 is a flowchart representation of our approach for CT interior image reconstruction. It is comprised of the following three main steps: (1) Applying Projection Extrapolation: Given the truncated projection y, projection extrapolation is performed by simulating an ellipsoid support. Then the FBP algorithm is used to obtain the initial guess µ 0 .
(2) Minimizing the Cost-Function of PWLS-TV: Given the current image estimation µ n , the steepest descent scheme is utilized to yield a new image estimation, i.e. µ P n , which can be expressed as follows: where ( ) µ ∇ TV n represents the gradient of ( ) µ TV n , η n represents the gradient step-size which can be calculated adaptively by the following estimator: (3) Performing Interior Feature Refinement: After each iteration of PWLS-TV, the reconstruction result µ P n is replaced by µ n IFR with a greater number of losing structure are extracted and recovered, which progressively improves the image quality. The feature refinement step is µ the residual image between two successive iterations in the process of the PWLS-TV minimization. f t n is the feature descriptor extracted from µ P n , and ν ⊗ f t n n is the detected structure image. Finally, update µ µ = + n n 1 IFR . Output µ + n 1 until the stopping criterion is satisfied. A modified gradient descent algorithm (O'Sullivan and Benac 2007) was adopted to solve the reconstruction problem. In the implementation, to find a stable convergent solution, the root mean square error (RMSE) metric was used to measure the reconstruction quality. In summary, the pseudo-code for the presented PWLS-TV-IFR approach can be described as follows: 1: Initialization: µ IFR 0 = FBP-EP { } y ; 2: While a stopping criterion is not met; 3: µ µ = − , n n IFR 1 4: calculate the gradient step-size η n using (8), 5: obtain a new image estimation µ P n using (7), 6: calculate the residual image ν µ µ = − n n P n , 7: calculate the feature descriptor f t n from µ P n using (6), 8: extract the structure information from the residual image and refine the image estimation µ µ ν = + ⊗ f n P n t n n IFR , n = n + 1 9: End In line 1, an initial estimate of the to-be-reconstructed image is set to be the preliminary image reconstructed by the FBP-EP method. Each loop (lines 3-8) is performed by two components: minimizing the cost-function of PWLS-TV (lines 4-5), and feature refinement step (lines 6-8).
Since the feature descriptor which is designed to distinguish structure from noise and artifacts plays an important role in our proposed framework, the selection on several scalar parameters is critical. For example, the image patch size has to be carefully selected because the large one leads to better detail-detection ability, while the small one can benefit the computational load. Practically, the patch size of 3 × 3 is a good choice for balancing the structuredetection ability and computational efficiency. Additionally, the parameter C 1 is included to avoid instability when σ σ + p q 2 2 is very close to zero. Therefore, it should be assigned a small constant (1 × 10 −5 in our work).

Experiments
A comprehensive numerical simulation was carried out to evaluate the performance of PWLS-TV-IFR. Additionally, in vivo Micro-CT datasets were acquired to demonstrate the feasibility of our proposed algorithm in practical applications. The TV based compressed sensing approach presented by  was carried out for comparison, which is referred to as 'ART-TV'. In this paper, the ART-TV approach was initialized with zero images similar to the algorithm by . In addition, the local FBP-EP approach and the TV regularized PWLS approach without an interior feature refinement step, which is referred to as 'PWLS-TV' was also carried out for comparison. For each experiment, a large range of different parameters were tried for each method and the ones with the best performance in terms of RMSE were selected for comparison. All algorithms used in this work were implemented in Matlab 7.0 (The MathWorkds, Inc.). The code was run on a PC with Intel(R) i3-2120, 3.30 GHz processor and 2 GB of memory.

Digital phantom simulations
In this experiment, a digital XCAT phantom (Segars et al 2010) was used for experimental data simulations (as shown in figure 2). We chose a geometry that is representative for a monoenergetic fan-beam CT scanner setup with a circular orbit to acquire 1160 views over 2π.
The distance from the rotation center to the curved detector is 570 mm and the distance from the x-ray source to the detector is 1140 mm. Each projection datum along an x-ray through the sectional image was calculated based on the known densities and intersection areas of the x-ray with the geometric shapes of the objects in the sectional image. An equi-spatial virtual detector was assumed, which is centered at the system origin and made perpendicular to the line from the system origin to the x-ray source. The number of channels per view is 168 with a total length of 235 mm. Each phantom is comprised of 512 × 512 square pixels with the intensity value from 0.004 to 0.019. The size of each pixel is 0.5 mm × 0.5 mm. A ROI is defined by the local scanning beam, and only the projection data through the ROI were obtained.
For the noisy truncated projection data, similar to the studies presented in (Ma et al 2012b, Huang et al 2013, Niu et al 2014a, after calculating the noise-free truncated line integral y as a direct projection operation based on the model (1), the noisy measurement b i at each bin i was generated according to the statistical model of pre-logarithm projection data: where I 0 denotes the incident x-ray intensity and σ e 2 is the background electronic noise variance, which was set to 10 in the experiment. The noisy measurements y i were calculated by the logarithm transform of b i . In the present study, the normal-and low-dose scanning protocols were simulated by setting the x-ray exposure level to 5.0 × 10 5 and 1.0 × 10 5 , respectively. In the case of fewer projection data, they were generated by under-sampling the 1160 projections of the normal-dose simulation to only 580 and 360 projections evenly over 360°.

In vivo Micro-CT data
In vivo imaging experiments were performed using a Micro-CT Skyscan1076 scanner. In this experiment, a CT scan of an anesthetized mouse was performed and only central slice projection was extracted, which is typical fan-beam geometry. The distance from the object to the curved detector is 165 mm and the distance from the object to the x-ray source is 121 mm. In our experiments, two scans were performed. One used a normal-dose (70 kVp, 225 mAs) protocol where 450 projections were uniformly collected over a 360º range, and the other used a low-dose (70 kVp, 27 mAs) protocol where 360 projections were uniformly collected over a 360° range. For each projection, 4000 detector elements were equi-angularly distributed with the total length of 50.12 mm. We first reconstructed the entire cross-section of the lungs in a 500 × 500 matrix covering a 12.5 mm × 12.5 mm region from each full-scan dataset. Then a circular region around the heart was chosen as an ROI. Finally, only the projection data through the ROI was kept to simulate an interior scan. To evaluate the performance of the proposed algorithm with fewer projections, the views were down sampled from 450 to 225 and 150 for the normal-dose protocol and from 360 to 180 for the low-dose protocol. Figure 3 is the normal-dose image reconstructed by using the FBP method with ramp filtering from the 450-views projection data and the corresponding interior ROI.

Performance evaluation
(1) Evaluation by reconstruction accuracy. To quantitatively evaluate the reconstruction accuracy, the root mean squared error (RMSE) metrics were calculated over a region of interest (ROI). The RMSE is defined as: where µ denotes the to-be-reconstructed image, µ xtrue denotes the normal-dose image or ground truth image. m indexes the pixels in the ROI and Q is the total number of pixels in the ROI.
(2) Evaluation by noise reduction. The following two metrics were utilized to evaluate the noise reduction for the comparisons among FBP-EP, ART-TV, PWLS-TV and PWLS-TV-IFR: (1) peak signal-to-noise ratio (PSNR); and (2) mean per cent absolute error (MPAE): where µ denotes the to-be-reconstructed image, µ xtrue denotes the ground truth image, and ( ) µ MAX xtrue represents the associated maximum intensity value. Q is the number of pixels in the ROI.
(3) Noise-resolution tradeoff. The noise-resolution tradeoff curves from PWLS-TV and PWLS-TV-IFR were generated from the simulated projection data using the XCAT phantom. The image resolution was analyzed by the edge spread function (ESF) along the vertical profile. Based on the strategy described in La Riviere and Billmire (2005), by assuming the broadening kernel is a Gaussian function with a standard deviation δ b , an error function (erf) can be used to represent the ESF function parameterized by δ b . Consequently, the parameter δ b can be calculated by fitting the vertical profiles to the error function, and the associated full-width at half-maximum (FWHM) of the Gaussian broadening kernel is denoted as 2.35δ b , which is used to indicate the resolution of the reconstructed image. In this study, the noise level of the reconstructed image was characterized by the standard deviation of a uniform region of size 40 × 40 in the background region. By varying the penalty parameter β, we obtained the noise-resolution tradeoff curves from the images reconstructed using PWLS-TV and PWLS-TV-IFR.  Figure 4 shows the images and zoomed ROIs reconstructed by FBP-EP, ART-TV, PWLS-TV and PWLS-TV-IFR using the simulated truncated normal dose projection data ( = × I 5 10 0 5 , 1160 views), respectively. It can be observed that the truncation artifacts are mostly corrected on the FBP-EP image. The use of ART-TV and PWLS-TV result in less artifacts but worse resolution and degraded details. The reconstruction using the proposed PWLS-TV-IFR method exhibits noticeable gains than those using ART-TV and PWLS-TV in preserving the small-scale and low-contrast structures, as indicated by the arrows in the zoomed four ROIs. Figure 5 depicts the horizontal profiles of the reconstructions in figure 4. The result reconstructed by ART-TV has an obvious bias while PWLS-TV achieves better performance. Moreover, it is noticeable that the profiles from PWLS-TV-IFR agree much better with the profiles from the reference. In other words, the presented PWLS-TV-IFR method can achieve higher low-contrast detestability as compared with other methods. Table 1 lists the RMSE, PSNR, and MPAE measures of the images reconstructed by FBP-EP, ART-TV, PWLS-TV, and PWLS-TV-IFR, as shown in figure 4. It can be seen that the results from ART-TV, PWLS-TV, and PWLS-TV-IFR exhibit an average of more than 20%  gains over that from the FBP-EP method in terms of the RMSE and MPAE, and of more than 8% gains in terms of the PSNR. The PWLS-TV-IFR method is found to perform better than both ART-TV and PWLS-TV with less noise and better reconstruction accuracy.  The appealing performance of the PWLS-TV-IFR method is mainly due to the contribution of the IFR step. Figure 6 shows the current image estimations µ n , the images µ P n after PWLS-TV minimization, the residual images ν n , the feature descriptors f t n and the detected feature images ν ⊗ f t n n (which is also the difference image before and after refinement) in the first, tenth and last iteration. It can be observed that in the first iteration, the image µ 0 exhibits many truncated artifacts and noise, most of which are filtered out in the residual image ν 0 . The feature descriptor f t n estimated fromµ P n could recognize the structure component from the artifacts and noise. As the iteration continues, the feature descriptor f t n gradually contains more and more detailed structures. Finally, more structural information is recovered and a better reconstruction is obtained. Meanwhile, the evolution of the figures along iterations can also provide an interpretation of the convergence property of our proposed framework.
To demonstrate the convergence property of our proposed method, figure 7 shows the values of the RMSE and PSNR against the number of iterations. It can be observed that, both PWLS-TV and PWLS-TV-IFR achieve convergence at the 80th iteration, since the convergence of the algorithm is guaranteed by the convexity of the objective function. Obviously, both PWLS-TV and PWLS-TV-IFR exhibit better image quality in the initial stage of the iteration. This could be due to the fact that PWLS-TV and PWLS-TV-IFR use the FBP-EP image as the initial guess, while the ART-TV method uses zero images. Meanwhile, PWLS-TV-IFR outperforms PWLS-TV in just first decades of iterations. It can be contributed to the FBP initialization which provides rich fine-feature information in initial iteration steps for IFR to work quickly, considering the nature of IFR for extracting structure information from the residual image between two successive iterations. The PWLS-TV-IFR method achieves a lower RMSE and higher PSNR after several iterations, implying that our algorithm converges to a stable and better solution. The parameter β is a hyper-parameter to balance the fidelity term and penalty term, which is determined by optimizing the trade-off between noise and image sharpness through a trialerror process. Figure 8 shows the noise-resolution curves obtained via both PWLS-TV and PWLS-TV-IFR. One vertical profile as indicated by the line in the image located at the bottomleft corner was selected. Additionally, one uniform region near the corresponding profile as indicated by the square in the background was selected for calculating the standard deviation  of the reconstructed image. When doing the noise-resolution experiment, the iteration number N was set as 100. By varying the penalty parameter β from 5.0 × 10 −3 to 5.0 × 10 −2 , we obtained the corresponding noise-resolution tradeoff curves. The noise-resolution curves from PWLS-TV and PWLS-TV-IFR have similar trends, but the PWLS-TV-IFR method shows a slightly higher resolution level in most noise range. This is probability due to the consideration of the feature refinement by the IFR step.

Reconstruction from simulated truncated low-mAs datasets.
This experiment is expected to evaluate the performance of the proposed PWLS-TV-IFR approach with simulated low-mAs and/or sparse-view projections ( = × I 1 10 0 5 , 1160, 580 and 360 views). Figures 9-11 show the reconstructed images by different reconstruction methods. It can be observed that the FBP-EP images have mitigated truncation artifacts, but are noisy and affected by undersampling. When the number of views was decreased from 1160 to 360, more streak artifacts appear in the FBP-EP images. The images reconstructed by the other three algorithms degenerated indistinctively, especially the PWLS-TV-IFR ones. In any case with different views, the FBP-EP method exhibits more noise and artifacts. ART-TV and PWLS-TV mitigate both cupping and streaks. However, the small-scale and low-contrast structures are particularly degraded. Comparing with the other three approaches, PWLS-TV-IFR always achieve the best performance in preserving small-scale and low-contrast structures. Table 2 lists the RMSE, PSNR, and MPAE measures of the images reconstructed by FBP-EP, ART-TV, PWLS-TV, and PWLS-TV-IFR with different level of views, as shown in figures 9-11, respectively. The quantitative accuracy of each algorithm increases as expected from the qualitative evaluation. For FBP-EP and ART-TV, however, on the low dose dataset, view down-sampling makes the results distinctly worse, especially in the case of 360 views. However, in any case with different views, the ROIs reconstructed by the PWLS-TV-IFR method are always the best one compared with the other methods. Although PWLS-TV does not work as well as PWLS-TV-IFR, it outperforms the ART-TV.  Figures 12-14 show the images and zoomed interior ROIs reconstructed by FBP-EP, PWLS-TV and PWLS-TV-IFR using the truncated normal-mAs projection data with different views (225 mAs, 450, 225 and 150 views), respectively. It is observed that the images reconstructed by FBP-EP contain severe noise which leads to a low visibility of details. The results reconstructed by the other two algorithms show a lower level of noise. For the FBP-EP method, the downsampling from 450 to 225 and 150 views results in some streak artifacts. The PWLS-TV results contain less artifacts but with degraded resolution and missing details. The reconstructions using the presented PWLS-TV-IFR method exhibit noticeable gains compared with those from FBP-EP and PWLS-TV in suppressing noise and preserving the smallscale and low-contrast structures, as indicated by the arrows in the zoomed ROIs.    respectively. Comparing with the normal-dose reconstructions, the images reconstructed by FBP-EP display a higher noise level, and contain some streak artifacts, especially in the 180-views case. On the other hand, the reconstructions using the PWLS-TV method exhibit less artifacts and noise but with worse resolution and detail loss. The less views that are used, the more the details become degraded. As shown in the zoomed ROIs, the reconstructions using the presented PWLS-TV-IFR method result in better performance than the PWLS-TV method in achieving higher small-scale and low-contrast detestability.

Conclusion and discussion
In this paper, a statistical interior tomography approach which includes feature refinement has been proposed. Two key ingredients were included in this approach. One is the framework of TV based interior tomography under the criteria of penalized weighed least-square, which takes into account the statistical nature of interior projections. The other is the feature refinement step which is inserted after each iteration of PWLS-TV to restore the missing fine features. The reconstruction results from both digital phantom and in vivo Micro-CT dataset show that the present PWLS-TV-IFR approach outperforms the existing methods, which demonstrates its useful potential for interior tomography under the low dose conditions.
Under the ideal noise-free conditions of the truncated projection, ART-TV can achieve accurate image reconstruction (Han et al 2009). However, in practice, several complex factors, such as data calibration and noise, will unavoidably violate the data piecewise constant assumption. This poor condition results in the failure of the ART-TV algorithm in noisy interior tomography. To improve interior tomography for further radiation dose reduction, it is natural to consider the statistical nature of the local projection data. Based on these observations, a maximization of a posteriori (MAP)  and a ML-EM based interior image reconstruction methods (Xu and Tsui 2011) have been proposed. Both of these two statistical related methods were developed by combining the two existing approaches, namely (1) when the full knowledge of a small region inside the interior ROI is known, and (2) when the complete interior ROI is piecewise constant. Contrast to above work, this paper is only based on the piecewise constant assumption, which has less limitations. A TV based interior tomography is performed under the criteria of penalized weighed least-square (PWLS), wherein the weighted least-square (WLS) term is built based on an accurate relationship between the variance and mean of projection data. The aim of the statistical image reconstruction process is to estimate the 'best' image, i.e. the image which has the highest probability to match the measured projection data and obtain exact interior reconstruction result from truncated noisy projection. It should be noted that the reconstruction results show that noise, truncated artifacts and streak artifacts can be well suppressed in the ultra-low-dose protocol combining the low-mAs and sparse-view simultaneously. A possible reason is that the constructed WLS fidelity term can suppress statistical noise and the TV-based penalty term can suppress both truncated and streak artifacts significantly. The hyper-parameter β is very important to achieve a better result, which can balance the fidelity term and penalty term. In this paper, the parameter β was selected by considering the noiseresolution tradeoff.
Another contribution worth mentioning is the feature refinement step. As pointed out by Osher et al ideally the denoising procedure would decompose the input image into a noisefree image and a residual image which only contains the additive noise. However, it is not fully attainable practically. In the noise and artifacts removing procedure, some features are often sacrificed (Tang et al 2009, Gao et al 2014. Local ROIs usually include rich details, which are important for radiologists in making diagnose. In this regard, recovering the losing features in interior tomography is highly desirable. To this end, describing the features is critical and can be accomplished by the objective and subjective schemes. For medical imagining applications in which images are ultimately to be viewed by human beings (radiologists), the only 'correct' method of quantifying diagnosis quality is through subjective evaluation. However, it is usually too inconvenient and time-consuming in practice. The same is true on feature description. Therefore, we chose an objective feature descriptor scheme that can automatically describe the features from the image. With the help of this descriptor, we propose to add a feature refinement step to recover the losing structure information and preserve details. It should be mentioned that the present feature refinement can be regarded as an enhancement filter, which is not only related to the TV-minimization but also related to the whole objective function.
Binary optimal reconstruction strategy was used for solving the objective function of our method. The initial guess for the statistical iteration used a projection data extrapolation method, which ensures continuity at the boundary between the measured and extrapolated regions. Minimizing the PWLS-TV the objective function was followed by the interior feature refinement step at each iteration. Just like many existing one-step-late (OSL) algorithms, the present algorithm also lacks strict global convergence proof (O'Sullivan and Benac 2007). However, the difference images between iterations as shown in figure 6 can provide an empirical interpretation of the convergence property of our proposed framework. Additionally, the RMSE and PSNR as a function of the number of iterations as shown in figure 7 demonstrate that our algorithm converges to a stable and better solution compared with other methods.
There are three limitations in our proposed framework. Firstly, the proposed PWLS-TV-IFR adds computational burden due to the added feature refinement step compared with the PWLS-TV approach. Secondly, in the protocol when the level of mAs and number of views are all very low, the image reconstructed by TV-minimization displays unphysical patchy artifacts and is different from the FBP reconstruction, especially in the in vivo datasets. It leads to degraded resolution of the reconstructed images. How to keep the resolution when reducing the noise is still an interesting topic for the ultra-low-dose image reconstruction (Wang et al 2008b, Gao et al 2014. Thirdly, it is difficult to combine the IFR filter into the cost-function with a unify formula.
Although the current work only focuses on 2D case, the proposed interior reconstruction strategy can be naturally extended to 3D case since the inherent methodology is the same. But it is worth mentioning that in practical 3D interior reconstruction, the associated projections would be severely transverse and/or axial truncated in circle or helical scans. Thus, an advanced 3D interior reconstruction method should be explored to yield a desired result. This is also an interesting topic for future research.