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

Authors: A. V. Narasimhadhan; Kasi Rajgopal
Journal: International Journal of Biomedical Imaging
 (2012)

---

## Abstract

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.


---

## Body

## 1. 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–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–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.


## 2. 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].


### 2.1. Parallel Beam Projections

Letf(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(1)p(ξ⃗,s)=∫R2f(x⃗)δ(x⃗⋅ξ⃗-s)dx⃗,
where s∈ℝ, ξ⃗=(-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.Figure 1
Parallel-beam geometry.The inversion off 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,(2)pR(ξ⃗,s)=(hR*p)(ξ⃗,s)=∫Rp(ξ⃗,s′)hR(s-s′)ds′,
where hR=∫-∞∞|σ|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, (3)pR(ξ⃗,s)=(hH*∂∂sp)(ξ⃗,s)=12π∂∂spH(ξ⃗,s),
where hH(s)=-1/2π∫-∞∞jsgn(σ)exp(j2πσs)dσ. The second step in the inversion is the backprojection of filtered projections resulting from the application of Hilbert filter. Mathematically,
(4)f(x⃗)=14π∫02π∂∂spH(ξ⃗,s)|s=x⃗⋅ξ⃗dθ.


### 2.2. 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 X-ray source along a circular trajectory as shown in Figure2. The source trajectory can be parameterized by an angular parameter β and is given by a⃗(β)=(R0cosβ,R0sinβ) where R0 is the radius of the circle. The fan-beam projections can be represented by
(5)g(β,α⃗)=∫0∞f(a⃗(β)+tα⃗)dt,α⃗∈S1,
where S1 is the set of all possible unit vectors in 2D space.Figure 2
Relation between the fan-beam and parallel-beam geometry.


### 2.3. Relation between the Parallel Beam Projections and Fan-Beam Projections

Recalling thatpH denotes the Hilbert transform of the 2D Radon transform p of f. The line M corresponds to parallel beam ray and fan-beam ray as shown in Figure 2. Mathematical representation of this line satisfies s=a⃗(β)·ξ⃗, where a⃗(β) is a point on this line. The value of pH corresponding to this line can be obtained from gH(β,ξ⃗), and the relation between pH and gH(β,ξ⃗) is given by
(6)pH(ξ⃗,s)|s=a⃗(β)⋅ξ⃗=gH(β,ξ⃗),
where
(7)gH(β,ξ⃗)=-∫S1hH(ξ⃗⋅α⃗)g(β,α⃗)dα⃗.
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.Figure 3
Fan-beam equispaced linear array detector scanning geometry.


### 2.4. First Algorithm with No Backprojection Weight for Equispace Linear Array Detector

Letg(β,u) be the fan-beam projections for equispace detector and p(θ,s) be the parallel-beam projections. Using change of variable γ=tan-1(u/R0) in (7), we define Hilbert transform of g(β,u)(8)gH(β,u′)=1πR02+u′2∫-umumg(β,u)u′-u1R02+u2du.
From the relation of Hamaker (6), θ=β-tan-1(u′/R0), and s=R0u′/R02+u′2, we have gH(β,u′)=pH(β-tan-1(u′/R0),R0u′/(R02+u′2)). The derivative of gH(β,u′) with respect to u′ and β using chain rule is given by
(9)∂∂u′gH(β,u′)=[∂s∂u′∂∂s+∂θ∂u′∂∂θ]pH(θ,s),∂∂βgH(β,u′)=[∂θ∂β∂∂θ]pH(θ,s).
Since the parallel beam parameters are functions of the fan-beam parameters,∂s/∂u′, ∂θ/∂u′, and ∂θ/∂β are given by
(10)∂s∂u′=R03R02+u′21R02+u′2,∂θ∂u′=-R0R02+u′2,∂θ∂β=1.
Substituting (10) in (9), we get(11)∂∂u′gH(β,u′)=[R03R02+u′21R02+u′2∂∂s+-R0R02+u′2∂∂θ]×pH(β-tan-1u′R0,R0u′R02+u′2),(12)∂∂βgH(β,u′)=∂∂θpH(β-tan-1u′R0,R0u′R02+u′2).Substituting (12) in (11), (11) can be written as
(13)∂∂u′gH(β,u′)=R0R02+u′2[R02R02+u′2∂∂spH(β-tan-1u′R0,R0u′R02+u′2)∂gH(β,u′)∂β∂gH(β,u′)∂β×R02R02+u′2-∂gH(β,u′)∂β].
In (4), (∂/∂s)pH is needed to reconstruct the function f, therefore (13) can be rewritten as(14)∂∂spH(β-tan-1u′R0,R0u′R02+u′2)=(∂∂u′gH(β,u′)R02+u′2R0+∂∂βgH(β,u′))R02+u′2R02.
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(15)f(x⃗)=14π∫02π[R02+u′2R02(∂∂u′gH(β,u′)R02+u′2R0+∂∂βgH(β,u′))R02+u′2R02×R02+u′2R02]dθ,
where β=θ-tan-1(u′/R0) and u′=x⃗·ξ⃗R0/R02-(x⃗·ξ⃗)2. The integration in this formula being over θ, it is a parallel-beam 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:(16)w(β,u′)=wi(1-0.5wa),
where wi is equal to 1 if β=(θ+tan-1(u′/R0))-floor((θ+tan-1(u′/R0))/2π)*2π≤π otherwise wi=0 and wa is equal to 1 if β=(θ-tan-1(u′/R0)+π)-floor((θ-tan-1(u′/R0)+π)/2π)*2π≤π otherwise wa is equal to 0. The reconstruction formula for this scanning range is given by(17)f(x⃗)=12π∫0πw(β,u′)[R02+u′2R02(∂∂u′gH(β,u′)R02+u′2R0+∂∂βgH(β,u′))R02+u′2R02]dθ.


### 2.5. 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
(18)f(x⃗)=12π∫02π1L[w(β,x⃗)gH(β,u′)]u′=u′(β,x⃗)dβ,
where
(19)gH(β,u′)=∫-umumhH(u′-u)R0R02+u2×(∂∂β+R02+u2R0∂∂u)g(β,u)du,(20)u′(x,y,β)=R0x⃗⋅e⃗uR0+x⃗⋅e⃗w,L=R0+x⃗⋅e⃗w.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,(21)w(β,x⃗)=‖x⃗-a⃗(β)‖2R0cos⁡γ′,
where γ′=arctan(u′/R0). 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
(22)w(β,x⃗)=LR02+u′22R03,
where L=∥x⃗-a⃗(β)∥cosγ′=R0+x⃗·e⃗w. The second fan-beam inversion formula with no backprojection weight for equispace linear array detector is obtained by substituting (22) in (18) and is given as
(23)f(x⃗)=12π∫02πR02+u′22R03gH(β,u′)u′=u′(β,x⃗)dβ,
where gH is given by (19) and u′=R0x⃗·e⃗u/(R0+x⃗·e⃗w).Figure 4
Illustration of a fan-beam scanning geometry showing redundantly measured line integral along lineK. The line integral along K corresponds to source positions a(β) and a(βc).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).


## 2.1. Parallel Beam Projections

Letf(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(1)p(ξ⃗,s)=∫R2f(x⃗)δ(x⃗⋅ξ⃗-s)dx⃗,
where s∈ℝ, ξ⃗=(-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.Figure 1
Parallel-beam geometry.The inversion off 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,(2)pR(ξ⃗,s)=(hR*p)(ξ⃗,s)=∫Rp(ξ⃗,s′)hR(s-s′)ds′,
where hR=∫-∞∞|σ|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, (3)pR(ξ⃗,s)=(hH*∂∂sp)(ξ⃗,s)=12π∂∂spH(ξ⃗,s),
where hH(s)=-1/2π∫-∞∞jsgn(σ)exp(j2πσs)dσ. The second step in the inversion is the backprojection of filtered projections resulting from the application of Hilbert filter. Mathematically,
(4)f(x⃗)=14π∫02π∂∂spH(ξ⃗,s)|s=x⃗⋅ξ⃗dθ.


## 2.2. 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 X-ray source along a circular trajectory as shown in Figure2. The source trajectory can be parameterized by an angular parameter β and is given by a⃗(β)=(R0cosβ,R0sinβ) where R0 is the radius of the circle. The fan-beam projections can be represented by
(5)g(β,α⃗)=∫0∞f(a⃗(β)+tα⃗)dt,α⃗∈S1,
where S1 is the set of all possible unit vectors in 2D space.Figure 2
Relation between the fan-beam and parallel-beam geometry.


## 2.3. Relation between the Parallel Beam Projections and Fan-Beam Projections

Recalling thatpH denotes the Hilbert transform of the 2D Radon transform p of f. The line M corresponds to parallel beam ray and fan-beam ray as shown in Figure 2. Mathematical representation of this line satisfies s=a⃗(β)·ξ⃗, where a⃗(β) is a point on this line. The value of pH corresponding to this line can be obtained from gH(β,ξ⃗), and the relation between pH and gH(β,ξ⃗) is given by
(6)pH(ξ⃗,s)|s=a⃗(β)⋅ξ⃗=gH(β,ξ⃗),
where
(7)gH(β,ξ⃗)=-∫S1hH(ξ⃗⋅α⃗)g(β,α⃗)dα⃗.
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.Figure 3
Fan-beam equispaced linear array detector scanning geometry.


## 2.4. First Algorithm with No Backprojection Weight for Equispace Linear Array Detector

Letg(β,u) be the fan-beam projections for equispace detector and p(θ,s) be the parallel-beam projections. Using change of variable γ=tan-1(u/R0) in (7), we define Hilbert transform of g(β,u)(8)gH(β,u′)=1πR02+u′2∫-umumg(β,u)u′-u1R02+u2du.
From the relation of Hamaker (6), θ=β-tan-1(u′/R0), and s=R0u′/R02+u′2, we have gH(β,u′)=pH(β-tan-1(u′/R0),R0u′/(R02+u′2)). The derivative of gH(β,u′) with respect to u′ and β using chain rule is given by
(9)∂∂u′gH(β,u′)=[∂s∂u′∂∂s+∂θ∂u′∂∂θ]pH(θ,s),∂∂βgH(β,u′)=[∂θ∂β∂∂θ]pH(θ,s).
Since the parallel beam parameters are functions of the fan-beam parameters,∂s/∂u′, ∂θ/∂u′, and ∂θ/∂β are given by
(10)∂s∂u′=R03R02+u′21R02+u′2,∂θ∂u′=-R0R02+u′2,∂θ∂β=1.
Substituting (10) in (9), we get(11)∂∂u′gH(β,u′)=[R03R02+u′21R02+u′2∂∂s+-R0R02+u′2∂∂θ]×pH(β-tan-1u′R0,R0u′R02+u′2),(12)∂∂βgH(β,u′)=∂∂θpH(β-tan-1u′R0,R0u′R02+u′2).Substituting (12) in (11), (11) can be written as
(13)∂∂u′gH(β,u′)=R0R02+u′2[R02R02+u′2∂∂spH(β-tan-1u′R0,R0u′R02+u′2)∂gH(β,u′)∂β∂gH(β,u′)∂β×R02R02+u′2-∂gH(β,u′)∂β].
In (4), (∂/∂s)pH is needed to reconstruct the function f, therefore (13) can be rewritten as(14)∂∂spH(β-tan-1u′R0,R0u′R02+u′2)=(∂∂u′gH(β,u′)R02+u′2R0+∂∂βgH(β,u′))R02+u′2R02.
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(15)f(x⃗)=14π∫02π[R02+u′2R02(∂∂u′gH(β,u′)R02+u′2R0+∂∂βgH(β,u′))R02+u′2R02×R02+u′2R02]dθ,
where β=θ-tan-1(u′/R0) and u′=x⃗·ξ⃗R0/R02-(x⃗·ξ⃗)2. The integration in this formula being over θ, it is a parallel-beam 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:(16)w(β,u′)=wi(1-0.5wa),
where wi is equal to 1 if β=(θ+tan-1(u′/R0))-floor((θ+tan-1(u′/R0))/2π)*2π≤π otherwise wi=0 and wa is equal to 1 if β=(θ-tan-1(u′/R0)+π)-floor((θ-tan-1(u′/R0)+π)/2π)*2π≤π otherwise wa is equal to 0. The reconstruction formula for this scanning range is given by(17)f(x⃗)=12π∫0πw(β,u′)[R02+u′2R02(∂∂u′gH(β,u′)R02+u′2R0+∂∂βgH(β,u′))R02+u′2R02]dθ.


## 2.5. 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
(18)f(x⃗)=12π∫02π1L[w(β,x⃗)gH(β,u′)]u′=u′(β,x⃗)dβ,
where
(19)gH(β,u′)=∫-umumhH(u′-u)R0R02+u2×(∂∂β+R02+u2R0∂∂u)g(β,u)du,(20)u′(x,y,β)=R0x⃗⋅e⃗uR0+x⃗⋅e⃗w,L=R0+x⃗⋅e⃗w.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,(21)w(β,x⃗)=‖x⃗-a⃗(β)‖2R0cos⁡γ′,
where γ′=arctan(u′/R0). 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
(22)w(β,x⃗)=LR02+u′22R03,
where L=∥x⃗-a⃗(β)∥cosγ′=R0+x⃗·e⃗w. The second fan-beam inversion formula with no backprojection weight for equispace linear array detector is obtained by substituting (22) in (18) and is given as
(23)f(x⃗)=12π∫02πR02+u′22R03gH(β,u′)u′=u′(β,x⃗)dβ,
where gH is given by (19) and u′=R0x⃗·e⃗u/(R0+x⃗·e⃗w).Figure 4
Illustration of a fan-beam scanning geometry showing redundantly measured line integral along lineK. The line integral along K corresponds to source positions a(β) and a(βc).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).


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


### 3.1. Cone-Beam Geometry

Letf(x,y,z) be the object density to be reconstructed and g(β,u,v) be the flat panel detector CB projection data collected with a helical scan trajectory given by a⃗(β)=[R0cosβ,R0sinβ,hβ]T where R0 is the radius of the source trajectory, h=P/2π, and P represents the helical pitch. To describe the CB data for a flat detector, we introduce the local detector coordinates (u,w,v) with unit vectors: e⃗u(β)=[-sinβ,cosβ,0]T, e⃗w(β)=[-cosβ,-sinβ,0]T, and e⃗v(β)=[0,0,1]T. Here e⃗w(β) points from the source point a⃗(β) to the center of the detector and the unit vectors e⃗u(β) and e⃗v(β) span the detector. The CB projection for flat panel detector geometry is given by(24)g(β,u,v)=∫0∞f(a⃗(β)+t(ue⃗u(β)+De⃗w(β)+ve⃗v(β)u2+D2+v2))dt,
where D is the source to detector distance.


### 3.2. 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 distanceR0 from a⃗(β) and hence contains the z-axis. Scaling all equations by substituting u→R0u/D and v→R0v/D is sufficient to describe the case of a flat panel array detector located at a distance D from the vertex.


#### 3.2.1. First FDK Formula with No Backprojection Weight for Helical Trajectory (HFDKW-1)

LethH(·) be Hilbert filter kernel in spatial domain. The HFDKW-1 formula is obtained from (15) by the standard FDK extension technique for a particular slice z,
(25)f̂(x,y,z)=14π∫(z/h)-π(z/h)+πR02+u′2R03(∫-∞∞ĝ(β,u,v′)hH(u′-u)du)dθ,
where
(26)ĝ(β,u,v)=R0R02+u2+v2(∂∂β+R02+u2R0∂∂u+uvR0∂∂v)g(β,u,v),β(x,y,θ)=θ+γ′(x,y,θ), and γ′(x,y,θ)=arcsin((-xsinθ+ycosθ)/R0). The variables u′ and v′ are computed by the following equations:(27)u′(x,y,θ)=R0(-xsin⁡β(x,y,θ)+ycos⁡β(x,y,θ))R0-xcos⁡β-ysin⁡β(x,y,θ),v′(x,y,θ)=R0(z-hβ(x,y,θ))R0-xcos⁡β(x,y,θ)-ysin⁡β(x,y,θ).


#### 3.2.2. 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(28)f̂(x,y,z)=14π∫(z/h)-π(z/h)+πR02+u′2R03∫-∞∞ĝ(β,u,v′)hH(u′-u)dudβ,
where(29)u′(x,y,β)=R0(-xsin⁡β+ycos⁡β)R0-xcos⁡β-ysin⁡β,v′(x,y,z,β)=R0(z-hβ)R0-xcos⁡β-ysin⁡β,
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 weightL-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 (R02+u2)/R03 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.


#### 3.2.3. 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 pitchP=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].


## 3.1. Cone-Beam Geometry

Letf(x,y,z) be the object density to be reconstructed and g(β,u,v) be the flat panel detector CB projection data collected with a helical scan trajectory given by a⃗(β)=[R0cosβ,R0sinβ,hβ]T where R0 is the radius of the source trajectory, h=P/2π, and P represents the helical pitch. To describe the CB data for a flat detector, we introduce the local detector coordinates (u,w,v) with unit vectors: e⃗u(β)=[-sinβ,cosβ,0]T, e⃗w(β)=[-cosβ,-sinβ,0]T, and e⃗v(β)=[0,0,1]T. Here e⃗w(β) points from the source point a⃗(β) to the center of the detector and the unit vectors e⃗u(β) and e⃗v(β) span the detector. The CB projection for flat panel detector geometry is given by(24)g(β,u,v)=∫0∞f(a⃗(β)+t(ue⃗u(β)+De⃗w(β)+ve⃗v(β)u2+D2+v2))dt,
where D is the source to detector distance.


## 3.2. 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 distanceR0 from a⃗(β) and hence contains the z-axis. Scaling all equations by substituting u→R0u/D and v→R0v/D is sufficient to describe the case of a flat panel array detector located at a distance D from the vertex.


### 3.2.1. First FDK Formula with No Backprojection Weight for Helical Trajectory (HFDKW-1)

LethH(·) be Hilbert filter kernel in spatial domain. The HFDKW-1 formula is obtained from (15) by the standard FDK extension technique for a particular slice z,
(25)f̂(x,y,z)=14π∫(z/h)-π(z/h)+πR02+u′2R03(∫-∞∞ĝ(β,u,v′)hH(u′-u)du)dθ,
where
(26)ĝ(β,u,v)=R0R02+u2+v2(∂∂β+R02+u2R0∂∂u+uvR0∂∂v)g(β,u,v),β(x,y,θ)=θ+γ′(x,y,θ), and γ′(x,y,θ)=arcsin((-xsinθ+ycosθ)/R0). The variables u′ and v′ are computed by the following equations:(27)u′(x,y,θ)=R0(-xsin⁡β(x,y,θ)+ycos⁡β(x,y,θ))R0-xcos⁡β-ysin⁡β(x,y,θ),v′(x,y,θ)=R0(z-hβ(x,y,θ))R0-xcos⁡β(x,y,θ)-ysin⁡β(x,y,θ).


### 3.2.2. 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(28)f̂(x,y,z)=14π∫(z/h)-π(z/h)+πR02+u′2R03∫-∞∞ĝ(β,u,v′)hH(u′-u)dudβ,
where(29)u′(x,y,β)=R0(-xsin⁡β+ycos⁡β)R0-xcos⁡β-ysin⁡β,v′(x,y,z,β)=R0(z-hβ)R0-xcos⁡β-ysin⁡β,
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 weightL-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 (R02+u2)/R03 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.


### 3.2.3. 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 pitchP=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].


## 3.2.1. First FDK Formula with No Backprojection Weight for Helical Trajectory (HFDKW-1)

LethH(·) be Hilbert filter kernel in spatial domain. The HFDKW-1 formula is obtained from (15) by the standard FDK extension technique for a particular slice z,
(25)f̂(x,y,z)=14π∫(z/h)-π(z/h)+πR02+u′2R03(∫-∞∞ĝ(β,u,v′)hH(u′-u)du)dθ,
where
(26)ĝ(β,u,v)=R0R02+u2+v2(∂∂β+R02+u2R0∂∂u+uvR0∂∂v)g(β,u,v),β(x,y,θ)=θ+γ′(x,y,θ), and γ′(x,y,θ)=arcsin((-xsinθ+ycosθ)/R0). The variables u′ and v′ are computed by the following equations:(27)u′(x,y,θ)=R0(-xsin⁡β(x,y,θ)+ycos⁡β(x,y,θ))R0-xcos⁡β-ysin⁡β(x,y,θ),v′(x,y,θ)=R0(z-hβ(x,y,θ))R0-xcos⁡β(x,y,θ)-ysin⁡β(x,y,θ).


## 3.2.2. 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(28)f̂(x,y,z)=14π∫(z/h)-π(z/h)+πR02+u′2R03∫-∞∞ĝ(β,u,v′)hH(u′-u)dudβ,
where(29)u′(x,y,β)=R0(-xsin⁡β+ycos⁡β)R0-xcos⁡β-ysin⁡β,v′(x,y,z,β)=R0(z-hβ)R0-xcos⁡β-ysin⁡β,
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 weightL-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 (R02+u2)/R03 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.


## 3.2.3. 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 pitchP=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].


## 4. 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 objectf can be written as the sum of three terms [19]:(30)f̂=f̂FDK+f̂H+f̂N,
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 termf̂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:(31)f̂=f̂CFDKW-1+f̂N,f̂=f̂CFDKW-2+f̂N.
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(32)f̂N(x⃗)=-14π2R02+z2R02(1-R02-z2R0)∫02πĝ(β,z)dβ,
where ĝ(β,z)=(∂/∂z2)∫-∞∞(R0/R02+u2+z2)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.


## 5. 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 sourceR0 is 2.4m and the object radius is 1m 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.125m (0.5m). 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 = π).Reconstructed slices atz=-0.25 of Shepp-Logan low-contrast phantom using HFDKW-1 formula for (a) short scan segment on a helical trajectory. (b) Less than short scan segment. The value of pitch used is 0.125m. Display window is [11.04].
(a)(b)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 thez direction. Each disc has a uniform attenuation coefficient of 0.7, and the ellipsoid has a diameter of 1.2m and a thickness of 0.12m, with a distance of 0.2m 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)).


