PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
IEEE Trans Med Imaging. Author manuscript; available in PMC 2010 June 15.
Published in final edited form as:
PMCID: PMC2885857
NIHMSID: NIHMS199632

Fast Exact/Quasi-Exact FBP Algorithms for Triple-Source Helical Cone-Beam CT

Yang Lu,* Alexander Katsevich, Jun Zhao, Member IEEE, Hengyong Yu, Senior Member, IEEE, and Ge Wang, Fellow, IEEE

Abstract

Cardiac computed tomography (CT) has been improved over past years, but it still needs improvement for higher temporal resolution in the cases of high or irregular cardiac rates. Given successful applications of dual-source cardiac CT scanners, triple-source cone-beam CT seems a promising mode for cardiac CT. In this paper, we propose two filtered-backprojection algorithms for triple-source helical cone-beam CT. The first algorithm utilizes two families of filtering lines. These lines are parallel to the tangent of the scanning trajectory and the so-called L lines. The second algorithm utilizes two families of filtering lines tangent to the boundaries of the Zhao window and L lines, respectively, but it eliminates the filtering paths along the tangent of the scanning trajectory, thus reducing the required detector size greatly. The first algorithm is theoretically exact for r < 0.265 R and quasi-exact for 0.265 Rr < 0.495 R, and the second algorithm is quasi-exact for r < 0.495 R, where r and R denote the object radius and the trajectory radius, respectively. Both algorithms are computationally efficient. Numerical results are presented to verify and showcase the proposed algorithms.

Index Terms: Cardiac CT, computed tomography (CT), cone-beam CT (CBCT), filtered backprojection (FBP), triple-source

I. Introduction

Thanks to the modern detection technique, projection data can now be collected on a 2-D array of detectors [1], which is much faster than with the conventional single- or multi-row detectors. The more detector rows we use, the larger the longitudinal coverage we have, and the faster the data acquisition is for volumetric scanning. However, the detector array cannot be too large in a cost and dose efficient way, and the mechanical rotation of an X-ray source has already approached its limit. One alternative solution is to use multi-source computed tomography (CT) architecture, which led to the emergence of Siemens' dual-source CT [2]. Given successful applications of dual-source cardiac CT scanners, triple-source cone-beam CT (CBCT) seems a promising mode for cardiac CT [3].

For the dual-source scanner, the minimum rotation interval is 90° plus the fan angle α compared to 180° + α in a third-generation geometry. Since the fan angle α is relatively small, there is > 40% reduction in acquisition time. As a result, the dual-source scanner reduces the data acquisition time to 0.35 s, yielding 83 ms temporal resolution with cardiac gating.

There are major problems with current cardiac CT including dual-source CT [2]. Actually, even in the favorable cases, retrospectively reconstructed cardiac images still suffer from substantial motion blurring for patients who have high and irregular heartbeats because in practice each projection sector covers a projection angular range of a substantial length. Within such an angular range, the heart moves appreciably, especially when it is not in a relative stationary phase. As a benchmark, we routinely achieve ~0.3 mm spatial resolution in spiral CT of the temporal bone where the motion magnitude is much less than that of the heart [4], [5]. On the other hand, the spatial resolution with cardiac CT is at best in the millimeter domain.

As a natural extension to dual-source CT, triple-source CT promises to bring temporal resolution down to 0.2 s and hopefully to be as fast as electron-beam CT. Interestingly, for helical cone-beam reconstruction [6], [7], the trinity (triple-source architecture) is better than the duality (dual-source architecture) because the triple-source helical scan allows a perfect mosaic of longitudinally truncated cone-beam data to satisfy the Orlov condition and yields better noise performance than the dual-source counterpart [8].

Since the 1-D Hilbert transform can be obtained on a generalized PI-line/chord by backprojecting the weighted differential projection data, which was first proposed by Gel'fand and Graev [9] and generalized in the 2D SPECT case [10] and rediscovered in the 2D and 3D CT cases [11], we previously formulated a backprojection-filtration (BPF) approach for triple-source helical CBCT [12] and a corresponding “slow” filtered-backprojection (FBP) algorithm [13]. In this paper, we present two “fast” exact/quasi-exact FBP algorithms in this important case. The algorithms are not specifically designed for cardiac gating which is widely used in cardiac CT. However, the triple-source scanner achieves a data acquisition speed three times faster than the conventional CT, thus the temporal resolution is significantly increased and may be applied without gating in multiple dynamic cases such as perfusion studies in which gating is generally not possible.

This paper is organized as follows. In Section II, we summarize the general reconstruction formulae and the notations of the inter-PI lines and inter-PI arcs. In Section III, we analyze the intersections of the plane through a reconstruction point and the inter-PI arcs associated with the point. In Section IV, we derive the filtering directions and backprojection coefficients for our first algorithm. In Section V, we present the second algorithm after defining a different weight function. In Section VI, we analyze the percentage of the incorrectly weighted Radon planes. In Section VII, we describe the numerical results to validate our algorithms. In Section VIII, we discuss the relevant issues and conclude the paper.

II. Background Materials

A. Geometry of the Triple-Source Helical CBCT

Let f(x) be an object function to be reconstructed. Assume that this function is smooth and vanishes outside the object cylinder

Ω={x=(x1,x2,x3)|x12+x22r2,x3minx3x3max,}0<r<R
(1)

where r is the radius of the object cylinder and R the radius of the scanning cylinder on which a scanning trajectory resides. In the Cartesian coordinate system (x1, x2, x3), the triple-helix trajectories can be expressed as

