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

**|**HHS Author Manuscripts**|**PMC2840998

Formats

Article sections

- Abstract
- I. INTRODUCTION
- II. CONE BEAM GEOMETRY AND PREVIOUS CALIBRATION METHODS
- III. THEORY
- IV. COMPUTER SIMULATION
- V. IMPLEMENTATION
- VI. DISCUSSION
- VII. CONCLUSIONS
- References

Authors

Related links

Med Phys. Author manuscript; available in PMC 2010 March 17.

Published in final edited form as:

Med Phys. 2006 June; 33(6): 1695–1706.

PMCID: PMC2840998

NIHMSID: NIHMS178801

Kai Yang, Department of Radiology, University of California, Davis Medical Center, 4701 X Street, Sacramento, California 95817 and Department of Biomedical Engineering, University of California, Davis, California 95616;

See other articles in PMC that cite the published article.

Cone beam CT systems are being deployed in large numbers for small animal imaging, dental imaging, and other specialty applications. A new high-precision method for cone beam CT system calibration is presented in this paper. It uses multiple projection images acquired from rotating point-like objects (metal ball bearings) and the angle information generated from the rotating gantry system is also used. It is assumed that the whole system has a mechanically stable rotation center and that the detector does not have severe out-of-plane rotation (< 2°). Simple geometrical relationships between the orbital paths of individual BBs and five system parameters were derived. Computer simulations were employed to validate the accuracy of this method in the presence of noise. Equal or higher accuracy was achieved compared with previous methods. This method was implemented for the geometrical calibration of both a micro CT scanner and a breast CT scanner. The reconstructed tomographic images demonstrated that the proposed method is robust and easy to implement with high precision.

Flat panel detector-based cone beam CT systems acquire projection images during the rotation of either the object itself, or the x-ray and detector systems around the object on a mechanically stable axis. The axis of rotation is referred to as the rotation center. A three-dimensional volume data set consisting of a number of tomographic images of the object is reconstructed from the projection data. Precise assessment of the CT system’s geometric parameters is crucial to achieve successful reconstruction with good spatial resolution and low artifact content in the reconstructed tomographic images. In this study, a method for estimating system geometric parameters necessary for cone beam reconstruction is discussed. An accurate determination of system geometric parameters results from the proposed method. The cone beam CT system calibration method utilizes a ball-bearing phantom (BB phantom) and angle information generated from the gantry encoder electronics. Because this method does not require precise information about the phantom (only a rough measurement of distance between two BBs is required), it is robust and easy to implement.

To describe the cone beam CT geometry, it is convenient to assume that the x-ray source and detector system are stationary and that the object rotates around the rotation center. As shown in Fig. 1, the rotation axis is defined as the *z* axis of the system. The axis which passes through the cone vertex (x-ray tube focal spot) and which is also perpendicular to the *z* axis is defined as the *x* axis. The axis perpendicular to the *x*-*z* plane which passes through the intersection of the *x* axis and *z* axis is defined as the *y* axis. In the detector plane, two more axes are defined along the detector, (*u*,υ) for horizontal and vertical, respectively. Thus, (*u*,υ)=(0,0) represents the top left (viewed from the x-ray source) detector pixel.

If the detector is aligned such that the υ axis is parallel to *z* axis and the *u* axis is parallel to *y* axis, the system geometric parameters include:

- Source to isocenter distance,
*R*, distance from the cone vertex to the rotation center, and thus the coordinate of the cone vertex (x-ray tube focal spot) is (−_{FI}*R*,0,0)._{FI} - Source to detector distance,
*R*, distance from the cone vertex to the detector plane, so the detector plane is located at_{FD}*x*=*R*−_{FD}*R*._{FI} *u*_{0}, horizontal location of the intersection of the*x*axis and detector plane.- υ
_{0}, vertical location of the intersection of the*x*axis and detector plane.The detector rotation is defined in three parameters, as shown in Fig. 2. ϕ and σ are out-of-plane rotation angles and η is the in-plane rotation angle. - ϕ, the rotation angle of the detector plane along the axis of
*u*=*u*_{0}, which is also the axis determined by*y*=0 and*x*=*R*−_{FD}*R*._{FI} - σ, the rotation angle of the detector plane along the axis of υ= υ
_{0}, which is also the axis determined by*z*=0 and*x*=*R*−_{FD}*R*._{FI} - η, the rotation angle of the detector plane along the point of (
*u*_{0},υ_{0}), which is also the point of (*R*−_{FD}*R*,0,0)._{FI}

For a position (*x _{i}*,

$$\begin{array}{c}{u}_{i}=\frac{{R}_{\mathit{\text{FD}}}}{{R}_{\mathit{\text{FI}}}+{x}_{i}}\xb7\frac{{y}_{i}}{\mathrm{\Delta}u}+{u}_{0},\hfill \\ {\upsilon}_{i}=\frac{{R}_{\mathit{\text{FD}}}}{{R}_{\mathit{\text{FI}}}+{x}_{i}}\xb7\frac{{z}_{i}}{\mathrm{\Delta}\upsilon}+{\upsilon}_{0}.\hfill \end{array}$$

(1)

For a non-ideally aligned detector with three rotation angles known as (ϕ,σ,η), a transform matrix *T* can be used to define the ideally aligned detector position. We use (*u*′,υ′,*w*′) to represent a point on the nonideal detector plane, and (*u*,υ,*w*) to represent the corresponding position on the ideal detector plane. The origin of the rotation in (*u*,υ,*w*) system is located at position (*u*_{0}·Δ*u*,υ_{0}·Δυ,*R _{FD}*−

Then we have

$$\left[\begin{array}{c}(u-{u}_{0})\xb7\mathrm{\Delta}u\\ (\upsilon -{\upsilon}_{0})\xb7\mathrm{\Delta}\upsilon \\ w-({R}_{\mathit{\text{FD}}}-{R}_{\mathit{\text{FI}}})\end{array}\right]=T\xb7\left[\begin{array}{c}(u\prime -{u}_{0})\xb7\mathrm{\Delta}u\\ (\upsilon \prime -{\upsilon}_{0})\xb7\mathrm{\Delta}\upsilon \\ w\prime -({R}_{\mathit{\text{FD}}}-{R}_{\mathit{\text{FI}}})\end{array}\right]\phantom{\rule{thinmathspace}{0ex}},$$

(2)

where

$$T=\left[\begin{array}{ccc}\text{cos}\phantom{\rule{thinmathspace}{0ex}}\sigma \phantom{\rule{thinmathspace}{0ex}}\text{cos}\phantom{\rule{thinmathspace}{0ex}}\eta -\text{sin}\phantom{\rule{thinmathspace}{0ex}}\sigma \phantom{\rule{thinmathspace}{0ex}}\text{sin}\varphi \phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}\eta & -\phantom{\rule{thinmathspace}{0ex}}\text{cos}\phantom{\rule{thinmathspace}{0ex}}\sigma \phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}\eta -\text{sin}\phantom{\rule{thinmathspace}{0ex}}\sigma \phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}\varphi \phantom{\rule{thinmathspace}{0ex}}\text{cos}\phantom{\rule{thinmathspace}{0ex}}\eta & -\phantom{\rule{thinmathspace}{0ex}}\mathrm{sin}\phantom{\rule{thinmathspace}{0ex}}\sigma \phantom{\rule{thinmathspace}{0ex}}\text{cos}\phantom{\rule{thinmathspace}{0ex}}\varphi \\ \text{cos}\phantom{\rule{thinmathspace}{0ex}}\varphi \phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}\eta & \text{cos}\phantom{\rule{thinmathspace}{0ex}}\varphi \phantom{\rule{thinmathspace}{0ex}}\text{cos}\phantom{\rule{thinmathspace}{0ex}}\eta & -\phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}\varphi \\ \text{sin}\phantom{\rule{thinmathspace}{0ex}}\sigma \phantom{\rule{thinmathspace}{0ex}}\text{cos}\phantom{\rule{thinmathspace}{0ex}}\eta \phantom{\rule{thinmathspace}{0ex}}+\phantom{\rule{thinmathspace}{0ex}}\text{cos}\phantom{\rule{thinmathspace}{0ex}}\sigma \phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}\varphi \phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}\eta & -\phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}\sigma \phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}\eta \phantom{\rule{thinmathspace}{0ex}}+\phantom{\rule{thinmathspace}{0ex}}\text{cos}\phantom{\rule{thinmathspace}{0ex}}\sigma \phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}\varphi \phantom{\rule{thinmathspace}{0ex}}\text{cos}\phantom{\rule{thinmathspace}{0ex}}\eta & \text{cos}\phantom{\rule{thinmathspace}{0ex}}\sigma \phantom{\rule{thinmathspace}{0ex}}\text{cos}\phantom{\rule{thinmathspace}{0ex}}\varphi \end{array}\right]\phantom{\rule{thinmathspace}{0ex}}.$$

(3)

Thus, the seven geometric parameters *R _{FD}*,

Many methods have been proposed for CT, PET, or SPECT system calibration, for both fan beam and cone beam geometries.^{1}^{–}^{8} Many of these methods used the coordinate information acquired from projections of “point-like” objects (small ball-bearings with higher or lower density compared to the background) to calculate the system geometric parameters.^{4}^{,}^{5}^{,}^{7}^{,}^{8} For early investigations, simultaneous nonlinear equations of these system parameters were constructed from the projection locations of point objects. To solve these equations, iterative methods, such as Levenberg-Marquard algorithm, were often applied.^{2} These methods required either precise knowledge of the point locations, using a special manufactured calibration phantom,^{8} or good initial values for the system parameters, which in some situations were difficult to assess. Problems involving parameter count restrictions, nonlinear optimization issues, algorithm convergence, and uniqueness of the solution, limited the effectiveness of these methods. To avoid those difficulties in solving the nonlinear equations, Noo *et al.*^{4} proposed a method by introducing intermediate parameters from fitting an ellipse to the projection-orbit data. This method only requires a small set of measurements of a simple phantom consisting of two BBs. This method assumes that one out-of-plane rotation angle σ=0. Also two BBs used in the calibration phantom need to be on the opposite side of the central plane, as an extra requirement for phantom placement. Recently, von Smekal *et al.*^{7} proposed an analytical method to solve six system parameters, except for the source to isocenter distance *R _{FI}*, based on Fourier analysis of the projection-orbit data. This method can explicitly solve these six parameters by canceling out

Based upon the results of previously reported methods and our own investigation, the following observations can be made:

- These two out-of-plane rotation angles, ϕ and σ, are quite difficult to determine with reasonable accuracy.
- These two angles have only a small influence on the image quality compared with other parameters.
- In practical implementation, these two angles can be kept small (less than 1°) by good mechanical design (CAD), high accuracy machining, and careful alignment of phantom.

In this paper, we present a simple method to estimate five system parameters with angle information from the rotation system. We assume that the two out-of-plane rotation angles, ϕ and σ, are zero, based on the observation of von Smekal *et al.*,^{7} the analytic validation in Sec. III A, and computer simulation results in Sec. IV. In addition, only linear regression techniques are used for data fitting, making the proposed method less sensitive to noise and relatively simple to implement.

In this calibration method, we assume that ϕ=0 and σ=0. To validate this assumption, the introduced errors in the calculation of the projected coordinates of the point-like object in the detector plane were calculated. As shown in Fig. 3, for an ideal point object with coordinates of (*x*_{obj},*y*_{obj},*z*_{obj}), if the detector is ideally aligned (i.e., ϕ=0 and σ=0), the point will be projected onto the detector plane at (*u*_{ideal},υ_{ideal}) and

$${u}_{\text{ideal}}={y}_{\text{obj}}\xb7\frac{{R}_{\mathit{\text{FD}}}}{{R}_{\mathit{\text{FI}}}+{x}_{\text{obj}}},\text{\upsilon ideal=zobj\xb7RFDRFI+xobj.}$$

(4)

We first consider out-of-plane rotation in one direction only, e.g., ϕ≠0 and σ=0. We can also assume −90° < ϕ < 90°. Then the position of this object on the detector plane with out-of-plane rotation is (*u*_{wr},υ_{wr}), and

