Accelerated Compressed Sensing Based CT Image Reconstruction

Hashemi, SayedMasoud; Beheshti, Soosan; Gill, Patrick R.; Paul, Narinder S.; Cobbold, Richard S. C.

doi:https://doi.org/10.1155/2015/161797

Computational and Mathematical Methods in Medicine

On this page

Abstract Introduction Results and Discussion Conclusion Appendix Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2015 | Article ID 161797 | https://doi.org/10.1155/2015/161797

Accelerated Compressed Sensing Based CT Image Reconstruction

SayedMasoud Hashemi,¹Soosan Beheshti,²Patrick R. Gill,³Narinder S. Paul,^1,4and Richard S. C. Cobbold¹

Academic Editor: Hugo Palmans

Received26 Nov 2014

Revised12 May 2015

Accepted20 May 2015

Published18 Jun 2015

Abstract

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.

1. 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–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 CS-based CT reconstructions are computationally expensive and the statistics of CT measurements are not usually incorporated in the problem formulation [11–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.

2. CS-Based CT Reconstruction and Its Challenges

Compressed sensing prescribes solving the optimization problem that can be represented byor other similar forms, for example, the norm of the image gradient, such asto 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, is the measurement matrix, is the column vector representation of the desired image (), is the measured data, is a sparsifying transform, , and denotes the total variation , 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 . Currently, in most available CS-based reconstruction methods, the measurement matrix 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, would be a 1080000 × 262144 matrix. Although this matrix is sparse, each iteration typically requires two multiplications by and , 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].

2.1. 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–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, samples are needed ( samples on 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 . 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)

(b)

A major objective of this paper is to accelerate the CS-based 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.

2.1.1. 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 , , , and are geometry parameters defined in Figure 2. Moreover, is the fan beam projected data and 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.

(a)

(b)

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

3. 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 , 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 ) and the photon counting noise (observed counts are drawn from a Poisson distribution of mean ). If we consider the discrete formulation in which denotes the vectorized , denotes the vectorized , and 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 is a function which depends upon measured data only and is a diagonal matrix. For the purpose of the MAP estimation may be ignored since it does not depend on . Ignoring this term, (6) describes a simplified CT model that can be written asin which is Gaussian distributed noise with a covariance matrix 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 . The th measured projection is given bywhere is noiseless and follows the Poisson distribution with . As a result the variance of projection data can be estimated from . Using as an unbiased estimation of , the diagonal elements of can be expressed asTo reconstruct the image from the projections, the MAP estimator can be used:Here 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.

3.1. 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, , can be modeled by generalized Gaussian distribution (GGD), that is, bywhere is the wavelet coefficients of the image , is inverse wavelet transform, and are the parameters of the GGD, and is the normalization parameter. It should be noted that when , the GGD is equivalent to Laplacian distribution and when , it describes a Gaussian distribution. Using (6), (10), and (11), the MAP model for CT images can be expressed aswhere , is a function of GGD parameters, and is typically in the range . 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 is a piecewise function spanned by the roof-top basis , the following class of probability distribution can be used to describe it:where , , is normalizing factor, and is a -valued random vector. When , this yields the total variation norm. Using 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, fromwhere and are regularization parameters.

