PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Tsinghua Sci Technol. Author manuscript; available in PMC 2010 June 16.
Published in final edited form as:
Tsinghua Sci Technol. 2010 February; 15(1): 56–61.
doi:  10.1016/S1007-0214(10)70009-2
PMCID: PMC2886312
NIHMSID: NIHMS183218

Consistency Conditions for Cone-Beam CT Data Acquired with a Straight-Line Source Trajectory*

Abstract

A consistency condition is developed for computed tomography (CT) projection data acquired from a straight-line X-ray source trajectory. The condition states that integrals of normalized projection data along detector lines parallel to the X-ray path must be equal. The projection data is required to be untruncated only along the detector lines parallel to the X-ray path, a less restrictive requirement compared to Fourier conditions that necessitate completely untruncated data. The condition is implemented numerically on simple image functions, a discretization error bound is estimated, and detection of motion inconsistencies is demonstrated. The results show that the consistency condition may be used to quantitatively compare the quality of projection data sets obtained from different scans of the same image object.

Keywords: data redundancy, John’s equation, consistency conditions

Introduction

Computed tomography (CT) is the extension of X-ray-based medical imaging from the familiar two-dimensional (2-D) shadowgrams to full three-dimensional (3-D) reconstructions of a patient’s internal structure. CT scanners measure the approximate line integrals of a patient’s linear attenuation coefficient along the X-ray path. Reconstructing an image from CT data requires solving the following inverse problem: given a set of measured line integrals, estimate the patient’s attenuation distribution function. In reality there are several setbacks to determining the exact image function. First, the data collected from a real CT system is not continuous, thus real CT reconstruction involves solving the discrete version of a continuous model equation which introduces unavoidable discretization errors. Second, the continuous model is an ideal model (idealization) that depends on perfect projection data. However real projection data is far from perfect; the data is wrought with inconsistencies due to practical considerations such as X-ray scatter, noise, beam hardening effects, patient motion, and machine calibration.

The aim of this article is to develop quantitative methods for evaluating the quality of projection data by deriving consistency conditions on the projection data so that inconsistencies may be systematically reduced prior to reconstruction. Conditions for 2-D parallel- and fan-beam CT geometries are well-known but those for cone-beam CT and other more practical configurations are less understood. In general redundancies in the projection data can be expressed by an ultrahyperbolic partial differential equation derived in 1938 by the mathematician Fritz John. Working years before the development of tomography, John observed that the X-ray integral naturally depends on four independent variables while the image function itself relies only on three, so that the reconstruction problem must contain redundancies[1]. Recent works have shown that John’s equation may integrated to determine cone-beam projections for certain point sources located near the X-ray trajectory[2,3]. Exact reconstruction may then be more readily achieved by first obtaining a new set of cone-beam projection data. John’s equation has been also applied to rebinning algorithms for helical cone-beam CT and has been extended beyond CT to include derivation of an exact rebinning algorithm for 3-D positron emission tomography (PET)[4]. However, the application of John’s equation to improve quality of known projection data has not been sufficiently explored.

1 Consistency Conditions and John’s Equation

Let fC02 be a real image function with compact support. The X-ray transform of f is defined as

g(ξ,η)=Rf(ξ+t(ηξ))dt
(1)

where ξ and η denote locations on the X-ray source trajectory and detector plane, respectively, and g is the measured projection data. The normalized projections

g~(ξ,η)=1ηξg(ξ,η)
(2)

are the measured projections scaled by [mid ] ηξ [mid ], the length of a line segment beginning at ξ on the source trajectory and ending at η in the detector plane.

We consider a mathematically convenient setup in which the X-ray trajectory lies in a plane parallel to and unit distance above the detector plane. As shown in Fig. 1 the projection data g (as well as the normalized data g~) depends on four independent variables, ζ1,ξ2,η1, and η2. The image function f can be at most a function of three independent spatial variables x, y, and z, so that Eq. (1) necessarily contains redundancies. The ultrahyperbolic partial differential equation,

(2η1ξ22η2ξ1)g~(ξ,η)=0
(3)

