Computational and Mathematical Methods in Medicine
Volume 2015 (2015), Article ID 691021, 9 pages
http://dx.doi.org/10.1155/2015/691021

Research Article

NUFFT-Based Iterative Image Reconstruction via Alternating Direction Total Variation Minimization for Sparse-View CT

Bin Yan, Zhao Jin, Hanming Zhang, Lei Li, and Ailong Cai

China National Digital Switching System Engineering and Technological Research Center, Zhengzhou 450002, China

Received 18 September 2014; Accepted 11 January 2015

Academic Editor: Yi Gao

Copyright © 2015 Bin Yan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Sparse-view imaging is a promising scanning method which can reduce the radiation dose in X-ray computed tomography (CT). Reconstruction algorithm for sparse-view imaging system is of significant importance. The adoption of the spatial iterative algorithm for CT image reconstruction has a low operation efficiency and high computation requirement. A novel Fourier-based iterative reconstruction technique that utilizes nonuniform fast Fourier transform is presented in this study along with the advanced total variation (TV) regularization for sparse-view CT. Combined with the alternating direction method, the proposed approach shows excellent efficiency and rapid convergence property. Numerical simulations and real data experiments are performed on a parallel beam CT. Experimental results validate that the proposed method has higher computational efficiency and better reconstruction quality than the conventional algorithms, such as simultaneous algebraic reconstruction technique using TV method and the alternating direction total variation minimization approach, with the same time duration. The proposed method appears to have extensive applications in X-ray CT imaging.

1. Introduction

X-ray computed tomography (CT) has been widely used for imaging applications in various fields, such as industrial nondestructive testing [1] and medical diagnosis [2], for its advantages of noninvasive and high spatial resolution. However, in many practical applications of X-ray computed tomography, complete projection data set cannot be obtained because of the limitation of scanning time, space, dose, and so on. Therefore, sparse angle scanning scheme is adopted to tackle these problems. On one hand, this scheme can speed up the scanning rate and decrease the X-ray radiation dose, such as breast and vascular imaging [3–6]. On the other hand, sparse data sampling can save much scanning time and it is of practice value when high reconstruction precision is not that urgent. To solve the sparse-view reconstruction problem, the classical methods should be upgraded, and a new algorithm needs to be developed.

Given the unsatisfactory Tuy-Smith condition [7, 8] in sparse-view, a CT image cannot be accurately reconstructed via analytic method. To solve the ill-posed problem [9, 10], numerous iterative algorithms [11–13] have been proposed based on spatial domain. However, these iterative algorithms are time-consuming and have a great demand for hardware resources. Despite applying hardware speedup technology, such as an ordinary graphics processing unit [14], these algorithms still consume a considerable amount of time.

Compressive sensing theory by Candés et al. [15–17] provided a new idea for the exact recovery of an image from the sparse samples of its discrete Fourier transform. The exact reconstruction relies on the assumption that there exists sparse representation for an image. A number of cases are known to have sparse gradient-magnitude images. In some cases, minimizing the total variation (TV) can generate accurate images from sparse samples [18–20]. Therefore, combining TV regularization with the iterations of the simultaneous algebraic reconstruction technique (SART), hereinafter called SART–TV [21], can improve reconstruction image quality while decreasing mean-squared error. Based on the projection onto the convex sets (POCS) algorithm, the adaptive steepest descent-POCS (ASD-POCS) algorithm [20] can effectively handle incomplete datasets and demonstrates excellent performance in sparse-view CT applications.

The rapid increase in the size of scanning data has highlighted the importance of reducing reconstruction time and improving reconstruction quality. It is known that processing the same signal in the frequency space is faster than that in the spatial domain by fast Fourier transform (FFT). Several algorithms for image reconstruction in the frequency space can also be developed on the basis of fast Fourier transform. Several studies have been conducted to achieve this goal. In 1981, Stark et al. [22] developed direct Fourier methods (DFM) using central slice theorem and obtained favorable results. In 2003, Seger and Danielsson [23] analyzed the missing projection data in the frequency domain and proposed a reconstruction method for the scanned timber data according to Fourier transform. In 2013, Fahimian et al. [24] presented a Fourier-based iterative reconstruction in medical X-ray CT, and numerical experiment results showed that this method required less computation time than other iterative algorithms. These achievements facilitated the development of an improved algorithm for solving the sparse-view reconstruction problem in the frequency domain.

