Iterative reconstruction for sparse-view X-ray CT using alpha-divergence constrained total generalized variation minimization

Abstract

BCKGROUND:

Accurate statistical model of the measured projection data is essential for computed tomography (CT) image reconstruction. The transmission data can be described by a compound Poisson distribution upon an electronic noise background. However, such a statistical distribution is numerically intractable for image reconstruction.

OBJECTIVE:

Although the sinogram data is easily manipulated, it lacks a statistical description for image reconstruction. To address this problem, we present an alpha-divergence constrained total generalized variation (AD-TGV) method for sparse-view x-ray CT image reconstruction.

METHODS:

The AD-TGV method is formulated as an optimization problem, which balances the alpha-divergence (AD) fidelity and total generalized variation (TGV) regularization in one framework. The alpha-divergence is used to measure the discrepancy between the measured and estimated projection data. The TGV regularization can effectively eliminate the staircase and patchy artifacts which is often observed in total variation (TV) regularization. A modified proximal forward-backward splitting algorithm was proposed to minimize the associated objective function.

RESULTS:

Qualitative and quantitative evaluations were carried out on both phantom and patient data. Compared with the original TV-based method, the evaluations clearly demonstrate that the AD-TGV method achieves higher accuracy and lower noise, while preserving structural details.

CONCLUSIONS:

The experimental results show that the presented AD-TGV method can achieve more gains over the AD-TV method in preserving structural details and suppressing image noise and undesired patchy artifacts. The authors can draw the conclusion that the presented AD-TGV method is potential for radiation dose reduction by lowering the milliampere-seconds (mAs) and/or reducing the number of projection views.

Keywords

Sparse-view X-ray CT iterative reconstruction alpha-divergence total generalized variation

1 Introduction

X-ray CT has been widely used in clinical diagnosis, image-guided interventions, and radiotherapy because of high-resolution anatomic imaging. However, the extra radiation dose in CT examination is harmful which may induce genetic, cancerous and other diseases [1, 2]. Therefore, minimizing radiation dose is one of the major topic in CT fields.

Over the years, two strategies have been widely used to reduce radiation dose: the first one is to lower the milliampere-seconds (mAs) in data acquisition protocols [3 –10] and the second one is to decrease the number of projections per rotation around the object [11 –18]. The first strategy will result in noisy sinogram data and the reconstructed image quality will be degraded [5 , 10]. The latter strategy will cause serious streaking artifacts in the reconstructed images because of insufficient angular sampling [14 –16].

To achieve high-quality images, various advanced CT image reconstruction methods have been reported [8 , 10–18]. Statistical iterative reconstruction (SIR) methods have shown great potential to reduce noise and artifacts as compared with filtered back-projection (FBP) algorithm [5 , 10]. SIR method is composed of data-fidelity and regularization terms. The fidelity term models the statistics of measured projection data and the regularization term reflects the prior information of the reconstructed image. The projection data can be modelled accurately using the compound poisson noise of polyenergetic x-ray spectrum plus Gaussian electronic noise [10]. However, such a complicated statistical distribution is numerically intractable for image reconstruction [19]. The projection data after linearity calibration on the transmission data can be easily manipulated, but lacks a statistical description for image reconstruction [20]. To address this problem, we present an alpha-divergence constrained total generalized variation (AD-TGV) method for CT image reconstruction. The novelty of the AD-TGV method is twofold. First, alpha-divergence fidelity and TGV regularization are incorporated into one framework for CT image reconstruction. Alpha-divergence is adapted to measure the discrepancy between projection data and corresponding estimator; the higher order derivative of the images is incorporated into TGV regularization could yield visually pleasant images with more continuous boundaries and less patchy artifacts. Second, a modified alternating optimization algorithm is proposed to minimize the objective function of the AD-TGV method with robust convergence.

The rest of this paper is organized as follows. The AD-TGV model for CT image reconstruction is presented in section 2. The experimental setup and evaluation metrics are also presented in this section. The results are reported in section 3. Finally, the discussion and conclusion are presented in section 4 and 5, respectively.

2 Materials and methods

2.1 Alpha-divergence

