PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Tsinghua Sci Technol. Author manuscript; available in PMC 2010 June 21.
Published in final edited form as:
Tsinghua Sci Technol. 2010 February; 15(1): 17–24.
doi:  10.1016/S1007-0214(10)70003-1
PMCID: PMC2888110
NIHMSID: NIHMS200124

Filtered Backprojection Reconstruction with Depth-Dependent Filtering*

Abstract

A direct filtered-backprojection (FBP) reconstruction algorithm is presented for circular cone-beam computed tomography (CB-CT) that allows the filter operation to be applied efficiently with shift-variant band-pass characteristics on the kernel function. Our algorithm is derived from the ramp-filter based FBP method of Feldkamp et al. and obtained by decomposing the ramp filtering into a convolution involving the Hilbert kernel (global operation) and a subsequent differentiation operation (local operation). The differentiation is implemented as a finite difference of two (Hilbert filtered) data samples and carried out as part of the backprojection step. The spacing between the two samples, which defines the low-pass characteristics of the filter operation, can thus be selected individually for each point in the image volume. We here define the sample spacing to follow the magnification of the divergent-beam geometry and thus obtain a novel, depth-dependent filtering algorithm for circular CB-CT. We evaluate this resulting algorithm using computer-simulated CB data and demonstrate that our algorithm yields results where spatial resolution and image noise are distributed much more uniformly over the field-of-view, compared to Feldkamp’s approach.

Keywords: computed tomography (CT), reconstruction, filtered backprojection, circular cone-beam CT

Introduction

Most X-ray scanners for volumetric computed tomography (CT) that are currently in clinical use apply a filtered backprojection (FBP) algorithm, such as in Refs. [1-4], to convert a series of acquired X-ray conebeam (CB) projections into the 3-D density distribution of the investigated object. In a conventional FBP algorithm, CB projections are first weighted and convolved with a non-local 1-D high-pass filter (very often the ramp filter) along a set of lines to create filtered CB data. Subsequently, each sample of the filtered data is transferred into the image domain by a weighted backprojection, where it additively contributes to the reconstruction at all points that are located on the ray connecting the sample position with the X-ray source point.

FBP approaches that operate directly in the fan-beam or CB geometry are very attractive in terms of data flow and efficiency, but they yield reconstruction results where image noise is nonstationary[5,6] and where spatial resolution is locally varying over the field-of-view. A fundamental reason for this phenomena is that these FBP methods neglect the fact that during data filtering, the kernel band-limitation would need to vary with the depth of the voxel of interest along the backprojection ray[6]. For efficiency reasons, however, conventional FBP methods use a kernel that is independent of this depth, as explained in the previous paragraph. They thus reconstruct slices of spatially-varying image characteristics.

In the last few years, several direct FBP reconstruction methods have been suggested that aim at achieving more homogenous image quality. Some of these methods involve special backprojection operators[7-10], while others aim at beneficially making use of the redundancies in the acquired data set[11-13].

In this paper, we propose a new direct FBP reconstruction algorithm that improves the uniformity of image noise and spatial resolution. Our approach involves depth-dependent filter kernels, where the band-pass characteristics of the filter may depend on the distance between the detector and the point to be reconstructed. The proposed algorithm is efficient, and can be used for both, full-scan and short-scan geometries.

1 Geometry and Notation

Assume that CB data is acquired in the geometry illustrated in Fig. 1. The X-ray linear attenuation coefficient of the object under investigation is described using the function f (x) with x = (x, y, z). The object is entirely contained in a cylindrical region with circular base of radius R0 that is centered on the z-axis. During the scan, the focal spot of the X-ray source describes the circular path a(λ) = (Rcosλ, Rsinλ, 0), where λ describes the angular position of the source, and R corresponds to the scan radius. Data is acquired with a 2-D flat detector that has distance D from the source and that is parallel to the vectors eu(λ)=(−sinλ, cosλ, 0) and ev(λ) = (0, 0, 1). The vector ew(λ)=(cosλ, sinλ, 0) is perpendicular to the detector and faces the source. CB data can then be expressed as

g(λ,u,v)=0f(a(λ)+tα(λ,u,v))dt
(1)

where u and v are the detector coordinates measured along eu(λ) and ev(λ), respectively, and where

α(λ,u,v)=ueu+vevDew(λ)u2+v2+D2
(2)

gives the unit direction of the ray that connects the source with the detector point (u, v).

Fig. 1
Top-view of the circular CB acquisition geometry

