International Journal of Biomedical Imaging
Volume 2008 (2008), Article ID 242841, 11 pages
doi:10.1155/2008/242841
Research Article

An Efficient Estimation Method for Reducing the Axial Intensity Drop in Circular Cone-Beam CT

Lei Zhu,1 Jared Starman,1,2 and Rebecca Fahrig1

1Department of Radiology, Stanford University, Stanford, CA 94305, USA
2Department of Electrical Engineering, Stanford University, Stanford, CA 94305, USA

Received 9 April 2008; Accepted 8 August 2008

Academic Editor: Jiang Hsieh

Copyright © 2008 Lei Zhu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Reconstruction algorithms for circular cone-beam (CB) scans have been extensively studied in the literature. Since insufficient data are measured, an exact reconstruction is impossible for such a geometry. If the reconstruction algorithm assumes zeros for the missing data, such as the standard FDK algorithm, a major type of resulting CB artifacts is the intensity drop along the axial direction. Many algorithms have been proposed to improve image quality when faced with this problem of data missing; however, development of an effective and computationally efficient algorithm remains a major challenge. In this work, we propose a novel method for estimating the unmeasured data and reducing the intensity drop artifacts. Each CB projection is analyzed in the Radon space via Grangeat_s first derivative. Assuming the CB projection is taken from a parallel beam geometry, we extract those data that reside in the unmeasured region of the Radon space. These data are then used as in a parallel beam geometry to calculate a correction term, which is added together with Hu’s correction term to the FDK result to form a final reconstruction. More approximations are then made on the calculation of the additional term, and the final formula is implemented very efficiently. The algorithm performance is evaluated using computer simulations on analytical phantoms. The reconstruction comparison with results using other existing algorithms shows that the proposed algorithm achieves a superior performance on the reduction of axial intensity drop artifacts with a high computation efficiency.

1. Introduction

Circular cone-beam (CB) scans are commonly used in X-ray CT. However, given Tuy's data sufficiency condition [1], the Radon space data of the scanned object cannot be completely measured in such an imaging geometry, and an exact reconstruction is possible only in the plane of the source trajectory (midplane). Many approximate reconstruction algorithms have been proposed in the literature. The FDK algorithm, developed by Feldkamp et al. [2], is by far the most popular mainly due to its structure of one-dimensional (1D) shift-invariant filtering and backprojection. Although originally derived as a heuristic extension of the exact fan-beam reconstruction, the FDK algorithm has been shown to be equivalent to an exact 3D reconstruction if the unmeasured Radon space data are assumed zeros [3, 4], except that a small correction term is also needed [5, 6]. Therefore, the FDK algorithm is exact for an object with uniform distribution in the longitudinal direction [2], whose Radon space data (the first derivative of the Radon transform) are exactly zeros in the unmeasured region of a circular trajectory. However, in general, zero is not a good approximation of the missing Radon space data, for the case when the scanned object is nonuniform and compactly supported, and has nonnegative attenuation coefficients. Consequently, the reconstructed images have CB artifacts, such as the well-known intensity drop in the axial direction [2, 5, 7–10]. The reconstruction can be improved by using an auxiliary trajectory in addition to the circular trajectory to measure the missing data [11–13]. In this work, we focus on using circular trajectories only, and develop an estimation-based method to reduce the intensity drop artifacts.

Reduction of artifacts in circular cone-beam CT (CBCT) can be achieved by estimating the unmeasured data using interpolation or extrapolation. Using Grangeat's formula [3], CB projection data can be analyzed as the first derivative of the Radon transform of the scanned object, and data estimation using interpolation can be performed in the Radon space to suppress CB artifacts [8]. However, this method requires three-dimensional (3D) data gridding, and hence it is computationally intensive; and in addition the backprojection structure is lost. Estimation methods of the Radon space data are also proposed in the space of reconstructed images, using multiple scans with different source-to-axis distances [7, 14, 15]. Zeng et al. developed improved algorithms to reduce the intensity drop using iterations [16]. These algorithms require either multiple reconstructions using different imaging geometric parameters or an iterative reconstruction that involves several computationally intense forward and backprojection steps. Other researchers have developed many improved algorithms in the framework of filtered-backprojection (FBP), without using Grangeat's framework [9, 10, 17, 18]. Since FBP has a compact structure, it is not easy to implement data interpolation/extrapolation in the Radon domain to further reduce the intensity drop artifacts.

Hu discovered that in a circular CB trajectory, the original object can be written as the summation of three terms [5]: 𝑓 𝑓 = F D K + 𝑓 𝐻 + 𝑓 𝑁 , ( 1 ) where 𝑓 is the original object, 𝑓 F D K is the FDK reconstruction, 𝑓 𝐻 is Hu's correction term which represents the information contained in a circular CB scan but not utilized in the FDK reconstruction, and 𝑓 𝑁 represents the information that is missing in the circular trajectory and cannot be reconstructed exactly. Hu proposed an algorithm that includes the first two measured terms, which shows reduced intensity drop in the reconstruction as compared to the FDK reconstruction [5]. Based on Grangeat's formula, Yang et al. proposed an algorithm that estimates the missing term 𝑓 𝑁 and effectively reduces the intensity drop artifacts [19]. However, the formula of the estimated term takes the form of shift-variant filtering and backprojection, two steps that both require intense computation.

The work presented in this paper is also based on Grangeat's formula and Hu's theory. However, the derivation of estimation formula for 𝑓 𝑁 is different from that of Yang, and the resulting implementation of the final formula is very efficient. We first analyze the CB projection data in the Radon domain via Grangeat's formula. Then, the unmeasured Radon space data are estimated from the CB projection by assuming that the projection is acquired in a parallel beam geometry. This approximation is equivalent to a data interpolation in the Radon domain. The estimated data are reconstructed as in the parallel beam geometry. More approximations are made to avoid the expensive steps of shift-variant filtering and backprojection in the calculation. The result is then added to the standard FDK reconstruction together with Hu's correction term to form the final reconstruction. Although the derivation assumes projections on a full scan, it can be readily extended to short-scan reconstructions using heuristic weighting schemes, such as Parker's weighting [20]. The performance of the proposed algorithm is verified using computer simulations on the 3D Shepp-Logan phantom, the FORBILD head phantom and the Defrise disc phantom. To fully evaluate the algorithm, the reconstructions are also compared with those using other existing algorithms.