The alpha-divergence is a generalization of KL divergence and χ²-divergence which is used to measure the difference between probability distributions [7 , 22]. For any tow distributions r (x) and t (x), the alpha-divergence is defined as follows: $D_{α} (r, t) = \frac{1}{α (1 - α)} \int [α r + (1 - α) t - r^{α} t^{1 - α}] dx$ (1) where α ∈ (- ∞ , + ∞). We denote $D_{0} (r, t) = lim_{α \to 0} D_{α} (r, t)$ and $D_{1} (r, t) = lim_{α \to 1} D_{α} (r, t)$ to ensure the continuity of alpha-divergence on α.

2.2 TGV regularization

TGV was first proposed in [23] for image denoising, which is capable to measure image characteristics up to a certain order of differentiation. The second-order TGV of a given image f is defined as follows: ${TGV}_{α}^{2} (f) = min_{w \in BD (Ω)} α_{1} \int_{Ω} | \nabla f - w | dx + α_{0} \int_{Ω} | ɛ (w) | dx$ (2) where BD (Ω) denotes the space of vector fields of bounded deformation, ∇ represents the gradient operator, and $ɛ (w) = \frac{1}{2} (\nabla w + \nabla w^{T})$ denotes the symmetrized derivative. TGV provides a way to balance the first and second derivatives through the ratio of positive weights α₀ and α₁. In contrast, TV only takes the first derivative into account. For compatibility with TV, we set α₁ = 1 throughout the study. Moreover, based on our previous work [15], the value α₀ = 3 is suitable for CT image reconstruction.

2.3 Reconstruction model

The CT measurement can be formulated as a discrete linear system under a monochromatic beam: $y = Hf$ (3) where y = (y₁, y₂, ⋯ , y_M) ^T represents the obtained projection data, f = (f₁, f₂, ⋯ , f_N) ^T is the vector of attenuation coefficients to be reconstructed, T denotes the matrix transpose, and H represents the system or projection matrix with the size of M × N. The element h_ij of H is the length of the intersection of projection ray I with pixel j. CT image reconstruction problem can be formulated as to retrieve the unknown image f from the projection y according to Equation (3).

To find an approximation of Equation (3), alpha-divergence was used to measure discrepancy between y and Hf: $D_{α} (y, Hf) = \frac{1}{α (1 - α)} \sum_{i = 1}^{M} [α y_{i} + (1 - α) (Hf)_{i} - y_{i}^{α} (Hf)_{i}^{1 - α}]$ (4)

The objective image is the minimum of D_α (y, Hf) under the constraints f ≥ 0 (f_j ≥ 0 for all j). As opposed to minimizing the Equation (4) directly, we propose the following AD-TGV minimization for CT image reconstruction $min_{f \geq 0} D_{α} (y, Hf) + β TGV (f)$ (5) where β > 0 is a hyper-parameter to balance the alpha-divergence fidelity and TGV regularization.

2.4 Numerical algorithm

Since the objective function in the reconstruction problem (5) is convex, it can be solved by iteratively performing the following two steps: $\begin{matrix} (P 1) : u^{k + 1} = f^{k} - \frac{1}{μ} \nabla D_{α} (y, {Hf}^{k}) \\ (P 2) : f^{k + 1} = arg min_{f} \frac{μ}{2 β} {∥ f - u^{k + 1} ∥}_{2}^{2} + {TGV}_{α}^{2} (f) \end{matrix}$ (6) where k is the index for iteration steps, μ > 0 is a constant and u is an auxiliary vector. This two-step iteration is actually the so-called proximal forward-backward splitting algorithm, whose convergence result has been given in Combettes and Wajs’s work [24].

