Few-view CT image reconstruction using improved total variation regularization

Abstract

X-ray radiation is harmful to human health. Thus, obtaining a better reconstructed image with few projection view constraints is a major challenge in the computed tomography (CT) field to reduce radiation dose. In this study, we proposed and tested a new algorithm that combines penalized weighted least-squares using total generalized variation (PWLS-TGV) and dictionary learning (DL), named PWLS-TGV-DL to address this challenge. We first presented and tested this new algorithm and evaluated it through both data simulation and physical experiments. We then analyzed experimental data in terms of image qualitative and quantitative measures, such as the structural similarity index (SSIM) and the root mean square error (RMSE). The experiments and data analysis indicated that applying the new algorithm to CT data recovered images more efficiently and yielded better results than the traditional CT image reconstruction approaches.

Keywords

Dictionary learning total generalized variation regularization CT image reconstruction few-view

1 Introduction

CT is widely used in hospitals for diagnosis and intervention. However, X-ray radiation is harmful to human health, and many clinical studies have indicated that CT radiation doses beyond the normal range can cause diseases such as metabolic abnormalities and cancer. Because CT radiation doses accumulate over time, repeated CT scans increase the probability of cancer [1], highlighting the importance of low-dose CT. However, reducing the X-ray dose degrades the reconstructed image quality. Therefore, obtaining high-quality diagnostic CT images from either sparse view acquisition or low-dose radiation is an important research topic in the field of CT.

Two strategies exist to reduce the radiation dose: minimizing the milliampere-seconds (mAs) or peak kilovoltage (kVp) [2, 3] and reducing the number of projection views per rotation around the body [4, 5]. Various techniques, including better image reconstruction methods [5 –7] and scanning schemes [8 –10], have been extensively studied to reduce radiation doses in CT examinations. Among these, a method referred to as statistical iterative reconstruction (SIR) has been proposed. SIR optimizes the maximum-likelihood or penalized-likelihood function formulated according to the statistical characteristics of the projection data, which reduces the radiation dose compared to the filtered back-projection (FBP) reconstruction algorithm [11, 12]. In general, however, the conventional regularization term tends to produce unfavorable oversmoothing effects at the edge regions by penalizing the differences between local neighboring pixels. To solve this shortcoming, several edge-preserving regularization terms have been proposed [13, 14]. There are two typical examples: total variation (TV) regularization with the piecewise constant assumption (PCA) [15 –17] and total generalized variation (TGV) regularization [18, 19] adapted under the penalized weighted least-squares (PWLS) criterion. Unfortunately, neither images reconstructed by TV regularization nor by TGV regularization achieve satisfactory results from few-view projections.

Dictionary learning (DL) has received increasing attention in recent years and has been incorporated into iterative reconstruction frameworks by various algorithms [20 –23], such as adaptive dictionary-based statistical iterative reconstruction (ADSIR) and 3D dictionary learning (3D-DL) [24]. The combination of dictionary learning and iterative reconstruction in these methods can effectively improve the imaging performance. Other methods treat dictionary learning as postprocessing to reduce artifacts in images reconstructed with FBP, such as artifact suppression dictionary learning (ASDL) algorithms [25 –31].

To reconstruct a better image from a few-view projection (116-view projections had been used), we propose an improved CT statistical iterative reconstruction method in this paper. Our approach is to reconstruct intermediate images using TGV under the PWLS standard [18] and then postprocess the results using sparse coding dictionary learning to eliminate residual noise and produce clinically acceptable CT images. The proposed method is abbreviated as PWLS-TGV-DL for simplicity. The novelty of our approach is that PWLS-TGV-DL can produce better images with less noise and fewer patchy artifacts. Both qualitative and quantitative assessments of digital and physical models are presented and evaluated in terms of accuracy and resolution.

The remainder of this paper is organized as follows. Section 2 reviews the PWLS image reconstruction model, TGV model, dictionary learning, the PWLS-TGV-DL algorithm, its optimization algorithm, and its workflow. In Section 3, we present the experimental setup, the evaluation metrics, and the results of the PWLS-TGV-DL algorithm along with those of comparison methods using both simulation and physical experiments. Finally, Section 4 provides a discussion and a conclusion.

2 Methods

2.1 PWLS image reconstruction

The penalized weighted least-squares (PWLS) approach for iterative reconstruction of X-ray CT images was studied previously [32]. Based on the noise properties of CT projection data, the PWLS criterion for CT image reconstruction can be written as follows [33]: $m i n_{f \geq 0} {(y - H f)}^{T} \sum^{- 1} (y - H f) + R (f),$ (1) where f is a vector of the attenuation coefficients to be reconstructed (i.e., f = (f₁, f₂, ... , f_N) ^T ), and y represents the obtained sinogram data (projections after system calibration and logarithm transformation) (i.e., y = (y₁, y₂, ... , y_M) ^T ). Here, T denotes the matrix transpose, and the operator H represents the system or projection matrix with a size of M×N. The element h_ij is the length of the intersection of projection ray i with pixel j. ∑is a diagonal matrix with the ith element of $σ_{i}^{2}$ , which is the variance of the sinogram data y_i. R(f) represents the regularization term, and β is the parameter used to control the weight.