known as John’s equation, expresses these redundancies as a consistency condition, and the solution of Eq. (3) under appropriate boundary conditions yields a solution to the integral problem (1).

Fig. 1
Geometrical illustration of redundancies in the X-ray transform Eq. (1). The path from a source ξ = (ξ1,ξ2,1) to a point on the detector η = (η1,η2,0) is parametrized by four independent variables while ...

Beginning with Eq. (3) we take the Fourier transform of g~ with respect to the detector variable η,

R2(2η1ξ22η2ξ1)g~(ξ,η)e2πiηκdη=0
(4)

where κ is the wave vector associated with the detector variable η. Integrating Eq. (4) by parts gives

κξg^(ξ,κ)=0
(5)

a wave equation with solution

g^(ξ,κ)=g^(ξ+τκ,κ)τR
(6)

The known boundary data are the normalized projections g~(ξ,η) acquired from source points ξ on the X-ray source trajectory. This data may be used together with Eq. (6) to compute certain transformed normalized projections g~(ξ,κ) from source points not on the X-ray source trajectory. As shown in Fig. 2 when the X-ray source moves along a circular trajectory [mid ]ξ[mid ]2=1, one can compute g~(ξ,κ) for all [mid ]ξ [mid ]2<1 and build up the two-dimensional partial Fourier space with the one-dimensional transformed projections. Then the unmeasured normalized projections inside the source trajectory can be recovered by inverting the two-dimensional transform[2].

Fig. 2
The partial Fourier transform of the normalized projections with respect to the detector variable can be computed for all points inside the circular X-ray source trajectory[2].

In this article our goal is not to determine unmeasured projection data; rather, we aim to use Eq. (6) to construct a metric to evaluate the quality of known projection data. We consider the situation where the source moves along a line; without loss of generality, the source ξ = (ξ,0) moves in the κ[perpendicular] = (κ,0) direction. Equation (6) then reads

g^(ξ,κ)=g^(ξ+t,κ)tR
(7)

where t = κτ, and after inverting the partial Fourier transform, the consistency condition becomes

G=g~(ξ,η)dη=g~(ξ+t,η)dηtR
(8)

We therefore conclude from Eq. (8) that for any two point sources on a linear X-ray trajectory, the integrals of the normalized projections acquired from the two sources over the detector must be equal. The advantage condition (8) has over (6) is that (8) only requires untruncated projection data for lines on the detector parallel to the linear X-ray trajectory; condition (6) requires completely untruncated data everywhere to compute the Fourier transform. These conclusions are illustrated in Fig. 3. Condition (8) is somewhat reminiscent of Grangeat’s relation between cone-beam data and the radial derivative of the 3-D Radon transform. In fact, Grangeat’s relation can also be used to provide complimentary information to the condition developed here[5]. It is also worth noting that Eq. (8) is not equivalent to the parallel beam condition equating the integrals of projection data over the detector for each view. The latter is an equal area condition which is not the case here, as the data is weighted by the X-ray path length.

Fig. 3
For any two sources ξ and ξ + t on a linear X-ray source trajectory, the integral of the normalized measured projections g~ over a parallel line in the detector plane must be equal. Equation (8) may be used when the projection data is ...

Condition (8) was derived by assuming that two sources lie on a linear X-ray trajectory. This condition must be extended to general X-ray trajectories before it can be applied to practical circular and helical cone-beam scanning configurations. To do so, we let the linear trajectory in the earlier derivation be the intersource line joining two sources on a general trajectory. The integration in condition (8) is then performed along any projection of the intersource line onto the detector. The projection line is defined by a point x0 and a direction vector n^1 and is determined from the intersection of a plane containing the two sources and a virtual detector plane, a plane parallel to the intersource line direction vector. For a given pair of sources ξ1 and ξ2 lying on a general X-ray trajectory, there exist an infinite number of planes containing both. As shown in Fig. 4, a particular source plane may be chosen by specifying the plane normal vector n^s. Consider one of these source planes, and a virtual detector plane with normal vector n^d and which contains the point d0 . Provided they are not parallel (n^dkn^s, k[set membership]R), the source and virtual detector planes intersect in a line with direction vector n^1=n^d×n^s. The point x0 on the intersection line is found by solving linear system Mx0 = b where

