Low Dose CT Image Reconstruction Based on Structure Tensor Total Variation Using Accelerated Fast Iterative Shrinkage Thresholding Algorithm

Abstract

Low dose computed tomography (CT) has drawn much attention in the medical imaging field because of its ability to reduce the radiation dose. Recently, statistical iterative reconstruction (SIR) with total variation (TV) penalty has been developed to low dose CT image reconstruction. Nevertheless, the TV penalty has the drawback of creating blocky effects in the reconstructed images. To overcome the limitations of TV, in this paper we firstly introduce the structure tensor total variation (STV₁) penalty into SIR framework for low dose CT image reconstruction. Then, an accelerated fast iterative shrinkage thresholding algorithm (AFISTA) is developed to minimize the objective function. The proposed AFISTA reconstruction algorithm was evaluated using numerical simulated low dose projection based on two CT images and realistic low dose projection data of a sheep lung CT perfusion. The experimental results demonstrated that our proposed STV₁-based algorithm outperform FBP and TV-based algorithm in terms of removing noise and restraining blocky effects.

1. Introduction

The X-ray computer tomography (CT) has been extensively utilized in industry nondestructive testing and medical diagnosis. However, repeated use of conventional CT could significantly increase the chance of developing cancer and other diseases due to the high radiation exposure [1,2,3]. Hence, low dose CT which was firstly proposed by Naidich [4], has become one of the research hot spots in the CT field. Lowering the tube current values (milliampere (mA) or milliampere second (mAs)) or the voltage values (kilovolt (kV)) is the most straightforward way to reduce the radiation dose because it does not need to change the scanning structure of existing CT equipment. Nonetheless, this method will result in insufficient number of x-ray photons detected at the detector and hence, upgrade the quantum noise level on the sinogram. In this situation, for most current commercial CT scanners, the often used Feldkamp-Davis-Kress algorithm (or its variants) will lead to severe image quality degradation due to noisy projection. Hence, it is highly desirable to develop a new method to reconstruct the high-quality image for LDCT imaging.

Due to its advantages in effective physical noise modeling and possibilities of incorporating priors of the reconstructed image, the statistical iterative reconstruction (SIR) method [5,6] had been shown to be superior in removing image noise and streak artifacts. Recently, inspired by compressed sensing (CS) theory [7], in 2008, Sidky et al. [8] first introduced the total variation (TV) into a low dose CT reconstruction field and proposed an adaptive-steepest-descent-projection onto convex sets (ASD-POCS) algorithm. Subsequently, Tang et al. [9] introduced the TV regularization term into SIR framework for low dose CT reconstruction and further developed a Gauss-Siedel iterative algorithm to minimize the objective function. Choi et al. [10] investigated a primal-dual first-order method to minimize the TV-based SIR framework for CBCT image reconstruction with sparse and noisy low-dose projection data. However, the TV regularization term has the drawback of creating staircase artifacts, particularly at low dose levels. To further eliminate the block effects, Xu et al. [11] proposed a dictionary learning-based approach for low dose CT reconstruction, in which the sparse constraint in terms of a redundant dictionary was incorporated into an objective function in a SIR reconstruction framework, and further proposed an alternating minimization reconstruction algorithm. Sun et al. [12] proposed a Hessian penalty and developed an effective algorithm to minimize the objective function using a maximization-minimization (MM) approach. Zhang et al. [13] described two iterative reconstruction algorithms for low dose CT image reconstruction based on a Gamma regularization. Shangguan et al. [14] introduced a modified Markov random field regularization term into the SIR framework and utilized a modified alternation iterative algorithm to optimize the cost function. Very recently, Shi et al. [15] combined the weighted TV and Hessian penalties for low dose CT reconstruction in a structure-adaptive way, in which a MM approach was designed to optimize the objection function. Xu et al. [16] introduced the TV and wavelet-based L₁ penalties into SIR framework and solved the objective function using accelerated variants of the fast iterative soft-thresholding algorithm. Zhang et al. [17] introduced an adaptive fractional order TV penalty and used a separable paraboloid surrogate (SPS) method to minimize the objection function. Liu et al. [18] discussed and compared the behaviors of several convex Hessian Schatten penalties with orders 1,2 and $\infty$ for low dose CT reconstruction. At the same time, K. Kim et al. [19] proposed low dose CT reconstruction using spatially encoded nonlocal penalty, in which an ordered subset SQS method for log-likelihood is developed and the patch-based similarity constraint with a spatially variant factor is developed to reduce the noise effectively and preserve features simultaneously. Very recently, Cai et al. [20] investigated block-marching sparsity regularization and developed a practical reconstruction algorithm using hard thresholding and projection onto convex set methods for low dose CT reconstruction.