4. 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 -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 isUsing (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 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 models the interpolation distance. Consequently, if the dose at each projection is high enough to ignore the electric noise , the EAW can be rewritten asUsing this definition, our proposed CS formulation to reconstruct the CT images can be expressed byThis 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, in this formulation acts similar to a Jacobi preconditioner and therefore could accelerate the convergence rate [40].

4.1. 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 . 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 . This process is illustrated in Figure 3.

In practice, 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.

5. Solving the Proposed CS Formulation

To solve the proposed formulation, a fast composite splitting algorithm (FCSA) [41–43] is used to decompose (19) into two simpler subproblems given byin which is the measured data interpolated/rebinned on the equally sloped lines, is the PPFT-based Radon transform, and is its adjoint. By calculating and , 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,in which 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 . To find , the optimization problem in step (2) of this algorithm can be solved by a wavelet soft thresholding algorithm [36]. Moreover, to calculate in step (3), the split Bregman TV based denoising algorithm as proposed in [45] was used. Finally, to estimate in the th iteration, (21) was used.

Initialize: , , , , , maxiter, tol
while > tol or < maxiter do
(1)
(2)
(3)
(4)
(5)
(6)
(7)
(8)
endwhile

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

6. 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 scanner (Toronto General Hospital, Canada). The scanner gathers data from 900 projection angles in each 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.

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

7.1. 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 adaptive-steepest-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 problem in ART step, a randomized Kaczmarz algorithm is used [4]. If is a linear system of equations and 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 , with probability proportional to , is the th row of , 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.

Initialize: , inner, outer, , tol
while > tol or < outer do
(1)
(2)
(3)
(4)
endwhile

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) 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 is 2 × 10⁻².

(a)

(b)

7.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). Figure 7 compares the reconstructed 512 × 512 images with the original image (1) using the inverse pseudopolar Fourier transform (least squares method), (2) 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.

(a)

(b)

(c)

(d)

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

(a)

(b)

(c)

(d)

(e)

(f)

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.

(a)

(b)

(c)

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

Appendix

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 fan-beam 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 , is the cone beam projections at helical position, is the distance between the source and the origin of the detector, and , , and are geometry parameters defined in Figure 12. Each interpolated fan beam projection is weighted by where , , is the pitch of the helical trajectory, and is the maximum fan angle. The parallel beams can then be estimated from the weighted fan beams using (4), from which 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 and in this test pitch factor . 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)

(b)

(a)

(b)

