A new adaptive-weighted total variation sparse-view computed tomography image reconstruction with local improved gradient information

Abstract

Inspired by the compressed sensing (CS) theory, introducing priori information of sparse image into sparse-view reconstruction algorithm of computed tomography (CT) can improve image quality. In recent years, as a special case of CS, total variation (TV) reconstruction algorithm that uses both image sparsity and prior information of edge direction have attracted much research interest in sparse-view image reconstruction due to its ability to preserve image edges. In this paper, we propose a new adaptive-weighted total variation (NAWTV) algorithm for CT image reconstruction, which is derived by considering local gradient direction continuity and the anisotropic edge property. The anisotropic edge property is used to consolidate the image sparsity, where the associated weights are expressed as a combination of exponential and cosine function. The weights can also be adjusted adaptively according the local image intensity gradient. The NAWTV algorithm is numerically implemented with gradient descent method. The typical Shepp-Logan phantom and FORBILD head phantom are employed to perform image reconstruction simulation. To evaluate performance of NAWTV algorithm, we compared it with TV and AwTV reconstruction algorithms in experiments. Numerical experimental results verified the effectiveness and feasibility of the proposed algorithm. Comparison results also showed that the NAWTV algorithm achieved a satisfactory performance in suppressing artifacts and preserving the edge structure details information of the reconstructed image.

Keywords

Compressed sensing computed tomography sparse-view reconstruction total variation

1 Introduction

X-ray computed tomography (CT) is one of the widely used technologies in the field of medical imaging. However, the extra radiation dose in clinical examination can increase the health risk, which may affect the physiological function of human body and increase the chance of carcinogenesis, especially for children [1]. Therefore, it is very important to control radiation dose in clinical examination. It is well known that both lowering the x-ray intensity and reducing the number of projections when the patient is scanned to get the projection data are effective ways to reduce radiation dose [2 –4]. When the intensity of the X-ray is reduced by adjusting the parameters of the equipment, the acquisition projection data will unavoidably increase the unstable noise, which causes the reconstructed image quality to fall [5]. The acquisition of reducing the number of projections is a problem of the incomplete projection data reconstruction of the under-sampled, which is an ill-posed inverse problem in mathematical theory [6].

When the analytic algorithms such as filtered back-projection (FBP) [7] is used for under-sampled sparse projection reconstruction, image quality can be greatly affected by the aliasing artifacts due to insufficient projection information [6, 8]. However, the application of traditional iterative algorithms such as algebraic reconstruction technique (ART) [9] and the maximum likelihood expectation maximization (MLEM) [10] will not be disturbed by aliasing artifacts. The iterative algorithm can reconstruct a clearer image, but the reconstruction image will also show some unnecessary artifacts which can’t meet the clinical needs [11]. In recent years, the development of compressed sensing (CS) theory [12 –14] has provided new ideas for solving this problem. The key point of CS is that CT images can be accurately reconstructed by mining its L1 norm or the L1 norm of its sparse transformation from under-sampled projection if the original image is sparse or can find a suitable transform domain to make it sparse [15]. The sparse transform commonly used in the image domain mainly includes discrete gradient transform [12], wavelet transform [16], curvelet transform [17] and shearlet transform [18] et al. In CT reconstruction field, the iterative reconstruction algorithm based on CS employs the total variation (TV, i.e. the L1 norm of image gradient) [19] as regularization item in constructing the optimization problems. It has been widely studied in the reconstruction from under-sampled projection data [20 –22].

In 2006, Sidky et al. [23] proposed the TV-POCS (projection onto convex sets) reconstruction algorithm, which shows that this kind of algorithm can reconstruct the CT image in the sparse-view projection data. However, because of the piecewise constant assumption for the image, the edges of reconstructed images tend to be frequently over-smoothed and sometimes produce staircase effects, which leads to the loss of low-contrast information. In order to mitigate this drawback, some scholars have improved this algorithm and many related algorithms have been put forward. These efficient and accurate algorithms have been widely investigated: 1) considering the anisotropic edge property of image by introducing penalty weights [25 –27]; 2) increasing more constraint on TV regular terms [28 –30]; and 3) changing the order of TV, considering the first order difference operator and applying other order difference operators instead of using any fixed combination [31 –33]. These algorithms have achieved a good performance in suppressing artifacts and preserving the edge structure details information of the reconstructed image.

In this paper, considering both local gradient direction continuity information and the anisotropic edge property of an image, we propose a new adaptive-weighted TV reconstruction algorithm (NAWTV). The rest of this paper is organized as follows. In Section 2, the imaging model, TV model and the process of building NAWTV model for CT image reconstruction is presented. In Section 3, the proposed algorithm is applied to solve sparse-view problems with numerical simulation experiments. In Section 4, we conclude the paper and some related issues discussion is presented.