$$\begin{array}{c}{u}_{\text{wr}}={y}_{\text{obj}}\xb7\frac{{R}_{\mathit{\text{FD}}}+{u}_{\text{wr}}\xb7\text{sin}\phantom{\rule{thinmathspace}{0ex}}\varphi}{{R}_{\mathit{\text{FI}}}+{x}_{\text{obj}}}\xb7\frac{1}{\text{cos}\phantom{\rule{thinmathspace}{0ex}}\varphi},\hfill \\ {\upsilon}_{\text{wr}}={z}_{\text{obj}}\xb7\frac{{R}_{\mathit{\text{FD}}}+{u}_{\text{wr}}\xb7\text{sin}\phantom{\rule{thinmathspace}{0ex}}\varphi}{{R}_{\mathit{\text{FI}}}+{x}_{\text{obj}}}.\hfill \end{array}$$

(5)

The relationship between the coordinates with and without rotation can be expressed as

$${u}_{\text{wr}}=\frac{{u}_{\text{ideal}}}{\text{cos}\phantom{\rule{thinmathspace}{0ex}}\varphi -\frac{{u}_{\text{ideal}}}{{R}_{\mathit{\text{FD}}}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}\varphi},\text{\upsilon wr=uideal1\u2212uidealRFDtan\varphi .}$$

(6)

If we define 2α as the fan angle, then

$$\frac{{u}_{\text{ideal}}}{{R}_{\mathit{\text{FD}}}}=\text{tan}\phantom{\rule{thinmathspace}{0ex}}\alpha .$$

(7)

The percentage errors in υ and *u* (for *u*_{ideal}≠0,υ_{ideal}≠0) will be a function of α and ϕ as

$$\begin{array}{c}{\text{error}}_{u}=\frac{{u}_{\text{wr}}-{u}_{\text{ideal}}}{{u}_{\text{ideal}}}=\frac{1}{\text{cos}\phantom{\rule{thinmathspace}{0ex}}\varphi -\text{tan}\phantom{\rule{thinmathspace}{0ex}}\alpha \phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}\varphi}-1,\hfill \\ {\text{error}}_{\upsilon}=\frac{{\upsilon}_{\text{wr}}-{\upsilon}_{\text{ideal}}}{{\upsilon}_{\text{ideal}}}=\frac{\text{tan}\phantom{\rule{thinmathspace}{0ex}}\alpha \phantom{\rule{thinmathspace}{0ex}}\text{tan}\phantom{\rule{thinmathspace}{0ex}}\varphi}{\text{tan}\phantom{\rule{thinmathspace}{0ex}}\alpha \phantom{\rule{thinmathspace}{0ex}}\text{tan}\phantom{\rule{thinmathspace}{0ex}}\varphi -1}.\hfill \end{array}$$

(8)

If we consider the effect of the other out-of-plane rotation angle σ, the results will be identical to those discussed earlier due to symmetry, which means the effect of ϕ on *u* will be exactly the same as the effect of σ on υ. Therefore, Fig. 4 only shows the results of ϕ.

Percentage error in (a) *u*, (b) υ, as a function of out-of-plane rotation angle ϕ, and fan angle α. Fan angle α=5°, 10°, 15°, 20°, to avoid redundant curves, only positive fan angle curves **...**

As shown in Fig. 4, errors increase with increasing fan angle and rotation angle. For α 15°, −2° ϕ 2°, the maximal errors in both υ and *u* are less than 1%, and this 1% error only occurs to the pixels at the edge of the detector with a 15° fan angle. Pixels corresponding to smaller fan angles have errors smaller than 1%. As shown in Fig. 5, the maximal percentage error in both *u* and υ drops to 0.68% for α=10° and 0.37% for α=5°, with ϕ=2°.

From the above discussion, we can assume that these two out-of-plane rotation angles ϕ=0 and σ=0, because they only have minor (normally less than 1% for the worst case on the edge of the detector) effects on system calibration if the detector does not have severe out-of-plane rotations (<2°). Normally with the advantage of computer aided design (CAD), less than 1° out-of-plane rotations can be achieved, which means the maximal error in pixel position will be even smaller than 0.5% on the edge. For example, careful measurement of the detector rotation angle, ϕ, on the breast CT system in our laboratory demonstrated that ϕ < 0.5°. Therefore, with the assumptions that the detector can be physically mounted to yield negligible errors, coupled with the observation that small alignment errors have minimal effect on calibration performance, the system calibration is simplified to a problem of estimating the remaining five system parameters, *R _{FD}*,

For a cone beam CT system, if we assume that only the object rotates and the x-ray tube and detector remain stationary, the orbit of a point in the object (such as a BB) during the scan is a circle in a plane parallel to *x*-*y* plane. The projection of this circle on the detector plane will be an ellipse. Individual points on this ellipse correspond to the BB’s angular position on the circle, relative to an initial reference position. We refer to two points on the circular orbit that are exactly 180° out of phase as a “radial pair,” as illustrated in Fig. 6(a). We can calculate ρ, the distance between a radial pair of points on the detector plane as described in Appendix A.

In Appendix A, we prove that ρ will have its maximum and minimum values when the point object is on the *x* axis or *y* axis (respectively) for cone beam CT systems with a fan angle less than 60° and a cone angle less than 30°.