Problem (P1) is a gradient descent method with a step-size 1/μ of for the minimization problem $min_{f} D_{α} (y, Hf)$ . The introduction of μ is to ensure the numerical stability of the gradient descent method. The efficiency of the algorithm can be improved by substituting this one step gradient descent with a surrogate function (SF) algorithm in [21, 22], the overall convergence can be enhanced. The detailed implementation of this SF algorithm is performed as follows: $u_{j}^{k + 1} = f_{j}^{k} {\sum_{i = 1}^{M} {h_{ij} (\frac{y_{i}}{{({Hf}^{k})}_{i}})}^{α}}^{1 / α}, j = 1, 2, \dots, N$ (7) where the superscript k represents the iteration steps. The expression of the update for α = 0 can be performed as follows: $u_{j}^{k + 1} = f_{j}^{k} \prod_{i = 1}^{M} {(\frac{y_{i}}{{({Hf}^{k})}_{i}})}^{h_{ij}}, j = 1, 2, \dots, N .$ (8)

Prime-dual algorithm described in [25] was employed for solving the problem (P2), which can be rewritten as the following convex-concave saddle-point problem: $min_{f \in F, w \in W} max_{p \in P, q \in Q} \frac{1}{2 λ} {∥ f - u^{k + 1} ∥}_{2}^{2} + p^{T} (\nabla f - w) + q^{T} ɛ (w)$ (9) where λ = β/μ, $F = ℝ^{NN}$ , $W = ℝ^{2 NN}$ , _{P = {p ∈R^2NN|_∥p ∥∞ ≤α₁ }}, $Q = {q \in ℝ^{3 NN} | {∥ q ∥}_{\infty} \leq α_{0}}$ . We assume that proj_P (p) and proj_Q (q) denote the Euclidean projectors onto the convex sets P and Q, respectively. By the pointwise operation, the associated projections can be expressed as: ${proj}_{P} (p) = \frac{p}{max (1, | p | / α_{1})}, {proj}_{Q} (q) = \frac{q}{max (1, | q | / α_{0})} .$ (10)

The proximal map prox₁ (f) is defined as ${prox}_{1}^{τ} (f) = \underset{\tilde{f} \in F}{arg min} \frac{{∥ \tilde{f} - τ u ∥}_{2}^{2}}{2 λ} + \frac{{∥ \tilde{f} - f ∥}_{2}^{2}}{2 τ} = \frac{λ f + τ u}{λ + τ} .$ (11)

In summary, the implementation of the AD-TGV method can be described as follows: $\begin{matrix} 1 : initial f^{0}, w^{0}, {\bar{w}}^{0}, p^{0} and q^{0}; \\ 2 : initial λ, ρ, τ, α_{0}, α_{1} and k : = 0; \\ 3 : while stop criterion is not met; \\ 4 : for j = 1, 2, \dots, N; \\ 5 : if α = 0; \\ 6 : u_{j}^{k + 1} = f_{j}^{k} \prod_{i = 1}^{M} {(\frac{y_{i}}{{({Hf}^{k})}_{i}})}^{h_{ij}}, \\ 7 : else \\ 8 : u_{j}^{k + 1} = f_{j}^{k} {\sum_{i = 1}^{M} {h_{ij} (\frac{y_{i}}{{({Hf}^{k})}_{i}})}^{α}}^{1 / α}, \\ 9 : end if; \\ 10 : end for; \\ 11 : z^{0} = u^{k + 1}, {\bar{z}}^{0} = u^{k + 1} \\ 12 : for l = 0, \dots, L - 1; \\ 13 : p^{l + 1} = {proj}_{P} (p^{l} + ρ (\nabla {\bar{z}}^{l} - {\bar{w}}^{l})), \\ 14 : q^{l + 1} {= proj}_{Q} (q^{l} + ρ ɛ ({\bar{w}}^{l})), \\ 15 : z^{l + 1} {= prox}_{1}^{τ} (z^{l} + τ {div}_{1} p^{l + 1}), \\ 16 : {\bar{z}}^{k + 1} = 2 z^{l + 1} - z^{l}, \\ 17 : w^{l + 1} = w^{l} + τ (p^{l + 1} {+ div}_{2} q^{l + 1}), \\ 18 : {\bar{w}}^{l + 1} = 2 w^{l + 1} - w^{l}, \\ 19 : end for; \\ 20 : w^{0} = w^{L}, {\bar{w}}^{0} = {\bar{w}}^{L}, p^{0} = p^{L}, {\bar{q}}^{0} = {\bar{q}}^{L}; \\ 21 : f^{k + 1} = \max {z^{L}, 0}; \\ 22 : end if stop criterion satisfies . \end{matrix}$

