Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Phys Med Biol. Author manuscript; available in PMC 2010 June 21.
Published in final edited form as:
PMCID: PMC2888295

Exact image reconstruction with triple-source saddle-curve cone-beam scanning

Yang Lu,1,2 Jun Zhao,1,2,3 and Ge Wang2


In this paper, we propose an exact shift-invariant filtered backprojection (FBP) algorithm for triple-source saddle-curve cone-beam CT. In this imaging geometry, the x-ray sources are symmetrically positioned along a circle, and the trajectory of each source is a saddle curve. Then, we extend Yang’s formula from the single-source case to the triple-source case. The saddle curves can be divided into four parts to yield four datasets. Each of them contains three data segments associated with different saddle curves, respectively. Images can be reconstructed on the planes orthogonal to the z-axis. Each plane intersects the trajectories at six points (or three points at the two ends) which can be used to define the filtering directions. Then, we discuss the properties of these curves and study the case of 2N + 1 sources (N ≥ 2). A necessary condition and a sufficient condition are given to find efficient curves. Finally, we perform numerical simulations to demonstrate the feasibility of our triple-source saddle-curve approach. The results show that the triple-source geometry is advantageous for high temporal resolution imaging, especially important for cardiac imaging and small animal imaging.

1. Introduction

Since its introduction in 1973, x-ray CT has revolutionized clinical imaging and become a cornerstone of radiology departments. Closely correlated to the development of x-ray CT, the research for better image quality has been pursued for biomedical applications. The first dynamic CT system was the Dynamic Spatial Reconstructor (DSR) built at the Mayo Clinic in 1979 (Robb et al 1983). In a 1991 SPIE conference, we presented a spiral cone-beam scanning mode and methods to solve the long object problem (Wang et al 1991). In 1990s, single-slice spiral CT became the standard scanning mode of clinical CT (Kalender 1995). In 1998, multi-slice spiral CT entered the market (Taguchi and Aradate 1998). With the fast evolution of the source, detector and reconstruction technology, cone-beam CT is recognized as a central issue for development of the next generation clinical CT systems (Wang et al 2000). Moreover, just as there have been strong needs for clinical CT, there are strong demands and exciting results for pre-clinical CT, cone-beam micro-CT of mice and rats in particular.

Despite the impressive advancement of the CT technology, there are critical and immediate needs for better image quality in many biomedical investigations, of which we single out cardiac CT as a primary example. Earlier CT applications are in relatively static body parts such as the head, chest, abdomen and extremities, but the dynamic body parts such as the heart and specifically the coronary arteries remain challenging. Multi-row CT has made imaging of the lumen of the coronary arteries possible with excellent diagnostic accuracy. Nevertheless, imaging the lumen of the coronary arteries provides only a partial understanding of coronary atherosclerosis. In the early 1990s, it became clear that a majority of myocardial infarctions did not occur at stenotic locations. Clinical trials of thrombolytic agents in acute myocardial infarction showed that the lesion underlying the clot was often not stenotic. More recent results from the ‘Courage’ trial indicated that fixing the stenosis did not reduce death and other major cardiac events. The current research emphasis is to improve our understanding of the pathobiology and genetics of coronary artery diseases. Preliminary results demonstrated that the CT technology has good correlation with intracoronary ultrasound of plaques and can classify them into lipid containing fibrotic and calcified but with significant overlaps. Hence, cardiac CT has been identified as an important tool.

However, even with the state-of-the-art Siemens dual-source CT scanner (Flohr et al 2006), temporal resolution of 83 ms is only achievable with the aid of ECG gating around stable phases (30% or 70% of the R–R period). Major problems still exist, including (1) the incapability of offering similar image quality at an arbitrary cardiac phase for truly dynamic cardiac imaging; (2) prerequisite for breath-holding during the entire cardiac CT scan of 7–13 s, which is not feasible in the cases of pediatric and emergent scans; (3) reconstruction difficulties in the cases of high/irregular heart rates; (4) limitations with circular and helical scanning—while with a circular scan there are more image artifacts further away from the mid-plane, with a helical scan the cardiac motion is only scanned once targeting a single cardiac phase and yet a substantial number of x-ray photons emitted near the two ends of a helical trajectory segment are under-utilized; (5) failure to reach often preferred higher temporal resolution—since the mean velocity of the right coronary artery is 69.5 mm s−1, a temporal resolution of <20–80 ms is needed. In fact, most demanding cardiac imaging used to be done by electron-beam CT (EBCT) with a temporal resolution of 50–100 ms, which is expensive and rarely available, while CT of animals as human cardiac disease models may require much higher temporal resolution; for example, the heart rates of mice/rats are twice as fast or more than that of humans.

The multi-source CT architecture is an effective strategy for the improvement of temporal resolution. To date, most commercial CT scanners are based on the so-called third-generation geometry, in which a single-source detector assembly is rotated around a patient. For the dual-source scanner, the minimum rotation interval is 90° + α, compared to 180° + α in a third-generation geometry. Since the fan angle α is relatively small, there is a >40% reduction in acquisition time. As a result, the dual-source scanner reduces the data acquisition time to 0.35 s, which leads to 83 ms temporal resolution with cardiac gating. Interestingly, for exact helical cone-beam reconstruction, the trinity is better than the duality because the triple-source arrangement allows a perfect mosaic of longitudinally truncated cone-beam data to satisfy the Orlov condition (Orlov 1975) and yield better noise performance than the dual-source counterpart. In 2006, we presented the first paper on exact image reconstruction with a triple-source helical cone-beam scan (Zhao et al 2006, 2009). In the 2008 SPIE conference, we reported our work on exact image reconstruction with a triple-source saddle-curve cone-beam scan (Lu et al 2008). In the long run, we believe that the triple-source cone-beam CT scanner will play an important role for medical imaging due to its unique merits relative to the dual-source counterpart despite increased hardware cost.

This paper is organized as follows. In section 2, we overview the key formulas obtained by Pack and Noo (2005) and Yang et al (2006a, 2006b). In section 3, we propose the concept of efficient curves that leads to shift-invariant FBP reconstruction. In section 4, we study the (2N + 1)-source geometry and associated temporal resolution. In section 5, we present numerical results. Finally, in section 6, we discuss relevant issues and conclude the paper.

2. Preliminaries

In 2002, Katsevich (2002) proved the first exact and efficient reconstruction formula for helical cone-beam CT. Since then, exact reconstruction algorithms for cone-beam CT have been actively studied (Wang et al 2008). In 2005, Pack and Noo derived a general formula (Pack and Noo 2005). It can be summarized as follows.

Let x be an arbitrary point in an object, e be the filtering direction, a(λ) be a point on a scanning trajectory, [omega with right arrow above] be a unit vector in S2, λ and λ+ be the angular parameters for x, and a) and a+) be collinear. Then,


gF(λ,x,e)=ππdγ1sinγ[partial differential]g(λ,θ^(λ,x,e,γ))[partial differential]λ,




β^(λ,x,e)=e(e[center dot]a^(λ,x))a^(λ,x)e(e[center dot]a^(λ,x))a^(λ,x).

Pack and Noo (2005) proved that

κ(x,e,λ,λ+)=18π2s2dω(Rf)(ω,x[center dot]ω)σ(x,ω,e,λ,λ+),


σ(x,ω,e,λ,λ+)=12sgn(ω[center dot]e)[sgn(ω[center dot]a^(λ,x))sgn(ω[center dot]a^(λ+,x))].

This formula is attractive because R″ f is the second-order derivative of the 3D Radon transform. This means that if σ(x, [omega with right arrow above], e, λ, λ+) = 1, we have f(x) = κ(x, e, λ, λ+). If we can find trajectories with appropriate filtering directions to make σ(x, [omega with right arrow above], e, λ, λ+) a constant, the object function can be exactly reconstructed. Note that this finding is closely linked to Katsevich’s general scheme (Katsevich 2003) and the finding by Ye and Wang (2005) as a special case of Katsevich’s general scheme (Katsevich 2003).

In 2006, Yang et al (2006b) proved that an object function f can be reconstructed along a rectangle on a plane perpendicular to the rotation axis. Then, they extended the formula to a more general case (Yang et al 2006a):

Theorem 1

An object function can be reconstructed by

f(x)=12κ(x,eN,1,λN,λ1)+12i=1N1κ(x,ei,i+1,λi,λi+1),    N=3,4,5,

where ei,j=a(λj)a(λi)a(λj)a(λi) are the filtering directions and λi are points where the scanning trajectory intersects with the reconstruction plane, and the reconstructed region formed by λi is a convex polygon.

3. Efficient curves

Using equation (8), we can immediately reconstruct the image on any plane which has at least three non-collinear intersection points with the trajectory.

This method is convenient but not efficient. For different planes, we have to recalculate the intersection points, determine the filtering directions and then perform the filtration. As shown in figure 1(a), if we directly use equation (8) to reconstruct the plane Π0, although the lines la and lb on Π0 are parallel, at the vertex point a(λ), the filtration needs to be performed twice, respectively, along the lines lha and lhb on the detector plane. On the other hand, as shown in figure 1(b) if we directly use equation (8) to reconstruct the planes Π1 and Π2, at the vertex point a(λ), the filtration will only be performed along the line lh on the detector plane, where lh is the projection of both l1 and l2 (l1 and l2 are on Π1 and Π2, respectively). Thus, unless the trajectory and the reconstructed plane are appropriately selected, the algorithm will not be efficient.

Figure 1
Analysis on the computational efficiency. The detector coordinates v and u are used with (u, v) = (0, 0) being the orthogonal projection of a(λ) onto the detector plane. (a) la and lb are parallel lines on the reconstructed plane ...

In this section, we first introduce a family of curves which is called efficient, which means that shift-invariant filtration can be performed as in the classic Feldkamp-type reconstruction or Katsevich-type reconstruction. With such an efficient scanning curve, an efficient FBP algorithm can be obtained.

Our question is whether this algorithm can be efficiently applied along a general closed curve. The key is to use a shift-invariant FBP algorithm. The conditions need to be found so that the filtration can be done in a shift-invariant fashion.

Let us first introduce necessary notation. Let C be a continuous closed curve in the space and S be a family of parallel planes each of which intersects C at three points or more. S can be written as a function S(z) = {z = z′ [set membership] [zs, ze]}, zs, ze [set membership] R. When z varies from zi to zi+1, zszi < zi+1ze, if the polygons formed by the intersection points are similar to each other and the corresponding edges of the polygons are parallel to each other, we consider the curve C between zi and zi+1 efficient. The volume bounded by an efficient part of C or the corresponding portion of S(z) is denoted as V(zi, zi+1) that can be exactly and efficiently reconstructed in a shift-invariant fashion. If an ROI of the object is in the space V=i=1nV(zi,zi+1), it can be reconstructed in the shift-invariant FBP format.