We focus on the problem of recovering f (x) from CB data g(λ, u, v) that is known over at least a short-scan, i.e., for λ[sm epsilon][λ1, λ2] with λ2 - λ1 >= π + 2arcsin(R0/R). Moreover, we assume that no truncation occurs in u, i.e., g(λ,u, v) is known for u [sm epsilon] [−um,um] with um=R0DR2R02.

2 Feldkamp FBP Approach

Following the notation of Section 1, the Feldkamp formula for CB reconstruction[1] can be expressed as

fe(x)=λ1λ2RD(Rxe(λ))2gF(λ,u,v)dλ
(3)

where u* and v* are the coordinates of the CB projection of x onto the detector, and where the subscript e emphasizes that the outcome is an estimate of the true object density function f. The function gF denotes weighted and filtered data,

gF(λ,u,v)=umumhr(uu)Dw(λ,u)g(λ,u,v)u2+v2+D2du
(4)

In this equation, w(λ,u) is a smooth weighting function that approximately accounts for redundancies in the data[14] and hr (u) is the ramp filter kernel[15]. Note that reconstruction from sampled data requires regularization of the reconstruction process[16]. Typically, this is achieved by replacing the ramp filter kernel in Eq. (4) by its band-limited version,

hrb(u)=Φ(uu)hr(u)du
(5)

The function [var phi] (u) describes a low-pass filter that can be used to globally balance noise and spatial resolution in the reconstruction results.

3 Depth-Dependent Filtering Approach

We now derive a new method for direct CB FBP reconstruction. One important feature of this method is that points that are on the same backprojection ray may receive different filter values, dependent on their depth along the ray. In this feature, our method differs from conventional FBP algorithms. The derivation is based on decomposing the ramp filter operation into the combination of a local component and a remaining global component by making use of the identity

hr(u)=12πuhh(u)
(6)

where hh (u) denotes the kernel of the Hilbert transform as shown in Ref. [17]. This means that data filtering according to Eq. (4) can alternatively be achieved by first applying the Hilbert transform on the weighted data to obtain an intermediate function

gH(λ,u,v)=umumhh(uu)Dw(λ,u)g(λ,u,v)u2+v2+D2du
(7)

and by subsequently computing

gF(λ,u,v)=12πδ(uu)gH(λ,u,v)du
(8)

where δ’ is the kernel of the differentiation operator. The benefit of this decomposition is that the differentiation can be computed very efficiently compared to the convolution with the non-local ramp kernel in Eq. (4). That means that the differentiation component of the filtering process can be carried out individually for each voxel without significantly increasing the computational effort. In order to obtain voxel-dependency in Eq. (8), we apply the change of variable ul (at fixed x and fixed λ) that is defined by

u(l)=DlRxew(λ)
(9)

and that has Jacobian |D (R - x·ew(λ))|. Note that u and l are related to each other via a scaling factor that depends on the distance between the point x and the detector. Picture l as the coordinate that describes locations on a virtual detector, i.e., a detector that is parallel to the physical detector but shifted along ew(λ) such that it contains the point x. For a geometric illustration as shown in Fig. 2. The right hand side of Eq. (8) can then be developed into

12πDRxew(λ)δ(D(ll)Rxew(λ))gH(λ,u(l),v)dl=12πRxew(λ)Dδ(ll)gH(λ,u(l),v)dl=Rxew(λ)2πDlimε0gH(λ,u(l+ε),v)gH(λ,u(lε),v)2ε
(10)
Fig. 2
Top-view of the physical detector as well as of the virtual detector through x and of some of the associated quantities

For numerical implementation, we suggest approximating the limit in Eq. (10) by a finite difference, defined by a parameter Δl that affects the distance between the required samples of gH, namely by

gH(λ,u(l+Δl),v)gH(λ,u(lΔl),v)2Δl
(11)

The actual selection of Δl has an impact on the band-pass characteristics imposed on the differentiation operation. Combining Eqs. (11), (8), and (3) yields the new CB FBP formula for estimating the object density,

fe(x)=R4πΔlλ1λ21Rxew(λ).(gH(λ,u+u(Δl),v)gH(λ,uu(Δl),v))dλ
(12)

with gH defined in Eq. (7). Image reconstruction according to our new approach is thus achieved by filtering each CB projection according to Eq. (7) and by backprojecting the filter result into the image domain, such that each image point receives the difference of two filter values, separated by 2ul). One important feature of our algorithm is that the sample spacing 2ul) depends on the distance between x and the X-ray source so that it is closely related to the magnification of the acquisition geometry. In other words, points in the image domain may receive different filter values (computed using a different realization of the finite difference) even if they project onto the same detector coordinate u*. Consequently, our algorithm carries out a filter operation with a depth-dependent band-pass characteristics, without losing the general algorithmic structure of filtering and backprojection.