In this study, an initial estimate of f⁰ was set to be the preliminary image reconstructed by the FBP method with ramp filter. The initial estimation of $w^{0}, {\bar{w}}^{0}, p^{0} and q^{0}$ were set to be zeros. In line 2, ρ and τ are step-size, which can be determined according to the work in [25]. The stop criterion is $\frac{| | f^{k + 1} - f^{k} | |_{2}}{| | f^{k} | |_{2}} \leq 10^{- 6}$ .

2.5 Parameter selections

2.5.1 Selection of the value of α

Determining an optimum α to describe the statistical distribution of projection data is very important for AD-TGV image reconstruction. It is known that if α = 2, then the α-divergence measure has the following form: $D_{2} (r, t) = \frac{1}{2} \sum_{i = 1}^{N} \frac{{(r_{i} - t_{i})}^{2}}{t_{i}}$ (12) where r = (r₁, r₂, ⋯ , r_N) ^T, t = (t₁, t₂, ⋯ , t_N) ^T. It can be observed that the case of is the weighted least-squares (WLS) model in [26, 27]. Meanwhile, based on Minka’s fixed-point scheme [28], KL-distance is equivalent to alpha-divergence at α = 1. With above discussion, in this study, we focus on the value of α between 1.0 and 2.0.

2.5.2 Selection of the value of λ

The penalty parameter λ balances the fidelity and regularization terms in Equation (9). The penalty parameter determines the smoothness of the reconstructed image. Practically, the penalty parameter is empirically determined through a visual inspection for eye-appealing result. In this work, the penalty parameter λ was determined case by case on the basis of the noise level of reconstructed image and projection.

2.6 Experimental data

2.6.1 XCAT phantom

A slice of digital XCAT phantom with a tumor lesion (Fig. 1(a)) was used in this study [29]. The imaging geometry is as follows: (1) the distance from the detector arrays to the source is 1040 mm; (2) the distance from the rotation center to the source is 570 mm; (3) the number of an arc detector bin is 672; and (4) the space of each bin is 1.407 mm. The size of the reconstructed images is512×512.

Fig.1

Illustration of a slice of digital XCAT phantom (a), an anthropomorphic torso phantom (b), and the gold standard (c) for the anthropomorphic torso phantom reconstructed by FBP algorithm with the averaged sinogram data of 150 times repeated scans at a protocol of 100 mAs and 120 kVp.

Similar to our previous studies [8, 15], the noisy transmission data I was simulated as follows: $I = Poisson (I_{0} exp (- \hat{y})) + Normal (0, σ_{e}^{2})$ (13) where I₀ is the incident photon intensity, $\hat{y}$ is the noise-free sinogram data, and $σ_{e}^{2}$ is the background electronic noise variance. In this study, I₀ and $σ_{e}^{2}$ were set to be 1.0×10⁶ and 10.0, respectively. Finally, the noisy sinogram data y was obtained by performing the logarithm transformation on the transmission data I.

2.6.2 Anthropomorphic torso phantom

An anthropomorphic torso phantom (Fig. 1(b)) was scanned by a Siemens SOMATOM Sensation 16 CT scanner at 40 mAs, 120 kVp. The sparse-view projection was evenly extracted from the acquired sinogram data to ensure the same noise level for each projection view.

2.6.3 Patient data

The patient data were obtained with a patient consent for a chest CT study. The experimental data were acquired by a Siemens SOMATOM Sensation 16 CT scanner. The scanning parameters were as follows: the time per gantry rotation is 0.5 s, the tube current is 200 mA, and the tube voltage is 120 kVp.

2.7 Performance evaluation

2.7.1 Image reconstruction accuracy

