FDK-Type Algorithms with No Backprojection Weight for Circular and Helical Scan CT

We develop two Feldkamp-type reconstruction algorithms with no backprojection weight for circular and helical trajectory with planar detector geometry. Advances in solid-state electronic detector technologies lend importance to CT systems with the equispaced linear array, the planar (flat panel) detectors, and the corresponding algorithms. We derive two exact Hilbert filtered backprojection (FBP) reconstruction algorithms with no backprojection weight for 2D fan-beam equispace linear array detector geometry (complement of the equi-angular curved array detector). Based on these algorithms, the Feldkamp-type algorithms with no backprojection weight for 3D reconstruction are developed using the standard heuristic extension of the divergent beam FBP algorithm. The simulation results show that the axial intensity drop in the reconstructed image using the FDK algorithms with no backprojection weight with circular trajectory is similar to that obtained by using Hu's and T-FDK, algorithms. Further, we present efficient algorithms to reduce the axial intensity drop encountered in the standard FDK reconstructions in circular cone-beam CT. The proposed algorithms consist of mainly two steps: reconstruction of the object using FDK algorithm with no backprojection weight and estimation of the missing term. The efficient algorithms are compared with the FDK algorithm, Hu's algorithm, T-FDK, and Zhu et al.'s algorithm in terms of axial intensity drop and noise. Simulation shows that the efficient algorithms give similar performance in axial intensity drop as that of Zhu et al.'s algorithm while one of the efficient algorithms outperforms Zhu et al.'s algorithm in terms of computational complexity.


Introduction
In divergent beam ramp FBP reconstruction algorithms, it is known that the properties such as noise, resolution and divergent beam artifacts of reconstructed image are influenced by the position-dependent weight in the backprojection. The backprojection weight (which is related inversely to the distance of the source to the point being reconstructed) causes spatially nonuniform distribution of noise and resolution, and the divergent beam artifacts are more pronounced away from the rotation axis [1][2][3]. Many algorithms using different approaches have been proposed to address this issue. Pan [2,4] used shift variant ramp filtering approach while Wang et al. [5] proposed spatially variant weighting in the backprojection. These approaches are computationally expensive. Recently, divergent beam FBP reconstruction algorithms have been proposed with inverse distance weight [6] and no backprojection weight [7,8]. The reconstructed images with these algorithms have better properties in terms of reduced noise level and uniformity in noise distribution and uniformity in resolution as compared to the classical ramp FBP reconstruction algorithm. The fan-beam reconstruction algorithms with inverse distance backprojection weight have been extended to obtain 3D cone-beam (CB) FDK-type approximate reconstruction algorithms with improved noise and resolution properties [9]. On the other hand, algorithms have been developed to improve the noise performance with the modified redundancy weight in the projection domain for circular and helical CB CT [10][11][12]. The fan-beam reconstruction algorithms with no backprojection weight can be extended to obtain 3D CB FDK-type approximate reconstruction algorithms to improve noise and resolution properties.
The most widely used practical method for CB reconstruction has been the approximate algorithm of Feldkamp, Davis and Kress (FDK) [13] originally proposed for circular scan CB reconstruction, and FDK has been extended to helical scan CB geometry by Wang et al. [14]. The circular scan trajectory geometry does not satisfy Tuy's data sufficiency condition on the projections for exact and stable reconstruction [15,16] permitting only approximate reconstruction and also leads to intensity drop in the reconstructed images. Many methods have been suggested to give improvements in the image quality when reconstructing from circular CB data. Some of them use the iteration schemes to reconstruct the image with circular scan FDK as an intermediate result [17,18]. These algorithms are computationally inefficient due to loss of FBP structure. Another approach to reduce the axial intensity drop is given by Hu [19] which adds, to FDK reconstruction, a correction term obtained from unused circular scan CB data. Zhu et al. [20] propose to further add a "missing term" to obtain a more accurate FDK reconstruction. The missing term is obtained by estimating the missing data, required by the data sufficiency condition, from the CB scan projection data.
In this paper, two exact Hilbert FBP reconstruction algorithms with no backprojection weight for 2D fan-beam with equispace linear array detector geometry (complement of the equi-angular curved array detector) have been derived. Using these fan-beam Hilbert FBP inversion algorithms, Feldkamp-type reconstruction algorithms with no backprojection weight for 3D reconstruction for circular and helical scan trajectories are developed. Further, we present efficient algorithms to reduce the axial intensity drop for circular CB CT using the proposed FDK reconstruction algorithms with no backprojection weight. Simulation studies demonstrate the properties of the proposed reconstruction algorithms with no backprojection weight and the efficient algorithms for reducing axial intensity drop in 3D CB CT. We compare the performance of these algorithms with other algorithms in the literature.