Based on previous works already available [34], in this study, the variance of $σ_{i}^{2}$ is defined by the following mean-variance relationship: $σ_{i}^{2} = \frac{1}{I_{0}} exp ({\bar{y}}_{ℓ}) (1 + \frac{1}{I_{0}} exp ({\bar{y}}_{ℓ}) (σ_{e}^{2} - 1.25)),$ (2) where I₀ is the incident X-ray intensity, ${\bar{y}}_{ℓ}$ is the mean of the sinogram data at bin i, and $σ_{e}^{2}$ is the background electronic noise variance.

2.2 Total generalized variation

The TGV algorithm was first proposed by Bredies et al. in the image denoising model [18], which is used to measure image characteristics up to a certain order of differentiation. Mathematically, TGV can be defined as follows: ${TGV}_{α}^{k} (u) = sup {\int_{Ω} u di v^{k} ν dx | ν \in C_{c}^{k} | (Ω, Sy m^{k} ({Ropf;}^{d})), {∥ di v^{l} ν ∥}_{\infty} ⩽ α_{l}},$ (3) where Ω is a bounded domain, div represents the divergence operator, v denotes the dual variable of the exact TGV definition, ℝ ;^d denotes the d-dimensional real space, l = 0, ... , k– 1, k∈ ℕ ; is an order of TGV, and α= (α₀, α₁, ... , α_k–1) represents the positive weights for TGV. Sym^k( ℝ ;^d) denotes the space of the symmetric k-tensors. The l-divergence of a symmetric k-tensor field is given by ${(di v^{l} ν)}_{η} = \sum_{γ \in M_{l}} \frac{l!}{γ!} \cdot \frac{\partial^{l} υ_{η + γ}}{\partial x^{γ}}$ for each component η∈M_k_–1. Here, M_k is the multi-index of order k (i.e., $M_{k} = {η \in {ℕ;}^{d} | \sum_{i = 1}^{d} η_{i} = k}$ ). The ∞-norm for symmetric k-vector fields is defined as follows: ${∥ ν ∥}_{\infty} = sup_{x \in Ω} {{(\sum_{η \in M_{k}} \frac{k!}{η!} ν_{η} {(x)}^{2})}^{\frac{1}{2}}}$ (4)

In this work, we focus only on the second-order TGV, which can be mathematically defined as ${TGV}_{α}^{2} (u) = sup {\int_{Ω} u di v^{2} ν dx | ν \in C_{c}^{2} (Ω, S^{d \times d}), {∥ ν ∥}_{\infty} ⩽ α_{0}, {∥ div ν ∥}_{\infty} ⩽ α_{1}},$ (5) where S^d^×d denotes the space of the symmetric d×d matrices, and α= (α₀, α₁) is assumed to be positive. Furthermore, the first and second divergences of the symmetric matrix can be calculated as follows: ${(div ν)}_{i} = \sum_{j = 1}^{d} \frac{\partial ν_{ij}}{\partial x_{j}}, di v^{2} ν = \sum_{i = 1}^{d} \frac{\partial^{2} ν_{ii}}{\partial x_{i}^{2}} + \sum_{i < j} \frac{\partial ν_{ij}}{\partial x_{i} \partial x_{j}} .$ (6)

The ∞-norm in Eq. 5 can be calculated by ${∥ ν ∥}_{\infty} = sup_{x \in Ω} {(\sum_{i = 1}^{d} {| ν_{ii} (x) |}^{2} + 2 \sum_{i < j} {| ν_{ij} (x) |}^{2})}^{\frac{1}{2}}$ (7) and the second-order TGV can be calculated as follows: ${TGV}_{α}^{2} (u) = min_{w} α_{1} \int_{Ω} | \nabla u - w | dx + α_{0} \int_{Ω} | ɛ (w) | dx .$ (8)

The minimum is taken over all vector fields in Ω, and the weak symmetrized derivative $ɛ (w) = \frac{1}{2} (\nabla w + \nabla w^{T})$ is a matrix-valued Radon measure. Eq. 8 balances the first and second derivatives by using the weights α₀ and α₁.

2.3 Dictionary learning

Let N and K be integers and ℝ ; be the real space. Given a 2D image f of size $\sqrt{N} \times \sqrt{N}$ can be expressed as a N dimensional vector f∈ ℝ ; ^N ^×1 and a vector f_j_,k representation of an image patch of size $\sqrt{n} \times \sqrt{n}$ , which can be expressed as a n dimensional vector f_j_,k∈ ℝ ; ⁿ ^×1, N > > n, each image patch f_j_,k is uniquely indexed by the location of its top-left corner pixel (j, k) in image f. A matrix D ∈ ℝ ; ⁿ ^×K exists whose column d_k =∈ ℝ ; ⁿ ^×1(k = 1, ... ,K), K > > n is an n-dimensional vector as a dictionary, which is also called an atom. The dictionary is usually redundant or overcomplete. Assume that each patch f_j_,k can be exactly or approximately represented as a sparse linear combination Dα _j _,k of dictionary atoms; that is, ${∥ R_{j, k} f - D α_{j, k} ∥}_{2}^{2} ⩽ ɛ,$ (9) where ɛ≥0 is a small error bound, and that the representation vector α _j _,k∈ ℝ ;^K ×1 has few nonzero entries, ∥α_j,k ∥ ₀ ⪡ n ⪡ K, where ∥ · ∥ ₀ is the l₀-norm. The matrix R_j_,k ∈ ℝ ; ⁿ ^×N denotes the operator that extracts a patch f_j_,k from f as f_j_,k = R_j_,kf. The dictionary redundancy implies that the number of atoms is larger than the length of an atom.