The relative root mean square error (rRMSE), which calculates the accuracy of the resulting image, was utilized to measure the quality of the reconstructed image. The rRMSE is defined as follows: $rRMSE = \sqrt{\frac{\sum_{m = 1}^{N} {(f (m) - f_{xtrue} (m))}^{2}}{\sum_{m = 1}^{N} {(f_{xtrue} (m))}^{2}}}$ (14) where f (m) and f_xtrue (m) represents the reconstructed and true attenuation coefficient at pixel m, respectively. N is the total number of pixels of the desired image.

2.7.2 Image similarity

To explore the performance of the proposed method, the universal quality index (UQI) [30] was used to conduct an ROI based analysis by evaluating the similarity between the reconstructed image and gold standard. UQI ranges from 0 to 1, and its value closer to 1 suggests a better similarity between the reconstructed and gold standard images.

2.7.3 Comparison method

To validate and evaluate the performance of AD-TGV method, the AD-TV method described in our previous work [7] was adopted for comparison. The cost function of AD-TV method is written as follows: $min_{f \geq 0} D_{α} (y, Hf) + β TV (f)$ (15) where β > 0 is a hyper-parameter to balance the fidelity term (i.e., the first term of Equation (15)) and TV regularization term.

3 Results

3.1 XCAT phantom study

3.1.1 Determination of the value of α

The influence of the parameter α on the reconstruction accuracy should be considered carefully because the performance of the present AD-TGV method heavily depends on the α setting. In this study, the AD-TGV method was validated quantitatively from different values of α with the range of α ∈ [1, 2]. The rRMSE versus the value of α for the AD-TGV method from different sparse-views (30-, 40-, and 60-view projections) data were shown in Fig. 2. We can observe that rRMSE increases with the value of alpha increases. Hence, we just focus on the case of α = 1, which exhibits the lowest rRMSE.

Fig.2

rRMSE versus the value of alpha for the AD-TGV method from different projections: (a) 30-view projection; (b) 40-view projection; (c) 60-view projection.

3.1.2 Convergence analysis

The rRMSE versus the iteration steps for the present AD-TGV method from different projection views (30-, 40- and 60-view) are displayed in Fig. 3. It can be observed that the proposed AD-TGV method can converge to a steady solution after enough iteration. In addition, the convergence speed is accelerated as the projection views increased with the same penalty parameter. The results show that the present AD-TGV method can yield a steadily convergence solution in terms of the rRMSE measurement.

Fig.3

rRMSE versus iteration steps for the AD-TGV method from different projections: (a) 30-view projection; (b) 40-view projection; and (c) 60-view projection.

3.1.3 Visual evaluation

The results reconstructed by different methods from the 30-, 40- and 60-view projections are shown in Fig. 4. Serious streak artifacts can be observed in the images reconstructed by the FBP method. The steak artifacts have been mostly suppressed in the images reconstructed by AD-TV method, but undesired artifacts can still be observed. The images reconstructed by AD-TGV method achieve more gains over the AD-TV method in terms of successful artifacts suppression and edges detail preservation.

Fig.4

The images reconstructed by the FBP with ramp filter, AD-TV and AD-TGV methods from 30-, 40- and 60-vews projections, respectively.

3.1.4 UQI study

To quantitatively evaluate the present AD-TGV method, the UQI were studied in this subsection. The ROI as indicated by a red rectangular window in Fig. 1(a) was selected to calculate the UQI values. The corresponding UQI values are shown in Fig. 5. We can observed that the UQI values of AD-TGV results are higher than the AD-TV results.

Fig.5

UQI values of images reconstructed by FBP, AD-TV, and AD-TGV methods from 30, 40, 60, 116, 232, 290, 387, 580, and 1160-projection views.

3.2 Anthropomorphic torso phantom study

3.2.1 Visual evaluation

The results reconstructed by different methods from the sinogram data acquired at 40 mAs with different projection views are shown in Fig. 6. We can observe that AD-TGV method achieves more gains than AD-TV method in terms of successful patchy artifacts suppression and structure detail preservation. Figure 7 depicts the associated vertical profiles of the reconstructed images across the 245th column, from the 233th row to the 283th row. The results suggest that the AD-TGV method achieved better matching results than the AD-TV method.

Fig.6

Images reconstructed by the FBP, AD-TV, and AD-TGV methods from 58-, 116-, 232-, 290- view projection at 40 mAs.