Fan-Beam Algorithms with No Backprojection Weight for Equispace Linear Detector
In this section, we revisit the relation between the parallel beam projections and the fan-beam projections through the Hilbert kernel. Using this relation, we derive the fan-beam reconstruction formulae with no backprojection weight for equispace linear array detector, which are counterpart of the algorithms for the curved array detector [7,8].

Parallel Beam
Projections. Let f ( x) denote the 2D object density function to be reconstructed where x = (x, y) T . The parallel-beam projections of f are given by the expression where s ∈ R, ξ = (− sin θ, cos θ) is any unit vector and δ is Dirac delta function. The delta function in (1) selects the value of x that satisfies x · ξ = s, so that p( ξ, s) is the integral of f on the line perpendicular to ξ at a distance s from the origin x = (0, 0) as shown in Figure 1.
The inversion of f from its projection p using the FBP algorithm involves two steps. The first step is to filter the projection p( ξ, s) by applying 1D ramp filter in s, where h R = ∞ −∞ |σ| exp( j2πσs)dσ is the impulse response of |σ| filter. The ramp filter is seen as equivalent to the successive application of a derivative in s and Hilbert transform, where h H (s) = −1/2π ∞ −∞ j sgn(σ) exp( j2πσs)dσ. The second step in the inversion is the backprojection of filtered projections resulting from the application of Hilbert filter. Mathematically,

Fan-Beam
Projections. The fan-beam projections are obtained by using the fan-beam sampling by an equi-angular curved or equispaced collinear detector arrays. The following discussion makes no assumption about detector geometry. The fan-beam projection are measured by moving the Xray source along a circular trajectory as shown in Figure 2.
The source trajectory can be parameterized by an angular parameter β and is given by a(β) = (R 0 cos β, R 0 sin β) where R 0 is the radius of the circle. The fan-beam projections can be represented by where S 1 is the set of all possible unit vectors in 2D space.
International Journal of Biomedical Imaging where This relation is given by Hamaker et al. [21]. The use of this relation has been explained through the parameterization of α = cos γ e w + sin γ e u and ξ = − sin θ e w + cos θ e u [6], where e u and e w are unit vectors and these vectors are illustrated in Figure 3. The proof of (6) can be found in [6,21]. Using (6), we derive two reconstruction algorithms with no backprojection weight for fan-beam equispace linear array detector.
R 0 u /(R 2 0 + u 2 )). The derivative of g H (β, u ) with respect to u and β using chain rule is given by Since the parallel beam parameters are functions of the fanbeam parameters, ∂s/∂u , ∂θ/∂u , and ∂θ/∂β are given by Substituting (10) in (9), we get Substituting (12) in (11), (11) can be written as 4 International Journal of Biomedical Imaging In (4), (∂/∂s)p H is needed to reconstruct the function f , therefore (13) can be rewritten as This result is used in (4) to get the first fan-beam inversion formula with no backprojection weight for equispace linear array detector. The inversion formula is given by where The integration in this formula being over θ, it is a parallelbeam backprojection. Equation (15) is applicable for full scan, short scan, and very short scan. For the short scan and over scan, (15) degenerates to the 180 • parallel-beam reconstruction. In the case of super short scan, You and Zeng [8] suggested the redundancy weight similar to Noo's redundancy weighting. If the scanning range is [0 π], the redundancy weight can be calculated as follows: where wi is equal to The reconstruction formula for this scanning range is given by