The rest of this paper is organized as follows. Section 2 reviews the Radon transform and Grangeat's formula. The main algorithm is then derived and a reconstruction scheme is also proposed. Section 3 presents the results of computer simulations. Finally, Section 4 summarizes the paper.

2. Method

2.1. The System Geometry

The system geometry is shown in Figure 1. In this paper, we use an equally spaced flat panel detector with a finite size. Algorithms when other types of detectors are used can be derived similarly. During data acquisition, the X-ray source 𝑆 rotates about the 𝑧 axis in the 𝑥 - 𝑦 plane, with a fixed distance 𝐷 to the center of rotation 𝑂 . Angle 𝛾 𝑚 is the full fan angle determined by the size of the detector and the focal spot to detector distance. We derive the algorithm assuming that the range of the view angle 𝛽 is 3 6 0 degrees in a full-scan mode. The detector is placed perpendicular to → 𝑆 𝑂 for each projection. The object to be reconstructed is described by a compactly supported nonnegative function 𝑓 ( → 𝑟 ) , where → 𝑟 = ( 𝑥 , 𝑦 , 𝑧 ) is the Cartesian coordinate. In the derivation, it is also assumed that there is no truncation of the projection data; this condition, however, will be relaxed based on the final formula.

242841.fig.001
Figure 1: CB data acquisition geometry and coordinate system.

Denoting the distance from 𝑂 to detector as 𝐷 𝑜 𝑑 , the relationship between 𝑝 ( 𝑢 , 𝑣 , 𝛽 ) , the real projection image, and 𝑝 𝑣 ( 𝑢 , 𝑣 , 𝛽 ) , the image on a virtual detector that is parallel to the real detector and passes through 𝑂 , is as follows: 𝑝 𝑣 𝑢 ( 𝑢 , 𝑣 , 𝛽 ) = 𝑝 𝐷 + 𝐷 𝑜 𝑑 𝐷 , 𝑣 𝐷 + 𝐷 𝑜 𝑑 𝐷 , 𝛽 . ( 2 ) For simplicity, we use 𝑝 𝑣 in the reconstruction hereafter, and the parameter 𝛽 is dropped in the algorithm derivation if it is not used.

Main variables used in this paper are listed in Table 1 for clarity.

tab1
Table 1: Variable glossary.
2.2. The Radon Inversion and Grangeat's Formula

Reconstruction of the original object from its projections can be solved using the 3D Radon inverse formula [3, 21]. However, it is not efficient to directly use the Radon inversion on the X-ray projection data. To reduce the computational complexity, Grangeat established a fundamental relationship between the X-ray projection image and the first derivative of the Radon transform of the scanned object.

Figure 2 shows the geometric parameters of Grangeat's formula. One line 𝐿 on the projection image 𝑝 𝑣 can be specified by two parameters, its distance to the origin ‖ → 𝑂 𝑀 ‖ = 𝜌 , and the vector → 𝑚 in the image plane and perpendicular to 𝐿 . The line 𝐿 and the focal spot 𝑆 determine a plane 𝑃 , which can also be specified by two parameters, its normal vector → 𝑛 and its displacement to the origin, ‖ → √ 𝑂 𝑁 ‖ = 𝐷 𝜌 / 𝐷 2 + 𝜌 2 . Define an intermediate function 𝑆 𝑝 as a weighted sinogram of the 2D projection image on the virtual detector 𝑝 𝑣 : 𝑆 𝑝 𝜌 , → 𝑚 = → 𝑡 ∈ 𝐿 ( 𝜌 , → 𝑚 ) 𝐷 √ 𝐷 2 + 𝑢 2 + 𝑣 2 𝑝 𝑣 → 𝑡 𝑑 → 𝑡 , ( 3 ) where → 𝑡 = ( 𝑢 , 𝑣 ) is a vector in the plane of the projection image, and √ 𝐷 / 𝐷 2 + 𝑢 2 + 𝑣 2 is the cosine weighting for an oblique incident angle.

242841.fig.002
Figure 2: Illustration of the geometry of Grangeat's formula.

Denote 𝑅 𝑓 as the Radon transform of the scanned object 𝑓 , which is specified by a unit vector → 𝑛 0 and a scalar 𝜌 0 : 𝜌 𝑅 𝑓 0 , → 𝑛 0 = → 𝑟 ∈ 𝑃 0 ( 𝜌 0 , → 𝑛 0 ) 𝑓 → 𝑟 𝑑 → 𝑟 , ( 4 ) where 𝑃 0 is a 2D plane, with a normal vector → 𝑛 0 and a displacement 𝜌 0 from the origin.

The relationship between the first derivative of 𝑆 𝑝 and the first derivative of 𝑅 𝑓 can be found: ‖ ‖ → ‖ ‖ 𝑂 𝑆 2 ‖ ‖ ⇒ ‖ ‖ 𝑁 𝑆 2 𝑆 𝑝 ‖ ‖ → ‖ ‖ , 𝑂 𝑀 → 𝑚 = 𝑅 𝑓 ‖ ‖ → ‖ ‖ , 𝑂 𝑁 → 𝑛 ⇒ 𝜌 2 + 𝐷 2 𝐷 2 𝑆 𝑝 𝜌 , → 𝑚 = 𝑅 𝑓 ( 𝐷 𝜌 √ 𝐷 2 + 𝜌 2 , → 𝑛 ) , ( 5 ) where both of the first derivatives are with regard to the first parameter.