Because of the limitation of FFT, that is, its unsuitability for application to nonuniform samples, this technique requires further enhancement to improve its universality. To this end, nonuniform FFT (NUFFT) [25] has been recently developed to overcome this limitation without increasing the computation complexity of FFT. NUFFT is also basis of the proposed reconstruction algorithm in this study. Motivated by the feature of NUFFT for data distribution, some approaches have been proposed to reconstruct a CT image to deal with frequency data. In 2004, Matej et al. [26] proposed an iterative tomographic image reconstruction method using NUFFT and obtained better results with this technique than with the filtered back-projection (FBP) algorithm. In 2006, Zhang-O’Connor and Fessler [27] proposed Fourier-based forward- and back-projectors for fan-beam tomographic image reconstruction. However, these proposed NUFFT cannot effectively solve the problems in sparse-view image reconstruction. The NUFFT just was especially applied as a transition during iteration in spatial domain, which in turn burdened computation consumption.

In this paper, our study aims to present a promising contribution to the task of image reconstruction from sparse-view by combining the alternating direction total variation minimization (ADTVM) technique with NUFFT to establish a new method which is suitable for large-scale reconstruction because of its low computational requirement. The algorithm is developed under the framework of alternating direction method (ADM) which shows high efficiency and stability. The advantages of the proposed algorithm are verified by the results of several groups of experiments.

The organization of this paper is organized as follows. Section 1 concisely reviews the basic CT reconstruction and the state of the art of sparse-view image reconstruction. Section 2 describes the basic principles of the proposed method, including the reconstruction model and the corresponding algorithm based on NUFFT and ADTVM. Section 3 demonstrates groups of typical experiments results, including numerical simulation and real data ones. Section 4 discusses the findings of the experiments and concludes the paper.

2. Method and Material

2.1. Image Reconstruction Model

In this work, we consider temporarily parallel geometry. In parallel geometry, 2D function (objection function) is defined in a compact support of spatial domain, which means that it vanishes outside a finite region of the plane. In the plane, the general formation for the line integral, known as the Radon transform of , is

In Figure 1, projection consists of a collection of line integrals (1) taken along straight parallel lines in the plane that means a collection of with constant and .

Figure 1: An object and its parallel projection .

CT image reconstruction is an inverse problem, and the observed projections should be converted into images which reflect the distribution of the attenuation coefficient of the interested physical object. A conventional reconstruction method, that is, the direct Fourier methods (DFM), is established based on the Fourier slice theorem. Basically, the steps of DFM can be summarized as follows:

(a)1D discrete Fourier transform of the parallel projections taken at different angles;(b)polar to Cartesian grid interpolations;(c)2D inverse Fourier transform.

From the perspective of DFM, the observation equation in the Fourier domain can be expressed as follows:

where is observed by the Fourier transform of the measured projection data and is the frequency variable of Fourier transform. The process proceeds with the use of FFT and is characterized by high accuracy. According to Fourier central slice theory, the 1D FFT of the projection is equal to , which is derived from the 2D FFT of the reconstructed image in a certain angle. Therefore, the relationship is described as follows:

where . The equation shows an obvious equivalence corresponding to frequency projection with in polar coordinates.

To avoid interpolation errors, such as DFM in the image and frequency domains, this study introduces NUFFT, which can translate polar coordinates into the image space without interpolation. This technique can significantly improve the accuracy of reconstruction. Let represent the NUFFT operator, such that the following equation can be derived:

The reconstruction module can be discretely shown as follows:

where the (observed) constant and the variant are the vector forms of Fourier sampling and objection function, respectively. Matrix stands for the NUFFT of . The Fourier transform in (5) and its adjoint can be implemented by using FFT to generate a fast and accurate evaluation.

In sparse-view reconstruction, (5) is ill-posed, and the projection data are insufficient for exact reconstruction. Mathematically, the problem that we consider here involves insufficient data, such that (5) is underdetermined. To solve this linear and underdetermined equation, we specify a TV minimization algorithm that considers the reconstruction to be a task of finding the best solution to the following optimization problem:

where denotes the discretization of the TV term and . By applying the directional gradients operators [20, 28], model (6) can also be written as follows:

where is the fidelity parameter to control the data consistency in the object function.

Therefore, the overall reconstruction flowchart can be summarized as Figure 2.

Figure 2: The flowchart of sparse-view image reconstruction for the model.

2.2. NUFFT with ADM for the Model

The above constrained optimization is addressed by converting the equation into its unconstrained form by applying the augmented Lagrange function:

where is a scalar that denotes the penalty coefficient and denotes the multipliers. The minimization processes with respect to variables and cannot be easily realized simultaneously by directly performing the optimization. Moreover, decomposing the variables by using ADM has a low computation cost. The ADM approach decouples the augmented Lagrange function into two subproblems, namely, the -subproblem and the -subproblem [29].

The -subproblem can be written as follows:

The -subproblem is separable with respect to , and problem (9) can be efficiently solved by using the shrinkage operator [30], which is expressed as follows:

In addition, with the aid of , the optimization of -subproblem can be achieved by solving the following:

is clearly a quadratic function, the gradient of which is expressed as follows:

Force and the exact solution for is presented as follows:

where denotes the Moore-Penrose inverse of matrix . Theoretically, the exact minimizer can be used to solve the -subproblem. However, the inverse or pseudo-inverse is too costly to compute numerically at each iteration. The augmented Lagrangian function (8) is expected to be minimized by solving the -subproblem and the -subproblem alternately. Therefore, solving the -subproblem accurately at each sweep may be unnecessary. A robust and efficient nonmonotone alternating direction algorithm [31] is used to solve problem (13).

By using the solutions of and , the multipliers are updated as follows:

The optimized solution for (8) is attained by circularly applying (10) and (13) until is minimized jointly with respect to .

2.3. Algorithm of the Overall Framework

All issues in handling the subproblems have been addressed in Section 2.2. In light of all derivations presented above, the new algorithm for solving the reconstruction problem can be stated as follows.

Algorithm 1. Input projection data . Initialize and starting points , for all . Set .

(1)make 1D FFT of with respect to

 while “not achieved maximum iteration loops,” Do(2)compute frequency domain via NUFFT;(3)compute by

(4)compute by

(5)update by

(6) End Do

In this study, the proposed NUFFT reconstruction technique is developed on the basis of ADTVM. This technique is called NUFFT-ADM. According to the above algorithm, the proposed method demonstrates fast convergence and effective iteration through ADM. This method can be effectively implemented in large-scale reconstruction in sparse-view because of its low computational cost, thus making this technique promising in practical applications.

3. Experiments

To verify the performance of the proposed algorithm, both numerical simulations and real CT scan data experiments are conducted. All experiments are performed on an AMAX Tesla workstation with Intel Xeon E5520 dual-core CPU 2.27 GHz and 24 GB memory. All programs are performed using MATLAB 2011a. In all experiments, the parameter of TV is that primary penalty parameter and secondary penalty parameter are 1024 and 32, respectively.

3.1. Numerical Phantom Simulation

The above algorithm is applied to validate its high efficiency. A group of 2D image reconstruction experiments are performed using a 2D Shepp-Logan phantom with a size of . This phantom is generated according to the definition of the ellipse phantom image. The scanning and reconstruction parameters in the experiment are listed in Table 1. The detector elements are equidistantly spaced 0.127 mm from one another.

Table 1: Parameters in the simulation of a sparse-view scan.

Parameters
Configuration

Detection elements
512

Source to axis distance
300 mm