Second Fan-Beam Algorithm with No Backprojection
Weight for Equispace Linear Array Detector. Noo et al. [6] have proposed Hilbert FBP fan-beam inversion formula which can also reconstruct with less than short scan. The inversion formula is given by where The ray defined by (β, u ( x, β)) contains the point x. L is the distance from source position obtained by projecting the point x on the ray passing through the origin. Illustration of u and L are shown Figure 3. Since the redundancy weighting w(β, u ) in (18) applies after the Hilbert filtering, it can be applied in the image domain as part of the backprojection step. This is not feasible with the standard FBP formula because redundancy weight is applied before filtering.
Dennerlein et al. [7] have suggested the valid redundancy weight for full scan, where γ = arctan(u /R 0 ). From the geometry shown in Figures 3 and 4, it can be proved that w(β, γ ) + w(β c , γ c ) = 1, where w(β, γ ) and w(β c , γ c ) are the redundancy weights corresponding to the redundant line K which passes through the point x. Equation (21), in terms of equispace array geometry parameters, can be written as where The second fanbeam inversion formula with no backprojection weight for equispace linear array detector is obtained by substituting (22) in (18) and is given as where g H is given by (19) and Thus, the backprojection weight (L) in (18) can be eliminated for each ray that corresponds to a line integral that is measured twice by weighting each ray with function w(β, x) of (21) during the backprojection for short scan or super short scan. However, it does not allow elimination of the weight for singly measured rays (i.e, in case of short scan or a super short scan).

New FDK-Type Formulae
Feldkamp et al. [13] and Wang et al. [14] have developed 3D CB CT reconstruction algorithms for circular and helical trajectories, respectively, from the conventional fan-beam reconstruction formula. Similarly, in this section, we propose International Journal of Biomedical Imaging two new FDK-type inversion formulae for 3D CB reconstruction with the flat panel array detector CB geometry which are obtained by heuristic extension of the fan-beam inversion formulae with no backprojection weight. These formulae are the counterpart of the formulae for cylindrical array detector given by Narasimhadhan et al. [22,23] for circular and helical trajectory CB CT reconstructions.
where D is the source to detector distance.

New FDK Formulae for Helical and Circular CB Geometry.
To simplify the formulae given below, we assume that the flat panel virtual detector is parallel to the actual detector panel at a distance R 0 from a(β) and hence contains the zaxis. Scaling all equations by substituting u → R 0 u/D and v → R 0 v/D is sufficient to describe the case of a flat panel array detector located at a distance D from the vertex.

Second FDK Formula with No Backprojection Weight for
Helical Trajectory (HFDKW-2). The HFDKW-2 formula is obtained from the (23) by standard FDK extension technique f x, y, z where and g(β, u, v) is given by (26). Note that the integration in second formula in (28) is over β and therefore it is a fan-beam backprojection while the integration in first HFDKW formula (25) being over θ results in a parallel-beam backprojection. The first HFDKW formula (25) can be modified to include redundancy weight for partial scan data as given by You and Zeng [8] (25).
The conventional FDK formula involves a weighted backprojection. The weight L −2 in the backprojection, where L is the distance from source position obtained by projecting the point x on the ray passing through the origin, is therefore position-dependent. The weight (R 2 0 + u 2 )/R 3 0 required in the HFDKW-2 can be evaluated and multiplied with filtered projection data before the backprojection. Since there is no position dependent weight in the backprojection step, HFDKW-2 is computationally more efficient than the conventional FDK algorithm. HFDKW-1 algorithm performs an additional interpolation in the backprojection, hence the computational complexity of the HFDKW-1 is more than the conventional HFDK and HFDKW-2 algorithms.

FDK Formulae with No Backprojection Weight for
Circular Trajectory. The HFDKW-1 and HFDKW-2 algorithms can be reduced to circular FDK algorithms with no backprojection weight by setting the helical pitch P = 0. The resulting two algorithms are named as CFDKW-1 and CFDKW-2, respectively. It is seen from the Figure 7 that the CFDKW algorithms, Hu's algorithm, and T-FDK algorithm give similar performance in terms of axial intensity drop. This motivates us to develop an efficient approach using the estimate of missing term as given by Zhu et al. [20].

