A Simple and Precise Procedure for a Complete Characterization of a Cone-Beam Computed Tomography System

1. Introduction

Existing geometric calibrations rely on the decision of a parameters-set that accurately characterize the imaging system’s geometry, which is setting the axis of rotation (AOR) normal to the sample platform. The well-known calibration phantoms containing steel beads set in some special pattern are used in most geometric calibration methods. The most types of extensively used phantom includes beads in circular [1,2,3], helical [4,5], a line [6,7,8,9,10], or other special arrangements [11,12].

Shih et al. [6] proposed and validated an analytic geometric scheme to calibrate, find and correct any misalignment in the imaging detector relative to the X-ray source. Here, we modify that scheme to include rectifying images taken with a tilt in the sample platform (or holder), thus calibrating the complete cone-beam computed tomography (CBCT) system including its tilt of sample platform. We validate this new scheme by using a helical-beads phantom and a line-beads phantom, allowing a tilt between the AOR and the central axis of the phantom (AOP). Sec. 2.1 will show the dental CBCT geometry and phantoms of calibration. Sec. 2.2 will show how the beads in the two chosen phantoms can have varied rotational radii. The fundamental theorem of an ideal CBCT system, where the beads in the phantoms may have varied rotating radii will be defined in Sec. 2.3. The precise relations between minor misalignment and perfect alignment will be delineated in Sec. 2.4–the relations that provide the necessary variables to depict the misalignment to constitute the foundation of the new precision alignment loop (PAL) generalized from our earlier version [6]. Sec. 3.1 will show the criterion for validating PAL’s analytical simulations. Secs. 3.2 and 3.3 will validate the new PAL by detecting the misaligned parameter-set in the simulations for the helical-beads phantom (Sec. 3.2) and the line-beads phantom (Sec. 3.3) having a slightly tilted AOP with respect to the AOR. In our analytical simulations, we set the beads and the source of X-ray as points, and each bead was projected on a misaligned detector by stretching a line from the source to the bead to intersect the misaligned detector. Sec. 4 will present the calibrating results of a dental CBCT system using our new scheme.

2. Materials and Methods

2.1. The Dental CBCT Geometry

The dental CBCT system under study exposed in Figure 1a has a flat detector and an X-ray source on the opposite sides of a AOR mounted on a rotating gantry. The central ray from the X-ray focal point S, passes O on the AOR perpendicularly, and intersects the detector at the geometric center G, or simply “the center G”. The three points, S, O, and G lie on the system’s principal axis, the z-axis, originating from O pointing to S. The y-axis is the AOR. The x-axis with its origin at O follows the rule of right-hand screw. The distance between S and O is source-to-rotational-axis distance (SRD), and the distance between S and G is source-to-detector distance (SDD). The coordinate system that defines the detector's orientation is illustrated in Figure 1b.

In an ideally aligned system, the U’-V’-W’ coordinate, a local system, is a translation of the x-y-z coordinate, the global system, along z-axis to the detector. The position (U_i, V_j) of the (i, j)^th pixel is given by U-V system, the body coordinate of detector, co-planar with the U’-V’ system. If the system is not aligned, then the misalignment (θ, φ, η) can be characterized by a rotation of the axes at the center G: slant angle η about the W’-axis, skew angle φ about the V’-axis, and tilt angle θ about the U’-axis. In a misaligned system, as θ or φ is not zero, the z-axis is not overlapping with the W’-axis. The set (U_G, V_G, SRD, SDD, θ, φ, η) parameterizes the CBCT system. In an ideal aligned system, (θ, φ, η) = (0, 0, 0).

A misaligned system produces distorted images. Shih et al. [6] showed to find angles (θ, φ, η) accurately with precision can revive the fidelity of image. Our current approach varies from the previous works (cited in Sec. 1) by using PAL’s iterative process that starts with a data set either from a simulation or an actual CT image to acquire the angles (θ, φ, η). The essence of PAL is exposed in Figure 2, but we modify it to include finding the tilt angle ζ between the AOP and the AOR in this study.