Finding a sparse representation α _j _,k∈ ℝ ; ^K ^×1 of an image patch f_j_,k∈ ℝ ; ^N ^×1 with respect to a given dictionary is equivalent to solving the following optimization problem: $min_{α_{j, k}} {∥ α_{j, k} ∥}_{0} s . t . {∥ R_{j, k} f - D α_{jk} ∥}_{2}^{2} ⩽ ɛ .$ (10)

Equation 10 can be rewritten in the following unconstrained form by the Lagrange method $min_{α_{j, k}} {∥ R_{j, k} f - D α_{j, k} ∥}_{2}^{2} + ν_{j, k} {∥ α_{j, k} ∥}_{0},$ (11) where v_j, k is the Lagrange multiplier. Equations 10 and 11 are equivalent when a suitable parameter has been chosen. However, because solving Equation 11 directly is NP-hard, an approximate alternative strategy is desirable. Several greedy algorithms are available to accomplish this, such as the matching pursuit (MP) [35] and orthogonal matching pursuit (OMP) algorithms [36, 37]. In order to facilitate calculation, the l₀-norm is usually represented by the l₁-norm to make the problem convex and manageable by a basis pursuit (BP) algorithm [38].

2.4 Proposed new model

We formulate the CT reconstruction problem using DL and discuss the strategy for dictionary construction. Inspired by the studies of TGV in image reconstruction [19], we provide the following cost function of our algorithm: $mi n_{f ⩾ 0} {(y - Hf)}^{T} G^{- 1} (y - Hf) + β_{2} {TGV}_{α}^{2} (f) + β_{3} (\sum {∥ R_{j, k} f - D α_{jk} ∥}_{2}^{2} + \sum ν_{j, k} {∥ α_{j, k} ∥}_{0}),$ (12) where $G = \frac{1}{β_{1}} H H^{T} + \sum$ , and β₁, β₂ and β₃ are hyperparameters to balance the fidelity and regularization terms, respectively. We use the modified weighted least-squares approach in Eq. 12 instead of the weighted least-squares technique in Equation 1. By introducing a vector μ, we have $min_{μ} {(y - H μ)}^{T} \sum^{- 1} (y - H μ) + β_{1} {∥ μ - f ∥}_{2}^{2} = {(y - Hf)}^{T} G^{- 1} (y - Hf) .$ (13)

Therefore, the minimization problem in Equation 12 can be rewritten as follows: $min_{μ, f} \begin{matrix} {(y - H μ)}^{T} \sum^{- 1} (y - H μ) + β_{1} {∥ μ - f ∥}_{2}^{2} + β_{2} {TGV}_{α}^{2} (f) \\ + β_{3} (\sum {∥ R_{j, k} f - D α_{jk} ∥}_{2}^{2} + \sum ν_{j, k} {∥ α_{j, k} ∥}_{0}) \end{matrix} .$ (14)

2.5 Algorithm to solve the proposed model

Because Equation 14 includes several variables, we propose an improved alternative optimization method that involves three subproblems, which are mathematically defined as follows:

Subproblem A: $μ = arg mi n_{μ} {(y - H μ)}^{T} \sum^{- 1} (y - H μ) + β_{1} {∥ μ - f ∥}_{2}^{2}$ (15)

This subproblem is a classic statistical image reconstruction problem. The solution to this problem in the implementation has been studied by Elbakri and Fessler et al. [11], and this solution is a separable paraboloidal surrogate algorithm written as follows: $μ_{j}^{k + 1} = μ_{j}^{k} - \frac{\sum_{i = 1}^{M} ((1 / - σ_{i}^{2}) h_{ij} ({[H μ^{k}]}_{i} - y_{i})) + β_{1} (μ_{j}^{k} - f_{j}^{k})}{\sum_{i = 1}^{M} ((1 / - σ_{i}^{2}) h_{ij} \sum_{t = 1}^{N} h_{it}) + β_{1}},$ (16) where the superscript k = 1, 2, ... , K represents the iteration index.

Subproblem B: $f = arg min_{f} β_{1} {∥ μ - f ∥}_{2}^{2} + β_{2} {TGV}_{α}^{2} (f)$ (17)

