Analysis of Cone-Beam Artifacts in off-Centered Circular CT for Four Reconstruction Methods

Cone-beam (CB) acquisition is increasingly used for truly three-dimensional X-ray computerized tomography (CT). However, tomographic reconstruction from data collected along a circular trajectory with the popular Feldkamp algorithm is known to produce the so-called CB artifacts. These artifacts result from the incompleteness of the source trajectory and the resulting missing data in the Radon space increasing with the distance to the plane containing the source orbit. In the context of the development of integrated PET/CT microscanners, we introduced a novel off-centered circular CT cone-beam geometry. We proposed a generalized Feldkamp formula (α-FDK) adapted to this geometry, but reconstructions suffer from increased CB artifacts. In this paper, we evaluate and compare four different reconstruction methods for correcting CB artifacts in off-centered geometry. We consider the α-FDK algorithm, the shift-variant FBP method derived from the T-FDK, an FBP method based on the Grangeat formula, and an iterative algebraic method (SART). The results show that the low contrast artifacts can be efficiently corrected by the shift-variant method and the SART method to achieve good quality images at the expense of increased computation time, but the geometrical deformations are still not compensated for by these techniques.


INTRODUCTION
Tomographic reconstruction has been a very active research field for over twenty years. Early works were rather focused on the inverse problem of tomography [1], the theoretical basis of which was founded on the Radon inversion formula, while later contributions have mostly been devoted to practical reconstruction methods for specific acquisition geometries. Tomographic reconstruction methods are generally classified into two classes of reconstruction algorithms (see, e.g., [2,3]): analytical and algebraic methods. The discretization of continuous inversion formulae leads to analytical methods while algebraic methods are based on the resolution of a linear system modeling the projection process. In algebraic reconstruction, the choice of the resolution technique results in a given algorithm (ART [4], SIRT [5], SART [6], etc.) characterized by an iteration formula which is completely independent of the acquisition geometry. In contrast, analytical methods are strongly dependent on the scanner geometry. Successive generations of scanners led to two-dimensional (2D) tomography with parallel-beam (PB) geometry, fan-beam (FB) geometry [7,8], and today 3D cone-beam (CB) tomography [9][10][11]. 3D tomography enables the reconstruction of larger (multislice) volumes with reduced acquisition duration and irradiation.
The trajectory of the X-ray source is an important feature in 3D CB tomography. In addition to simple circle, numerous trajectories such as circle plus line [12], two circles [13], set of lines [14], helix, and so forth, have been investigated but currently, most contributions deal with the circular and the helical trajectories [15,16]. The completeness condition on the source trajectory establishes whether or not exact reconstruction can be achieved from the projection data. According to the condition due to Kirillov [17] and formulated by Tuy [18] and Smith [19], a trajectory is complete if any plane crossing the support of the object intersects the vertex path at least once. No exact reconstruction can be expected if the source trajectory does not fulfill this condition owing to the incomplete projection data set.
In regular circular CB geometry, the most commonly used reconstruction algorithm is the FDK [20] which is a generalization to 3D CB of the standard filtered backprojection (FBP) algorithm. It is clear that the circular trajectory is incomplete since a plane parallel to the source path might intersect the support of the reconstructed volume while not crossing the source vertex. Reconstruction is therefore approximate except in the mid-plane (where it is equivalent to FB reconstruction). If reconstructions of good quality can be obtained with FDK for limited CB angles, artifacts rapidly increase with the distance to the mid-plane for a given focus length. Figure 1 represents the mid-plane and a vertical slice of a 3D Shepp-Logan phantom reconstructed with a total CB angle equal to 36 • . The reconstruction suffers from two kinds of CB artifacts: drop of low contrast intensity and geometrical deformations which appear at contrast jumps along the rotation axis (here: vertical axis). In order to reduce CB artifacts, different approaches have been reported in the past few years. Grass et al. [21] proposed the T-FDK algorithm which relies on changing the ramp filtering direction through CB to parallel fan-beam (PFB) rebinning. Yu et al. [22] have improved the reconstruction resolution via shift-variant filtering. Another class of methods is based on the Grangeat formula which relates the CB projections to the derivative of the 3D Radon transform [23]. The Radon transform reconstituted from the CB projection collected along a circular trajectory forms a torus instead of a sphere. This difference, referred to as the shadow zone, corresponds to the missing data. The reconstruction obtained with zero padding of the shadow was presumed to be equivalent to FDK and therefore suffered from CB artifacts. Lee et al. showed that filling in the shadow zone by interpolation could improve reconstruction [24]. Moreover, Hu made an important contribution providing a link between the Grangeat and the FDK formulae. He established that the FDK formula only deals with the inner part of the Radon torus and proposed a method which added the contribution of the Radon data on the torus shell [25]. Yang et al. completed this approach with information on the shadow zone [26]. More empirical approaches were suggested such as a 3D weighting of the projection data before the backprojection step [23] or error-correction-based methods [27][28][29] which associate analytical and algebraic approaches.
We recently introduced a novel acquisition geometry for CT data which enables simultaneous acquisition of PET and CT data for small animal [30]. This concept is based on photon counting pixel X-ray detectors expected to have a very high counting rate at X-rays energy (around 50 keV) while not stopping gamma rays (511 keV) [31]. This feature permits their placement inside a full micro-PET ring, in front of the PET crystals. However, this configuration imposes the X-ray source to be located outside the PET ring on an independently rotating system and focused on the center of the PET field of view (FOV). In this design, the X-ray source is off-centered with respect to the volume of interest (VOI) which results in the rotation of the source and the detector in two distinct planes. We first proposed a generalized formula called α-FDK [30] allowing for this topology, which comes down to standard FDK when applied to the simple circular geometry. However, the off-centered geometry is clearly exposed to CB artifacts since the entire volume of interest might be outside the exactly reconstructed mid-plane. The α-FDK performs the reconstruction of a shifted (or not) VOI defined by the angle α and partially compensates for the lowcontrast attenuation via the addition of an offset to obtain the correct value in the central plane. Nevertheless, the formula produces geometrical deformations and does not suppress the ramp of the reconstructed intensity. Reconstruction results obtained on simulations with the α-FDK therefore suffered from severe CB artifacts with large values of α and multislice volumes.
The purpose of this paper is to evaluate four different reconstruction techniques in our specific off-centered geometry in terms of CB artifacts reduction. We consider three analytical reconstruction methods, the α-FDK [30], the adaptation to our geometry of the shift-variant filtering (SVF) technique [22], and an FBP Grangeat-based (GB) formula [25]. We also evaluate an algebraic method, SART [6].
The paper is organized as follows. Section 2 gives a detailed description of the acquisition geometry. The four reconstruction methods are presented in Section 3. The reconstruction results obtained by computer simulations are presented and discussed in Section 4.