4 Numerical Evaluation

A numerical evaluation of the depth-dependent filtering (DDF) approach introduced in Section 3 was carried out and its performance was compared to that of Feldkamp’s algorithm (Section 2). We first focused on the reconstruction within the plane of the source trajectory, to avoid the issue of CB artifacts. Inside this plane, we investigated image noise and spatial resolution in x-y and focused on how uniformly these image properties are distributed over the reconstructed region. Afterwards, we briefly investigated 3-D image quality by visually assessing CB artifacts in the reconstructions of the FORBILD head phantom.

4.1 Implementation details

To allow an accurate numerical comparison, the two reconstruction formulas were implemented using identical discretization schemes, wherever possible. Backprojection, e.g., was applied with bilinear interpolation in u and v to obtain the filter value at coordinates u* and v*. The low-pass filter [var phi] used for regularizing the ramp kernel in the Feldkamp algorithm according to Eq. (5) was set to a Gaussian function with standard deviation σ. In the DDF algorithm, we applied the Hilbert transform with a half-pixel shift, as described in Ref. [18].

4.2 Image quality in the scan plane

Spatial x-y resolution and noise at z = 0 mm were both evaluated from CB data that was simulated for a full-scan (i.e., λ2 = λ1 + 2π) using the parameters listed in Table 1. A virtual detector containing the z-axis was used and the data sample at any pixel of this detector was computed using an upsampling strategy. More specifically, we first subdivided the pixel, then computed the object density integrals along the lines that connect the X-ray source point with every subdivided patch and finally averaged these integral values to obtain a sample of the function g. For the noise studies, each detector pixel along the u-axis was uniformly subdivided by a factor of 9, and for the resolution studies, a factor of 25 was used.

Table 1
Geometry parameters for the noise and resolution study

(1) Spatial resolution

We first studied spatial resolution by estimating and evaluating transaxial point-spread-functions (PSFs), using an approach similar to that of Ref. [11]. The redundancy weighting function was set to w(λ,u) = 0.5, so that spatial resolution in a transaxial slice depended only on the distance to the z-axis. Resolution at distance d was estimated using a long cylindrical object of tiny circular base (radius 0.2 mm ), which was shifted to (x, y) = (d,0). Reconstruction was carried out at z = 0 mm, on a fine Cartesian grid with Δx = Δy = 0.03 mm. The reconstruction result, which was assumed to represent the transaxial PSF, was evaluated on 360 radial profiles diverging from the PSF center. We determined the full-width at half maximum (FWHM) for each profile and used the average of these 360 FWHM values to quantify the x-y resolution at the current cylinder location.

The process described in the previous paragraph was repeated for 10 cylinder locations defined by d ={5, 15, 25, …, 95} mm, yielding resolution estimates along a radial profile within the plane of the circular scan. We then computed the mean and the standard deviation of these 10 values and used these two figures-of-merit to globally characterize spatial resolution in the reconstruction result.

Figure 3 shows these two figures-of-merit for both algorithms, and for various selections of Δl and σ. In these plots, the average resolution is indicated with a small gray circle, while the length of the errorbars is related to the spatial variation of x-y resolution. These results show that in the Feldkamp approach, resolution varies noticeably over the investigated slice, and this variation becomes even more pronounced when using smoother filter kernels. Resolution achievable with the DDF approach, in contrast, is more uniformly distributed. The advantage of the new approach becomes particularly visible for larger selections of Δl: x-y resolution then tends to become almost constant over the field-of-view, at the cost of a slight increase in resolution anisotropy far away from the isocenter as shown in Fig. 3b. Clearly, Δl and σ play a similar role in the algorithms, and it is possible to determine a pair of these parameters, for which some specific image property, such as the mean x-y resolution along the profile, is matched. The results of Fig. 3 indicate that such a matching is for instance achieved for Δl = 0.9Δu and σ = Δu; Figure 4 shows some PSFs obtained for this specific setting.

Fig. 3
Evaluation of spatial x-y resolution in the plane of the scan. Each figure displays the mean (circle) and the variance (error-bar) in the FWHM value along a radial line in the field-of-view, as a function of the smoothing parameters σ or Δ ...
Fig. 4
Illustration of 5 PSFs computed at z = 0 mm and at (from left to right) 10 mm, 20 mm, 30 mm, 40 mm and 50 mm from the z-axis. (Top) Feldkamp with σ = Δu, (bottom) the new approach for Δl = 0.9Δu. Each PSF was normalized ...

(2) Noise