We calibrate the CBCT system using a cylindrical phantom of helical-beads (Figure 3a), and a cylindrical phantom of line-beads (Figure 3c).The helical-beads phantom is a cylinder of 40-mm radius having eleven 2-mm diameter steel beads evenly separated by 36° in azimuth and 10 mm in height. The line-beads phantom is a cylinder of 40-mm radius having eleven 2-mm diameter steel beads by 10 mm even spacing on a line parallel to the AOP.

We acquired a projection images-set of a calibration phantom, inverted the images, and then used a method of threshold segmentation to remove each projection’s background around beads (Figure 3b,d), then applied the center-of-mass idea to the images to calculate out the exact beads’ coordinates precisely.

2.2. Varied Rotational Radii of Beads in Phantoms

The beads on the two cylindrical phantoms shown in Figure 4 can have systematically varying radii of rotation under two conditions: 1) When the helical-beads phantom’s AOP is parallel to the AOR, but displaced by a distance c in the x-z plane as exposed in Figure 4a; the beads’ rotating radii of this phantom are within the range, r_p-c to r_p+c, where r_p is the radius of the phantom. 2) When there is a small tilt angle ζ between the line-beads phantom’s AOP and the AOR, the beads’ rotating radii of this phantom vary as exposed in Figure 3b. The tilt angle ζ can be deliberate by Eq. (1):

$$\sin \zeta ={\displaystyle \frac{(r_N-r_1)}{(N-1)d}}$$

(1)

where d is two adjacent beads’ vertical separation, N is the number of beads, and r₁, r₂, r₃, …, and r_N are the respective beads’ rotating radii of the line-beads phantom. Eq. (1) is also applied to the helical-beads phantom, since the first bead and the last bead are located on a line parallel to the AOP. For detecting the tilt angle accurately, it is important to setup the AOR, AOP, and the line of beads in the same plane.

The actual vertical separation d’ must be modified according to Eq. (2):