M=[n^sn^dn^l]andb=[ξin^sd0n^dd0n^l]
(9)

Here M is always invertible. After finding x0 and n^1 for a given source plane, the consistency condition (8) may then be applied. Note that although the data g~(ξ,η) is acquired on a real detector, the integration in Eq. (8) is performed along lines on the virtual detector.

Fig. 4
The intersection of a plane containing two sources ξ1 and ξ2 and the virtual detector is a line with direction n^s×n^d , parallel to the intersource line ξ2ξ1 (solid lines). The source plane is not unique. ...

2 Numerical Implementation

The integral condition (8) is a quantitative metric for evaluating the quality of CT measured projection data. For perfect CT projection data,

Q=g~(ξ,η)dηg~(ξ+t,η)dη
(10)

is identically zero, but in practical situations the data contain inconsistencies and Q is nonzero. Furthermore, Eq. (8) is a continuous condition and must be discretized before implementation as real data is discrete and fixed by the number of bins on the detector. Discretization errors are thus unavoidable and are introduced independent from the quality of the projection data. Condition (8) is then only useful if one can distinguish between data inconsistencies (noise, systematic drift, etc.) and discretization errors.

Let δ be a bound for the discretization errors that result from numerical implementation of Eq. (8). Then Qδ if the projection data is perfect or contain inconsistencies with sufficiently small total error. An estimate for the error bound δ can be determined quantitatively with perfect projection data by computing Q numerically for a realistic number of equally sized detector bins. Consider for example an object composed of two squares with respective constant attenuations f0 and f1, side lengths 2a0 and 2a1, and perpendicular distances h0 and h1 from the intersource line as shown in Fig. 5. The analytical value of Eq. (8) for an individual square is computed by rewriting

g~(ξ,η)dη=01f[ξ+t(ηξ)]dtdη01f(tη,td)dtdη=0df(x,y)ydydx=2af0lnh0+a0h0a0
(11)

using transformation x = , z = td, and for an object composed of n squares we then have

g~(ξ,η)dη=2i=0naifilnhi+aihiai
(12)
Fig. 5
Geometrical setup employed to rewrite Eq. (8) in an analytically tractable form. The result for a square with constant image function f (x, y) = f0, side 2a0, and center distance h0 from the intersource line is given by Eq. (11).

For the two sources we numerically integrate the corresponding normalized projections over the detector and plot the results as a function of the number of detector bins N. One can see in Fig. 6 that as the number of bins increases, the numerical integration converges to the analytical value of the integral indicated by the red line. In particular, when there are a realistic number of bins ( N = 672 ), the percent difference Q is on the order of 10−3 and thus δ ~ 0.001%, shown in Fig. 7. To measure the quality of the object’s real CT projection data, one would compute Q for this real data and compare the result to δ ~ 0.001%. If Q > δ we conclude that the real projection data has measurable inconsistencies and we rank the quality of projection data from different scans by the amount that Q deviates from δ. However when Qδ we cannot conclude that the projection data is free of inconsistencies; we can only conclude that the errors from the inconsistencies cannot be resolved due to the presence of larger discretization errors.

Fig. 6
Numerical values of G defined by Eq. (8) as a function of detector bin number N for a source (top) lying on a circular trajectory and separated from a second source (bottom) by polar angle 90°. The integration is for midplane projections of a ...
Fig. 7
Top: percent difference Q as a function of detector bin number N for midplane projections for two sources separated by polar angle 90°. Bottom: discretization error is on the order of 10−3% for N =672.

3 Detection of Motion Inconsistencies

One application of the consistency condition is in detecting motion inconsistencies that occur in many tomographic systems such as image-guided radiation therapy (IGRT), C-arm CT, breast tomosynthesis and CT, and cardiac CT. Motivated by the cardiac CT application, we studied motion detection for a circular cone-beam scan by comparing inconsistencies between views from two simulated heart phantoms, one stationary and one beating. The projection data from the stationary phantom was used as a baseline reference for which to compare the projection data from the dynamic phantom.

