Accelerated Compressed Sensing Based CT Image Reconstruction

In X-ray computed tomography (CT) an important objective is to reduce the radiation dose without significantly degrading the image quality. Compressed sensing (CS) enables the radiation dose to be reduced by producing diagnostic images from a limited number of projections. However, conventional CS-based algorithms are computationally intensive and time-consuming. We propose a new algorithm that accelerates the CS-based reconstruction by using a fast pseudopolar Fourier based Radon transform and rebinning the diverging fan beams to parallel beams. The reconstruction process is analyzed using a maximum-a-posterior approach, which is transformed into a weighted CS problem. The weights involved in the proposed model are calculated based on the statistical characteristics of the reconstruction process, which is formulated in terms of the measurement noise and rebinning interpolation error. Therefore, the proposed method not only accelerates the reconstruction, but also removes the rebinning and interpolation errors. Simulation results are shown for phantoms and a patient. For example, a 512 × 512 Shepp-Logan phantom when reconstructed from 128 rebinned projections using a conventional CS method had 10% error, whereas with the proposed method the reconstruction error was less than 1%. Moreover, computation times of less than 30 sec were obtained using a standard desktop computer without numerical optimization.


Introduction
Compared to conventional radiography, CT results in a relatively large radiation dose to patients, which is of serious long-term concern in its potential for increasing the risk of developing cancer [1,2]. As a result, low dose CT imaging that maintains the resolution and achieves good contrast to noise ratio has been the goal of many CT developments over the past decade. However, low dose CT images reconstructed with conventional filtered back projection (FBP), which directly calculates the image in a single reconstruction step, suffer from low contrast to noise ratios. Iterative reconstruction approaches, namely, Algebraic Reconstruction Technique (ART) [3][4][5] and statistical iterative reconstruction (SIR) [6,7], have been proposed to improve the reconstruction quality and to decrease image artifacts. The iterative algorithms improve the quality by considering more accurate models for the CT images and geometries. However, they significantly increase the computational complexity, compared to the FBP based methods.
Iterative reconstruction methods have progressed with the introduction of compressed sensing (CS) [8,9]. Such methods are capable of reconstructing high quality images from a substantially smaller number of views than those needed in FBP [10], thereby permitting the use of a much lower dose scanning protocol than that needed in conventional reconstruction methods. However, conventional CSbased CT reconstructions are computationally expensive and the statistics of CT measurements are not usually incorporated in the problem formulation [11][12][13][14][15][16].
In this paper, we propose a fast weighted CS-based CT reconstruction algorithm, the weights of which are direct consequences of the geometry and the CT statistics. The first part of this paper leads to the proposed weighted CS formulation, which is solved by a computationally efficient method discussed in the second part.

CS-Based CT Reconstruction and Its Challenges
Compressed sensing prescribes solving the ℓ 1 optimization problem that can be represented bŷ or other similar forms, for example, the ℓ 1 norm of the image gradient, such aŝ to recover a sparse signal from few samples. In these equations, acts as a regularization parameter specifying a trade-off between the image prior model and the fidelity to observations, A is the measurement matrix, x is the column vector representation of the desired image ( ), y is the measured data, is a sparsifying transform, ‖ ‖ = (∑ | | ) 1/ , and TV denotes the total variation TV( ) = ∑ √(∇ ) 2 + (∇ ) 2 , where ∇ and ∇ are the first derivatives in the and directions of the desired image.
The main challenge in solving these optimization problems within a reasonable amount of time arises from the size of the measurement matrix A. Currently, in most available CS-based reconstruction methods, the measurement matrix A is a projection matrix which models the rays going through the patient. To reconstruct a 512 × 512 pixel image from 900 sensors and 1200 projection angles, A would be a 1080000 × 262144 matrix. Although this matrix is sparse, each iteration typically requires two multiplications by A and A T , resulting in a very significant increase in the computation burden for reconstructing a 512 × 512 image [11,12] as compared to FBP based methods. To enable the CS-based CT reconstruction to be done in a reasonable computation time, GPU based algorithms have been proposed [17].

