

ORIGINAL ARTICLE 

Year : 2022  Volume
: 12
 Issue : 1  Page : 824 

Lowdose conebeam computed tomography reconstruction through a fast threedimensional compressed sensing method based on the threedimensional pseudopolar fourier transform
N Teyfouri^{1}, Hossein Rabbani^{1}, I Jabbari^{2}
^{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  30Dec2019 
Date of Decision  24Apr2021 
Date of Acceptance  20Aug2021 
Date of Web Publication  28Dec2021 
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 PPbased compressed sensingHSW. In the reconstruction range of 120 projections around the 360° rotation, the image quality is visually similar to reconstructed images by FeldkampDavisKress 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 highresolution, lownoise reconstruction results in lowdose level are also shown.
Keywords: Threedimensional compressed sensing, threedimensional pseudopolar Fourier transform, cone beam computed tomography reconstruction, hessian
How to cite this article: Teyfouri N, Rabbani H, Jabbari I. Lowdose conebeam computed tomography reconstruction through a fast threedimensional compressed sensing method based on the threedimensional pseudopolar fourier transform. J Med Signals Sens 2022;12:824 
How to cite this URL: Teyfouri N, Rabbani H, Jabbari I. Lowdose conebeam computed tomography reconstruction through a fast threedimensional compressed sensing method based on the threedimensional pseudopolar fourier transform. J Med Signals Sens [serial online] 2022 [cited 2022 Jan 26];12:824. Available from: https://www.jmssjournal.net/text.asp?2022/12/1/8/334132 
Introduction   
Today, Xray 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 Xray 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 Xray imaging, the signal to noise ratio (SNR) has quadratic dependency on the dose of Xray. As other conditions are identical, any reduction in the Xray 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 Xray 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 Xray 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 fewview, limitedangle, 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 maximumlikelihood or penalizedlikelihood 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 highquality images from roughly onetenth of the number of views required by filtered back projection (FBP) in twodimensional (2D) images, thereby allowing a much lower dose scanning protocol than required in applying conventional reconstruction methods.^{[5]} However, a major drawback with CSbased reconstruction algorithms is their being prohibitively computationally intensive for clinical use.^{[11],[21]}
In this paper, Fourierbased reconstruction algorithms (FRAs), potentially enjoying highspeed 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 parallelray 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 backprojection.
Interpolation errors are defined as limitation of Fourierbased methods of reconstruction. Our method, pseudopolar based compressed sensing (PPCS), relies on as threedimensional (3D) pseudopolar Fourier transform (PPFT), which is a new form of FFT. In PPFT, the grid points in Fourier domain lie on equallysloped lines rather than on equally angled lines. PPFT removes the errors associated with iterative regridding 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 CSbased 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 Xray 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 PPFTbased 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 PPCSHSW.
The rest of the article is organized as follows: Section II outlines different stages of the method including stepbystep 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 l_{1}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 Xray 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 xrays and voxels in the image f, i_{T} and j_{T} are the total number of projection measurements and image voxels, respectively. represents the vector of attenuation coefficients of the object to be reconstructed and y_{i}=[ϕ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 I_{i} that are related to the line integral projections by Beer's law of attenuation. They represent the detected xray intensity after attenuation by the scanned object, and follow a Poisson distribution. m_{i} 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. I_{T} denotes the impinging Xray photon intensity emitted from the Xray.^{[29]}
Modeling the electronic noise is significant in lowdose 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,. d_{i} (y^{3}) 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}. d_{i} 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 I_{T} is the number of radiated photons from the Xray source. I_{i} 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}, I_{i},^{[9]} d_{i} can be defined as:^{[36]}
The tomography operator, ϕ in y = ϕf+w, is illposed 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 CSbased reconstruction methods, projection operator or 3D Radon transform, ϕ consist of integrals over planes through image data points and is a (N_{I}×N_{J}×N_{K}) 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 threedimensional radon space. (b) Sinogram Sampling on a conventional threedimensional 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 FeldkampDavisKress (FDK)based methods. To do the CSbased 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 pseudopolar 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 pseudopolar grid. (a) Threedimensional conventional Fourier slice theorem on spherical coordinate, sinogram sampling in the Radon space. (b) Threedimensional Fourier slice theorem on pseudopolar 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 pseudopolar (PP) grid with complexity of O(N^{3}logN) [Figure 2]b. PP grid shown in [Figure 2]b would be of help in bypassing the interpolation stage.^{[28],[43],[44]} The PPFTbased 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) Twodimensional cartesian grid. (b) Twodimensional PP grid. (c) Threedimensional pseudopolar grid
Click here to view 
In reconstructing undersampled radial data with CS, we would need regridding and inverseregridding to transfer data between the frequency and image domains. In each CS iteration, 3D interpolations are implemented twice in the regridding and inverseregridding. That would give rise to errors and affect the reconstruction quality. To resolve these problems, a radiallike 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 PPbased FFT and its inverse, setting the gridding error to zero, (2) PPFT has a fastforward and a fastbackward calculation algorithm,^{[23],[28],[45]} enabling our proposed algorithm to avoid the regridding step used in iterative nonCartesian Fourierbased 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 l_{1} norm term, effectively. FCSA is motivated by FISTA^{[46]} which minimizes the following relation:
Where P denotes a smooth convex function with Lipschitz constant L_{P}, 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 f^{k} = prox_{ρ}(Q)(f_{Q}) 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 l_{1} regularization problem as the second step f^{k}= prox_{ρ}(β‖Ψf‖_{1})(f_{Q}) has a close form solution.
The HS regularization problem, f^{k}= prox_{ρ}(α‖Ψf‖_{1})(f_{Q}) is solved by the fast alternating minimization (FAM) algorithm.^{[4],[47],[48],[49]}
However, FISTA does not solve the composite l_{1} 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).
Threedimensional 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 firstorder derivatives results in the unwanted staircase artifact that would often make reconstructed images oversharpened and unnatural^{[4],[48],[49],[50]} In this article, we suggest a secondorder 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 l_{1} norm of the gradient magnitude,^{[49]} specified as:
Isotropic Higher degree TV (HDTV) regularization is specified by (21).
where f_{u,n} is the nth degree directional derivative defined as:
Given a specific nth degree differential operator D=∑_{α=2}c_{α} ∂^{α}, (α is a multiindex and the c_{α} are constants), we define the generalized HDTV penalty as:
where and D_{u} 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.
Threedimensional 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) Threedimensional 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 D_{t}:(s_{t},α_{t},ψ_{t}), where (s_{t},α_{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 CSbased 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(N^{3}logN).
To generate initial image for Algorithm 1, use of PPFT based method makes it possible to reduce the arithmetic operations from O(N^{4}) with the conventional FDK method^{[44]} to O(N^{3}log(N))^{[23],[28]} for an image with N × N × N voxels. This result in a faster convergence compared with setting f^{1}=FDK.
FCSA solver in the PPCSHSW 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(N^{3}).
The cost of FISTA solver is O(N^{3} log(N^{3})) in each iteration. The step 4, 6 and 7 of “for” loop of algorithm only involve adding vectors or scalars, thus cost only O(N^{3}). In the step 1, since in this case. Thus, this step only costs O(N^{3} log(N^{3})). As introduced above, both of two steps f^{k}= prox_{ρ}(2α‖f‖_{Hess})(f_{g}) and f^{k} = prox_{ρ}(2β‖Ψf‖_{1})(f_{g}) has a close form solution and can be computed with cost O(N^{3} log(N^{3})) 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 {f^{k}} is the sequence generated by the , then, f^{k} will converge to .
Summary of pseudopolar based compressed sensingHSW algorithm
After combining the entire process together, Algorithm 1 will elaborate the proposed method.
Algorithm 1: Summary of pseudopolar based compressed sensingHSW 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 SheppLogan phantom. [Figure 5]a demonstrates a representative slice of a noisefree 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 contrasttonoise ratio (CNR), as well as a ball with linear gradual transitioning intensity resembling the octahedron.
To demonstrate the performance of PPCSHSW 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 Xray intensity, I_{0}, and the background electronic noise variance, σ, are set to 0.6 × 10^{5} 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. C_{1} and C_{1} are two constants to avoid instability. In this study, C_{1} and C_{1} were chosen as C_{2}=(0.01μmax)^{2} and C_{1}=(0.01μmax)^{2}. While the two windows move pixelbypixel 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 PPCSHSW, we have compared it with some other algorithms. First orderedsubset simultaneous algebraic reconstruction techniques (OSSART)^{[57]} was implemented. second the adaptive steepest descent projections onto convex set (ASDPOCS)^{[58],[59]} was implemented. Third and fourth, PPCSTVW and PPCSHS method were presented to demonstrate the effect of Hessian and l_{1} 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 SARTtype 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. OSSART updates the image with a subset of the projections.^{[57]}
ASDPOCS 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 Xrays and voxels in the image.
In PPCSTVW we have:
Where,
We solved (32) by FCSA including two FISTA solvers.^{[46]}
In PPCSHSW without l1 norm named PPCSHS, 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 PPCSHSW 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, PPCSTVW, PPCSHS, and PPCSHSW 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) Pseudopolar based compressed sensingTVW. (d) Pseudopolar based compressed sensingHS. (e) Pseudopolar based compressed sensing HSW. (f) iRadon. (g) Adaptive steepest descent projections onto convex set. (h) Pseudopolar based compressed sensingTVW. (i) Pseudopolar based compressed sensingHS. (j) Pseudopolar based compressed sensingHSW. (k) iRadon. (l) Adaptive steepest descent projections onto convex set. (m) Pseudopolar based compressed sensingTVW. (n) Pseudopolar based compressed sensingHS. (o) Pseudopolar based compressed sensingHSW
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 2^{nd} and 3^{rd} 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 PPCSHS look more natural and smoother, emphasizing on ability of PPCSHS in preserving smooth regions as well as suppressing the noise.
Compared to PPCSHS, using a sparsity regularization in an orthogonal basis was made to better reconstruct edges in the PPCSHSW 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 PPCSHSW with sparse views and the fixed incident Xray intensity I_{0} of 0.6 × 10^{5} 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, PPCSHSW outperformed others. PPCSHSW 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 PPCSHSW 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) Pseudopolar based compressed sensingTVW. (e) Pseudopolar based compressed sensingHS. (f) Pseudopolar based compressed sensingHSW
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 PPCSHSW was whiter, indicating a higher ability in preserving structures.  Figure 8: Structural similarity index map of the head phantom using different methods. (a) Pseudopolar based compressed sensing TVW (structural similarity index = 0.65). (b) Pseudopolar based compressed sensing HSW (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 edgepreservation, 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 PPCSTVW and PPCSHSW were 2.35 and 1.2, respectively. However, PPCSHSW 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, ASDPOCS and PPCSTVW.  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, Pseudopolar based compressed sensingTVW [Figure 11]c and Pseudopolar based compressed sensingHSW [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) Pseudopolar based compressed sensingTVW. (d) Pseudopolar based compressed sensingHSW
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: Contrasttonoise 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 PPCSHSW 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 PPCSHSW. 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) Pseudopolar based compressed sensingTVW. (d) Pseudopolar based compressed sensingHSWv
Click here to view 
Ability to remove blocky artifact
[Figure 11]c and [Figure 11]d shows the ability of PPCSHSW method compared to PPCSTVW 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 TVbased reconstruction produced several staircase artifacts, while the reconstructed results using HS penalties were far smoother [Figure 13].
The regions of smoothlychanged intensity, such as octahedron in the upperright corner and the sphere in the lowerleft 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 SheppLogan 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 frequencydomain methods in increasing the speed of iterative regularization problem solving, we made a comparison between PPCSTVW and PPCSHSW and two algorithms in the spatial domain, OSSART^{[57]} and ASDPOCS.^{[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 PPCSHSW >PPCSTVW >ASDPOCS >OSSART. This finding proves to be correct at all levels of iterations, as shown in [Figure 14].  Figure 14: The reconstructed images of the SheppLogan phantom, using different algorithms, as a function of 10, 30, and 50 iterations. A total of 36 projections in conebeam geometry were used for the reconstructions
Click here to view 
As [Figure 15] indicates, the PPCSTVW and PPCSHSW algorithms show convergence at about 20 and 15 iterations, respectively. That is while the OSSART and ASDPOCS algorithms need further convergence at 50 iterations.  Figure 15: Meansquared 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 f_{i,j,k} denotes the voxel values in the reconstructed volume f and represents the original values of the SheppLogan 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 OSSART and ASDPOCS 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 PPCSTVW, solved by FSCA in,^{[9]} by Bregman iterative regularization in^{[11]} and by nonlocal TV regularization in^{[21]} has been shown that PPCSTVW is faster than FBP and algebraic reconstruction techniqueTV based method for fan beam images reconstruction. In addition, Huang^{[14]} used Fourier based CSTVW 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 PPCSTVW and PPCSHSW 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 realtime 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 fewview problem. To evaluate the proposed algorithm for fewview tomography, the number of projections around 360° was downsampled from 720 to 120 and 36, respectively. The PPCSHSW method was then applied. Also, the FDK, iRadon, OSSART and ASDPOCS 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 2^{nd} column is related to 22^{th} vertical slice of the image, unlike the others that are 22^{th} 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 PPCSHSW were far better than the rest of the reconstructions. The PPCSHSW 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 OSSART and ASDPOCS have patchy artifacts and wash out the fine structures in lowdose setting, while the higherdegree scheme, PPCSHSW, provided very similar results with more precise reconstruction (see the regions indicated by the red arrows).
As shown in [Figure 16], the PPCSHSW 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/20^{th} (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 PPCSHSW 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/inverseregridding 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(N^{3}log(N^{3})) 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 Xray 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 PPCSHSW 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:
D_{u}f denotes the multidimensional convolution of f with D_{u} discrete filter.^{[49]}
Given the Huber function and its halfquadratic dual [49], (36) can be transformed to (37).
Where Z_{u} is an auxiliary function. In order to solve (37), we use the alternative minimization method, which produces two efficiently solved subproblems: the zsubproblem that can be solved in one shrinkage step; and the fsubproblem 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 zSubproblem: 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 softshrinkage formula provides the exact componentwise solution to the above problem:
 The fSubproblem: 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, D_{u}= ES. It should be note that S is steering function and C= S^{T}S. Finally, q, the projected shrinkage defined as where the samples u_{i}∈S and weights w_{i} 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 Bspline operators.
 nfilters = size of G.
 D = Build derivatve operators from filters G.
 D_{t} = 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} su_{t}su and identity matrix.
 Dvec = resize D (im0) to siz× nfilters.
 CD = resize (Dvec×C) to siz×nfilters.
 D_{t}CD = FFT (D_{t} (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 = fy.
 e=diff^{2}+αΣDD: objective function
 ind = ind + 1, e (ind) = e.
 If  e (ind)e (ind1) <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 + αλD_{t} (DF)).
 F2 = αλ × D_{t}CD + 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 Xray CT reconstruction through sparsitybased MAP modeling. arXiv Prepr 2014;abs/1402.1801:110. 
2.  Zhang C, Zhang T, Li M, Peng C, Liu Z, Zheng J. Lowdose CT reconstruction via L1 dictionary learning regularization using iteratively reweighted leastsquares. 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. 2916. 
4.  Liu L, Li X, Xiang K, Wang J, Tan S. Lowdose CBCT reconstruction using hessian schatten penalties. IEEE Trans Med Imaging 2017;36:258899. 
5.  Xu Q, Yu H, Mou X, Zhang L, Hsieh J, Wang G. Lowdose Xray CT reconstruction via dictionary learning. IEEE Trans Med Imaging 2012;31:168297. 
6.  Bai T, Yan H, Jia X, Jiang S, Wang G, Mou X. Zindex parameterization for volumetric CT image reconstruction via 3D dictionary learning. IEEE Trans Med Imaging 2017;36:246678. 
7.  Zhang H, Wang J, Ma J, Lu H, Liang Z. Statistical models and regularization strategies in statistical image reconstruction of lowdose Xray CT : A survey. arXiv Prepr 2014;1412.1732 :167. 
8.  Wang AS, Stayman JW, Otake Y, Kleinszig G, Vogt S, Siewerdsen JH. Nesterov's method for accelerated penalizedlikelihood statistical reconstruction for Carm conebeam CT. In: Image Formation in XRay Computed Tomography. Curran Associates, Inc. 57 Morehouse Lane Red Hook, NY 12571; 2014. p. 40913. 
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 leastsquares image reconstruction for dual energy Xray transmission tomography. IEEE Trans Med Imaging 2000;19:107581. 
11.  Mao Y, Fahimian BP, Osher SJ, Miao J. Development and optimization of regularized tomographic reconstruction algorithms utilizing equallysloped tomography. IEEE Trans Image Process 2010;19:125968. 
12.  Donoho DL, Donoho DL. Compressed sensing. IEEE Trans Inf Theory 2016;52:1289306. 
13.  Yang Y, Liu F, Li M, Jin J, Weber E, Liu Q, et al. Pseudopolar Fourier transformbased compressed sensing MRI. IEEE Trans Biomed Eng 2017;64:81625. 
14.  Huang J, Zhang S, Metaxas D. Efficient MR image reconstruction for compressed MR imaging. Med Image Anal 2011;15:6709. 
15.  Lee MS, Kim HJ, Cho HM, Hong DK, Je UK. Compressedsensing (CS)based 3D image reconstruction in conebeam CT (CBCT) for lowdose, highquality dental Xray imaging. J Korean Phys Soc 2013;63:106671. 
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. 14. 
17.  Hashemi S, Beheshti S, Cobbold RS, Paul NS. Subbanddependent compressed sensing in local CT reconstruction. Signal Image Video Process 2016;10:100915. 
18.  Xu D, Huang Y, Kang JU. Volumetric (3D) compressive sensing spectral domain optical coherence tomography. Biomed Opt Express 2014;5:392134. 
19.  Ragab M, Omer OA, AbdelNasser M. Compressive sensing MRI reconstruction using empirical wavelet transform and grey wolf optimizer. Neural Comput Appl 2018;7:2705724. 
20.  Bonnet S, Peyrin F, Turjman F, Prost R. Nonseparable waveletbased conebeam reconstruction in 3D rotational angiography. IEEE Trans Med Imaging 2003;22:3607. 
21.  Fahimian BP, Zhao Y, Huang Z, Fung R, Mao Y, Zhu C, et al. Radiation dose reduction in medical xray CT via Fourierbased iterative reconstruction. Med Phys 2013;40:031914. 
22.  Matej S, Fessler JA, Kazantsev IG. Iterative tomographic image reconstruction using Fourierbased forward and backprojectors. IEEE Trans Med Imaging 2004;23:40112. 
23.  Averbuch A, Shkolnisky Y. 3D Fourier based discrete Radon transform. Appl Comput Harmon Anal 2003;15:3369. 
24.  Averbuch A, Shkolnisky Y. 3D discrete Xray transform. Applied and Computational Harmonic Analysis 2004;17:25976. 
25.  Averbuch A, Coifman RR, Donoho DL, Elad M, Israeli M. Analysis, “Fast and accurate Polar Fourier transform”. Appl Comput Harmon Anal 2006;21:14567. 
26.  Averbuch A, Sedelnikov I, Shkolnisky Y. CT reconstruction from parallel and fanbeam projections by a 2D discrete Radon transform. IEEE Trans Image Process 2012;21:73341. 
27.  Keller Y, Shkolnisky Y, Averbuch A. Volume registration using the 3D Pseudopolar Fourier transform. IEEE Trans Signal Process 2006;54:432331. 
28.  Teyfouri N, Rabbani H, Kafieh R, Jabbari I. “An Exact and Fast CBCT Reconstruction via PseudoPolar Fourier TransformBased Discrete Grangeat's Formula.” IEEE Transactions on Image Processing 2020;29: 58325847. 
29.  Qiaoqiao DD, Yong L, Zhang X, Fessler JA. Statistical image reconstruction using mixed poissongaussian noise model for Xray CT. arXiv 2018;1801.09533:111. 
30.  Thibault JB, Sauer KD, Bouman CA, Hsieh J. A threedimensional statistical approach to improved image quality for multislice helical CT. Med Phys 2007;34:452644. 
31.  Wang G, Zhou J, Yu Z, Wang W, Qi J. Hybrid prelog and postlog image reconstruction for computed tomography. IEEE Trans Med Imaging 2017;36:245765. 
32.  Thibault JB, Bouman CA, Sauer KD, Hsieh J. A recursive filter for noise reduction in statistical iterative tomographic imaging. Proc SPIE 2006;6065:110. 
33.  Bouman CA, Sauer K. A unified approach to statistical tomography using coordinate descent optimization. IEEE Trans Image Process 1996;5:48092. 
34.  Sauer K, Bouman C. A local update strategy for iterative reconstruction from projections. IEEE Trans Signal Process1993;41:53448. 
35.  
36.  Zhang H, Wang J, Zeng D, Tao X, Ma J. Regularization strategies in statistical image reconstruction of lowdose xray CT: A review. Med Phys 2018;45:e886907. 
37.  
38.  Flores L, Vidal V, Verdú G. Iterative reconstruction from fewview projections. ProcediaProcedia Comput Sci 2015;51:70312. 
39.  Ha S, Matej S, Ispiryan M, Mueller K. “GPUAccelerated Forward and BackProjections with Spatially Varying Kernels for 3D DIRECT TOF PET Reconstruction,” IEEE Transactions on Nuclear Science 2013;60:16673. 
40.  Riabkov D, Xue X, Tubbs D, Cheryauka A. Accelerated conebeam backprojection using GPUCPU hardware, in Proc. 9th Int. Meeting Fully ThreeDimensional Image Reconstruction in Radiology and Nuclear Medicine (Fully3D '07), China; 2007. p. 6871. 
41.  Gottlieb D. On the direct Fourier method for computer tomography. IEEE Trans Med Imaging 2000;19:22332. 
42.  Miao J, Förster F, Levi O, Förster F, Levi O. Equally sloped tomography with oversampling reconstruction. Phys Rev B 2005;72:36. 
43.  Dusaussoy NJ. VOIR: A volumetric image reconstruction algorithm based on Fourier techniques for inversion of the 3D radon transform. IEEE Trans Image Process 1996;5:12131. 
44.  Axelsson C, Danielsson PE. Threedimensional reconstruction from conebeam data in O (N3 log N) time. Phys Med Biol 1994;39:47791. 
45.  Averbuch A, Shabat G, Shkolnisky Y. Direct inversion of the 3D pseudopolar Fourier transform. SIAM J Sci Comput 2016;38:A1100A1120. 
46.  Beck A, Teboulle M. A Fast iterative shrinkagethresholding algorithm for linear inverse problems. SIAM J Imaging Sci 2009;2:183202. 
47.  Sun T, Sun N, Wang J, Tan S. Iterative CBCT reconstruction using Hessian penalty. Phys Med Biol 2015;60:196587. 
48.  Hu Y, Jacob M. Higher degree total variation (HDTV) regularization for image recovery. IEEE Trans Image Process 2012;21:255971. 
49.  Hu Y, Ongie G, Ramani S, Jacob M. Generalized higher degree total variation (HDTV) regularization. IEEE Trans Image Process 2014;23:242335. 
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:151121. 
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. Dualenergy subtraction Xray digital tomosynthesis: Basic physical evaluation. J Med Imaging 2012;2:11117. 
54.  Wang Z, Bovik AC, Sheikh HR, Simoncelli EP. Image quality assessment: From error visibility to structural similarity. IEEE Trans Image Process 2004;13:60012. 
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. Orderedsubset simultaneous algebraic reconstruction techniques (OSSART). XRay Science and Technology 2003;12:169177. 
58.  Sidky EY, Pan X. Image reconstruction in circular conebeam computed tomography by constrained, totalvariation minimization. Phys Med Biol 2008;53:4777807. 
59.  Biguri A, Dosanjh M, Hancock S, Soleimani M. TIGRE: A MATLABGPU 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. CONRADa software framework for conebeam 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 opensource conebeam 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:3547. 
[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]