In 2015, S. Lefkimmiatis et al. [21] proposed a novel structure tensor total variation (STV) penalty for image denoising and image deblurring. In contrast to TV penalty, STV penalized the image variation at every point of the domain by taking into account the information in a local neighborhood. Hence, it supplies a richer and more robust measure of image variation which translates to improve reconstruction performance. Inspired by this work, Zeng et al. [22] firstly introduced STV penalty into medical imaging field for multi-energy CT reconstruction. Later, Zeng et al. [23] developed a robust perfusion deconvolution approach via STV regularization for estimating an accurate residue function in cerebral perfusion CT with low mAs data acquisition. In 2018, Gong et al. [24] developed a SIR approach by incorporating a precontrast normal-dose CT scan of robust dynamic myocardial perfusion CT. Inspired by the above research, in this paper we introduce the STV penalty into the SIR framework for low dose CT reconstruction, and then develop an accelerated fast iterative shrinkage thresholding algorithm (AFISTA) to minimize the associated objective function. At last, the resulting method is evaluated using simulated low dose projection data and real CT projection data

The remainder of this paper is organized as follows. Section 2 describes the SIR, SIR framework with STV₁ penalty, and AFISTA for the proposed reconstruction model. Section 3 illustrates the performance of the proposed reconstruction approach via two numerical simulated experiments and realistic CT projection. Finally, we will discuss relevant issues and conclude this paper in Section 4.

2. Materials and Methods

2.1. SIR

The SIR algorithm [5] for reconstructing a CT image $u$ can be expressed as the following minimization problem:

$$u^{\ast }=\underset{u}{\arg \min }\left\{\Phi (u)\triangleq \frac{1}{2}{\Vert Au-p\Vert }_W^2+\beta R(u)\right\}$$

(1)

where $u={\left(u_1,u_2,\dots ,u_N\right)}^T$ is the image vector representing the attenuation coefficients of the imaged object, $p={\left(p_1,p_2,\dots ,p_m\right)}^T$ is the vector representing the raw detector measurements, $A$ is the $m\times N$ system matrix that models the imaging system, the diagonal matrix $W$ provides statistical weighting that accounts for the ray-dependent variance of the noise, $R(u)$ is the function of the image $u$ which is called regularization term, $\beta$ is the regularization parameter that balances the data-fitting term and the regularization term, and $\Phi (u)$ is the objective function.

2.2. SIR with STV₁ Penalty for Low Dose CT Reconstruction

Discrete STV₁ penalty [21] can be expressed as:

$${\mathrm{STV}}_1(u)={\displaystyle \sum _{n=1}^N{\Vert {\left[J_{K}u\right]}_n\Vert }_{S_1}}$$

(2)

where $u$ denotes an image, N denotes the number of pixels in the image $u$, ${\Vert \cdot \Vert }_{S_1}$ denotes the Schatten nuclear norm, corresponding to the ${\mathrm{S}}_1$ norm of the singular values, $J_K$ is the patch-based Jacobian operator of the image $u$, which is defined as the linear mapping $J_K:{\mathbb{R}}^N\mapsto 𝒳$, where $𝒳\triangleq {\mathbb{R}}^{N\times L\times 2}$, $K$ is a non-negative, rotationally symmetric convolution kernel, ${[J_{K}u]}_n$ denotes a matrix with the patch-based Jacobian applied on the nth pixel of the $u$, which is defined as

$${[J_{K}u]}_n={\left({\left[T_{s_1,w}\circ \nabla u\right]}_n,\dots ,{\left[T_{s_L,w}\circ \nabla u\right]}_n\right)}^T$$

(3)

where ${(\cdot )}^T$ is the transpose operator, $\nabla$ denotes the discrete gradient operator, $\circ$ denotes the composition of operator. The shift vectors $s_l(l=1,\dots ,L)$ are the elements in the neighborhood of the nth pixel, where $L={(2L_k+1)}^2$, $L_k$ is the radius of the convolution kernel, $T_{S_l,w}$ is a weighted translation operator, which is defined as:

$${\left[T_{{}_{S_l,w}}\circ \nabla u\right]}_n=w(s_l)\nabla u[x_n-s_l]$$

(4)

where $x_n$ represents the coordinates of image $u$ and $w(s)=\sqrt{K[s]}$ denotes the window function of the convolution kernel [21]. In this work, we chose Gaussian kernel to design the structure tensor.

To overcome the block effects of TV penalty and inspired by the successful application of STV₁ penalty in image processing field [21], a new SIR-STV₁ framework for low dose CT reconstruction is proposed as follows:

$$u^{\ast }=\arg \underset{u}{\min }\left\{\Phi (u)\triangleq \frac{1}{2}{\Vert Au-p\Vert }_W^2+\beta {\mathrm{STV}}_1(u)\right\}$$

(5)

2.3. AFISTA Algorithm for Solving SIR-STV₁

2.3.1. General FISTA for Solving SIR-STV₁

FISTA has been widely used for image denoising and deblurring in image processing field due to its ability to minimize a cost function that is specified by the sum of a smooth and convex data fidelity term and a convex but possibly nonsmooth penalty [25]. Since the STV₁ penalty is not differentiable everywhere, the traditional gradient-based algorithm cannot be directly applied to minimize the cost function (5). Therefore, in this paper we propose to adapt the FISTA algorithm to solve problem (5).

Let

$$\mathrm{F}(u)=\frac{1}{2}{\Vert Au-p\Vert }_w^2$$

(6)

the FISTA algorithm decoupled the minimization problem (5) into two steps. The first step is to minimize the data fitting term $\mathrm{F}(y^{(i)})$ using a gradient descent method, which can be expressed as:

$$v^{(i)}=y^{(i)}-\frac{1}{Q}\nabla \mathrm{F}(y^{(i)})=y^{(i)}-\frac{1}{Q}{A}^{T}W(A{y}^{(i)}-p)$$

(7)

where i is the iteration number, Q is the Lipschitz constant of $\nabla F(u)$. The second step is STV₁-based denoising problem, which can be represented as:

$$u^{(i)}=\underset{u}{\arg \min }\left\{\frac{1}{2}{\Vert u-v^{(i)}\Vert }_2^2+\frac{\beta }{Q}{\mathrm{STV}}_1(u)\right\}$$

(8)

Let $\lambda =\beta /Q$, Equation (8) can be solved efficiently by a progressive proximal map algorithm [21]. In all, the detailed workflow of the general FISTA is listed in Algorithm 1.

Algorithm 1 Workflow of the general fast iterative shrinkage thresholding algorithm (FISTA).

Input: system matrix $A$, projection data $p$, Q is the Lipschitz constant of $\nabla F(u)$

Initial Step:$y^{(1)}=u^{(0)}=0$, $t_1=1$, maximum iteration number I, regularization parameter $\lambda$, convolution kernel $K$

fori = 1,2,...,I

Update intermediate image $v^{(i)}$ with (7)

Update image $u^{(i)}$ with (8)

$t_{i+1}=\frac{1}{2}(1+\sqrt{1+4t_i^2})$

$y^{(i+1)}=u^{(i)}+((t_i-1)/t_{i+1})/(u^{(i}{}^{)}-u^{(i-1)})$

end for

output:$u^{(I)}$

2.3.2. AFISTA for Solving SIR-STV₁

Although the general FISTA algorithm can be applicable to solve problem (5), it poses challenges for optimization. Lipschitz constant of the gradient data-fitting term is very large for low dose CT reconstruction [6,26,27], which results in small gradient steps leading to slow convergence. To solve this problem, we proposed to use the separable quadratic surrogate (SQS) method [6,28], that replaces the data fitting term $\mathrm{F}(y^{(i)})$ by a surrogate function $\varphi (v;y^{(i)})$ at ith iteration.

$$\varphi (v;y^{(i)})\triangleq F(y^{(i)})+\nabla F(y^{(i)})(v-y^{(i)})+\frac{D}{2}{(v-y^{(i)})}^2$$

(9)

where $D\triangleq diag\{[A^{T}WA]1\}$ is a diagonal Hessian (second order derivatives) matrix. $\varphi (v;y^{(i)})$ can be easily minimized by zeroing the first derivative, this leads to the following updating algorithm:

$$v^{(i+1)}=y^{(i)}-D^{-1}\nabla \mathrm{F}(y^{(i)})=y^{(i)}-D^{-1}{A}^{T}W(A{y}^{(i)}-p)$$

(10)

To further accelerate updating Equation (10), we adopt order subset (OS) methods [29], by grouping the projection views into M subsets evenly and using only the subset of measured data to approximate the exact gradient of the cost function.

Furthermore, to find the solution of Equation (8), we first introduce a dual norm of ${\mathrm{STV}}_1(u)$ and write it as follows [21]:

$${\displaystyle \sum _{n=1}^N{{\displaystyle \Vert {[J_{K}u]}_n\Vert }}_{S_1}}=\underset{\Omega \in B_{\infty ,\infty }}{\max }{\langle \Omega ,J_{K}u\rangle }_{𝒳}$$

(11)

where $𝒳\triangleq {ℝ}^{N\times L\times 2}$ is the target space of $J_K$, $\Omega$ denotes the variable in the target space $𝒳$, $B_{\infty ,\infty }$ denotes the $l_{\infty }-S_{\infty }$ unit-norm ball, which can be expressed as:

$$B_{\infty ,\infty }=\{\Omega \in 𝒳:\Vert {\Omega }_n\Vert {}_{S_n}\le 1,\forall n=1,\dots ,N\}$$

(12)

Then, the proximal map of STV₁ can be rewritten as follows:

$$u^{\ast }=\underset{u\in C}{\arg \min }\frac{1}{2}{\Vert u-v\Vert }_2^2+\lambda \underset{\Omega \in B_{\infty ,\infty }}{\max }{\langle J_K^{*}\Omega ,u\rangle }_2$$

(13)

where convex set $C={ℝ}^N$, ${\langle J_K^{*}\Omega ,u\rangle }_2={\langle \Omega ,J_{K}u\rangle }_{𝒳}$ and $J_K^{*}$ denotes the adjoint of the patch-based Jacobian operator. This formulation naturally leads us to the following minmax problem:

$$\underset{u\in C}{\min }\underset{\Omega \in B_{\infty ,\infty }}{\max }ℒ(u,\Omega )=\frac{1}{2}{\Vert u-v\Vert }_2^2+\lambda {\langle J_K^{*}\Omega ,u\rangle }_2$$

(14)

Since $ℒ$ is a strictly convex w.r.t. $u$ and concave w.r.t. $\Omega$, a saddle point of $ℒ$ can be obtained. Therefore, the order of the minimum and the maximum in Equation (13) does not affect the solution. This means that there exists a common saddle point $(u^{\ast },{\Omega }^{\ast })$ when the minimum and the maximum are interchanged, i.e.,

$$\underset{u\in C}{\min }\underset{\Omega \in B_{\infty ,\infty }}{\max }ℒ(u,\Omega )=ℒ(u^{\ast },{\Omega }^{\ast })=\underset{\Omega \in B_{\infty ,\infty }}{\max }\underset{u\in C}{\min }ℒ(u,\Omega )$$

(15)

Based on Equation (15), two optimization problems, the primal and the dual one can be defined. This can be accomplished by identifying the primal and dual objective functions, respectively, with the following minmax problem [21]:

$$p(u)=\underset{\Omega \in B_{\infty ,\infty }}{\max }ℒ(u,\Omega )=\frac{1}{2}{\Vert u-v\Vert }_2^2+\lambda {\displaystyle \sum _{n=1}^N{\Vert {[J_{K}u]}_n\Vert }_{S_1}}$$

(16)

$$d(\Omega )=\underset{u\in C}{\min }ℒ(u,\Omega )=\frac{1}{2}{\Vert z-{\displaystyle {\prod }_C(z)}\Vert }_2^2+\frac{1}{2}\left({\Vert v\Vert }_2^2-{\Vert z\Vert }_2^2\right)$$

(17)

where $z=v-\lambda J_K^{*}\Omega$ and ${\prod }_C$ is the orthogonal projection operator on the convex set C. According to the conclusion in Ref. 21, the minimizer $u^{\ast }$ of the primal objective function can be obtained from the maximizer ${\Omega }^{\ast }$ of the dual objective function. It can be expressed as follows:

$$u^{\ast }={\prod }_C(v-\lambda J_K^{\ast }{\Omega }^{\ast })$$

(18)

where ${\Omega }^{\ast }$ can be derived as the optimization of the corresponding dual objective function in Equation (17). For minimizing the dual objective function in Equation (17), the progressive proximal map algorithm [21] can be used and listed as follows (Algorithm 2):

Algorithm 2 Workflow of the progressive proximal map algorithm

Input:$v$, $\lambda$, ${\prod }_C$

Initial Step:${\Psi }^{\left(1\right)}={\Omega }^{(0)}=0$, $t_1=1$

fori = 1,2,...,I

${\Omega }^{(i)}={\prod }_{B\infty ,\infty }({\Psi }^{(i)}+\frac{1}{8\sqrt{2}}{J}_K{\prod }_C(v-\lambda J_{{}_K}^{\ast }{\Psi }^{(i)}))$

