Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2886312

Formats

Article sections

- Abstract
- Introduction
- 1 Consistency Conditions and John’s Equation
- 2 Numerical Implementation
- 3 Detection of Motion Inconsistencies
- 4 Conclusions
- References

Authors

Related links

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-2PMCID: PMC2886312

NIHMSID: NIHMS183218

Department of Radiology, University of Chicago, 5841 S. Maryland Ave., Chicago, IL 60637, USA

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.

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.

Let $f\in {C}_{0}^{2}$ be a real image function with compact support. The X-ray transform of *f* is defined as

$$g(\xi ,\eta )={\int}_{\mathbf{R}}f(\xi +t(\eta -\xi \left)\right)\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}t$$

(1)

where ** ξ** and

$$\stackrel{~}{g}(\xi ,\eta )=\frac{1}{\mid \eta -\xi \mid}g(\xi ,\eta )$$

(2)

are the measured projections scaled by ** η** –

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 $\stackrel{~}{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,

$$\left(\frac{{\partial}^{2}}{\partial {\eta}_{1}\partial {\xi}_{2}}-\frac{{\partial}^{2}}{\partial {\eta}_{2}\partial {\xi}_{1}}\right)\stackrel{~}{g}(\xi ,\eta )=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).

Beginning with Eq. (3) we take the Fourier transform of $\stackrel{~}{g}$ with respect to the detector variable ** η**,

$${\int}_{{\mathbf{R}}^{2}}\left(\frac{{\partial}^{2}}{\partial {\eta}_{1}\partial {\xi}_{2}}-\frac{{\partial}^{2}}{\partial {\eta}_{2}\partial {\xi}_{1}}\right)\stackrel{~}{g}(\xi ,\eta ){\mathrm{e}}^{-2\pi \mathrm{i}\eta \cdot \kappa}\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}\eta =0$$

(4)

where ** κ** is the wave vector associated with the detector variable

$${\kappa}^{\perp}\cdot {\nabla}_{\xi}\phantom{\rule{thinmathspace}{0ex}}\widehat{g}(\xi ,\kappa )=0$$

(5)

a wave equation with solution

$$\widehat{g}(\xi ,\kappa )=\widehat{g}(\xi +{{\tau}_{\kappa}}^{\perp},\kappa )\phantom{\rule{1em}{0ex}}\forall \tau \in \mathbf{R}$$

(6)

The known boundary data are the normalized projections $\stackrel{~}{g}(\xi ,\eta )$ 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 $\stackrel{~}{g}(\xi ,\kappa )$ 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

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 ** ξ** = (

$$\widehat{g}(\xi ,\kappa )=\widehat{g}(\xi +t,\kappa )\phantom{\rule{1em}{0ex}}\forall t\in \mathbf{R}$$

(7)

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

$$G={\int}_{-\infty}^{\infty}\stackrel{~}{g}(\xi ,\eta )\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}\eta ={\int}_{-\infty}^{\infty}\stackrel{~}{g}(\xi +t,\eta )\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}\eta \phantom{\rule{1em}{0ex}}\forall t\in \mathbf{R}$$

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

For any two sources *ξ* and *ξ* + *t* on a linear X-ray source trajectory, the integral of the normalized measured projections $\stackrel{~}{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 *x*_{0} and a direction vector ${\widehat{\mathit{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 ${\widehat{n}}_{\mathrm{s}}$. Consider one of these source planes, and a virtual detector plane with normal vector ${\widehat{\mathit{n}}}_{\mathrm{d}}$ and which contains the point *d*_{0} . Provided they are not parallel (${\widehat{\mathit{n}}}_{\mathrm{d}}\ne k{\widehat{\mathit{n}}}_{\mathrm{s}}$, *k***R**), the source and virtual detector planes intersect in a line with direction vector ${\widehat{\mathit{n}}}_{1}={\widehat{\mathit{n}}}_{\mathrm{d}}\times {\widehat{\mathit{n}}}_{\mathrm{s}}$. The point *x*_{0} on the intersection line is found by solving linear system *Mx*_{0} = ** b** where

$$\mathit{M}=\left[\begin{array}{c}\hfill {\widehat{\mathit{n}}}_{\mathrm{s}}\hfill \\ \hfill {\widehat{\mathit{n}}}_{\mathrm{d}}\hfill \\ \hfill {\widehat{\mathit{n}}}_{\mathrm{l}}\hfill \end{array}\right]\text{and}\phantom{\rule{1em}{0ex}}\mathit{b}=\left[\begin{array}{c}\hfill {\xi}_{i}\cdot {\widehat{\mathit{n}}}_{\mathrm{s}}\hfill \\ \hfill {\mathit{d}}_{0}\cdot {\widehat{\mathit{n}}}_{\mathrm{d}}\hfill \\ \hfill {\mathit{d}}_{0}\cdot {\widehat{\mathit{n}}}_{\mathrm{l}}\hfill \end{array}\right]$$

(9)

Here ** M** is always invertible. After finding

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

$$Q=\mid {\int}_{-\infty}^{\infty}\stackrel{~}{g}(\xi ,\eta )\mathrm{d}\eta -{\int}_{-\infty}^{\infty}\stackrel{~}{g}(\xi +t,\eta )\mathrm{d}\eta \mid $$

(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 *f*_{0} and *f*_{1}, side lengths 2*a*_{0} and 2*a*_{1}, and perpendicular distances *h*_{0} and *h*_{1} from the intersource line as shown in Fig. 5. The analytical value of Eq. (8) for an individual square is computed by rewriting

$$\begin{array}{c}\hfill {\int}_{-\infty}^{\infty}\stackrel{~}{g}(\xi ,\eta )\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}\eta ={\int}_{-\infty}^{\infty}{\int}_{0}^{1}f[\xi +t(\eta -\xi \left)\right]\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}t\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}\eta \hfill \\ \hfill {\int}_{-\infty}^{\infty}{\int}_{0}^{1}f(t\eta ,td)\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}t\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}\eta ={\int}_{-\infty}^{\infty}{\int}_{0}^{d}\frac{f(x,y)}{y}\mathrm{d}y\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}x=\hfill \\ \hfill 2a{f}_{0}\phantom{\rule{thinmathspace}{0ex}}\mathrm{ln}\mid \frac{{h}_{0}+{a}_{0}}{{h}_{0}-{a}_{0}}\mid \hfill \end{array}$$

(11)

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

$${\int}_{-\infty}^{\infty}\stackrel{~}{g}(\xi ,\eta )\phantom{\rule{thinmathspace}{0ex}}\mathrm{d}\eta =2\sum _{i=0}^{n}{a}_{i}{f}_{i}\phantom{\rule{thinmathspace}{0ex}}\mathrm{ln}\mid \frac{{h}_{i}+{a}_{i}}{{h}_{i}-{a}_{i}}\mid $$

(12)

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.

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 *a _{x}* =20,

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 *x*_{0} = −2 **...**

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.

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.

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.

^{*}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)

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

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |