Exact Weighted-FBP Algorithm for Three-Orthogonal-Circular Scanning Reconstruction

Recently, 3D image fusion reconstruction using a FDK algorithm along three-orthogonal circular isocentric orbits has been proposed. On the other hand, we know that 3D image reconstruction based on three-orthogonal circular isocentric orbits is sufficient in the sense of Tuy data sufficiency condition. Therefore the datum obtained from three-orthogonal circular isocentric orbits can derive an exact reconstruction algorithm. In this paper, an exact weighted-FBP algorithm with three-orthogonal circular isocentric orbits is derived by means of Katsevich's equations of filtering lines based on a circle trajectory and a modified weighted form of Tuy's reconstruction scheme.


Introduction
The cone-beam scanning configuration with a circular trajectory remains one of the most popular scanning configuration and has been widely employed for data acquisition in 3-D X-ray computed tomography (CT), because it allows for operating at a high rotating speed due to its symmetry, avoiding the need to axially translate the patient, such as in helical or step-and-shoot CT [1,2]. Unfortunately, a circular trajectory does not satisfy Tuy's data sufficiency condition [3] which requires every plane that passes through the reconstruction region (ROI) must also cut the trajectory at least once and therefore only approximate reconstructions are possible. In the pursuit of the data sufficiency condition, various scanning trajectories have been proposed, such as circle and line [4,5], circle plus arc [6,7], double orthogonal circles [8,9] and dual ellipses [10].
Recently, the FDK algorithm [11] has been implemented for a cone-beam vertex trajectory consisting of three-orthogonal isocentric circles [12]. On the other hand, we know that 3D image reconstruction based on three-orthogonal circular isocentric orbits is sufficient in the sense of Tuy data sufficiency condition. Therefore the datum obtained from three-orthogonal circular isocentric orbits can derive an exact reconstruction algorithm. In this paper, another theoretically exact and general inversion formula which was proposed in Katsevich [13] will be implemented for three-orthogonal circular isocentric orbits. Compared with the aforementioned algorithms, the distinctive features of Katsevich's algorithm can be summarized as the choice of a more general weight and a novel way to dealing with discontinuous in weight. Furthermore, Katsevich's algorithm has been implemented for various scanning trajectories, such as circle and line [14,15], circle plus arc [16], several circular segment [17], two orthogonal circles [13] and helix [18,19].
In contrast to a singular circular trajectory, there are several advantages in using weighted reconstruction with the three-orthogonal-circular trajectory. First, for the reconstruction noise due to the algorithm is 'trajectory' dependent. Using three-orthogonal geometry with weighting function, the noise is suppressed. Second, the three-orthogonal setting yields better image quality in comparison with classical scheme for larger cone-angles. Last, the scan method for three-orthogonal geometry can be implemented with minimal requirements and can provide a better 3D reconstruction for small animals or objects. At the same time, compared to the two-circular trajectory, artifacts of the boundaries of the ROI where artifacts are likely to appear are suppressed by the additional filtering lines and according weighting functions.
In our implementation, a new weighting scheme is developed so that all measurements are used by accurately averaging over multiply measured projections of the three-orthogonal-circular scanning. Moreover, this new weighted reconstruction is differs from the weighted FDK [11] reconstruction with the three-orthogonal-circular trajectory [12] considerably, in that it is a theoretically exact formulation of a general shift invariant filtered back-projection (FBP) reconstruction framework. In this paper, we show the derivation of our cone-beam FBP reconstruction algorithm starting with Tuy's inversion formula. In other words, we discuss how to derive a FBP cone-beam reconstruction formula from the Tuy's classical inversion formula. To construct the reconstruction formula, several operators are needed. First, using some properties of the cone-beam transform, the classical Tuy's inversion formula can be written as (the first-order derivative of the Radon transform in R 3 ) Grangeat's inversion formula [20] in a weighted form. Then, using the methodology in [13], we can get the final inversion formula. The features of our inversion formula can be summarized as follow: using some properties of the cone-beam transform to derive the inversion formula, a proper choice of the weighting function, deriving equations of the filtering lines and describing geometric properties of filtering lines in the planar detector plane.
The organization of this paper is as follows. In Section 2, we derive Katsevich's inversion formulation by Tuy's weighted form. In Section 3, we derive equations of the filtering planes and filtering lines. In Section 4, we show the geometric properties of filtering lines in the planar detector plane and according weighting functions.

Weighting Function for Tuy's Formula
Let a trajectory C be a differentiable curve in R 3 described by y(s), s ∈ R. The object density function is f (x), where x is a vector in the (y 1 , y 2 , y 3 ) coordinate system, and f (x) is an infinitely differentiable real integrable function with a compact support Ω ⊂ R 3 \C. The modified cone-beam projection of f (x) along the direction of α/ α at the focal point location y(s) is defined as: Then, let Df be extended from C × S 2 to R 3 × R 3 as: First consider the three-dimensional Fourier transform of D y f (z) for a fixed f (x) as: It can be easily verified that: In the following, we will show how to derive a new reconstruction formula from Tuy's reconstruction scheme in a weighted form.