2 Method

2.1 CT imaging model

The CT imaging system is usually described by linear equations. The discrete form of the reconstruction model can be expressed as: $p = Wf$ (1)

Where f denotes the image to be reconstructed, p is the projection data, W is the system matrix which associates the image with the projection data. The CT reconstruction problem can be expressed as a solution to the unknown image f based on the system matrix W and projection matrix p. The traditional iterative reconstruction algorithms such as ART and MLEM can reconstruct an image with under-sampled or limited projection data, but fail to meet the clinical needs. However, the image quality can be further improved by introducing a priori knowledge to regularized constraints in the reconstruction model [35].

2.2 TV minimization reconstruction algorithm

For most natural images are non-sparse, but if we can find a suitable sparse transform domain and transform it, we can take the sparsely transformed signal as sparse. Rapid variations often occur at the edges of some structures in medical images. Therefore, the L1 norm of image gradient is often used in sparse-view image reconstruction. As mentioned in Section 1, the objective function of the TV reconstruction algorithm can be expressed as: $f^{*} = \underset{f}{arg min} {∥ f ∥}_{TV}, s . t . Wf = p, f ⩾ 0$ (2) and the TV of image is defined as ${∥ f ∥}_{TV} = \sum_{s, t} \sqrt{{(f_{s, t} - f_{s - 1, t})}^{2} + {(f_{s, t} - f_{s, t - 1})}^{2}}$ (3)

Where s and t represent the coordinates of the current pixel. The constraint term is used to constrain the consistency between the measured projection data and the reconstructed model data. The TV reconstruction algorithm can effectively use the information of image sparsity and edge direction, and can reconstruct the image with high quality under the sparse-view projection data. However, this algorithm can’t fully adjust the local information of the image, which lead to the edges of the reconstructed image tend to be over-smoothed.

2.3 NAWTV minimization reconstruction algorithm

To further improve the performance of the TV reconstruction algorithm, we propose a new adaptive-weighted TV reconstruction algorithm (NAWTV) for sparse-view CT image reconstruction. This algorithm enhances the local gradient continuity by increasing the gradient directionality information and consolidates the sparsity of the image by improved weight, while TV reconstruction algorithm only considers the sparsity of gradient in horizontal and vertical directions. Therefore, the sharp edges of CT images can be better preserved and the image reconstruction results are more accurate for NAWTV reconstruction algorithm.

2.3.1 NTV model