Source to detection distance
600 mm

Views of projection data
18

Projection data
512 18

Reconstruction size
256 256 pixels

Pixel size
0.127 mm 0.127 mm

To demonstrate the reconstruction accuracy quantitatively, the root-mean-squared error (RMSE) is used as a measurement of the reconstruction error. RMSE is defined as follows:

where and denote the ideal phantom and the reconstruction image, respectively; and denotes the total number of pixels of the image. The image was reconstructed using SART-TV, ADTVM, the proposed method, respectively, and their results are presented in Figure 3. Two hundred iterations are performed for each algorithm. The profiles of these images along the central horizontal and vertical rows are shown in Figure 4 for the different methods.

(a) Phantom

(b) SART-TV

(c) ADTVM

(d) New method

Figure 3: Image reconstructions of the Shepp-Logan phantom in 60-view scan. Display Window (0.1 0.5). (a) shows the original image. (b) shows the image after applying the SART-TV algorithm at 200 iterations. (c) gives the image after vua the ADTVM algorithm at 200 iterations. (d) presents the reconstruction image after using the NUFFT-ADM algorithm at 200 iterations.

(a) SART-TV of the vertical rows

(b) SART-TV of the horizontal rows

(c) ADTVM of the vertical rows

(d) ADTVM of the horizontal rows

(e) New method of the vertical rows

(f) New method of the horizontal rows

Figure 4: The image profile in Figure 1 shows (a) vertical profiles along the center of the SART-TV result, (b) horizontal profiles along the center of the SART-TV result, (c) vertical profiles along the center of the ADTVM result, and (d) horizontal profiles along the center of the ADTVM result. (e) Vertical profiles along the center of the NUFFT-ADM result and (f) horizontal profiles along the center of the NUFFT-ADM result.

RMSE is used as an evaluation criterion for different iteration times. The result is shown in Figure 5. The accuracy and running time of each reconstruction method at different iterations are presented in Table 2 for comparison.

Table 2: Accuracy and the running time of different method.

Iteration number
SART-TV
ADTVM
NUFFT-ADM

RMSE
Running time
RMSE
Running time
RMSE
Running time

100
0.0377
177.786 s
0.0165
107.439 s
0.0079
8.986 s

200
0.0302
349.544 s
0.0015
199.133 s
0.0012
16.752 s

500
0.0174
856.158 s

508.376 s

51.863 s

Figure 5: RMSEs as functions of iterations for three different algorithms.

The RMSEs, as well as the accuracy and running time of different methods, show that NUFFT-ADM significantly outperforms SART-TV and ADTVM. On one hand, the convergence of the new method is faster than that of SART-TV because of the use of the optimal solution with ADM. On the other hand, by taking advantage of the frequency NUFFT operator instead of the projection and back-projection in the spatial domain, which consumes the greatest amount of time among all components, NUFFT-ADM has a higher computation capability than the general algorithm in the spatial domain.

3.2. Reconstruction Using Real Data

To further verify the performance of the proposed algorithm, several experiments are performed to reconstruct a head model from real data using the new method. We compare the proposed algorithm with SART-TV and ADTVM. Table 3 lists the scanning and reconstruction parameters. The detector elements are equidistantly spaced 0.635 mm from one another. The number of iterations for NUFFT-ADM, ADTVM, and SART-TV is 200.

Table 3: Parameters in the simulation of a sparse-view scan.

Parameters
Configuration

Detection elements
640

Source to axis distance
678 mm

Source to detection distance
1610 mm

Views of projection data
60

Projection data
60 72

Reconstruction size
256 256 pixels

Pixel size
0.582 0.582 mm²

The reconstruction results are presented in Figure 6.

(a) SART-TV (at 100 iterations)

(b) ADTVM (at 100 iterations)

(c) New method (at 100 iterations)

(d) SART-TV (at 300 iterations)

(e) ADTVM (at 300 iterations)

(f) New method (at 300 iterations)

Figure 6: Reconstruction results of SART-TV, ADTVM, and NUFFT-ADM.