As illustrated in Fig. 6(b), after the maximum and minimum distances ρ between a radial pair are determined, four corresponding benchmark points can be defined on the detector plane as *A _{ij}*(

For every individual BB indexed by *i*, we define

$${Y}_{i}=\frac{{\upsilon}_{i1}-{\upsilon}_{i2}}{{A}_{i3}{A}_{i4}},\text{Xi=\upsilon i1+\upsilon i22.}$$

(9)

As proven in Appendix B, a linear function *X*=*a*_{1} + *b*_{1}*Y* can fit all the (*Y _{i}*,

$$\begin{array}{c}{\upsilon}_{0}={a}_{1},\hfill \\ {R}_{\mathit{\text{FD}}}={b}_{1}.\hfill \end{array}$$

(10)

To locate the projection line of the rotation center, we first need to locate the projection center of each ellipse *A*_{i0}(*u*_{i0},υ_{i0}), which can be determined by four benchmark points and every other four points (positions) which are 90° apart from each other, as shown in Fig. 7. For any group of four points, we have

$$\begin{array}{c}{u}_{i0}=\frac{\left|\begin{array}{c}\begin{array}{cc}\left|\begin{array}{cc}{u}_{i1}\hfill & {\upsilon}_{i1}\hfill \\ {u}_{i2}\hfill & {\upsilon}_{i2}\hfill \end{array}\right|& {u}_{i1}-{u}_{i2}\end{array}\hfill \\ \begin{array}{cc}\left|\begin{array}{cc}{u}_{i3}\hfill & {\upsilon}_{i3}\hfill \\ {u}_{i4}\hfill & {\upsilon}_{i4}\hfill \end{array}\right|& {u}_{i3}-{u}_{i4}\end{array}\hfill \end{array}\right|}{\left|\begin{array}{cc}\begin{array}{c}{u}_{i1}-{u}_{i2}\\ {u}_{i3}-{u}_{i4}\end{array}\hfill & \begin{array}{c}{\upsilon}_{i1}-{\upsilon}_{i2}\\ {\upsilon}_{i3}-{\upsilon}_{i4}\end{array}\hfill \end{array}\right|},\hfill \\ {\upsilon}_{i0}=\frac{\left|\begin{array}{c}\begin{array}{cc}\left|\begin{array}{cc}{u}_{i1}\hfill & {\upsilon}_{i1}\hfill \\ {u}_{i2}\hfill & {\upsilon}_{i2}\hfill \end{array}\right|& {\upsilon}_{i1}-{\upsilon}_{i2}\end{array}\hfill \\ \begin{array}{cc}\left|\begin{array}{cc}{u}_{i3}\hfill & {\upsilon}_{i3}\hfill \\ {u}_{i4}\hfill & {\upsilon}_{i4}\hfill \end{array}\right|& {\upsilon}_{i3}-{\upsilon}_{i4}\end{array}\hfill \end{array}\right|}{\left|\begin{array}{cc}\begin{array}{c}{u}_{i1}-{u}_{i2}\\ {u}_{i3}-{u}_{i4}\end{array}\hfill & \begin{array}{c}{\upsilon}_{i1}-{\upsilon}_{i2}\\ {\upsilon}_{i3}-{\upsilon}_{i4}\end{array}\hfill \end{array}\right|}.\end{array}$$

(11)

We can go through all the different four-point groups and average the results to get a more accurate estimation of the projection center of an individual ellipse. All these projection centers of the different BB trajectories will form the projection line of the rotation center, and we can use a linear function *U*=*a*_{2} + *b*_{2}*V* to fit all the coordinates (*u*_{i0},υ_{i0}). Since we already have the value of υ_{0}:

$$\begin{array}{c}{u}_{0}={a}_{2}+{b}_{2}{\upsilon}_{0},\hfill \\ \eta ={\text{tan}}^{-1}{b}_{2}.\hfill \end{array}$$

(12)

To calculate the last parameter *R _{FI}*, we need the distance between two BBs,

If we label these two BBs as *i*=1 and *i*=2, as shown in Fig. 8, *h* is defined as the vertical distance between two orbits, *r*_{1} and *r*_{2} are radius for each orbit, and we have

$${l}^{2}={h}^{2}+{r}_{1}^{2}+{r}_{2}^{2}-2{r}_{1}{r}_{2}\phantom{\rule{thinmathspace}{0ex}}\text{cos}({\alpha}_{1}-{\alpha}_{2}).$$

(13)

From previous steps, we already have

$$h=\frac{{A}_{10}{A}_{20}}{{R}_{\mathit{\text{FD}}}}\xb7{R}_{\mathit{\text{FI}}},\text{}{r}_{1}=\frac{{A}_{13}{A}_{14}}{2\xb7{R}_{\mathit{\text{FD}}}}\xb7{R}_{\mathit{\text{FI}}},\text{}{r}_{2}=\frac{{A}_{23}{A}_{24}}{2\xb7{R}_{\mathit{\text{FD}}}}\xb7{R}_{\mathit{\text{FI}}}.$$

(14)

If we define the original angular positions (before the system starts rotating) of BBs 1 and 2 as θ_{10} and θ_{20}, respectively, then

$$\begin{array}{c}{\alpha}_{1}={\theta}_{10}-{\theta}_{14},\hfill \\ {\alpha}_{2}={\theta}_{20}-{\theta}_{24}.\hfill \end{array}$$

(15)

We thus solve *R _{FI}* using

$${R}_{\mathit{\text{FI}}}=\frac{l\xb7{R}_{\mathit{\text{FD}}}}{\sqrt{({A}_{10}{A}_{20}{)}^{2}+(\frac{{A}_{13}{A}_{14}}{2}{)}^{2}+(\frac{{A}_{23}{A}_{24}}{2}{)}^{2}-\frac{{A}_{13}{A}_{14}\xb7{A}_{23}{A}_{24}\xb7\text{cos}\left({\alpha}_{1}-{\alpha}_{2}\right)}{2}}}.$$

(16)

To validate the accuracy of the calibration method described in Sec. III, the projection coordinates of BBs on the detector plane were simulated. The calibration method was then applied on the simulated data and the calculated system parameters were compared with their known values. Gaussian noise of different standard deviation was added to the projection coordinates to simulate the effect of noise in the practical measurement environment and to evaluate the consistency of this method in the presence of noise.

In this simulation, the detector matrix size was set to 2048 × 1024 with a pixel size of 48 µm × 48 µm, similar to the mouse CT scanner assembled in our laboratory. The five system parameters were set as: *R _{FD}*=400.0000 mm,

Eight BBs were simulated and the distance between the first and the last BB was *l*=14.39 mm. Five hundred points were simulated for each projection orbit, corresponding to a 500-view CT scan in 360°. One example of the simulation is shown in Fig. 9.

Simulated projection orbits with Gaussian noise. R_{FD}=400 mm, R_{FI}=150 mm, u_{0}=1005 pixel, v_{0}=480 pixel, η=−1.0°, standard deviation of Gaussian noise=0.4 pixel.

Simulations were performed in three different groups:

- No preset out-of-plane rotation, with different noise levels, results are shown in Table I and Table II;
- With different out-of-plane rotation angles, without noise, results are shown in Table III;
- With specific preset out-of-plane rotations (ϕ=1.5000°,σ=1.2000°), with different noise levels, results are shown in Table IV and Table V.

Linear regression was used in two steps of the calibration method, and the *R*-square value or the coefficient of determination, which is the indicator of how well the fitting works, was shown in Table I and Table IV. Gaussian noise with standard deviations *s* of 1%, 20%, 40%, and 100% of the pixel size were independently added into the horizontal and vertical coordinates of every projection position. Most *R*-square values calculated were close to 1.0, corresponding to an ideal linear fitting. Especially in calculation of *R _{FD}* and υ

The system parameters calibrated from computer simulation are listed in Table II, Table III, and Table V. Though the results degraded slightly with increasing noise, the calibrated values were nevertheless accurate in comparison with their true values.

From the results shown in Table III–Table V, when the out-of-plane rotation angles are relatively small (less than 2°), there will be minor effects caused by neglecting these angles compared with the effects caused by the noise in the measurement of BB positions. And normal engineering design and machinery can satisfy this loose requirement by limiting out-of-plane rotations to within 1°. These results demonstrate that good calibration accuracy can be achieved assuming ϕ=0 and σ=0.

Computer simulations were performed with various system parameters including the number of views, number of BBs used, and different combinations of system parameters. All the calibration results under various conditions validated the accuracy of this calibration method.

To compare with previous results, simulations with same rotation angles and noise level as in Smekal’s^{7} work were performed and calibration results are given in Table VI. In Table VI, relative errors were calculated by the ratio between uncertainties and mean values. The method presented here has equal or better accuracy as the results reported by Smekal *et al.*^{7}

For practical implementation of the calibration method, small ball bearings (with diameters of 0.25 or 2.3 mm in our study for two different scanners) were used as point-like objects to provide projection coordinates. Because of its finite size and the cone beam system magnification, the projection of an individual BB covers more than one pixel in the detector plane, as shown in Fig. 10. The projection images were first processed to reverse the gray scale and thresholded to reduce the background noise. Then the projection coordinates were calculated as the mass center of the x-ray shadow area by

$$\begin{array}{c}{u}_{m}=\frac{{\mathrm{\Sigma}}_{u}{\mathrm{\Sigma}}_{\upsilon}\text{Image}(u,\upsilon )\xb7u}{{\mathrm{\Sigma}}_{u}{\mathrm{\Sigma}}_{\upsilon}\text{Image}(u,\upsilon )},\hfill \\ {\upsilon}_{m}=\frac{{\mathrm{\Sigma}}_{u}{\mathrm{\Sigma}}_{\upsilon}\text{Image}(u,\upsilon )\xb7\upsilon}{{\mathrm{\Sigma}}_{u}{\mathrm{\Sigma}}_{\upsilon}\text{Image}(u,\upsilon )},\hfill \end{array}$$

(17)

while Image(*u*,υ) is the gray scale value of detector pixel (*u*,υ) after processing. The mass center was calculated iteratively three times for each view (Please refer to Appendix C). Computer simulation was performed to demonstrate that the calculated mass center has the error of less than 0.05 pixel from the true projection center of the BB. A computer program was developed to automatically track BB positions along the trajectory continuously through all the views. Detailed simulation and tracking procedures are also provided in Appendix C.

The calibration method was applied to a micro CT system developed in our lab. In this system, an x-ray tube (XTF5011, Oxford Instruments, Scotts Valley, CA) with a 70 µm focal spot was used as the x-ray source, and was capable of 4–50 kVp and 0–1.0 mA operation. The CMOS detector (Shad-o-Box 2048, Rad-icon Imaging Corp., Santa Clara, CA) used in this system has an active area of 50 × 100 mm^{2} and pixel size of 48 µm, and produces 1024 × 2048 projection images. A rotary stage driven by a stepping motor (MDrive 23, Intelligent Motion Systems, Inc., Marlborough, CT) was used to rotate the object. The whole system is shown in Fig. 11(a). Eight stainless steel ball bearings with a diameter of 254 µm were taped on the surface of a 19-mm-diameter test tube to make the BB phantom. The distance between BB 1 and BB 8 was measured as: *l*=17.2 mm (Fig. 11(b)). Every scan was performed with 30.0 kVp, 0.5 mA, and 500 views. The system parameters were measured 10 times under exactly the same conditions. The measured mean and standard deviation values were calculated and shown in Table VII.

The measured results demonstrated excellent consistency of the proposed calibration method. Example CT images of a mouse reconstructed with these system parameters are shown in Figs. 12(a)–12(d). These images depict excellent detail and sharp edges of the bone structures, qualitatively illustrating the performance of the calibration method.

In contrast to the micro CT system, in which the object rotates, a breast CT system constructed in our lab was designed such that the x-ray tube and detector rotate around the object. This system uses a 40 × 30 cm amorphous silicon detector (PaxScan 4030CB, Varian Medical Systems, Salt Lake City, UT) with 194 µm pixel size, corresponding to a 2048 × 1536 image size. In this study, the detector was working in the mode with 2 × 4 pixel binning, corresponding to a 1024 × 384 image size. A Pantak HF160 x-ray generator with a 0.4 mm focal spot x-ray tube (Comet, MXR-160/20), capable of 160 kVp and 6 mA (continuous) operation was used. The system design is shown in Fig. 13(a). A set of lead ball bearings with a diameter of 2.3 mm was mounted into a foam board to make the BB phantom. The distance between BB 1 and BB 6 was measured as: *l*=100.00 mm (Fig. 13(b)). Every scan was performed with 80.0 kVp, 7.0 mA, and 500 views. Ten calibration scans were acquired in a period of about 2 months with the detector normally positioned (<1° in-plane rotation), and the calibration results are shown in Table VIII. The results demonstrated excellent reproducibility of the proposed method over a relative long time period.

To test the performance of this method under a large inplane rotation angle, the detector was deliberately tilted with an in-plane rotation angle of ~7°. In this configuration, only one calibration was obtained. The system parameters were measured as: *u*_{0}=383.5924 pixel, υ_{0}=49.7066 pixel, η=7.5995°, *R _{FI}*=465.1768 mm,

The R-square values for the linear fitting were 0.999 920 4 and 0.999 992 9, respectively.

CT images of a bone phantom reconstructed with these system parameters are shown in Fig. 12(e). As compared with Fig. 12(f), excellent details and sharp edges of the bone structures were demonstrated with correct system parameters. The calibration method works consistently even with a relative large detector in-plane rotation angle.

In this work, a new and simple method for cone beam CT calibration was presented. The previous seven-parameter problem was simplified to a five-parameter problem. The maximal error introduced by this assumption was found to be less than 1% for typical cone beam CT systems. High resolution angle information from the rotation system, basic information for CT systems, is leveraged to achieve higher precision than with previous methods. Linear regression performance validated the accuracy and consistency of this method in the presence of noise. Only a rough measurement of the distance between two BBs was required to calculate the source to isocenter distance *R _{FI}*. The error in

From the results of computer simulation and practical system implementation in a micro-CT scanner and a breast CT scanner, equal or higher accuracy was achieved with a simpler implementation. The proposed method can be applied to any cone beam CT system with a mechanically stable rotation center and without severe detector out-of-plane rotations (<2°).

This work was supported in part by National Institutes of Health Grant Nos. R21 EB04643 and R01 EB02138, and by grant from the California Breast Cancer Research Program (CA BCRP 11IB-0114).

As shown in Fig. 14, the distance between a radial pair of points on the detector plane, ρ, will be a function of θ(0° θ 90°). This distance does not depend on the in-plane rotation of the detector, so we can assume η=0. For a circular trajectory located at plane *z*=*z _{i}*, with radius

$$\begin{array}{c}{u}_{1}={u}_{0}+\frac{{R}_{\mathit{\text{FD}}}\xb7r\phantom{\rule{thinmathspace}{0ex}}\text{cos}\phantom{\rule{thinmathspace}{0ex}}\theta}{{R}_{\mathit{\text{FI}}}+r\phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}\theta},\text{}{u}_{2}={u}_{0}-\frac{{R}_{\mathit{\text{FD}}}\xb7r\phantom{\rule{thinmathspace}{0ex}}\text{cos}\phantom{\rule{thinmathspace}{0ex}}\theta}{{R}_{\mathit{\text{FI}}}-r\phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}\theta},\hfill \\ {\upsilon}_{1}={\upsilon}_{0}+\frac{{R}_{\mathit{\text{FD}}}\xb7{z}_{i}}{{R}_{\mathit{\text{FI}}}+r\phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}\theta},\text{}{\upsilon}_{2}={\upsilon}_{0}+\frac{{R}_{\mathit{\text{FD}}}\xb7{z}_{i}}{{R}_{\mathit{\text{FI}}}-r\phantom{\rule{thinmathspace}{0ex}}\text{sin}\phantom{\rule{thinmathspace}{0ex}}\theta},\hfill \end{array}$$