$t_{i+1}=\frac{1}{2}(1+\sqrt{1+4t_i^2})$

${\Psi }^{(i+1)}={\Omega }^{(i)}+((t_i-1)/t_{i+1})/({\Omega }^{(i)}-{\Omega }^{(i-1)})$

end for

Output$u={\prod }_C(v-\lambda J_K^{\ast }{\Omega }^{(I-1)})$

In all, the proposed AFISTA can be listed as follows (Algorithm 3):

Algorithm 3 Workflow of the proposed accelerated (A)FISTA

Input: System matrix $A$, projection data $p$, $v$, $\lambda$, ${\prod }_C$

Initial Step:$y^{\left(1\right)}=u^{(0)}=0$, $t_1=1$, maximum iteration number I, regularization parameter $\lambda$, convolution kernel K, the number of order subset M

fori = 1,2...I

for m = 1,...M

$k=\left(i-1\right)\times I+m$

Update intermediate image using Equation (10):

$v^{(k)}=y^{(k)}-D^{-1}M{A}_m^T{W}_m(A_m{y}^{\left(k\right)}-p_m)$

$A_m$, $W_m$, $p_m$ are submatrices of $A$, $W$, $p$ corresponding to the mth subset

Update image $u^{(k)}$ with progressive proximal map algorithm using Algorithm 2

$t_{k+1}=\frac{1}{2}(1+\sqrt{1+4t_k^2})$

$y^{(k+1)}=u^{(k)}+((t_k-1)/t_{k+1})/(u^{(k}{}^{)}-u^{(k-1)})$

end for

end for

Output$u^{(I\times M)}$

3. Experimental Results

In this subsection, a series of numerical simulation experiments were designed to evaluate the performance of the proposed method in CT image reconstruction from a low dose situation. In addition, low dose and under sampled raw projections of a sheep lung perfusion were acquired to validate the feasibility of our proposed method in practical applications. Meanwhile, the TV-based SIR (denoted by SIR-TV) and filtered backprojection(FBP) methods were presented for comparison.

3.1. Brain Image Numerical Simulation

A human brain CT image, downloaded from website [30], was used to validate the performance of the proposed algorithm. This image was a 2 mm thick corrected form of a female patient. During scanning, the X-ray tube voltage and tube current was 120 KVp and 200 mAs, respectively, a monochromatic spectrum was adopted, scattered photons were not considered, and the equi-angularly detector was used. The image is shown in Figure 1a. The spatial dimension is 512 × 512 pixels and the size of each pixel is 0.4883 × 0.4883 mm.

In this simulation, we adopt the fan-beam geometry with an equi-angular detector to simulate projection. The distance from the X-ray source to the detector is 1140 mm and the distance from the rotation center to the curved detector is 570 mm. Uniformly collected in [0, 2π] were 1160 projections. For each projection, 672 detector elements spanned a field-of-view (FOV) of 25 cm in radius. To obtain noisy projection, we firstly compute noiseless projection using siddons’s ray-driven algorithm. Then, Poisson noise assuming 5 × 10³, 1 × 10⁴, and 5 × 10⁴ photons per detector element was superimposed on the obtained noiseless projection to simulated three different noise levels, respectively.

3.1.1. Convergence Analysis

To examine the convergence of our proposed algorithm, Figure 2 shows the value of objective function in terms of iteration steps for the brain image numerical simulation with 5 × 10⁴ incident photon numbers. The curve indicates that our proposed algorithm can converge to a steady solution after 20 iterations. Furthermore, it can be seen that the proposed AFISTA method converges much faster than FISTA.

3.1.2. Visual Quality Comparison

Figure 3 shows the reconstructed brain images and zoomed regions of interest(ROI) (indicated by the red square in Figure 1a) reconstructed from low dose projections. Images in the left, middle, and right column are reconstructed by the FBP, SIR-TV, and SIR-STV₁ method, respectively. Images in the first, second, and third row are reconstructed from noisy projections with incident photon numbers 5 × 10³, 1 × 10⁴, and 5 × 10⁴, respectively. In the SIR-STV₁ method, parameters $\sigma$ and $L_K$ are selected according to the method described in Ref. [21], i.e.,$\sigma =0.5$, $L_K=3$, the number of subset were set to 40. For the reconstruction from noisy projections with incident photon numbers 5 × 10³, 1 × 10⁴, and 5 × 10⁴, the parameter $\lambda$ were set as 6 × 10⁻⁵, 2 × 10⁻⁵, and 9 × 10⁻⁶, respectively. It can be seen that the FBP results contains serious noise, and became worse and worse when the number of photon decreased from 5 × 10⁴ to 5 × 10³. The SIR-TV method removes noise but the reconstructed image suffers from over smoothing effect. The reconstructed results using the SIR-STV₁ perform better in restraining the blocky effect and removing image noise than other methods, and are almost visually the same as the original image.