ACQUISITION GEOMETRY
The off-centered geometry parameterized by angle α is depicted in Figure 2. Let f (x, y, z) be a 3D attenuation function whose support is a sphere of radius r centered on the origin O of a Cartesian coordinate system. The plane (O, x, y) is called the central plane. The X-ray source is located at a S. Valton et al. distance R from the origin where R > r and rotates in the plane z = c where c = 0, called the mid-plane. The X-ray tube is oriented so that the central ray passes through O. Its position on the vertex path is parameterized by angle β. The acquisition geometry is characterized by α, defined as the angle between the central ray and the central plane. Detector plane is perpendicular to the central ray and faces the source at a distance D. For the sake of simplicity and without loss of generality, we will consider a virtual detector whose center is on the z axis, that is, D = R. Let (u, ν) be the coordinate system centered on the orthogonal projection of the source on the detector, where u refers to the lines of the detector (z constant) and ν refers to the columns. Let p α,β (u, ν) denote the 2D off-centered CB projection at angular position β in the off-centered geometry defined by angle α.

α-FDK
We derived a modified FDK expression adapted to the parameterization of the off-centered geometry. We recall the generalized inversion α-FDK formula [30] where the weight U is given by and (t , s , r ) is the coordinates system obtained by rotation of angle β around z axis followed by rotation of angle −α around x axis so that the s axis points toward the source.
Coordinates t , s , and r are given by The preweighting and filtering operations are given by where h is the ramp filter. This formula is equivalent to the standard FDK when α = 0. From a practical point of view, the α-FDK offers the same advantages as FDK for implementation: filtering is made on the lines of the detector, and each projection can be processed independently making the algorithm easily parallelizable.