A corresponding discrete version to solve this problem was studied by Bredies et al. [39], which can be rewritten as follows: $min_{f \in F, w \in W} \frac{1}{2 λ} {∥ μ - f ∥}_{2}^{2} + α_{1} {∥ \nabla f - w ∥}_{1} + α_{0} {∥ ɛ (w) ∥}_{1},$ (18) where λ = β₂/ - 2β₁, F = ℝ ; ^NN , and W = ℝ ;^2NN, and the differential operators div, ɛ, and ∇ are approximated by using first-order finite differences. These operators are chosen adjoint to each other, i.e., ${(div)}_{2}^{*} = - ɛ$ , ${(div)}_{1}^{*} = - \nabla$ . A Chambolle– Pock first-order primal dual algorithm was provided by Chambolle and Pock et al. as a solution to Equation (18), which is written as follows: $min_{f \in F, w \in W} max_{p \in P, q \in Q} \frac{1}{2 λ} {∥ μ - f ∥}_{2}^{2} + p^{T} (\nabla f - w) + q^{T} ɛ (w),$ (19) where p and q are bivariate, and the set associated with these variables is given byP ={ p ∈ ℝ ; ^2NN | ∥ p ∥ _∞ ⩽ α₁ },Q ={ q ∈ ℝ ; ^3NN | ∥ q ∥ _∞ ⩽ α₀ }. Let ${proj}_{P} (\tilde{p})$ and ${proj}_{Q} (\tilde{q})$ denote the values representing the Euclidean projector on the convex sets P and Q, respectively. Then, the associated projection can be calculated as follows: $pro j_{P} (\tilde{p}) = \frac{\tilde{p}}{max (1, | \tilde{p} | / - α_{1})}, pro j_{Q} (\tilde{q}) = \frac{\tilde{q}}{max (1, | \tilde{q} | / - α_{0})} .$ (20)

Likewise, another ${prox}_{1} (\tilde{f})$ can be written as $pro x_{1} (\tilde{f}) = arg min_{f \in F} \frac{{∥ f - τ μ ∥}_{2}^{2}}{2 λ} + \frac{{∥ f - \tilde{f} ∥}_{2}^{2}}{2 λ} = \frac{λ f + τ μ}{λ + τ} .$ (21)

Subproblem C: $f = arg min_{α_{j, k}} {∥ R_{j, k} f - D α_{j, k} ∥}_{2}^{2} + ν_{j, k} {∥ α_{j, k} ∥}_{0}$ (22)

This subproblem aims at finding the sparsest representation α _j _,k with respect to the dictionary D and updating an image with respect to an adaptive dictionary. DL-based image reconstruction in the few-view case was studied by Xu et al. [40] and Liao and Sapiro [41]. These studies used the orthogonal matching pursuit (OMP) algorithm [42] and the K-singular value decomposition (K-SVD) algorithm [23] to solve this problem. Then, the DL step was incorporated in the iterative reconstruction procedure to update an intermediate image.

The entire PWLS-TGV-DL algorithm can be summarized as follows:

Initialization: $μ^{0}, f^{0}, {\bar{f}}^{0}, w^{0}, {\bar{w}}^{0}, p^{0} and q^{0}$ ;

Initialization: β₁, β₂, β₃, ρ, τ, α₀, α₁, and k = 0;

While the termination criterion is not met;

For j = 1, 2, ... , N;

$μ_{j}^{k + 1} = μ_{j}^{k} - \frac{\sum_{i = 1}^{M} ((1 / - σ_{i}^{2}) h_{ij} ({[H μ^{k}]}_{i} - y_{i})) + β_{1} (μ_{j}^{k} - f_{j}^{k})}{\sum_{i = 1}^{M} ((1 / - σ_{i}^{2}) h_{ij} \sum_{t = 1}^{N} h_{it}) + β_{1}},$

End For;

For n = 0, 1, ... , J – 1;

$p^{n + 1} = pro j_{P} (p^{n} + ρ (\nabla {\bar{f}}^{n} - {\bar{w}}^{n}))$ , $q^{n + 1} = pro j_{Q} (q^{n} + ρ ɛ ({\bar{w}}^{n}))$

f^old = fⁿ, $f^{n + 1} = {prox}_{1}^{τ} (f^{n} + τ di v_{1} p^{n + 1})$ , ${\bar{f}}^{n + 1} = 2 f^{n + 1} - f^{old}$ ,

w^old = wⁿ, wⁿ⁺¹ = wⁿ +τ(pⁿ⁺¹ + div₂ qⁿ⁺¹), ${\bar{w}}^{n + 1} = 2 w^{n + 1} - w^{old}$ ,

: End For;

${\bar{f}}^{0} = {\bar{f}}^{J}$ , w⁰ = w^J, ${\bar{w}}^{0} = {\bar{w}}^{J}$ , p⁰ = p^J, q⁰ = q^J

Using the OMP and SVD algorithms, update ${\bar{f}}^{k + 1}$ from images ${\bar{f}}^{k}$ .

If $f_{j}^{J} > 0$ ; then $f_{j}^{k + 1} = f_{j}^{k + 1}$ ; j = 1, 2, ... , N;

Else $f_{j}^{k + 1} = 0$ ; j = 1, 2, ... , N;

End If;

End if termination criterion is satisfied.

2.6 Parameter selections

2.6.1 Selection of β₁, β₂ and β₃

β₁, β₂ and β₃ are hyperparameters to balance the data fidelity and regularization terms. How to choose appropriate parameters is a tricky problem in the field of CT image reconstruction. Practically, selecting them was usually judged by experience, and then we compared the reconstructed image obtained from the normal dose images using the selected parameters via visual examination. In our work, we found that the reconstructed images were less sensitive to the value of β₁within the range 1×10^–3≤β₁≤1×10^–2. The two parameters β₂ and β₃ are especially important in controlling the smoothness of the reconstructed image; selection of these variables was optimized case by case on the basis of the noise level of the image and sparse-view projections.