### 5.1. Axial Intensity Drop

The axial intensity drop and variation become significant in conventional helical FDK algorithms due to the approximate nature of the algorithm when large helical pitch is used. We reconstructed the slices of 3D Shepp-Logan low-contrast phantom with helical pitch0.5m 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.Slices of Shepp-Logan low-contrast phantom atx=0. Illustration of axial intensity drop due to large value of pitch 0.5m in full-scan reconstruction (a) standard helical FDK. (b) HFDKW-1 algorithm and (c) HFDKW-2 algorithm. Display window is [11.04].
(a)(b)(c)From left to right column, reconstructed slices of Shepp-Logan low contrast, ForBild head and Defrise disk phantom, respectively, atx=0. From top row to bottom row, reconstruction with FDK algorithm, Hu’s algorithm, T-FDK algorithm, CFDKW-1, and CFDKW-2 respectively. Display windows are [0.981.05],[1.011.09] and [0.50.7] for Shepp-Logan, ForBild head and Defrise disk phantom, respectively.
(a)(b)(c)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 Figure7 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.Figure 8
Illustration of axial intensity drop in full-scan reconstructions of Shepp-Logan low-contrast phantom, vertical intensity profile of the Shepp-Logan low-contrast phantom atz=0 in Figure 7.Figure 9
Illustration of axial intensity drop in full-scan reconstructions of Defrise disk phantom, vertical intensity profile of the Defrise disk phantom atz=0 in Figure 7.We reconstruct the three phantoms using the EM algorithms and compare them with the reconstructions obtained by Zhu et al.’s algorithm. Figure10 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 that Zhu et al.’s, EM-1, and EM-2 reconstruction algorithms give the same performance in terms of the axial intensity drop and significant improvement over reconstructions shown in 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.From left to right column, reconstructed slices of the Shepp-Logan low contrast, ForBild head and Defrise disk phantom, respectively, atx=0. From top row to bottom row, reconstruction with Zhu et al.’s algorithm, EM-1 (CFDKW-1 + FN), and EM-2 (CFDKW-2 + FN), respectively. Display windows are [0.981.05],[1.011.09] and [0.50.7] for Shepp-Logan, ForBild head, and Defrise disk phantom, respectively.
(a)(b)(c)Figure 11
Illustration of axial intensity drop in full-scan reconstructions of Shepp-Logan low contrast phantom, vertical profile of the Shepp-Logan low contrast phantom atz=0 in Figure 10.Figure 12
Illustration of axial intensity drop in full-scan reconstructions of Defrise disk phantom, vertical profile of the Defrise disk phantom atz=0 in Figure 10.