For an evaluation of image noise, we simulated CB data of a cylindrical phantom of radius 100 mm that has the density of water (μ =0.01836 / mm). The phantom contains small circular inserts of different sizes and of densities that are slightly higher than that of the phantom background (the maximum contrast is 45 HU). Prior to reconstruction, we added Poisson noise to the CB data, with a noise level that corresponds to an emission of 200 000 X-ray photons per ray. We did not include any bowtie filter in our simulations. Figure 5 shows the reconstructions at z=0 mm from one noisy data set, obtained with the DDF and the Feldkamp approach. The reconstruction for the DDF approach was carried out using Δl=0.4Δu. For the Feldkamp algorithm, however, we selected two different smoothing parameters: σ = 0.35Δu, for which resolution was matched to the DDF results close to the rotation axis (at d = 5 mm), and σ = 0.5Δu, which yields a mean resolution along the profile that is identical to that of the DDF approach as shown in Fig. 3.

Fig. 5
The slice z = 0 mm through the reconstruction result obtained from one noisy full-scan fan-beam data set, in a grayscale window of [−100, 100] HU. (Left) The DDF approach with Δl = 0.4Δu, (center) the Feldkamp result with a mean ...

Figure 6 displays the noise profile along the y-axis for each of the three results, which is computed as the pixel-wise standard deviation across the reconstructions obtained from 1000 distinct noise realizations. From these results, we observe that the DDF approach can effectively avoid the strong noise amplification at the boundaries of the field-of-view that is visible in the Feldkamp results. The DDF method thus clearly improves noise uniformity over the field-of-view for the selected smoothing parameters.

Fig. 6
Noise profiles along the y-axis for the images displayed in Fig. 5, estimated using 1000 distinct noise realizations

Figure 7 displays the mean and homogeneity of image noise (computed as the mean and the standard deviation of the noise profile) achievable with both methods, for various selections of Δl and σ. Clearly, the DDF approach improves noise uniformity at identical mean noise, contributing to more uniform image characteristics over the entire FOV.

Fig. 7
Evaluation of image noise in (a) the Feldkamp approach and (b) the new approach, for various selections of the algorithmic parameters. The gray circles in each figure indicate the mean noise level along the y-axis and the error-bars indicate the uniformity ...

4.3 3-D image quality

It is well-known that CB data acquired along a planar trajectory is insufficient for numerically stable and exact reconstruction. There is no indication that the DDF method should react noticeably differently to these data insufficiencies than the classical Feldkamp algorithm. In order to support this statement, we simulated CB data of the FORBILD head phantom without ears and investigated the CB artifact in two off-center trans-axial slices. The phantom was shifted by 40 mm in positive z-direction. Data was simulated in a typical C-arm geometry, using a short-scan with λ1=0° and λ2=216°, a scan radius of R=750 mm, and again a virtual detector that contains the z-axis and has pixels with Δu = Δv = 0.5 mm.

The slices at z = 25 mm and z = 45 mm through the reconstructions obtained with the DDF and the Feldkamp method are shown in Fig. 8. These results, which are displayed in a narrow grayscale window of 100 HU width, show no differences other than some slight variations in discretization artifacts. CB artifacts are thus essentially identical.

Fig. 8
The slices (top) z = 25 mm and (bottom) z = 45 mm through the reconstructions of the head phantom, displayed in [0, 100] HU. The images on the left were obtained with Feldkamp and σ= 0.5Δu, and those on the right with the new DDF method ...

5 Conclusions

In this paper, we introduced a new algorithm for CT reconstruction from circular CB data. Our approach was derived from the classical CB algorithm of Feldkamp et al. A nice property of our approach is that it efficiently allows the band-pass characteristics of the filter operation to vary with the depth of the point to be reconstructed along the backprojection ray, without loosing the general algorithmic structure of filtering and backprojection. Numerical evaluation in the full-scan scenario and inside the plane of the scan demonstrated that this property allows for more uniformly distributed spatial x-y resolution and image noise, compared to the Feldkamp algorithm. This improvement is particularly visible when enforcing smooth image characteristics, i.e., for larger selections of Δl and σ. In terms of CB artifacts, our new depth-dependent filtering method shows no noticeable differences to the Feldkamp method.

In the last years, other direct reconstruction methods that come with a reduction in the spatial variations in image noise and resolution have been suggested[7-13]. It would be interesting to compare our method to these using figures-of-merit similar to those described above. However, this task is not necessarily easy, because for some methods, the result depends on the size of a voxel in the image domain, unlike in the method presented above. It is not clear if PSF resolution measurements then yield meaningful results, because they typically involve the reconstruction on voxel grids that are much finer than the ones used for final reconstruction.

