An Approximate Cone Beam Reconstruction Algorithm for Gantry-Tilted CT Using Tangential Filtering

FDK algorithm is a well-known 3D (three-dimensional) approximate algorithm for CT (computed tomography) image reconstruction and is also known to suffer from considerable artifacts when the scanning cone angle is large. Recently, it has been improved by performing the ramp filtering along the tangential direction of the X-ray source helix for dealing with the large cone angle problem. In this paper, we present an FDK-type approximate reconstruction algorithm for gantry-tilted CT imaging. The proposed method improves the image reconstruction by filtering the projection data along a proper direction which is determined by CT parameters and gantry-tilted angle. As a result, the proposed algorithm for gantry-tilted CT reconstruction can provide more scanning flexibilities in clinical CT scanning and is efficient in computation. The performance of the proposed algorithm is evaluated with turbell clock phantom and thorax phantom and compared with FDK algorithm and a popular 2D (two-dimensional) approximate algorithm. The results show that the proposed algorithm can achieve better image quality for gantry-tilted CT image reconstruction.


INTRODUCTION
In some applications of clinical CT scanning, it is required that the CT gantry be tilted. For example, in order to avoid exposure of eyes to X-rays, gantry is tilted during the head scanning procedure. To meet such special requirements, a number of algorithms have been developed for gantry-tilted helical MSCT (multislice computed tomography) image reconstruction [1][2][3][4]. Among these algorithms, Kacherließ et al. [2] developed a gantry-tilted reconstruction algorithm based on the earlier developed 2D algorithm ASSR [5]. Hein et al. [3] developed a gantry-tilted reconstruction algorithm based on the 3D FDK algorithm [6,7]. In the gantry-tilted FDK algorithm in [3], the reconstruction plane is perpendicular to the rotating axis of the CT scanning. Thus the scanning cone angle increases with the pitch value and slice number which can lead to unavoidable artifacts. Recently, Noo et al. [4] developed a general framework which can be applied to the gantry-tilted CT for exact [8,9] and approximate image reconstruction.
The FDK-type algorithms for the approximate reconstruction using the projection data filtering in the horizontal direction suffer from considerable artifacts due to the large scanning cone angle for both normal and gantry-tilted helical MSCT imaging. Recently, some of the improvements have been made for the FDK-type algorithms. The tilted plane technique is combined with the FDK algorithm to reduce the cone angle and further reduce artifacts caused by the large cone angle [10]. Several methods reduce the artifacts by filtering projection data along the tangential direction of the helix [11][12][13]. Such a technique was earlier proposed by Yan and Leahy [11] and further improved by Sourbelle and Kalender [12] for short scan FDK-type algorithms. The improved FDK-type algorithm can achieve better image quality than that of the conventional FDK algorithm.
Motivated by the observation that artifacts for FDK-type algorithms can be effectively reduced by filtering projection data along the helix tangential direction, this paper extends the existing tangential filtering technique to present a 3D FDK-type approximation algorithm for gantry-tilted helical MSCT. Taking into account the gantry-tilted geometry, we provide a general formula for gantry-tilted CT reconstruction for different gantry-tilted angles. As a special case, the reconstruction formula reduces to the standard tangential filtering for conventional CT with zero gantry-tilted angle.
To deal with the complicated gantry-tilted geometry and large scanning cone angle, our proposed algorithm first 2 International Journal of Biomedical Imaging Centre plane applies the ramp filtering to reconstruct a sequence of image planes and these planes are perpendicular to the rotating axis. The final horizontal image planes are then obtained by interpolating the tilted image planes. We present simulation results to show that the proposed gantry-tilted algorithm can provide considerably improved imaging for large scanning cone angle.
The key technique of the proposed algorithm is to reconstruct images on planes perpendicular to the rotation axis using the ramp filtering along the helix tangential direction. This is essentially different from the tilted plane reconstruction technique, such as that in [10], which is based on geometrically optimized reconstruction plane for reducing the effective cone angle.
It is noted that the tangential filtering technique can be implemented as a special case of the general framework for gantry-tilted CT proposed by Noo et al. [4]. Following its procedure, the gantry-tilted geometry is first transformed into a conventional CT scanning geometry with zero gantrytilted angle and a projection data set of the transformed scanning geometry is formed. Then the tangential filtering technique can be applied on the data set for the image reconstruction. The procedure of the proposed algorithm in this paper is different from Noo et al. 's general framework in that it computes the helix tangential direction and applies the filtering directly without rebinning the projection data set for the transformed scanning geometry. This leads to more efficient computing and image reconstruction.
The rest of this paper is organized as follows. Section 2 presents the geometric scheme of the gantry-tilted CT. Section 3 is on the projection data set formation for tangential filtering. The proposed approximate reconstruction algorithm presented in Sections 4 and 5 is on simulation and evaluation of the proposed algorithm.