Fig.7

Profiles (245th column) of the images reconstructed by the AD-TV and AD-TGV methods from different projections at 40 mAs: (a) 58-view projection; (b) 116-view projection; (c) 232-view projection; (d) 290-view projection.

To further show the gains of the present AD-TGV method, the zoomed ROIs (as indicated by the squares in Fig. 1(c) are shown in Fig. 8. Serious artifacts can be found in the FBP results for all the cases, and AD-TGV method can achieve noticeable gains compared with AD-TV method in terms of staircase and patchy artifacts suppression.

Fig.8

The zoom-in views of images in Fig. 7.

3.2.2 UQI study

The gold standard image (Fig. 1(c)) was reconstructed by the FBP method with ramp filter from the averaged full view sinogram data of 150 times repeated scans at the protocol of 100 mAs and 120 kVp. The curve of UQI values versus the number of projection views is shown in Fig. 9. The UQI values of AD-TGV results are higher than that of the AD-TV results. Thus, the AD-TGV method achieves better matching results than the AD-TV method.

Fig.9

UQI value versus the numbers of projection views.

3.3 Patient study

The chest CT images reconstructed by different methods are displayed in Fig. 10. The image reconstructed by the FBP method with ramp filter from full view projection is shown in Fig. 10(a). The image reconstructed by the FBP method with ramp filter from 387 projection views is shown in Fig. 10 (b). The image reconstructed by the AD-TV method from 387 projection views is shown in Fig. 10 (c). The image reconstructed by the AD-TGV method from 387 projection views is shown in Fig. 10 (d). To visualize the difference between the AD-TV and AD-TGV methods, the zoomed ROIs of the images (Fig. 10) are shown in Fig. 11. The results show that the AD-TGV method achieves better noise reduction, artifacts suppression, and structure preservation than the AD-TV method.

Fig.10

Clinical chest CT image reconstruction results: (a) image reconstructed by the FBP method with ramp filter from full view projection at 100 mAs, 120 kVp; (b) image reconstructed by the FBP method with ramp filter from 387-view projection at 100 mAs, 120 kVp; (c) image reconstructed by the AD-TV method from 387-view projection at 100 mAs, 120 kVp; and (d) image reconstructed by AD-TGV method from 387-view projection at 100 mAs, 120 kVp.

Fig.11

The zoomed ROI in Fig. 10. (a) image reconstructed by the FBP method with ramp filter from the full view projection at 100 mAs, 120 kVp; (b) image reconstructed by the FBP method with ramp filter from 387-view projection at 100 mAs, 120 kVp; (c) image reconstructed by the AD-TV method from 387-view projection at 100 mAs, 120 kVp; and (d) image reconstructed by the AD-TGV method from 387-view projection at 100 mAs, 120 kVp.

4 Discussion

In this study, we proposed an alpha-divergence constrained total generalized variation minimization for sparse-view CT image reconstruction. The alpha-divergence is a generalization of KL divergence and χ²- divergence [30] which is usually used to measure the difference between two probability distributions. Some special cases in D_α (y, Hf) listed below have been frequently applied to iteration reconstruction: $D_{α = 1} (y, Hf) = \sum_{i = 1}^{M} [y_{i} ln \frac{y_{i}}{(Hf)_{i}} - y_{i} + (Hf)_{i}]$ (16) $D_{α = 2} (y, Hf) = \frac{1}{2} \sum_{i = 1}^{M} \frac{[y_{i} - (Hf)_{i}]^{2}}{(Hf)_{i}}$ (17) $D_{α = - 1} (y, Hf) = \frac{1}{2} \sum_{i = 1}^{M} \frac{[y_{i} - (Hf)_{i}]^{2}}{y_{i}}$ (18)

Ignoring the constants in Equation (16) it coincides with opposite log-likelihood function. Equation (17) and Equation (18) are the corresponding WLS models in [26, 27]. An accurate statistical distribution of the sinogram data is very important for high-quality CT image reconstruction. By the determined adequate α, the corresponding alpha-divergence can accurately describe the statistical distribution ofsinogram data.

TV regularization has been widely used in CT image reconstruction with noticeable gains [11, 17]. However, TV regularization often leads to noticeable staircasing effect and patchy artifacts because of the piecewise constant assumption [15, 32]. In this work, the TGV regularization was used to eliminate the staircasing effect and patchy artifacts. The performance of the presented AD-TGV method has been extensively validated and evaluated in terms of the noise reduction and structure detail preservation. As shown in the study of digital XCAT phantom, the performance of the AD-TGV method is impressive by quantitative measure and visual inspection. Furthermore, qualitative and quantitative studies were conducted using using physical phantom and patient data in realistic CT imaging. The results show that the AD-TGV method outperforms the AD-TV method in terms of visual inspection and UQI evaluation.

For the implementation of the AD-TGV method, several parameters should be optimally selected. This is an open problem for the iterative reconstruction methods. In this study, we just focus on the value of α = 1. Furthermore, the rRMSE versus the vale of α in the range of [10^-6, 2] for AD-TGV method with 60-view projection (XCAT phantom) was shown in Fig. 12. We note that the optimal value of α is around 0.15, and rRMSE can reach its minimum at α = 1 in the range of [1, 2]. Determining the optimal value of α in a greater interval is an interesting topic for future research. The parameters ρ and τ can be selected according to the scheme introduced in [25]. Similar to most SIR methods, the computational time of the proposed AD-TGV could be a challenge for its practical use because of multiple re-projection and back-projection operations. Nevertheless, with the development of fast computers and dedicated hardware [33], iterative reconstruction algorithms including the present AD-TGV method may be used for clinical CT reconstruction in the near future.

Fig.12

rRMSE versus the value of alpha for the AD-TGV method with 60-view projection (XCAT phantom).

5 Conclusion

In this work, we presented an iterative reconstruction method for sparse-view X-ray CT using alpha-divergence constrained TGV. The alpha-divergence is used to measure the discrepancy between the measured projection and estimation. The TGV regularization is proposed to eliminate the staircase and patchy artifacts which is often observed in TV regularization. Experimental results show thatAD-TGV method is more beneficial than AD-TV method in terms of noise reduction, artifact suppression and structure preservation. In other words, the presented AD-TGV method can potentially reduce the radiation dose by lowering the mAs and reducing the number of projection views.

Footnotes

Acknowledgments

This work was partially supported by the National Natural Science Foundation of China (No.81371544, No.61571214, No.11661007 and No.11661008), the National Science and Technology Major Project of the Ministry of Science and Technology of China (No. 2014BAI17B02), Science and Technology Program of Guangzhou (No.201510010039), Natural Science Foundation of Jiangxi Province(No.20161BAB212055), Science and Technology Program of Jiangxi Education Committee(GJJ150994), and NIH/NCI(# CA143111 and # 082402).

References

Brenner

D.J.

and Hall

E.J.

, Computed tomography–an increasing source of radiation exposure, N Engl J Med 357 (2007), 2277–2284.

Hall

E.J.

and Brenner

D.J.

, Cancer risks from diagnostic radiology, Br J Radiol 81 (2008), 362–378.

Bian

, et al., SR-NLM: A sinogram restoration induced non-local means image filtering for low-dose computed tomography, Comput Med Imaging Graph 37 (2013), 293–303.

, et al., An efficient augmented lagrangian method for statistical X-ray CT image reconstruction, PloS One 10 (2015), e0140579.

La Rivière

P.J.

and Billmire

D.M.

, Reduction of noise-induced streak artifacts in X-ray computed tomography through spline-based penalized-likelihood sinogram smoothing. IEEE Transactions on Medical Imaging 24 (2005), 105–111.

, et al., Low-dose computed tomography image restoration using previous normal-dose scan, Med Phys 38 (2011), 5713–5731.

, et al., Low-dose computed tomography image reconstruction by α-divergence constrained total variation minimization, Fully three-dimensional image reconstruction in radiology and nuclear medicine Proceedings Germany, 2011.

, et al., Iterative image reconstruction for cerebral perfusion CT using a pre-contrast scan induced edge-preserving prior, Phys Med Biol 57 (2012), 7519–7542.

Niu

, et al., Low-dose cerebral perfusion computed tomography image restoration via low-rank and total variation regularizations, Neurocomputing 197 (2016), 143–160.

10.

Wang

, et al., Penalized weighted least-squares approach to sinogram noise reduction and image reconstruction for low-dose X-ray computed tomography, IEEE Trans Med Imaging 25 (2006), 1272–1283.

11.

Chen

G.H.

, Tang

and Leng

, Prior image constrained compressed sensing (PICCS): A method to accurately reconstruct dynamic CT images from highly undersampled projection data sets, Med Phys 35 (2008), 660–663.

12.

Huang

, et al., Sparse angular CT reconstruction using non-local means based iterative-correction POCS. Comput Biol Med 41 (2011), 195–205.

13.

Liu

, et al., Adaptive-weighted total variation minimization for sparse data toward low-dose x-ray computed tomography image reconstruction, Phys Med Biol 57 (2012), 7923–7956.

14.

Bendory

and Feuer

, Sparse sampling in helical cone-beam CT perfect reconstruction algorithms, Journal of X-ray Science and Technology 24 (2016), 389–405.

15.

Niu

, et al., Sparse-view x-ray CT reconstruction via total generalized variation regularization, Phys Med Biol 59 (2014), 2997–3017.

16.

Sidky

E.Y.

, Kao

C.M.

and Pan

X.H.

, Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT, Journal of X-Ray Science and Technology 14 (2006), 119–139.

17.

Sidky

E.Y.

and Pan

, Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization, Phys Med Biol 53 (2008), 4777–4807.

18.

Tian

, et al., Low-dose CT reconstruction via edge-preserving total variation regularization, Phys Med Biol 56 (2011), 5949–5967.

19.

Lasio

G.M.

, Whiting

B.R.

and Williamson

J.F.

, Statistical reconstruction for x-ray computed tomography using energy-integrating detectors, Physics in Medicine and Biology 52 (2007), 2247–2266.

20.

Macovski

, Medical imaging systems, Prentice Hall 1983.

21.

Teng

and Zhang

, Iterative reconstruction algorithms with alpha-divergence for PET imaging, Comput Med Imaging Graph 35 (2011), 294–301.

22.

Teng

and Zhang

, Generalized EM-type reconstruction algorithms for emission tomography, IEEE Trans Med Imaging 31 (2012), 1724–1733.

23.

Bredies

, Kunisch

and Pock

, Total Generalized Variation, Siam Journal on Imaging Sciences 3 (2010), 492–526.

24.

Combettes

P.L.

and Pesquet

J.C.

, A proximal decomposition method for solving convex variational inverse problems, Inverse Problems 24 (2008), 065014.

25.

Knoll

, et al., Second order total generalized variation (TGV) for MRI, Magn Reson Med 65 (2011), 480–491.

26.

Anderson

J.M.

, et al., Weighted least-squares reconstruction methods for positron emission tomography, IEEE Trans Med Imaging 16 (1997), 159–165.

27.

Fessler

J.A.

, Penalized weighted least-squares image reconstruction for positron emission tomography, IEEE Trans Med Imaging 13 (1994), 290–300.

28.

Minka

, Divergence measures and message passin, Technical report, Microsoft Research, 2005.

29.

Segars

W.P.

, et al., 4D XCAT phantom for multimodality imaging research, Med Phys 37 (2010), 4902–4915.

30.

Wang

and Bovik

A.C.

, A universal image quality index, IEEE Signal Processing Letters 9 (2002), 81–84.

31.

Shun-ichi

, Differential-geometrical methods in statistics. Springer Science & Business Media, 2012.

32.

Tang

, Nett

B.E.

and Chen

G.H.

, Performance comparison between total variation (TV)-based compressed sensing and statistical iterative reconstruction algorithms, Phys Med Biol 54 (2009), 5781–5804.

33.

and Mueller

, Accelerating popular tomographic reconstruction algorithms on commodity PC graphics hardware, IEEE Transactions on Nuclear Science 52 (2005), 654–663.