Based on Grangeat's formula (5), Figure 3 shows the data supports in the domain of 𝑅 𝑓 for a circular CB trajectory. For simplicity, hereafter we refer to the domain of the first derivative of the Radon transform as the Radon domain, and the data in this domain as the Radon space data. For one CB projection, the surface of a sphere is measured in the Radon domain; as the projection rotates, the sphere rotates as well, and after a full scan, a torus of data are measured. It can be seen that not all of the Radon space is measured, and the missing data problem is clear. The diameter of the sphere of one projection is 𝐷 , the distance from the focal spot to the center of rotation. As this distance increases, the region of missing data gets smaller. In a parallel beam geometry (with an infinite 𝐷 ), this sphere surface becomes a plane, and the missing data problem disappears.

242841.fig.003
Figure 3: Data supports in the Radon domain.

This formalism provides a clearer understanding of the three terms in (1), which is consistent with Hu's original arguments [5]. The FDK reconstruction 𝑓 F D K represents all the information inside the torus of measured data and partial information on the surface of the torus; Hu's correction term 𝑓 𝐻 compensates for the unused data on the surface of the torus; the term 𝑓 𝑁 represents the missing information outside the torus, and it is our goal to estimate these data.

2.3. The Algorithm Derivation

The missing data represented by the term 𝑓 𝑁 result in CB artifacts in the reconstructed images. One solution to this problem is estimating this term using interpolation in the Radon domain. Inspired by the fact that a parallel beam geometry is free of the missing data problem, the missing data are estimated from measured CB projections, assuming that they are acquired in a parallel geometry. Then reconstruction is carried out using only these estimated data as though acquired in a parallel beam geometry, and the result is added to the standard FDK term and Hu's term as a correction to form a final reconstruction. Using this method, partial data on the spherical surface of one CB projection as shown in Figure 3 are filled in the region of missing data, and a data interpolation/extrapolation is carried out implicitly.

The concept of implicit data interpolation/extrapolation is illustrated in the Radon space in Figure 4. As similar to Figure 3, the sphere surfaces in the figure are measured using CB projections. If we assume the CB projections are acquired in a parallel beam geometry and reconstruct the object using a parallel beam geometry as well, the measured sphere surfaces become planes, that is, data points on the measured sphere surfaces are relocated onto the planes. Note that measured sphere surfaces associated with CB projections from the opposite projection angles are approximated as the same plane in the Radon space. As shown in Figure 4, data point 𝐶 on the plane maps to two points 𝐴 and 𝐵 associated with two different CB projections. After the CB to parallel approximation, both data points 𝐴 and 𝐵 are “moved” to data point 𝐶 ; equivalently, the approximated data point 𝐶 is a weighted summation of data points 𝐴 and 𝐵 . Therefore, this process can be considered as a linear interpolation of data point 𝐶 from measured data points 𝐴 and 𝐵 .

242841.fig.004
Figure 4: Illustration of the concept of implicit data interpolation/extrapolation. Data point 𝐶 on the approximated plane is a weighted summation of measured data points 𝐴 and 𝐵 .

Using the idea of implicit interpolation/extrapolation, we can estimate Radon space data on a plane from two CB projections. As the projection angle changes, all the Radon space data can be estimated. However, we only need to estimate the data in the missing region of the Radon space in a CB geometry, and use them in the correction of the artifacts. Figure 5 is a 2D plot of the vertical plane in the Radon domain in Figure 3. If a parallel beam geometry is used in the data acquisition, the data of this plane can be measured completely using one projection from the direction of the normal of the plane. In a CB geometry, however, the measured data in the plane consists of two discs (shown as shaded regions), and the data outside are missing. If the CB projection from the direction of the normal of the plane is assumed to be a parallel projection, the measured data can be calculated in the whole plane; then only the data outside the shaded region are used in the reconstruction of the correction term 𝑓 𝑁 . The data of 𝑅 𝑓 in the missing region can be separated from the measured data by multiplying by 𝑤 ( 𝑠 , 𝑡 ) , a map function that is zero in the shaded region of Figure 5, and one outside. Denote 𝑝 𝑝 as the projection image in a parallel beam geometry from the direction of the normal of the plane, and 𝑆 𝑝 𝑝 as its sinogram. The relationship between 𝑅 𝑓 and 𝑆 𝑝 𝑝 can be found using Grangeat's formula, by letting 𝐷 go to infinity in (5): 𝑆 𝑝 𝑝 𝜌 , → 𝑚 = 𝑅 𝑓 𝜌 , → 𝑛 . ( 6 ) From Figure 2, it is also seen that the unit vector → 𝑚 becomes identical to → 𝑛 in a parallel geometry. Let → 𝑚 or → 𝑛 = ( c o s 𝛼 , s i n 𝛼 ) . For simplicity, the functions 𝑆 𝑝 𝑝 ( 𝜌 , → 𝑚 ) and 𝑅 𝑓 ( 𝜌 , → 𝑛 ) are rewritten as 𝑆 𝑝 𝑝 ( 𝜌 , 𝛼 ) and 𝑅 𝑓 ( 𝜌 , 𝛼 ) in the polar coordinate system. Based on Figure 5, the map function 𝑤 in the ( 𝜌 , 𝛼 ) coordinate system can be found as | | 𝜌 | | > | | | | , 𝑤 ( 𝜌 , 𝛼 ) = 1 , 𝐷 c o s 𝛼 0 , o t h e r w i s e . ( 7 )

242841.fig.005
Figure 5: A vertical plane in the Radon domain.

Now the multiplication of 𝑤 ( 𝜌 , 𝛼 ) on 𝑅 𝑓 ( 𝜌 , 𝛼 ) can be done directly on 𝑆 𝑝 𝑝 , according to (6): 𝑆 𝑝 𝑝 ( 𝜌 , 𝛼 ) 𝑤 ( 𝜌 , 𝛼 ) = 𝑅 𝑓 ( 𝜌 , 𝛼 ) 𝑤 ( 𝜌 , 𝛼 ) . ( 8 )

