Nowadays, x-ray computed tomography (CT) has been widely used in clinical and industrial applications.^[1] The involved x-ray radiation has also attracted increasing concern. Low-dose imaging has become a major challenge in the CT field.^[2–4] The radiation dose in the data acquisition process can be reduced in two ways: by decreasing the current of the x-ray source and by reducing the number of projections. However, the details of an imaging object may be ignored with decreased scanning current. Therefore, the sparse scanning imaging process has become popular on the premise that it ensures reconstruction quality.

Research on sparse scanning has been conducted over decades. Given the low number of projections for the Tuy–Smith condition,^[5,6] which prescribes a number limit for an accurate reconstruction. Hence, to obtain ideal results through analytic reconstruction methods is difficult, such as the filtered back-projection (FBP)^[7] algorithms, which creates numerous artifacts because of the data inconsistency of sparse-view CT measurements. In fact, to the ill-posed problem, iterative reconstructions are currently the primary, widely used method of solving sparse projections. However, tolerating these methods is difficult because of their time-consuming and high-memory allocation in the spatial domain.

By contrast, the Fourier-based reconstruction method has advantages, such as high computational efficiency and low memory usage. These excellent properties are mainly due to fast Fourier transform (FFT). Various methods have been developed in accordance with central slice theorem.^[8] In 1987, Peng et al.^[9] derived fan-beam direct Fourier method by rebinning a fan to a parallel projection; this method has fewer operations than filtered convolution back-projection. In 2006, Fessler et al.^[10] presented a projection method for fan-beam tomography based on non-uniform FFT (NUFFT), which significantly improves the speed of iterative reconstruction. In 2014, Zhao et al.^[11] developed the generalized Fourier slice theorem and avoided the interpolation for fan beam rebinning.

However, the aforementioned Fourier methods are aimed at all angle projection data, and sparse-view reconstruction remains a challenge in the Fourier domain. Fortunately, CS-based iterative reconstruction^[12–14] draws significant interest because of its beneficial potential in under-sampled CT reconstruction. This relatively recent innovation in signal processing allows image recovery from projections fewer than those required by the Nyquist sampling theorem. Reconstruction research on compressed sensing mainly focuses on the development in the spatial domain. In 2008, Pan et al.^[15] proposed adaptive decent methods via a projection onto convex sets (POCS), which is conventionally used for algebraic reconstruction techniques (ART). Afterwards, the simultaneous ART (SART)^[16] that is combined with total variation (TV) regularization, namely SART–TV,^[17] improves reconstruction image quality and accelerates convergence. In 2013, the idea of alternating direction method (ADM) to solve constrained TV minimization was introduced by Zhang et al.,^[18] which is superior to other spatial algorithms in image quality. In 2014, Wang et al.^[19] reviewed basic conclusions and classical algorithms in sparse optimization in few-view reconstruction and predicted the future research direction of sparse optimization based few-view reconstruction for CT image reconstruction.

The theory of compressed sensing also offers new perspectives for Fourier-based iterative reconstruction. The related studies have also inspired the flourishing research in the compressive sensing area. In 2014, Choi et al.^[20] presented a Fourier-based CS technique, which provides a fast convergence rate for parallel beam CT reconstruction in sparse view. In 2015, Yan et al.^[21] proposed a reconstruction method of NUFFT-based iterative image reconstruction via alternating direction TV minimization (ADTVM) for sparse-view parallel beam CT. The above work concentrated on sparse view for parallel beam due to the central slice theorem which directly derived from parallel beam geometry.

For Fourier-based fan reconstruction, the process of rebinning is necessary because of the limit of central slice theorem. In general, because the number of angular sampling is smaller than full scanning, the missing data cannot be interpolated without aliasing.^[22] Therefore, the rebinning of fan sparse views is not straightforward suitable, and it needs to combine with other methods.

Inspired by previous work, we aim to present a promising contribution to the task of image reconstruction from fan sparse views in Fourier domain. To conveniently record the rebinning data between Fourier domain and image space, the selective matrix is presented to keep the known data constant. Combining the alternating direction total variation minimization (ADTVM) technique with NUFFT, a new method is proposed which is suitable for large-scale reconstruction because of high computational efficiency of FFT. The advantages of the proposed algorithm are identified based on the results of several group experiments.

This paper is organized as follows. Section 1 briefly reviews the basic CT reconstruction and the current research state of fan sparse-view image reconstruction. Section 2 describes the basic principles of the proposed method, which mainly include a selective matrix and a constrained TV minimization model. Section 3 demonstrates groups of typical experiment results, including numerical and real-data simulations. Finally, Section 4 presents the discussion and conclusion.