$$d\text{'}=d\times \cos \zeta$$

(2)

2.3. The Ideal System

As shown in our previous work [6], the projections of the rotating beads on the detector are ellipses of different eccentricities. In an ideal detector, θ = φ = η = 0. Figure 5 is a ray diagram showing points H, I, J, K, L, and M on the ideal detector (in yellow) as projections of the respective positions A, B, C, D, E, and F on a circular track traced out by a bead on the rotating phantom. The six critical projected points are sufficient to derive the equations to yield the necessary set of elliptical parameters: (U_L, V_L), the elliptical center; a and b, the semi-major and semi-minor axes; and η*, the slant angle of the ellipse. In an ideal detector, η* = 0.

As exposed in Figure 5a, two elliptical co-vertices J and M on the detector are the respective projected positions of C and F, therefore we have

$$S_{\mathrm{J}}=h{\displaystyle \frac{SDD}{SRD+r}};S_{\mathrm{M}}=h{\displaystyle \frac{SDD}{SRD-r}}$$

(3)

where S_J and S_M are the distances from M and J (the co-vertices) to the center G (Figure 5a), respectively, r the bead’s rotating radius, and h the distance from the rotating bead to the z-axis.

The system’s magnification factor at D, the center of bead’s circular track, M is related to the position of K, the projection of D on the detector:

$$M={\displaystyle \frac{S_{\mathrm{K}}}{h}}={\displaystyle \frac{SDD}{SRD}};S_{\mathrm{K}}={\displaystyle \frac{SDD}{SRD}}h=M\times h$$

(4)

where S_K is the distance from K to the center G, as exposed in Figure 5a.

The magnification factor M’, at the critical positions of A, B, and E, is related to their respective projections at H, I, and L on the detector, as exposed in Figure 5b:

$$M\text{'}={\displaystyle \frac{SDD}{SRD-\delta }}={\displaystyle \frac{S_{\mathrm{L}}}{h}}={\displaystyle \frac{a}{\sqrt{r^2-{\delta }^2}}}={\displaystyle \frac{SR{D}^2}{SR{D}^2-r^2}}{\displaystyle \frac{SDD}{SRD}}$$

(5)

where δ =r²/SRD, and S_L is the distance from the elliptical center L to the center G.

$$S_{\mathrm{L}}={\displaystyle \frac{(S_{\mathrm{M}}+S_{\mathrm{J}})}{2}}=\left({\displaystyle \frac{SR{D}^2}{SR{D}^2-r^2}}{\displaystyle \frac{SDD}{SRD}}\right)h=M\text{'}\times h$$

(6)

The semi-major axis a and the rotating radius rare related according to Eq. (7):

$$a=r{\displaystyle \frac{SDD}{\sqrt{SR{D}^2-r^2}}};r=a{\displaystyle \frac{SRD}{\sqrt{SD{D}^2+a^2}}}$$

(7)

Since the beads’ rotating radii vary, those projected ellipses on the detector have different semi-major axes.

For the projected ellipse on the detector, the relations between ‘b’, the semi-minor axis, ‘h’, the rotating bead’s height, ‘r’, this bead’s rotating radius, and ‘S_L’, the distance between L, the center of the ellipse and the center G, are described by Eq. (8).

$$b={\displaystyle \frac{(S_{\mathrm{M}}-S{}_{\mathrm{J}})}{2}}=\left({\displaystyle \frac{SR{D}^2}{SR{D}^2-r^2}}{\displaystyle \frac{SDD}{SRD}}{\displaystyle \frac{r}{SRD}}\right)h=S_{\mathrm{L}}{\displaystyle \frac{r}{SRD}};S_{\mathrm{L}}={\displaystyle \frac{b}{r}}SRD$$

(8)

Eq. (8) gives S_L, the beads’ projected elliptical centers as a function of b/r, semi-minor axis divides by respective bead’s rotating radius, and S_L = S_L(b/r). The S_L vs. b/r plot for the elliptical centers above and below the center G, yields V_G, since all the projected ellipses’ centers constitute a vertical line, then we acquire the center G, (U_G, V_G), from which. After finding the center G, S_K can be calculated out by Eq. (4).

If a rotating bead is on the system’s principal axis (z-axis), i.e., h = 0, then on the detector, b =0, and S_L = 0; and projection is a horizontal line of length 2a passing through the center G, a crucial point on the detector. Without a bead located at h = 0, however, we can still find G by considering the following. If we plot all the projected ellipses on the U-V plane, we see the ellipses flatten as the corresponding circular track’s h approaches 0, as expected. This suggests that in a V_L vs. S_L plot, the ellipses will provide points to which two straight lines with opposite slopes can be fitted and they intersect at the V for G. This will be shown when we present the calibration results in Sec. 4.

As the rotating bead is closer to the z-axis, then b→0, on the detector. The elliptical parameters (U_L, V_L, a, b, η*) will be acquired by the linear interpolation rather than fitting an ellipse.

In an ideal detector, the vertical projected distance, S_d, between two adjacent beads’ rotating centers is the same, i.e.,

$$S_d=Md\text{'}$$

(9)

2.4. The Misaligned System

We consider the condition as η, φ, and θ ≤ 10°, then sinη～η, sinφ～φ, and sinθ～θ, cosη～cosφ～cosθ～1, and after neglecting the terms of second-order, then the simplified matrix of transformation T for a misaligned system is as:

$$T=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \mathrm{\theta }& -\sin \mathrm{\theta }\\ 0& \sin \mathrm{\theta }& \cos \mathrm{\theta }\end{array}\right]\left[\begin{array}{ccc}\cos \mathrm{\phi }& 0& \sin \mathrm{\phi }\\ 0& 1& 0\\ -\sin \mathrm{\phi }& 0& \cos \mathrm{\phi }\end{array}\right]\left[\begin{array}{ccc}\cos \mathrm{\eta }& \sin \mathrm{\eta }& 0\\ -\sin \mathrm{\eta }& \cos \mathrm{\eta }& 0\\ 0& 0& 1\end{array}\right]\approx \left[\begin{array}{ccc}1& \mathrm{\eta }& \mathrm{\phi }\\ -\mathrm{\eta }& 1& -\mathrm{\theta }\\ -\mathrm{\phi }& \mathrm{\theta }& 1\end{array}\right]$$

(10)

The point (U₀, V₀, 0) on an aligned U’–V’–W’ local system can be translated to (U₀+ V₀η, –U₀η + V₀, –U₀φ + V₀θ) for a misaligned detector, with respect to the original coordinate system. (U₀ + V₀η, –U₀η + V₀, –RDD –U₀φ + V₀θ) are their respective global Cartesian coordinates, where RDD = SDD – SRD. The small z coordinate’s variation cause the varied magnification factor of four quadrants on a misaligned detector.

M_l and M_u are the respective magnification factors of the beads below and above z-axis:

$$M_l=M(1-{\displaystyle \frac{h_l\times \mathrm{\theta }}{SRD}});M_u=M(1-{\displaystyle \frac{h_u\times \mathrm{\theta }}{SRD}})$$

(11)

where h_l (h_l< 0) and h_u are the bead farthest below and farthest above z-axis, respectively. Therefore, Eq. (9), the average difference ΔS_d of the vertical projected separation between two adjacent beads’ rotating centers below and above z-axis is:

$$\Delta S_d=(M_u-M_l)\times d\prime ={\displaystyle \frac{M(h_u-h_l)\times d\prime \times \theta }{SRD}}={\displaystyle \frac{S_d\times (phantom\prime sheight)\times \theta }{SRD}}$$

(12)

$$\theta ={\displaystyle \frac{\Delta S_d}{S_d}}\times {\displaystyle \frac{SRD}{(\mathrm{phantom}\prime \mathrm{s} \mathrm{height})}}$$

(13)

where N is total bead number, (N-1)*d’ ≤ (h_u–h_l) ≤ ( N+1)*d’, and (h_u–h_l) ≈ phantom’s height. The tilt angle θ of detector is acquired from Eq. (13). N*d’ is set as the phantom’s height, since the calculated θ value and the setting value of the model are closer. Eq. (13) then acquires the tilt angle θ.

In a misaligned detector, the bead’s positions A and B have distinct magnification factors caused by the tilt angle θ and skew angle φ. The positions H and I are their respective projection on the detector. Let M_A and M_B be the magnification factors for the bead’s respective positions of A and B, then

$$M_{\mathrm{A}}={\displaystyle \frac{[SDD+a\mathrm{\phi }-S_{\mathrm{L}}\theta ]}{SRD-\delta }}; M_{\mathrm{B}}={\displaystyle \frac{[SDD-a\mathrm{\phi }-S_{\mathrm{L}}\theta ]}{SRD-\delta }}$$

(14)

In the case η = 0, θ ≠ 0, and φ ≠ 0, then the two projected elliptical vertex points H and I on the detector are as follows:

$$(U_{\mathrm{H}}, V_{\mathrm{H}})=(M_{\mathrm{A}}\sqrt{r^2-{\delta }^2},M_{\mathrm{A}}h);\hspace{1em}(U_{\mathrm{I}}, V_{\mathrm{I}})=(-M_{\mathrm{B}}\sqrt{r^2-{\delta }^2},M_{\mathrm{B}}h)$$

(15)

$$\mathrm{\eta }*\approx \tan \mathrm{\eta }*={\displaystyle \frac{(V_{\mathrm{H}}-V_{\mathrm{I}})}{(U_{\mathrm{H}}-U_{\mathrm{I}})}}={\displaystyle \frac{(M_{\mathrm{A}}-M_{\mathrm{B}})h}{(M_{\mathrm{A}}+M{}_{\mathrm{B}})\sqrt{r^2-{\delta }^2}}}\approx {\displaystyle \frac{h}{SRD-\delta }}\mathrm{\phi }={\displaystyle \frac{S_{\mathrm{L}}}{SDD}}\mathrm{\phi }$$

(16)

A skew angle φ in a misaligned detector will cause a slant angle η* = φ×S_L/SDD of this projected ellipse. As this rotating bead is on the system’s principal axis, a horizontal line on the detector is traced out by the rotating bead, therefore η* = 0. As η ≠ 0, then η is the rotational angle about “W’-axis”, which is normal to the plane of detector, therefore η* = η at the center G. The slant angle η* in Eq. (16) is necessary to add “η” to the original η* in each ellipse, then:

$$\mathrm{\eta }*\approx \mathrm{\eta }+{\displaystyle \frac{S_{\mathrm{L}}}{SDD}}\mathrm{\phi }$$

(17)

which shows η* is linear to S_L and detector’s skew angle φ. The slant angles η* for all projected ellipses are plotted against ‘S_L’ of each respective bead to acquire the detector’s skew angle φ from the plot’s slope given by Eq. (17), and the slant angle η at the center G.

Since the first order angular approximation is used, there is larger difference between the calculated misalignment and the actual data values at larger misalignment ($>$10°). The PAL effectively reduces those differences.

The misaligned system is calibrated by the resulting angles (θ, φ, η) to acquire a new set of the near ideal detector’s images. After that, we acquire the center G, (U_G, V_G), and M, the CBCT system’s magnification factor, and SRD is re-acquired from Eq. (4), but we measure SDD from the CBCT system directly.

3. Validations

3.1. The Criterion for the Validation of Analytical Simulation

PAL’s ability to detect misalignment and to restore alignment was tested on a detector of a dental CBCT having 2176×1792 pixels, with a pixel size of 0.139 mm in both dimensions. The SRD and SDD are respectively set at 430 mm and 620 mm, and SRD is re-acquired by Eq. (4). The inputs and results of these tests are presented in Secs. 3.2 and 3.3, and summarized in Tables 1 and 2, respectively. The input setting values of calibrated parameters are in bold shown in the row labeled “Model”. In the row labeled “Found” are the misalignment found by the PAL’s iterative process. PAL can accurately yield all the misalignment angles, SRD, and the geometric center (U_G, V_G) at once. Alignment is reached when (θ, φ, η) < 0.002° (3×10^-4 rad) after calibration, as listed in the row labeled “Aligned”. For this setup, all misaligned angles will be < 0.002° after calibration.

3.2. The Helical-Beads Phantom’s AOP Has a Tilt Angle ζ from the AOR

First, the helical-beads phantom (Figure 3a) has its AOP set parallel to the AOR but offset by c = 10.0 mm (Figure 4a), after that, the AOP is tilted by a small angle ζ. Thus, the rotating radii of the beads are set from 30.0 mm to (50.0+100.0*sinζ) for all tests. The results of nine analytical simulations are exposed in Table 1. The errors in SRD are within 0.1 %, and that in (U_G, V_G) is within 0.5 pixels of the setting values. The errors of the misaligned angles (θ, φ, η) are within ±0.03° of the relevant setting values.

3.3. The Line-Beads Phantom’s AOP Has a Tilt Angle ζ from the AOR

The line-beads phantom is exposed in Figure 3c, and its AOP has a small tilt angle ζ from the AOR. The rotating radii of the beads (Figure 4b) are set from 40.0 mm to (40.0 + 100.0*sinζ) mm for all simulations. The results of nine analytical simulations are exposed in Table 2. The errors in SRD are within 0.1%, and that in (U_G, V_G) is within 0.5 pixels of the setting values. The errors of the misaligned angles (θ, φ, η) are within ±0.03° of the relevant setting values.

From the simulated tests of Table 1 and Table 2, we find the two sets of results are almost the same with negligible difference. We conclude the calibrating results are independent of the choice of those two phantoms. The error of the tilt angle ζ between the AOP and the AOR is within ±0.004°of the setting values as ζ ≤ 10°. The error of angle ζ is increasing as ζ increasing, and this error is around 0.04%.

4. Results

4.1. Our Dental CT System

Having validated the efficacy of PAL, we use it calibrate our home-made dental CBCT system with the same geometry in Sec. 3. Ellipses on the detector are traced out by phantom’s rotating beads, which form the PAL’s input data, and the output of final result is the system to remove all misalignment after calibration, i.e., (θ, φ, η) ≈ (0, 0, 0). With variances < 0.002° (3×10^-5 rad).

4.2. Calibrating Our Dental CBCT System

Table 3 lists the elliptical center, (U_L, V_L), semi-axes a and b (in pixels), slant angle η* for each ellipse, r is each bead’s rotating radius, the value of (b/r)×SRD is the projected distance of elliptical center to geometric center G, and V_K the projected distance of beads to the principal axis, all decided from the image CBCT by PAL handling directly from the calibration phantom, thus “before calibration”. In Table 3 and Table 4, each bead’s rotating radius r, varies, as we expected when the AOP has a small tilt from the AOR.

Figure 6 is derived from Table 3. Figure 6a,b yield (U_G, V_G); Figure 6c and Eq. (13) yield θ = -0.277°; Figure 6d and Eq. (17) yield φ = 1.169° and η = -0.292°; thus (θ, φ, η) = (-0.277°, 1.169°, -0.292°) is the misalignment from the primary calculation before calibration. Through PAL, (θ, φ, η) = (-0.231°, 1.152°, -0.291°) is the misalignment from the final calculation. After calibration, we converted the misaligned parameters into "ideal aligned parameters", and final results are exposed in Table 4 and Figure 7. From Table 4, the tilt angle ζ between AOR and the AOP can be determined:

ζ = sin^-1{(51.821 - 53.646)/[(11 - 1) * 10)] } = -1.046° (±0.001°).

After calibration, the geometric center G, (U_G, V_G) = (1028.706, 885.237) in pixels, and the CBCT magnification factor M is 1.4328 from Eq. (4), thence the SRD is 432.474 mm; Eq. (13) gives θ=(-5×10^-4)°; Eq. (17) gives φ=(1×10^-8)° and η=(-3×10^-8)°. The misalignment of the detector are effectively reduced from (θ, φ, η) = (-0.231°, 1.152°, -0.291°) before calibration to (θ, φ, η) = ((-5×10^-4)°, (1×10^-8)°, (-3×10^-8)°) after calibration by PAL, and the slant angles η* of all ellipses are less than 4×10^-5 rad (0.003°). The CBCT system reached within ±0.001° in all three directions of the ideal alignment after calibration.

4.3. Calibrating a Dental CBCT System with Analytical Simulation

An imitation CBCT system, similar to our dental CBCT, is calibrated to validate. The results of simulations is exposed in Table 5. The errors in SRD are within 0.01%, and that in (U_G, V_G) is within 0.02 pixels of the setting values. The errors of the misaligned angles (θ, φ, η) are within ±0.001° of the relevant setting values. The error of the tilt angle ζ between the AOP and the AOR is within ±0.001°of this setting value. The errors of measurement uncertainty in a real CBCT system are not included.

4.4. The Reconstructed Images of the HA Phantom

The HA phantom is designed to analyze the CT reconstruction quality. The reconstructed image of the HA phantom is exposed on Figure 8. There is blur circular edge for the reconstructed image without calibration. We do SNR and CNR analysis for image in CS region (slice number 365). Those two values are improved after calibration. The results are exposed in Table 6 and Table 7.

5. Discussions and Conclusions

Yang et al. [3] gave the errors of misaligned angles (θ, φ, η) within 0.04° and 0.1% error in SDD, but need precisely set beads’ locations and configurations, and the AOP and the AOR overlap. Our previous work [6] required only the two axes be parallel but not overlap, and the beads’ rotating radius > 10mm.

In this paper, we have revised and validated PAL, our precise calibrating algorithm, by using analytical simulations for the setting models and to calibrate a CBCT system by PAL with cylindrical phantoms of a helical-beads (Sec. 3.2) and a line-beads (Sec. 3.3), under three conditions: 1) the detector’s misaligned angles (θ, φ, η) be each<10°, to fit the first-order approximation of the trigonometric functions in Eq. (10), the transformation matrix T; 2) all critical knowledge of misalignment comes from the parameters describing the ellipses projected on the detector by the beads with varied radii of rotation about the AOR; and 3) the AOP needs not be parallel to the AOR, but any small tilt angle ζ between the two axes can be determined.

We found our new procedure to: 1) calibrate the CBCT system more precisely than our previous work [6] in the errors of misaligned angles (θ, φ, η); 2) determine the tilt angle between the AOP and AOR, when not 0°; and 3) work equally well with different phantoms. We conclude that our generalized procedure presented above can give a complete, precise characterization of any CBCT system including a tilting sample platform.