Another popular strategy for fan-beam reconstruction is to apply a parallel-rebinning of the data and to obtain object density using a parallel-beam reconstruction formula[15]. This indirect approach comes with additional computational costs due to the rebinning step and with additional interpolations, compared to the direct methods, but it also contributes to a homogenization of image quality. A comparison between the rebinning method and our direct approach could be an exciting topic of further investigations.

Note that Δl does not need to be constant, but might for instance be defined as a function of λ and x. We expect further improvements in image quality when making use of these additional degrees-of- freedom, but still need to quantify in detail these improvements. Note finally that the idea presented in this paper is not restricted to circular fan-beam or CB-CT, but can be applied in a straight-forward manner to exact or approximate FBP reconstruction in many other geometries, involving non-planar and/or non-circular source trajectories.

Footnotes

*Supported in part by the US National Institute of Health (Nos. R01EB007236 and R21EB009168). The concepts presented in this paper are based on research and are not commercially available.

Contributor Information

Frank Dennerlein, Siemens AG, Healthcare Sector, Forchheim 91301, Germany.

Holger Kunze, Siemens AG, Healthcare Sector, Forchheim 91301, Germany.

References

[1] Feldkamp LA, Davis LC, Kress JW. Practical conebeam algorithm. J. Opt. Soc. Am. A. 1984;1(6):612–619.
[2] Katsevich A. Image reconstruction for the circle and line trajectory. Phys. Med. Biol. 2004;49(22):5059–5072. [PubMed]
[3] Katsevich A. Analysis of an exact inversion algorithm for spiral conebeam CT. Phys. Med. Biol. 2002;47(15):2583–2597. [PubMed]
[4] Stierstorfer K, Rauscher A, Boese J, et al. Weighted fbp–A simple approximate 3D fbp algorithm for multislice spiral ct with good dose usage for arbitrary pitch. Phys. Med. Biol. 2004;49(11):2209–2218. [PubMed]
[5] Pan X. Optimal noise control in and fast reconstruction of fan-beam computed tomography image. Med. Phys. 1999;26(5):689–697. [PubMed]
[6] Zeng GL. Nonuniform noise propagation by using the ramp filter in fan-beam computed tomography. IEEE Trans. Med. Imag. 2004;23(6):690–695. [PubMed]
[7] Man BD, Basu S. Distance-driven projection and backprojection in three dimensions. Phys. Med. Biol. 2004;49(11):2463–2475. [PubMed]
[8] De Man B, Basu S. Distance-driven projection and backprojection. Proceedings of IEEE Nuclear Science Symposium Conference; Norfolk, Virginia. 2002. pp. 1477–1480.
[9] Wang J, Lu H, Li T, et al. An alternative solution to the nonuniform noise propagation problem in fan-beam fbp image reconstruction. Med. Phys. 2005;32(11):3389–3394. [PMC free article] [PubMed]
[10] Sunnegardh J, Danielsson PE. A new anti-aliased projection operator for iterative CT reconstruction. Proceedings of Fully 3D Meeting and HPIR Workshop; Lindau, Germany. 2007.
[11] Dennerlein F, Noo F, Lauritsch G, et al. Fanbeam filtered-backprojection reconstruction without backprojection weight. Phys. Med. Biol. 2007;52(11):3227–3240. [PubMed]
[12] Xia D, Cho S, Pan X. Backprojection filtration algorithm with improved noise property for a circular cone-beam CT. Proceedings of Fully 3D Meeting and HPIR Workshop; Beijing, China. 2009. pp. 170–173.
[13] Narasimhadhan AV, Anoop KP, Rajgopal K. FDK algorithms with no backprojection weight. Proceedings of Fully 3D Meeting and HPIR Workshop; Beijing, China. 2009. pp. 158–161.
[14] Parker DL. Optimal short scan convolution reconstruction for fan-beam CT. Med. Phys. 1982;9(2):254–257. [PubMed]
[15] Kak AC, Slaney M. Principles of Computerized Tomographic Imaging. IEEE Press; New York: 1988.
[16] Shepp LA, Logan BF. The fourier reconstruction of a head section. IEEE Trans. Nucl. Sci. 1974;21:21–43.
[17] Noo F, Defrise M, Clackdoyle R, et al. Image reconstruction from fan-beam projections on less than a short scan. Phys. Med. Biol. 2002;47(14):2525–2546. [PubMed]
[18] Noo F, Pack JD, Heuscher D. Exact helical reconstruction using native cone-beam geometries. Phys. Med. Biol. 2003;48(23):3787–3818. [PubMed]