In this section, the new method is described in detail, which mainly includes the methods of central slice theorem, fan projection data rebinning, selective matrix, and constrained TV minimization reconstruction via NUFFT.

2.1. Central slice theorem

Central slice theorem (also called projection theorem) is outlined to review the foundation of Fourier reconstruction in the following. The basic geometry is shown in Fig. 1. In a compact support of spatial domain, the general formation for the line integral, known as the Radon transform of f(x,y) (objection function), is

Then, in the Fourier domain, the integral is expressed as follows:

where P_θ(ρ) is observed with the Fourier transform of the measured projection data P_θ (ρ), and ρ is the frequency variable of the Fourier transform. On the other hand, the two-dimensional (2D) Fourier transform of objection is denoted as

Combined with Eqs. (2) and (3), the following relationship can be obtained for the projection and frequency spectrum of objection:

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. Related Fourier transform of objection and projection.

Namely, the one-dimensional (1D) FFT P_θ (ρ) of the parallel projection is equal to 2D FFT F (u,v) of reconstruction objection in the direction perpendicular to the projection. In addition, the polar coordinate data cannot be directly transformed to spatial domain by IFFT and it needs transform to Cartesian coordinate which make plenty of streak artifacts due to the introduction of the interpolation. The adjoint of NUFFT can realize reconstruction from polar coordinate data and avoid interpolation form polar coordinate to Cartesian coordinate. Based on the latest NUFFT, we can design the following efficiency reconstruction.

2.2. Fan projection data rebinning

Given the limitation of parallel projection, other imaging geometries, including fan beam imaging geometry, needs a rebinning process before developing Fourier-based reconstruction methods. According to their geometric relationship (see Fig. 2), a straightforward approach would be to rebin every fan-beam ray into a parallel-beam ray.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. Geometry relationship of a fan beam and a parallel beam geometry for equal spacing.

From Fig. 2, the following coordinate relationship can be satisfied for rebinning:

where β, D, and L denote the rotation angle, the distance from the source to rotation axis, and the distance from the source to the linear detector, respectively. Moreover, s is the signed distance from the isoray along the linear detector. When these variables satisfy the preceding formula, p(t,θ) is equal to R(s,β). However, this rebinning approach is not preferred because rebinning requires data interpolation with involving coordinate conversion, which introduces error. The idea of rebinning is feasible, but it requires a combination of other methods if the reconstruction result is sufficiently accurate.

2.4. Constrained TV minimization reconstruction via NUFFT

Let f denote the two-dimensional image, then sparse model can be linearly expressed as follows:

where F_N and s represent the NUFFT operator and selective matrix, respectively. Then, is the 1D Fourier transform of the sparse parallel P, which is a result of the rebinning of sparse fan P_f.

The preceding Fourier model can be concisely written as follows by letting F_s = s F_N:

To solve this linear and underdetermined equation, we specify a TV minimization algorithm that considers the reconstruction a task of finding the best solution to the following optimization problem:

where denotes the discrete gradient of f at pixel i, and D_i represents the differential operator along the direction i of the image. Instead of minimizing the mentioned TV model directly, we consider an equivalent variant of Eq. (9),

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

Its corresponding augmented Lagrangian function of Eq. (10) is

Given that equation (10) remains a convex problem, the global convergence theorem can guarantee convergence, while applying the augmented Lagrangian method to deal with it. By letting and f* represent the true minimizers of Eq. (11), the following update formula of multipliers^[23,24] is obtained:

The task of solving Eq. (9) is equal to minimizing L(w_i, f) efficiently at each iteration. In this study, the ADM is used to solve the minimization of Eq. (11). The ADM approach decouples the augmented Lagrange function into two subproblems, namely, w_i-subproblem and f-subproblem.

Suppose that f_k and w_i,k, respectively, denote the approximate minimizers at the k-th iteration which refers to the inner iteration, while solving the subproblem. If f_j and w_i,j are available for all j = 0,1,…,k, then w_i,k+1 can be attained by

The w-subproblem is separable with respect to w_i, and equation (13) can be efficiently solved with the shrinkage operator.^[25] The closed form solution is expressed as follows:

In addition, with the aid of w_i,k+1, the f-subproblem can be achieved by solving the following equation:

where Q_k(f) is clearly a quadratic function, and its gradient is

Force d_k(f) = 0 and the exact solution for Q_k(f) is

where M⁺ stands for the Moore–Penrose pseudo inverse of matrix M. In theory, accepting the exact minimizer as the solution of the f-subproblem is ideal. However, computing the inverse or pseudoinverse at each iteration is too costly to implement numerically. In general, we use a robust and efficient algorithm, namely, nonmonotone alternating direction algorithm (NADA),^[26] to solve Eq. (17).