(c)

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 aswhere 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:where is the detector value at a distance of from the center of detector along the line perpendicular to (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.

(a)

(b)

(c)

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 :where , , 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.

Conflict of Interests

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

Acknowledgments

This study was partially supported by Toshiba Medical Systems Group, Canada, and NSERC Canada Grant no. 3247-2012.

References

D. J. Brenner and E. J. Hall, “Computed tomography—an increasing source of radiation exposure,” The New England Journal of Medicine, vol. 357, no. 22, pp. 2277–2284, 2007.
View at: Publisher Site | Google Scholar
M. S. Pearce, J. A. Salotti, M. P. Little et al., “Radiation exposure from CT scans in childhood and subsequent risk of leukaemia and brain tumours: a retrospective cohort study,” The Lancet, vol. 380, no. 9840, pp. 499–505, 2012.
View at: Publisher Site | Google Scholar
M. Jiang and G. Wang, “Convergence of the simultaneous algebraic reconstruction technique (SART),” IEEE Transactions on Image Processing, vol. 12, no. 8, pp. 957–961, 2003.
View at: Publisher Site | Google Scholar | MathSciNet
T. Strohmer and R. Vershynin, “A randomized Kaczmarz algorithm with exponential convergence,” Journal of Fourier Analysis and Applications, vol. 15, no. 2, pp. 262–278, 2009.
View at: Publisher Site | Google Scholar | MathSciNet
T. G. Herman and B. L. Meyer, “Algebraic reconstruction techniques can be made computationally efficient,” IEEE Transactions on Medical Imaging, vol. 12, no. 3, pp. 600–609, 1993.
View at: Publisher Site | Google Scholar
J. Fessler, “Statistical image reconstruction methods for transmission tomography,” in Handbook of Medical Imaging, Volume 2: Medical Image Processing and Analysis, M. Sonka and J. M. Fitzpatrick, Eds., SPIE, Bellingham, Wash, USA, 2000.
View at: Google Scholar
J. Browne and A. B. de Pierru, “A row-action alternative to the EM algorithm for maximizing likelihood in emission tomography,” IEEE Transactions on Medical Imaging, vol. 15, no. 5, pp. 687–699, 1996.
View at: Publisher Site | Google Scholar
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.
View at: Publisher Site | Google Scholar | MathSciNet
D. L. Donoho, “Compressed sensing,” IEEE Transactions on Information Theory, vol. 52, no. 4, pp. 1289–1306, 2006.
View at: Publisher Site | Google Scholar | MathSciNet
E. Y. Sidky, C.-M. Kao, and X. Pan, “Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT,” Journal of X-Ray Science and Technology, vol. 14, no. 2, pp. 119–139, 2006.
View at: Google Scholar
G.-H. Chen, J. Tang, and S. Leng, “Prior image constrained compressed sensing (PICCS): a method to accurately reconstruct dynamic CT images from highly undersampled projection data sets,” Medical Physics, vol. 35, no. 2, pp. 660–663, 2008.
View at: Publisher Site | Google Scholar
H. Lee, L. Xing, R. Davidi, R. Li, J. Qian, and R. Lee, “Improved compressed sensing-based cone-beam CT reconstruction using adaptive prior image constraints,” Physics in Medicine and Biology, vol. 57, no. 8, pp. 2287–2307, 2012.
View at: Publisher Site | Google Scholar
K. Choi, J. Wang, L. Zhu, T.-S. Suh, S. Boyd, and L. Xing, “Compressed sensing based cone-beam computed tomography reconstruction with a first-order method,” Medical Physics, vol. 37, no. 9, pp. 5113–5125, 2010.
View at: Publisher Site | Google Scholar
Z. Zhu, K. Wahid, P. Babyn, D. Cooper, I. Pratt, and Y. Carter, “Improved compressed sensing-based algorithm for sparse-view CT image reconstruction,” Computational and Mathematical Methods in Medicine, vol. 2013, Article ID 185750, 15 pages, 2013.
View at: Publisher Site | Google Scholar
J. C. Park, B. Song, J. S. Kim et al., “Fast compressed sensing-based CBCT reconstruction using Barzilai-Borwein formulation for application to on-line IGRT,” Medical Physics, vol. 39, no. 3, pp. 1207–1217, 2012.
View at: Publisher Site | Google Scholar
M. Chen, D. Mi, P. He, L. Deng, and B. Wei, “A CT reconstruction algorithm based on L_1/2 regularization,” Computational and Mathematical Methods in Medicine, vol. 2014, Article ID 862910, 8 pages, 2014.
View at: Publisher Site | Google Scholar
J. Chen, J. Cong, L. A. Vese, J. Villasenor, M. Yan, and Y. Zou, “A hybrid architecture for compressive sensing 3-D CT reconstruction,” IEEE Journal on Emerging and Selected Topics in Circuits and Systems, vol. 2, no. 3, pp. 616–625, 2012.
View at: Publisher Site | Google Scholar
D. Gottleib, B. Gustafsson, and P. Forssen, “On the direct fourier method for computer tomography,” IEEE Transactions on Medical Imaging, vol. 19, no. 3, pp. 223–232, 2000.
View at: Publisher Site | Google Scholar
T. Knopp, S. Kunis, and D. Potts, “A note on the iterative MRI reconstruction from nonuniform k-space data,” International Journal of Biomedical Imaging, vol. 2007, Article ID 24727, 9 pages, 2007.
View at: Publisher Site | Google Scholar
Y. Mao, B. P. Fahimian, S. J. Osher, and J. Miao, “Development and optimization of regularized tomographic reconstruction algorithms utilizing equally-sloped tomography,” IEEE Transactions on Image Processing, vol. 19, no. 5, pp. 1259–1268, 2010.
View at: Publisher Site | Google Scholar | MathSciNet
J. Miao, F. Förster, and O. Levi, “Equally sloped tomography with oversampling reconstruction,” Physical Review B, vol. 72, no. 5, Article ID 052103, 2005.
View at: Publisher Site | Google Scholar
E. Lee, B. P. Fahimian, C. V. Iancu et al., “Radiation dose reduction and image enhancement in biological imaging through equally-sloped tomography,” Journal of Structural Biology, vol. 164, no. 2, pp. 221–227, 2008.
View at: Publisher Site | Google Scholar
A. Averbuch, R. R. Coifman, D. L. Donoho, M. Elad, and M. Israeli, “Fast and accurate polar fourier transform,” Applied and Computational Harmonic Analysis, vol. 21, no. 2, pp. 145–167, 2006.
View at: Publisher Site | Google Scholar | MathSciNet
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.
View at: Publisher Site | Google Scholar
S. Hashemi, S. Beheshti, P. R. Gill, N. S. Paul, and R. S. C. Cobbold, “Fast fan/parallel beam CS-based low-dose CT reconstruction,” in Proceedings of the 38th IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP '13), pp. 1099–1103, IEEE, Vancouver, Canada, May 2013.
View at: Publisher Site | Google Scholar
A. Averbuch, I. Sedelnikov, and Y. Shkolnisky, “CT reconstruction from parallel and fan-beam projections by a 2-d discrete radon transform,” IEEE Transactions on Image Processing, vol. 21, no. 2, pp. 733–741, 2012.
View at: Publisher Site | Google Scholar | MathSciNet
J. Hsieh, Computed Tomography: Principles, Design, Artifacts, and Recent Advances, vol. PM188, SPIE Press Book, Bellingham, Wash, USA, 2nd edition, 2009.
Q. Xu, H. Y. Yu, X. Q. Mou, L. Zhang, J. Hsieh, and G. Wang, “Low-dose X-ray CT reconstruction via dictionary learning,” IEEE Transactions on Medical Imaging, vol. 31, no. 9, pp. 1682–1697, 2012.
View at: Publisher Site | Google Scholar
P. He, B. Wei, S. Wang, S. R. Stock, H. Yu, and G. Wang, “Piecewise-constant-model-based interior tomography applied to dentin tubules,” Computational and Mathematical Methods in Medicine, vol. 2013, Article ID 892451, 8 pages, 2013.
View at: Publisher Site | Google Scholar | MathSciNet
S. Ramani and J. A. Fessler, “A splitting-based iterative algorithm for accelerated statistical X-ray CT reconstruction,” IEEE Transactions on Medical Imaging, vol. 31, no. 3, pp. 677–688, 2012.
View at: Publisher Site | Google Scholar
J. Tang, B. E. Nett, and G. H. Chen, “Performance comparison between total variation (TV)-based compressed sensing and statistical iterative reconstruction algorithms,” Physics in Medicine and Biology, vol. 54, no. 19, pp. 5781–5804, 2009.
View at: Publisher Site | Google Scholar
E. J. Candès, M. B. Wakin, and S. P. Boyd, “Enhancing sparsity by reweighted $l_{1}$ minimization,” The Journal of Fourier Analysis and Applications, vol. 14, no. 5-6, pp. 877–905, 2008.
View at: Publisher Site | Google Scholar | MathSciNet
S. R. Deans, The Radon Transform and Some of Its Applications, John Wiley & Sons, New York, NY, USA, 1983.
View at: MathSciNet
J.-B. Thibault, K. D. Sauer, C. A. Bouman, and J. Hsieh, “A three-dimensional statistical approach to improved image quality for multislice helical CT,” Medical Physics, vol. 34, no. 11, pp. 4526–4544, 2007.
View at: Publisher Site | Google Scholar
S.-I. Maeda, W. Fukuda, A. Kanemura, and S. Ishii, “Maximum a posteriori X-ray computed tomography using graph cuts,” Journal of Physics: Conference Series, vol. 233, no. 1, Article ID 012023, 2010.
View at: Publisher Site | Google Scholar
M. Hashemi and S. Beheshti, “Adaptive Bayesian denoising for general Gaussian distributed signals,” IEEE Transactions on Signal Processing, vol. 62, no. 5, pp. 1147–1156, 2014.
View at: Publisher Site | Google Scholar | MathSciNet
M. Basseville, A. Benveniste, K. C. Chou, S. A. Golden, R. Nikoukhah, and A. S. Willsky, “Modeling and estimation of multiresolution stochastic processes,” IEEE Transactions on Information Theory, vol. 38, no. 2, pp. 766–784, 1992.
View at: Publisher Site | Google Scholar
M. Lassas and S. Siltanen, “Can one use total variation prior for edge-preserving Bayesian inversion?” Inverse Problems, vol. 20, no. 5, pp. 1537–1563, 2004.
View at: Publisher Site | Google Scholar | MathSciNet
J. Yang, Y. Zhang, and W. Yin, “A fast alternating direction method for TVL1-L2 signal reconstruction from partial Fourier data,” IEEE Journal on Selected Topics in Signal Processing, vol. 4, no. 2, pp. 288–297, 2010.
View at: Publisher Site | Google Scholar
H. Song, “Preconditioning techniques analysis for cg method,” ECS 231 Large-Scale Scientific Computation Course, 2013.
View at: Google Scholar
A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM Journal on Imaging Sciences, vol. 2, no. 1, pp. 183–202, 2009.
View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
J. Huang, S. Zhang, H. Li, and D. Metaxas, “Composite splitting algorithms for convex optimization,” Computer Vision and Image Understanding, vol. 115, no. 12, pp. 1610–1622, 2011.
View at: Publisher Site | Google Scholar
M. Jiang, J. Jin, F. Liu et al., “Sparsity-constrained SENSE reconstruction: an efficient implementation using a fast composite splitting algorithm,” Magnetic Resonance Imaging, vol. 31, no. 7, pp. 1218–1227, 2013.
View at: Publisher Site | Google Scholar
P.-E. Maingé, “Strong convergence of projected subgradient methods for nonsmooth and nonstrictly convex minimization,” Set-Valued Analysis, vol. 16, no. 7-8, pp. 899–912, 2008.
View at: Publisher Site | Google Scholar | MathSciNet
T. Goldstein and S. Osher, “The split Bregman method for $L1$ -regularized problems,” SIAM Journal on Imaging Sciences, vol. 2, no. 2, pp. 323–343, 2009.
View at: Publisher Site | Google Scholar | MathSciNet
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.
View at: Publisher Site | Google Scholar
L. Ritschl, F. Bergner, C. Fleischmann, and M. Kachelrieß, “Improved total variation-based CT image reconstruction applied to clinical data,” Physics in Medicine and Biology, vol. 56, no. 6, pp. 1545–1561, 2011.
View at: Publisher Site | Google Scholar
J. Bian, J. H. Siewerdsen, X. Han et al., “Evaluation of sparse-view reconstruction from flat-panel-detector cone-beam CT,” Physics in Medicine and Biology, vol. 55, no. 22, pp. 6575–6599, 2010.
View at: Publisher Site | Google Scholar
J. M. Bioucas-Dias and M. A. T. Figueiredo, “A new TwIST: two-step iterative shrinkage/thresholding algorithms for image restoration,” IEEE Transactions on Image Processing, vol. 16, no. 12, pp. 2992–3004, 2007.
View at: Publisher Site | Google Scholar | MathSciNet
K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, “Image denoising by sparse 3-D transform-domain collaborative filtering,” IEEE Transactions on Image Processing, vol. 16, no. 8, pp. 2080–2095, 2007.
View at: Publisher Site | Google Scholar | MathSciNet
S. W. Lee, G. Cho, and G. Wang, “Artifacts associated with implementation of the Grangeat formula,” Medical Physics, vol. 29, no. 12, pp. 2871–2880, 2002.
View at: Publisher Site | Google Scholar
F. Noo, M. Defrise, and R. Clackdoyle, “Single-slice rebinning method for helical cone-beam CT,” Physics in Medicine and Biology, vol. 44, no. 2, pp. 561–570, 1999.
View at: Publisher Site | Google Scholar
Y. Keller, Y. Shkolnisky, and A. Averbuch, “Volume registration using the 3-D pseudopolar Fourier transform,” IEEE Transactions on Signal Processing, vol. 54, no. 11, pp. 4323–4331, 2006.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2015 SayedMasoud Hashemi 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.

PDF Download Citation

Download other formats

Order printed copies

Views

2668

Downloads

1798

Citations