GEOMETRY FOR GANTRY-TILTED HELICAL MSCT
The helical MSCT scanning set up consists of an X-ray source and a detector array forming a source-detector framework.
The geometry of the source-detector framework in a global cartesian coordinate system x − y − z is shown in Figure 1(a), where C s denotes the X-ray source, the detector array is a rectangular surface with a geometric centre C d . The sourcedetector framework defines a rotating cartesian coordinate system s − l − u with the X-ray source C s as the origin, the straight line C s − C d being the s axis, the l axis being parallel to the horizontal lines, and the u axis being parallel to the vertical lines of the rectangular detector surface.
The s − l plane containing C s and C d is called the centre plane and its geometry is shown in Figure 1(b), where β is the projection angle, the distances from the rotating centre point o to the detector and the X-ray source are R d and R f , respectively. In the gantry-tilted scanning process, the sourcedetector framework rotates around the axis u which is parallel to the u axis and intersects the z axis at o as shown in Figure 1(a). The gantry-tilted angle is represented by μ between the z and u axes. And the projection angle β is defined in the s − l plane as shown in Figure 1 For gantry-tilted CT, the X-ray source trajectory can be considered as a combination of two movements: the X-ray source rotates on the tilted circular trajectory with a radius R f and the centre of the circle moves straightforward along the table feeding direction. As a result, the X-ray source trajectory is the sum of the circular rotation trajectory on the centre plane and the linear translation in the z direction in the following: Ming Yan et al. where β 0 is the initial projection angle, z 0 is the initial z position of the centre of the gantry, p is the pitch value of the helical cone beam scanning, S is the slice thickness, and M is the number of detector slice.
With the X-ray source trajectory given in (1), we take the derivative with respect to β to obtain its tangential direction at C s in the global coordinate system represented in the following form: where h = pSM/2π. The X-ray source trajectory and its tangential direction at C s are shown in Figure 2. It follows from the geometric relation between the global x − y − z and the rotating s − l − u coordinate system, as shown in Figures 1(a) and 1(b), that the transformation of the tangential direction (x t , y t , z t ) T as given in (2) to the rotating s − l − u coordinate system can be obtained by first rotating this vector by an angle μ around the x axis followed a rotation angle of β − π around the u axis. This results in the trajectory tangential direction in the rotating s − l − u coordinate system, denoted by (s t , l t , u t ) T , in the following: To implement the tangential filtering for the gantry-tilted scanning projection data, we further use the tangential direction (s t , l t , u t ) T to set up a tangential cartesian coordinate system ξ − η − ζ. The origin of this coordinate system is the X-ray source C s , the η axis is defined by the tangential direction of the X-ray source trajectory at C s , as given in (3), the ζ axis is parallel to the detector plane, and the ξ − ζ plane is orthogonal to the η axis.
Using the tangential direction (s t , l t , u t ) T in (3), we introduce two angular quantities σ and γ as follows: Following from its definition, the ξ −η−ζ coordinate system can be obtained by first rotating the s − l − u coordinate system around the u axis by an angle σ to obtain an intermediate coordinate system ξ − η − u followed by rotating the ξ − η − u coordinate system around the ξ axis by an angle γ. Such coordinate rotations and transformation are illustrated in Figure 3. It follows that the transformation between the s − l − u and ξ − η − ζ coordinate systems is given by where the transformation matrix T is

DATA REFORMATION FOR TANGENTIAL FILTERING
In conventional FDK-type algorithms, the ramp filtering is performed on detector rows and the filtering direction is parallel to the l axis on the detector surface. Motivated by the tangential filtering technique for CT without gantry tilting, we propose in this paper that the ramp filtering is performed in the tangential direction of the X-ray source trajectory of the gantry-tilted CT as shown in Figure 2. For this purpose, the projection data are reformed such that rows of the reformed projection data set are parallel to the tangential direction of the X-ray source trajectory. The reformed projection data set enables the tangential filtering technique for the conventional CT reconstruction being applied. The reformed data set is built up by introducing a virtual detector array surface in the tangential coordinate system ξ − η − ζ and transforming the projection data from the real detector array surface in the rotating s − l − u coordinate system to the virtual detector array surface. The virtual detector array surface is placed on the plane ξ = (R f + R d ), as shown in Figure 4, with the rows being parallel to the η axis.
As sown in Figure 4, each projection ray radiates on the real detector array at the point ((R f +R d ), l p , u p ) T in the s−l− u coordinate also radiates on a point ((R f + R d ), η p , ζ p ) T on the virtual detector array. Let the ξ − η − ζ coordinate of the point ((R f + R d ), l p , u p ) T on the real detector array surface be (ξ p ,η p ,ζ p ) T . The transformation from ((R f + R d ), l p , u p ) T to (ξ p ,η p ,ζ p ) T is determined by the coordinate transformation matrix T in (7), that is, In view of Figure 4, the coordinate (η p , ζ p ) T of the virtual detector cell can be expressed, in terms of that of the real detector cell in the ξ − η − ζ coordinate system, as This, together with (8), can determine the coordinate of the virtual detector cell (η p , ζ p ) T in the tangential coordinate system ξ − η − ζ from a given real detector cell at ((R f + R d ), l p , u p ) T in the rotating coordinate system. Let the projection datum collected from the real detector cell at ((R f + R d ), l p , u p ) T in the s − l − u coordinate at the projection angle β be denoted by D(β, l p , u p ) and the corresponding projection datum on the virtual detector cell at ((R f + R d ), η p , ζ p ) T in the ξ − η − ζ coordinate by the same projection ray be denoted by D v (β, η p , ζ p ). Since D(β, l p , u p ) and D v (β, η p , ζ p ) are due to the same projection ray, we have Actual detector Virtual detector Figure 4: Illustration of the relationship between (ζ p ,η p ,ξ p ) and This can be used to obtain the reformed projection data set following from the coordinate transform from each real detector cell at ((R f + R d ), l p , u p ) T in the rotating coordinate system to the corresponding virtual detector cell at ((R f + R d ), η p , ζ p ) T in the tangential coordinate system.