3.1.3. Quantitative Comparison

To quantify the reconstruction accuracy of the low dose reconstruction algorithm, signal-to-noise ratio (PSNR) [31], relative reconstruction error (RRE), and structure similarity (SSIM) [32] indexes are employed in this subsection.

$$\mathrm{PSNR}=10{\log }_{10}\frac{{(\max ({\mu }_n))}^2}{{\displaystyle \sum _n{({\mu }_n-{\mu }_n^{\ast })}^2/N}}$$

(19)

$$\mathrm{RRE}=\frac{{\Vert \mu -{\mu }^{\ast }\Vert }_2^2}{{\Vert \mu \Vert }_2^2}$$

(20)

$$\mathrm{SSIM}=\frac{(2\overline{\mu }{\overline{\mu }}^{\ast }+c_1)(2{\sigma }_{\mu {\mu }^{\ast }}+c_2)}{({\overline{\mu }}^2-{\overline{\mu }}^{\ast }{}^2+c_1)({\sigma }_{\mu }^2+{\sigma }_{{\mu }^{\ast }}^2+c_2)}$$

(21)

where ${\mu }_n$ is the reconstructed value, ${\mu }_n^{\ast }$ is the golden reference value,$J$ is the number of the pixels in the reconstructed image, $\overline{\mu }$ and ${\overline{\mu }}^{\ast }$ are the mean value of $\mu$ and ${\mu }_{}^{\ast }$, ${\sigma }_{\mu }^{}$ and ${\sigma }_{{\mu }^{\ast }}^{}$ are the variation of $\mu$ and ${\mu }_{}^{\ast }$, ${\sigma }_{\mu {\mu }^{\ast }}$ is the covariance of $\mu$ and ${\mu }_{}^{\ast }$, $c_1$ and $c_1$ are the constants that we choose according to Reference [32].

Table 1. Image quality metrics.

Incident Photon Number	Algorithm	PSNR (dB)	RRE	SSIM
5 × 103	FBP	26.6131	0.2049	0.5526
	SIR-TV	33.8416	0.0892	0.7539
	SIR-STV1	35.7597	0.0715	0.7639
1 × 104	FBP	29.8176	0.1417	0.6135
	SIR-TV	35.7454	0.0716	0.8162
	SIR-STV1	38.2683	0.0536	0.8248
5 × 104	FBP	36.7103	0.0641	0.7318
	SIR-TV	40.2409	0.0427	0.9016
	SIR-STV1	42.8079	0.0307	0.9275

gives the PSNR, RRE, and SSIM values of the reconstructed whole brain images in Figure 3. It can be seen that the SIR-STV₁ method has the best performance in all evaluation metrics in all noise levels.

3.2. Thorax Image Numerical Simulation

To further validate the performance of the SIR-STV₁ method, a human Thorax image is downloaded from website [33], which is shown in Figure 1b. The spatial dimensional of the image is 512 × 512 pixels, the real size of each pixel and the imaging geometry are exactly the same as that used in the brain numerical simulation. The Poisson noise assuming 5 × 10⁴ photons per detector element was added to simulate noisy projection.

3.2.1. Visual Quality Comparison

Figure 4 shows the reconstructed images from noisy projections with incident photon numbers 5 × 10⁴ by the FBP, SIR-TV, and SIR- STV₁ methods. In the proposed SIR-STV₁ method,$\sigma$ is set to 0.5, $L_K$ is set to be 3, and $\lambda$ is set to $2\times {10}^{-6}$. It can be seen that SIR-TV and SIR-STV₁ can remove noise effectively and the SIR-STV₁ method performs better in eliminating blocky effect than SIR-TV.

3.2.2. Quantitative Comparison

To quantitatively evaluate the performance of the SIR-STV₁ method, we compare the performance of the three methods on the reconstruction of ROIs with detailed structures, which were marked by red rectangles in Figure 1b. The quantitative results are given in

Table 2. A summary of the evaluation indexes of the reconstructed image at different noise levels in the Thorax image simulated studies.

Different ROI	Algorithm	PSNR (dB)	RRE	SSIM
ROI 1	FBP	30.7537	0. 0447	0.7403
	SIR-TV	35.0948	0.0271	0.9206
	SIR-STV1	35.7648	0.0251	0.9370
