|
 
 |
| ORIGINAL ARTICLE |
|
| Year : 2022 | Volume
: 12
| Issue : 1 | Page : 8-24 |
|
Low-dose cone-beam computed tomography reconstruction through a fast three-dimensional compressed sensing method based on the three-dimensional pseudo-polar fourier transform
N Teyfouri1, Hossein Rabbani1, I Jabbari2
1 Medical Image and Signal Processing Research Center, School of Advanced Technologies in Medicine, Isfahan University of Medical Sciences, Isfahan, Iran
2 Department of Nuclear Engineering, Faculty of Advanced Sciences and Technologies, University of Isfahan, Isfahan, Iran
| Date of Submission | 30-Dec-2019 |
| Date of Decision | 24-Apr-2021 |
| Date of Acceptance | 20-Aug-2021 |
| Date of Web Publication | 28-Dec-2021 |
Correspondence Address:
Hossein Rabbani
Medical Image and Signal Processing Research Center, School of Advanced Technologies in Medicine, Isfahan University of Medical Sciences, Isfahan 8174673461
Iran
 Source of Support: None, Conflict of Interest: None
DOI: 10.4103/jmss.jmss_114_21
Background: Reconstruction of high quality two dimensional images from fan beam computed tomography (CT) with a limited number of projections is already feasible through Fourier based iterative reconstruction method. However, this article is focused on a more complicated reconstruction of three dimensional (3D) images in a sparse view cone beam computed tomography (CBCT) by utilizing Compressive Sensing (CS) based on 3D pseudo polar Fourier transform (PPFT). Method: In comparison with the prevalent Cartesian grid, PPFT re gridding is potent to remove rebinning and interpolation errors. Furthermore, using PPFT based radon transform as the measurement matrix, reduced the computational complexity. Results: In order to show the computational efficiency of the proposed method, we compare it with an algebraic reconstruction technique and a CS type algorithm. We observed convergence in <20 iterations in our algorithm while others would need at least 50 iterations for reconstructing a qualified phantom image. Furthermore, using a fast composite splitting algorithm solver in each iteration makes it a fast CBCT reconstruction algorithm. The algorithm will minimize a linear combination of three terms corresponding to a least square data fitting, Hessian (HS) Penalty and l1 norm wavelet regularization. We named it PP-based compressed sensing-HS-W. In the reconstruction range of 120 projections around the 360° rotation, the image quality is visually similar to reconstructed images by Feldkamp-Davis-Kress algorithm using 720 projections. This represents a high dose reduction. Conclusion: The main achievements of this work are to reduce the radiation dose without degrading the image quality. Its ability in removing the staircase effect, preserving edges and regions with smooth intensity transition, and producing high-resolution, low-noise reconstruction results in low-dose level are also shown.
Keywords: Three-dimensional compressed sensing, three-dimensional pseudo-polar Fourier transform, cone beam computed tomography reconstruction, hessian
How to cite this article:
Teyfouri N, Rabbani H, Jabbari I. Low-dose cone-beam computed tomography reconstruction through a fast three-dimensional compressed sensing method based on the three-dimensional pseudo-polar fourier transform. J Med Signals Sens 2022;12:8-24 |
How to cite this URL:
Teyfouri N, Rabbani H, Jabbari I. Low-dose cone-beam computed tomography reconstruction through a fast three-dimensional compressed sensing method based on the three-dimensional pseudo-polar fourier transform. J Med Signals Sens [serial online] 2022 [cited 2023 Mar 23];12:8-24. Available from: https://www.jmssjournal.net/text.asp?2022/12/1/8/334132 |
| Introduction | |  |
Today, X-ray computed tomography (CT) is widely used in hospitals and clinics with the purposes of diagnosis and intervention.[1],[2],[3],[4],[5],[6] There is no doubt about harmful side effects of X-ray radiation such as genetic problems, cancerous tissues, and other diseases. That justifies increased attention to radiation risk, yielding to famous slogan: As low as reasonably achievable principle is employed.[5],[6]
Considering quantum accumulation process in X-ray imaging, the signal to noise ratio (SNR) has quadratic dependency on the dose of X-ray. As other conditions are identical, any reduction in the X-ray dose will degrade the quality of image dramatically.[7],[8]
A hot topic in the field of CT is reconstructing adequate CT images at minimum radiation dose level. Two strategies are already developed to address this problem: Reducing the X-ray flux toward each detector element and decreasing the number of projections across the object to be reconstructed.[2] The former is often carried out by adjusting the operating current, the exposure time of an X-ray tube and the operating potential, which would result in noisy projections. The latter, which we utilize in this study, produces inadequate projection data, flawed with few-view, limited-angle, interior scan, and other deficiencies. Such deficiencies encourage algorithmic solutions in this area.
The first algorithmic strategy is to preprocess the projection data (with optimization variables) before image reconstruction. The second algorithmic strategy, which is also known as statistical iterative reconstruction (SIR), optimizes the maximum-likelihood or penalized-likelihood function formulated based on the statistical characteristics of projection data.[7],[8],[9],[10],[11] The image pixels or voxels would be accounted as optimization variables.
The compressive sensing (CS) theory,[12],[13],[14],[15],[16],[17],[18],[19] as a subcategory from SIR methods, has shot to prominence. It is instrumental in reconstructing a reliable and clean image from a limited number of noisy projection measurements.[9] The CS theory allows for accurate reconstruction of a sparse signal from samples much less than required by the Shannon/Nyquist sampling theorem. CS owes its success to the sparsity of a signal being studied. Although generally speaking an object is not sparse, a sparsifying transform may be used to convert it into a domain in which the signal has a sparse representation.[14],[17],[19],[20] This method would be able to reconstruct high-quality images from roughly one-tenth of the number of views required by filtered back projection (FBP) in two-dimensional (2D) images, thereby allowing a much lower dose scanning protocol than required in applying conventional reconstruction methods.[5] However, a major drawback with CS-based reconstruction algorithms is their being prohibitively computationally intensive for clinical use.[11],[21]
In this paper, Fourier-based reconstruction algorithms (FRAs), potentially enjoying high-speed implementation in the frequency domain[9],[11],[21],[22] are incorporated to reduce the computational complexity of image reconstruction. The Fourier reconstruction approaches are based on the relationship between the fast Fourier transform (FFT) of the image and FFT of the parallel-ray projections. The critical two steps are the estimations of the samples of the projection transform, on the central section through the origin of Fourier space, from the samples of the transform of the image, and vice versa for back-projection.
Interpolation errors are defined as limitation of Fourier-based methods of reconstruction. Our method, pseudo-polar based compressed sensing (PPCS), relies on as three-dimensional (3D) pseudo-polar Fourier transform (PPFT), which is a new form of FFT. In PPFT, the grid points in Fourier domain lie on equally-sloped lines rather than on equally angled lines. PPFT removes the errors associated with iterative re-gridding and the PP trajectory keeps the details untouched. PPFT has been proven to be a an algebraically precise and fast method in addition to being geometrically faithful and invertible.[22],[23],[24],[25],[26],[27],[28]
Our CS-based method has three steps: encoding, sensing, and decoding. In encoding, f is encoded into a smaller vector y = ϕf, including projections, by a linear transform ϕ. It's clear that y contains less information than f; therefore, it is a compression of f. In many CS applications, the linear transform ϕ is not calculated by a computer; rather it has been obtained by certain physical or digital means. In cone beam computed tomography (CBCT), for example, ϕ represents a PPFT, where it is “partial” as ϕf would only yield the Fourier coefficients corresponding to an incomplete set of frequencies. It has to be noted that because f is unknown during this step, ϕ will only be chosen independently of f. In sensing, y is acquired by X-ray projections on the objects that are measured on the detectors before being sent to a computer. Decoding requires recovering f from y. As f is sparse, it may be found as the sparsest solution to the underdetermined equations y = ϕf unless another even sparser solution is found to these equations.
The main advantage with this article in comparison with the published compressed sensing techniques is the acceleration of the CBCT reconstruction by lowering the CS complexity. To that effect, 3D PPFT-based Radon transform as suggested in[28] and fast composite splitting algorithm (FCSA) are used for solution. Furthermore, the simultaneous use of Hessian (HS) penalty and 3D wavelet transform in the reconstruction of CBCT images will remove blocky and staircase artifacts while the edges are retained. We refer to our method as PPCS-HS-W.
The rest of the article is organized as follows: Section II outlines different stages of the method including step-by-step implementation and verification of the optimization problem. Section III expresses the criteria and data used to evaluate the algorithm. Section IV is dedicated to experiments. We provide the results obtained from computer simulations to illustrate the performance of the method. Discussions and conclusion are provided in Sections V.
| Methods | | |
In this section, we elaborate a fast method developed to solve a l1-HS constrained regularization problem for a 3D image reconstruction of conventional geometry of single circular source trajectory CBCT with a limited number of projections.
Problem formulation
In the case of normal clinical exposures, assuming the source to be monochromatic, the X-ray CT measurements are modeled mainly as the sum of a Poisson distribution representing photon counting data and an independent Gaussian distribution denoting additive electronic noise,[29]
where  is a matrix describing the intersections of x-rays and voxels in the image f, iT and jT are the total number of projection measurements and image voxels, respectively.  represents the vector of attenuation coefficients of the object to be reconstructed and yi=[ϕf]i represents the line integrals through the image for a variety of positions and projection angles,[30] where  denotes the vector of projection measurements. The CT transmission scan does not provide the projection measurements directly, but rather is formed of a collection of recorded detector measurements Ii that are related to the line integral projections by Beer's law of attenuation. They represent the detected x-ray intensity after attenuation by the scanned object, and follow a Poisson distribution. mi and σi represents the mean and standard deviation of electronic noise that has been converted to photon units, The offset mean of m background signals such as dark current would be estimated using blank measurements before each scan and subtracted from the measured intensity. Therefore, we assume m=0 hereafter. IT denotes the impinging X-ray photon intensity emitted from the X-ray.[29]
Modeling the electronic noise is significant in low-dose CBCT. While there is no simple analytical form for the likelihood function of the combined Poisson and Gaussian model in (1), the sum, after adding a constant, can be approximated by a Poisson distribution.[31]
The reconstruction problem may be formulated in the Bayesian framework as the maximum a posteriori (MAP) estimate.
where P(.) denotes the probability, which is equivalent to:
Using the second order Taylor series expansion of the Poisson distribution,[30] the log likelihood of the measurements is given by:
where  is a diagonal weighting matrix with the ith diagonal element,. di (y3) may be ignored since it does not depend on f. Ignoring this term, (5) describes a shortened CT model that can be written as y = ϕf + w, where w denotes an additive Gaussian noise with a covariance matrix D-1. di is proportional to the detector counts, corresponding to the maximum likelihood of the inverse of the variance of the projection data, that is, to  .[32],[33],[34]
Where IT is the number of radiated photons from the X-ray source. Ii follows the Poisson distribution with  , where Īi is the average detector count.
Based on (6) we have:
If h(X) is a function of X, then Taylor expansions for the moments functions of random variables establish that:[35]
Where is E [X] the expected value of X. h(X) = log(X) So if then
Since electronic noise is an additive process with standard deviation σn, the variance in the total detector counts is  . By means of an unbiased estimation of Īi, Ii,[9] di can be defined as:[36]
The tomography operator, ϕ in y = ϕf+w, is ill-posed and therefore cannot be inverted. We apply variational reconstruction methods wherein a solution is found through convex optimization problem.[37]
where J(f) is a prior energy, whose choice will be explained in the following sections.
The main challenge in solving (11) within a reasonable amount of time arises from the size of the measurement matrix ϕ. Currently, in most available CS-based reconstruction methods, projection operator or 3D Radon transform, ϕ consist of integrals over planes through image data points and is a (NI×NJ×NK) large scale matrix with entries ϕijk
where  is a vector with length of {ρi}(0≤i≤I) in the direction of {θj=2π/J}(0≤j≤J) and {φk=2π/K}(0≤k≤K) in spherical coordinates [Figure 1], that describes the length and direction of the normal from origin to the integration plane. | Figure 1: Sphere coordinate. (a) Coordinate of a point on three-dimensional radon space. (b) Sinogram Sampling on a conventional three-dimensional Radon diameter
Click here to view |
To solve (11) iteratively, each iteration would typically need two multiplications by ɸ and ɸT,[38] which would result in a very significant increase in the computation burden for image reconstruction, as compared to Feldkamp-Davis-Kress (FDK)-based methods. To do the CS-based CT reconstruction within a reasonable computation time, graphics processing units (GPU)-based algorithms are suggested.[39],[40]
Reduction of complexity and elimination of interpolation errors using the pseudo-polar Fourier transform
In order to alleviate the computational burden on the Radon transform, we use the Fourier slice theorem (FST), which builds up a 3D FFT of the object through 1D FFT of the projections.[9],[11],[22],[26],[41],[42]
As depicted in [Figure 2]a the 3D conventional FST[28] shows that (Ŷφk,θj)(ω), the 1D radial FFT of (ɸφ,θ f)(ρ) along a line with direction of φk,θj, is equivalent to a line of the 3D FFT of f in the same direction,  φk,θj,[37] where ω is frequency variable. | Figure 2: Comparing conventional Fourier slice theorem with Fourier slice theorem on pseudo-polar grid. (a) Three-dimensional conventional Fourier slice theorem on spherical coordinate, sinogram sampling in the Radon space. (b) Three-dimensional Fourier slice theorem on pseudo-polar grid, linogram sampling in the Radon space[28]
Click here to view |
According to FST, we explore a FFT version of ɸ in (11), since it corresponds to the sampled FFT of the image along discretized rays.
The major challenge with conventional FST is the provision of the 3D Radon and Fourier space on spherical coordinates as sinogram sampling [Figure 2]a. In the calculation of 3D inverse FFT, the data needs to be on Cartesian coordinate. Meanwhile, accurate 3D interpolation is needed in the frequency domain to exchange between spherical and Cartesian lattice. 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.
Rebinning is transforming the diverging cone beam projections to parallel beam projections. Using the 2D and 3D FST and debt recovery tribunals[23] in the stage of the rebinning, instead of sinogram sampling in Radon space, an accurate and fast linogram sampling has been proposed by Teyfour et al.,[28] This algorithm applies 3D PPFT and its inverse on Radon transform (T and ɸT) and in pseudo-polar (PP) grid with complexity of O(N3logN) [Figure 2]b. PP grid  shown in [Figure 2]b would be of help in bypassing the interpolation stage.[28],[43],[44] The PPFT-based Radon space is ideally acquired at angles corresponding to the lines of the PP grid, which consists of concentric squares with a horizontal and a vertical groups of lines [Figure 3]. | Figure 3: Different grids. (a) Two-dimensional cartesian grid. (b) Two-dimensional PP grid. (c) Three-dimensional pseudo-polar grid
Click here to view |
In reconstructing under-sampled radial data with CS, we would need regridding and inverse-regridding to transfer data between the frequency and image domains. In each CS iteration, 3D interpolations are implemented twice in the regridding and inverse-regridding. That would give rise to errors and affect the reconstruction quality. To resolve these problems, a radial-like PP trajectory for the CS CBCT applications would allow for an exact image reconstruction with PPFT rather than interpolations.
PPFT has three significant properties making it an ideal substitute to conventional Fourier methods: (1) PPFT provides an exact PP-based FFT and its inverse, setting the gridding error to zero, (2) PPFT has a fast-forward and a fast-backward calculation algorithm,[23],[28],[45] enabling our proposed algorithm to avoid the regridding step used in iterative non-Cartesian Fourier-based reconstruction methods, and (3) PPFT has an analytical inverse function.[45]
Designing the optimization problem
It has been shown that each projection is capable of reconstructing a limited number of Radon space diameters.[28] In this study, due to a limited number of projections and using of PPFT version of ɸ, there are only a few nonzero diameters in 3D PPFT Radon space.
Hence, recovering from tomography measurements, y = ϕf is equivalent to inpainting the missing Fourier frequencies. We assume the partial noisy Fourier measures as
Where ŷ is 1D radial FFT of the 3D Radon space,  denote 3D PPFT of 3D unknown image and Ŵ[ω] is the measured noise on projections in the frequency domain.
Iterative reconstruction algorithms have demonstrated excellent performance in improving the quality of the image. We therefore the use an improved algorithm by manipulating the penalty term in objective function to increase the quality of the reconstructed images. For this purpose, HS penalty term (‖f‖HS) and wavelet transform (‖Ψf‖1) are proposed to be added to the optimization problem. Justifications for such choices are elaborated in section II. C. After combining all terms, we will reach (15) instead of (11).
Solving the optimization problem
FCSA[14] is selected to solve (15) by solving a composite regularization problem including HS and l1 norm term, effectively. FCSA is motivated by FISTA[46] which minimizes the following relation:
Where P denotes a smooth convex function with Lipschitz constant LP, and Q denotes a convex function that can be nonsmooth. The efficiency of FISTA largely depends on its ability to quickly solve its second step fk = proxρ(Q)(fQ) With a continuous convex function Q(x) and any scalar P>0, the proximal map associated with function Q is defined as:
FISTA easily solves the l1 regularization problem as the second step fk= proxρ(β‖Ψf‖1)(fQ) has a close form solution.
The HS regularization problem, fk= proxρ(α‖Ψf‖1)(fQ) is solved by the fast alternating minimization (FAM) algorithm.[4],[47],[48],[49]
However, FISTA does not solve the composite l1 and HS regularization in (15) efficiently since there would be no efficient algorithm to solve (18).
FCSA,[14] solves the problem (18) by combining the composite splitting denoising, FISTA and FAM.
In this section, we describe two subproblems related (18).
Three-dimensional hessian penalty
The total variation denoising problem is as follows:
Where u is an observed noisy image. This is particularly suitable for medical imaging from body organs with relatively constant gray value, thereby resembling the cartoon image model. That indicates the ability of total variation in recovering sharp features while inpainting Fourier measures.
However, the TV penalty calculated from first-order derivatives results in the unwanted staircase artifact that would often make reconstructed images over-sharpened and unnatural[4],[48],[49],[50] In this article, we suggest a second-order derivative penalty involving the Frobenius norm of the HS matrix from an image for CBCT reconstruction.[49] When using the second order penalty, some of the most favorable properties of the TV penalty, including homogeneity, convexity as well as rotation and translation invariance, are retained. The performance of the second order penalty is better in preserving gradual transition structures in the reconstructed images.[47] HS penalties can reflect the smooth intensity transitions of the underlying image without causing any staircase effect.
For 3D signals f: Ω C, we know the standard isotropic TV regularization penalty is the l1 norm of the gradient magnitude,[49] specified as:
Isotropic Higher degree TV (HDTV) regularization is specified by (21).
where fu,n is the nth degree directional derivative defined as:
Given a specific nth degree differential operator D=∑|α|=2cα ∂α, (α is a multi-index and the cα are constants), we define the generalized HDTV penalty as:
where  and Du is discrete derivative operators for u∈S that penalizing all rotations in 3D.
One choice of D is so that the expression of G[D p](f) will be as (24).[49]
That H is HS matrix and ‖Hf‖Frob is Frobenius norm of the HS in 3D, equivalent to the isotropic generalized HDTV penalty.[49] Thus, instead of (19), the image recovery problem with the HS penalty in this work is as (25).
To solve (25), we used the effective algorithm of FAM whose speed is ten times higher than the iteratively reweighted majorize minimize algorithm.[49] In order to facilitate understanding of the FAM optimization parameters given in the section IV, the details of this algorithm for special case of denoising is given in the appendix.
Three-dimensional sparsifying transform
To provide further improvement in the results, 3D Daubechies wavelet transform, Ψ is proposed in (26) to serve as a sparsifying function.
We use a cycle spinning hard threshold denoising by taking average from the denoising results of the translated version of the signal which is expected to lighten blocking artifacts due to thresholding in orthogonal wavelets.[51] The threshold was adjusted to T, since we achieved the best result by trial and error.
Mask generation
In the rebinning step (redistribution of the cone beam projections to equally sloped lines in a PP grid), each characteristic point of the 3D Radon space [Figure 2]b, C:(ρ,θ,φ) is mapped to a point D:(s,α,ψ) on the detectors based on Grangeat's formula [Figure 4].[28] | Figure 4: (a) Three-dimensional Radon space on sphere coordinate. (b) Detector on polar coordinate.[28] The OC line is perpendicular to the integration plane and Xψ is the virtual detector at the angle Ψ of the projection
Click here to view |
Assuming limited number of projection angles, ψt, an interpolation step is needed to calculate all of C:(ρ,θ,φ) using Dt:(st,αt,ψt), where (st,αt) show the polar coordinate of pixels on each detector.
The larger the angular distance between Ψ and the nearest ψt would be, the greater the error of interpolating will be. To compensate interpolation errors and obtain a high quality image, we consider matrix ξ as the mask operator. ξ is as large as ϕ. In this matrix, the points located far from original projection angles (ψt) are set to zero. So all points on some diameters of Radon space with the same angular direction of will be zero. Their values will be computed by the CS algorithm. The mask operator is then constructed as follows:
Where ψ1 and ψ2 are the two nearest ψs to ψt, and empirically, we found that considering M=(ψ1-ψ2)/4 give the best results. Note that ξ is computed only once before the iterative algorithm.
In short, let f be the 3D image to be reconstructed in the object domain and  be the 3D PPFT of f. At first, we compute Radon transform on PP grid in space domain (the rebinning process) by Teyfouri's method.[28] After applying 1D radial FFT on each diameter (ŷ), what we obtain from tomography data is:
Where ξ denotes a linear operator or mask matrix selecting the entry of approachable frequency values and ŷ denotes an array in the Fourier domain containing the known frequency values. According to “frequency constraint,”  is restricted in the subspace  , and only changeable elements are unknown (equal to zero) elements. We solve (15) inside the space  .
It must be also pointed out that the rebinning process has to be applied only once, before the initiation of the iterative process. Since, the parallel projections are not calculated at all the lines of the PP grid, rather only limited portion of the PP grid in Fourier space are filled with ŷ, the rest of the lines are filled in by the iterative CS algorithm.
Time complexity and convergence analysis
Applying CS-based CBCT reconstruction algorithms to medical imaging has proven difficult as it requires a prohibitive time to compute. A single, complete iteration requires at least one forward and one backward projection calculations that are computationally expensive. Thus, the efficiency of an algorithm requires (1) a minimal number of forward and backward projection calculations per iteration, in addition to the necessity of (2) converging in a minimal number of total iterations.
We now discuss the computational complexity of Algorithm 1. In each iteration, there are three operators with the highest computational cost: the PPFT  , its inverse operator,  T and the FCSA solver. The computational complexity of  [23],[28] and  T[28],[45] has been proven to be O(N3logN).
To generate initial image for Algorithm 1, use of PPFT based method makes it possible to reduce the arithmetic operations from O(N4) with the conventional FDK method[44] to O(N3log(N))[23],[28] for an image with N × N × N voxels. This result in a faster convergence compared with setting f1=FDK.
FCSA solver in the PPCS-HS-W method includes FAM and FISTA solvers. We used an efficient FAM solver has been proposed In[49] which is faster than TV regularization problem with cost O(N3).
The cost of FISTA solver is O(N3 log(N3)) in each iteration. The step 4, 6 and 7 of “for” loop of algorithm only involve adding vectors or scalars, thus cost only O(N3). In the step 1,  since  in this case. Thus, this step only costs O(N3 log(N3)). As introduced above, both of two steps fk= proxρ(2α‖f‖Hess)(fg) and fk = proxρ(2β‖Ψf‖1)(fg) has a close form solution and can be computed with cost O(N3 log(N3)) in each iteration.
Another important specification is its fast convergence performance borrowed from FISTA. FISTA is able to obtain an ∈-optimal solution in  iterations,[14] implying that the required number of the iterations for reaching  .
In[14] it has been proved that if {fk} is the sequence generated by the  , then, fk will converge to  .
Summary of pseudo-polar based compressed sensing-HS-W algorithm
After combining the entire process together, Algorithm 1 will elaborate the proposed method.
Algorithm 1: Summary of pseudo-polar based compressed sensing-HS-W algorithm for solving[15]
- Input:
 , Zero matrix and CBCT parameters. - R = Compute 3D Radon space from cone beam projections on PP grid by Teyfouri's method[28].