Efficient Approaches to Reduce the Axial Intensity Drop in Circular CB CT
It has been shown by reformulating Grangeat's algorithm for circular CB trajectory that the reconstruction of the true object f can be written as the sum of three terms [19]: where f FDK is the reconstructed object using FDK algorithm, f H is Hu's correction term, and f N is the missing term. Hu's term represents the information contained in a circular CB scan but not utilized in the FDK reconstruction. First two terms of the equation are obtained using the circular CB data [19]. The third term can be obtained in different ways. One approach is to estimate the missing term from the available circular CB scan data as proposed by Zhu et al. [20], which includes all the three terms in (30).
In this section, we propose to use the circular FDK inversion with no backprojection weight (CFDKW algorithms) and the estimate of missing term f N using Zhu et al.'s to improve the axial intensity profile of the reconstructed image. Thus, the reconstructed object f of the true object f , with the two algorithms, CFDKW-1 and CFDKW-2, respectively, can be written as: These two efficient algorithms are called EM-1 and EM-2, respectively, for further reference. The estimation of missing term is based on Hu's theory [19] and Grangeat's theory [24]. The formula for estimating the missing term [20] is given by where g(β, z) = (∂/∂z 2 ) ∞ −∞ (R 0 / R 2 0 + u 2 + z 2 )g(β, u, z)du. We show through simulations that the efficient algorithms and Zhu et al.'s [20] algorithm give qualitatively same intensity profile and improvement in intensity drop.

Simulation Results
Cone-beam circular and helical X-ray source trajectory scanner is used in the computer simulations. The radius of the circular and helical trajectory of the source R 0 is 2.4 m and the object radius is 1 m for the simulations. Circular CB projection data is simulated with 450 uniformly spaced source positions over 2π interval, using 283 detector elements, and 283 detector rows. Helical CB projection data is simulated with 450 uniformly spaced source positions over 2π interval, using 283 detector elements, and 29 (111) detector rows with pitch value P = 0.125 m (0.5 m). Image matrix for a slice of the 3D phantom is 256 × 256.
Simulations are carried out on full, short, and very short scan data obtained from 3D Shepp-Logan low contrast phantom [14] for HFDKW algorithms. HFDKW-2 algorithm is based on fan-beam backprojection over 2π scan and therefore it can only handle full-scan reconstruction. HFDKW-1 is a more general algorithm which can accommodate both short scan as well as very short scan data. Figure 5(a) shows the reconstructed slice at z = −0.25 of Shepp-Logan low-contrast phantom using HFDKW-1 for short scan (scan angle of π + fan angle) and Figure 5(b) shows the reconstruction with very short scan case (scan angle = π). The 3D Shepp-Logan low-contrast phantom, ForBild head phantom and Defrise disk phantom, are used to illustrate the performance of the proposed algorithms for 3D reconstruction with circular scan trajectories under different conditions. Defrise disk phantom consists of nine ellipsoidal discs stacked along the z direction. Each disc has a uniform attenuation coefficient of 0.7, and the ellipsoid has a diameter of 1.2 m and a thickness of 0.12 m, with a distance of 0.2 m between discs. We evaluate the performance of the CFDKW algorithms in terms of longitudinal or axial intensity drop. We compare the CFDKW algorithms with three different algorithms: the FDK algorithm which is only the first term of (30), Hu's algorithm (consists of first two terms of (30)) [19], and T-FDK algorithm which has been developed heuristically with a structure of shift-invariant filtering [25]. We also illustrate the performance of the EM-1 and EM-2 in terms of reduction in the axial intensity drop, noise level, and computational complexity. We compare the EM-1 and EM-2 with Zhu et al.'s [20] algorithm (i.e, all the three terms in (30)). of 3D Shepp-Logan low-contrast phantom with helical pitch 0.5 m using FDK, HFDKW-1, and HFDKW-2 algorithms. It is observed from Figures 6(a), 6(b), and 6(c) that the reconstructed image using the helical FDK shows the axial intensity drop and variation compared to the reconstructed images using HFDKW algorithms. Thus, HFDKW algorithms can be applicable for large value of pitch. We study the axial intensity drop associated with the circular FDK algorithms due to data incompleteness of the circular trajectory in 3D reconstruction. The Shepp-Logan low-contrast phantom, ForBild phantom, and the Defrise disk phantom are used to illustrate the performance of the CFDKW algorithms, in comparison with FDK, Hu's, and T-FDK algorithms. The first column of the Figure 7 shows the reconstruction of the Shepp-Logan low-contrast phantom with full scan circular trajectory using FDK, Hu's, T-FDK and CFDKW algorithms. It is observed that the drop in density in the case of Hu's, T-FDK, CFDKW-1, and CFDKW-2 in Figure 7 is same and significantly less than the circular FDK. This can be clearly seen in the intensity profiles in the Figure 8. The second column of the Figure 7 shows the reconstruction of ForBild phantom which has high-density object in the central region leading to artifacts. It is seen in Figure 7 that the reconstructed images of ForBild head phantom have white and black streaks using the FDK and Hu's algorithms. These streaks are significantly reduced in the reconstructed images using the CFDKW and T-FDK algorithms. The Defrise disk phantom is used to show the performance of the algorithms to reconstruct sharp change in the axial intensity. The strong high-frequency components in the axial direction due to sharp change in density represents the most challenging aspect of the reconstruction using circular CB data. Figure 9 shows the 1D axial intensity profile of Defrise disk phantom illustrating the observations made earlier in Figure 7. It is seen from the axial intensity profile in Figure 9 that only the central disk is reconstructed properly by all the algorithms.