{y1(s)=(Rcoss,Rsins,h2πs)y2(s)=(Rcos(s+23π),Rsin(s+23π),h2πs)y3(s)=(Rcos(s+43π),Rsin(s+43π),h2πs).
(2)

where h > 0 is the pitch of each helix, and s [set membership] R is the rotation angle. Fig. 1 illustrates the triple-source helical CBCT geometry.

Fig. 1
Geometry of triple-source helical CBCT. Three X-ray sources are rotated around the x3-axis along the helices y1(s), y2(s), and y3(s), respectively. The helices y1(s), y2(s), and y3(s) are on a cylinder of radius R. An object to be reconstructed is confined ...

Previously, we defined inter-helix PI-lines (for simplicity, we call them inter-PI lines thereafter) and extended the traditional Tam–Danielsson window to the Zhao window in the case of triple helices [14]. Specifically, for each source position yj(s), j [set membership] {1, 2, 3}, the corresponding Zhao window is the region on the surface of the scanning cylinder bounded by the nearest helix turn of yj mod 3+1(s) and the nearest helix turn of y(j+i) mod 3+1(s), j [set membership] {1, 2, 3}. In Fig. 2, we use Γ±1 and Γ±2 to denote the boundaries of the Zhao window and the Tam–Danielsson window on the detector plane, respectively. In this paper, our algorithms are implemented on flat-panel detectors. However, they can also be applied to curved detectors. In these cases, the algorithms will use almost the same steps as described, for example, in Noo's “Native geometry” paper [15]. In particular, they can be applied for both flat and curved detectors without any rebinning.

Fig. 2
Illustration of the Zhao window bounded by solid lines Γ±1 and the Tam–Danielsson window bounded by dashed lines Γ±2. The detector plane is represented by the Cartesian coordinate system (u, v).

B. Properties of the Inter-PI Lines and Inter-PI Arcs

Recall that an inter-PI line for yj(s) and yj mod 3+1(s), j [set membership] {1, 2, 3}, is the line that (1) intersects yj(s) at one point and yj mod 3+1 (s) at another point; and (2) the absolute difference between the angular parameter values at the two intersection points is less than 2π [14]. We already proved the existence and uniqueness of the inter-PI line in the following theorem [14].

Theorem 1: Through any fixed x [set membership] Ω, there exists one and only one inter-PI line for any pair of the three helices defined by (2).

In the triple-helix case, there are three inter-PI lines for a fixed x [set membership] Ω and corresponding inter-helix PI-arcs whose end points are along the corresponding helices and share the intersection points of the inter-PI lines. The three inter-PI arcs represent the source trajectory arcs along which the sources illuminate the point x (Fig. 3).

Fig. 3
Inter-PI arcs (thick solid curve-arcs) s1ss1e, s2ss2e, and s3ss3e for x.

III. General Approach

A. Katsevich Scheme

In 2003, Katsevich proposed a general scheme for constructing inversion algorithms for CBCT [16]. It can be stated as follows:

f(x)=14π2Imcm(s,x)|xy(s)|×02πqDf(y(q),cosγβ(s,x)+sinγα(s,x,θm))|q=sdγsinγds
(3)
β(s,x):=xy(s)|xy(s)|
(4)
α(s,x,θ):=β(s,x)×α(s,x,θ)
(5)
cm(s,x):=limε0+(ϕ(s,x,θm+ε)ϕ(s,x,θmε))
(6)
ϕ(s,x,θ):=sgn(αy˙(s))n(s,x,α)
(7)
I:=l=1Lc[al,bl]3,Isy(s)3
|y˙(s)|0onI

where Df(y, β) is the cone-beam transform of f, θ the polar angle in the plane perpendicular to β(s, x), α(s, x, θ) a unit vector perpendicular to β(s, x), θm a point where [var phi](s, x, θ) is discontinuous, n(s, x, α) a weight function, C a finite union of C curves in R3, −∞ < al< bl< ∞, and y(s) := dy/ds.

The aforementioned general inversion formula can be applied to any trajectory that satisfies Tuy's condition, but only when the weight function n(s, x, α) is well designed can the inversion formula have a shift-invariant filtering structure. To derive fast exact FBP algorithms for triple-source helical CBCT, our general approach involves the following key concepts of and analyses on the inflection line, A−, T−, L-, and Bs-curves.

B. Inflection Line

On the detector plane, the boundaries of the Zhao window are expressed as

u(s)=Dsins1coss,υ(s)=Dh2πR(s+Δs)(1coss).
(8)

where D is the distance between the detector and the source, s is the angular parameter relative to the corresponding source position, Δs = −2/3π and Δs = −4/ are for the top and bottom boundaries respectively. Then

u˙(s)=Dcoss1
(9)
u¨(s)=Dsins(coss1)2
(10)
v˙(s)=Dh2πR[1coss(s+Δs)sins](1coss)2
(11)
v¨(s)=Dh2πR(1coss)3[(coss1)(s+Δs)coss+2sins(coss+ssins+Δssins1)]
(12)
d2vdu2=v¨(s)u˙(s)u¨(s)v˙(s)[u˙(s)]3=h2πD(s+Δssins).
(13)

The inflection point exists when d2v/du2 = 0. Thus, we obtain su = 2.6053 and sd = 3.6779 when Δs = −2/3π and −4/3π. The slope of the tangent line at s is computed as

dvdu=v˙(s)u˙(s)=h[1coss(s+Δs)sins]2πR(coss1)=h2πRcoss
(14)

Because cos su = cos sd = −0.8596, the slope is the same (−0.1368h/R) at both inflection points. For practical medical applications, it is common that rFOV≤ 0.5R, and we can include a boundary limitation x12+x22r2,r=0.495R, which is shown as the vertical lines Γl and Γr in Fig. 4. Now, the inflection lines (the tangent lines at s^d and s^u, where s^d and s^u are the projection of y(j+1) mod 3+1(sd) and yj mod 3+1(su), j [set membership] {1, 2, 3} on the detector plane) and the boundary lines split the Zhao window into the following three regions: G1, G2 and G3. Only the points in G1, and G3 can have tangent lines with Γ±1.

Fig. 4
Decomposition of the Zhao window into the regions G1,G2, and G3. Lit, and Lib are the inflection lines at s^u and s^d, respectively.

C. A-Curve and T-Curve

To construct an appropriate weight function, at first we must know how Radon planes intersect with the trajectories. Recall that the number of intersection points only changes when a Radon plane is tangent to the trajectory or contains one PI line/inter-PI line. Hence, if we find all such Radon planes, we can determine the distribution of the intersection points. Since each plane is uniquely determined by its normal vector, in the following sections we use unit vectors instead of the Radon planes. Let us define several key notions introduced in [17]. An A-curve consists of all unit vectors orthogonal to an inter-PI line. A T-curve consists of all unit vectors

α(s)=±(xy)(s))×y˙(s)|(xy(s))×y˙(s)|
(15)

where s belongs to an inter-PI arc. Actually, the A-curve represents all Radon planes containing one inter-PI line, and the T-curve represents all Radon planes tangent to the trajectory. Since there are three inter-PI lines and three inter-PI arcs for a fixed x, there are accordingly three A-curves and three T-curves. Similar to [17], let us use spherical coordinates (θ1, θ2) to describe these curves on the unit sphere

α=(cosθ1sinθ2,sinθ1sinθ2,cosθ2),πθ1π,0θ2π.
(16)

With the identification (θ1, θ2) [congruent with] ((θ1+π) mod 2π, πθ2), each α corresponds to a unique plane through x with the normal vector α.

As an example, the A-curves and T-curves of point x = (0.1, 0, 0) are illustrated in Fig. 5, where R = 1, h = 2π. T1, T2, and T3 stand for the T-curves corresponding to the inter-PI arcs s1ss1e, s2ss2e, and s3ss3e, respectively. Similarly, A1, A2, and A3 are for the A-curves corresponding to the inter-PI lines s1ss2e¯, s2ss3e¯, and s3ss1e¯, respectively.

Fig. 5
Visualization of the domains delimited by the A-curves and T-curves on the surface of the unit sphere. (a) A-curves and T-curves in spherical coordinates for x = (0.1, 0, 0); (b), (c), and (d) the zoom-in versions of the areas bounded by the bottom left, ...

The A-curves and T-curves divide the surface of the unit sphere into several connected domains, in each of which all the planes through x have the same number of intersection points (IPs) with the inter-PI arcs of x. Given an object point and one source trajectory, the number of IPs changes only when a Radon plane is tangent to the trajectory or contains an endpoint of the trajectory. The A-curve represents all planes containing the endpoints of the trajectory, and the T-curve represents all planes tangent to the trajectory. If we pick any Radon plane and rotate it around one direction, the normal vector of this plane forms a curve on the unit sphere. Clearly, only when this curve intersects with the A-curve or T-curve does the number of IPs change. Thus, the A -curve and T-curve define the boundaries of different domains in which the number of IPs is constant. The distribution of IPs over the inter-PI arcs is listed in Table I.

TABLE I
Distribution of IPs on Each of the Domains Delimited by the A-Curves and T-Curves on the Surface of the Unit Sphere. Symbol “-” Indicates no IP on the Corresponding Inter-PI Arc

To determine the distribution of IPs, we first pick a vector α(θ1, θ2) in each domain, and then generate the plane through x and perpendicular to α(θ1, θ2), and compute numerically the number of IPs.