RECONSTRUCTION AND INTERPOLATION
Because the gantry is tilted and final images should be on horizontal planes in the global x − y − z coordinate system, we first reconstruct images on a sequences of intermediate planes parallel to the centre plane and then interpolate them to obtain the horizontal images. The introduction of the intermediate planes effectively reduces the projection cone angle in comparison with constructing directly the horizontal planes. Let ψ = {P i : i = 1, 2, 3, . . . , n} be the set of the intermediate planes and let the intersection of each tilted plane P i and the z-axis be at point o i = (0, 0, z i ) T , with z i = z 0 + iΔz, i = 1, 2, 3, . . . , n, and Δz > 0. The projection angle for the X-ray source being at z i in the z direction is And the equation for the intermediate plane P i is The ramp filtering is along the η direction which is defined as the tangential direction of the source trajectory. Let g(·) denote the ramp filter. Applying it to the reformed data set on the virtual detector array in the tangential direction of the source trajectory and performing standard cone beam weighting yields  Then the reconstruction formula is derived from Yan and Leahy's [11] paper for tangential filtering reconstruction.
where · denotes Euclidean norm and β m is the the length of the projection angle interval. Given the reconstructed images of the tilted planes P i , for i = 1, 2, 3, . . . , n, and a point (x, y, z) T on a plane parallel to the x − y plane where the image is to be obtained by interpolation, the z-positions of the intermediate planes P i , i = 1, 2, 3, . . . , n, at (x, y) T and denoted by z pi (x, y) can be determined by (12). Thus there exist two points (x, y, z Pj ) T and (x, y, z Pk ) T on two tilted planes P j and P k , respectively, at the upper and lower sides of (x, y, z) T , respectively, which are closest to (x, y, z) T . Using the obtained attenuation functions f j (x, y, z Pj ) and f k (x, y, z Pk ) for planes P j and P k , respectively, the interpolated attenuation function can be obtained using the following interpolation formula: We can now summarise the proposed reconstruction procedure in the following.
(1) Determine the z position z i , centre of the optimal plane P i , its centre projection angle β i , and obtain the reconstruction plane using (12). For half-scan, β m = π + θ fan . (2) Obtain the reformed projection data set on the virtual detector array using (9).  Table 1. Results obtained from the simulation of Turbell clock phantom are shown in Figure 5, the left column is reconstructed images and the right column is the centre vertical line for reconstructed images. In Figure 5, Figure 5(a) is the image reconstructed by our proposed algorithm, Figure 5(c) is reconstructed by the gantrytilted FDK algorithm and Figure 5(e) shows the image reconstructed by the gantry-tilted ASSR. The results demonstrate that the image reconstructed by our proposed algorithm contains fewer artifacts than that of the other two existing algorithms for gantry-tilted CT.
In the simulation of Turbell clock phantom, the centre vertical lines (y = 0 mm) for the reconstructed images are also constructed and displayed by the proposed algorithm ( Figure 5(b)), gantry-tilted FDK algorithm ( Figure 5(d)), and gantry-tilted ASSR ( Figure 5(f)), repectively. In Figure 5, dashed lines represent the original image and solid lines represent reconstructed profiles. The results show that the profile of the proposed algorithm as shown in Figure 5(b) has smaller variance than the other two results.
A thorax phantom is further simulated which is designed by referring to human thorax consisting of many important organs and often scanned in CT examination. These organs include lungs, heart, aorta, ribs, spine, sternum, and shoulders. The phantom definitions are obtained from a 6 International Journal of Biomedical Imaging  world phantom database FORBILD. The simulation results are shown in Figure 6, where Figure 6(a) is the original thorax phantom, Figure 6(b) shows the image reconstructed by the proposed algorithm, Figure 6(c) is the image reconstructed by the gantry-tilted FDK algorithm, and Figure 6(d) is reconstructed by the gantry-tilted ASSR algorithm. The right-hand column shows zoomed reconstruction images. It is shown that there are obvious artifacts around ribs in the images reconstructed by conventional FDK algorithm and ASSR algorithm. In contrast, the proposed algorithm provides better image quality and more accurate reconstruction.

CONCLUSION
This paper presents an approximate algorithm for gantrytilted helical MSCT image reconstruction. It is based on