Complexity Reduction Using the Pseudopolar Fourier Transform (PPFT).
To reduce the computational burden on the Radon transform, the central slice theorem (CST) or direct Fourier reconstruction (DFR) has been used [18]. This relates the 1D Fourier transform of the projections to the 2D Fourier transform of the image. Such a method requires the interpolation of polar data onto a Cartesian grid followed by an inverse FFT on the same grid to reconstruct the CT image. Since interpolation does not have a known analytical adjoint, its use in iterative algorithms is not a practical option. In addition, inclusion of a gridding and regridding step at each iteration increases the overhead computational complexity. This problem has been extensively studied in non-Cartesian magnetic resonance imaging reconstruction algorithms [19].
An equally sloped tomography (EST) method was originally proposed for electron beam tomography [20][21][22] to improve the DFR-based algorithms. EST is an iterative method that makes use of the pseudopolar Fourier transform (PPFT) [23]. It calculates the Fourier coefficients of an image directly on pseudopolar grids, which contain two types of samples: basically horizontal (BH) and basically vertical (BV), as can be seen in Figure 1. To reconstruct an × image from its PPFT coefficients, 4 2 samples are needed (2 samples on 2 equally sloped radial lines). A fast algorithm has been proposed by Averbuch et al. [23] to calculate the PPFT and its adjoint with complexity of ( 2 log ). This algorithm can then be used to implement a fast and efficient 2D Radon transform on the equally sloped radial lines. The PPFT has three important properties which makes it a good alternative to conventional DFR methods: (1) it is closer to a polar (equiangular line) grid than to a Cartesian grid, which significantly decreases the gridding error, (2) it has both a fast forward and a fast backward calculation algorithm [23], which enables our proposed algorithm to avoid the regridding step used in iterative non-Cartesian Fourier based reconstruction methods, and (3) it has an analytical adjoint function. As a result, it can efficiently be used in iterative algorithms, including compressed sensing [24,25]. However, it should be noted that Fourier-based reconstruction algorithms, for example, ESR and DFR, are only valid for parallel X-ray projections.
A major objective of this paper is to accelerate the CSbased CT reconstruction by decreasing the CS complexity using PPFT-based Radon transform proposed in [26]. The application of the proposed method is extended to equiangular parallel and nonparallel geometries using rebinning.

Rebinning Process.
To enable use of the PPFT-based Radon transform for nonparallel geometries, the projected rays must first be transformed to parallel beams [27]. This requires two interpolation steps. At first, projections are interpolated on equally sloped radial lines on the following angles: This step makes use of the following relationships between fan and parallel beams: where , , = BH ∪ BV , and are geometry parameters defined in Figure 2. Moreover, R( , ) is the fan beam projected data and ( sin , + ) is the corresponding rebinned parallel ray. In the second interpolation step, the  radial samples are interpolated on the pseudopolar grids shown in Figure 1(b). To reduce the interpolation error, these radial lines are zero-padded and the 1D Fourier transforms of the zero-padded radial lines are interpolated. This is equivalent to oversampling in the Fourier domain, which makes the interpolation error manageably small.

Measurement Noise and CS-Based CT Reconstruction.
The measurement noise in CT scanners can best be modeled by a Poisson distribution [27], while the noise that is considered in classical CS formulations, such as those given by (1) and (2), is white additive Gaussian noise [8,9]. Therefore, to enable a more accurate low dose CS-based CT reconstruction, the classical CS formulations should be modified. To address this problem in prior studies different approaches have been used. For instance, in [14,28,29] to account for the statistical properties associated with low-dose measurements, an iterative SIR based technique followed by TV denoising was used. The penalized weighted least squares (PWLS) formulation of statistical CT reconstruction was used in [30,31] to improve the quality of reconstructed images. This PWLS formulation can be characterized as a weighted ℓ 1 minimization problem as proposed by Candès et al. [32], who showed that by using appropriate weights the quality of the recovered signal can be improved. In this modified CS formulation, the weights could account for the statistical characteristics of the signals.
In this paper, a fast CS-based CT reconstruction is proposed using a pseudopolar based Radon transform to decrease the complexity of the CS recovery. The method requires the projections to be rebinned on equally sloped lines. We propose a weighted CS formulation in the framework of statistical CT image reconstruction algorithms. The weight, denoted by error adaptation weight (EAW), is a function of the rebinning interpolation error [25] and the Poisson noise of the CT projections [31] calculated from a maximum a posteriori (MAP) model of CT reconstruction, described in next section.