The reconstruction acquired results using real data clearly show that the quality of the reconstructed image is improved as the number of iterations is increased. Under the same number of iterations, the reconstruction results of NUFFT-ADM are superior to those of SART-TV, especially in terms of the high-frequency information that shows the image detail or volatile part. Compared to ADTVM, the detail of the image by the new method is nearly the same.

4. Conclusions

An optimal algorithm based on NUFFT for CT image reconstruction is presented in this work. The validity of the NUFFT-ADM algorithm is verified by conducting numerical simulations and real data experiments. The reconstruction results show that the proposed reconstruction algorithm improves reconstruction quality, accelerates convergence speed, and demonstrates lower computation complexity than other iterative algorithms. That is, the NUFFT-ADM algorithm can practically deal with fast image reconstruction from sparse projection measurements to reduce the radiation dose in X-ray CT. In principle, the proposed method can be extended to fan-beam geometry via the rebinning method to broad its application.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was supported in part by the National High Technology Research and Development Program of China under the Grant no. 2012AA011603 and the National Natural Science Foundation of China under Grant no. 61372172.

References

H. E. Martz, S. G. Azevedo, J. M. Brase, K. E. Waltjen, and D. J. Schneberk, “Computed tomography systems and their industrial applications,” International Journal of Radiation Applications and Instrumentation, Part A: Applied Radiation and Isotopes, vol. 41, no. 10-11, pp. 943–961, 1990.

S. Hu, E. A. Hoffman, and J. M. Reinhardt, “Automatic lung segmentation for accurate quantitation of volumetric X-ray CT images,” IEEE Transactions on Medical Imaging, vol. 20, no. 6, pp. 490–498, 2001.

X. Gong, S. J. Glick, B. Liu, A. A. Vedula, and S. Thacker, “A computer simulation study comparing lesion detection accuracy with digital mammography, breast tomosynthesis, and cone-beam CT breast imaging,” Medical Physics, vol. 33, no. 4, pp. 1041–1052, 2006.

S. J. Glick, S. Thacker, X. Gong, and B. Liu, “Evaluating the impact of X-ray spectral shape on image quality in flat-panel CT breast imaging,” Medical Physics, vol. 34, no. 1, pp. 5–24, 2007.

C. J. Henri and T. M. Peters, “Three-dimensional reconstruction of vascular trees. Theory and methodology,” Medical Physics, vol. 23, no. 2, pp. 197–204, 1996.

S. M. Jorgensen, O. Demirkaya, and E. L. Ritman, “Three-dimensional imaging of vasculature and parenchyma in intact rodent organs with X-ray micro-CT,” American Journal of Physiology—Heart and Circulatory Physiology, vol. 275, no. 3, pp. H1103–H1114, 1998.

H. K. Tuy, “An inversion formula for cone-beam reconstruction,” SIAM Journal on Applied Mathematics, vol. 43, no. 3, pp. 546–552, 1983.

B. D. Smith, “Image reconstruction from cone-beam projections: necessary and sufficient conditions and reconstruction methods,” IEEE Transactions on Medical Imaging, vol. 4, no. 1, pp. 14–25, 1985.

F. Natterer, The Mathematics of Computerized Tomography, B. G. Teubner, Stuttgart, Germany, 1986.

F. Natterer, Mathematical Methods in Image Reconstruction, SIAM, 2001.

L.-Y. Wang, L. Li, B. Yan, C.-S. Jiang, H.-Y. Wang, and S.-L. Bao, “An algorithm for computed tomography image reconstruction from limited-view projections,” Chinese Physics B, vol. 19, no. 8, Article ID 088106, 2010.

M. Beister, D. Kolditz, and W. A. Kalender, “Iterative reconstruction methods in X-ray CT,” Physica Medica, vol. 28, no. 2, pp. 94–108, 2012.