Inspired by the wok in [35], here a new TV (NTV) model is presented. The NTV model not only uses the gradient information in horizontal and vertical directions, but also enforces the gradient continuity. In Fig. 1 (a), it is seen that the pixels in the smooth region have similar values to the surrounding pixels, and any direction of the partial gradient can be seen as continuous. And when the partial gradient is continuous in addition to its own direction, the edge of the image can be obtained. Similarly, the continuity of partial gradient in the neighborhood (f_s,t, f_s-1,t-1, f_s-1,t, f_s,t-1) can be considered. In Fig. 1 (b), it is seen that the gradient continuity depends on the direction of the image edge, if it exists in the neighborhood. For different directions, the gradient continuity constraints are: ${\begin{matrix} {∥ {Gxx}_{s, t} ∥}_{1} = {∥ {Gx}_{s, t} - {Gx}_{s - 1, t} ∥}_{1} = 0, & horizontal direction \\ {∥ {Gyy}_{s, t} ∥}_{1} = {∥ {Gy}_{s, t} - {Gy}_{s, t - 1} ∥}_{1} = 0, & vertical direction \\ {∥ {Gxy}_{s, t} ∥}_{1} = {∥ {Gx}_{s, t} + {Gy}_{s - 1, t} ∥}_{1} = 0, & left - lower direction \\ {∥ {Gyx}_{s, t} ∥}_{1} = {∥ {Gx}_{s, t} - {Gy}_{s, t} ∥}_{1} = 0, & right - upper direction \end{matrix}$ (4)

Fig. 1

(a) pixel distribution; (b) partial gradient distribution.

Thus, we can enforce the directional continuity of Gx and Gy by ∥Gxx_s,t ∥ ₁,∥Gyy_s,t ∥ ₁, ∥Gxy_s,t ∥ ₁, ∥Gyx_s,t ∥ ₁, where Gxx_s,t is the derivative of Gx_s,t along the horizontal direction and Gyy_s,t is the derivative of Gy_s,t along the vertical direction. Actually, ∥Gxx_s,t ∥ ₁ = ∥ f_s,t + f_s-1,t-1 - f_s,t-1 - f_s-1,t ∥ ₁ = ∥ Gyy_s,t ∥ ₁. When we calculate the sparsest gradient, ∥Gxx_s,t ∥ ₁ = ∥ Gyy_s,t ∥ ₁ can be regarded as the redundancy of the entire equation. We define NTV model: $NTV (f) = \sum_{s, t} \sqrt{{Gx}^{2} + {Gy}^{2}} + γ \sum_{s, t} \sqrt{{Gxy}^{2} + {Gyx}^{2}}$ (5) where Gx = f_s,t - f_s,t-1, Gy = f_s,t - f_s-1,t, Gxy = f_s,t - f_s-1,t-1 and Gyx = f_s-1,t - f_s,t-1 are respectively horizontal, vertical and two diagonal partial gradient operators. γ plays a crucial role in the whole process since it measures the continuity of gradient direction. It can be understood that the NTV model is an extension of TV. When the weight γ → 0, the NTV model is equal with TV model.

2.3.2 NAWTV model

In 2012, Liu et al. [25] proposed an adaptive-weighted total variation (AwTV) reconstruction algorithm. This algorithm improved the preservation of edge details by an adaptive-weighted which brings the local information into the above conventional TV model. We apply the adaptive-weighted technique in AwTV algorithm into the above NTV algorithm, we called it NWTV: $\begin{matrix} NWTV (f) = \sum_{s, t} \sqrt{ω_{11} * {Gx}^{2} + ω_{12} * {Gy}^{2}} + γ \sum_{s, t} \sqrt{ω_{21} * {Gxy}^{2} + ω_{22} * {Gyx}^{2}} \\ ω_{11} = e^{- \frac{{Gx}^{2}}{h^{2}}}, ω_{12} = e^{- \frac{{Gy}^{2}}{h^{2}}}, ω_{21} = e^{- \frac{{Gxy}^{2}}{h^{2}}}, ω_{22} = e^{- \frac{{Gyx}^{2}}{h^{2}}} \end{matrix}$ (6) where h is a scale factor which controls the strength of the diffusion during each iteration, and the scale factor h can be set as some fixed value according to some experimental trials or histogram of the gradient magnitude throughout the image [25, 36]. The weight kernel function serves as an essential role in the smoothness of the edges and details of the reconstructed image. The main purpose of adding weight in image reconstruction is to make the reconstruction algorithm clearly distinguish the edges and smooth regions of images. This is to say the weight is closely related to the neighborhood similarity, which increases as the similarity increases.

In 2008, Tian et al. [37] compared the performance of several kernel functions in denoising. The idea of these kernel functions is in line with our requirements for weight functions. Several weight functions are shown in Fig. 2, the Euclidean distance is used to represent the neighborhood pixel similarity, and the neighborhood pixel similarity is inversely proportional to the Euclidean distance.

Fig. 2

Diagram of the weighted function curves.

From Fig. 2, it can be seen the exponential function (a) and cosine function (a) both get a larger weight when the similarity of neighborhood pixels is relatively high. However, with the decrease of similarity, the exponential function (a) decrease rapidly and cosine function (a) decreases slowly. When the relative parameters are adjusted to change the trend of exponential function, the range of weights will be smaller, such as the exponential function (b); the relative parameters are adjusted to change the trend of the cosine function, so that some small weights of the local image will be lost, such as the cosine function (b). Therefore, the exponential function has under-weighted shortcoming while the cosine function has over-weighted shortcoming.

To balance both under-weighted and over-weighted case, a new exponential cosine kernel function is proposed by combining exponential function and cosine function. From the exponential cosine function (a) and (b) in Fig. 2, it can be seen this new weight function can easily adjust the trend of weight change in the range of 0–1. This new weight is defined as follows: $ω (s, t) = {\begin{matrix} e^{- \frac{d {(s, t)}^{2}}{h_{1}^{2}}} cos (\frac{π d (s, t)}{2 h_{2}}), & d (s, t) ⩽ h_{2} \\ 0, & d (s, t) > h_{2} \end{matrix}$ (7)

where the Euclidean distance d (s, t) is related to the direction gradient, h₂ is a scale factor which controls the range of weight, the value of h₂ can be determined by maximum Euclidean distance to ensure the minimum value of the weight is 0. Replacing the weight function in formula (6) with the weight function defined by formula (7), a new adaptive-weighted TV (NAWTV) imaging model is proposed as follows: $\begin{matrix} f^{*} = \underset{f}{arg min (} {∥ f ∥}_{wtv} + γ {∥ f ∥}_{wtv 1}) s . t . Wf = p, f ⩾ 0 . \\ {∥ f ∥}_{wtv} = \sum_{s, t} \sqrt{ω_{11} * {Gx}^{2} + ω_{12} * {Gy}^{2}} \\ {∥ f ∥}_{wtv 1} = γ \sum_{s, t} \sqrt{ω_{21} * {Gxy}^{2} + ω_{22} * {Gyx}^{2}} \end{matrix}$ (8) where the weight function is expressed as: $\begin{matrix} ω_{11} = e^{- \frac{{Gx}^{2}}{h_{1}^{2}}} cos (\frac{π Gx}{2 h_{2}}), ω_{12} = e^{- \frac{{Gy}^{2}}{h_{1}^{2}}} cos (\frac{π Gy}{2 h_{2}}), \\ ω_{21} = e^{- \frac{{Gxy}^{2}}{h_{1}^{2}}} cos (\frac{π Gxy}{2 h_{2}}), ω_{22} = e^{- \frac{{Gyx}^{2}}{h_{1}^{2}}} cos (\frac{π Gyx}{2 h_{2}}) \end{matrix}$ (9)

2.3.3 Implementation of NAWTV algorithm

Inspired by the TV implementation strategy, the steepest descent method and FISTA algorithm [38] are employed in NAWTV implementation for solving the optimization problem. In summary, the implementation of the NAWTV method can be described as follows:

(A) Initializing the beginning reconstructed image: $n = 1 and f^{ART - DATA} (n, 0) = 0$ where n is the total iteration number and f^ART-DATA is the reconstructed image by ART algorithm.

(B) The ART algorithm is used to reconstruct the image and calculate the median image satisfying the projection data consistency, for m = 1, …, N_data: $f_{j}^{ART - DATA} (n, m) = f_{j}^{ART - DATA} (n, m - 1) + λ \frac{W_{ij}}{\sum_{j = 1}^{N} W_{ij}} (p_{i} - \sum_{j = 1}^{N} W_{ij} f_{j}^{ART - DATA} (n, m - 1))$ where N_data is the total number of rays.

(C) Positivity Constraint: $f^{pos} (n) = {\begin{matrix} f^{ART - DATA} (n) & , & f^{ART - DATA} (n) ⩾ 0 \\ 0 & , & f^{ART - DATA} (n) < 0 \end{matrix}$

(D) Calculating and updating weights: $\begin{matrix} ω_{11} = e^{- \frac{{Gx}^{2}}{h_{1}^{2}}} cos (\frac{π Gx}{2 h_{2}}), ω_{12} = e^{- \frac{{Gy}^{2}}{h_{1}^{2}}} cos (\frac{π Gy}{2 h_{2}}), \\ ω_{21} = e^{- \frac{{Gxy}^{2}}{h_{1}^{2}}} cos (\frac{π Gxy}{2 h_{2}}), ω_{22} = e^{- \frac{{Gyx}^{2}}{h_{1}^{2}}} cos (\frac{π Gyx}{2 h_{2}}) \end{matrix}$

(E1) TV gradient descent initialization: $f^{TV - GRAD} (n, 1) = f^{pos} (n)$

(E2) According equation (8) and (9), TV gradient descent method is applied, for m = 1, …, N_grad $\begin{matrix} v (n, m) = \frac{\partial f_{wtv}}{\partial f_{s, t}} + γ \frac{\partial f_{wtv 1}}{\partial f_{s, t}} \\ d = ∥ f^{TV - GRAD} (n, 1) - f^{ART - DATA} (n) ∥ \\ f^{ART - DATA} (n) = f^{TV - GRAD} (n, 1) \\ f^{TV - GRAD} (n, m + 1) = f^{TV - GRAD} (n, m) - bd \frac{v}{∥ v ∥} \end{matrix}$ where N_grad is the iteration number of TV term.

(F) Fasten the convergence with FISTA algorithm: $\begin{matrix} t_{n + 1} = \frac{1 + \sqrt{1 + 4 t_{n}^{2}}}{2} \\ f^{ART - DATA} (n + 1, 0) = f^{TV - GRAD} (n) + \frac{t_{n} - 1}{t_{n + 1}} (f^{TV - GRAD} (n + 1) - f^{TV - GRAD} (n)) \end{matrix}$

(G) Return to (B) until the stopping criterion is satisfied.

3 Numerical results

To investigate the feasibility of proposed NAWTV algorithm for sparse-view CT reconstruction, we implemented it in MATLAB2016b under a PC with 3.60 GHz Intel Core i7-7700 CPU, 8.0 G memory, Windows 8 operating system environment. In experiments investigation, the Shepp-Logan phantom [39] and the low-contrast FORBILD head phantom [40] were used to simulate the fan-beam projection data based on circular orbit, and the size of simulation data is 256*256. In the simulated environmental, we use an equi-view virtual detector, the edge ray just covers the target. And the distance between source and center of rotation is 655.36 mm, the number of detector channels was set to 256. In the following numerical experiments, we use TV, AwTV, NTV, NWTV and NAWTV algorithms to reconstruct the simulation data and compare their experimental results, the maximum iteration number is set 1000.

3.1 Evaluation measures

To verify the effectiveness of the proposed method, the performance of the reconstructed image is evaluated by the root mean square error (RMSE) and the universal quality index (UQI). The definitions of RMSE and UQI are as follows: $RMSE = \sqrt{\frac{1}{m \times n} \sum_{s = 1, t = 1}^{m, n} {(f_{s, t} - f_{s, t}^{*})}^{2}}$ (10) $UQI (f, f^{*}) = \frac{4 σ_{{ff}^{*}} \bar{f} \times {\bar{f}}^{*}}{({\bar{f}}^{2} + {\bar{f}}^{* 2}) (σ_{f}^{2} + σ_{f *}^{2})}$ (11) where f_s,t and $f_{s, t}^{*}$ are the original image and the reconstructed image respectively, $\bar{f}$ and ${\bar{f}}^{*}$ are respectively the mean of f and f^*, $σ_{f}^{2}$ and $σ_{f^{*}}^{2}$ are respectively the variance of f and f^*, σ_{ff
^*} are the covariance of f and f^*. RMSE is used to evaluate the reconstruction accuracy, which ranges from 0 to 1, and the value closer to 0 suggests the reconstruction results better. UQI ranges from –1 to 1, and its value closer to 1 suggests a better similarity between the reconstructed and original images.

3.2 Shepp-Logan phantom based experiment

For Shepp-Logan phantom based experiment, the phantom is uniformly sampled from 0 to 2π, and 24 projections are obtained. The two weight parameters in the implementation of the NAWTV algorithm are set to h₁ = 0.21 and h₂ = 1, and the gradient’s continuity weight is set to γ = 0.2. The weight parameter in the implementation of AwTV algorithm and NWTV algorithm is set to h = 0.21. The relaxation factor of the ART algorithm is set to λ = 1 and it is repeated once per total cycle. The outer iteration number for all methods are 1000 and the inner iteration number of TV terms are all N_data = 20. The reconstructed images with TV, AwTV, NTV, NWTV, NAWTV algorithms and performance evaluation measures are shown in Fig. 3 and Table 1 respectively.

Fig. 3

Reconstructed results of the Shepp-Logan phantom from 24 noise-free projection data. (a) Original image. (b-f) The reconstructed results of NAWTV, NWTV, NTV, AwTV and TV.

Table 1

Performance evaluation of the reconstructed results by different algorithms from noise-free projection data

	TV	AwTV	NTV	NWTV	NAWTV
UQI	0.994147	0.997835	0.997357	0.999677	0.999685
RMSE	0.015276	0.006901	0.009927	0.002164	0.002139

From Fig. 3 we can see that using 24 projection data with 1000 times iterations, the reconstructed image using NAWTV and NWTV method almost have no difference with original image in visual effect. While the reconstructed image using NTV and AwTV method contains less artifacts, and the inner distribution near edge is more uniform than the image reconstructed with TV. Combined with the quantitative evaluation in Table 1, we can conclude that the NAWTV algorithm performs best among five of them.

To further investigate whether the performance of proposed algorithm is better, the comparison of the midline profiles of horizontal and vertical between original Shepp-Logan phantom and reconstructed images using five different reconstruction algorithms shown in Fig. 4. It can be seen from the local zoom in figure (a’) and (b’) that the horizontal and vertical profile of the image reconstructed by the NAWTV and NWTV algorithm are the closest to the original profile. In addition, the algorithm is well in protecting the edge property of the original image. At last, the change curves of the UQI and RMSE of this experiment are presented in Fig. 5. From the figure, we can see that the NTV algorithm converges faster than the TV algorithm, while NWTV and NAWTV converge faster than AwTV. The convergence speed difference is not very big between NWTV and NAWTV in Fig. 5. To facilitate observation, a local zoom of this two change curves is also shown in Fig. 5. In conclusion, we can get that the convergence speed of NAWTV algorithm is the fastest among five of them.

Fig. 4

The comparison of the midline profiles between reconstructed images using five different reconstruction algorithms and Shepp-Logan phantom. (a) is the comparison of horizontal midline, (b) is the comparison of vertical midline. (a’) and (b’) is the corresponding zoomed in part of (a) and (b).

Fig. 5

The RMSE and UQI change curves of reconstructed images with different reconstruction algorithms at 24 projection data for Shepp-Logan phantom. The center rectangle region is a local zoomed in results of the change curves of NAWTV and NWTV algorithms in the dashed frame.

3.3 FORBILD head phantom based simulation data reconstruction

In order to further validate the robustness of the proposed algorithm in sparse-view reconstruction, simulation experiments were performed on FORBILD phantoms with more detail. The simulation model is uniformly sampled from 0 to 2π, and 36 projection data are obtained. The two weight parameters in the implementation of the NAWTV algorithm are set to h₁ = 0.21 and h₂ = 1.8, and the other reconstruction parameters are consistent with Section 3.2. The reconstructed images and performance evaluations are respectively shown in Fig. 6 and Table 2. The comparison of the midline profiles of horizontal and vertical between original FORBILD head phantom and reconstructed images using five different reconstruction algorithms shown in Fig. 7.

Fig. 6

Reconstructed results of the FORBILD head phantom from 36 noise-free projection data. (a) Original image. (b-f) The Reconstructed results of NAWTV, NWTV, NTV, AwTV and TV.

Fig. 7

The comparison of the midline profiles between reconstructed images using five different reconstruction algorithms and FORBILD head phantom. (a) is the comparison of horizontal midline, (b) is the comparison of vertical midline. (a’) and (b’) is the corresponding zoomed in part of (a) and (b).

Table 2

Performance evaluation of the reconstructed results by different algorithms from noise-free projection data

	TV	AwTV	NTV	NWTV	NAWTV
UQI	0.983815	0.993146	0.991191	0.999054	0.999070
RMSE	0.055058	0.032853	0.042934	0.011592	0.011481

From Fig. 6 we can see that using 36 projection data with 1000 times iterations, the reconstructed image using NAWTV and NWTV method have little difference in visual effect. Although there are some pleated artifacts, the reconstructed image of NAWTV and NWTV are the most similar to the original image in several reconstruction algorithms. In addition, the reconstructed image using NTV and AwTV method contains slightly more artifacts, but the visual effect of the reconstructed image is still better than the image reconstructed with TV. Combined with the quantitative evaluation in Table 2, we can conclude that the NAWTV algorithm is the best method for FORBILD phantom reconstruction.

From the local zoomed in results of the horizontal and vertical profile in Fig. 7 (a’) and (b’), we can see that the profile of the reconstructed image fluctuates more volatile, but the profile of NAWTV and NWTV algorithm are still the closest to the original profile in several reconstructed results, especially in the corner region. The change curves of the UQI and RMSE of the reconstructed image are shown in Fig. 8. Among them, the local zoomed in results of the change curves of NAWTV and NWTV algorithm are also shown in Fig. 8. In the reconstruction of phantom data with more details, the convergence speed of NAWTV algorithm is still the best, and the result is consistent with the simulation result of Shepp-Logan phantom.

Fig. 8

The RMSE and UQI change curves of reconstructed images with different reconstruction algorithms at 36 projection data for FORBILD head phantom. The center rectangle region is a local zoomed in results of the change curves of NAWTV and NWTV algorithms in the dashed frame.

3.4 Validation of the algorithm’s anti-noise performance

To validate the anti-noise performance of proposed algorithm, 24 projection data are uniformly sampled from the Shepp-Logan model, then random noise is added in projection data. Five reconstruction algorithms’ performance on noise projection data are compared in detail. The reconstruction parameters in the implementation of each algorithm are consistent with Section 3.2. The reconstructed images and performance evaluations are respectively shown in Fig. 9 and Table 3. From the reconstructed results and the evaluation measures, the edge of the reconstructed image on noise data is a little blurred comparing with the noise-free case, and the evaluation index is slightly reduced. For this noise projection data reconstruction experiment, the proposed NAWTV algorithm performs the best among five of them.

Fig. 9

Reconstructed results of the Shepp-Logan phantom from 24 noise projection data. (a) Original image. (b-f) The reconstructed results of NAWTV, NWTV, NTV, AwTV and TV.

Table 3

Performance evaluation of the reconstructed results by different algorithms from noise projection data

	TV	AwTV	NTV	NWTV	NAWTV
UQI	0.988243	0.995928	0.989167	0.996598	0.996634
RMSE	0.023691	0.009069	0.022907	0.007445	0.007347

The comparison of the midline profiles of horizontal and vertical between Shepp-Logan phantom and reconstructed images on noise projection data are shown in Fig. 10. It can be seen from the local zoomed in views in Fig. 10 (a’) and (b’), the reconstructed image profile of NTV is closer to the original image than the reconstructed image profile of TV, while the reconstructed image profile of NAWTV and NWTV are closer to the original image than the reconstructed image profile of AwTV. The change curves of the UQI and RMSE of the reconstructed images from noise projection data and local zoom in results of the change curves are shown in Fig. 11. It can be clearly seen that the NTV algorithm converges faster than the TV algorithm, the NWTV algorithm converges faster than the AwTV algorithm, while the NAWTV algorithm converges the fastest among five of them. This experimental result is consistent with the simulation result of the noise-free projection data reconstruction. In all, whether for noise-free or noise projection data reconstruction, the proposed algorithm NAWTV gives a better reconstruction performance and faster convergence than other algorithms.

Fig. 10

The comparison of the midline profiles between reconstructed images from projection data with noise on Shepp-Logan phantom. (a) is the comparison of horizontal midline, (b) is the comparison of vertical midline. (a’) and (b’) is the corresponding zoomed in part of (a) and (b).

Fig. 11

The RMSE and UQI change curves of reconstructed images with different reconstruction algorithms at 24 noise projection data for Shepp-Logan phantom. The center rectangle region is a local zoomed in results of the change curves of NAWTV and NWTV algorithms in the dashed frame.

4 Discussion and conclusion

In this paper, the simulation experiments are implemented with Shepp-Logan phantom and FORBILD head phantom. The reconstructed image and numerical results of TV, AwTV, NTV, NWTV and NAWTV were compared qualitatively and quantitatively. It can be seen in terms of visual effect or performance evaluations, the algorithm of NAWTV has the biggest advantages in reconstructing image quality. The NAWTV reconstruction algorithm has a higher reconstruction accuracy, which can suppress artifacts and noise while preserving more edge structure information.

Since the lack of clinical CT projection data, this paper only carried out on simulation experiments, the experiment on real clinical projection data will be studied in the near future. In addition to the reconstruction algorithm is solved by gradient descent method, which takes a long time. The proposed algorithm can be further improved with state of the art technique, such as GPU based acceleration technique and new solution algorithms. In the further, we will try to explore the limited-view projection reconstruction problem to reduce radiation dose.

In summary, we proposed a new adaptive-weighted TV reconstruction algorithm for sparse-view projection data and demonstrated its promising performance. First, the proposed method introduces local gradient direction continuity to the TV reconstruction algorithm and the NTV model is proposed to improve the local gradient continuity. Second, the adaptive-weighted technique of AwTV algorithm is introduced into the NTV model to enhance the gradient sparsity. Last, the NAWTV model is proposed by improving weight function makes the weight function can balance the under-weighted and over-weighted cases. The NAWTV is used as a regularization term to constrain gradient sparsity and gradient continuity for sparse-view CT image reconstruction. In addition, it also fastens the convergence by incorporating the FISTA step into ours methods.

Footnotes

Acknowledgments

This work is supported by the Natural Science Foundation of Shanghai, China (Grant No. 18ZR1426900) and the National Natural Science Foundation of China (Grant No. 61201067).

References

Brenner

D.J.

and Hall

E.J.

, Computed tomography–An increasing source of radiation exposure, New England Journal of Medicine 357 (2007), 2277–2284.

Gao

, et al., Low-dose X-ray computed tomography image reconstruction with a combined low-mAs and sparse-view protocol, Optics Express 22 (2014), 15190–15210.

Dong

, et al., A Simple Low-dose X-ray CT Simulation from High-dose Scan, IEEE Transactions on Nuclear Science 62 (2015), 2226–2233.

Zhang

, et al., Noise correlation in CBCT projection data and its application for noise reduction in low-dose CBCT, Medical Physics 41 (2014), 031906.

Huang

, et al., Projection data restoration guided non-local means for low-dose computed tomography reconstruction, IEEE International Symposium on Biomedical Imaging 48 (2011), 1167–1170.

Niu

, et al., Iterative reconstruction for sparse-view x-ray CT using alpha-divergence constrained total generalized variation minimization, Journal of X-ray Science and Technology 25 (2017), 673–688.

Feldkamp

L.A.

, Davis

L.C.

and Kress

J.W.

, Practical cone-beam algorithm, Josa A 1 (1984), 612–619.

Zeng

, et al., Spectral CT Image Restoration via an Average Image-Induced Nonlocal Means Filter, IEEE Transactions on Biomedical Engineering 63 (2016), 1044–1057.

Gordon

, Bender

and Herman

G.T.

, Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and x-ray photography, Journal of Theoretical Biology 29 (1970), 471–476.

10.

Dempster

A.P.

, Laird

N.M.

and Rubin

D.B.

, Maximum Likelihood from Incomplete Data via the EM Algorithm, Journal of the Royal Statistical Society 39 (1977), 1–38.

11.

Yan

and Li

, “CT image reconstruction algorithm,” Science Press, (2014), pp. 97–99.

12.

Candes

E.J.

, Romberg

and Tao

, Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information, IEEE Trans Inf Theory 52 (2006), 489–509.

13.

Donoho

D.L.

, Compressed sensing, IEEE Trans Inf Theory 52 (2006), 1289–1306.

14.

Needell

and Vershynin

, Signal Recovery From Incomplete and Inaccurate Measurements Via Regularized Orthogonal Matching Pursuit, IEEE Journal of Selected Topics in Signal Processing 4 (2010), 310–316.

15.

Niu

and Zhu

, Accelerated barrier optimization compressed sensing (ABOCS) reconstruction for cone-beam CT: Phantom studies, Medical Physics 39 (2012), 3892–3893.

16.

Rantala

, et al., Wavelet-based reconstruction for limited-angle X-ray tomography, IEEE Transactions on Medical Imaging 25 (2006), 210–217.

17.

Starck

J.L.

, Candès

E.J.

and Donoho

D.L.

, The curvelet transform for image denoising, IEEE Transactions on Image Processing 11 (2002), 670–684.

18.

Vandeghinste

, et al., Iterative CT Reconstruction Using Shearlet-Based Regularization, IEEE Transactions on Nuclear Science 60 (2013), 3305–3317.

19.

Rudin

L.I.

, Osher

and Fatemi

, Nonlinear total variation based noise removal algorithms, Eleventh International Conference of the Center for Nonlinear Studies on Experimental Mathematics: Computational Issues in Nonlinear Science: Elsevier North-Holland 60 (1992), 259–268.

20.

, Chen

, Wang

, et al., Sparse CT reconstruction based on multi-direction anisotropic total variation (MDATV), Biomedical Engineering Online 13 (2014), 92–119.

21.

Kim

, et al., Non-local total-variation (NLTV) minimization combined with reweighted L1-norm for compressed sensing CT reconstruction, Physics in Medicine & Biology 61 (2016), 6878–6891.

22.

Chen

, et al., A new Mumford—Shah total variation minimization based model for sparse-view x-ray computed tomography image reconstruction, Neurocomputing 285 (2018), 74–81.

23.

Sidky

E.Y.

and Pan

, Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization, Physics in Medicine & Biology 53 (2008), 4777–4807.

24.

Tian

, Jia

, Yuan

, Pan

and Jiang

S.B.

, Low-dose CT reconstruction via edge-preserving total variation regularization, Physics in Medicine & Biology 56 (2011), 5949–5967.

25.

Liu

, Ma

, Fan

, Liang, Z, Adaptive-weighted Total Variation Minimization for Sparse Data toward Low-dose X-ray Computed Tomography Image Reconstruction, Physics in Medicine & Biology 57 (2012), 7923–7956.

26.

Zhang

, Zhang

and Zhou

, Accurate sparse-projection image reconstruction via nonlocal TV regularization, The Scientific World Journal 22 (2014), 458496.

27.

Chen

M.Y.

, et al., Computed Tomography Image Reconstruction from Few-Views Data by Multi-Directional Total Variation, Journal of Medical Imaging & Health Informatics 5 (2015), 309–316.

28.

Yazdanpanah

A.P.

and Regentova

E.E.

, Sparse-View CT Reconstruction Using Curvelet and TV-Based Regularization, International Conference Image Analysis and Recognition 2016, pp. 672–677.

29.

Ertas

, Yildirim

, Kamasak

and Akan

, Iterative image reconstruction using non-local means with total variation from insufficient projection data, Journal of X-ray Science and Technology 24 (2016), 1–8.

30.

Liu

, et al., Median prior constrained TV algorithm for sparse view low-dose CT reconstruction, Computers in Biology & Medicine 60 (2015), 117–131.

31.

Debatin

and Hesser

, Accurate low-dose iterative CT reconstruction from few projections by Generalized Anisotropic Total Variation minimization for industrial CT, Journal of X-ray Science and Technology 23 (2015), 701–726.

32.

Zhang

H.M.

, et al., Constrained Total Generalized p-Variation Minimization for Few-View X-Ray Computed Tomography Image Reconstruction, Plos One 11 (2016), e0149899.

33.

, et al., A weighted difference of L1 and L2 on the gradient minimization based on alternating direction method for circular computed tomography, Journal of X-ray Science and Technology 25 (2017), 813–829.

34.

Zhang , et al., Few-view cone-beam CT reconstruction with deformed prior image, Medical Physics 41 (2014), 527–527.

35.

Shu

and Ahuja

, Hybrid compressive sampling via a new total variation TVL1, European Conference on Computer Vision 6316 (2010), 393–404.

36.

Perona

, Scale-space and edge detection using anisotropic diffusion, IEEE Transactions on Pattern Analysis & Machine Intelligence 12 (2002), 629–639.

37.

Tian

, Yu

W.Y.

and

S.L.

Xie , On the kernel function selection of nonlocal filtering for image denoising, International Conference on Machine Learning and Cybernetics 5 (2008), 2964–2969.

38.

Beck

and Teboulle

, A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems, Siam Journal on Imaging Sciences 2 (2009), 183–202.

39.

Shepp

L.A.

and Logan

B.F.

, The Fourier reconstruction of a head section, IEEE Transactions on Nuclear Science 21 (1974), 21–43.

40.

, et al., Simulation tools for two-dimensional experiments in x-ray computed tomography using the FORBILD head phantom, Physics in Medicine & Biology 57 (2012), 237–252.