Maximum a Posteriori (MAP) Model of CT
X-ray projections of the parallel beam CT can be expressed as the Radon transform of the object. The Radon transform is defined as [33] ( , ) which is the integral along a ray at angle and at the distance from the origin, ( , ) is Dirac delta function, and ( , ) is the object attenuation at ( , ). However, this is not what the scanners directly measure. Scanner detectors measure the number of photons that hit the detector, ( , ), which is usually modeled by Poisson distribution with expected value of ( , ) [6,27]. The relation between ( , ) and ( , ) is ( , ) = −log( ( , )/ ), where is the number of radiated photons from the X-ray source. It should be noted that ( , ) is usually corrupted with two kinds of noise: electrical noise of the detectors (with variance of 2 ) and the photon counting noise (observed counts are drawn from a Poisson distribution of mean ). If we consider the discrete formulation in which y denotes the vectorized ( , ), x denotes the vectorized ( , ), and A is the projection matrix, using the second order Taylor series expansion of the Poisson distribution, the log likelihood of the measurements is given by [34,35] in which (y 3 ) is a function which depends upon measured data only and is a diagonal matrix. For the purpose of the MAP estimation (y 3 ) may be ignored since it does not depend on x. Ignoring this term, (6) describes a simplified CT model that can be written as in which n is Gaussian distributed noise with a covariance matrix −1 and , the th diagonal element of , is proportional to the detector counts, corresponding to the maximum likelihood of the inverse of the variance of the projection measurements, that is, to 1/ 2 . The th measured projection y is given by where y is noiseless and follows the Poisson distribution with 2 = . As a result the variance of projection data can be estimated from 2 Using as an unbiased estimation of , the diagonal elements of can be expressed as To reconstruct the image from the projections, the MAP estimator can be used: Here ℎ( ) = log (x) acts as a penalty function, which is used to statistically model the wavelet coefficients distribution and the piecewise constant (locally constant) nature of CT images.

Piecewise Constant and Sparsity of the CT Images in MAP
Model. Many studies [36,37] have shown that the wavelet transform of a variety of images, = x, can be modeled by generalized Gaussian distribution (GGD), that is, by where is the wavelet coefficients of the image x = , is inverse wavelet transform, and are the parameters of the GGD, and ( , ) is the normalization parameter. It should be noted that when = 1, the GGD is equivalent to Laplacian distribution and when = 2, it describes a Gaussian distribution. Using (6), (10), and (11), the MAP model for CT images can be expressed aŝ where ‖y−Ax‖ 2 = (y−Ax) (y−Ax), is a function of GGD parameters, and is typically in the range 0 < ≤ 1. Another prior on ( ) is that the objects being imaged are piecewise constant. A variation distribution can be used to describe the piecewise constant functions [38]. If ( ) = ∑ =1 Ψ ( ) is a piecewise function spanned by the roof-top basis Ψ ( ), Computational and Mathematical Methods in Medicine 5 the following class of probability distribution can be used to describe it: where > 0, 0 = +1 = 0, , is normalizing factor, and [ 1 , . . . , ] is a R -valued random vector. When = 1, this yields the total variation norm. Using = 1 in (13), (10) becomes the following MAP model for CT images: As can be seen, (14) and (12) are generalized forms of the CS models given by (2) and (1), respectively. It has been shown that the quality of the reconstructed image can be improved by combining the sparsity and total variation penalty terms [39]. Therefore, both penalty functions were used in the proposed method to reconstruct CT images from the undersampled data, that is, from where 1 and 2 are regularization parameters.