The optimized solution for Eq. (11) is attained by circularly applying Eqs. (14) and (17) until L(w_i, f) is minimized jointly with respect to (w_i, f). The following equations show the (k + 1)-th iteration of the proposed method:

2.5. Algorithm of the overall framework

Thus, all issues in the process of handling the subproblems have been settled in the preceding sections. In light of all the preceding derivations, the new algorithm for fan sparse-view reconstruction is stated as follows.

Algorithm 1: υi = vi,0, w = wi,0 is initialized for all i, fan sparse projections Pf, λ, βi > 0 are inputted, and k = 0 is set. 1) Data rebinning is performed, then the selective matrix is built. 2) 1D FFT of P(t, θ) is performed with respect to t by While “not achieved maximum iteration loops,” Do 3) Compute f with 4) Compute w with 5) Update υi with 6) k ← k + 1 End Do

It is noticeable that this property of iteration mainly depends on alternating direction method. Therefore, the convergence of algorithm is inevitable for the proof convergence of ADM.^[27] The Fourier-based algorithm here is validated that the complexity decreases from O(N³) to O(N²logN) for NUFFT and its adjoint. The proposed non-uniform Fourier reconstruction based on TV regularization algorithm for sparse fan projection is called NUTVM in this paper.

We implement the proposed NUTVM in MATLAB environment to demonstrate its excellent performance. It is also evaluated with both digital simulation data and real CT data. All experiments are performed on an HP Z820 workstation with Intel Xeon E5-2620 dual-core CPU (2 GHz) and 24 GB memory.

3.1. Digital simulation study

We perform the simulation studies with a 2D modified Shepp–Logan phantom, which is generated according to the definition of the ellipse phantom image, to validate the performance of our proposed method. Circular equally spaced fan beam geometry is adopted. In the experiment, the detector elements are linearly spaced within the line detector where the size of one pitch is assumed to be 0.127 mm × 0.127 mm. The scanning and reconstruction parameters in the experiment are listed in Table 1. In the following experiments, we perform a series of comparisons, including ASD-POCS, ADTVM, and the proposed NUTVM, under the preceding configuration.

Table 1.

Table 1.

Table 1. Parameters in the simulation of a sparse fan scan. . Parameters Configuration Detection elements 512 Source to axis distance 400 mm Source to detection distance 600 mm Views of projection data 60 Projection data 512×60 Reconstruction size 256×256 pixels Pixel size 0.127 mm× 0.127 mm

Table 1. Parameters in the simulation of a sparse fan scan. .

To further evaluate the reconstruction results, the root mean square error (RMSE) is introduced, which can be written as

where f_t and f₀ denote the reconstruction and reference image of N pixels. We performed the experiments to further check the capability of the proposed algorithm by reconstructing the images from noisy projections in all tests. In MATLAB, noise is generated with a Poisson model^[28]

where N₀ denotes the x-ray initial intensity, and P is the normalized projection in real space.

The image was reconstructed with ASD-POCS, ADTVM, and the proposed method, and their results are presented in Fig. 4 with N = 5.9×10⁵. Then, 400 iterations are performed for each algorithm.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. Image reconstructions of the Shepp–Logan phantom in 60 view scanning. Display window is (0.1, 0.5). (a) The original image. (b) The image after applying the SART–TV algorithm at 400 iterations. (c) The reconstructed image via the ADTVM algorithm at 400 iterations. (d) The reconstruction image after using the NUTVM algorithm at 400 iterations.

In Fig. 4, the corresponding profiles of these images along the central horizontal and vertical rows are shown for different algorithms.

RMSE is used as an evaluation criterion for three different algorithms with iteration time. The result is shown in Fig. 6. It is obvious that the RMSE of ASD-POCS algorithm is zero in the first twenty seconds because the time of each iteration is twenty seconds. So are ADTVM and NUTVM algorithm. Moreover, the accuracy and running time of each reconstruction method at the same iteration number are presented in detail in Table 2 for comparison.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. (a) ASD-POCS of the vertical rows, (b) ADTVM of the vertical rows, (c) new method of the vertical rows, (d) ASD-POCS of the horizontal rows, (e) ADTVM of the horizontal rows, (f) new method of the horizontal rows.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. RMSEs with iteration time for three different algorithms.

Table 2.

Table 2.

Table 2. Accuracy and the running time and memory allocation of different methods. . Iteration number 200 400 800 ASD-POCS RMSE 0.0036 0.0045 0.0082 Running time/s 3455.97 7092.54 19732.58 Memory allocation/GB 15.425 ADTVM RMSE 0.0018 0.0024 0.0035 Running time/s 354.95 843.96 2200.48 Memory allocation/GB 14.889 NUTVM RMSE 0.0250 0.0020 0.0025 Running time/s 42.54 84.23 186.75 Memory allocation/GB 0.134