By construction, a T-curve always starts from an A-curve and ends on another A-curve. It can be seen from Fig. 5 that a T-curve is possibly not smooth at some point ac, but the limits of unit tangent vectors at ac and ac+ are equal. Such a point ac is called a “cusp.” The cusp indicates that the two vectors determine the same plane, and ac is the normal vector to that plane. It has been proved in [17] that the cusp is equivalent to the osculating plane Πc(x) which goes through yi(sic(x)), i [set membership] {1, 2, 3}, is parallel to y˙i(sic(x)), y¨i(sic(x)), and contains x (see Fig. 6). On the detector plane, this corresponds to a point where the projected boundary has a point of zero curvature, i.e., the point of inflection [18], [19].

Fig. 6
Illustration of the osculating plane Πc.

The diagram we plot in spherical coordinates deforms smoothly as a function of x. The new diagram is equivalent to the old one if the distribution of IPs remains the same. An essential change could happen when three boundaries intersect each other at one point, which is called “critical event.” A “critical event” may happen in the following seven cases.

  1. Three A-curves intersect at one point.
  2. Three T-curves intersect at one point.
  3. Two A-curves and one T-curve intersect at one point.
  4. One A-curve and two T-curves intersect at one point.
  5. A T-curve becomes tangent to an A-curve at a point of nonsmoothness (i.e., cusp).
  6. When the order of tangency (i.e., the zero derivatives of this order) at the beginning of the T-curve is increased, the T-curve reemerges on the other side of the A-curve.
  7. T-curve develops a smooth dent and becomes tangent to an A-curve.

From Lemmas 2–3, we know Case 1, 2, 4, 6, 7 do not occur for r < 0.495R and Case 5 does not occur for r < 0.265R. Case 3 is possible. If Case 3 takes place, the tangency of T-curve and A-curve will move across another A-curve, then one domain disappears. For example, when x = (0.1, 0, 0) gradually changes to x = (0, −0.15, 0), in Fig. 5(d) the tangency of T3 and A2 will move across A3, and domain D10 disappears [see Fig. 7(a)]. Case 5 is possible for r ≥ 0.265R. If it takes place, a T-curve will only intersect A-curves at the endpoints. That is, the cusp of that T-curve and one domain disappear [see Fig. 7(b)].

Fig. 7
(a) Close-up view of the diagram for x = (0, −0.15, 0). (b) Close-up view of the diagram for x = (0, −0.3, 0).

D. L-Curve

The L-curve is used to split the domain D4 into subdomains, making the weight function n continuous across all the A-curves. This is the key requirement, which will allow us to develop efficient reconstruction algorithms (see Section IV below). Thus, the L-curve should not go across an A-curve. Let us fix x and run s over the three inter-PI arcs, x^ forms a trajectory on the detector plane. Because at the endpoint of the inter-PI arc the line connecting yi(s) and x happens to be an inter-PI line, x^ always starts from one endpoint of the inter-PI arc on the boundary of the Zhao window, and ends at the other one. Hence, whatever the trajectory of x^ is, part of the trajectory is in G2. In other words, x^ will run across one inflection line, then move in G2, and finally cross the other inflection line. Note that x^ on the inflection line indicates a plane containing the inflection line, i.e., a cusp in one T-curve. From Lemma 3, the cusps always belong to the boundary of domain D4. Thus, they can be used as the endpoints of the L-curve. A family of L-curves is formed as follows. Run s over the three inter-PI arcs of x. If x^ is in G2 and above s^u where Γ+1 intersects Lit (Fig. 4), find the plane through x^ and s^u. If x^ is in G2 and below s^d where Γ−1 intersects Lib, find the plane through x^ and s^d. If x^ is in G2 and between s^d and s^u, find the plane through x^ and parallel to the u-axis. Then, we can plot all the normal unit vectors of these planes in the spherical coordinates (θ1, θ2). This gives us three L-curves. The corresponding lines on the detector plane are called L -lines. Fig. 8 shows the L-curves on the diagram in spherical coordinates (θ1, θ2), where L1, L2, and L3 denote the L-curves corresponding to the inter-PI arcs s1ss1e, s2ss2e, and s3ss3e, respectively. As is seen from the above construction, the L-curve always starts and ends on the cusps, and not defined for those parameter values when x^ is not in G2.

Fig. 8
Illustration of L-curves in the spherical coordinates (θ1, θ2).

For r ≥ 0.265R, one or more cusps will disappear if “critical event Case 5” occurs, then the L-curve may start from the intersection of T-curve and A-curve, and end at one A-curve. Also, the L-curve may start from one cusp and end at one A-curve or start from the intersection of T-curve and A-curve, and end at one A-curve. For example, see Fig. 9, L1 starts from the intersection of T3 and A2, and ends at the intersection of T2 and A2; L2 starts from one point on A2, and ends at the cusp of T1; L3 starts from the cusp of T1, and ends at one point on A2.

Fig. 9
L-curves for x = (0.2, −0.3, 0). (a) is the full diagram shows different regions spit by A–, T- and L-curves; (b), (c), (d) are the zoom-in versions of the regions bounded by the up, bottom right, and bottom left circles in (a).

Whatever the endpoints of L-curves are, the L-curves intersects at one point, i.e., θ2 = π or 0 in the spherical coordinates (which corresponds to the plane containing x and parallel to the x1x2 plane). Then D4 is split into several sub-domains. If the endpoints of L-curves are cusps, by Lemma 5, each sub-domain contains only one A-curve. If not, small “line segments” on A-curves may appear and the sub-domains may contain more than one A-curve [see Fig. 9(b)].

E. Bs-Curve

A Bs-curve consists of all unit vectors perpendicular to xyi(s), i [set membership] {1, 2, 3}. Each intersection of Bs- and A-curves corresponds to a plane containing an inter-PI line and yi(s). Each intersection of Bs and T-curves corresponds to a plane tangent to an inter-PI arc and containing yi(s). For example, pick x [set membership] Ω with x^ [set membership] G2, x^ is above L0, where L0 is the projection of the helical tangent at the current source position. Denote L(θ) := DP (s) ∩ Π(x, α(s, x, θ)), where DP(s) is the detector plane corresponding to the source position s, and L(θ) is the projection of the plane through x with the normal vector α(s, x, θ) (Fig. 10). As θ increases, α(s, x, θ) [set membership] β[perpendicular](s, x) rotates clockwise on DP(s), and the following sequence of events takes places. First, Π(x, α(s, x, θ)) intersects s^1ss^2e¯, and a pair of IPs is born. On the unit sphere, this is seen as an intersection of Bs and A1, after which Bs enters D5 [Fig. 11(a)]. Second, Π(x, α(s, x, θ)) intersects s^3ss^1e¯ and another pair of IPs is born. On the unit sphere, this is seen as an intersection of Bs and A3, after which Bs enters D11. Third, a swap of two IPs takes place. On DP(s) this happens when θ = θ0, L(θ0) is parallel to the helical tangent. On the unit sphere, this means that Bs is tangent to T1 at α0 = α(s, x, θ0). Fourth, Bs exits D11 by intersecting T1. On DP(s) this takes place when L(θ0) is tangent to Γ+2. Finally, Π(x, α(s, x, θ)) intersects the L-line. This will not change the number of IPs but it will be useful for construction of the weight function. On the unit sphere, this is seen as an intersection of Bs and L1.

Fig. 10
Domains on the detector plane.
Fig. 11
(a) Bs-curve is tangent to a T-curve in D11 for the source on y1(s). (b) Bs-curve is tangent to a T -curve in D5 for the source on y1(s). (c) Bs-curve passes across the second T-curve when the source on y1(s).