2.6.2 Selection of parameters in solving (P2 and P3)

The iteration number of the sub-iteration step is an important factor for obtaining a successful result. Fewer iterations may fail to achieve the expected convergence effect, and more iterations may cause the related computational load to be very heavy. In practice, assuming that the trial solution of (P1) is a reasonable initial value for (P2) and that the trial solution of (P2) is a reasonable initial value for (P3), the total number of iterations could be greatly reduced. Furthermore, the step variables ρ and τ that control the step lengths of the updating procedure are also crucial; we must choose suitable values in the experiment. In our studies, ρ and τ were optimized by using α₀ = 2, α₁ = 1. The dictionary redundancy improves the sparsity of representation. Similarly, the parameters of dictionary learning should also be in an appropriate range. A patch size of 5×5 pixels was used in this paper. In the dictionary learning process, we solve the optimization problem (22) by minimizing the representation error with ∈= 2.5×10^–5.

3 Experiments and results

3.1 Experimental setup

To evaluate the performance of the PWLS-TGV-DL algorithm for CT image reconstruction, we conducted experiments on the digital extended cardiac-torso (XCAT) phantom and the head physical phantom.

3.1.1 Digital XCAT phantom

Figure 1(a) shows a slice of the XCAT phantom. We chose a geometry representative of a monoenergetic fan-beam CT scanner setup with a circular orbit to acquire 1160 projection views over 2π. The number of channels per view was 672. The distance from the detector arrays to the X-ray source was 1040 mm, and the distance from the rotation center to the X-ray source was 570 mm. The reconstructed images were composed of 512×512 square pixels. Each projection datum along an X-ray through the sectional image was calculated based on the known densities and intersection areas of the ray with the geometric shapes of the objects in the sectional image.

Fig.1

Digital and physical phantoms used in the studies: (a) a slice of the digital XCAT phantom; (b) the standard image reconstructed by the FBP method of the head phantom.

Similar to previous studies by Niu et al. [19], we first simulated the noise-free sonogram data ŷ and then generated the noisy transmission measurement I according to the statistical model of the prelogarithm projection data, that is, $y_{i} = Poisson (b_{i} exp (- \hat{y})) + Normal (0, σ_{e}^{2}),$ (23) where b_i is the incident X-ray intensity, and $σ_{e}^{2}$ is the background electronic noise variance. In the simulation, b_i and $σ_{e}^{2}$ were set to 1.0×10⁶ and 10.0, respectively, to simulate a low-dose scan. Finally, the noisy sinogram data y were calculated by performing the logarithm transformation on the transmission data y_i. For the digital XCAT phantom experiment, the sparse-view projections were generated by undersampling the 1,160 views of normal-dose simulation to only 116 views evenly over 2π.

3.1.2 Head physical phantom

Figure 1(b) shows the standard image of the head phantom reconstructed by the FBP method. This experiment was performed on an in-house CT imaging bench in our lab. The system had a rotating-anode tungsten target diagnostic level X-ray tube (Varex G-242, Varex Imaging Corporation, UT, USA) and was operated in 120.00 kV continuous fluoroscopy mode with a 0.40 mm nominal focal spot. The X-ray tube current was set to 11.00 mA. The X-ray detector was an energy-resolving photon-counting detector (XC-Hydra FX50, XCounter AB, Sweden) made from cadmium telluride (CdTe). The projection datasets were acquired with a 1×1 detector binning mode and were rebinned to 6×6 during the postprocessing procedures. Because there was no need to discriminate the photon energies, the detector was operated and calibrated with only a single energy threshold (10.00 keV). The source-to-detector distance was 1500.00 mm, and the source-to-rotation center was 1000.00 mm. A dental and diagnostic head phantom (Atom Max 711-HN, CIRS Inc., VA, USA) was imaged in this work. The phantom was rotated through 360 degrees at an 1.0 degree angular interval.

3.1.3 Performance evaluation

