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

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.


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][4][5][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][12][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][16][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][19][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, 2 Computational and Mathematical Methods in Medicine 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 sparseview 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 sparseview 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.

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 ∈ [0, /2] and ∈ [− /2, /2]. 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. Medicine   3 From the perspective of DFM, the observation equation in the Fourier domain can be expressed as follows:

Computational and Mathematical Methods in
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 tô( , V), which is derived from the 2D FFT of the reconstructed image in a certain angle. Therefore, the relationship is described as follows: . The equation shows an obvious equivalence corresponding to frequency projection ( ) witĥ( , V) 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 F 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 F and the variant f are the vector forms of Fourier sampling and objection function, respectively. Matrix F stands for the NUFFT of f. The Fourier transform F f in (5) and its adjoint F f 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 ‖f‖ TV denotes the discretization of the TV term and ‖f‖ TV = ∑ ‖ f‖ 1 . 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. 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 f and w 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 w-subproblem and the f-subproblem [29].
The w-subproblem can be written as follows: The w-subproblem is separable with respect to w, and problem (9) can be efficiently solved by using the shrinkage operator [30], which is expressed as follows: In addition, with the aid of w , the optimization of fsubproblem can be achieved by solving the following: w (f) is clearly a quadratic function, the gradient of which is expressed as follows: where + denotes the Moore-Penrose inverse of matrix . Theoretically, the exact minimizer can be used to solve the f-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 w-subproblem and the f-subproblem alternately. Therefore, solving the f-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 w * and f * , the multipliers are updated as follows: The optimized solution for (8) is attained by circularly applying (10) and (13) until (f, w ) is minimized jointly with respect to (f, w ).

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. (1) make 1D FFT of ( ) with respect to while "not achieved maximum iteration loops, " Do (3) compute f by (4) compute w by 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 sparseview because of its low computational cost, thus making this technique promising in practical applications.

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.   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. 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.
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.

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.
The reconstruction results are presented in Figure 6. 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.

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.