The jumps across an A-curve can only be of two types: from a 1-IP domain to a 3-IP domain and from a 3-IP domain to a 5-IP domain. Note that the Bs-curve is tangent from the inside to T1, which means a swap of two intersection points at α = α0 where sgn(α0 · y(s)) = 0 (see [17]). For a fixed s, if x is allowed to change slightly inside the Zhao window, the tangency point will move from D11 to D5 across A3 (or from D11 to D6 across A1) [Fig. 11(a) and (b)]. If x projects into G1 or G3, the Bs-curve will pass not only through D5 (or D6) and D11 but also through D7 and D12 [Fig. 11(c)]. The similar results can be obtained if the source is on y2(s) or y3(s).

IV. First Fast FBP Algorithm

Now, we are ready to design our first fast FBP algorithm for triple-source helical CBCT. We will start with specifying the weight function n(s, x, α). Then, we will determine the filtering directions by finding the discontinuities of [var phi](s, x, θ) := sgn(α · y(s))n(s, x, α). Then, the backprojection coefficients can be calculated according to (6). Once the filtering lines and the backprojection coefficients are determined, we can use (3) to reconstruct the object.

A. Construction of the Weight Function n(s, x, α)

To have an efficient FBP structure, the weight function n(s, x, α) should be continuous across all A-curves. Thus, the weight function can be defined as shown in Table II. The values in Table II are the weights assigned to IPs. For example, in the D1 domain the Radon plane has only one IP on the inter-PI segment s1ss1e. Accordingly, we assign weight 1 to this IP and use “-” to indicate that there is no IP on the inter-PI segments s1ss1e, s2ss2e, and s3ss3e. In the D11 domain the Radon plane has three IPs on s1ss1e, one IP on s2ss2e, and one IP on s3ss3e. Thus, we assign weight 1 to one IP on s1ss1e and weight 0 to all other IPs.

TABLE II
Weight Function for the First Fast FBP Algorithm

B. Discontinuity of sgn(α · y(s))

To find the backprojection coefficients, let us pick a representative point in each area, determine the discontinuities of [var phi](s, x, θ) := sgn(α · y(s))n(s, x, α), and extend the results by continuity to the entire area. A discontinuity of sgn occurs only when a Bs-curve is tangent to a T-curve from the inside. In other words, the Bs-curve does not go across the T-curve, but stays on one side in a neighborhood of the point of tangency. On the detector plane, L(θ) is parallel to the helical tangent for all x^ [set membership] G1 [union or logical sum] G2 [union or logical sum] G3. Hence, this gives a family of filtering lines parallel to L0, where L0 is the projection of the helical tangent at y(s). The swap of two IPs changes the weight from [var phi]+ = 0 to [var phi] = 1. The backprojection coefficient is computed as c0 = [var phi]+[var phi] = (+1)(0) − (−1)(1) = 1 [Fig. 12(a)].

Fig. 12
(a) Determination of c0. (b) Determination of cl.

C. Discontinuity of n(s, x, α)

By construction, the weight function n is continuous across all inter-PI lines. A discontinuity of n occurs only when a Bs-curve intersects a T-curve or an L-curve. Without loss of generality, pick y1(s0) on s1ss1e. For x^ [set membership] G2, after the swap mentioned in the above paragraph the weight at the current position is zero. Hence, when the Bs-curve passes through a T-curve, i.e., from D11 to D4, n is continuous. Possible jumps of n may only occur when a Bs-curve passes through an L-curve, i.e., from D43 to D41 or from D42 to D41 in Fig. 8. On the detector plane, this occurs when L(θ) overlaps the L-line of x^. Then, the backprojection coefficients are computed as cl = [var phi]+[var phi] = (+1)(1) − (+1)(0) = 1 [Fig. 12(b)]. For x^ [set membership] G1 [union or logical sum] G3, the Bs-curve will not enter D41. Instead, it passes through a second T-curve twice, i.e., from D43 to D12 and from D7 to D1. From Table II, the jumps of n may only occur in the latter intersection. On the detector plane, this happens when L(θ) overlaps the line tangent to Γ±1. Then, the backprojection coefficients are computed as ct = [var phi]+[var phi] = (+1)(1) − (+1)(0) = 1.

D. Required Detector Area

Fig. 13(a) and (b) summarizes the filtering lines and the backprojection coefficients discussed above. In these figures, L0' is the line parallel to L0 and Ll denotes the L-line. To implement the proposed algorithm, the filtering lines cannot be truncated. Thus, the detector size should be large enough to cover the area bounded by Γr, Γl, Lmax, and Lmin, where Lmax and Lmin are the lines across the intersections of (1) Γl and Γ+1 and (2) Γr and Γ−1, respectively, and parallel to L0 (Fig. 14). The required detector area is determined by two factors: 1) the ratio of the pitch h and the scanning radius R and 2) the ratio of the object support radius r and the scanning radius R. If R is fixed, the required detector area grows as h or r increases.

Fig. 13
(a) Filtering lines in the case of x^ [set membership] G1 [union or logical sum] G3 for the first fast FBP algorithm. (b) Filtering lines in the case of x^ [set membership] G2 for the first fast FBP algorithm.
Fig. 14
Required detector area is bounded by Γr, Γl, Lmax, and Lmin for the first algorithm, and by Γr, Γl, Lmax', and Lmin' for the second algorithm.

V. Second Fast FBP Algorithm

Again, the design of our second fast FBP algorithm starts with specifying new weights (Table III). By construction, n(s, x, α) is continuous across all inter-PI lines. More importantly, a swap of two IPs takes place when a Bs-curve becomes tangent to a T-curve [17], and n(s, x, α) changes from +1 to −1. Recall that the discontinuity of sgn(α · y(s)) appears only when a Bs-curve is tangent to a T-curve from inside. Since both n(s, x, α) and sgn(α · y(s)) are discontinuous at that point, the function [var phi](s, x, θ) := sgn(α · y(s))n(s, x, α) is continuous. Thus, the filtering operation along the tangent of the scanning trajectory is eliminated.

TABLE III
Weight Function for the Second Fast FBP Algorithm

A discontinuity of n can occur only when a Bs-curve intersects a T-curve or an L-curve. Following the discussion in Section IV, jumps of n may occur when 1) a Bs-curve passes through a T-curve, i.e., from D11 to D4 or from D4 to D12 in Fig. 11(a) and (c), and 2) a Bs-curve passes through an L-curve, i.e., from D42 to D41 or from D43 to D41 in Fig. 8. On the detector plane, this gives two families of filtering lines: the lines tangent to Γ±2 or Γ±1 and the L-lines. Note that the filtering lines tangent to Γ±1 are different from those for our first fast FBP algorithm (Fig. 15), because the discontinuity of n(s, x, α) occurs on the different side of the cusp. Then, the backprojection coefficients can be calculated as ct = [var phi]+[var phi] = (+1)(1) − (+1)(0) = 1 and cl = [var phi]+[var phi] = (+1)(1) − (+1)(0) = 1.

Fig. 15
Filtering lines for our two fast FBP algorithms when x^ is above Lit · Lt2 and Lt1 are for the first and second algorithms, respectively.

The reconstruction formula for the second algorithm is very similar to that for the first algorithm. The only difference lies in the selection of the filtering lines. For clarity, our second fast FBP algorithm is illustrated in Fig. 16. Because the filtering paths along the tangent of the scanning trajectory are eliminated, the required detector area is reduced by at least 30% (Fig. 14).