We applied numeric metrics to perform an objective assessment. The first metric is the structural similarity index (SSIM), which measures the similarity between two images by considering the luminance l, contrast c, and structural information s, where f is the reconstructed image, and y is the reference image. SSIM is defined mathematically as follows: ${\begin{matrix} l (f, y) = \frac{2 μ f μ_{y} + C_{1}}{μ_{f}^{2} + μ_{y}^{2} + C_{1}} \\ c (f, y) = \frac{2 σ_{f} σ_{y} + C_{2}}{σ_{f}^{2} + σ_{y}^{2} + C_{2}} \\ s (f, y) = \frac{σ_{fy} + C_{s}}{σ_{f} σ_{y} + C_{s}} \end{matrix}$ (24) $SSIM (f, y) = l (f, y) \cdot c (f, y) \cdot s (f, y),$ where μ _f and σ _f are the mean and variance of image f, respectively, μ _y and σ _y are the mean and variance of image y, respectively, σ _fy is the covariance of x and y, and C₁, C₂, and C₃ are constants. The SSIM index yields a value from 0 to 1; values closer to 1 denote greater image similarity.

The second metric is root mean square error (RMSE), which indicates the difference between the reconstructed image and the ground truth image and characterizes the reconstruction accuracy: $RMSE = \sqrt{\frac{\sum_{x = 1}^{X} {(y_{x} - f_{x}^{*})}^{2}}{X}},$ (25) where y_x is the true image, $f_{x}^{*}$ is the reconstructed image, and X is the number of image pixels.

To validate and evaluate the performance of PWLS-TGV-DL, TV regularization and TGV regularization were also conducted using the PWLS criterion for comparison. The methods being compared are referred to simply as PWLS-TV and PWLS-TGV. All the experiments were implemented using MATLAB 2016a and executed on a PC equipped with an Intel Core i7-6700 CPU @ 3.40 GHz and 56.0 GB of memory.

3.2 Performance evaluation on digital phantom

In the digital phantom study, the original phantom data were directly used as the ground-truth image. As mentioned above, low-dose CT can be implemented by limiting the number of projection views. To achieve a more comprehensive study, we tested our method using a few-view case.

In the few-view case, the CT scan projection views were set to 116 views over 360°. The reconstructed images and the zoomed ROIs of the tested methods are shown in Fig. 2. For simplicity, the two ROIs shown in the figure are referred to herein as ROI1 and ROI2, respectively. The PWLS-TV result in Fig. 2(a) contains substantial noise throughout the entire reconstructed region; some small structures are almost entirely obscured by noise, which can cause misdiagnosis in the clinic. Figure 2(b) shows the image reconstructed by PWLS-TGV, which is better than that of PWLS; the image is much cleaner after being denoised by TV, but it still contains some noise. Figure 2(c) shows the reconstruction by our proposed new model, which is to eliminate most noise and artifacts, indicating that PWLS-TGV-DL is effective at reducing noise. These images demonstrate that the proposed approach is better than the traditional CT reconstruction for few-view reconstruction.

Fig.2

Image reconstructed by (a) PWLS-TV, (b) PWLS-TGV, and (c) our proposed new model.

We also provide RMSE values and SSIM values in Table 1 for the reconstructed images of the three algorithms. The PWLS-TGV-DL method achieves the lowest RMSE value, and its SSIM value is closest to 1. Perhaps the RMSE and SSIM data from the full image are not significantly better. In order to better represent the advantages of our algorithm, we also provide the RMSE values and SSIM values from ROI1 and ROI2 in Table 2.

Table 1

Numeric results from three algorithms (PWLS-TV, PWLS-TGV and Proposed method)

Comparison	PWLS-TV	PWLS-TGV	PWLS-TGV-DL
RMSE	0.7566	0.6624	0.6021
SSIM	0.8964	0.9558	0.9856

Table 2

Numerical results of (ROI1 and ROI2) from three algorithms (PWLS-TV, PWLS-TGV and the proposed method)

Comparison	PWLS-TV	PWLS-TGV	PWLS-TGV-DL
RMSE (ROI1)	0.8725	0.6754	0.5206
SSIM (ROI1)	0.8424	0.9545	0.9912
RMSE (ROI2)	0.8975	0.6872	0.5564
SSIM (ROI2)	0.8342	0.9454	0.9884

The profile images and residual images are compared in Figs. 3 and 4, respectively. Clearly, the proposed method curve is closer to the Phantom curve than are those of the other two methods. The results show that compared to the other methods, our proposed method achieves images with superior quality.

Fig.3

The profiles were located at the pixel positions on the x-axis from 160 to 400 and y = 320. The Phantom curve represents the profile of the true image in Fig. 1(a). The PWLS-TGV curve represents the profile of the reconstructed image using the PWLS-TGV method in Fig. 2(b). The curve in Fig. 2(c) represents the profile of the reconstructed image using the proposed PWLS-TGV-DL method.

Fig.4

Residual images of the reconstructed results from the simulation XCAT data based on the (a) PWLS-TV method, (b) PWLS-TGV method, and (c) our proposed method. All the images are displayed in the same window.

3.3 Performance evaluation on physical phantom

In the physical phantom study, an image reconstructed with 720 views by the FBP method was used as the true image. To further evaluate our proposed method, we also performed an experiment in which the views were reduced to 90 and 180.

Figure 5 shows the images reconstructed by the different methods at 90 and 180 views. To further display the advantages of PWLS-TGV-DL, zoomed ROIs (as indicated by the squares in Fig. 5(a)) are shown in Fig. 6. Serious artifacts existed in the PWLS-TV results in all cases. Our proposed method yielded more noticeable gains than the PWLS-TGV method in terms of patchy artifact suppression. These results suggest that compared with the PWLS-TGV method, PWLS-TGV-DL achieves profiles that better match the gold standard.

Fig.5

Experimental results of the head phantom. The images in (a)– (c) are reconstructed with 90 views, and the images in (d)– (f) are reconstructed with 180 views. The images in (a) and (d) were reconstructed by the PWLS-TV algorithm; the images in (b) and (e) were reconstructed by the PWLS-TGV algorithm; and the images in (c) and (f) were reconstructed by our proposed algorithm.

Fig.6

Zoomed-in views of images reconstructed by PWLS-TV ((a) and (d)), PWLS-TGV((b) and (e)), and our proposed method ((c) and (f)) from 90-view (1st row) and 180-view (2nd row) projections.

Table 3 compares the numeric results of these methods. Note that PWLS-TGV-DL obtains the best results on all the numeric metrics. We can conclude that our proposed method has a clear performance advantage over the compared methods.

Table 3

Numeric results from three algorithms on head phantom

Comparison	PWLS-TV	PWLS-TGV	PWLS-TGV-DL
RMSE (90 views)	0.1392	0.0856	0.0262
SSIM (90 views)	0.9264	0.9493	0.9675
RMSE (180 views)	0.1049	0.0746	0.0190
SSIM (180 views)	0.9324	0.9604	0.9710

4 Discussion and conclusion

In this paper, we proposed a new few-view CT reconstruction solution by combining TV minimization and sparse dictionary learning based on the PWLS standard, and then we discussed the results in terms of image qualitative and quantitative studies such as SSIM and RMSE. Whether in terms of visual effects or performance evaluations, our proposed method has considerable advantages in terms of reconstructed image quality relative to traditional algorithms such as PWLS-TV and PWLS-TGV. The proposed method has a higher reconstruction accuracy, and it can suppress artifacts and noise while preserving more edge structure information.

In summary, we proposed a new algorithm that combines penalized weighted least-squares using total generalized variation and dictionary learning for few-view projection data and demonstrated its promising performance. First, an intermediate image is reconstructed using TGV minimization; then, it is postprocessed using dictionary learning to remove residual noise and produce a clinically acceptable CT image. The proposed PWLS-TGV-DL method is targeted at efficiently eliminating noise and reducing the artifacts of the TGV-based method for few-view CT image reconstruction. From the simulation and physical experimental results presented in Section 3, in which PWLS-TGV-DL was compared with reconstruction methods such as PWLS-TV and PWLS-TGV, the proposed method improves the quality of reconstructed images.

The better reconstruction performance of the proposed method opens a new idea to study how dictionary learning can improve the reconstructed image quality of traditional algorithms. The proposed method outperforms the conventional denoising scheme and the edge-preserving regularization in terms of image quality at a similar level of noise suppression. The approach is convenient and attractive for clinical applications. However, the main shortcoming of PWLS-TGV-DL is that the matrix update process during dictionary learning increases the computational burden and requires a long running time. To solve this problem, a fast computer and dedicated hardware are required. Future research work will include evaluations of the approach in more realistic situations or for clinical raw data reconstruction. In the future, we believe that iterative-based image reconstructions such as the PWLS-TGV-DL algorithm will be widely used in medical clinics.

Footnotes

Acknowledgments

The authors would like to thank Prof. Jianhua Ma and Zhaoying Bian at Southern Medical University for providing the total generalized variation code. This work was supported by the National Natural Science Foundation of China (81871441), Guangdong Special Support Program of China (2017TQ04R395), the Natural Science Foundation of Guangdong Province in China (2017A030313743), the Guangdong International Science and Technology Cooperation Project of China (2018A050506064), the National Natural Science Foundation of China (61601426). The authors would like to thank the editor and anonymous reviewers for their constructive comments and suggestions.

References

et al., Radiation dose reduction in computed tomography: Techniques and future perspective, Imaging Med 1(1) (2009), 65–84.

T. F.

et al., Nonlinear sinogram smoothing for low-dose X-ray CT, IEEE Transactions on Nuclear Science 51(5) (2004), 2505–2513.

Wang

, Li

, Xing

, Iterative image reconstruction for CBCT using edge-preserving prior, Medical Physics 36(1) (2009), 252–260.

Sidky

E.Y.

, Kao

C.M.

, Pan

X.H.

, Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT, J Xray Sci Technol 14(2) (2006), 119–139.

Sidky

E.Y.

, Pan

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

, Yang

, Jiang

, Wang

, Supplemental analysis on compressed sensing based interior tomography,N, Phys Med Biol 54(18) (2009), 425–32.

Xia

X.Q.

, Hu

X.D.

, Han

Z.Y.

, Zhang

, Xu

, Zou

, Scatter Correction in Cone Beam CT for Metal Additive Manufacturing Components, JOM 71(3) (2019), 1082–1087.

Kalra

M.K.

et al., Techniques and applications of automatic tube current modulation for CT, Radiology 233(3) (2004), 649–657.

McCollough

C.H.

, Bruesewitz

M.R.

, KoflerJ.M.> Jr, CT dose reduction and dose management tools: Overview of available options, Radiographics 26(2) (2006), 503–512.

10.

, Bruesewitz

M.R.

, Thomas

K.B.

, Fletcher

J.G.

, Kofler

J.M.

, McCollough

C.H.

, Optimal tube potential for radiation dose reduction in pediatric CT: Principles, Clinical Implementations, and Pitfalls, Radiographics 31(3) (2011), 835–848.

11.

Elbakri

I.A.

, Fessler

J.A.

, Statistical image reconstruction for polyenergetic X-ray computed tomography, IEEE Trans Med Imaging 21(2) (2002), 89–99.

12.

Lauzier

P.T.

, Chen

G.H.

, Characterization of statistical prior image constrained compressed sensing (PICCS): II. Application to dose reduction, Med Phys 40(2) (2013), 021902.

13.

, Feng

, Huang

, Chen

, Generalized Gibbs priors based positron emission tomography reconstruction, Comput Biol Med 40(6) (2010), 565–571.

14.

Liu

et al., Total variation-stokes strategy for sparse-view X-ray CT image reconstruction, IEEE Trans Med Imaging 33(3) (2014), 749–763.

15.

Wang

, Qi

, A new adaptive-weighted total variation sparse-view computed tomography image reconstruction with local improved gradient information, J Xray Sci Technol 26(6) (2018), 957–975.

16.

Panin

, Zeng

, Gullberg

, Total variation regulated EM algorithm, IEEE Nuclear Science Symposium 1998 3 (1998), 1562–1566.

17.

Rong

J.Y.

et al., CT reconstruction from few-views with anisotropic edge -guided total variance, (in English), Nuclear Instruments & Methods in Physics Research Section a-Accelerators Spectrometers Detectors and Associated Equipment 820 (2016), 54–64.

18.

Bredies

, Kunisch

, Pock

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

19.

Niu

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

20.

et al., Image reconstruction from few– view CT data by gradient– domain dictionary learning, J Xray Sci Technol 24(4) (2016), 627–638.

21.

, Yu

, Mou

, Zhang

, Hsieh

, Wang

, Low-dose X-ray CT reconstruction via dictionary learning, IEEE Trans Med Imaging 31(9) (2012), 1682–1697.

22.

Chen

et al., Artifact suppressed dictionary learning for low-dose CT image processing, IEEE Trans Med Imaging 33(12) (2014), 2271–2292.

23.

Aharon

, Elad

, Bruckstein

, K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation, IEEE Trans on Signal Processing 54(11) (2006), 4311–4322.

24.

Liu

, Chen

, Hu

, Luo

, Low-dose CBCT reconstruction via 3D dictionary learning, in pp, 2016 IEEE 13th International Symposium on (2016), 735–738.

25.

Z.L.

, Liang

, Xia

, Zheng

H.R.

, Compressive sampling in computed tomography: Method and application, Nuclear Instruments & Methods in Physics Research Section a-Accelerators Spectrometers Detectors and Associated Equipment 748 (2014), 26–32.

26.

Trinca

, Libin

, Performance of the sinogram-based iterative reconstruction in sparse view X-ray computed tomography, J Xray Sci Technol 27(1) (2019), 37–49.

27.

Singh

et al., Adaptive statistical iterative reconstruction technique for radiation dose reduction in chest CT: A pilot study, Radiology 259(2) (2011), 565–573.

28.

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

29.

Bian

et al., Evaluation of sparse-view reconstruction from flat-panel-detector cone-beam CT, Phys Med Biol 55(22) (2010), 6575–9659.

30.

Huang

, Ma

, Liu

, Feng

, Chen

, Projection data restoration guided non-local means for low-dose computed tomography reconstruction, in pp, 2011 IEEE International Symposium on (2011), 1167–1170.

31.

Han

, Bian

, Ritman

E.L.

, Sidky

E.Y.

, Pan

, Optimization-based reconstruction of sparse images from few-view projections, Phys Med Biol 57(16) (2012), 5245–5273.

32.

Ginat

D.T.

, Gupta

, Advances in computed tomography imaging technology, Annu Rev Biomed Eng 16 (2014), 431–453.

33.

Wang

, Li

, Lu

, Liang

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

34.

Liu

, Ma

, Fan

, Liang

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

35.

Mallat

, Zhang

, Matching pursuit with time-frequency dictionaries, Courant Institute of Mathematical Sciences New York United States 1993.

36.

Chen

, Billings

S.A.

, Luo

, Orthogonal least-squares methods and their application to non-linear system-identification, International Journal of Control 50(5) (1989), 1873–1896.

37.

Pati

Y.C.

, Rezaiifar

, Krishnaprasad

P.S.

, Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition, IEEE Proc of 27th Asilomar Conference on in Signals, Systems and Computers (1993), 40–44.

38.

Chen

S.S.B.

, Donoho

D.L.

, Saunders

M.A.

, Atomic decomposition by basis pursuit, (in English), Siam Review 43(1) (2001), 129–159.

39.

Bredies

, Dong

Y.Q.

, Hintermuller

, Spatially dependent regularization parameter selection in total generalized variation models for image restoration, (in English), International Journal of Computer Mathematics 90(1) (2013), 109–123.

40.

, Yu

, Wang

, Mou

Frontiers of Medical Imaging (2015), 99–119 Dictionary learning based low-dose X-ray CT reconstruction, World Scientific.

41.

Liao

H.Y.

, Sapiro

, Sparse representations for limited data tomography, in pp, 2008. ISBI 2008. 5th IEEE International Symposium on (2008), 1375–1378.

42.

Tropp

J.A.

, Greed is good: Algorithmic results for sparse approximation, IEEE Trans Information Theory 50(10) (2004), 2231–2242.

Few-view CT image reconstruction using improved total variation regularization

Abstract

Keywords

1 Introduction

2 Methods

2.1 PWLS image reconstruction

2.6.1 Selection of β1, β2 and β3

2.6.2 Selection of parameters in solving (P2 and P3)

3 Experiments and results

3.1 Experimental setup

3.1.1 Digital XCAT phantom

3.1.3 Performance evaluation

Footnotes

Acknowledgments

References

2.6.1 Selection of β₁, β₂ and β₃