T-FDK with shift-variant filtration
The algorithm proposed in [22] by Yu et al. for CB reconstruction is an FBP applied to parallel fan-beam data. The CB-to-PFB rebinning is performed via the Fourier space. The rebinned PFB projections are given by pr ϕ(u, ν) where angle ϕ satisfies The projections obtained after the CB-to-PFB rebinning operation lie on a curved virtual detector where the native coordinates are no longer equispaced. This step is thus followed by a vertical rebinning to obtain vertically equispaced pr ϕ (u, ve) projection data where Then, a similar horizontal rebinning is performed in [21] while Yu et al. alternatively propose to modify the 1D ramp filter. A shift-variant filtering of the nonrebinned data is applied instead of the ramp filtering of the equispaced data. Let ξ be the equispaced coordinate along the detector lines given by Then one obtains via a changing variable the following reconstruction formula: 4 International Journal of Biomedical Imaging This formula has two drawbacks compared to the standard FDK: (i) the CB-to-PFB rebinning requires all the CB projections which is memory and time consuming, (ii) the shift-variant filter kernel prevents from implementing the filtering in the Fourier space.

Radon-based method (RB)
The 3D Radon transform f (ρ, n) of the function f corresponds to its integral along the plane of normal vector n at a distance ρ from the origin: Hu showed in [25] that the function f RB reconstructed with the Grangeat formula based on the 3D Radon inversion could be expressed as the sum of three terms: The first term corresponds to the FDK formula, or the inversion of the 3D Radon data contained into the torus. The second term can be interpreted as the inversion of the 3D Radon data situated on the torus shell and was expressed as a filtered backprojection. The last term is a contribution from the shadow zone which is set to 0 in Hu's method. The first term is thus given by The second term is expressed as where the Radon data on the circle which describes the radon shell is obtained by where the preponderate projection p β is given by (4) with p β = p α=0,β . It is important to note that both the RB and the SVF methods are given for standard trajectory. We therefore performed reconstructions with these methods by considering a shifted VOI to account for the off-centered geometry without modifying the algorithms.

SART
Algebraic methods are based on a discrete modeling of the reconstruction problem. The 3D attenuation function is represented by an N-element vector F and for each β, the projections by an M-element vector P αβ obtained by the following matrices product: where R αβ is the M ×N matrix describing the projection process. The purpose is to find F, given projections P αβ for β describing a discrete set of angles. The solution F is then approached by iterative corrections of the volume F k . The algorithm we implemented is the SART (simultaneous algebraic reconstruction technique) and adopts the following scheme. For each cycle, (i) for each projection, the following is required: (1) estimation of the projection of the reconstructed volume F k , (2) difference between real projection and estimated projection, (3) normalization of the error, (4) backprojection of the error. This algorithm is summarized by the iteration formula where N is the number of voxels, M the number of pixels on the projections, and λ is a relaxation parameter. The SART algorithm was implemented on the basis of the α-FDK, with minimum modifications. The projection matrices used for the backprojection were used reversely to compute the reprojections. A bilinear interpolation on the detector was used both for back-and reprojection, and the relaxation parameter λ was set to 1. This implementation produces numerical artifacts which were smoothed by the application of a mean filter to the final reconstructed volume.

Simulation parameters
Simulations were performed with the 3D Shepp-Logan phantom [32]. The acquisition system was defined with the features given in Table 1 where all distances are given in centimeters.
We tested the reconstruction with six values of angle α between 0 and 0.5 radian which corresponds to 0 • . to 28.5 • . The resulting half-CB angle varies from 9.5 • to 50 • which is very large compared to current values in standard CB geometry.
The four methods developed in Section 3 were implemented and results were studied on simulations. The peakto-peak signal-to-noise ratio (PPSNR) is used for the quantitative evaluation of the reconstructions. It was computed  using the following expression: where MSE is the mean square error between the reconstructed and the reference phantom volume. In the following subsections, we present and discuss reconstruction results obtained with the three analytical methods and the SART algorithm after ten cycles.  Figure 4. As expected, the volumes reconstructed with α-FDK (Figure 3(a)) show two kinds of asymmetrical artifacts when α = 0: vertical geometrical deformations and drop of lowcontrast intensity. The vertical profiles displayed in Figure 4 include the result obtained with the standard FDK and reveal an offset between the two formulae. The RB algorithm (Figure 3(b)) partially compensates for the intensity drop but since no offset is used, the reconstructed intensity stays below the result of the α-FDK (Figure 4). In Figure 3, the result with α = 28.5 • appears dark because the same grey level window was used for comparison between the different reconstructions. The SVF algorithm (Figure 3(c)) better compensates for the intensity drop but the vertical elongation remains. The results obtained with standard SART after 10 cycles displayed in Figure 3(d) show that the intensity drop is only partially compensated for while the elongation is slightly corrected. The PPSNRs reported in Table 2 indicate that, quantitatively, the best reconstruction is obtained with the algebraic algorithm.