Fig. 16
Illustration of the second fast FBP algorithm.

VI. Analysis on the Quasi-Exactness

By Lemma 3, there are two types of “line segments” according to different critical events. First, let us consider the “line segment” related to a critical event in Case 3. Recall that before entering D4 the Bs-curve will be tangent to a T-curve. For the first algorithm, at the tangency the weight n changes from 1 to 0, then it does not change whether the Bs-curve enters D4 across an A-curve or a T-curve. For the second algorithm, if the weight n changes from 1 to −1 at the tangency, then it will jump from −1 to 0 when the Bs-curve enters D4. If the Bs-curve enters D4 across an A-curve [i.e., the “line segment”; see Fig. 7(a)], the FBP structure is ruined. Thus, the critical event in Case 3 will only affect the second algorithm, without damaging the FBP structure of the first algorithm. Then, let us consider the “line segment” related to a critical event in Case 5, which only occurs when r ≥ 0.265R. From the discussion in Section III.D such “a line segment” is the boundary of the one-IP and three-IP domains. Hence, for both algorithms the weight n will jump from 1 to 0 if the Bs-curve enters D4 across the “line segment,” and the FBP structure is ruined (if the exactness is to be maintained. By sacrificing the exactness we retain the FBP structure). Consequently, the first algorithm is theoretically exact for r < 0.265R and not exact for 0.495Rr ≥ 0.265R, and the second algorithm is not exact for r < 0.495R.

Since the algorithms are not always exact, we are interested in estimating what percentage of the Radon planes is incorrectly calculated. If the Radon planes with approximate weighting only have a small percentage, the algorithms can be considered quasi-exact, and we can still reconstruct with high image quality.

First, let us consider the incorrectly weighted planes caused by the critical event in Case 3. It appears in the area r < 0.495R. Let us fix x for r < 0.495R, denote the intersection of the Radon plane and the detector plane as L(θ), θ [set membership] [0, 2π], run s over the three inter-PI arcs, and see what happens with x^ and L(θ). Based on the discussion in Section III-E, for x^ in G2, if the critical event occurs, the Bs-curve will first intersect a T-curve, and then go cross an A-curve. For example, in Fig. 11(a) the Bs-curve will first intersect T1, then enter D4 across A3. On the detector plane, this corresponds to that L(θ) intersects the tangent of Γ+2 before the inter-PI line s^1es^3s¯ while L(θ) is rotated clockwise. Therefore, the Radon planes between the tangent of Γ+2 and s^1es^3s¯ are not exactly weighted. Because the slope of s^1es^3s¯ is positive and the slope of the tangent of Γ+2 is less than h/2πR, the percentage of the incorrectly weighted Radon planes is less than ρ = β/π, where α<β<2arctan((2)/(3)tan(α/(2)), α = arctan(h)/(2πR) (see Appendix). It is common that h/R < 0.2 in practical applications, hence ρ < 1.17%. For x^ in G1 (G3), if the critical event occurs, the Bs-curve will first go across an A-curve, and then over a T-curve. On the detector plane, this corresponds to the case when L(θ) intersects the tangent to Γ+1−1) before the inter-PI line s^3es^2s¯ while L(θ) is rotated clockwise. Hence, the Radon planes between the tangent of Γ+1−1) and s^3es^2s¯ are not exactly weighted. Because the slope of s^3es^2s¯ is negative and the slope of the tangent of Γ+1−1) is more than −0.35h/R, the percentage of the incorrectly weighted Radon planes is less than ρ = 2.57% for h/R = 0.2.

Then, let us consider the incorrectly weighted planes caused by the critical event in Case 5. Recall that the L-curves are used to split the domain D4 into sub-domains, making the weight function n continuous across all the A-curves, and the cusps are the starting and ending points of the L-curves. If the endpoints of the L-curves are not the cusps, there will be small fractions (or “line segments”) on the A-curves, making the weight function n discontinuous across them and ruining the FBP structure of our algorithms. It possibly occurs for r ≥ 0.265R. As discussed above, the Bs -curve will first enter a 1-IP domain from a 3-IP domain across the “line segment,” and then pass through an L-curve. On the detector plane, this corresponds to that L(θ) intersects the inter-PI line s^3es^2s¯ before the L-line while L(θ) is rotated clockwise. Recall that if the cusp is not in D4, s^2s is possibly to the left of s^ u or s^3e is to the right of s^ d. Thus, the slope of s^3es^2s¯ is always more than −0.35h/R. Because the slope of the L-curve is never positive, the percentage of the incorrectly weighted Radon planes is less than ρ = 2.57% for h/R = 0.2. On the other hand, based on the discussion on Lemma 3, one or more cusps may possibly remain even when r ≥ 0.265R, which means that less “line segments” related to critical events will appear in Case 5, and in fact more Radon planes may be correctly weighted.

VII. Numerical Results

To verify and showcase our proposed fast FBP algorithms, numerical tests were performed using the Clock phantom. This phantom consists of ellipses, as parameterized in Tables IV. In the simulations, the origin of the reconstruction coordinate system was set to the center of each phantom. The spherical phantom support was of 375 mm for the experiment. Three sources were arranged uniformly along a circle with their corresponding detectors on the opposite side. The source-detector distance was 1000 mm. Projections were generated from 1000 view angles while the sources and the detectors were constantly moved along three helixes in one turn. The helix was of 750 mm in radius and 100 mm in pitch. The detector plane consisted of 1300 × 200 detection elements of 1.0 × 1.0 mm2.

TABLE IV
Parameters of the clock phantom

The implementation of our algorithms consists of the following steps.

  • Step 1) Differentiate each projection with respect to variable s.
  • Step 2) For each yi(s),i [set membership] {1,2,3} perform the Hilbert transform of derivative data along the given filtering directions on the corresponding detector plane.
  • Step 3) Backproject the filtered data on the inter-PI segments to reconstruct the object point.

A more detailed description of the proposed algorithm is similar to [20] except for the following differences: 1) in triple-helices geometry the filtered data are backprojected on inter-PI segments; 2) there are two families of filtering lines for each algorithm so that each point on the detector plane will be filtered twice. Since our algorithms allow shift-invariant filtration, all results are in Cartesian coordinates directly, and there is no coordinate transform like what we have used in the slow-FBP algorithm [13] or BPF algorithm [12].

The algorithms were coded in MATLAB and executed on a regular PC (Intel Core2 Duo CPU 3.06 GHz, 4 GB RAM). Reconstructed images are shown in Fig. 17. Our numerical results show that in the case of r = 0.495R both two algorithms produced high quality images.

Fig. 17
Reconstructed images of the Clock phantom with r = 375 mm using (a) the first fast FBP algorithm and (b) the second fast FBP algorithm. Images were reconstructed at x3 = 0 mm and displayed in [0.95, 1.05]. Additionally, (c) and (d) show the differences ...

VIII. Discussions and Conclusion

Although our previously published BPF algorithm for triple-source helical CBCT can indeed produce excellent image quality, FBP algorithms (either “slow” or “fast”) are computationally desirable for several reasons, such as being amendable for parallel processing as discussed in [13]. In particular, while the computational structures of our BPF algorithm and FBP algorithms are quite similar, the FBP algorithms avoid densely sampled intermediate reconstruction in the PI-line-based coordinate system, and more importantly they can reconstruct a region of interest (ROI) or volume of interest (VOI) much more efficiently than the BPF counterpart. Note that ROI/VOI reconstruction is very common in medical imaging. A related technology called “interior tomography” is being actively developed to target this type of problems [21], [22]. Then, an interesting possibility would be to develop tripe-source interior CBCT.