Proposed CS Formulation: Measurement Noise and Interpolation Error
In Section 3, it was shown that the MAP estimator of CT is a form of the weighted CS problem given by (15), in which the weight is a function of noise variance and is denoted by in (9). To avoid the computational burden associated with the huge projection matrix of CS-based CT, our proposed algorithm makes use of the fast PPFT-based Radon transform described in Section 2.1. This not only accelerates the computations by reducing the computational complexity, but also substantially reduces the gridding error and eliminates the regridding step (since it has a fast backward calculation algorithm, there is no need to regrid the updated image to pseudopolar grids at each iteration).
As described in Section 2.1.1, the proposed method is generalized for use with nonparallel geometries by rebinning the X-ray beams onto equally sloped radial lines as given by (3). The rebinning step induces interpolation error to the measured data, which propagates in each iteration of CS-based CT reconstruction. This problem has not received much attention in the literature. Fahimian et al. [24] proposed an EST method for reconstructing fan beam and helical cone beam images, in which they overcome the rebinning interpolation problem at each iteration by using a nonlocal total variation minimization smoothing step. In the method proposed by Hashemi et al. [25], an ℓ 2 -TV optimization scheme was used to reconstruct the CT images from fan beam projections. To compensate for the interpolation error, a confidence matrix was added to the CS scheme, enabling control of the error propagation in successive iterations.
To control the interpolation error, we make use of a MAP model of the CT reconstruction process. Denoting the variance of interpolation error by , the variance of the measurements is Using (8), the effects of both the noise variance and interpolation error can be lumped together into the form of an error adaptation weight (EAW), denoted by a diagonal matrix C with diagonal elements: Assuming a nearest neighbor interpolation, the interpolation error is linearly dependent on the distance between the desired and the measured grids; that is, = × in which ∈ [0, ∞) models the interpolation distance. Consequently, if the dose at each projection is high enough to ignore the electric noise , the EAW can be rewritten as Using this definition, our proposed CS formulation to reconstruct the CT images can be expressed bŷ This implicitly models the effect of polar and nonparallel projections to pseudopolar gridding, which in turn affects the noise in the data. This error is considered to be linearly dependent on the interpolation distance, while it is typically smaller using more accurate interpolation methods, that is, in Kaiser-Bessel based interpolation [19]. Therefore, the proposed CS formulation can be thought of as a minimization of an upper bound of the pseudopolar rebinning error. In addition, C in this formulation acts similar to a Jacobi preconditioner and therefore could accelerate the convergence rate [40].

Calculation of .
The value of represents the error of the interpolated samples. If the interpolated sample is close to the original measurements, the value of is small and the confidence about the interpolated value is high. If the angular distance of the measured data from the interpolated line is more than the angular difference of the equally sloped lines, the interpolation error is considered to be high ( → ∞): this follows from the fact that the distance of points on the line from the true measured values is maximal and therefore the error is maximal. Using (18), this condition corresponds to → 0. The closer the equally sloped lines are to the rays on which the measurements are made, the smaller the interpolation error will be, so that s on that line get closer to zero. Finally, if the desired equally sloped rays are exactly on the polar lines, the interpolation error is zero, which is equivalent to = 0. This process is illustrated in Figure 3.
In practice, (1/(1 + )) can be estimated by rebinning an all-ones matrix with the same size as the measured data onto  the equally sloped radial angles followed by an interpolation on the pseudopolar grid. Note that this has to be calculated only once before the reconstruction.

Solving the Proposed CS Formulation
To solve the proposed formulation, a fast composite splitting algorithm (FCSA) [41][42][43] is used to decompose (19) into two simpler subproblems given bŷ in which y is the measured data interpolated/rebinned on the equally sloped lines, A is the PPFT-based Radon transform, and A is its adjoint. By calculatingx 1 andx 2 , the FCSA method proposes that the solution to the problem can be obtained by a linear combination of the solutions of the two subproblems; that is,x in which Δ = 2 /( 1 + 2 ) is a function of the values of the objective functions of the two subproblems. Each of these subproblems can be solved by a subgradient-projection based method [44]. The pseudocode of the proposed recovery is shown in Algorithm 1, in which prox{ ( ), } = arg min ( ) + (1/2)‖ − ‖ 2 2 . To findx 1 , the optimization problem in step (2) of this algorithm can be solved by a wavelet soft thresholding algorithm [36]. Moreover, to calculatex 2 in step (3), the split Bregman TV based denoising algorithm as proposed in [45] was used. Finally, to estimatê x in the th iteration, (21) was used.
In the proposed CT reconstruction algorithm, summarized in Figure 4, the Daubechies wavelets with four vanishing moments in 5 levels are used as the sparsifying transform . The regularization parameters are manually tuned to 1 = 0.05‖ (A y)‖ ∞ and 2 = 1 × 10 −3 TV(A y).