### 5.2. Noise Performance

We simulate noisy projection data using the Shepp-Logan low-contrast phantom with Poisson noise based on an emission ofN0=3×105 photons. The number of photons detected N was calculated by -ln(N/N0). 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/N0). 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.Table 1
Variance of the noisy images in Figure13 based on the noise-free images in Figure 7.
Algorithmvalue×10-4FDK1.3986Hu’s1.3999T-FDK0.27443CFDKW-10.24365CFDKW-20.26443Zhu et al.’s1.3999EM-10.24365EM-20.26443Slices of Shepp-Logan low contrast phantom atx=0, reconstruction with (a) FDK algorithm, Hu’s algorithm, T-FDK algorithm, (b) CFDKW-1, CFDKW-2, Zhu et al.’s algorithm, (c) EM-1 and EM-2 with added noise. Display window is [0.981.05].
(a)(b)(c)


### 5.3. Computational Complexity

We compute the execution time for FDK algorithm, Hu’s term, missing term, Hu’s algorithm, Zhu et al.’s algorithm, CFDKW-2, and EM-2 algorithms in MATLAB 7.0 on a Intel P4 desktop computer with a dual core2.4GHz processor platform. The algorithm modules considered for estimating the execution time are (1) 1D filtering of the CB projections and (2) backprojection of a single slice. The execution time of the reconstruction for FDK, Hu’s term, missing term, Zhu et al.’s algorithm, CFDK-2, and EM-2 algorithms are computed using these modules. The execution times are given in the Table 2. CFDKW-2 does not have position-dependent weight in the backprojection, hence CFDKW-2 algorithm has less computational complexity than FDK and Hu’s algorithms. The computational complexity of EM-2 algorithm is less than Zhu et al.’s algorithm. The EM-2 algorithm consists of the reconstruction of the image using CFDKW-2 algorithm and estimating the missing term, whereas Zhu et al.’s algorithm consists of the reconstruction of the image using conventional FDK, estimation of Hu’s term, and estimation of missing term. Hence the computational complexity of EM-2 is less than Zhu et al.’s algorithm.Table 2
Execution time (in seconds).
Single sliceTotal volume (256 × 256 × 256)BackprojectionFilteringBackprojectionTotalFDK11.876.430203096.4Hu’s term5.62.31433.61435.9Missing term5.32.41356.81359.2Hu’s algorithm17.478.74453.64542.3Zhu et al.’s algorithm22.781.15811.25892.3CFDKW-2119028162906EM-216.392.44172.84265.2