(A1)

$$\begin{array}{cc}{\rho}^{2}\hfill & =({u}_{1}-{u}_{2}{)}^{2}+({\upsilon}_{1}-{\upsilon}_{2}{)}^{2}\hfill \\ & =4\xb7{R}_{\mathit{\text{FD}}}^{2}\xb7\left(\frac{{R}_{\mathit{\text{FI}}}^{2}-{z}_{i}^{2}}{{R}_{\mathit{\text{FI}}}^{2}-{r}^{2}\phantom{\rule{thinmathspace}{0ex}}{\text{sin}}^{2}\theta}+\frac{{R}_{\mathit{\text{FI}}}^{2}({z}_{i}^{2}+{r}^{2}-{R}_{\mathit{\text{FI}}}^{2})}{{({R}_{\mathit{\text{FI}}}^{2}-{r}^{2}\phantom{\rule{thinmathspace}{0ex}}{\text{sin}}^{2}\theta )}^{2}}\right)\phantom{\rule{thinmathspace}{0ex}}.\hfill \end{array}$$

(A2)

Let

$$x=\frac{2\xb7{R}_{\mathit{\text{FD}}}}{{R}_{\mathit{\text{FI}}}^{2}-{r}^{2}\phantom{\rule{thinmathspace}{0ex}}{\text{sin}}^{2}\phantom{\rule{thinmathspace}{0ex}}\theta},$$