The proposed two fast exact/quasi-exact FBP algorithms for triple-source helical CBCT have their advantages and disadvantages. From the perspective of exact reconstruction, the first algorithm is more desirable than the second algorithm because it is not affected by critical events in Case 3. However, in terms of efficient data acquisition, it requires a larger detector area than the second algorithm. In the medical CT field, the rectangular detector shape is most popular, and the helical pitch may be varied case by case. Therefore, it is practically possible to have projection data for reconstruction using either or both of our two fast FBP algorithms. A detailed comparison of these two algorithms should be done in terms of image quality and radiation dosage. Other designs for fast FBP in the triple-source helical CBCT are also possible.

One concern with triple-source CT is the increment in scattering artifacts, which was addressed in our earlier paper for the BPF formulation [13]. Although the height of the Zhao window for triple-source helical CT is one third of the Tam-Danielsson window, our proposed fast FBP formulation requires additional data beyond the Zhao window. In this case, the total dose delivered to the patient can be controlled to be comparable to that in the single-source case, and the total amount of scattering would be similar as well. Another way is to utilize source multiplexing so that not all the sources are on simultaneously. Other scattering correction means, such as interior tomography [21], [22], may be applied as well.

In conclusion, based on our previous triple-source helical CBCT work and guided by the Katsevich general reconstruction scheme, we have proposed two fast FBP algorithms for triple-source helical CBCT. The proposed algorithms utilize the inter-PI line and inter-PI arcs, and have a shift-invariant filtering structure. Unlike our slow-FBP algorithm performing filtration spatial-variantly line by line, the proposed fast-FBP algorithms filter projection data spatial-invariantly view by view, representing a significant computational benefit. Since triple-source helical CBCT may triple temporal resolution, it seems a promising mode for cardiac CT and other applications, and our proposed algorithms may find applications in this context.

Acknowledgments

This work was supported in part by the National Science Foundation of China (30570511 and 30770589), the National High Technology Research and Development Program of China (2007AA02Z452), the National Basic Research Program of China (2010CB834300), and NIH/NIBIB (EB002667, EB004287, and EB007288). The work of A. Katsevich was supported in part by the National Science Foundation under Grant DMS-0806304.

The authors would like to thank the anonymous reviewers for their comments and suggestions.

Appendix

A. Auxiliary Lemmas

Let us fix a point x [set membership] Ω and find its three associated inter-PI lines as shown in Fig. 3. Then, we can select a source position s(sjs,sje), j [set membership] {1, 2, 3} and determine how the inter-PI lines project onto the corresponding detector plane. For simplicity, here and below the projection of yj(s), j [set membership] {1, 2,3} on a detector plane is denoted by S^ j.

Lemma 1: On a detector plane, the slopes of the projected inter-PI lines s^jss^jmod3+1e¯ and s^(j+1)mod3+1ss^je¯ are always positive, and that of the inter-PI line s^jmod3+1ss^(j+1)mod3+1e¯, j [set membership] {1, 2, 3} is always negative.

Proof: Without loss of generality, suppose the source position we select is y1(s0), s0(s1s,s1e). By construction, s2es1s<2π, s1es3s<2π, and s3es2s>0. Hence, the projections of s1s, and s3e and s3e are always to the left of those of s2e, s2s, and s1e respectively (Fig. 18). When s(s1s,s1e) changes, the point x^, i.e., the projection of x onto the detector plane, moves inside the region G := G1 [union or logical sum] G2 [union or logical sum] G3. Clearly, x^ reaches its highest (resp., lowest) position in the vertical direction when x^ is at the intersection of Γl and Γ+1 (resp., of Γr and Γ−1). Also, the vertical coordinates at these points are υmax = 1.3794(Dh)/(2πR) and υmin = −1.3794(Dh)/(2πR), respectively. Moreover, the lowest point on Γ+2 and the highest point on Γ−2 are υ'min=1.3801(Dh)/(2πR) and υ'min=1.3801(Dh)/(2πR), respectively. Evidently, υmaxυ'min and υ'maxυmin. Since s^1s and s^3s are to the left of s^2e and s^1e, the slopes of the inter-PI lines s^2es^1s¯ and s^1es^3s¯ are positive for all x^ in G.

Fig. 18
Projected inter-PI lines s^1ss^2e¯, s^2ss^3e¯ and s^3ss^1e¯ on the detector plane. Thick curve segments denote the inter-PI arcs.

As follows from the definition in [14], the inter-PI line s2ss3e¯ satisfies:

{x1=Rtcoss2s+R(1t)coss3ex2=Rtsins2s+R(1t)sins3ex3=h2πt(s2s2π3)+h2π(1t)(s3e4π3)
(17)

where G2 and s2s, s3e(s0,s0+2π). Let

x1=r0cosμ0,x2=r0sinμ0
(18)

where r0 [set membership] [0, 0.495R] and μ0 [set membership] [0, 2π]. We have

coss3es2s2=r0sin(μ0s2s)R2+r022Rr0cos(μ0s2s).
(19)

Then, (19) can be rewritten as

coss3es2s2=sin(μ0s2s)(Rr0+cos(μ0s2s))2+sin2(μ0s2s).
(20)

When μ0s2s is fixed and r0 is reduced, |cos(1)/(2)(s3es2s)| decreases. Therefore, the right side of (20) reaches its maximum or minimum when r0 is maximized, i.e., at r0 = 0.495R. Those maximum and minimum values can be numerically calculated [Fig. 19(a)], and we have

Fig. 19
(a) Plot of (20) with r0 = 0.4951R. (b) Plot of Φ(s) over the range s [set membership] (0, 2π); (c) Plot of Ψ (s2s) over the range 0s2ss04.1773.
0.4949cos12(s3es2e)0.4949

or

2.1062s3es2s4.1770.
(21)

Next we prove that 0s2ss04.1773 implies

υ1(s3es0)υ2(s2ss0)>0
(22)

where υ1(s) = (Dh)/(2πR)(s – 4π/3)/(1 – cos s) and υ2(s) = (Dh)/(2πR)(s–2π/3)/(l–cos s). Fig. 19(b) shows the function Φ(s) = 1 – cos s – (s – 4/3π)sin s in the range s [set membership] (0,2π), demonstrating that υ̇1(s) is always positive. Hence υ1(s3es0) is monotonically increasing. Let us fix s3e=s2s+2.1062 and plot the function Ψ(s2s)=υ1(s3es0)υ2(s2ss0) [Fig. 19(c)]. Clearly, this function is always positive in the range 0s2ss04.1773. Note that s = 4.1773 + s0 is the intersection of Γl and Γ+1 and s^2s cannot be to the left of this point (otherwise, x^ is outside G).

Equation (22) indicates that s^2s is always lower than s^3e in the vertical direction. Since s^2s is to the right of s^3e, the slope of the inter-PI line s^2ss^3e¯ is always negative. Due to symmetry, the other two cases s(s2s,s2e) and s(s3s,s3e) can be handled similarly. This finishes the proof.

Lemma 2: A T-curve cannot be tangent to an A-curve at an interior point of T.

Proof: The interior point of T-curve is any point of the T-curve except an endpoint. It has been proved in [17] that a T-curve is smooth everywhere, except possibly at a cusp. Let sic, i [set membership] {1, 2,3} be the point where the cusp occurs. Suppose a T-curve is tangent to an A-curve at α(s). If s=sic(sis,sie), where sis, sie are the endpoints of the T-curve, then the osculating plane Πc(x) intersects the helix yi(s) at only one point and it contains one inter-PI line. By construction, ΠC(x) intersects the detector plane at the asymptote of the Tam-Danielson window boundary and x^ belongs to the asymptote. Connecting x^ and s^is, x^ and s^ie we get two inter-PI lines. Clearly, ΠC(x) will not contain any of them. By Lemma 1, the third inter-PI line has a negative slope, thus it will not overlap the asymptote. Hence, ΠC(x) will not contain it. Consequently, ssic and T-curve is smooth in a neighborhood of α(s). Let θ be the polar angle for the great circle (xyi(sa)),sa{sis,sie},i{1,2,3}. Then, the A-curve consists of all the unit vectors α1(θ) [set membership] (xyi(sa))[perpendicular]. Clearly, [alpha with dot above]1(θ) is perpendicular to α1(θ) and (xyi(sa)). By construction, the T-curve is tangent to the A-curve at α(s). Hence, [alpha with dot above](s) is parallel to [alpha with dot above](θ). That is, [alpha with dot above](s) is perpendicular to α1(θ) and (xyi(sa)). Because [alpha with dot above](s) is also perpendicular to (xyi(st)),st(sis,sie),(xyi(st)) is parallel to (xyi(sa)) and st = sa, which contradicts the assumption that st is an inner point of the T-curve. This finishes the proof.

Lemma 3: Case 1, 2, 4, 6, 7 do not occur for r < 0.495R and Case 5 does not occur for r < 0.265R.

Proof: Cases 1 and 2 are impossible because they mean that there will be a plane containing three inter-PI lines or tangent to three inter-PI arcs.

In Case 4, there will be one plane Π containing one inter-PI line and tangent to two inter-PI arcs. Assume this inter-PI line is s1ss2e¯, let us pick a point s=s1s on yi(s) and denote L = Π [union or logical sum] DP (s0). By construction, on the detector plane x^ is on Γ+1 and it overlaps s^2e. Then, L may be tangent to Γ+1 and Γ+2 or Γ+1 and Γ−1, see Fig. 20. Because the points of tangency are on the inter-PI arcs, the endpoint s^1e is to the left of the tangency for case A and s^3sis to the right of the tangency for case B. Connecting s^1e and x^ (or s^3s and x^) we find that the slope of the inter-PI line s^3ss^1e¯ is negative. By Lemma 1, these two cases are impossible.

Fig. 20
Possible locations of the “critical event” for Case 4.

In Case 5, there will be one plane containing one inter-PI line and tangent to one inter-PI arc at the inflection point. Thus, the inter-PI line on the detector overlaps the inflection line, i.e., s^2ss^3e¯ overlaps Lit when, s2s=su and s3e=2su, where su is difference between s2s and s3e. Let us look at (20). When μ0s2s is fixed, the absolute value of cos(s3es2s)/(2) is monotonically decreasing when ro is reduced. Then, the range of s3es2s is narrowed. In other words, the difference between s3e and s2s becomes closer to π. Case 5 occurs when s3es2s=su. If we want to exclude Case 5, the range of s3es2s cannot cover the value su = 2.6053. Hence, the minimum range of s3es2s is 2.6053<s3es2s<2π2.6053. That is, 0.2649<cos(s3es2s)/(2)<0.2649 and cos(s3es2s)/(2) reaches its extreme when r0 = 0.265R. From (20) and (21) we have s3es2s=su only for r ≥ 0.265R, which contradicts our condition. Hence, Case 5 is impossible for r < 0.265R.

In Case 6, a T-curve will intersect one A-curve twice before meeting a cusp. Suppose this takes place at inter-PI arc s1ss1e.^. Let us pick a point s=s1s on y1(s) and see what happens on the detector plane when s moves. By construction, the plane II containing y1(s) and x intersects the detector plane at the line L which is parallel to the helix tangent across x^. At s=s1s, x^ is on Γ+1 and Π contains inter-PI line s^1ss^2c¯. As s moves along y1(s), x^ moves downwards. Notice that L is parallel to the asymptote of the Tam–Dannielson window, so it will not intersect Γ−2 provided that x^ moves across the asymptote, at where the cusp occurs. Hence, n will not contain the inter-PI line s^1ss^2e¯ and Case 6 is impossible.

By Lemma 2, Case 7 is impossible. This finishes the proof. Let G21 be the area bounded by Lit, Γl' and the line υ = Dh(su – (2)/(3)π)/[2π R(l – cos su)] (i.e., the line across s^u and parallel to the u -axis) and G22 the area bounded by Lib, Γr' and the line υ = Dh(sd – (4)/(3)π)/[2πR(l – cos sd)] (i.e., the line across s^ d and parallel to the u-axis, see Fig. 21). Γl' and Γr' are the vertical lines across S^ d and s^ u. They correspond to the boundary condition r0 = 0.265R.

Fig. 21
Illustration of the regions G21 and G22.

Lemma 4: The inflection point s^ u (s^ d) is inside the inter-PI arc when x^ is in G21(G22).

Proof: By Lemma 3, any point in the area r < 0.265R has three cusps in the diagram. Note that there is one IP in each inter-PI arc within D4. Since all three cusps are in D4, an osculating plane of one inter-PI arc intersects two other inter-PI arcs exactly once at one point. Assume that this osculating plane Πc contains x. Consider Πc of the second inter-PI arc (i.e., of y2 (s)). Let s10 be the point where it intersects the first inter-PI arc (i.e., on y1 (s)). Let us move s along the first inter-PI arc and observe what happens with x^ on the detector. When s=s1s, x^ enters the Zhao window through Γ±1. When s=s10, x^ belongs to Lit. As follows from the diagram, the point su must be inside the second inter-PI arc, i.e., between s2s and s2e. As the point s moves further, the difference s2ss becomes smaller, and the point s2s moves to the right of s^ u along Γ+1. The inter-PI line s^1ss^2e¯ has a positive slope. Thus, as long as x^ is inside G21, the point s^2e is always to the left of s^ u. The case where x^ is in G22 can be similarly treated. This proves Lemma 4.

Lemma 5: An L -curve never intersects an A -curve for r < 0.265R.

Proof: If an L- and A -curves intersect, an L -line through x^ overlaps the inter-PI line. Let us pick one point on the inter-PI arc s3ss3e and consider the slope of the inter-PI line on the corresponding detector plane. By construction, an L-curve always starts from a cusp of a T-curve and ends on a cusp of another T-curve. For the osculating plane ΠC, its intersection with the detector plane is the line tangent to Γ+1−1) at s^ u(s^ d). By Lemma 4, the endpoints of the inter-PI arc on Γ+1−1) are on different sides of s^ u(s^ d), and one of them on the right (left) side is also an endpoint of the inter-PI line for helices y2(s) and y3(s), and denoted as s^2s(s^3e) in Fig. 22.

Fig. 22
Relationship among the inter-PI line La, L-line Ll and inflection line Lit(Lib). (a) x^ in G2 and above s^ u, and (b) x^ in G2 and below s^d.

If x^ is in G2 and above s^ u, we can have the L-line by connecting x^,s^ u, and the inter-PI line by connecting x^, s^2s. If x^ is in G2 and below S^ d, we can have the L-line by connecting x^,s^ d, and the inter-PI line by connecting x^, s^3e. Clearly, in any case the slope of L-line is between zero and the slope of the inter-PI line. That is, the L-line will not overlap the inter-PI line. For other x^, the L-line is parallel to the u axis. By Lemma 1, it is always between two inter-PI lines and will not overlap with any of them.