## 5.1. Axial Intensity Drop

The axial intensity drop and variation become significant in conventional helical FDK algorithms due to the approximate nature of the algorithm when large helical pitch is used. We reconstructed the slices of 3D Shepp-Logan low-contrast phantom with helical pitch0.5m 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.Slices of Shepp-Logan low-contrast phantom atx=0. Illustration of axial intensity drop due to large value of pitch 0.5m in full-scan reconstruction (a) standard helical FDK. (b) HFDKW-1 algorithm and (c) HFDKW-2 algorithm. Display window is [11.04].
(a)(b)(c)From left to right column, reconstructed slices of Shepp-Logan low contrast, ForBild head and Defrise disk phantom, respectively, atx=0. From top row to bottom row, reconstruction with FDK algorithm, Hu’s algorithm, T-FDK algorithm, CFDKW-1, and CFDKW-2 respectively. Display windows are [0.981.05],[1.011.09] and [0.50.7] for Shepp-Logan, ForBild head and Defrise disk phantom, respectively.
(a)(b)(c)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 Figure7 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.Figure 8
Illustration of axial intensity drop in full-scan reconstructions of Shepp-Logan low-contrast phantom, vertical intensity profile of the Shepp-Logan low-contrast phantom atz=0 in Figure 7.Figure 9
Illustration of axial intensity drop in full-scan reconstructions of Defrise disk phantom, vertical intensity profile of the Defrise disk phantom atz=0 in Figure 7.We reconstruct the three phantoms using the EM algorithms and compare them with the reconstructions obtained by Zhu et al.’s algorithm. Figure10 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 that Zhu et al.’s, EM-1, and EM-2 reconstruction algorithms give the same performance in terms of the axial intensity drop and significant improvement over reconstructions shown in 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.From left to right column, reconstructed slices of the Shepp-Logan low contrast, ForBild head and Defrise disk phantom, respectively, atx=0. From top row to bottom row, reconstruction with Zhu et al.’s algorithm, EM-1 (CFDKW-1 + FN), and EM-2 (CFDKW-2 + FN), respectively. Display windows are [0.981.05],[1.011.09] and [0.50.7] for Shepp-Logan, ForBild head, and Defrise disk phantom, respectively.
(a)(b)(c)Figure 11
Illustration of axial intensity drop in full-scan reconstructions of Shepp-Logan low contrast phantom, vertical profile of the Shepp-Logan low contrast phantom atz=0 in Figure 10.Figure 12
Illustration of axial intensity drop in full-scan reconstructions of Defrise disk phantom, vertical profile of the Defrise disk phantom atz=0 in Figure 10.


## 5.2. Noise Performance

We simulate noisy projection data using the Shepp-Logan low-contrast phantom with Poisson noise based on an emission ofN0=3×105 photons. The number of photons detected N was calculated by -ln(N/N0). 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/N0). 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.Table 1
Variance of the noisy images in Figure13 based on the noise-free images in Figure 7.
Algorithmvalue×10-4FDK1.3986Hu’s1.3999T-FDK0.27443CFDKW-10.24365CFDKW-20.26443Zhu et al.’s1.3999EM-10.24365EM-20.26443Slices of Shepp-Logan low contrast phantom atx=0, reconstruction with (a) FDK algorithm, Hu’s algorithm, T-FDK algorithm, (b) CFDKW-1, CFDKW-2, Zhu et al.’s algorithm, (c) EM-1 and EM-2 with added noise. Display window is [0.981.05].
(a)(b)(c)


## 5.3. Computational Complexity

We compute the execution time for FDK algorithm, Hu’s term, missing term, Hu’s algorithm, Zhu et al.’s algorithm, CFDKW-2, and EM-2 algorithms in MATLAB 7.0 on a Intel P4 desktop computer with a dual core2.4GHz processor platform. The algorithm modules considered for estimating the execution time are (1) 1D filtering of the CB projections and (2) backprojection of a single slice. The execution time of the reconstruction for FDK, Hu’s term, missing term, Zhu et al.’s algorithm, CFDK-2, and EM-2 algorithms are computed using these modules. The execution times are given in the Table 2. CFDKW-2 does not have position-dependent weight in the backprojection, hence CFDKW-2 algorithm has less computational complexity than FDK and Hu’s algorithms. The computational complexity of EM-2 algorithm is less than Zhu et al.’s algorithm. The EM-2 algorithm consists of the reconstruction of the image using CFDKW-2 algorithm and estimating the missing term, whereas Zhu et al.’s algorithm consists of the reconstruction of the image using conventional FDK, estimation of Hu’s term, and estimation of missing term. Hence the computational complexity of EM-2 is less than Zhu et al.’s algorithm.Table 2
Execution time (in seconds).
Single sliceTotal volume (256 × 256 × 256)BackprojectionFilteringBackprojectionTotalFDK11.876.430203096.4Hu’s term5.62.31433.61435.9Missing term5.32.41356.81359.2Hu’s algorithm17.478.74453.64542.3Zhu et al.’s algorithm22.781.15811.25892.3CFDKW-2119028162906EM-216.392.44172.84265.2


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



---

*Source: 969432-2012-02-16.xml*