then ρ^{2} will be a quadratic function of *x*,

$${\rho}^{2}={R}_{\mathit{\text{FI}}}^{2}({z}_{i}^{2}+{r}^{2}-{R}_{\mathit{\text{FI}}}^{2}){x}^{2}+2\xb7{R}_{\mathit{\text{FD}}}\xb7({R}_{\mathit{\text{FI}}}^{2}-{z}_{i}^{2})x.$$

(A3)

Because 0° θ 90°, then

$$\frac{2\xb7{R}_{\mathit{\text{FD}}}}{{R}_{\mathit{\text{FI}}}^{2}}x\frac{2\xb7{R}_{\mathit{\text{FD}}}}{{R}_{\mathit{\text{FI}}}^{2}-{r}^{2}},\text{and this implies}\phantom{\rule{thinmathspace}{0ex}}{R}_{\mathit{\text{FI}}}r.$$

*Case 1*: If
${z}_{i}^{2}+{r}^{2}-{R}_{\mathit{\text{FI}}}^{2}=0$.

${\rho}^{2}=2\xb7{R}_{\mathit{\text{FD}}}\xb7({R}_{\mathit{\text{FI}}}^{2}-{z}_{i}^{2})x$, ρ^{2} will be a linear function of *x* and have its maximum value when θ=90° and minimum value when θ=0°.

*Case 2*: If
${z}_{i}^{2}+{r}^{2}-{R}_{\mathit{\text{FI}}}^{2}<0$.