Both phantoms were composed of an ellipsoidal chest with constant attenuation f =1 and semi-principal axes ax =20, ay = az =10. In the stationary phantom, the heart was modelled by a sphere with constant attenuation f = 0.05 and radius r = 5, and the chest and heart were centered at the origin of the coordinate system. As shown in Fig. 8, the heart in the dynamic phantom varied sinusoidally in both radius and location with respect to the chest center. All simulations were performed using an N = 512 bin detector. Figure 9 shows Q as a function of the polar angle between views for both midplane and off-mid-plane projection data. For both projection sets the motion inconsistencies are clearly detected and vary sinusoidally with frequency ω. It can also be seen that the inconsistency measure Q approaches a minimum for views at the same place in the simplied cardic cycle. The nonzero Q for the consistent static data is a result of detector discretization since our detector only had 512 bins. If we increased the number of bins to 1024, the measure Q would be closer to zero, and only in the continuous limit N → ∞ would we have Q = 0. The important point here is that motion inconsistencies in the dynamic heart are distinguishable despite using a detector with only 512 bins. This measure may prove particularly useful for irregular heartbeats which have the tendency to complicate the standard echocardiographic (ECG) gated-CT.

Fig. 8
Dynamic spherical heart phantom embedded in a stationary ellipsoidal chest. The heart attenuation was 5% larger than that of the background chest. Both the heart radius and center varied sinusoidally with frequency ω = 4 from r = 4 and x0 = −2 ...
Fig. 9
Q as a function of polar angle between views Δθ for the midplane (top) and an off-midplane (bottom) projection for both static (blue) and dynamic (red) heart phantoms. The angle between the two planes is 3°. The nonzero Q for the ...

Condition (8) is the only method currently available for testing motion inconsistencies between views. For a stationary subject it is possible to rebin the projection data to parallel beam and then apply conjugate ray symmetry and moment conditions to assess inconsistencies between views. However, this process becomes problematic for dynamic subjects. In this case the rays in each rebinned parallel view are measured at different times and one could not expect the parallel consistency conditions to hold. Furthermore applying conjugate ray and moment conditions after rebinning only works for the midplane projections, whereas Eq. (8) applies to other planes such as the one represented in the bottom of Fig. 8 where the angle is 3° between this plane and the midplane. Thus Eq. (8) allows for the possibility of combining the Q from different planes.

4 Conclusions

This study shows that consistency conditions based on John’s equation may be used to quantitatively compare the quality of projection data sets obtained from different scans of the same image object. The projection data may be partially truncated; the only restriction is for the data to be untruncated on detector lines parallel to the X-ray path. The condition may prove useful for system calibration as well as in assessing quality when the data contains inconsistencies due to the presence of noise, X-ray beam polychromaticity and scatter, and metal. It is applicable to practical circular and helical cone-beam scanning configurations, and we have demonstrated the utility of the condition for detecting motion inconsistencies by implementing it on static and dynamic heart phantom data.

Acknowledgements

The authors would like to thank Thomas Köhler for discussion of Grangeat’s relation as a consistency condition at the 2009 fully 3-D meeting.

Footnotes

*Supported in part by NIH R01 (Nos. CA120540 and EB000225) and the Illinois Department of Public Health Ticket for the Cure Grant. E.Y. Sidky was supported in part by a Career Development Award from NIH SPORE (No. CA125183-03)

References

[1] John F. The ultrahyperbolic equation with 4 independent variables. Duke Math. J. 1938;4:300–322.
[2] Patch SK. Consistency conditions upon 3D CT data and the wave equation. Phys. Med. Biol. 2002;47:2637–2650. [PubMed]
[3] Finch D. Cone beam reconstruction with sources on a curve. SIAM J. Appl. Math. 1985;45:665–673.
[4] Defrise M, Liu X. Fast rebinning algorithms for 3D positron emission tomography using John’s equation. Inverse Probl. 1999;15:1047–1065.
[5] Grangeat P. Mathematical framework of cone-beam 3D reconstruction via the first derivative of the radon transform. Mathematical Methods in Tomography. 1991;1497:66–97.