Axial Intensity
We reconstruct the three phantoms using the EM algorithms and compare them with the reconstructions obtained by Zhu et al.'s algorithm. Figure 10 shows the reconstruction using Zhu et al.'s algorithm, EM-1, and EM-2 algorithms of Shepp-Logan phantom in the first column, ForBild phantom in the second column and Defrise disk in the third column, respectively. It is seen qualitatively from the Figure 10 Figure 7. This is also evident quantitatively from the axial intensity profiles shown in Figures 11 and 12 of the Shepp-Logan and Defrise disk phantoms in Figure 10, respectively.

Noise Performance.
We simulate noisy projection data using the Shepp-Logan low-contrast phantom with Poisson noise based on an emission of N 0 = 3 × 10 5 photons. The Table 1: Variance of the noisy images in Figure 13 based on the noise-free images in Figure 7.  N/N 0 ). With this as the mean to a Poisson random process, N for the case of noisy projection is computed, corresponding noisy projection is then given by − ln(N/N 0 ). The image is reconstructed using the eight different algorithms to make a qualitative and quantitative evaluation of the CFDKW and the EM algorithms in comparison to other FDK algorithms as given in Table 1. Figure 13 shows the reconstructed noisy images for the eight algorithms in the same order as given in Table 1. To compute the error variances of the images, the error images are obtained by taking the difference of reconstructed noisy images ( Figure 13) and corresponding noise-free images as shown in Figures 7 and 10. The variances of error images are measured and the variances are given in the Table 1. It is seen from Table 1 that the noise level is reduced in the CFDKW and EM algorithms. This is due to the additional numerical differentiation step, which is usually implemented as two-or three-point difference, and these algorithms are based on Hilbert filtering which avoids the approximation introduced by the nonuniform cutoff frequency required in the ramp filter-based FBP algorithm. Since CFDKW-1 and EM-1 algorithms perform additional interpolation in the source direction in the backprojection, these algorithms give better noise performance among all algorithms.

Conclusion
We have given two fan-beam algorithms with no backprojection weight for equispace linear array detector, which are counterparts of the algorithms with no backprojection weight for equi-angular detector geometry. Two new variants of FDK algorithms for circular and helical CB scan geometry with planar detector geometry are proposed. The CFDKW-1 and HFDKW-1 are applicable to full, short and very short scan data. The CFDKW-2 and HFDKW-2 deal only with the full scan data. The CFDKW and HFDKW algorithms are based on Hilbert filtering and have no position-dependent weights in the backprojection step and hence reduction in the axial intensity drop and variation when compared with the classical FDK algorithm. The CFDKW, T-FDK, and Hu's algorithms give identical performance in the reduction of axial intensity drop in the circular CB CT. CFDKW-2 algorithm is more efficient than the FDK and Hu's algorithm in terms of computational complexity. We have given efficient algorithms to reduce the axial intensity drop in the circular scan CB CT. These efficient algorithms and Zhu et al.'s algorithm give a similar performance in terms of the axial intensity drop. The EM-2 algorithm outperforms the other algorithms in terms of computational complexity. Noise performance is better in case of CFDKW and EM algorithms than conventional FDK, Hu's algorithm, T-FDK and Zhu et al.'s algorithms.