${\rho}^{2}={R}_{\mathit{\text{FI}}}^{2}({z}_{i}^{2}+{r}^{2}-{R}_{\mathit{\text{FI}}}^{2}){x}^{2}+2\xb7{R}_{\mathit{\text{FD}}}\xb7({R}_{\mathit{\text{FI}}}^{2}-{z}_{i}^{2})x$ will be a quadratic equation and have maximum value at,

$${x}_{\text{max}}=\frac{{R}_{\mathit{\text{FD}}}({R}_{\mathit{\text{FI}}}^{2}-{z}_{i}^{2})}{{R}_{\mathit{\text{FI}}}^{2}({R}_{\mathit{\text{FI}}}^{2}-{z}_{i}^{2}-{r}^{2})}.$$

(A4)

We can easily prove that if ${z}_{i}^{2}+2{r}^{2}<{R}_{\mathit{\text{FI}}}^{2}$, then

$${x}_{\text{max}}<\frac{2{R}_{\mathit{\text{FD}}}}{{R}_{\mathit{\text{FI}}}^{2}}<\frac{2{R}_{\mathit{\text{FD}}}}{{R}_{\mathit{\text{FI}}}^{2}-{r}^{2}},$$

as shown in Fig. 15(a).

So ρ^{2} will have its maximum value when θ=0° and minimum value when θ=90°.

*Case 3*: If
${z}_{i}^{2}+{r}^{2}-{R}_{\mathit{\text{FI}}}^{2}>0$.

${\rho}^{2}={R}_{\mathit{\text{FI}}}^{2}({z}_{i}^{2}+{r}^{2}-{R}_{\mathit{\text{FI}}}^{2}){x}^{2}+2\xb7{R}_{\mathit{\text{FD}}}\xb7({R}_{\mathit{\text{FI}}}^{2}-{z}_{i}^{2})x$ wil1 be another quadratic equation. From Figs. 15(b) and 15(c), we can see that, in either situation, ρ^{2} will have its maximum value when θ=90° and minimum value when θ=0°.