Simulation Methods
Fan beam simulations were performed using a Shepp-Logan phantom available in MATLAB (MathWorks, Massachusetts, USA), a custom made phantom that mimics different cardiac plaques and a clinical patient. This study was approved by our institutional (Toronto General Hospital, Toronto, ON, Canada) review board and individual patient consent was waived. X-ray projections of the phantoms and the patient were taken using a Toshiba Aquilion ONE © scanner (Toronto General Hospital, Canada). The scanner gathers data from 900 projection angles in each 360 ∘ rotation. To be compatible with the available hardware, when the images were reconstructed from fewer than the 900 projections, the projection views were selected equiangularly. For all the scan protocols, the X-ray tube current-exposure time product was 50 mAs and the peak voltage was 120 kV. This current/voltage is high enough to ignore the electric noise in the simulations. Data from the central row of a volumetric scan on one single rotation served as the fan beam data.

Results and Discussion
In this section, the performance of the proposed method is evaluated in terms of the reconstruction time and accuracy. In addition, the effect of the EAW is evaluated in reducing the influence of the interpolation error and the Poisson measurement noise.

Reconstruction Time Acceleration.
Compared to the other CS-based reconstruction techniques, the primary improvement of our proposed method is the major reduction in the computational burden. The reconstruction time of the proposed method is compared with a conventional ART-TV based CT reconstruction algorithm. The convergence of this method is justified by projection on convex set (POCS) algorithm. As an example in [46], an ART-TV based CT reconstruction method is proposed denoted by adaptivesteepest-decent POCS (ASD-POCS), which has been used by many other researchers [14,29,47,48]. Here, we use a simple method in which the updates are calculated using an ART based method described below. This step is followed by a TV minimization step to project the updated image on a piecewise constant space. To solve the y = Ax problem in ART step, a randomized Kaczmarz algorithm is used [4]. If Initialize: 1 , 2 , C, 1 Algorithm 1: Algorithm used to solve (19). y = Ax is a linear system of equations and x 0 is an arbitrary initial approximation to the solution, randomized Kaczmarz applies the following updating step at each iteration: where ( ) is chosen randomly from the set {1, 2, . . . , }, with probability proportional to ‖ ( ) ‖ 2 2 , is the th row of A, and ⟨⋅, ⋅⟩ is the inner product of two vectors. Using this algorithm followed by the split Bregman TV minimization (denoted by SBROF in Algorithm 2), the ART-TV based method is described in Algorithm 2. In our simulations, the number of inner iterations used in Kaczmarz algorithm is 10 and the number of outer iterations is 50. Figure 5 compares the recovery time using (1) filtered back projection (FBP), (2) the proposed method, and (3) the ART-TV based method described in Algorithm 2 (the ART method implementation is based on the codes provided in algebraic iterative reconstruction methods (AIR Tools) and Tomobox packages). It can be seen that for a 512 × 512 image the recovery time for the proposed method is approximatively 10-30 sec, using MATLAB on an Intel(R) Initialize: 1 = 0, inner, outer, , tol Core(TM) i5 (3.2 GHz) CPU desktop PC with 16 GB of RAM. Figure 6 shows the phantom reconstructed by ART-TV described in Algorithm 2 from 200 projections. The normalized reconstruction error ‖x −x‖ 2 /‖x‖ 2 is 2 × 10 −2 .