We then approximate the parallel projection data 𝑆 𝑝 𝑝 using the CB projection data 𝑆 𝑝 . Mathematically, the following approximation is made: 𝜌 2 + 𝐷 2 𝐷 2 𝑆 𝑝 ( 𝜌 , 𝛼 ) 𝑤 ( 𝜌 , 𝛼 ) ≈ 𝑆 𝑝 𝑝 ( 𝜌 , 𝛼 ) 𝑤 ( 𝜌 , 𝛼 ) . ( 9 ) Note that 𝐷 goes to infinity in a parallel beam geometry, and the weight ( 𝜌 2 + 𝐷 2 ) / 𝐷 2 equals one on the right-hand side. Equations (8) and (9) show that the missing Radon space data are approximated using the measured CB projection data.

We first derive the correction term that compensates for the missing data in terms of 𝑆 𝑝 𝑝 , and then apply the approximation shown in (9). To calculate the correction term, reconstruction must be carried out using only partial projection data in the domain of 𝑆 𝑝 𝑝 . This can be done by using a 2D reconstruction of the projection data, followed by a 3D parallel reconstruction. Note that these intermediate steps are used only in the derivation, and the derived final formula is implemented efficiently without these steps.

Since 𝑆 𝑝 𝑝 is the sinogram of the projection image, we can compute the parallel projection image 𝑝 𝑝 using FBP in a 2D parallel geometry: 𝑝 𝑝 1 ( 𝑢 , 𝑣 ) = 2 0 2 𝜋 ∞ − ∞ 𝑆 𝑝 𝑝 ( 𝜌 , 𝛼 ) 𝑔 0 = 1 ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 − 𝜌 ) 𝑑 𝜌 𝑑 𝛼 ( 1 0 ) 4 𝜋 0 2 𝜋 ∞ − ∞ 𝑆 𝑝 𝑝 ( 𝜌 , 𝛼 ) 𝑔 ℎ ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 − 𝜌 ) 𝑑 𝜌 𝑑 𝛼 , ( 1 1 ) where 𝑔 0 ( 𝑢 ) is the ramp-filter kernel and 𝑔 ℎ ( 𝑢 ) = 1 / 𝜋 𝑢 is the Hilbert kernel. In (11), we use a Hilbert transform after a first derivative operation to substitute for the ramp-filtering, so that the calculation can be directly applied on 𝑆 𝑝 𝑝 .

Denote 𝑝 𝑚 as the projection image reconstructed from partial data of 𝑆 𝑝 𝑝 . Using (11), we have 𝑝 𝑚 = 1 ( 𝑢 , 𝑣 ) 4 𝜋 0 2 𝜋 ∞ − ∞ 𝑆 𝑝 𝑝 ( 𝜌 , 𝛼 ) 𝑤 ( 𝜌 , 𝛼 ) × 𝑔 ℎ ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 − 𝜌 ) 𝑑 𝜌 𝑑 𝛼 ( 1 2 )

Denote 𝑓 𝑐 as an estimate of the missing data 𝑓 𝑁 in (1). We can compute the correction term 𝑓 𝑐 using the 3D parallel reconstruction of 𝑝 𝑚 . This computation can be done as slices of 2D parallel reconstructions. The FBP-based reconstruction is 𝑓 𝑐 1 ( 𝑥 , 𝑦 , 𝑧 ) = 2 0 2 𝜋 𝑝 𝐹 ( 𝑥 c o s 𝛽 + 𝑦 s i n 𝛽 , 𝑧 , 𝛽 ) 𝑑 𝛽 , ( 1 3 ) where 𝑝 𝐹 is the ramp-filtered parallel projection at view angle 𝛽 : 𝑝 𝐹 ( 𝑢 , 𝑣 ) = ∞ − ∞ 𝑝 𝑚 𝑔 𝑢 , 𝑣 0 𝑢 − 𝑑 𝑢 , 𝑣 𝑢 = 1 ( 1 4 ) 2 𝜋 ∞ − ∞ 𝑝 𝑚 𝑔 𝑢 , 𝑣 ℎ 𝑢 − 𝑑 𝑢 , 𝑣 𝑢 . ( 1 5 )

Calculation of 𝑓 𝑐 using (12), (13), and (15) can be simplified. Based on the two-step method developed by Noo et al. [22], 2D parallel reconstruction can be carried out using a derivative backprojection followed by a 1D Hilbert transform on the reconstructed image. Therefore, the Hilbert transform is taken out of the integral: 𝑝 𝑚 = 1 ( 𝑢 , 𝑣 ) 4 𝜋 ∞ − ∞ 𝑔 ℎ ( 𝑢 ) 0 2 𝜋 s g n ( c o s 𝛼 ) × 𝑆 𝑝 𝑝 𝑢 − 𝑢 ) c o s 𝛼 + 𝑣 s i n 𝛼 , 𝛼 × 𝑤 𝑢 − 𝑢 ) c o s 𝛼 + 𝑣 s i n 𝛼 , 𝛼 𝑑 𝛼 𝑑 𝑢 , ( 1 6 ) where sgn is the signum function.

Define the function 𝑝 ℎ as 𝑝 ℎ ( 𝑢 , 𝑣 ) = 0 2 𝜋 s g n ( c o s 𝛼 ) 𝑆 𝑝 𝑝 ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 , 𝛼 ) × 𝑤 ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 , 𝛼 ) 𝑑 𝛼 . ( 1 7 )

Then, we have 𝑝 𝑚 1 ( 𝑢 , 𝑣 ) = 4 𝜋 ∞ − ∞ 𝑔 ℎ ( 𝑢 ) 𝑝 ℎ ( 𝑢 − 𝑢 , 𝑣 ) 𝑑 𝑢 . ( 1 8 ) Both first derivatives are with regard to the first parameter.

Insert (18) into (15), and the two Hilbert transforms become − 1 since the Hilbert transform of a Hilbert transform of a function equals the negative of the function: 𝑝 𝐹 1 ( 𝑢 , 𝑣 ) = − 8 𝜋 2 𝑝 ℎ ( 𝑢 , 𝑣 ) . ( 1 9 )

The problem is simplified to the calculation of 𝑝 ℎ . Applying the approximation in (9), we have 𝑝 ℎ ≈ 𝜕 ( 𝑢 , 𝑣 ) ( 𝜕 𝑢 0 2 𝜋 s g n ( c o s 𝛼 ) ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 ) 2 + 𝐷 2 𝐷 2 × 𝑆 = 𝑝 ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 , 𝛼 ) 𝑤 ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 , 𝛼 ) 𝑑 𝛼 ) 0 2 𝜋 s g n ( c o s 𝛼 ) ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 ) 2 + 𝐷 2 𝐷 2 × 𝜕 𝑆 𝜕 𝑢 + 𝑝 ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 , 𝛼 ) 𝑤 ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 , 𝛼 ) 𝑑 𝛼 0 2 𝜋 𝜕 s g n ( c o s 𝛼 ) ( 𝜕 𝑢 ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 ) 2 + 𝐷 2 𝐷 2 ) × 𝑆 = 𝑝 ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 , 𝛼 ) 𝑤 ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 , 𝛼 ) 𝑑 𝛼 ( 2 0 ) 0 2 𝜋 s g n ( c o s 𝛼 ) ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 ) 2 + 𝐷 2 𝐷 2 c o s 𝛼 × 𝑆 + 𝑝 ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 , 𝛼 ) 𝑤 ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 , 𝛼 ) 𝑑 𝛼 0 2 𝜋 s g n ( c o s 𝛼 ) 2 ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 ) 𝐷 2 c o s 𝛼 × 𝑆 = 𝑝 ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 , 𝛼 ) 𝑤 ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 , 𝛼 ) 𝑑 𝛼 0 2 𝜋 | | | | c o s 𝛼 𝑤 ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 , 𝛼 ) × ( ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 ) 2 + 𝐷 2 𝐷 2 𝑆 + 𝑝 ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 , 𝛼 ) 2 ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 ) 𝐷 2 𝑆 𝑝 ( 𝑢 c o s 𝛼 + 𝑣 s i n 𝛼 , 𝛼 ) ) 𝑑 𝛼 . ( 2 1 ) The discontinuous positions of the function 𝑤 correspond to the points on the surface of the torus of measured data in Figure 3, which are compensated for by Hu's correction term. Therefore, in (20), the derivatives on the discontinuities of 𝑤 are not included.

Equation (21) has a structure of shift-variant filtering (due to the shift-variance of the multplication by the weighting function) and backprojection, therefore the implementation is not very efficient. Furthermore, since a weighted sinogram of each 2D projection image is needed in the calculation, it is required that no truncation is present in the projection. In particular, the method suffers from projection truncation in the longitudinal direction, so-called long object problem.

We simplify (21) using further approximations. Since the weighting of 𝑆 𝑝 is usually much larger than that of 𝑆 𝑝 , that is, ( 𝜌 2 + 𝐷 2 ) / 𝐷 2 versus 2 𝜌 / 𝐷 2 , the second term associated with 𝑆 𝑝 is ignored. For a circular trajectory with a not very large cone angle, we have 𝐷 ≫ 𝜌 . The weighting function 𝑤 is nonzero only when | c o s 𝛼 | < | 𝜌 / 𝐷 | , that is, in small neighborhoods of 𝛼 at 𝜋 / 2 and 3 𝜋 / 2 . In these neighborhoods, the following approximations can be made: 𝑆 𝑣 s i n 𝛼 ≈ s g n ( s i n 𝛼 ) 𝑣 , 𝑢 c o s 𝛼 ≈ 0 , 𝑝 s g n ( s i n 𝛼 ) 𝑣 , 𝛼 ≈ 𝑆 𝑝 𝜋 𝑣 , 2 . ( 2 2 )

Now the calculation of 𝑝 ℎ can be simplified as 𝑝 ℎ ( 𝑢 , 𝑣 ) ≈ 2 𝜋 0 | | | | 𝑣 c o s 𝛼 𝑤 ( 𝑣 , 𝛼 ) 2 + 𝐷 2 𝐷 2 𝑆 𝑣 𝑝 ( 𝑣 , 𝛼 ) 𝑑 𝛼 ≈ 2 2 + 𝐷 2 𝐷 2 𝑆 𝑝 𝜋 𝑣 , 2 𝜋 / 2 + a r c c o s ( | | 𝑣 | | / 𝐷 ) 𝜋 / 2 − a r c c o s ( | | 𝑣 | | / 𝐷 ) | | | | 𝑣 c o s 𝛼 𝑑 𝛼 = 4 2 + 𝐷 2 𝐷 2 √ ( 1 − 𝐷 2 − 𝑣 2 𝐷 ) 𝑆 𝑝 𝜋 𝑣 , 2 𝑣 = 4 2 + 𝐷 2 𝐷 2 √ ( 1 − 𝐷 2 − 𝑣 2 𝐷 ) × 𝜕 2 𝜕 𝑣 2 ( ∞ − ∞ 𝐷 √ 𝑢 2 + 𝑣 2 + 𝐷 2 𝑝 𝑣 ( 𝑢 , 𝑣 , 𝛽 ) 𝑑 𝑢 ) . ( 2 3 ) The correction term 𝑓 𝑐 can be calculated by combining (13), (19), and (23).

Note that the approximate 𝑝 ℎ (or 𝑝 𝐹 ) is a function independent of parameter 𝑢 . Backprojection of 𝑝 𝐹 shown in (13) can be implemented efficiently by simply changing variable 𝑣 in the projection space to 𝑧 in the reconstruction space.

2.4. Practical Reconstruction Scheme

The final reconstruction is the summation of the FDK reconstruction 𝑓 F D K , Hu's term 𝑓 𝐻 , and the correction term 𝑓 𝑐 : 𝑓 𝑓 = F D K + 𝑓 𝐻 + 𝑓 𝑐 . ( 2 4 ) The practical implementation of this formula is summarized below. The derivations of the FDK reconstruction and Hu's term can be found in [2, 5], respectively, and the formulae are presented here as a reference: 𝑓 F D K 1 ( 𝑥 , 𝑦 , 𝑧 ) = 2 0 2 𝜋 𝐷 𝑥 + 𝐷 2 𝛽 𝑞 𝐹 ( 𝐷 𝑦 𝑥 + 𝐷 𝛽 , 𝐷 𝑧 𝑥 + 𝐷 𝛽 𝑞 , 𝛽 ) 𝑑 𝛽 , 𝐹 ( 𝑢 , 𝑣 , 𝛽 ) = ∞ − ∞ 𝐷 √ 𝑢 2 + 𝑣 2 + 𝐷 2 𝑝 𝑣 𝑔 𝑢 , 𝑣 , 𝛽 0 𝑢 − 𝑢 𝑑 𝑢 , ( 2 5 ) where the subscript 𝛽 stands for the coordinate transformation of rotation about the 𝑧 axis by 𝛽 ; and 𝑓 𝐻 1 ( 𝑥 , 𝑦 , 𝑧 ) = − 2 𝜋 0 2 𝜋 ( 𝑧 ( 𝑥 + 𝐷 ) 2 ) 𝛽 𝑞 𝐷 ( 𝐷 𝑧 𝑥 + 𝐷 𝛽 𝑞 , 𝛽 ) 𝑑 𝛽 , 𝐷 1 ( 𝑣 , 𝛽 ) = 𝜕 2 𝜋 𝜕 𝑣 ∞ − ∞ 𝐷 √ 𝑢 2 + 𝑣 2 + 𝐷 2 𝑝 𝑣 𝑑 𝑢 , 𝑣 , 𝛽 𝑢 . ( 2 6 )

𝑓 𝑐 is calculated by using (13), (19), and (23): 𝑓 𝑐 1 ( 𝑥 , 𝑦 , 𝑧 ) = − 4 𝜋 2 𝑧 2 + 𝐷 2 𝐷 2 √ ( 1 − 𝐷 2 − 𝑧 2 𝐷 ) × 0 2 𝜋 𝜕 2 𝜕 𝑧 2 ( ∞ − ∞ 𝐷 √ 𝑢 2 + 𝑧 2 + 𝐷 2 𝑝 𝑣 ( 𝑢 , 𝑧 , 𝛽 ) 𝑑 𝑢 ) 𝑑 𝛽 . ( 2 7 ) This estimation formula of the missing term 𝑓 𝑁 in (1) is the main result of the paper. The equation shows a simple structure of calculation. Note that, since the second derivative operation is very sensitive to high-frequency errors and the intensity drop artifacts are mostly low-frequency signals in the longitudinal direction, filtering techniques are used to suppress the errors in the calculation. A practical implementation can be divided into the following steps.

Step 1. Take each projection 𝑝 𝑣 ( 𝑢 , 𝑣 , 𝛽 ) , integrate along the direction of 𝑢 .Step 2. Take the second derivative with respect to variable 𝑣 .Step 3. Filter the 1D profile obtained from Step 2 using a median filter and a window filter.Step 4. Integrate the processed 1D profile along the projection angle 𝛽 .Step 5. Change coordinate variable from the projection space to the reconstruction space (from 𝑣 to 𝑧 ).Step 6. Weight by ( − 1 / 4 𝜋 2 ) ( ( 𝑧 2 + 𝐷 2 ) / 𝐷 2 √ ) ( 1 − 𝐷 2 − 𝑧 2 / 𝐷 ) .Step 7. Replicate the 1D profile of 𝑧 in the directions of 𝑥 and 𝑦 to generate a 3D volume. As discussed earlier, Step 3 is important to suppress high-frequency misestimation and remove the streak artifacts that are otherwise present in 𝑓 𝑐 . The median filter is able to remove high spikes caused by object boundaries, and the window filter is able to smooth out small fluctuations. In all the implementations presented in this paper, we used a median filter with a width of 10 pixels and a Hamming window filter. It is worth mentioning that since the filtering is applied only on 𝑓 𝑐 to enforce low-frequency estimation, it will not affect the resolution of the reconstructed images obtained by the first two terms in (24).

In (27), we take the second derivative of the projection images along the vertical direction, and therefore the proposed algorithm can survive the long object problem. The calculation of (27) is also very efficient, since neither a shift-variant filtering step nor a backprojection step is used. This feature makes the proposed algorithm distinct from other existing algorithms, such as Yang's method [19]. As will be shown in the section of numerical results, in our implementations, the proposed method is typically 7-8 times more efficient than Yang's method. Note that the calculation of Hu's term (26) has the same FBP structure as the FDK reconstruction, and the cone-beam backprojection steps of these two calculations can be combined to reduce the computation cost. Since in FBP reconstructions, the backprojection step takes the majority of the computation time, the computation complexity of the proposed reconstruction (24) is close to that of the FDK reconstruction only.

3. Numerical Results

3.1. Simulation Details

The algorithm performance was evaluated using computer simulations. Table 2 summarizes the system parameters used in the simulations. Three computer phantoms were used in this study. The first was the 3D Shepp-Logan phantom as defined in [23], which contains low-contrast objects. The second was the FORBILD head phantom (http://www.imp.uni-erlangen.de/forbild/). This phantom contains high-contrast objects, and therefore it results in more missing data in circular CBCT geometry. To further verify the algorithm, the Defrise disc phantom was also used. The Defrise phantom consists of seven ellipsoidal discs stacked in the 𝑧 direction. Each disc has a uniform attenuation coefficient of 1 m m − 1 , and the ellipsoid has a diameter of 1 4 0  mm and a thickness of 1 4 mm, with a distance of 2 5  mm between discs. This phantom has strong high-frequency components in the 𝑧 direction, and therefore has high values in its first derivative of the Radon transform in the region where the data are unmeasured in a circular CB scan. It represents the most challenging case of reconstruction using the circular CB data.

tab2
Table 2: Simulation parameters.

Simulations were carried out on full-scan data. To test the stability of the algorithm, reconstructions on noisy data of the Shepp-Logan phantom were also investigated. In the simulation, we used 3 0 0 0 0 0 photons per ray, and the base material of the Shepp-Logan phantom is modeled as water at 80 keV, with a linear attenuation coefficient of 0 . 0 1 8 3 7 m m − 1 .

To demonstrate the merit of the proposed algorithm, we compared the reconstructions using five different algorithms: the FDK algorithm [2], Hu's algorithm (only the first two terms of (24) are included) [5], the T-FDK algorithm [10], Yang's algorithm [19], and the proposed algorithm (24). All the five algorithms are in the category of analytical reconstruction algorithms for circular CBCT. As discussed earlier in the introduction section, the T-FDK algorithm was developed heuristically with a structure of shift-invariant FBP, and Yang's algorithm was based on interpolation in the Radon space, with a structure of shift-variant FBP.

3.2. Reconstruction Results

Figures 6, 9, and 10 show the reconstructed images on a full scan of the Shepp-Logan phantom, the FORBILD head phantom and the Defrise disc phantom, respectively. Comparisons of 1D vertical profiles of these images are also shown in Figures 8 and 11. The algorithm performances on the FORBILD head phantom are close to those on the Shepp-Logan phantom, and the 1D profile comparison of the reconstructions of the FORBILD head phantom is omitted here.

242841.fig.006
Figure 6: Reconstructions of the Shepp-Logan phantom. Top row: 𝑥 - 𝑧 views; bottom row: 𝑦 - 𝑧 views. Display window: [ 0 . 9 8 1 . 0 5 ] .
242841.fig.007
Figure 7: Images of the second term (Hu's term) and the third term (the proposed correction term) of (24), using the projection data on the Shepp-Logan phantom. Top row: 𝑥 - 𝑧 views; bottom row: 𝑦 - 𝑧 views. Display window: [ 0 0 . 0 8 ] .
fig8
Figure 8: 1D central vertical profile comparison of Figure 6. (a) The reduction of the axial intensity drop. (b) The reconstructed object edge using different algorithms.
242841.fig.009
Figure 9: Reconstructions of the FORBILD head phantom. Top row: 𝑥 - 𝑧 views; bottom row: 𝑦 - 𝑧 views. Display window: [ 1 . 0 2 1 . 0 8 ] .
242841.fig.010
Figure 10: Reconstructions of the Defrise phantom, 𝑥 - 𝑧 views. Display window: [ 0 . 7 1 . 1 ] .
242841.fig.011
Figure 11: 1D central vertical profile comparison of Figure 10.

From the comparison, it is seen that the CB artifacts of axial intensity drop are apparent in the reconstructions using the FDK algorithm or Hu's algorithm, and these artifacts are effectively suppressed using the proposed method.

The performance of the T-FDK algorithm on the reduction of intensity drop artifacts is also inferior to that of the proposed algorithm, as shown in the comparison of reconstructions. It is worth mentioning that the T-FDK algorithm is slightly more efficient than the FDK algorithm [10], although it causes resolution loss [18]. Recall that the proposed algorithm requires computation close to that of the FDK algorithm. Therefore, the T-FDK algorithm is more efficient than the proposed algorithm. Note that the reconstruction results using Hu's algorithm and the T-FDK algorithm are very similar. This similarity holds under certain conditions, and readers can refer to [24] for detailed discussion.

Compared to Yang's algorithm, our proposed algorithm has an advantage of high computation efficiency. As discussed in Section 2.4, the correction term in our algorithm only involves integration and derivative operations, and it is computed very efficiently. Yang's correction term has a structure of shift-variant filtering and backprojection; both steps require intense computation. In our implementations, Yang's algorithm is typically 7-8 times slower than our proposed algorithm.

Both the proposed algorithm and Yang's algorithm achieve similar reduction of the intensity drop away from the object edge. Nonetheless, estimation of high-frequency data in Yang's algorithm causes relatively large errors. The resulting artifacts are around the object edges in the axial direction, typically streaks with negative magnitudes. This problem can be seen in Figure 8(b), and is more obvious in Figure 11, where data estimation is more challenging due to the high-frequency missing data along the axial direction. Our algorithm does not estimate the high-frequency missing data, and the negative streaks do not appear around the object edges.

Reconstructions on noisy projection data are shown in Figure 12 to demonstrate the stability of the proposed algorithm in the presence of noise. The algorithm performance is similar to that in previous comparisons. Based on the noise-free reconstructions as shown in Figure 6, the noise variances in the images are also measured. The noise in the image using the proposed algorithm remains in the same level as that in the image using the FDK algorithm or Hu's algorithm.

242841.fig.012
Figure 12: Reconstructions of the Shepp-Logan phantom, using noisy projection data, 𝑥 - 𝑧 views. Display window: [ 0 . 9 8 1 . 0 5 ] . Based on the noise-free reconstructions as shown in Figure 6, from left to right, the noise variances in the images are measured as 5 . 8 7 2 × 1 0 − 5 , 5 . 8 7 2 × 1 0 − 5 , and 5 . 8 7 4 × 1 0 − 5 .

4. Conclusions and Discussion

In this work, we propose an efficient estimation method to reduce the intensity drop in the CB reconstruction on circular scans. The algorithms are derived using data analysis in the Radon domain via Grangeat's formula. Assuming the CB projections are measured in a parallel beam geometry, we estimate the unmeasured data from the CB projections. These data are then reconstructed to form a correction term to improve the FDK reconstruction with Hu's term included. Equivalently, an implicit data interpolation/extrapolation is carried out in the Radon domain. It is interesting to note that Hu's term takes the first derivative of the projection data along the axial direction, while our correction term takes the second derivative. Although our algorithm is derived for circular CBCT on a full scan, it can be easily extended to short-scan reconstructions using weighting schemes, such as Parker's weighting.

The algorithm performances are evaluated using computer simulations on the 3D Shepp-Logan phantom, the FORBILD head phantom, and the Defrise disc phantom. The results show that the proposed method greatly suppresses the axial intensity drop in the FDK reconstructions and its performance improves on Hu's algorithm. Residual artifacts are mainly due to the high-frequency Radon space data missing in a circular CB geometry, which cannot be estimated accurately using interpolation or extrapolation in general. As demonstrated in the results of the Defrise disc phantom, relatively large reconstruction errors are expected around high intensity objects, such as bones in a clinical case, in the longitudinal direction.

Our algorithm also outperforms the T-FDK algorithm on the reduction of intensity drop artifacts. As compared to other existing algorithms, such as Yang's algorithm, that are based on interpolation in the Radon space, our algorithm has an advantage of high efficiency. The calculation of the correction term requires only simple integration, 1D derivative and multiplication operations, and the total computation of the proposed algorithm is close to that of the FDK algorithm. In our implementations, the proposed algorithm is 7-8 times faster than Yang's algorithm.

Acknowledgments

This project is supported in part by NIH R01 EB003524 and the Lucas Foundation.

References

  1. H. K. Tuy, “An inversion formula for cone-beam reconstruction,” SIAM Journal on Applied Mathematics, vol. 43, no. 3, pp. 546–552, 1983.
  2. I. A. Feldkamp, L. C. Davis, and J. W. Kress, “Practical cone-beam algorithm,” Journal of the Optical Society of America A, vol. 1, no. 6, pp. 612–619, 1984.
  3. P. Grangeat, “Mathematical framework of cone-beam 3D reconstruction via the first derivativeof the Radon transform,” in Mathematical Methods in Tomography, vol. 1497 of Lecture Notes in Mathematics, pp. 66–97, Springer, Berlin, Germany, 1991.
  4. M. Defrise and R. Clack, “Cone-beam reconstruction algorithm using shift-variant filtering and cone-beam backprojection,” IEEE Transactions on Medical Imaging, vol. 13, no. 1, pp. 186–195, 1994.
  5. H. Hu, “An improved cone-beam reconstruction algorithm for the circular orbit,” Scanning, vol. 18, pp. 572–581, 1996.
  6. T.-L. Zhuang, B. E. Nett, X. Tang, and G.-H. Chen, “A cone-beam FBP reconstruction algorithm for short-scan and super-short-scan source trajectories,” in Developments in X-Ray Tomography IV, vol. 5535 of Proceedings of SPIE, pp. 566–576, Denver, Colo, USA, August 2004.
  7. F. Noo and D. J. Heuscher, “Image reconstruction from cone-beam data on a circular short-scan,” in Medical Imaging 2002: Image Processing, vol. 4684 of Proceedings of SPIE, pp. 50–59, San Diego, Calif, USA, February 2002.
  8. S. W. Lee and G. Wang, “A Grangeat-type half-scan algorithm for cone-beam CT,” Medical Physics, vol. 30, no. 4, pp. 689–700, 2003.
  9. H. Turbell, Cone-beam reconstruction using filtered backprojection, Ph.D. thesis, Department of Electrical Engineering, Linkopings Universitet, Linkoping, Sweden, February 2001.
  10. M. Grass, Th. Köhler, and R. Proksa, “3D cone-beam CT reconstruction for circular trajectories,” Physics in Medicine and Biology, vol. 45, no. 2, pp. 329–347, 2000.
  11. X. Wang and R. Ning, “A cone-beam reconstruction algorithm for circle-plus-arcdata-acquisition geometry,” IEEE Transactions on Medical Imaging, vol. 18, no. 9, pp. 815–824, 1999.
  12. G. L. Zeng and G. T. Gullberg, “A cone-beam tomography algorithm for orthogonal circle-and-line orbit,” Physics in Medicine and Biology, vol. 37, no. 3, pp. 563–577, 1992.
  13. X. Tang and R. Ning, “A cone beam filtered backprojection (CB-FBP) reconstruction algorithm for a circle-plus-two-arc orbit,” Medical Physics, vol. 28, no. 6, pp. 1042–1055, 2001.
  14. H. Kudo and T. Saito, “Derivation and implementation of a cone-beam reconstruction algorithm for nonplanar orbits,” IEEE Transactions on Medical Imaging, vol. 13, no. 1, pp. 196–211, 1994.
  15. K. Zeng, Z. Chen, L. Zhang, and G. Wang, “A half-scan error reduction based algorithm for cone-beam CT,” Journal of X-Ray Science and Technology, vol. 12, no. 2, pp. 73–82, 2004.
  16. K. Zeng, Z. Chen, L. Zhang, and G. Wang, “An error-reduction-based algorithm for cone-beam computed tomography,” Medical Physics, vol. 31, no. 12, pp. 3206–3212, 2004.
  17. X. Tang, J. Hsieh, A. Hagiwara, R. Nilsen, C. Shaughnessy, and E. Drapkin, “Axial reconstruction algorithms for cone-beam volumetric CT,” in Proceedings of the Radiological Society of North America (RSNA '05), Chicago, Ill, USA, November-December 2005.
  18. L. Yu, X. Pan, and C. A. Pelizzari, “Image reconstruction with a shift-variant filtration in circular cone-beam CT,” International Journal of Imaging Systems and Technology, vol. 14, no. 5, pp. 213–221, 2005.
  19. H. Yang, M. Li, K. Koizumi, and H. Kudo, “FBP-type cone-beam reconstruction algorithm with Radon space interpolationcapabilities for axially truncated data from a circular orbit,” Medical Imaging Techonology, vol. 24, no. 3, pp. 201–208, 2006.
  20. D. L. Parker, “Optimal short scan convolution reconstruction for fanbeam CT,” Medical Physics, vol. 9, no. 2, pp. 254–257, 1982.
  21. F. Natterer, The Mathematics of Computerized Tomography, Wiley-Teubner, New York, NY, USA, 1986.
  22. F. Noo, R. Clackdoyle, and J. D. Pack, “A two-step Hilbert transform method for 2D image reconstruction,” Physics in Medicine and Biology, vol. 49, no. 17, pp. 3903–3923, 2004.
  23. A. C. Kak and M. Slaney, Principles of Computerized Tomography Imaging, IEEE Press, Piscataway, NJ, USA, 1987.
  24. L. Zhu, J. Starman, and R. Fahrig, “The relationship between the T-FDK algorithm and the Hu-FDK algorithm,” in Proceedings of the 9th International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine, pp. 366–369, Lindau, Germany, July 2007.