In summary, we only need
${z}_{i}^{2}+2{r}^{2}<{R}_{\mathit{\text{FI}}}^{2}$ or (*z _{i}*/

$$(\frac{r}{{R}_{\mathit{\text{FI}}}}{)}_{\text{max}}=\text{tan}\phantom{\rule{thinmathspace}{0ex}}\alpha ,\text{(ziRFI)max=tan\beta ,}$$

(A5)

If a system has a cone angle less than 30° and a fan angle less than 60°. (tan β)^{2}+2(tan α)^{2}<1. This limit is relatively loose for a practical system. For example, our micro-CT system has a *R _{FI}*=200 mm, with a detector active area of 10 × 5 cm

As shown in Fig. 16,

$${A}_{i1}{A}_{0}=\frac{{z}_{i}}{{R}_{\mathit{\text{FI}}}-{r}_{i}}\xb7{R}_{\mathit{\text{FD}}},\text{}{A}_{i2}{A}_{0}=\frac{{z}_{i}}{{R}_{\mathit{\text{FI}}}+{r}_{i}}\xb7{R}_{\mathit{\text{FD}}}.$$

(B1)

Then

$$\begin{array}{c}\frac{{A}_{i1}{A}_{0}+{A}_{i2}{A}_{0}}{2}=\frac{{z}_{i}\xb7{R}_{\mathit{\text{FD}}}\xb7{R}_{\mathit{\text{FI}}}}{{R}_{\mathit{\text{FI}}}^{2}-{r}_{i}^{2}},\hfill \\ \frac{{A}_{i1}{A}_{0}-{A}_{i2}{A}_{0}}{2}=\frac{{z}_{i}\xb7{R}_{\mathit{\text{FD}}}\xb7{r}_{i}}{{R}_{\mathit{\text{FI}}}^{2}-{r}_{i}^{2}}.\hfill \end{array}$$

(B2)

So

$$\frac{{A}_{i1}{A}_{0}-{A}_{i2}{A}_{0}}{2}=\frac{{r}_{i}}{{R}_{\mathit{\text{FI}}}}\xb7\frac{{A}_{i1}{A}_{0}+{A}_{i2}{A}_{0}}{2}.$$

(B3)

And we also have

$$\begin{array}{c}\frac{{r}_{i}}{{R}_{\mathit{\text{FI}}}}=\frac{{A}_{i3}{A}_{i4}}{2\xb7{R}_{\mathit{\text{FD}}}},\hfill \\ {A}_{i1}{A}_{0}=\frac{({\upsilon}_{i1}-{\upsilon}_{0})}{\text{cos}\phantom{\rule{thinmathspace}{0ex}}\eta},\text{}{A}_{i2}{A}_{0}=\frac{({\upsilon}_{i2}-{\upsilon}_{0})}{\text{cos}\phantom{\rule{thinmathspace}{0ex}}\eta}.\hfill \end{array}$$

(B4)

Insert these into Eq. (B3), and we have

$$\frac{({\upsilon}_{i1}+{\upsilon}_{i2})}{2}={\upsilon}_{0}+{R}_{\mathit{\text{FD}}}\frac{{\upsilon}_{i1}-{\upsilon}_{i2}}{{A}_{i3}{A}_{i4}}.$$

(B5)

Define:

$${Y}_{i}=\frac{{\upsilon}_{i1}-{\upsilon}_{i2}}{{A}_{i3}{A}_{i4}},\text{}{X}_{i}=\frac{{\upsilon}_{i1}+{\upsilon}_{i2}}{2}.$$

(B6)

Then use a linear function *X*=*a*_{1} + *b*_{1}*Y* to fit all the (*Y _{i}*,

$$\begin{array}{c}{\upsilon}_{0}={a}_{1},\hfill \\ {R}_{\mathit{\text{FD}}}={b}_{1}.\hfill \end{array}$$

(B7)

Appendix B. Calculation of υ_{0} and *R*_{FD}.

The cone beam projection image of a spherical object was simulated by analytical calculation of the attenuation of an x-ray beam. In this simulation, monoenergetic photons were emitted from an isotropic point source. Output photons were uniformly distributed in 4π space. Gaussian noise was added to the incident photons and an example of the simulated projection image is shown in Fig. 17. This image was then processed to reverse the gray scale and thresholded to reduce the interference of background noise. A initial position
$({u}_{0}^{0},{\upsilon}_{0}^{0})$ and a window size *k* were set manually from observation. Thus the mass center for individual view with index *n* was calculated within the window
$({u}_{n}^{0}\pm k,{\upsilon}_{n}^{0}\pm k)$. This calculation was repeated iteratively as

$$\begin{array}{c}{u}_{n}^{i}=\frac{{\displaystyle {\sum}_{u={u}_{n}^{i-1}-k}^{u={u}_{n}^{i-1}+k}{\displaystyle {\sum}_{\upsilon ={\upsilon}_{n}^{i-1}-k}^{\upsilon ={\upsilon}_{n}^{i-1}+k}I(u,\upsilon )\xb7u}}}{{\displaystyle {\sum}_{u={u}_{n}^{i-1}-k}^{u={u}_{n}^{i-1}+k}{\displaystyle {\sum}_{\upsilon ={\upsilon}_{n}^{i-1}-k}^{\upsilon ={\upsilon}_{n}^{i-1}+k}I(u,\upsilon )}}},\hfill \\ {\upsilon}_{n}^{i}=\frac{{\displaystyle {\sum}_{u={u}_{n}^{i-1}-k}^{u={u}_{n}^{i-1}+k}{\displaystyle {\sum}_{\upsilon ={\upsilon}_{n}^{i-1}-k}^{\upsilon ={\upsilon}_{n}^{i-1}+k}I(u,\upsilon )\xb7\upsilon}}}{{\displaystyle {\sum}_{u={u}_{n}^{i-1}-k}^{u={u}_{n}^{i-1}+k}{\displaystyle {\sum}_{\upsilon ={\upsilon}_{n}^{i-1}-k}^{\upsilon ={\upsilon}_{n}^{i-1}+k}I(u,\upsilon )}}},\text{i=1,2,\u2026,Ni.}\hfill \end{array}$$

(C1)

*I*(*u*,υ) is gray scale value of the projection image at (*u*,υ) after processing. Ni is the total number of iterations.

Once the projection center position was determined, it was used as the initial position for the next view and the same procedure was repeated for the next projection image. The complete trajectory was determined through all the views over 360°. The trajectory tracking diagram is shown in Fig. 18.

Diagram of trajectory tracking. *I*_{n}(*u*,υ) is the gray scale value of position (*u*,υ) in view # *n*, after processing. *N* is the total number of views, *k* is the window size set from observation so that a 2k*2k window completely covers one BB **...**

Computer simulation was performed to validate the accuracy of the mass center calculation and determine the optimal iteration numbers. For a properly chosen window size, the results were always convergent in less than three iterations, even when the initial position was set at the edge of the x-ray shadow area.

In the computer simulation, 1 × 10^{6} photons with a noise level of 1 × 10^{4} were incident into each detector pixel. The x-ray shadow covered an area of 25 × 25 pixel^{2}. True values of the projection center coordinates were calculated as: *u*_{0}=415.7407 pixel, υ_{0}=242.1296 pixel.

With three iterations for each simulation, a 30 × 30 window, initial position at *u*_{0}=402 pixel, υ_{0}=230 pixel, the results of 10 runs were: *u*_{0}=415.75±0.05 pixel, υ_{0}=242.12±0.03 pixel.

Simulations using different shadow areas (BB size), noise levels, and initial positions were studied and the calculated results all had an error less than 0.05 pixel from the true values.

Kai Yang, Department of Radiology, University of California, Davis Medical Center, 4701 X Street, Sacramento, California 95817 and Department of Biomedical Engineering, University of California, Davis, California 95616.

Alexander L. C. Kwan, Department of Radiology, University of California, Davis Medical Center, 4701 X Street, Sacramento, California 95817.

DeWitt F. Miller, Department of Radiology, University of California, Davis Medical Center, 4701 X Street, Sacramento, California 95817 and Department of Biomedical Engineering, University of California, Davis, California 95616.

John M. Boone, Department of Radiology, University of California, Davis Medical Center, 4701 X Street, Sacramento, California 95817 and Department of Biomedical Engineering, University of California, Davis, California 95616.

1. Gullberg GT, Tsui BM, Crawford CR, Edgerton ER. Estimation of geometrical parameters for fan beam tomography. Phys. Med. Biol. 1987;32:1581–1594.

2. Gullberg GT, Tsui BM, Crawford CR, Ballard JG, Hagius JT. Estimation of geometrical parameters and collimator evaluation for cone beam tomography. Med. Phys. 1990;17:264–272. [PubMed]

3. Hsieh J. Three-dimensional artifact induced by projection weighting and misalignment. IEEE Trans. Med. Imaging. 1999;18:364–368. [PubMed]

4. Noo F, Clackdoyle R, Mennessier C, White TA, Roney TJ. Analytic method based on identification of ellipse parameters for scanner calibration in cone-beam tomography. Phys. Med. Biol. 2000;45:3489–3508. [PubMed]

5. Karolczak M, Schaller S, Engelke K, Lutz A, Taubenreuther U, Wiesent K, Kalender W. Implementation of a cone-beam reconstruction algorithm for the single-circle source orbit with embedded misalignment correction using homogeneous coordinates. Med. Phys. 2001;28:2050–2069. [PubMed]

6. Stevens GM, Saunders R, Pelc NJ. Alignment of a volumetric tomography system. Med. Phys. 2001;28:1472–1481. [PubMed]

7. von Smekal L, Kachelriess M, Stepina E, Kalender WA. Geometric misalignment and calibration in cone-beam tomography. Med. Phys. 2004;31:3242–3266. [PubMed]

8. Cho Y, Moseley DJ, Siewerdsen JH, Jaffray DA. Accurate technique for complete geometric calibration of cone-beam computed tomography systems. Med. Phys. 2005;32:968–983. [PubMed]

9. Feldkamp LA, Davis LC, Kress JW. Practical cone-beam algorithm. J. Opt. Soc. Am. A. 1984;1:612–619.

10. Kak AC, Slaney M. Principles of Computerized Tomographic Imaging. Philadelphia: SIAM; 1988.

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