ROI 2	FBP	29.9316	0.0540	0.7730
	SIR-TV	33.1043	0.0374	0.9145
	SIR-STV1	34.2960	0.0326	0.9348
ROI 3	FBP	30.0721	0.0446	0.7021
	SIR-TV	33.8564	0.0288	0.8927
	SIR-STV1	34.9798	0.0253	0.9183

. It can be seen that the SIR- STV₁ method has the lowest RRE and highest SSIM for all of the ROIs.

3.2.3. Profile-Based Comparison

In order to further visualize the differences among the three methods in this Thorax numerical simulation experiment, the profiles of the resulting images corresponding to the white line in Figure 1b were plotted in Figure 5. It can be demonstrated that the profiles of the Thorax images reconstructed by the SIR-STV₁ method achieve the best performance in terms of noise suppression and fine structure preservation.

3.2.4. Analysis of the Parameter

(1) Regularization parameter $\lambda$: To sense the sensitivity of $\lambda$, experiments are performed with various $\lambda =5\times {10}^{-7}$, $1\times {10}^{-6}$, $2\times {10}^{-6}$, $3\times {10}^{-6}$, $4\times {10}^{-6}$, and $5\times {10}^{-6}$ with $\sigma =0.5$, $L_K=3$. The reconstruction images are shown in Figure 6. The PSNR, RRE, and SSIM curves corresponding to different $\lambda$ are plotted in Figure 7. It is clear that the value of PSNR and SSIM are highest and the value of RRE is lowest when $\lambda =2\times {10}^{-6}$.

(2) Parameters $\sigma$ and L_K

When $\sigma$ was set to 0.5, Gaussian kernel with 99.7% of its energy will be within three pixels, so L_K and σ should be changed together instead of separately. To sense the impact of parameter $\sigma$ and L_K, six different sets of parameters including $\sigma =0.5,L_K=3$, $\sigma =0.8,L_K=5$, $\sigma =1.2,L_K=7$, $\sigma =1.5$, $L_K=9$ $\sigma =1.8,$ $L_K=11$, and $\sigma =2.1,L_K=13$ are tested. The results are given in Figure 8. In addition, the quantitative results are given in

Table 3. A summary of the evaluation indexes of the reconstructed image in Figure 8.

Parameter	PSNR (dB)	RRE	SSIM
$\sigma =0.5$, LK = 3	37.7983	0.0290	0.9449
$\sigma =0.8$, LK = 5	37.9762	0.0283	0.9477
$\sigma =1.2$, LK = 7	37.9653	0.0284	0.9482
$\sigma =1.5$, LK = 9	37.8765	0.0287	0.9480
$\sigma =1.8$, LK = 11	37.7748	0.0291	0.9474
$\sigma =2.1$, LK = 13	37.6704	0.0294	0.9466

. The results show that the performance of our proposed method is not quite sensitive to parameter L_K and σ and that there are no noticeable difference for different parameters. To balance the performance and computational time, the $\sigma$ is set to 0.5, L_K is set to 3 in this paper.

(3) Number of projections:

To evaluate the impact of the number of projections, simulated experiments with 580, 290, 145, and 116 views were performed. For all experiments, $\sigma$ were set to 0.5, $L_K$ were set to be 3. When the number of projection view decreases, the regularization parameter $\lambda$ should be increased accordingly to obtain high quality reconstructed results. For the reconstructions from 580, 290, 145, and 116 views, $\lambda$ were set to 3 × 10⁻⁶, 4 × 10⁻⁶, 6 × 10⁻⁶, and 8 × 10⁻⁶, the number of subset were set to 20, 10, 5, and 4, respectively. The reconstruction results are shown in Figure 9. It can be observed that our proposed method can suppress the noise effectively. In the case of 580 views, the reconstruction result was almost as good as that reconstructed from 1160 views in Figure 4c. In the cases of 290, 145, and 116 views, our proposed method still obtains high quality reconstruction results.

3.3. Realistic Sheep Lung Experiments