Interpolation Error and Noise Correction Using EAW.
Equiangular fan beam projections of the Shepp-Logan phantom were computed on 128 projection views. This data was then rebinned to parallel rays on equally sloped angles as described by (3). without using EAW weights and they were solved using an iterative soft threshold-based method [49], and (3) using the proposed method with EAW. Note that the lower row consists of an expanded view of the marked central region, from which it is clear that the proposed method yields results very close to that shown in Figure 7(a). Based on the same phantom, Figure 8 compares the accuracy of the reconstruction error for all three methods as the number of projections is varied from 50 to 1024. Both of these figures show that the recovery accuracy is improved significantly by the inclusion of EAW to correct for rebinning errors.
In particular, Figure 8 shows that the use of more than 300 projections for a 512 × 512 image does not significantly affect the reconstruction accuracy. Since the purpose of this particular study was to examine effects of EAW inclusion on rebinning interpolation error, noise was not included in the simulations.
To examine the effects of EAW on measurement noise, 128 equiangular projections through the Shepp-Logan phantom were computed and Poisson noise was added to the projections. Figure 9 shows the effect of EAW inclusion for different input peak signal-to-noise ratios (PSNR) on the PSNR of the reconstructed images. Images were reconstructed with the proposed method once with including EAW that is calculated by (9) and once without EAW. As can be seen, the PSNR is improved when the input noise is larger (small input PSNR) and its effect is less when the noise is low (larger input PSNR). The input PSNR was measured from the FBP reconstructed images.
Results obtained from the custom fabricated cardiac plaque phantom are shown in Figures 10(a)-10(c), which provides a comparison of reconstructions using FBP from 900 projections with the results obtained from 200 equiangular projections. It should be noted that the number of projection is chosen based on Figure 8, which shows that the error of the images reconstructed from 200 projections is in an acceptable range. Qualitative image evaluation was performed using a continuous linear scale of 1-5 (Excellent). Statistical analysis was performed using the paired Student -test, which showed no statistical difference in image quality between the cardiac phantom images reconstructed from 200 projections using the proposed method and the full projection FBP images ( = 0.05).
Finally, Figures 10(d)-10(f) show reconstructions of a chest CT scan from a hospital patient using FBP from 900 projections and the proposed method with 200 projections. The image reconstructed with the proposed method has almost the same quality as the FBP. Some small details in bony structures and flat lung regions are removed or have decreased contrast in the image reconstructed by the proposed method. However, all the important details are preserved. Thus, in comparison with FBP, the proposed method uses about four times less projections; that is, the radiation dose is decreased by more than a factor of four (78% dose reduction).
Current commercial CT scanners are unable to switch their X-ray sources on and off fast enough to achieve the proposed equiangular simulations. To overcome this problem, the mask used in Figure 11(a) is used that addresses this concern by turning the X-ray source off in the black areas over a range of angles and then turning it on in the white areas. To reconstruct high quality CT images scanned by this protocol, the reconstruction algorithm is modified by stacking the similar patches into 3D stacks [50]. Applying 3D wavelet thresholding/shrinkage on the 3D stacks of the similar patches increases the sparsity of the wavelet coefficients, which in turn improves the image reconstruction. The similar patches are selected from overlapped 15 × 15 neighborhoods and the patches are 6 × 6.

Conclusion
It has been shown that CT reconstruction can be statistically modeled as a weighted compressed sensing optimization problem. Our proposed weighted CS-based CT reconstruction algorithm was derived from the MAP model of CT imaging, considering the sparsity of the wavelet coefficients and piecewise constant nature of the CT images. Subsequently, a fast CS recovery method was proposed in which the pseudopolar based Radon transform was used as the measurement function to reduce the computational complexity. Moreover, to reconstruct CT images from nonparallel projections, rebinning to parallel beams was used. To remove the interpolation error caused by rebinning and the measurement noise, a weighting approach (EAW) was proposed. This enabled CT images to be reconstructed from a reduced number of projections. It was shown that using EAW improves the reconstruction quality substantially. For instance, a lung CT image was reconstructed with 78% lower dose but the same diagnostic quality as the image reconstructed by FBP from full data. The greatly reduced computational complexity of the proposed algorithm enabled a 512 × 512 image to be reconstructed in less than 30 sec on a desktop computer without numerical optimizations. Thus, our proposed method may be among the first CS-based CT reconstruction methods whose computational complexity is sufficiently small to enable low-dose image reconstructions to be performed without using either time-consuming computations or a complex computational system. Finally, it should be noted that the proposed method can be extended to 3D geometries by using the approaches described in Appendix A.