- ŷ =Apply radial 1D FFT on diameters of R.
- ξ = Compute mask matrix.
- ŷ = ξŷ.
- For k = 1 to K do{
radical FFT of Randon transform.
Algorithm 1 is executed easily due to its concise structure. However, we need to address some details of the process of implementation. Two things need to be considered. We want first to have good performance and fast convergence. Then, we intend to minimize the number of parameters needed to tune up in practice. Some parameters are available for tuning the optimization problem of (15). The guidelines for parameters selection are summarized in [Table 1].
| Materials and Evaluation Indexes | | |
The experiments were carried out on three computer simulation phantoms: A compressed sensing (CS) phantom,[52] a Head phantom and a modified Shepp-Logan phantom. [Figure 5]a demonstrates a representative slice of a noise-free image out of the CS phantom with size of 64 × 64 × 64. In the phase of initial Radon generation, due to the memory constraint of Matlab, in none of the experiments we were able to consider the image size larger than 64 × 64 × 64 voxels. | Figure 5: One representative slice of CS phantom, (b) horizontal profile of center of octahedron at the shown slice
Click here to view |
The phantom with a uniform background contains an octahedron with a linearly gradual transitioning intensity. [Figure 5]b plots the horizontal profile of the center of octahedron along the shown slice. Furthermore, the phantom contains a set of line objects to measure the resolution, a set of circular cylinders of different diameters and different intensities to measure the contrast-to-noise ratio (CNR), as well as a ball with linear gradual transitioning intensity resembling the octahedron.
To demonstrate the performance of PPCS-HS-W method in the low dose configurations, the number of projections in our framework was selected to be 720 in high, 120 in medium and 360° in low dose level in each rotation, and zero mean Gaussian noises was added to the projections. The incident X-ray intensity, I0, and the background electronic noise variance, σ, are set to 0.6 × 105 and 10, respectively.
The parameter values, customized for different phantoms used in the experiments, are summarized in [Table 2].
To measure the quality of reconstructed images using a variety of penalty terms, we measured peak signal to noise ratio (PSNR), improvement signal to noise ratio (ISNR), CNR,[53] structural similarity index (SSIM),[47],[54],[55] and full width at half maximum (FWHM).[56]
The SSIM metric is defined as:[47]
Where a and b are two local windows of size 8 × 8 pixels in two images. The two windows have the same position, and μa and σa and μb and σb are their mean and standard deviation, respectively. σa,b is the covariance between the two windows. C1 and C1 are two constants to avoid instability. In this study, C1 and C1 were chosen as C2=(0.01μmax)2 and C1=(0.01μmax)2. While the two windows move pixel-by-pixel over the reference image and the reconstructed image, a SSIM map will be obtained.[47]
| Results | | |
In this section, the performance of the proposed method is evaluated in terms of the accuracy (section), convergency and reconstruction time (section IV.) and reduction of radiation dose (section).
To evaluate the performance of PPCS-HS-W, we have compared it with some other algorithms. First ordered-subset simultaneous algebraic reconstruction techniques (OS-SART)[57] was implemented. second the adaptive steepest descent projections onto convex set (ASD-POCS)[58],[59] was implemented. Third and fourth, PPCS-TV-W and PPCS-HS method were presented to demonstrate the effect of Hessian and l1 norm regularization terms in (15). In addition, the inverse of Radon transform (iRadon)[28] as an initial images of our iterative method and FDK reconstruction as a conventional method are compared.
A SART-type family is a set of algorithms derived originally from the Kaczmarz method. It is also adapted to work projection by projection rather than row by row. This group of algorithms follows the equation:
Where V and W denote weight matrices based on ray length. OS-SART updates the image with a subset of the projections.[57]
ASD-POCS is an algorithm where the total variation norm requires the image to be piecewise smooth. It is expressed as:
Where ϕ is a matrix in space domain or 3D Radon space describing the intersections of X-rays and voxels in the image.
In PPCS-TV-W we have:
Where,
We solved (32) by FCSA including two FISTA solvers.[46]
In PPCS-HS-W without l1 norm named PPCS-HS, we have:
All notations are similar to (15). We solved (34) by FAM.[49]
In each case, the regularization parameters are optimized to obtain the optimized SNR that would ensure fair comparisons between different schemes. In our method, to achieve good quality reconstructions, the regularization parameters α and β are set as 0.001 and 0.035. In all of cases presented in this work, we stopped the computation at 50 iterations at which good convergence was seen through the observation of the reconstructed image at each iteration. In our experiment, the step size is set as 1. We set the max inner and outer iteration time of Hessian regularization problem as 7 and 10, respectively, and choose 20 and 2 for continuation parameter λ initialization and λfac Error tolerance is as 10-7 inner loop convergence tolerance. T, the threshold of wavelet coefficient is 0.18.
Accuracy of the reconstruction
The quality performance of PPCS-HS-W was assessed in subsections to A.7.
Visual assessment
We chose the red rectangles in [Figure 5]a as region of interest (ROI) to compare the quality of different methods. The enlarged rhombus ROI of reconstruction results acquired through iRadon, PPCS-TV-W, PPCS-HS, and PPCS-HS-W method with 36, 120, and 720 projections as low, medium and high dose levels are shown in [Figure 6]a,[Figure 6]b,[Figure 6]c,[Figure 6]d,[Figure 6]e,[Figure 6]f,[Figure 6]g,[Figure 6]h,[Figure 6]i,[Figure 6]j,[Figure 6]k,[Figure 6]l,[Figure 6]m,[Figure 6]n,[Figure 6]o respectively. | Figure 6: Visual assessment of region of interest (rhombus) in reconstructed images with different algorithms and different projection numbers. (a) iRadon. (b) Adaptive steepest descent projections onto convex set. (c) Pseudo-polar based compressed sensing-TV-W. (d) Pseudo-polar based compressed sensing-HS. (e) Pseudo-polar based compressed sensing HS-W. (f) iRadon. (g) Adaptive steepest descent projections onto convex set. (h) Pseudo-polar based compressed sensing-TV-W. (i) Pseudo-polar based compressed sensing-HS. (j) Pseudo-polar based compressed sensing-HS-W. (k) iRadon. (l) Adaptive steepest descent projections onto convex set. (m) Pseudo-polar based compressed sensing-TV-W. (n) Pseudo-polar based compressed sensing-HS. (o) Pseudo-polar based compressed sensing-HS-W
Click here to view |
As the reconstruction results indicate, all IR methods worked well in high dose level, while the iRadon algorithm could not reduce noise and caused significant artifacts even in high dose level. The TV penalty in both 2nd and 3rd column of [Figure 6] caused a serious staircase effect and gave rise to several artificial piecewise constant regions in low dose level. The reconstructed images with PPCS-HS look more natural and smoother, emphasizing on ability of PPCS-HS in preserving smooth regions as well as suppressing the noise.
Compared to PPCS-HS, using a sparsity regularization in an orthogonal basis was made to better reconstruct edges in the PPCS-HS-W method [red arrow in [Figure 6]e], while denoizing artifacts are reduced by using a translation invariant wavelet transform [orange arrow in [Figure 6]d].
Quantitative comparison
We first studied the performance of PPCS-HS-W with sparse views and the fixed incident X-ray intensity I0 of 0.6 × 105 photon counts per ray. To assess quality of the reconstructed image with different methods, the mean PSNR, ISNR, and SSIM of the whole slice obtained in ten repetitions for CS phantom are listed in [Table 3]. Quantitatively, PPCS-HS-W outperformed others. PPCS-HS-W would produce less noise and more structures than two others. | Table 3: Peak signal to noise ratio, improvement signal to noise ratio, and structural similarity Index values of images with different reconstruction algorithms for compressed sensing phantom
Click here to view |
Ability to preserve local structure
[Figure 7] illustrates a representative slice of the results acquired through applying different reconstruction methods for the head phantom from 120 projections. The ability of PPCS-HS-W in preserving local structure is clear. | Figure 7: A representative slice reconstructed in low dose condition. (a) Original head phantom. (b) iRadon. (c) Adaptive steepest descent projections onto convex set. (d) Pseudo-polar based compressed sensing-TV-W. (e) Pseudo-polar based compressed sensing-HS. (f) Pseudo-polar based compressed sensing-HS-W
Click here to view |
We used SSIM (the closer to 1, the higher structural similarity) in our study to measure the degree of similarity in local structures between reconstructed images and the original image. The SSIM map of the anthropomorphic head phantom [Figure 8] in PPCS-HS-W was whiter, indicating a higher ability in preserving structures. | Figure 8: Structural similarity index map of the head phantom using different methods. (a) Pseudo-polar based compressed sensing -TV-W (structural similarity index = 0.65). (b) Pseudo-polar based compressed sensing -HS-W (structural similarity index = 0.93)
Click here to view |
Comparison of edge preservation ability
In order to provide a quantitative analysis on performance of penalty terms in edge-preservation, we focused our study on the FWHM of reconstructed images along the red line in [Figure 9]. | Figure 9: Representative slices of original image in order to explore the ability to preserve the edges along the red line
Click here to view |
The penalties showed clear differences in sharpness along the red line. The FWHM value for PPCS-TV-W and PPCS-HS-W were 2.35 and 1.2, respectively. However, PPCS-HS-W conduced to the smallest FWHM value; it means that reconstructed image using HS penalty has the sharpest edges, while TV penalties produced reconstructed images with less sharp edges.
Furthermore, we selected the strip area [Barcode indicate by the pink square in [Figure 5]a. Line profile of the reconstructed barcode in [Figure 10] indicates that the proposed method preserved the edges better than iRadon, ASD-POCS and PPCS-TV-W. | Figure 10: Line profile of the pink region of interest in Figure 5 to explore the ability of our proposed method to preserve the edges
Click here to view |
Ability to provide high spatial resolution and deblurring
[Figure 11] shows a representative slice of the reconstructed Shepp Logan phantom images generated using different methods of iRadon [Figure 11]b, Pseudo-polar based compressed sensing-TV-W [Figure 11]c and Pseudo-polar based compressed sensing-HS-W [Figure 11]d. The noise level was calculated by using the red rectangle area. Two ROIs are selected to make a quantitative comparison of CNR in different reconstruction methods [shown by arrows in [Figure 11]a. | Figure 11: A representative slice of different reconstructed algorithms with low dose view. (a) Original phantom. (b) iRadon. (c) Pseudo-polar based compressed sensing-TV-W. (d) Pseudo-polar based compressed sensing-HS-W
Click here to view |
[Table 4] lists CNRs corresponding to ROIs in [Figure 11]. It is clear that by changing projection number, the noise level shows variation. In comparison, HS was better than the TV penalty for denoising with different projection numbers. | Table 4: Contrast-to-noise ratio for different regions of interest of Shepp - Logan phantom in [Figure 8]
Click here to view |
Since CNR is dependent on background noise level, this value should not be ignored in the evaluation of the quality of reconstructed images.
To underscore the significant contribution of PPCS-HS-W in deblurring task, we calculate FWHM in blue line indicated in [Figure 9] FWHM in the TV penalty is 2.69 and that of HS stands at 1.75, showing the higher image resolution in PPCS-HS-W. These numerical observations are consistent with visual inspection in [Figure 12], in which the iRadon caused significant artifacts [Figure 12]b, the TV penalty slightly blurred the image edges [Figure 12]c and the HS penalty preserved better edges [Figure 12]d. | Figure 12: A representative slice of image reconstructed in low dose view. (a) Original CS phantom. (b) iRadon. (c) Pseudo-polar based compressed sensing-TV-W. (d) Pseudo-polar based compressed sensing-HS-Wv
Click here to view |
Ability to remove blocky artifact
[Figure 11]c and [Figure 11]d shows the ability of PPCS-HS-W method compared to PPCS-TV-W in removing blocky artifacts. The image is reconstructed from 36 projections around 360° and the power of the proposed method in removing blocky artifacts is quite clear through the visual comparison.
Ability to remove staircase effect and preserve smooth images
The intensity curves along the red line in [Figure 12]a show that TV-based reconstruction produced several staircase artifacts, while the reconstructed results using HS penalties were far smoother [Figure 13].
The regions of smoothly-changed intensity, such as octahedron in the upper-right corner and the sphere in the lower-left corner of [Figure 12], exhibited numerous small and unnatural constant intensity artifacts in the reconstructed image, where the TV penalty was used. It stemmed from the staircase effect [Figure 12]c. The results yielded by HS penalty were more natural [Figure 12]d, indicating its capacity in suppressing the staircase effect.
Convergence analysis and central processing unit time
We implement our method on central processing unit (CPU) and apply it on a 3D Shepp-Logan phantom of size 64 × 64 × 64 voxels with 36 projections around 360° to compare its speed with other algorithms. Our MATLAB codes ran on a PC with one i7–6800K CPU@ 3.40GHz and a 32 GB RAM.
A rebinning step was performed before initiating the iterative algorithm to transform the cone beam projections to parallel projections along equally sloped lines of the PP gird.[28]
To prove the ability of the frequency-domain methods in increasing the speed of iterative regularization problem solving, we made a comparison between PPCS-TV-W and PPCS-HS-W and two algorithms in the spatial domain, OS-SART[57] and ASD-POCS.[58],[59] [Figure 14] shows the reconstructed 3D images using these algorithms in three different iteration numbers. The level of quality agreement goes in the order of PPCS-HS-W >PPCS-TV-W >ASD-POCS >OS-SART. This finding proves to be correct at all levels of iterations, as shown in [Figure 14]. | Figure 14: The reconstructed images of the Shepp-Logan phantom, using different algorithms, as a function of 10, 30, and 50 iterations. A total of 36 projections in cone-beam geometry were used for the reconstructions
Click here to view |
As [Figure 15] indicates, the PPCS-TV-W and PPCS-HS-W algorithms show convergence at about 20 and 15 iterations, respectively. That is while the OS-SART and ASD-POCS algorithms need further convergence at 50 iterations. | Figure 15: Mean-squared relative error as a function of the number of iterations, for the respective four algorithms
Click here to view |
The relative error is the mean squared percent error from the original image pixel values.
Where fi,j,k denotes the voxel values in the reconstructed volume f and  represents the original values of the Shepp-Logan phantom used.
A fair comparison study on the recovery time vs CPU time (s) would require the MATLAB code of the other CBCT reconstruction algorithms with CPU implementation. Since the codes of OS-SART and ASD-POCS and some available others are in GPU,[59],[60],[61],[62] we didn't report any comparative study on the speed of our algorithms vs time. 2D version of PPCS-TV-W, solved by FSCA in,[9] by Bregman iterative regularization in[11] and by nonlocal TV regularization in[21] has been shown that PPCS-TV-W is faster than FBP and algebraic reconstruction technique-TV based method for fan beam images reconstruction. In addition, Huang[14] used Fourier based CS-TV-W in 2D magnetic resonance imaging (MRI) images reconstruction, solved by FSCA too, and showed that it is faster than the other Fourier based CS algorithms, such as TVCMRI, RecPF, Sparse MRI.
The average computation times of the PPCS-TV-W and PPCS-HS-W method were 1336 and 1535 seconds with 50 iterations in the given studies with the phantom of size 64 × 64 × 64. For image applications involving real-time reconstruction, the computational performance can be enhanced by implementing the reconstruction code in more effective formats, such as C language, rather than the commonly used MATLAB scripts.
Ability to recover images with dose reduction
Reducing the number of projections is an important strategy to reduce image time and radiation dose, giving the few-view problem. To evaluate the proposed algorithm for few-view tomography, the number of projections around 360° was down-sampled from 720 to 120 and 36, respectively. The PPCS-HS-W method was then applied. Also, the FDK, iRadon, OS-SART and ASD-POCS methods were performed for comparison. The results are in [Figure 16]. The number of iterations was set to 50 in all three iterative methods. | Figure 16: The reconstructed images of the Head phantom (with size of 64 × 64 × 6, using different algorithms, as a function of 36, 120, and 720 projections. The 2nd column is related to 22th vertical slice of the image, unlike the others that are 22th horizontal slices
Click here to view |
It is seen that the performance of the FDK and iRadon reconstruction increased and improved as the number of views rose from 36 to 720. The findings of PPCS-HS-W were far better than the rest of the reconstructions. The PPCS-HS-W result was nearly as good as that reconstructed from 720 views in the case of 36 views [column two of [Figure 16]]. We note that the OS-SART and ASD-POCS have patchy artifacts and wash out the fine structures in low-dose setting, while the higher-degree scheme, PPCS-HS-W, provided very similar results with more precise reconstruction (see the regions indicated by the red arrows).
As shown in [Figure 16], the PPCS-HS-W algorithm achieved a reasonable image quality comparable to the image reconstructed by the FDK using all 720 projections, even with the dose reduction to 1/20th (36/720 projections). But, for a clinical patient case, about 120 projections or more were necessary to generate a reasonable quality images.
This achievement still represents a significant reduction in dose (120/720≅83%), but any further reduction in dose (i.e., less projections) is usually not recommended due to a rapid deterioration of image quality, as some structures in the PPCS-HS-W result may have been obscure or invisible with 36 projections. The possible reason that more projections are required in patients than in phantoms is that the internal anatomy of humans is much less detailed and therefore needs more information to depict it properly.
We presented a new rapid CBCT image reconstruction method which combines 3D PPFT and CS. Using PPFT would avoid repeated regridding/inverse-regridding operations that may give rise to interpolation errors in reconstruction based on Cartesian FFT. Conducted in three simulated phantoms, our studies demonstrated that HS penalty can wipeout the staircase effect while preserving both edges and regions with smooth intensity transition or piecewise constant. A significant advantage with using the PP based Radon transform is reducing the computational complexity. In addition, solving CS problem with FCSA for an image of size N×N×N voxels, requires a number of O(N3log(N3)) operations in each iteration and can obtain an ∈-optimal solution in  iterations and is therefore faster than ordinary iterative methods.[14]
Not only did the proposed method minimize the computational complexity, but it also decreased the dose of X-ray radiation by reducing the number of predictions. Nonetheless, approximately 120 calculations are required for a medical patient case to produce images of reasonable quality. This achievement represents a significant dose reduction of 83 percent, but any further dose reduction (i.e., lower number of projections) is usually not recommended due to a rapid deterioration in some structures of the image contrast and performance.
The future work is implementing of PPCS-HS-W in C# with the CUDA programming environment (NVIDIA, Santa Clara, CA) and utilizing the massive parallel computational capability of the GPU hardware to make it usable in clinical application and the large size 3D images.
Appendix
Fast alternating minimization algorithm for hessian regularized inverse problems
We change (25) as an image recovery problem with the Hessian penalty like:
Duf denotes the multi-dimensional convolution of f with Du discrete filter.[49]
Given the Huber function and its half-quadratic dual [49], (36) can be transformed to (37).
Where Zu is an auxiliary function. In order to solve (37), we use the alternative minimization method, which produces two efficiently solved subproblems: the z-subproblem that can be solved in one shrinkage step; and the f-subproblem involving the reversal of a linear process that can often be solved in one step using DFTs. Below are the specifics of these two sub problems.
- The z-Subproblem: Minimization With Respect to z, Assuming f to be Fixed
With the assumption of fixed f, we can minimize the cost function in (37) with respect to z in order to have
The soft-shrinkage formula provides the exact component-wise solution to the above problem:
- The f-Subproblem: Minimization With Respect to f, Assuming z to be Fixed
With the assumption of fixed z, we minimize (37) with respect to f, which leads to
The objective described above is quadratic in f, and from the normal equations its minimizer satisfies the following:
In general, a simple matrix solver such as the conjugate gradient algorithm can be used to solve the linear Eq. 41. It can be simplified by applying the DFT on both sides.
Where I is the identity matrix, E is circulant matrix corresponding to convolution with discretized partial derivatives, Du= ES. It should be note that S is steering function and C= STS. Finally, q, the projected shrinkage defined as  where the samples ui∈S and weights wi are determined by the choice of quadrature; more details on the choice and performance of specific quadratures are given at.[49]
The FAM pseudocode is shown below.
Algorithm 2: Pseudo code of fast alternating minimization to solve the hessian regularized problem(25)[49]
Initialization: {Input N_outer, N_inner, λ, λfac, siz = size (y), tol, nsamples = number of samples in Lebedev quadrature scheme.
Define derivative operators and steering functions:{
- G = Compute 3D filters cooresponding to discrete derivative operators defined by the tensor product of derivative of 1D B-spline operators.
- nfilters = size of G.
- D = Build derivatve operators from filters G.
- Dt = Adjoint of D (t is adjoint operator).
- (pts, weights) = Get quadrature samples/weights from nsamples.
- su = Compute steering polynomials of pts for Hessian penalties.
- su = su× weights.}
Precompute quantites:{
- C = tensor product between ∫S sutsu and identity matrix.
- Dvec = resize D (im0) to siz× nfilters.
- CD = resize (Dvec×C) to siz×nfilters.
- DtCD = FFT (Dt (CD)).}
- f = y: initialize f
Begin alternating minimization algorithm:
- e = [0]N_outer×N_inner, Ind = 0.
- for i = 1:N_outer
- for ii = 1:N_inner
Calculate error: {
- DF = Compute 4D array of all partial derivatives of f.
- DFvec = resize DF to siz× nfilters.
- DD = DFvec×su: Compute directional derivatives of f.
- diff = f-y.
- e=|diff2+αΣ|DD||: objective function
- ind = ind + 1, e (ind) = e.
- If | e (ind)-e (ind-1) |<tol, then break: Convergence test.}
Shrinkage step:{
- shrinkDD = shrink (DD,1/λ).
- Compute projected shrinkage step: DD = shrinkDD×DD.
- DF = resize DD×su to siz×nfilters.
Inversion step: {
- F1 = FFT3 ((2y + αλDt (DF)).
- F2 = αλ × DtCD + 2.
- f = iFFT3(F1/F2).}
- λ = λ ×λfac: Increase continuation parameter}
Financial support and sponsorship
None.
Conflicts of interest
There are no conflicts of interest.
| References | | |
| 1. | Hashemi S, Beheshti S, Gill PR, Paul NS, Cobbold RS. Efficient low dose X-ray CT reconstruction through sparsity-based MAP modeling. arXiv Prepr 2014;abs/1402.1801:1-10.  |
| 2. | Zhang C, Zhang T, Li M, Peng C, Liu Z, Zheng J. Low-dose CT reconstruction via L1 dictionary learning regularization using iteratively reweighted least-squares. Biomed Eng Online 2016;15:66. |
| 3. | Hou W, Zhang C. Analysis of compressed sensing based CT reconstruction with low radiation. In: 2014 International Symposium on Intelligent Signal Processing and Communication Systems (ISPACS). Malaysia. IEEE; 2014. p. 291-6. |
| 4. | Liu L, Li X, Xiang K, Wang J, Tan S. Low-dose CBCT reconstruction using hessian schatten penalties. IEEE Trans Med Imaging 2017;36:2588-99. |
| 5. | Xu Q, Yu H, Mou X, Zhang L, Hsieh J, Wang G. Low-dose X-ray CT reconstruction via dictionary learning. IEEE Trans Med Imaging 2012;31:1682-97. |
| 6. | Bai T, Yan H, Jia X, Jiang S, Wang G, Mou X. Z-index parameterization for volumetric CT image reconstruction via 3-D dictionary learning. IEEE Trans Med Imaging 2017;36:2466-78. |
| 7. | Zhang H, Wang J, Ma J, Lu H, Liang Z. Statistical models and regularization strategies in statistical image reconstruction of low-dose X-ray CT : A survey. arXiv Prepr 2014;1412.1732 :1-67. |
| 8. | Wang AS, Stayman JW, Otake Y, Kleinszig G, Vogt S, Siewerdsen JH. Nesterov's method for accelerated penalized-likelihood statistical reconstruction for C-arm cone-beam CT. In: Image Formation in X-Ray Computed Tomography. Curran Associates, Inc. 57 Morehouse Lane Red Hook, NY 12571; 2014. p. 409-13. |
| 9. | Hashemi S, Beheshti S, Gill PR, Paul NS, Cobbold RS. Accelerated compressed sensing based CT image reconstruction. Comput Math Methods Med 2015;2015:161797. |
| 10. | Sukovic P, Clinthorne NH. Penalized weighted least-squares image reconstruction for dual energy X-ray transmission tomography. IEEE Trans Med Imaging 2000;19:1075-81. |
| 11. | Mao Y, Fahimian BP, Osher SJ, Miao J. Development and optimization of regularized tomographic reconstruction algorithms utilizing equally-sloped tomography. IEEE Trans Image Process 2010;19:1259-68. |
| 12. | Donoho DL, Donoho DL. Compressed sensing. IEEE Trans Inf Theory 2016;52:1289-306. |
| 13. | Yang Y, Liu F, Li M, Jin J, Weber E, Liu Q, et al. Pseudo-polar Fourier transform-based compressed sensing MRI. IEEE Trans Biomed Eng 2017;64:816-25. |
| 14. | Huang J, Zhang S, Metaxas D. Efficient MR image reconstruction for compressed MR imaging. Med Image Anal 2011;15:670-9. |
| 15. | Lee MS, Kim HJ, Cho HM, Hong DK, Je UK. Compressed-sensing (CS)-based 3D image reconstruction in cone-beam CT (CBCT) for low-dose, high-quality dental X-ray imaging. J Korean Phys Soc 2013;63:1066-71. |
| 16. | Kaldate A, Patre BM, Harsh R, Verma D. MR Image Reconstruction Based on Compressed Sensing using Poisson sampling pattern. In: 2016 Second International Conference on Cognitive Computing and Information Processing (CCIP). India: IEEE; 2016. p. 1-4. |
| 17. | Hashemi S, Beheshti S, Cobbold RS, Paul NS. Subband-dependent compressed sensing in local CT reconstruction. Signal Image Video Process 2016;10:1009-15. |
| 18. | Xu D, Huang Y, Kang JU. Volumetric (3D) compressive sensing spectral domain optical coherence tomography. Biomed Opt Express 2014;5:3921-34. |
| 19. | Ragab M, Omer OA, Abdel-Nasser M. Compressive sensing MRI reconstruction using empirical wavelet transform and grey wolf optimizer. Neural Comput Appl 2018;7:2705-724. |
| 20. | Bonnet S, Peyrin F, Turjman F, Prost R. Nonseparable wavelet-based cone-beam reconstruction in 3-D rotational angiography. IEEE Trans Med Imaging 2003;22:360-7. |
| 21. | Fahimian BP, Zhao Y, Huang Z, Fung R, Mao Y, Zhu C, et al. Radiation dose reduction in medical x-ray CT via Fourier-based iterative reconstruction. Med Phys 2013;40:031914. |
| 22. | Matej S, Fessler JA, Kazantsev IG. Iterative tomographic image reconstruction using Fourier-based forward and back-projectors. IEEE Trans Med Imaging 2004;23:401-12. |
| 23. | Averbuch A, Shkolnisky Y. 3D Fourier based discrete Radon transform. Appl Comput Harmon Anal 2003;15:33-69. |
| 24. | Averbuch A, Shkolnisky Y. 3D discrete X-ray transform. Applied and Computational Harmonic Analysis 2004;17:259-76. |
| 25. | Averbuch A, Coifman RR, Donoho DL, Elad M, Israeli M. Analysis, “Fast and accurate Polar Fourier transform”. Appl Comput Harmon Anal 2006;21:145-67. |
| 26. | Averbuch A, Sedelnikov I, Shkolnisky Y. CT reconstruction from parallel and fan-beam projections by a 2-D discrete Radon transform. IEEE Trans Image Process 2012;21:733-41. |
| 27. | Keller Y, Shkolnisky Y, Averbuch A. Volume registration using the 3-D Pseudopolar Fourier transform. IEEE Trans Signal Process 2006;54:4323-31. |
| 28. | Teyfouri N, Rabbani H, Kafieh R, Jabbari I. “An Exact and Fast CBCT Reconstruction via Pseudo-Polar Fourier Transform-Based Discrete Grangeat's Formula.” IEEE Transactions on Image Processing 2020;29: 5832-5847. |
| 29. | Qiaoqiao DD, Yong L, Zhang X, Fessler JA. Statistical image reconstruction using mixed poisson-gaussian noise model for X-ray CT. arXiv 2018;1801.09533:1-11. |
| 30. | Thibault JB, Sauer KD, Bouman CA, Hsieh J. A three-dimensional statistical approach to improved image quality for multislice helical CT. Med Phys 2007;34:4526-44. |
| 31. | Wang G, Zhou J, Yu Z, Wang W, Qi J. Hybrid pre-log and post-log image reconstruction for computed tomography. IEEE Trans Med Imaging 2017;36:2457-65. |
| 32. | Thibault JB, Bouman CA, Sauer KD, Hsieh J. A recursive filter for noise reduction in statistical iterative tomographic imaging. Proc SPIE 2006;6065:1-10. |
| 33. | Bouman CA, Sauer K. A unified approach to statistical tomography using coordinate descent optimization. IEEE Trans Image Process 1996;5:480-92. |
| 34. | Sauer K, Bouman C. A local update strategy for iterative reconstruction from projections. IEEE Trans Signal Process1993;41:534-48. |
| 35. | |
| 36. | Zhang H, Wang J, Zeng D, Tao X, Ma J. Regularization strategies in statistical image reconstruction of low-dose x-ray CT: A review. Med Phys 2018;45:e886-907. |
| 37. | |
| 38. | Flores L, Vidal V, Verdú G. Iterative reconstruction from few-view projections. Procedia-Procedia Comput Sci 2015;51:703-12. |
| 39. | Ha S, Matej S, Ispiryan M, Mueller K. “GPU-Accelerated Forward and Back-Projections with Spatially Varying Kernels for 3D DIRECT TOF PET Reconstruction,” IEEE Transactions on Nuclear Science 2013;60:166-73. |
| 40. | Riabkov D, Xue X, Tubbs D, Cheryauka A. Accelerated cone-beam back-projection using GPU-CPU hardware, in Proc. 9th Int. Meeting Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine (Fully3D '07), China; 2007. p. 68-71. |
| 41. | Gottlieb D. On the direct Fourier method for computer tomography. IEEE Trans Med Imaging 2000;19:223-32. |
| 42. | Miao J, Förster F, Levi O, Förster F, Levi O. Equally sloped tomography with oversampling reconstruction. Phys Rev B 2005;72:3-6. |
| 43. | Dusaussoy NJ. VOIR: A volumetric image reconstruction algorithm based on Fourier techniques for inversion of the 3-D radon transform. IEEE Trans Image Process 1996;5:121-31. |
| 44. | Axelsson C, Danielsson PE. Three-dimensional reconstruction from cone-beam data in O (N3 log N) time. Phys Med Biol 1994;39:477-91. |
| 45. | Averbuch A, Shabat G, Shkolnisky Y. Direct inversion of the 3D pseudo-polar Fourier transform. SIAM J Sci Comput 2016;38:A1100-A1120. |
| 46. | Beck A, Teboulle M. A Fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM J Imaging Sci 2009;2:183-202. |
| 47. | Sun T, Sun N, Wang J, Tan S. Iterative CBCT reconstruction using Hessian penalty. Phys Med Biol 2015;60:1965-87. |
| 48. | Hu Y, Jacob M. Higher degree total variation (HDTV) regularization for image recovery. IEEE Trans Image Process 2012;21:2559-71. |
| 49. | Hu Y, Ongie G, Ramani S, Jacob M. Generalized higher degree total variation (HDTV) regularization. IEEE Trans Image Process 2014;23:2423-35. |
| 50. | Chen B, Xiang K, Gong Z, Wang J, Tan S. Statistical Iterative CBCT Reconstruction Based on Neural Network. IEEE Trans Med Imaging 2018;37:1511-21. |
| 51. | Fletcher AK, Ramchandran K, Goyal VK. Wavelet denoising by recursive cycle spinning. In: International Conference on Image Processing. USA: IEEE; 2002. p. 2. |
| 52. | |
| 53. | Tsutomu G, Nakajima M. Dual-energy subtraction X-ray digital tomosynthesis: Basic physical evaluation. J Med Imaging 2012;2:111-17. |
| 54. | Wang Z, Bovik AC, Sheikh HR, Simoncelli EP. Image quality assessment: From error visibility to structural similarity. IEEE Trans Image Process 2004;13:600-12. |
| 55. | |
| 56. | Smith SW. The Scientist and Engineer's Guide to Digital Signal Processing. San Diego, CA, USA: California Technical Publishing; 1997. |
| 57. | Wang G, Jiang M. Ordered-subset simultaneous algebraic reconstruction techniques (OS-SART). X-Ray Science and Technology 2003;12:169-177. |
| 58. | Sidky EY, Pan X. Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization. Phys Med Biol 2008;53:4777-807. |
| 59. | Biguri A, Dosanjh M, Hancock S, Soleimani M. TIGRE: A MATLAB-GPU toolbox for CBCT image reconstruction. Biomed Phys Eng Express 2016;2:055010. |
| 60. | Maier A, Hofmann HG, Berger M, Fischer P, Schwemmer C, Wu H, et al. CONRAD--a software framework for cone-beam imaging in radiology. Med Phys 2013;40:111914. |
| 61. | Rit S, Oliva MV, Brousmiche S, Labarbe R, Sarrut D, Sharp GC. The reconstruction toolkit (RTK), an open-source cone-beam CT reconstruction toolkit based on the insight toolkit (ITK). J Phys Conf Ser 2014;489:12079. |
| 62. | van Aarle W, Palenstijn WJ, De Beenhouwer J, Altantzis T, Bals S, Batenburg KJ, et al. The ASTRA toolbox: A platform for advanced algorithm development in electron tomography. Ultramicroscopy 2015;157:35-47. |
[Figure 1], [Figure 2], [Figure 3], [Figure 4], [Figure 5], [Figure 6], [Figure 7], [Figure 8], [Figure 9], [Figure 10], [Figure 11], [Figure 12], [Figure 13], [Figure 14], [Figure 15], [Figure 16]
[Table 1], [Table 2], [Table 3], [Table 4]
|
Search Pubmed for