Next Article in Journal
A New and Inexpensive Pyranometer for the Visible Spectral Range
Previous Article in Journal
Development of a Prototype Miniature Silicon Microgyroscope
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

College of Computer, Sichuan University, 610065, Chengdu, P. R. China
*
Author to whom correspondence should be addressed.
Submission received: 28 April 2009 / Revised: 21 May 2009 / Accepted: 12 June 2009 / Published: 12 June 2009
(This article belongs to the Section Chemical Sensors)

Abstract

:
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.

1. 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 ℝ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.

2. Weighting Function for Tuy's Formula

Let a trajectory C be a differentiable curve in ℝ3 described by y(s), s ∈ ℝ. The object density function is f(x), where x is a vector in the (y1, y2, y3) coordinate system, and f(x) is an infinitely differentiable real integrable function with a compact support Ω ⊂ ℝ3\C. The modified cone-beam projection of f(x) along the direction of α/‖α‖ at the focal point location y(s) is defined as:
D f ( y ( s ) , α ) = 0 f ( y + t α ) d t = 1 | | α | | D f ( y ( s ) , α | | α | | ) , ( y ( s ) , α ) C × S 2 .
Then, let Df be extended from C × Sensors 09 04606i12 to ℝ3 × ℝ3 as:
D y f ( z ) = D f ( y , z ) , y , z 3 .
First consider the three-dimensional Fourier transform of Dyf(z) for a fixed f(x) as:
D y f ( σ ) = 3 D f ( y , z ) e 2 π i z σ d z .
It can be easily verified that:
D y f ( t σ ) = 1 t 2 D y f ( σ ) , for t > 0 .
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 ℝ3 parameterized by a piecewise differentiable function y(s), s ∈ ℝ. For a fixed x ∈ Ω, we suppose that there exists a weighting function nx(s, σ) : ℝ × Sensors 09 04606i12 → ℝ such that nx(s, σ) is integrable with respect to the second variable σ Sensors 09 04606i12 for each s ∈ ℝ, and the set:
R ( x , σ ) = { s | x σ = y ( s ) σ , y ( s ) σ 0 }
is non-empty and finite for almost all σ Sensors 09 04606i12. Then:
f ( x ) = S 2 s R ( x , σ ) n x ( s , σ ) 2 π i y ( s ) σ q D y ( q ) f ( σ ) | q = s d σ
provided that nx(s, σ) fulfilled the completeness condition:
s R ( x , σ ) n x ( s , σ ) = 1 , a . e . i n σ S 2 .
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:
i D y f ( σ ) = i R 3 D f ( y ( q ) , z ) e 2 π i z σ d z = i S 2 0 1 r D f ( y ( q ) , θ ) e 2 π i z σ r 2 d r d θ ,
where r is the spherical radical coordinate, θ is a three dimensional unit vector such that z = and ‖θ‖ = 1. Next we introduce the quantity:
i D y f ( σ ) = 1 2 S 2 D f ( y ( q ) , θ ) d θ i r e 2 π i r θ σ d r + 1 2 S 2 D f ( y ( q ) , θ ) d θ i | r | e 2 π i r θ σ d r .
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:
i D y f ( σ ) = 1 4 π S 2 D f ( y ( q ) , θ ) δ ( θ σ ) d θ .
Thus, we rederived the weighted form of Tuy's inversion formula as follows:
f ( x ) = 1 8 π 2 S 2 s R ( x , σ ) n x ( s , σ ) y ( s ) σ S 2 q D f ( y ( q ) , θ ) | q = s δ ( θ σ ) d θ d σ .
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 D f ( y ( q ) , θ ) with respect to unit vector θ. Therefore Equation (10) can be written as follows:
f ( x ) = 1 8 π 2 S 2 s R ( x , σ ) n x ( s , σ ) y ( s ) σ q { S 2 θ , σ D f ( y ( q ) , θ ) δ ( θ σ ) } | q = s d θ d σ .
Using Grangeat's formula and change of variables ρq defined by ρ = y(q) · σ, we obtain:
1 y ( s ) σ q { S 2 θ , σ D f ( y ( q ) , θ ) d θ } | q = s = 2 ρ 2 f ^ ( ρ , σ ) | ρ = y ( s ) σ = x σ ,
where (ρ, σ) 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:
f ( x ) = 1 8 π 2 I d s | x y ( s ) | 0 2 π φ { sgn ( y ( s ) σ ) n x ( s , σ ) } d φ × 0 2 π q D f ( y ( q ) , cos γ β ( s , x ) + sin γ e x ( s , φ ) ) | q = s d γ sin γ
where β(s, x) = (xy))/‖ xy(s)‖ = (0,0,1),ex(s, φ) = β(s, x) × σx(s, φ), denote
κ x ( s , φ ) = sgn ( y ( s ) σ ) n x ( s , φ ) w m ( s , x ) = lim ɛ 0 ( k x ( s , φ m + ɛ ) k x ( s , φ m ɛ ) ) .
The values of wm(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 ex(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:
f ( x ) : = 1 8 π 2 I m w m ( s , x ) | x y ( s ) | × 0 2 π q D f ( y ( q ) , cos r β ( s , x ) + sin r e x ( s , φ m ) ) | q = s d r sin r d s .
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 cone-beam vertex trajectory consisting of three-orthogonal to each other circular. Furthermore, the filtering lines and the weighting functions will be given in detail.

3. 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 = C1C2C3, as shown in Figure 1.
where:
C 1 = { y ( s ) R 3 : y 1 ( s ) = R cos s , y 2 ( s ) = R sin s , y 3 ( s ) = 0 , s [ 0 , 2 π ] } ;
C 2 = { y ( s ) R 3 : y 1 ( s ) = R cos s , y 2 ( s ) = 0 , y 3 ( s ) = R sin s , s [ 2 π , 4 π ] } ;
C 3 = { y ( s ) R 3 : y 1 ( s ) = 0 , y 2 ( s ) = R cos s , y 3 ( s ) = R sin s , s [ 4 π , 6 π ] } .
Let f(x) is the density function to be constructed. Assume that the function is smooth and vanishes outside the ball:
Ω = { x = ( y 1 , y 2 , y 3 ) | y 1 2 + y 2 2 + y 3 2 r } , 0 < r < R ,
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 (y1, y2, y3) coordinate system, any point on the detector plane can be characterized completely by use of (u, v). It can be shown that:
y 1 = u sin s R cos s , y 2 = u cos s R sin s , y 3 = v , s [ 0 , 2 π ] ;
y 1 = u sin s R cos s , y 2 = v , y 3 = u cos s R sin s , s ( 2 π , 4 π ] ;
y 1 = v , y 2 = u sin s R cos s , y 3 = u cos s R sin s , s ( 4 π , 6 π ] .
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)C1 (the case y(s)C2,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:
σ x = ( sin s cos λ , 1 cos s cos λ , sin s sin λ ) .
The equation of tangent plane, which is a filter plane, is given by:
R ( x , σ ) = { x | ( x y ( s ) ) σ = 0 } .
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:
v ( u : s , λ ) = u ( cos s cos λ ) sin s sin λ + 2 R sin λ .
In the second case, the normal vector to the tangent plane is given by:
e = ( 1 cos s cos λ , cos s cos λ , cos s sin λ ) .
The equation of the filtering line on the planar-detector is given by:
v ( u : s , λ ) = u ( cos λ sin s ) sin λ cos s + 2 R sin λ .

4. 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 C1. We project stereographically the trajectories C2 and C3 onto the detector plane as shown in Figure 3. Let PC2 and PC3 denote these projections. The line PC1 is the projection of trajectory C1. 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 PC2 and PC3 split the entire ROI into three sub-ROIs: PRm, m = 1,2,3. The object point x is projected into the area PRm,m = 1,2,3. Let denotes this projection.
If PR1, 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 kx(s, φ) with equation (14) k is discontinuous if L(φ) is parallel to y'(s) or is tangent to C2. In the first case, such a plane has six points of intersection with C = C1C2C3, 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 L1(see the dot-dashed lines in Figure 3). The corresponding weighting function of filtering line L 1 is wx(s, φ1) = 1/3. In the second case, the number of intersection changes form four to six. So the filtering line is tangent to C2, which is denoted as L2 (see the dot-dashed lines in Figure 3). For the signum function in Equation (14) is unchagened, the corresponding weighting function of filtering line L2 is wx(s, φ2) = 1/12. If PR2, k is discontinuous only L(φ) is parallel to y'(s). So the filtering line L1 is parallel to y(s). The corresponding weight function is wx(s, φ1) = 1/3. If PR3, k is discontinuous if L(φ) is parallel to y'(s) or is tangent to C3. 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 L1 (as shown in Figure 3).The corresponding weighting function of filtering line L1 is wx(s, φ1) = 1/3. In the second case, the number of intersection changes form four to six. So the filtering line is tangent to C3, which is denoted as L3(see the dot-dashed lines in Figure 3). The corresponding weighting function of filtering line L3 is wx(s, φ3) = 1/12. Figure 3 summarizes the above information.

5. 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.

References and Notes

  1. Katsevich, A. Theoretically exact filtered backprojection-type inversion algorithm for spiral CT. SIAM J. Appl. Math. 2002, 62, 2012–2026. [Google Scholar]
  2. Zou, Y.; Pan, X.C. Image reconstruction on PI-lines by use of filtered backprojection in helical cone-beam CT. Phys. Med. Biol. 2004, 49, 2717–2731. [Google Scholar]
  3. Tuy, H.K. An inverse formula for cone-beam reconstruction. J. Appl. Math. 1983, 43, 546–552. [Google Scholar]
  4. Zeng, G.L.; Gullberg, G.T. A cone-beam tomography algorithm for orthogonal circle-and-line orbit. Phys. Med. Biol. 1992, 37, 563–577. [Google Scholar]
  5. Zeng, G.L.; Gullberg, G.T. Implementation of Tuy's cone-beam inversion formula. Phys. Med. Biol. 1994, 39, 493–507. [Google Scholar]
  6. Ning, R.; Tang, X.; Conover, D.; Yu, R. Flat panel detector-based cone beam computed tomography with a circle-plus-two-arcs data acquisition orbit: preliminary phantom study. Med. Phys. 2003, 30, 1694–1705. [Google Scholar]
  7. Wang, X.H.; Ning, R. A cone-beam reconstruction algorithm for circle-plus-arc data-acquisition geometry. IEEE Trans. Med. Imag. 1992, 18, 815–824. [Google Scholar]
  8. Finch, D.V. Cone beam reconstruction with sources on a curve. SIAM J. Appl. Math. 1985, 655–673. [Google Scholar]
  9. Chen, Z.K.; Ning, R. Volume fusion for two-circular-orbit cone-beam tomography. Appl. Opt. 2006, 45, 5960–5966. [Google Scholar]
  10. Noo, F.; Clack, R.; White, T.A.; Roney, T.J. The dual-ellipse cross vertex path for exact reconstruction of long objects in cone-beam tomography. Phys. Med. Biol. 1998, 43, 797–810. [Google Scholar]
  11. Feldkamp, L.A.; Davis, L.C.; Kress, J.W. Practical cone beam algorithm. J. Opt. Soc. Am. 1984, A1, 612–619. [Google Scholar]
  12. Delia, S.; Ivan, B.; Nicolas, P. Studies on circular isocentric cone-beam trajectories for 3D image reconstructions using FDK algorithm. Comput. Med. Imag. Grap. 2008, 32, 210–220. [Google Scholar]
  13. Katsevich, A. A general scheme for constructing inversion algorithms for cone beam CT Internat. J. Math. Math. Sci. 2003, 21, 1305–1321. [Google Scholar]
  14. Katsevich, A. Image reconstruction for the circle and line trajectory. Phys. Med. Biol. 2004, 49, 5059–5072. [Google Scholar]
  15. Zamyatin, A.A.; Katsevich, A.; Chiang, B.S. Exact image reconstruction for a circle and line trajectory with a gantry tilt. Phys. Med. Biol. 2008, 53, N423C–N435. [Google Scholar]
  16. Katsevich, A. Image reconstruction for the circle-and-arc trajectory. Phys. Med. Biol. 2005, 50, 2249–2265. [Google Scholar]
  17. Katsevich, A. Image reconstruction for a general circle-plus trajectory Inverse Problems. Inverse Probl. 2007, 23, 2223–2230. [Google Scholar]
  18. Katsevich, A. On two version of a 3PI algorithms for helical CT. Med. Phys. 2004, 49, 2129–2143. [Google Scholar]
  19. Katsevich, A. 3PI algorithms for helical computer tomography. Adv. Appl. Math. 2006, 36, 213–250. [Google Scholar]
  20. Grangeat, P. Mathematical framework of cone beam 3D reconstruction via the first derivative of the Radon transform. In Mathematical Methods in Tomography. Lect. Notes Math. 1991, 1497, 66–97. [Google Scholar]
Figure 1. Geometries for three-orthogonal circular trajectory and the planar-detector plane.
Figure 1. Geometries for three-orthogonal circular trajectory and the planar-detector plane.
Sensors 09 04606f1
Figure 2. The filtering plane is shown, which is tangent to one of circular trajectories and is perpendicular to σx.
Figure 2. The filtering plane is shown, which is tangent to one of circular trajectories and is perpendicular to σx.
Sensors 09 04606f2
Figure 3. Detector plane with various filtering lines covering projection of ROI shown.
Figure 3. Detector plane with various filtering lines covering projection of ROI shown.
Sensors 09 04606f3

Share and Cite

MDPI and ACS Style

Hu, H.; Zhang, J. Exact Weighted-FBP Algorithm for Three-Orthogonal-Circular Scanning Reconstruction. Sensors 2009, 9, 4606-4614. https://0-doi-org.brum.beds.ac.uk/10.3390/s90604606

AMA Style

Hu H, Zhang J. Exact Weighted-FBP Algorithm for Three-Orthogonal-Circular Scanning Reconstruction. Sensors. 2009; 9(6):4606-4614. https://0-doi-org.brum.beds.ac.uk/10.3390/s90604606

Chicago/Turabian Style

Hu, Hongli, and Jianzhou Zhang. 2009. "Exact Weighted-FBP Algorithm for Three-Orthogonal-Circular Scanning Reconstruction" Sensors 9, no. 6: 4606-4614. https://0-doi-org.brum.beds.ac.uk/10.3390/s90604606

Article Metrics

Back to TopTop