Discussion
In standard circular trajectory, CB artifacts correction is already a challenge since those artifacts come from the illconditioning of the inverse tomographic problem. In offcentered geometry, this problem is worse since the CB angle is more important and the VOI is asymmetrical with respect to the exactly reconstructed mid-plane.
The α-FDK provides a practical formula adapted to the parameterization by the angle α of the off-centered geometry while preserving the rapidity and simplicity of the standard formulation. It compensates for the low-contrast drop by adding an offset to the reconstructed value. However, since it does not suppress the increase of the drop with the distance to the mid-plane, this compensation is fully valuable only for a small volume around the central plane of the VOI.
The RB method reduces the CB artifacts in standard circular geometry but the correction is weak. In off-centered geometry, the α-FDK outperforms the RB method when α is large. In addition, the method presented by Hu takes 20% more time than the α-FDK.
Concerning the SVF method, the simulation results presented above show that the intensity drop is well corrected by this method which, however, does not deal with the geometrical deformations. We can notice in Table 2 that the PP-SNR is slightly inferior to that obtained with the α-FDK until α = 28.5 • . We presume that this difference is due to the numerical errors introduced by the different rebinning operations. Given the supplementary rebinning step and the shiftvariant filtering, this algorithm is more time and memory consuming than the α-FDK. Besides, the rebinning step prevents from starting the reconstruction before the end of the acquisitions. To conclude, reconstructions obtained with the SVF algorithm are better than those obtained with the α-FDK, but the computation time is twice as long. The choice of the optimal reconstruction method thus depends on the constraints of the application. We have seen that the reconstruction obtained with the α-FDK for a volume made of a limited number of slices is acceptable. Therefore, in the perspective of a bimodal scanner, it might be interesting to prefer this simple algorithm if the axial FOV of microPET scanners is limited to few slices, depending on the number of detector rings.
Our implementation of the SART algorithm produces numerical noise resulting from the model of forward projection and the value of the relaxation factor. This highfrequency noise appears very clearly in the picture if no mean filtering is applied because of the narrow grey level display window. In contrast, it has minor effects on the PPSNR thanks to its low amplitude and MSE. The high PPSNR compared to the other methods is related to the algebraic nature of the SART which minimizes a numerical error, namely, the difference between the real and the estimated projection. The profile plot in Figure 4 shows that the correction of intensity drop is less efficient than with the SVF method, but a close look at the right end of the graph points out that the elongation is slightly corrected by this algebraic reconstruction. Nevertheless, we think that this minor improvement is not worth the additional time expense needed to perform 10 cycles of the algorithm.
To our knowledge, no correction method compensates for the geometrical deformations due to the large CB angle. Actually, no clear mention of this kind of artifacts is made in the literature dealing with CB artifact since they clearly appear only with very large values of the CB angle. We believe that this aspect needs to be addressed in future works concerning CB artifacts reduction and that the correction of the elongation should, at least partially, compensate for the lowcontrast drop.

CONCLUSION
In this paper, a comparison between the results obtained with four different reconstruction techniques on data simulated with an off-centered CB circular acquisition geometry was presented. Since this geometry increases the proportion of missing data compared to the standard circular trajectory, the CB artifacts correction methods developed for the standard case need to be reinvestigated. On the one hand, the α-FDK reconstruction algorithm that we derived in a previous work produces strong CB artifacts, but a good correction is obtained in the planes closed to the VOI central plane, without additional time expense. On the other hand, the intensity drop is corrected by the RB algorithm, but not sufficiently enough, by the SART and even better by the SVF algorithm in return for a certain time expense. The preferable method therefore depends on the application, in terms of computation time and accuracy needs. Generally speaking, the SVF is preferable to the SART owing to the expensive computation time of the algebraic method, and to the α-FDK thanks to the attenuation correction. For limited number of slices and important value of α, the α-FDK gives better results than the RB method and can be used instead of SVF for its rapidity. Nonetheless, the geometrical deformations could not be efficiently corrected. Further research is therefore necessary to address this complicated problem inherent to CT.