Theorem 1. [Tuy's weighted form]
Let C be a curve in R 3 parameterized by a piecewise differentiable function y(s), s ∈ R. For a fixed x ∈ Ω, we suppose that there exists a weighting function n x (s, σ) : R × S 2 → R such that n x (s, σ) is integrable with respect to the second variable σ ∈ S 2 for each s ∈ R, and the set: is non-empty and finite for almost all σ ∈ S 2 . Then: provided that n x (s, σ) fulfilled the completeness condition: This inversion formula also tells us Tuy's data sufficiency conditions for an accurate reconstruction of a ROI from cone-beam projections. These conditions are: σ · y(s) = σ · x and σ · y (s) = 0. In the following, we will show how to derive a FBP reconstruction formula from Tuy's formula and show how Tuy's data sufficiency and nonended condition can be relaxed. In a first step we try to use Equation (1) to simplify Equation (4) as follows: where r is the spherical radical coordinate, θ is a three dimensional unit vector such that z = rθ and θ = 1. Next we introduce the quantity: In Equation (7) the first term is real and add, and the second term is imaginary and even if we take into account the fact that the factor y (s) · σ in the inversion formula is also odd in σ, we conclude that only the first term contributes in Tuy's inversion formula. The second term will vanish when we perform the integration over σ. Thus we have: Thus, we rederived the weighted form of Tuy's inversion formula as follows: Therefore, there is a derivative of delta function in the integration over θ. We can perform integration by part for the integration over θ in Equation (10). We thus obtain the directional derivative along the direction σ of cone-beam data Df (y(q), θ) with respect to unit vector θ. Therefore Equation (10) can be written as follows: Using Grangeat's formula and change of variables ρ → q defined by ρ = y(q) · σ, we obtain: wheref (ρ, σ) is the 3D Radon transform of f . Modifications from Tuy's original inversion (6) to Equation(12) with change of variables have been important progress. Accounting for the numerical implementation of formula Equation (11) is still complicated. We can change the two integrals over the unit vector θ = (− sin γ cos ϕ, sin γ cos ϕ, cos γ) and σ = (cos ϕ, sin ϕ, 0) into three single integrals over parameters s, ϕ and γ. Using the methodology in [13], we get: where β(s, x) = (x − y(s))/ x − y(s) = (0, 0, 1),e x (s, ϕ) = β(s, x) × σ x (s, ϕ), denote The values of w m (s, x) are determined by the definition of the weighting function and signum function near the discontinuous points ϕ m 's. Furthermore, in the FBP reconstruction scheme, the unit vector σ x (s, ϕ m ) is the normal vector of filtering plane and the vector e x (s, ϕ m ) denotes the direction of filtering lines which pass through the object point x and the source point y(s). After substituting Equation (14) into formula Equation (13), we obtain: , cos rβ(s, x) + sin re x (s, ϕ m ))| q=s dr sin r ds. (16) The final inversion formula is independent of the specific geometrical shape of the image object. Rather, it is determined by the scanning geometry. Thus, we need study the specific features of filtering lines and weighting functions. In the following, the inversion formula will be implemented for a conebeam vertex trajectory consisting of three-orthogonal to each other circular. Furthermore, the filtering lines and the weighting functions will be given in detail.

Equation of Filtering Lines for the Three-Orthogonal Circular Scanning
We propose the use of three independent, orthogonal to each other, circular isocentic scan-paths for X-ray projection acquisitions. Let a trajectory C = C 1 ∪ C 2 ∪ C 3 , as shown in Figure 1. where: Let f (x) is the density function to be constructed. Assume that the function is smooth and vanishes outside the ball: where r is the radius of the subject ball and R is the radius of the scanning circular. We assume that the physical detector is a planar-detector, which is denoted as DP (s), where s is the parameter of source point. Furthermore, let the planar-detector is at distance 2R from the source, as shown in Figure 1.
Denote by (u, v) the horizontal and vertical coordinates of the planar-detector plane, where u is parallel to the tangent vector of source, v is the normal vector of scanning trajectory and the origin is at the projection of the source point y(s) onto the detector plane. Therefore, in the (y 1 , y 2 , y 3 ) coordinate system, any point on the detector plane can be characterized completely by use of (u, v). It can be shown that: (21) Figure 2. The filtering plane is shown, which is tangent to one of circular trajectories and is perpendicular to σ x . Figure 2 shows a filtering plane R(x, σ)through the y(s) and is tangent to other trajectory at point y(λ). Without loss of generality we consider the source y(s) ∈ C 1 (the case y(s) ∈ C 2,3 is treated in a similar fashion). The point of tangency is given by either y(λ) = (R cos λ, 0, R sin λ) or y(λ) = (0, R cos λ, R sin λ). In the first case, the normal vector to the tangent plane is given by: s cos λ, 1 − cos s cos λ, sin s sin λ).
The equation of tangent plane, which is a filter plane, is given by: Filtering lines on the detector are obtained by intersecting the detector plane and filtering planes. Substituting Equation (20) and Equation (23) into Equation (24) and solving for v, we obtain the equation of the filtering line on the planar-detector: In the second case, the normal vector to the tangent plane is given by: cos s cos λ, cos s cos λ, cos s sin λ).
The equation of the filtering line on the planar-detector is given by:

Weighting Functions of Filtering Lines
In the following, we define adequate filtering lines for the three-orthogonal scanning case. In order to define the different filtering lines, which are necessary to perform the three-orthogonal circular reconstruction, we first, introduce certain curves, which separate the planar-detector into different regions. As an example consider the X-ray source moving on the trajectory C 1 . We project stereographically the trajectories C 2 and C 3 onto the detector plane as shown in Figure 3. Let P C 2 and P C 3 denote these projections. The line P C 1 is the projection of trajectory C 1 . Supp f (x) is supposed to be inside a ball centered as the origin and with sufficiently small radius r < R, so that the projection of ROI has a circular shape, as shown by the shaded region in figure 3. The curves P C 2 and P C 3 split the entire ROI into three sub-ROIs: P R m , m = 1, 2, 3. The object point x is projected into the area P R m , m = 1, 2, 3. Let x denotes this projection.
Ifx ∈ P R 1 , L(ϕ) is the projection of the plane through x with normal σ x (s, ϕ). As ϕ increase, σ x (s, ϕ) ∈ β ⊥ (s, x) rotates in the clockwise direction on DP (s). From definition of k x (s, ϕ) with equation (14), k is discontinuous if L(ϕ) is parallel to y (s) or is tangent to C 2 . In the first case, such a plane has six points of intersection with C = C 1 ∪ C 2 ∪ C 3 , so k = 1/6 on the side of the jump, where y (s) · σ > 0 and k = −1/6 on the other side. The filtering line is parallel to y (s), which is denoted as L 1 (see the dot-dashed lines in Figure 3). The corresponding weighting function of filtering line L 1 is w x (s, ϕ 1 ) = 1/3. In the second case, the number of intersection changes form four to six . So the filtering line is tangent to C 2 , which is denoted as L 2 ( see the dot-dashed lines in Figure 3). For the signum function in Equation (14) is unchagened, the corresponding weighting function of filtering line L 2 is w x (s, ϕ 2 ) = 1/12. Ifx ∈ P R 2 , k is discontinuous only L(ϕ) is parallel to y (s). So the filtering line L 1 is parallel to y (s). The corresponding weight function is w x (s, ϕ 1 ) = 1/3. Ifx ∈ P R 3 , k is discontinuous if L(ϕ) is parallel to y (s) or is tangent to C 3 . In the first case, such a plane has six points of intersection with trajectory, so k = 1/6 on the side of the jump, where y (s) · e > 0 and k = −1/6 on the other side. The filter line is parallel to y (s), which is denoted as L 1 (as shown in Figure 3).The corresponding weighting function of filtering line L 1 is w x (s, ϕ 1 ) = 1/3. In the second case, the number of intersection changes form four to six. So the filtering line is tangent to C 3 , which is denoted as L 3 ( see the dot-dashed lines in Figure 3). The corresponding weighting function of filtering line L 3 is w x (s, ϕ 3 ) = 1/12. Figure 3 summarizes the above information.

Conclusions
An exact shift invariant filtered back-projection (FBP) reconstruction algorithm for a cone-beam vertex trajectory consisting of three-orthogonal to each other circular was derived from Tuy's inversion formula. There are several major modifications to Tuy's formula. The first modification is not using the Fourier transform with the projection function to deduce the specific inversion formula, but instead using the inversion Radon transform. Second, we started with a Tuy-like inversion scheme. Weighting the redundant data and rewriting the weighted summation into an integral along the source trajectory resulted in a shift-invariant FBP reconstruction formula. Last, the "nontangential" condition in Tuy's original data sufficiency conditions is relaxed and the nonended condition is further relaxed. In implementing of the inversion formula for a cone-beam vertex trajectory consisting of three-orthogonal to each other circular, the concrete forms of the filtering lines and according weighting function are the most important and difficult segments. In this paper, we obtained above two segments base on analysis about the geometric properties of projections with the three-orthogonal circular trajectory and the radon planes which pass through the reconstruction point.