Table 2. Accuracy and the running time and memory allocation of different methods. .

From the reconstruction images, the corresponding profiles and RMSE, we can see that there is little difference between the reconstructed image using the NUTVM algorithm and the reference image, which demonstrate the correctness of the proposed algorithm. The reconstruction quality of NUTVM algorithm is an improvement over that of the ASD-POCS algorithm. In terms of reconstruction time, the new method is obviously superior to the other two contrastive algorithms. However, the reconstructed quality worsens for the existence of noise when the iteration number is over a certain number.

3.2. Reconstruction via real data

We performed the experiments to reconstruct a slice of the head model from real data to further check the capability of the NUTVM algorithm. The scanning and reconstruction parameters are listed in Table 3. The detector elements are equidistantly spaced (0.148 mm) from one another.

Table 3.

Table 3.

Table 3. Parameters in the sparse-view fan scanning. . Parameters Configuration Detection elements 3286 Source to axis distance 678 mm Source to detection distance 1610 mm Views of projection data 90 Projection data 3286 × 90 Reconstruction size 256×256 pixels Pixel size 0.582 mm×0.582 mm

Table 3. Parameters in the sparse-view fan scanning. .

Image reconstructions through ASD-POCS and ADTVM algorithms and the proposed algorithm are shown in Fig. 7. The numbers of iterations for the three algorithms are 200 and 400.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. Reconstruction results of ASD-POCS, ADTVM, and NUTVM. (a) ASD-POCS (at 200 iterations), (b) ADTVM (at 200 iterations), (c) New method (at 200 iterations), (d) ASD-POCS (at 400 iterations), (e) ADTVM (at 400 iterations), (f) New method (at 400 iterations).

Then the running time and memory allocation of each reconstruction method are presented in detail in Table 4 for real data.

Table 4.

Table 4.

Table 4. The running time and memory allocation of different methods. . Algorithm ASD-POCS ADTVM NUTVM Running time/s 6676.4 1119.6 169.1 Memory allocation/GB 8.46 6.87 0.16

Table 4. The running time and memory allocation of different methods. .

The acquired reconstruction results via real data clearly show that the quality of the reconstructed image improved as the number of iterations increased. Under the same number of iterations, the reconstruction results of ADTVM are superior to those of ASD-POCS, particularly in terms of the high-frequency information showing the image detail or volatile part. Compared with ADTVM, the results of the new method are approximate from the view of the image detail.

The simulated data agree well with the theoretical analysis for the 2D Fourier propriety to the reconstructed object. According to the preceding results, we examine the proposed iterative algorithm, which shows potential improvement. During the numerical simulation experiment, we determined that error from rebinning mainly depends on the choice of interpolation. In general, choosing an intricate interpolation method, such as polynomial function, helps in improving the quality of the reconstruction image to full angle scanning. However, the lack of projection data leads to an increased number of errors because the intricate interpolation method must use global data. Thus, excellent results must be realized by a simple interpolation, such as linear interpolation. According to the RMSE curve of the three different algorithms in the numerical simulation, the error stops decreasing when it arrives at the iteration number for the existing noise. This phenomenon is in accordance with the real reconstruction data because of the existing irregular noise. In the experiment, memory allocation by NUTVM is the smallest among the other contrastive algorithms. This preference mainly depends on FFT, encompassing low computation complexity and high computation time.

Due to the lack of sufficient measurements, the reconstruction procedure can be treated as an ill-posed inverse problem for sparse view, which generally has no unique solution. Based on the continuity of piecewise constant, the TV minimization model can uniquely reconstruct object by high probability. However, it is known that transformation between different data domains is the most computationally expensive part of the whole algorithm. Fortunately, NUFFT and adjoint NUFFT tackle with the transform, which makes this algorithm highly computationally efficient. In addition, unlike ART-based algorithms, the new method does not have to model the system matrix and the data consistency can be retained. The proposed method may have a greater capacity for reducing artifacts and reconstructing better quality images.

A chief contribution of this study is the proposal of a new, efficient TV minimization scheme with NUFFT for fan sparse projections. The validity of the proposed algorithm is verified by conducting numerical simulations and real-data experiments. The advantages of NUTVM indicate that the proposed algorithm is a good choice for fast, low-dose reconstruction in fan beam sparse scanning. Furthermore, it has potential in extending to the cone beam reconstruction in the Fourier domain because of fan basic properties.