H.-M. Zhang, L.-Y. Wang, B. Yan, L. Li, X.-Q. Xi, and L.-Z. Lu, “Image reconstruction based on total-variation minimization and alternating direction method in linear scan computed tomography,” Chinese Physics B, vol. 22, no. 7, Article ID 078701, 2013.

X. Jia, Y.-F. Lou, R.-J. Li, W. Y. Song, and S. B. Jiang, “GPU-based fast cone beam CT reconstruction from undersampled and noisy projection data via total variation,” Medical Physics, vol. 37, no. 4, pp. 1757–1760, 2010.

E. J. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Transactions on Information Theory, vol. 52, no. 2, pp. 489–509, 2006.

E. J. Candès, J. K. Romberg, and T. Tao, “Stable signal recovery from incomplete and inaccurate measurements,” Communications on Pure and Applied Mathematics, vol. 59, no. 8, pp. 1207–1223, 2006.

E. J. Candes and M. B. Wakin, “An introduction to compressive sampling,” IEEE Signal Processing Magazine, vol. 25, no. 2, pp. 21–30, 2008.

M. Chang, L. Li, Z. Chen, Y. Xiao, L. Zhang, and G. Wang, “A few-view reweighted sparsity hunting (FRESH) method for CT image reconstruction,” Journal of X-Ray Science and Technology, vol. 21, no. 2, pp. 161–176, 2013.

A.-L. Cai, L.-Y. Wang, H.-M. Zhang et al., “Edge guided image reconstruction in linear scan CT by weighted alternating direction TV minimization,” Journal of X-Ray Science and Technology, vol. 22, no. 3, pp. 335–349, 2014.

E. Y. Sidky and X. Pan, “Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization,” Physics in Medicine and Biology, vol. 53, no. 17, pp. 4777–4807, 2008.

H. Y. Yu and G. Wang, “Compressed sensing based interior tomography,” Physics in Medicine and Biology, vol. 54, no. 9, pp. 2791–2805, 2009.

H. Stark, J. W. Woods, I. Paul, and R. Hingorani, “Direct Fourier reconstruction in computer tomography,” IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 29, no. 2, pp. 237–245, 1981.

M. M. Seger and P.-E. Danielsson, “Scanning of logs with linear cone-beam tomography,” Computers and Electronics in Agriculture, vol. 41, no. 1-3, pp. 45–62, 2003.

B. P. Fahimian, Y. Zhao, Z. Huang et al., “Radiation dose reduction in medical X-ray CT via Fourier-based iterative reconstruction,” Medical Physics, vol. 40, no. 3, Article ID 031914, 2013.

A. Dutt and V. Rokhlin, “Fast Fourier transforms for nonequispaced data,” SIAM Journal on Scientific Computing, vol. 14, no. 6, pp. 1368–1393, 1993.

S. Matej, J. A. Fessler, and I. G. Kazantsev, “Iterative tomographic image reconstruction using Fourier-based forward and back-projectors,” IEEE Transactions on Medical Imaging, vol. 23, no. 4, pp. 401–412, 2004.

Y. Zhang-O'Connor and J. A. Fessler, “Fourier-based forward and back-projectors in iterative fan-beam tomographic image reconstruction,” IEEE Transactions on Medical Imaging, vol. 25, no. 5, pp. 582–589, 2006.

C. B. Li, An efficient algorithm for total variation regularization with applications to the single pixel camera and compressive sensing [M.S. thesis], Rice University, 2009.

Y. Xiao, J. Yang, and X. Yuan, “Alternating algorithms for total variation image reconstruction from random projections,” Inverse Problems and Imaging, vol. 6, no. 3, pp. 547–563, 2012.

Y. Wang, J. Yang, W. Yin, and Y. Zhang, “A new alternating minimizationalgorithm for total variation image reconstruction,” SIAM Journal on Imaging Sciences, vol. 1, no. 3, pp. 248–272, 2008.

C. Li, W. Yin, H. Jiang, and Y. Zhang, “An efficient augmented Lagrangian method with applications to total variation minimization,” Computational Optimization and Applications, vol. 56, no. 3, pp. 507–530, 2013.