For the point on other inter-PI arcs, the situation is the same. This finishes the proof.

B. Derivation of Formulae in Section VI

The angles in the detector plane are not equivalent to that in the spherical coordinate system because the detector plane is parallel to the x3-axis while the plane in spherical coordinates is perpendicular to yi(s) – x, i [set membership] {1,2,3}. To estimate the percentage of the incorrectly weighting planes, a transformation is necessary.

Let S = (0, 0, R) denote the source position, x one point in the sphere centered at o with radius r, A = (d, 0,0) the projection of x, B and C two points on the detector plane, and θ the angle between lines AB and AC. The plane perpendicular to the line SA intersects the lines SB and SC at P1 and P2, respectively. The angle [var phi] between lines AP1 and AP2 is what we need (see Fig. 23).

Fig. 23
Angular transformation from a spherical coordinate system to a detector.

By construction

{(P1A)(SA)=0(P1S)=c1(P1B)
(23)
{(P2A)(SA)=0(P2S)=c2(P2C)
(24)

where B = (d + d cos α, d sin α,0), C = (d + d cos(α + θ), d sin(α + θ), 0), c1 and c2 are two nonzero constants. Solving (23) and (24), we have

P1=(d(R2+d2)(1+cosα)R2+d2(1+cosα),d(R2+d2)sinαR2+d2(1+cosα),d2RcosαR2+d2(1+cosα))
(25)
P2=(d(R2+d2)(1+cos(α+θ))R2+d2(1+cos(α+θ)),d(R2+d2)sin(α+θ)R2+d2(1+cos(α+θ)),d2Rcos(α+θ)R2+d2(1+cos(α+θ)))
(26)
cosφ=(P1A)(P2A)(P1A)(P2A)=cosαcos(α+θ)+(1+k2)sinαsin(α+θ)1+k2sin2α1+k2sin2(α+θ).
(27)

where k = d/R.

Computing the derivative of (27) we find that [var phi] reaches its maximum at α = −θ/2 or α = −θ/2 + π. In this case, (27) becomes

cosφ=cos2θ2(1+k2)sin2θ2cos2θ2+(1+k2)sin2θ2.
(28)

For a fixed θ, cos [var phi] is monotonically decreasing while k is increased. Since r < 0.495 R, k reaches its maximum 1/3 when the line OX is perpendicular to the line SA and its minimum 0 when the line OX is parallel to the line OS. Hence,

θ<φ<arccos34tan2θ23+4tan2θ2

or

θ<φ<2arctan(23tanθ2).

Footnotes

Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org.

Contributor Information

Yang Lu, Department of Biomedical Engineering, Shanghai Jiao Tong University, Shanghai 200240, China (nc.ude.utjs@gnayvl)

Alexander Katsevich, Department of Mathematics, University of Central Florida, Orlando, FL 32816 USA (ude.fcu.liam@ivestaka)

Jun Zhao, Department of Biomedical Engineering, Shanghai Jiao Tong University, Shanghai 200240, China (nc.ude.utjs@oahznuj).

Hengyong Yu, SBES Division & ICTAS Center for Biomedical Imaging, Virginia Tech, Blacksburg, VA 24061 USA (gro.eeei@uy-gnoygneh)

Ge Wang, SBES Division & ICTAS Center for Biomedical Imaging, Virginia Tech, Blacksburg, VA 24061 USA (gro.eeei@gnaw-eg).

References

1. Wang G, Crawford CR, Kalender WA. Guest Editorial: Multirow detector and cone-beam spiral/helical CT. IEEE Trans Med Imag. 2000 Sep;19(9):817–821. [PubMed]
2. Flohr TG, McCollough CH, Bruder H, Petersilka M, Gruber K, Süß C, Grasruck M, Stierstorfer K, Krauss B, Raupach R. First performance evaluation of a dual-source CT (DSCT) system. Eur Radiol. 2006;16:256–268. [PubMed]
3. Kachelrieß M, Knaup M, Kalender WA. Multithreaded cardiac CT. Med Phys. 2006;33:2435–2447. [PubMed]
4. Vannier M, Wang G. Spiral CT refines imaging of temporal bone disorders. Diagnostic Imag. 1993;15:116–121. [PubMed]
5. Wang G, Zhao S, Yu H, Miller C, Abbas P, Gantz B, Lee S, Rubinstein J. Design, analysis and simulation for development of the first clinical micro-CT scanner1. Academic Radiol. 2005;12:511–525. [PubMed]
6. Wang G, Lin TH, Cheng PC, Shinozaki DM. A general cone-beam reconstruction algorithm. IEEE Trans Med Imag. 1993 Sep;12(3):486–496. [PubMed]
7. Katsevich A. Theoretically exact filtered backprojection-type inversion algorithm for spiral CT. SIAM J Appl Math. 2002;62:2012–2026.
8. Orlov SS. Theory of three-dimensional reconstruction. 1. Conditions of a complete set of projections. Sov Phys Crystallogr. 1975;20:312–314.
9. Gel'fand I, Graev M. Crofton's function and inversion formulas in real integral geometry. Functional Anal Its Appl. 1991;25:1–5.
10. Rullgård H. An explicit inversion formula for the exponential Radon transform using data from 180. Arkiv för Matematik. 2004;42:353–362.
11. Wang G, Ye Y, Yu H. Approximate and exact cone-beam reconstruction with standard and non-standard spiral scanning. Phys Med Biol. 2007;52:1–13. [PubMed]
12. Zhao J, Jiang M, Zhuang T, Wang G. An exact reconstruction algorithm for triple-source helical cone-beam CT. J X-ray Sci Technol. 2006;14:191–206.
13. Zhao J, Jin Y, Lu Y, Wang G. A filtered backprojection algorithm for triple-source helical cone-beam CT. IEEE Trans Med Imag. 2009 Mar;28(3):384–393. [PMC free article] [PubMed]
14. Zhao J, Jiang M, Zhuang T, Wang G. Minimum detection window and inter-helix PI-line with triple-source helical cone-beam scanning. J X-ray Sci Technol. 2006;14:95–107.
15. Noo F, Pack J, Heuscher D. Exact helical reconstruction using native cone-beam geometries. Phys Med Biol. 2003;48:3787–3818. [PubMed]
16. Katsevich A. A general scheme for constructing inversion algorithms for cone-beam CT. Int J Math Math Sci. 2003;21:1305–1321.
17. Katsevich A. 3PI algorithms for helical computer tomography. Adv Appl Math. 2006;36:213–250.
18. Kapralov M, Katsevich A. A one-PI algorithm for helical trajectories that violate the convexity condition. Inverse Problems. 2006;22:2123–2143.
19. Kapralov M, Katsevich A. On the study of 1PI algorithms for a general class of curves. SIAM J Imag Sci. 2008;1:418–459.
20. Yu H, Wang G. Studies on implementation of the Katsevich algorithm for spiral cone-beam CT. J X-ray Sci Technol. 2004;12:97–116.
21. Ye Y, Yu H, Wei Y, Wang G. A general local reconstruction approach based on a truncated Hilbert transform. Int J Biomed Imag. 2007;2007:63634. [PMC free article] [PubMed]
22. Kudo H, Courdurier M, Noo F, Defrise M. Tiny a priori knowledge solves the interior problem in computed tomography. Phys Med Biol. 2008;53:2207–2231. [PubMed]