Lemma 1

Let ai) (λi < λi+1, i = 1, 2, 3 …) be a series of intersection points at which S(z0) intersects C. If the curve C is efficient, we must have

a(λi)×a(λi+1)[center dot][a(λi+1)a(λi)]=0.

Proof. Consider another plane S(z0),z0z0, which intersects C at points a(λi),λi<λi+1,i=1,2,3. If C is efficient, from the definition of an efficient curve, the line connecting ai) and ai+1) is parallel to the line connecting a(λi)anda(λi+1). Moreover, points ai), ai + 1), a(λi)anda(λi+1) are coplanar. When z0 approaches z0, a(λi) moves to ai) and a(λi+1) to ai+1). These four points are still coplanar. Therefore, the line connecting a(λi) and ai), the line connecting a(λi+1) and ai+1), and the line connecting ai) and ai+1) are coplanar. Taking the limits of z0,a(λi)anda(λi+1), we have that the line connecting ai) and ai+1) and the tangent lines at ai) and ai+1) are coplanar, i.e., equation (9) is satisfied.

On the other hand, if one line does not satisfy equation (9), we can have the following four points on the curve: ai + β), ai − β), ai + 1 + β) and ai+1 − β), where β is some positive number. When β is sufficiently small, these four points are non-coplanar. Thus, for any δ1 and δ2 [set membership] [−β, β], the line through ai + δ1) and ai+1 + δ2) is not parallel to the line through ai) and ai+1). For any plane which intersects the trajectory at two points in a neighborhood of ai) and ai+1), its filtering direction (which is determined by ai + δ1) and ai+1 + δ2)) is always different from that of the plane intersecting the trajectory at ai) and ai+1). Therefore, these planes will not share the same filtering line on the detector. Since the set of these planes is of nonzero Lebesgue measure, it ruins the desirable structure for shift-invariant reconstruction. Thus, only efficient or partially efficient curves can be used in the aforementioned FBP algorithm.

Since it is not convenient to use equation (9) directly, in the following we provide two theorems to help us find the efficient curves. In doing so, we need the concept of Gaussian curvature, which is an intrinsic measure of curvature. It can be explained in terms of normal sectional curvature. Given a point on a surface and a direction in its tangent plane, the normal sectional curvature is computed by intersecting the surface with the plane determined by the point, the normal vector to the surface at that point and the direction. The normal sectional curvature is the signed curvature of this curve at that point. In this way, we can find the maximum and minimum normal sectional curvature values. Gaussian curvature is the product of these maximum and minimum values. A surface with zero Gaussian curvature is called a developable surface. It can be flattened onto a plane without any distortion.

Theorem 2

Any efficient curve must be on a series of developable surfaces.

Proof. For any continuous efficient curve, we can divide it into different parts according to the filtering directions. Suppose the curve has a filtering direction ei(λ), λ [set membership]i, λi+1]. Let a0, a1 be a pair of intersection points defining the filtering direction ei(λ), and η be the distance from a0 to a1. As shown in figure 2, we have


where λ varies from λi to λi+1, and all the filtering lines form a surface. Any point q on the surface can be denoted as

q=a0+eiμ,    μ[set membership]R.


qλ=[partial differential]q[partial differential]λ=a0,qμ=[partial differential]q(μ)[partial differential]μ=ei,n=a0×eia0×ei,qλλ=a0,qμμ=0,qλμ=0.
Figure 2
Filtering lines on an efficient curve.

Note that we write ei(λ) as ei since all filtering lines are parallel. Consequently, the Gaussian curvature equals zero at any point on this surface:

K=(qλλ[center dot]n)(qμμ[center dot]n)(qλμ[center dot]n)2(qλ[center dot]qλ)(qμ[center dot]qμ)(qλ[center dot]qμ)2=0.

Therefore, this part of the scanning trajectory must be on the developable surface (Toponogov 2006). Similarly, all parts of the curve should be on the developable surfaces. This finishes the proof.

In other words, theorem 2 gives a necessary condition for efficient curves. From the above discussion, we know that both the curve C and the set S are important. Usually, the choice of S may determine whether C is efficient or not. For the same C, various choices of S may lead to different computational efficiencies. Hence, how to choose an appropriate S becomes a major issue after the curve is known. However, we do not have a general procedure yet and have to find the optimal solution case by case.

An important family of curves with medical CT relevance can be described as


where λ [set membership] [0, 2π], H(s) is continuous and H(0) = H(2π). In this case, no matter how H(λ) changes, C(λ) is on the developable surface P(λ, z) = (R cos λ, R sin λ, z), z [set membership] R. Then, S can be specified to contain involved planes perpendicular to the z-axis.

Theorem 3

The curve C(λ) is efficient if for any minimum or maximum value of H0) at λ0, H0 − λ) = H0 + λ), λ [set membership] [0, Δλ], where Δλ is the smaller length of the two angular intervals reaching the adjacent extreme values of H0) and the range between H0 − Δλ) and H0) covers an ROI in an object in the z-direction.

Proof. Consider two points a0 − λ) and a0 + λ),



The derivatives of a0 − λ) and a0 + λ) are



Then, we have the triple product of a0 − λ), a0 + λ) and a0 + λ) − a0 − λ) as follows:

a(λ0λ)×a(λ0+λ)[center dot][a(λ0+λ)a(λ0λ)]=|Rsin(λ0λ)Rcos(λ0λ)H(λ0λ)Rsin(λ0+λ)Rcos(λ0+λ)H(λ0+λ)Rcos(λ0+λ)Rcos(λ0λ)Rsin(λ0+λ)Rsin(λ0λ)0|=R2[cos(2λ)1](H(λ0+λ)+H(λ0λ)).

Note that H(λ) = H(2λ0 − λ) and H′(λ) + H′(2λ0 − λ) = 0. Therefore, we obtain that H′0 − λ) + H′0 + λ) = 0, and the triple product always equals zero when λ [set membership] [0,Δλ]. That is, the curve C(λ) is efficient. This finishes the proof.

Theorem 3 provides a sufficient condition for the efficient curves. Using theorems 2 and 3, it can be shown that the saddle curve is efficient. Since the filtering directions are independent of the reconstruction point x, we can immediately obtain an efficient FBP algorithm characterized by the shift-invariant filtration of cone-beam projection data.

4. Analyses on FOV and temporal resolution

4.1. Field of view analysis

In the case of 2N + 1 sources, when the sources and the detector arrays are equally spaced, the radius of the field of view (FOV) is


which is the best choice to maximize the FOV (see figure 3). According to equation (20), the radius of FOV is decreased if the number of sources is increased. Hence, if we do not want to reduce FOV, we must enlarge the distance between the source and the object. This will decrease the effective x-ray flux and increase the gantry size, which is highly undesirable. With the Siemens SOMATOM volume zoom scanner as an example, the radii of the scanning circle and the FOV are 570 mm and 403 mm, respectively. Then, we can calculate the corresponding radii in a few representative cases by equation (20), as summarized in table 1. It is observed that in the triple-source geometry, the FOV (radius of 285 mm) is sufficiently large to cover a normal patient’s chest without any size increment of the scanning trajectory.

Figure 3
FOV determined by the number of the sources.
Table 1
Radii of the scanning circle and FOV for 2N + 1 sources.

4.2. Temporal resolution analysis in the case of triple sources

In this case, the three x-ray sources are symmetrically positioned along a circle, and the detectors are opposite the corresponding sources. Note that the detector area should be sufficiently large to cover all the projection lines through an object to be reconstructed. Cone-beam projections are collected when the sources and detectors are, respectively, moved along the three saddles defined by


where λ is the angular parameter, and R and h are the radius and the pitch of each saddle, respectively.

In one period of the saddles, the curves intersect at 12 points P1, P2, …, P12, as shown in figures 4 and and5.5. Let

Λ1={C1(8,10)}[union or logical sum]{C2(10,12)}[union or logical sum]{C3(12,8)},Λ2={C1(3,5)}[union or logical sum]{C2(5,1)}[union or logical sum]{C3(1,3)},Λ3={C1(11,7)}[union or logical sum]{C2(7,9)}[union or logical sum]{C3(9,11)},Λ4={C1(6,2)}[union or logical sum]{C2(2,4)}[union or logical sum]{C3(4,6)},

where Λ1, Λ2, Λ3 and Λ4 are the rebinned new trajectories, respectively, and for simplicity C1(8, 10) denotes {C1(λ)|λP8 < λ < λP10}, etc. All the rebinned curves are continuous and closed. Note that the scanning trajectories are periodic, and Λ1, Λ2, Λ3 and Λ4 will appear in turn.

Figure 4
Distribution of the 12 intersection points on the unraveled cylinder scanning surface. In the 2D view, a saddle is simply a sinusoid curve, and λT = λ mod (2π). The solid, dashed and dotted lines represent different saddles, respectively. ...
Figure 5
Rebinned curves in 3D views. The solid lines in (a), (b), (c) and (d) highlight the rebinned curves Λ1, Λ2, Λ3 and Λ4, respectively. The dotted, dash dotted with ‘+’ marks and dashed curves represent saddles ...

The intersection points can be grouped into two sets, {P1, P2, P3, P4, P5, P6} and {P7, P8, P9, P10, P11, P12}, each of which is characterized by an identical z-coordinate. The planes spanned by these sets divide the whole reconstructible volume into three parts, i.e. V1, V2 and V3, as illustrated in figure 4. Each rebinned curve covers two parts, as listed in table 2. In one period, V1 and V3 can be reconstructed twice while V2 can be reconstructed four times. Also, V2 appears in all four cases, and in one period the interval of the four phases of the object motion is the same, as shown in figure 6. In figure 6, suppose that the period is 1 s, i.e. the interval between two dashed lines is 1/12 s; then for V1 and V3, the time interval is 1/3 s with a break time of 1/6 s, and for V2, the time interval is 1/3 s without any break time. This means that with triple-source saddle CT, we can continuously perform dynamic reconstruction for V2 with three times better temporal resolution than with single-source saddle CT. Furthermore, with triple-source saddle CT, we can perform dynamic reconstruction in some specific time windows for V1 and V3 with three times better temporal resolution than with single-source saddle CT. More rebinned trajectories can be defined but the above conclusion on temporal resolution improvement remains the same.

Figure 6
Distributions of rebinned curves and reconstructible partial volumes in the case of triple sources. The horizontal axis is time. Vertical dashed lines divide each cycle into equal intervals. The gray boxes above the time axis indicate the minimum time ...
Table 2
Reconstructible partial volumes according to the different rebinned curves in the triple-source case.

Consider the plane Πz0= {(x, y, z) : z = z0, z0 [set membership] [zmin, zmax]}, which intersects any Λ [set membership]1, Λ2, Λ3, Λ4} at six points, where zmin, zmax indicate the minimum and maximum values of the saddle in the z-direction. We denote these intersection points as Fn, n = 1, 2, …, 6, and let Ωz0 be the convex polygon formed by them. By theorem 1, we have the reconstruction formula:


4.3. Temporal resolution analysis in the case of (2N + 1) sources

The main purpose of the multiple-source geometry is to achieve superior temporal resolution for dynamic imaging. Thus, it is important to analyze how the temporal resolution changes when the number of the sources increases. Recalling what we have done in the triple-source case, the original saddles are broken into many segments and then rebinned to new curves. The same strategies can be used in the case of (2N + 1) sources for N ≥ 2.

For N = 2, we have five original saddles C1, C2, C3, C4, C5 and eight rebinned continuous and closed curves Λ1, Λ2, Λ3, Λ4, Λ5, Λ6, Λ7, Λ8:

  • Λ1: {C1{(9, 30),(22, 3)}; C2{(1, 22),(24, 5)}; C3{(3, 24),(26, 7)}; C4{(5, 26),(28, 9)}; C5{(7, 28),(30, 1)}},
  • Λ2: {C1{(30, 22)}; C2{(22, 24)}; C3{(24, 26)}; C4{(26, 28)}; C5{(28, 30)}},
  • Λ3: {C1{(32, 13),(15, 36)}; C2{(34, 15),(17, 38)}; C3{(36, 17),(19, 40)}; C4{(38, 19), (11, 32)}; C5{(40, 11),(13, 34)},
  • Λ4: {C1{(13, 15)}; C2{(15, 17)}; C3{(17, 19)}; C4{(19, 11)}; C5{(11, 13)}},
  • Λ5: {C1{(4, 25),(27, 8)};C2{(6, 27),(29, 10)};C3{(8, 29),(21, 2)};C4{(10, 21),(33, 14)}; C5{(2, 23),(25, 6)}},
  • Λ6: {C1{(25, 27)}; C2{(27, 29)}; C3{(29, 21)}; C4{(21, 23)}; C5{(23, 25)}},
  • Λ7: {C1{(20, 31),(37, 18)}; C2{(12, 33),(39, 20)}; C3{(14, 35),(31, 12)}; C4{(16, 37), (33, 14)}; C5{(18, 39),(35, 16)},
  • Λ8: {C1{(18, 20)}; C2{(20, 12)}; C3{(12, 14)}; C4{(14, 16)}; C5{(16, 18)}},

where Λ1, Λ2, Λ3, Λ4, Λ5, Λ6, Λ7 and Λ8 are the rebinned new trajectories, respectively, and for simplicity, C1(9, 30) denotes {C1(λ)|λP9 < λ < λP30}, etc. Note that the scanning trajectories are periodic, and Λ123, Λ4, Λ5, Λ6, Λ7 and Λ8 will appear in turn.

The intersection points can be grouped into four sets, {P1, P2, P3, P4, P5, P6, P7, P8, P9, P10}, {P11, P12, P13, P14, P15, P16, P17, P18, P19, P20}, {P21, P22, P23, P24, P25, P26, P27, P28, P29, P30} and {P31, P32, P33, P34, P35, P36, P37, P38, P39, P40}, each of which is characterized by an identical z-coordinate. The planes spanned by these sets divide the whole reconstructible volume into five parts, i.e., V1, V2, V3, V4 and V5, as illustrated in figure 7. Similar to what we have done in the triple-source case, we can obtain the relationship between the rebinned curves and the reconstructible volume parts as listed in table 3.

Figure 7
Intersection points of the quintuple-source saddle trajectories on the unraveled cylinder scanning surface. In the 2D view, a saddle is simply a sinusoid curve, and λT = λ mod (2π). The lines with different patterns represent different ...
Table 3
Reconstructible partial volumes according to the rebinned curves in the quintuple-source case.

According to table 3, V1 and V5 can be reconstructed only twice while V2, V3 and V4 can reconstructed four times. The result is quite similar to that in the triple-source case. The only difference is that V1 here is smaller than that in the triple-source case, and V5 here is also smaller than V4 there. In fact, if we use seven sources or more, the conclusion remains the same. Thus, the increment of the source number does not improve the temporal resolution except for a small fraction of the volume.

Also, V2, V3 and V4 appear over the whole time axis, as shown in figure 8. Although in one period, all of these three parts have four phases, the time intervals for them are not the same. In figure 8, suppose that the interval between two dashed lines is 0.05 s; then for V1 and V5, the time interval is 0.2 s with break time 0.3 s, and for V2 and V4, the time intervals are 0.4 s and 0.2 s without any break time, and for V3, the time interval is 0.4 s without any break time. Remember that the shorter the time interval, the higher the temporal resolution will be. Interestingly, for V2 and V4, the time intervals are shorter in one phase but longer in the next phase than the counterparts in the triple-source case; for V3, the time intervals are longer than that in the triple-source case. In general, the temporal resolution does not become better in the quintuple-source case than in the triple-source case. In fact, V2 and V4 will become smaller and smaller when the number of source increases and V3 will be the major part of the object. Unlike in the triple-source case, the segment we select in a single curve is no longer continuous, i.e., C1{(9, 30), (22, 3)} in Λ1. We do this because the segment C1{(22, 24)} is redundant data for Λ1. However, when we calculate the time interval for a rebinned curve, such redundant segments must be taken into consideration since the x-ray source is still on for this segment. Thus, the time interval is N2N+1, which means if a single x-ray source finishes one saddle in 1 s (for a single source, one complete saddle is required to perform a fast FBP reconstruction), then it will take N2N+1S for 2N + 1 sources to acquire sufficient projections for a fast FBP reconstruction. As N2N+1 increases monotonically (with a limit 0.5), the minimum 1/3 is reached at N = 1.

Figure 8
Distributions of rebinned curves and reconstructible partial volumes in the quintuple-source case. The horizontal axis is for time. Vertical dashed lines divide each cycle into equal intervals. The gray boxes above the time axis indicate the minimum time ...

Recall that in the (2N + 1)-source helical CT case (Zhao et al 2006, 2009), the more sources, the higher the temporal resolution. Very interestingly, in the (2N + 1)-source saddle CT case, the best temporal resolution for continuous dynamic reconstruction of the central volume is achieved by triple-source saddle CT. As the number of sources is increased beyond three, the temporal resolution for continuous dynamic reconstruction of the central volume is decreased, although at the two end-parts of the object the temporal resolution within some time windows is indeed increased, and outside the time windows exact reconstruction cannot be achieved. For example, with quintuple-source saddle CT, we can continuously perform dynamic reconstruction for V2, V3, V4 with 2.5 times better temporal resolution than with single-source saddle CT, and for V1 and V5 within some time windows with five times better temporal resolution than with single-source saddle CT. Usually, we will not improve temporal resolution for a small portion of an object in an intermittent way. Therefore, the triple-source geometry is the best in our context.

5. Numerical simulation

5.1. Filtering lines



be a local coordinate system with cone-beam projection data measured on the detector plane parallel to eu(λ) and ev(λ) at a distance D from a(λ). Any position on the detector plane is specified by Cartesian coordinates (u, v), which are signed distances along eu(λ) and ev(λ), respectively, and (0, 0) denotes the orthogonal projection of a(λ) onto the detector plane. (u, v) can also be expressed as x = (xx, xy, xz). Let a = (R cos λ, R sin λ, H(λ)) be any point on the scanning trajectory and n0 = (0, 0, H(λ) + z) a point where the filtering plane Π(λ, z) intersects with the z-axis. Let n = (nx, ny, nz) be any vector in the space and e = (ex, ey, 0) be the filtering direction in the x–y plane.

Supposing that nn0 is the normal vector of the filtering plane Π(λ, z), we have

(nn0)[center dot](xn0)=0    [implies]    xxnx+xyny+(nzH(λ)z^)(xzH(λ)z^)=0.

Recall that the filtering plane Π(λ, z) passes through an0 and is parallel to the filtering direction e. Therefore,

(nn0)[center dot]e=0    [implies]    nx[center dot]ex+ny[center dot]ey=0,

(nn0)[center dot](an0)=0    [implies]    Rcosλ[center dot]nx+Rsinλ[center dot]ny=z^(nzH(λ)z^).

According to our detector setup, any point on the detector plane can be expressed as


Inserting equations (26)(28) into equation (25), we have

v(Rsinλ[center dot]exRcosλ[center dot]ey)=z^[(Dsinλucosλ)[center dot]ex(Dcosλ+usinλ)[center dot]ey],

which describes the filtering lines on the detector plane. Figure 9 shows the filtering lines on the detector plane for two source positions.

Figure 9
Filtering lines on the detector plane for two source positions. (a) λ = 0.83 π and (b) λ = 1.27 π in Λ4.

5.2. Algorithm description

In the triple-source case, the algorithm can be implemented as follows:

  • Step 1: Rearrange the saddles using equation (22).
  • Step 2: Choose one dataset and determine the filtering directions.
    For example, for the rebinned curve Λ1, the filtering directions are e1,2, e2,3, e3,4, e4,5, e5,6, e6,1, as shown in figure 10.
    Figure 10
    Illustration of the filtering directions.
  • Step 3: Compute the derivative data for every projection:
    g1(λ,u,v)=([partial differential][partial differential]λ+D(λ)u+u2+D2(λ)D(λ)[partial differential][partial differential]u+D(λ)v+uvD(λ)[partial differential][partial differential]v)gf(λ,u,v).
  • Step 4: Perform the filtration:
    • Step 4.1: Pre-weighting:
    • Step 4.2: Forward height-based rebinning:
      In table 4, λ2, λ4, λ6 are points where Λ1 reaches its maximum value along the z-axis while λ1, λ3, λ5 are points where Λ1 takes its minimum value along the z-axis.
      Table 4
      Selection of filtering directions for efficient FBP reconstruction.
    • Step 4.3: Hilbert transform in u:
      where hH is the kernel of the Hilbert transform.
    • Step 4.4: Backward height-based rebinning:
      v(Rsinλ[center dot]exRcosλ[center dot]ey)=z^[(Dsinλucosλ)[center dot]ex(Dcosλ+usinλ)[center dot]ey],
      ex and ey are the projections of the filtering direction e on the x-axis and y-axis, respectively. The filtering directions for different positions are listed in table 4.
    • Step 4.5: Post-weighting:
  • Step 5: Perform the backprojection:

For the rebinned curves Λ2, Λ3 and Λ4, the steps are essentially the same. If the object is stationary, one dataset is theoretically sufficient for exact image reconstruction. If the object changes over time, the four datasets should be used in turn for dynamic image reconstruction.

5.3. Numerical results

The performance of the proposed algorithm was numerically tested in the static and dynamic cases. As expected, the artifacts were significant when an object was dynamic. A conventional solution is to perform a motion correction in image and projection domains (Grangeat et al 2002, Li et al 2006, Yu and Wang 2007). On the other hand, the motion artifacts can be reduced if the scanning time is shortened (Flohr et al 2007, Yin and De Man 2007). Our triple-source geometry reduces the scanning time by 1/3 of that with the single-source geometry. In the dynamic case, we compared these two geometries in a numerical study.

The key simulation parameters are listed in table 5. In the static case, the Shepp–Logan phantom was used; the reconstructed slice and its representative profile in figure 11 show that our proposed algorithm performed very well. In the dynamic case, the clock phantom (Turbell and Danielsson 2000) was used. For comparison, two simulated reconstructions were performed using single-source cone-beam scanning along a saddle trajectory (Yang’s algorithm) and triple-source cone-beam scanning along helixes (Zhao’s algorithm, Zhao et al 2006), as shown in figures 12(a) and (b). The rotation speed of the x-ray source was 1 turn per second while the phantom was counterclockwise rotated by 10° per second. Figures 12(c), (d), (e) and (f) show the four reconstructed images corresponding to the time intervals 0–0.33 s, 0.25–0.58 s, 0.5–0.83 s and 0.75–1.08 s, respectively. To indicate the rotation direction of the phantom, figure 12(g) presents the difference between figures 12(f) and (c).

Figure 11
Reconstructed image of the Shepp–Logan phantom on the plane z = −15 mm and its profile along the line y = 4.9 mm. The display window is [1.0, 1.05].
Figure 12
Reconstructed images of the clock phantom on the plane z = 0. (a) The image reconstructed using Yang’s algorithm (Yang et al 2006b); (b) the image reconstructed using Zhao’s algorithm (Zhao et al 2006); (c), (d), (e) and (f) the images ...
Table 5
Parameters used in the numerical simulation.

Evidently, the images reconstructed in the triple-source geometry, whether along saddle curves or helical trajectories, were much better than those in the single-source geometry. Also, the images reconstructed using Zhao’s algorithm looked quite similar to those using our algorithm we have proposed here. Interestingly, in Zhao’s algorithm there were three backprojection segments for each point in an ROI, which were located on the different helixes and determined by three inter-helix PI-lines. In our experiment, it took less than 1/3 s to collect sufficient projections to reconstruct a single point; thus the motion artifacts were greatly reduced. However, in Zhao’s algorithm the combination of backprojection segments for each point was unique. To reconstruct the entire volume, we would need all the cone-beam projections collected along the trajectories. That is why we can only reconstruct one image from a full-scan dataset using Zhao’s algorithm. In contrast, from a full-scan dataset, four images in the four time intervals can be reconstructed using the algorithm we have proposed here.

6. Discussions and conclusions

While image artifacts are a main problem with circular cone-beam CT, saddle-curve scanning offers an attractive solution. For example, the Toshiba Aquiline One scanner utilizes 320 detector rows (0.5 mm width) to cover up to 16 cm longitudinally, which is sufficient to capture the entire heart or brain and show the physiological and pathological dynamics continuously (Rybicki et al 2008). Thus, it can greatly reduce the diagnostic imaging time and replace comprehensive exams and invasive procedures with a single scan. Unfortunately, well-known cone-beam artifacts are significant. The severity of these artifacts increases as a function of the distance from the mid-plane and may introduce substantial reconstruction errors and simulate/hide critical features. In addition to the above-discussed temporal resolution gain, our proposed cone-beam scanning along multi-saddle trajectories seems an ideal solution.

One concern with triple-source CT is that the improvement in temporal resolution may be offset by the increment in scattering artifacts. In fact, there were similar arguments against the dual-source system in comparison with the single-source counterpart prior to the introduction of the Siemens definition dual-source CT scanner. Actually, the Siemens dual-source CT scanner has received a very positive market response, which is encouraging for our proposed triple-source extension. The Siemens dual-source CT scanner has been claimed to have a dose benefit because the second detector is smaller, and the dose is reduced at the edge of the field of view. The dual-source CT scanner can double the power for cardiac applications without any compromise in image quality. It has been reported that cross scatter is smaller with the decreased collimation z-width or the decreased object size for dual-source CT (Kyriakou and Kalender 2007). Also, a potential strategy is to utilize shutters so that at any time instant the object is only exposed to one or two sources since the shutters can be rapidly opened in an alternative fashion. These guidelines should be valuable for optimization of triple-source CT. Furthermore, the recently developed interior tomography technology can significantly narrow the cone-beam aperture, thus reducing the scatter problem effectively (Ye et al 2007, Kudo et al 2008).

Another issue to be analyzed is the detector size for triple-source cone-beam CT. Briefly speaking, the necessary height of the detector remains the same as discussed by Yang et al (2006b). The width of the flat detector or the fan angle of the cylinder detector cannot be larger than 60° due to the angular limitation of the scanning mode. This may result in transverse data truncations when a patient is too large. Conventionally, the detector must cover all the projections through the patient or he/she cannot be exactly reconstructed. Again, interior tomography is being developed to allow projections to be truncated on one side or both sides (Ye et al 2007, Kudo et al 2008). Even when an ROI is entirely inside the chest, it can now be exactly reconstructed provided a priori knowledge is known on a sub-region in the ROI such as blood density in the aorta. With this new approach, triple-source CT can perform well with quite flexible data truncation to minimize radiation dose and scattering artifacts.

In conclusion, we have derived a necessary condition and a sufficient condition for construction of efficient curves so that the exact and efficient FBP algorithm can be easily designed for clinical applications. Also, we have derived a new exact and efficient algorithm for triple-source saddle-curve cone-beam CT. It reduces the scanning time, thus suppressing the motion artifacts. Finally, we have shown that triple-source saddle-curve cone-beam CT is the optimal solution in terms of FOV, hardware cost and temporal resolution among all multi-source saddle-curve cone-beam scanning modes. Therefore, we believe that triple-source saddle-curve cone-beam CT is particularly promising in reconstructing a localized dynamic volume with a beating heart as a primary example.


This work was partially supported by the National Natural Science Foundation of China (30570511 and 30770589), National High Technology Research and Development Program of China (863 Program) (2007AA02Z452), and NIH/NIBIB (EB002667 and EB004287).


  • 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]
  • Flohr T, Schoepf U, Ohnesorge B. Chasing the heart: new developments for cardiac CT. J. Thorac. Imaging. 2007;22:4. [PubMed]
  • Grangeat P, Koenig A, Rodet T, Bonnet S. Theoretical framework for a dynamic cone-beam reconstruction algorithm based on a dynamic particle model. Phys. Med. Biol. 2002;47:2611–2625. [PubMed]
  • Kalender WA. Thin-section three-dimensional spiral CT: is isotropic imaging possible? Radiology. 1995;197:578–580. [PubMed]
  • Katsevich A. Theoretically exact filtered backprojection-type inversion algorithm for spiral CT. SIAM J. Appl. Math. 2002;62:2012.
  • Katsevich A. A general scheme for constructing inversion algorithms for cone beam CT. Int. J. Math. Math. Sci. 2003;21:1305–1321.
  • 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]
  • Kyriakou Y, Kalender W. Intensity distribution and impact of scatter for dual-source CT. Phys. Med. Biol. 2007;52:6969. [PubMed]
  • Li T, Schreibmann E, Yang Y, Xing L. Motion correction for improved target localization with on-board cone-beam computed tomography. Phys. Med. Biol. 2006;51:253. [PubMed]
  • Lu Y, Zhao J, Wang G. Exact image reconstruction for triple-source cone-beam CT along saddle trajectories. Proc. SPIE. 2008 70780K.
  • Orlov SS. Theory of three-dimensional reconstruction: 1. Conditions of a complete set of projections. Sov. Phys. Crystallogr. 1975;20:312–314.
  • Pack JD, Noo F. Cone-beam reconstruction using 1D filtering along the projection of M-lines. Inverse Problems. 2005;21:1105.
  • Robb RA, Hoffman EA, Sinak LJ, Harris LD, Ritman EL. High-speed three-dimensional x-ray computed tomography: the dynamic spatial reconstructor. Proc. IEEE. 1983;71:308–319.
  • Rybicki FJ, Otero HJ, Steigner ML, Vorobiof G, Nallamshetty L, Mitsouras D, Ersoy H, Mather RT, Judy PF, Cai T. Initial evaluation of coronary images from 320-detector row computed tomography. Int. J. Cardiovasc. Imaging (formerly Card. Imaging) 2008;24:535–546. [PubMed]
  • Taguchi K, Aradate H. Algorithm for image reconstruction in multi-slice helical CT. Med. Phys. 1998;25:550. [PubMed]
  • Toponogov V. Differential Geometry of Curves and Surfaces: A Concise Guide. Boston, MA: Birkhauser; 2006.
  • Turbell H, Danielsson PE. Helical cone-beam tomography. Int. J. Imaging Syst. Technol. 2000;11:91–100.
  • Wang G, Crawford CR, Kalender WA. Guest editorial—Multirow detector and cone-beam spiral/helical CT. IEEE Trans. Med. Imaging. 2000;19:817–821. [PubMed]
  • Wang G, Lin TH, Cheng PC, Shinozaki DM, Kim HG. Scanning cone-beam reconstruction algorithms for x-ray microtomography. Proc. SPIE. 1991;1556:99.
  • Wang G, Yu H, De Man B. An outlook on x-ray CT research and development. Med. Phys. 2008;35:1051. [PubMed]
  • Yang H, Li M, Koizumi K, Kudo H. Application of pack and Noo’s cone-beam inversion formula to a wide class of trajectories; IEEE Nucl. Sci. Symp. Conf. Record; 2006a.
  • Yang H, Li M, Koizumi K, Kudo H. Exact cone beam reconstruction for a saddle trajectory. Phys. Med. Biol. 2006b;51:1157. [PubMed]
  • Yin Z, De Man B. 3D analytic cone-beam reconstruction for less than a full scan using a multi-source CT system; IEEE Nucl. Sci. Symp. Conf. Record; 2007.
  • Ye Y, Wang G. Filtered backprojection formula for exact image reconstruction from cone-beam data along a general scanning curve. Med. Phys. 2005;32:42. [PubMed]
  • Ye Y, Yu H, Wei Y, Wang G. A general local reconstruction approach based on a truncated Hilbert transform. Int. J. Biomed. Imaging. 2007;2007:63634. [PMC free article] [PubMed]
  • Yu H, Wang G. Data consistency based rigid motion artifact reduction in fan-beam CT. IEEE Trans. Med. Imaging. 2007;26:249. [PubMed]
  • 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.
  • Zhao J, Jin Y, Lu Y, Wang G. A filtered backprojection algorithm for triple-source helical cone-beam CT. IEEE Trans. Med. Imaging. 2009;28:384–393. [PMC free article] [PubMed]