To validate the effectiveness of our proposed method for real data, an anesthetized sheep lung was scanned at normal and low dose, respectively on a SIEMENS Somatom Sensation 64-Slice CT Scanner in a circular cone-beam scanning mode. A scan protocol was developed for low dose studies with ECG gating: Time point 1 for a normal X-ray dose scan (100 kV/150 mAs) before a contrast agent injection, and time points 2–21 for low dose scans (80 kV/17 mAs) after the contrast agent injection. All the sinograms of the central slice were extracted, which were in a fan-beam geometry. The radius of the trajectory was 57 cm. A total of 1160 projections were uniformly collected over a 360° range. For each projection, 672 detector elements were equi-angularly distributed to define a FOV of 25.05 cm in radius. In this experiment, the reconstructed images were 768 × 768 pixels with a physical size of 50 × 50 cm. To further demonstrate the performance of our proposed method for a few views of CT reconstruction, 580 and 290 views projection were uniformly exacted form 1160 views projection. In our proposed algorithm, $\sigma =0.5$, $L_K$ = 3. For the reconstruction from noisy projections with 1160, 580, and 290 views, the parameter $\lambda$ were set as 5 × 10⁻⁶, 7 × 10⁻⁶, and 2 × 10⁻⁵, the number of subset were set to 40, 20, and 10, respectively.

The reconstruction results are presented in Figure 10. It can be seen that the FBP results look very noisy and became worse and worse when the number of projection views decreased from 1160 to 290. Both the SIR-TV and our SIR-STV₁ method outperform the FBP algorithm in suppressing image noise. To clearly compare the reconstruction performance of all algorithms, ROIs are extracted from Figure 10 and magnified in Figure 11. From Figure 11, especially denoted by red arrow regions, we can observe that the SIR-TV method produces patchy artifacts obviously, while the SIR-STV₁ method could avoid the patchy artifacts effectively.

4. Discussion and Conclusions

With the development of STV₁ in recent years, STV₁ has been applied in many medical imaging problems, including multi-energy CT, dynamic myocardial perfusion CT, etc. Instead of penalizing the image at every pixel, STV₁ considers the available information from the neighborhood of every pixel by penalizing the eigenvalue of structure tensor at this point. In this paper, we introduced the STV₁ penalty into the SIR reconstruction framework for low dose CT reconstruction. Comparison studies among FBP and the SIR-TV method revealed that the proposed STV₁ penalty effectively removed image noise and restrained the staircase effect of the TV penalty.

There are three parameters, the standard deviation of Gaussian convolution kernel $\sigma$, the radius of Gaussian convolution kernel L_K, and the regularization parameter $\lambda$ that needs to be adjusted. Extensive experiments in Section 3.2.4 reveal that the performance of our proposed method is not quite sensitive to the parameters L_K and σ. Therefore, to balance the reconstruction performance and computational complexity, we chose σ = 0.5 and L_K = 3 in our paper. In addition, from our experience this choice is also suitable for a wide range of applications. For the selection of $\lambda$, a large number of $\lambda$ was manually tried to yield an optimal one in simulated experiments and realistic sheep lung data studies on the condition with fixed L_K and σ. It is known that this will require a large computation time and cost. In the future, we would study a regularization parameter selection method to obtain an approximate optimal regularization parameter for the presented AFISTA.

The OS-SQS method is massively parallelizable and well suited to modern computing architecture such as graphics processing unit (GPU), which is applied to replace the gradient descent method for minimizing the data fitting term in the original FISTA. Experimental results in Section 3.1.1 have demonstrated that the proposed AFISTA method converges much faster than FISTA. However, the convergence behavior of the AFISTA may be affected by the number and ordering of subsets, that is the proposed algorithms in some cases become unstable when it is too large. In the future, we would study how to balance the number of subsets and the convergence of the AFISTA. Meanwhile, stochastic gradient methods are investigated to realize the random ordering of the subsets.

Another important issue is the computation time. MATLAB codes are developed on a PC with a 3.20 GHz Intel core i7 processor and a 16 GB RAM. Our proposed AFISTA algorithm includes an image updating step and STV₁-based denoising step. In brain image numerical simulation with 1160 views, the image updating step and STV₁-based denoising step in each iteration takes 2.3463 s (on average) and 40 s (on average), respectively. GPU was employed in the projection and back-projection operation to speed up the proposed algorithm. In order to further improve the reconstruction speed of the AFISTA, we would use GPU in STV₁-based denoising step.

In conclusion, we introduced the STV₁ regularizer into SIR framework for low dose CT reconstruction and developed an effective algorithm to minimize the objective function using AFISTA. Both simulated data and realistic sheep data were used to validate the proposed method. From the experiments, we have seen that the proposed SIR-STV₁ have better performance in restraining the stair effect and noise suppression than FBP and SIR-TV method. Additionally, in the future, the proposed SIR-STV₁ framework will be refined and adapted to handle other topics in CT imaging, such as few-view reconstruction and interior CT.