A. Extending the Algorithm to 3-Dimensional Geometries
As discussed in Section 2.1.1, to be able to use our proposed method for nonparallel geometries, the projections should be transformed to the Radon transform data. Simple rebinnging algorithms can be used for fan beam and helical cone beam projections to redistribute the diverging beams into parallel beams, for which the Radon transform and the 2D projections are identical. However, this transformation is more complicated in 3D cases, since the 3D Radon is very different from the 3D X-ray projections.
In this appendix, we first describe a method, using which the helical cone-beam projections can be transformed to fan bean projections. This rebinning algorithm can be used in our proposed method to reconstruct 3D images from the helical cone beam projections.
In the second part, we describe a method that converts the 3D cone beam projections to 3D Radon transform. The 3D Radon transform of the object can be used in our proposed method along with 3D PPFT transform to reconstruct the cone beam images. This step is based on Grangeat's [51] formula, which can be used in our proposed method to reconstruct the 3D cone beam images.

A.1. Helical Cone Beam Reconstruction and Preliminary
Results. To reconstruct the helically scanned objects, the scanned cone beam data should first be converted to fanbeam data and then to parallel beams. This rebinning process is based on the method introduced in [52], called cone beam single slice rebinning (CB-SSRB). CB-SSRB consists of the following two steps: (1) For each source position in the helical trajectory, , fix the -sampling distances.
(2) For each -slice, calculate the complete fan-beam set, from which the image can be estimated. This step uses the following equation to interpolate the cone beam scanned data on the fan beam points of interest: where ( , ) is estimated fan beam projection at source angle and axial position , ( , V) is the cone beam projections at helical position, is the distance between the source and the origin of the detector, and , V, and are geometry parameters defined in Figure 12. Each interpolated fan beam projection is weighted by where = ( /2 + )(1 − Δ / ), = 0.5 ( /2 + )/(2 ), is the pitch of the helical trajectory, and 2 is the maximum fan angle. The parallel beams ( , ) can then be estimated from the weighted fan beams ( , ) using (4), from which y will be calculated by computing the 1D Fourier transform of ( , )s. Figure 13 shows a simple simulated phantom reconstructed by the proposed helical reconstruction method. The helix source position is defined as = [ cos( ), sin( ), ( /2 )] and in this test pitch factor = 0.5. As can be seen, aside from the start or end of the scan, the reconstruction is almost perfect. However, when the image is close to one of the endpoints, the error of rebinning increases and as a result the image reconstruction error increases.
A.2. Future Work: Volumetric Cone-Beam CT. While in 2D geometries the X-ray projection and the Radon transform can be used alternately, in higher dimensions (e.g., cone beam geometry) the Radon transform and X-ray projections are very different. This makes it impossible to use inverse Radon transform and Fourier slice theorem to reconstruct 3D images from the X-ray projections, directly. The 3D Radon transform of a 3D function (⋅), shown in Figure 14, is defined as where is the unit vector that passes through the origin and the point of interest, which is described by ( , , ) in spherical coordinate and = ( , , ). This means that the Radon transform of an object at a point is equal to the integral of the object on a plane that passes the point and is normal to the vector that connects the origin to that point. The Grangeat formula relates line integral of cone-beam data to the Radon data in the whole Radon space for exact image reconstruction [51]. The link between the 3D radon and the X-ray projections can be expressed as follows: (shown in Figure 15), denotes the source to origin distance, and is the distance between the source and an arbitrary point along line . Other parameters are defined in Figure 15.
To have an exact reconstruction, the Grangeat formula requires the sources on a curve to contain at least one source at each plane meeting support of the object, (⋅). This condition is obviously not satisfied for a plane circle for which the FDK approximation has been derived. As shown in Figure 16, some data needed for exact image reconstruction are missed in this geometry. The missed data can be estimated using compressed sensing and the error adaptation weights proposed in this paper.
To reconstruct the cone beam images, the following process is proposed. Using (A.4), the radial derivatives of the 3D Radon transform of the object are estimated on 3D pseudopolar grids. The 3D pseudopolar grids consist of three groups of grids ≡ 1 ∪ 2 ∪ 3 : where , = /2, . . . , /2, = 3 /2, . . . , 3 /2, and is the number of pixels of the 3D function in each dimension. Similar to the 2D case, a fast 3D pseudopolar Fourier transform is available [53] enabling a 3D object (⋅) to be reconstructed from the Fourier transform of  the derivatives of the estimated 3D Radon transform. The interpolation errors caused by numerical implementation of the Grangeat equation are tracked and included in the EAW to minimize their effect in the recovered images and to enable the optimization algorithm to recover the missed data more effectively.