X-ray CT geometrical calibration via locally linear embedding

Abstract

For X-ray computed tomography (CT), geometric calibration and rigid patient motion compensation are inter-related issues for optimization of image reconstruction quality. Non-calibrated system geometry and patient movement during a CT scan will result in streak-like, blurring and other artifacts in reconstructed images. In this paper, we propose a locally linear embedding based calibration approach to address this challenge under a rigid 2D object assumption and a more general way than what has been reported before. In this method, projections are linearly represented by up-sampled neighbors via locally linear embedding, and CT system parameters are iteratively estimated from projection data themselves. Numerical and experimental studies show that images reconstructed with calibrated parameters are in excellent agreement with the counterparts reconstructed with the true parameters.

Keywords

Geometric calibration rigid patient motion locally linear embedding (LLE)

1 Introduction

In the medical imaging field, X-ray computed tomography (CT) provides critical diagnostic information. Over the past decades, CT techniques have been greatly advanced for high-quality images at low radiation dose [1, 2]. One of the challenging problems is geometric calibration such as for C-arm CT [3] and ultra-high resolution CT [4]. This problem is also related to rigid patient motion compensation, since motion is relative between imaging components and a patient body. The Holy Grail for this type of technology is to extract system geometric parameters directly from projection data under practical conditions.

To perform geometric calibration and motion correction, a number of methods were proposed in the past decades. These methods can be analytic with a calibration phantom, and iterative with or without a calibration phantom. Analytic methods are widely used in industrial CT, one of which is based on the identification of elliptical parameters in cone-beam geometry [5]. On the other hand, more calibration methods are iterative, such as optimization-based calibration for cone-beam CT [6], self-calibration for cone-beam CT [7], and self-calibration for cone-beam micro-CT [8]. There is an overlap in the literature on geometric calibration and motion reduction. While some motion reduction methods utilize fast scanning, even with multi-source-detector systems [9, 10] or avoid motion affected data [11, 12], the other methods estimate patient motion and compensate for its effect [13 –17].

Inspired by the recently developed locally linear embedding (LLE) theory, here we propose a geometric calibration method by coupling LLE and re-projection [18, 19]. In our method, we start with an initial parameterization to reconstruct a CT image, calculate re-projected projection vectors that sample the parametric range densely. With the re-projected projection vectors and the original projection vectors, we calculate weight coefficients and neighbors by LLE, and update the parameter estimation. An image is iteratively reconstructed until a satisfactory quality is achieved. Being focused on fan-beam imaging geometry, the rest of this paper is organized as follows. In the next section, we introduce LLE and re-projection, and describe our calibration method. In the third section, we report numerical and experimental results. In the last section, we discuss relevant issues and conclude the paper.

2 Material and methods

2.1 Locally linear embedding

In 2000, Roweis and Saul proposed an unsupervised manifold learning algorithm called “locally linear embedding” (LLE) for dimensionality reduction. Related to classic dimensionality reduction methods such as principal component analysis (PCA) and multidimensional scaling (MDS), LLE computes eigenvectors but does so locally to preserve intrinsic geometry embedded in a high dimensional space. LLE is easy to implement and yet gives excellent performance when data are sufficiently dense. Specifically, the LLE consists of three steps: Step 1 is to find the K nearest neighbors for each high-dimensional point in terms of the Euclidean distance; Step 2 is to represent each point linearly with its neighbors and calculate weight coefficients for its neighbors; and Step 3 is to map high dimensional data with the weight coefficients to a low dimensional representation on an intrinsic manifold.

In this study, with the data vector b_i, the first step of LLE is to find its K nearest neighbors in its dataset vectors ${\tilde{b}}_{ij}$ according to the Euclidean distance $d_{ij} = {∥ b_{i} - {\tilde{b}}_{ij} ∥}_{2}^{2}$ (1)

With the K nearest vectors, the second step is to represent the original data vector linearly with its neighboring vectors: $b_{i} = \sum_{k = 1}^{K} w_{ik} {\tilde{b}}_{ik}$ (2) where ${\tilde{b}}_{ik}$ are the K nearest neighbors and w_ik are respectively weight coefficients, and ∑_kw_ik = 1. The weight matrix is solved by minimizing the following error: $min {∥ b_{i} - \sum_{k} w_{ik} {\tilde{b}}_{ik} ∥}_{2}^{2} \begin{matrix} s . t . \sum_{j} w_{ik} = 1 \end{matrix}$ (3)

Using the constraint ∑_kw_ik = 1, to solve Equation (3) is equivalent to solve the linear system $\sum_{k} c_{jk} w_{ik} = 1$ (4) where C = (c_jk) is the local covariance matrix calculated as $C = (b_{i} - {\tilde{b}}_{ik})^{T} (b_{i} - {\tilde{b}}_{ik})$ (5)

With the weight coefficients W = (w_ik), the third step is to calculate the global internal coordinate Y = (y_i) by solving the equation $\min_{Y} \sum_{i} {∥y_{i} - \sum_{k} w_{i k} y_{k}∥}_{2}^{2}$ (6)

By the Rayleitz-Ritz theorem [20], the solution of Equation (6) is given by the bottom d + 1 eigenvectors of the generalized eigenvalue problem $MY = λ Y$ (7) where λ is the eigenvalue, and M = (I - W) ^T (I - W).

2.2 Fan-beam geometry

Fan-beam geometry is popular with commercial spiral CT systems, as shown in Fig. 1 where the origin of the Cartesian coordinate system is the nominal center of an object to be scanned and reconstructed, and the X-ray source is intended to be on a circular trajectory whose center is the origin. The source-object distance (SOD) is the distance between the X-ray source and the system origin which is the radius of the scanning circle. The X-ray source emits fan-beam X-rays covering the object. The central ray in the fan-beam is perpendicular to a linear detector array of length L which receives X-rays after being attenuated by the object, the angle between a ray of interest and the central ray is γ. The object-detector distance (ODD) is the distance between the system origin and the detector midpoint. θ denotes the x-ray scanning angle, and β = θ - γ.

Because the X-ray source is on a circular trajectory, it is convenient to represent fan-beam geometry in a polar coordinate system. Various practical factors will lead to inaccurate geometric parameters in a view-dependent fashion. Specifically, inaccurate projection angles and other parameters are ${\tilde{θ}}_{i} = θ_{i} + δ_{i}, i = 1, 2, \dots, N$ (8)

${\tilde{p}}_{j} = p_{j} + ɛ_{j}, j = 1, 2, \dots, M$ (9) where θ_i are the accurate projection angles, p_j are the other accurate geometric parameters including source-object distance, object-detector distance, detector offset, detector tilt angle, and so on, which do not depend on θ_i. N is the number of views, M is the number of other parameters, δ_i and ɛ_j are respectively the angular error and the parametric error. Generally speaking, more accurate geometric parameters will help to reconstruct a better CT image.

To be general, a rigid patient motion problem is included in the geometric calibration problem. In fan-beam geometry as shown in Fig. 2, patient motion can be represented with an object x coordinate offset x_i, an object y coordinate offset y_i, and a projection angular error, while the projection angle is still expressed by Equation (8).

2.3 Calibration with LLE

Typically, an iterative CT reconstruction approach can be modeled to solve the following linear system of equations: $Au = b_{i}$ (10) where u = (u₁, u₂, …, u_J) is an image represented as a J dimensional vector, J is the number of pixels, b_i = (b₁, b₂, …, b_L) is projection data, L is the number of vector elements, and A = (a_jk) is a projection matrix related to the geometric parameters. If both the system matrix and projection data are known, an iterative algorithm [20] reconstructs a CT image with fixed projection matrix by reducing the difference between original and re-projected data. In the reconstruction algorithm, an important step is to calculate the projection matrix which is affected by the geometric parameters, such as projection angle and source-object distance; that is, $A = A (P)$ (11) where P is an estimated parameter vector, which include source-object distance, object-detector distance and detector offset distance, detector tilt angle, projection angle, and so on.

In this study, we calculate projections using the distance-driven method [21]. Based on the re-projection approach in the iterative reconstruction, we present a new iterative method to estimate the geometric parameters by minimizing the mean squared error between the projection data and re-projected projection data, which can be formulated as $P = arg min {∥ b_{i} - {\tilde{b}}_{ij} ∥}_{2}^{2} s . t . A (P) u = b_{i}$ (12) where b_i is the projection vector obtained from the measurement along different projection views, ${\tilde{b}}_{ij}$ is the corresponding re-projected projection vector from a reconstructed image with sampled parameters, and P is the updated vector of parameters. The re-projected projection vector is calculated by Equations (10) and (11) within a densely sampled parametric range: ${\tilde{P}}_{j} = (p_{j 1}, p_{j 2}, \dots, p_{jn})$ (13)

If the parametric sampling interval is sufficiently small, a true parameter vector is close to neighboring sampled parameter vectors, and the measured projection vector can be linearly expressed by the K nearest re-projected projection vectors associated with the sampled parameter vectors. That is, $b_{i} = \sum_{k = 1}^{K} w_{ik} {\tilde{b}}_{ik}$ (14) $P = \sum_{k = 1}^{K} w_{ik} {\tilde{P}}_{k}$ (15) where ${\tilde{b}}_{ik}$ are the K nearest re-projection vectors associated with the corresponding K vectors of parameters ${\tilde{P}}_{k}$ and w_ik are weight coefficients. The key relationship is that the weight coefficients for Equations (14) and (15) are the same. Therefore, the real parameter estimation can be refined by searching for the K nearest re-projected projection vectors and updating the parameter vector with the weight coefficients and the corresponding sampled parameters. Consequently, the geometric calibration problem is solved by dimensionality reduction via locally linear embedding.

With a densely sampled parametric domain and correspondingly re-projected projection vectors, we can find the K nearest re-projected projection vectors of an original projection vector and calculate the weight coefficients with LLE by Equation (1) and (4), (5). With the sampled parameters and corresponding weight coefficients, we can perform a parametric update according to Equation (15). Specifically, our calibration method is flowcharted in Fig. 3.

3 Results

First, the proposed LLE-based calibration approach was evaluated using an abdomen image phantom from clinic as shown in Fig. 4. The size of the abdomen phantom is 512 × 512. With this phantom, projection data with angular and other geometric errors were generated. Then, we applied our algorithm to the projection data with geometric parameter errors. Furthermore, we tested our algorithm in a rigid patient motion problem. The computational environment is MatLab 2010a on a computer with the CPU Intel Core 2 Duo E7600 @3.06 GHz, RAM 4.00GB, 64-bit OS.

To evaluate reconstructed images, we utilized the universal quality index (UQI) [22]. Also, we quantified geometric calibration results with the average angular error ${\bar{θ}}_{erroe}$ , parametric error P_Error, and average object coordinate offsets $\bar{x}$ and $\bar{y}$ .The UQI evaluates an image by integrating three factors including correlation distortion, brightness distortion and contrast distortion. The range of UQI is between [–1 1]. The closer to 1 the UQI, the better a reconstructed image will be. Given a reconstructed image $u_{ij}^{*}$ and the ground truth image u_ij with size of S × T, the UQI is defined as $UQI = \frac{4 σ_{{uu}^{*}} \bar{u} \times {\bar{u}}^{*}}{(σ_{u}^{2} + σ_{u^{*}}^{2}) [(\bar{u})^{2} + ({\bar{u}}^{*})^{2}]}$ (16) where $\bar{u} = \frac{1}{S \times T} \sum_{0 \leq i < S} \sum_{0 \leq j < T} u_{ij}$ (17) ${\bar{u}}^{*} = \frac{1}{S \times T} \sum_{0 \leq i < S} \sum_{0 \leq j < T} u_{ij}^{*}$ (18) $σ_{u}^{2} = \frac{1}{S \times T - 1} \sum_{0 \leq i < S} {\sum_{0 \leq j < T} (u_{ij} - \bar{u})}^{2}$ (19) $σ_{u^{*}}^{2} = \frac{1}{S \times T - 1} \sum_{0 \leq i < S} {\sum_{0 \leq j < T} (u_{ij}^{*} - {\bar{u}}^{*})}^{2}$ (20) $σ_{{uu}^{*}} = \frac{1}{S \times T - 1} \sum_{0 \leq i < S} \sum_{0 \leq j < T} (u_{ij}^{*} - {\bar{u}}^{*}) (u_{ij} - \bar{u})$ (21)

The average angular error is ${\bar{θ}}_{Error} = \frac{\sum_{i = 1}^{N} | \tilde{θ_{i}} - θ_{i} |}{N}$ (22) where ${\tilde{θ}}_{i}$ is the projection angle after calibration, and θ_i is the original projection angle. The parametric error is $P_{Error} = | P - P_{Ori} |$ (23) where P is the parameter vector after calibration, and P_Ori is the real parameter vector. The average object x and y coordinate offsets are $\bar{x} = \frac{\sum_{i = 1}^{N} | x_{i} - {\tilde{x}}_{i} |}{N}$ (24) $\bar{y} = \frac{\sum_{i = 1}^{N} | y_{i} - {\tilde{y}}_{i} |}{N}$ (25) where x_i and y_i are the original object x and y coordinate offsets, while ${\tilde{x}}_{i}$ and ${\tilde{y}}_{i}$ are the calibrated offsets respectively.

The image reconstruction process was performed using the Ordered Subset Simultaneous Algebraic Reconstruction Technique [23] (OS-SART) with Total Variation (TV) regularization [24]. As the reconstruction and calibration iteration stop criteria, we monitored for satisfactory UQI values in both the image and projection domains.

3.1 Geometric calibration

In this aspect, we evaluated the utility of our calibration approach in reducing projection angular and other geometric errors. The number of projection angles N = 360. The projection angles were randomly perturbed as follows: $θ_{i} = 360 \times (i - 1) / N + δ_{i}, i = 1, 2, \dots, N$ (26)

The other parameters are in Table 1, including the real values, initial values and final estimates. The length of the detector array is 80 cm. With the geometric parameters, projection data were generated with Poisson noise, assuming that the number of incoming x-ray photons was 10⁵.

To calibrate the parameters efficiently, we first calibrated the geometric parameters one by one rather than calibrating all the parameters together. The calibration process proceeded in the following sequence: detector offset, source object distance, object detector distance, detector tilt and projection angle. The angular sampling step for each parameter were 0.04 cm, 0.1 cm, 0.1 cm, 0.01° and 0.02° respectively, the sampling ranges were [0 cm 0.4 cm], [46 cm 54 cm], [46 cm 54 cm], [0° 2°] and [-1° 1°] accordingly, and the number of nearest neighbors was 2. The UQI threshold for stopping was 0.9 for reconstruction and 0.6 for calibration.

After 4 calibration iterations, the intermediate and final reconstructed images are in Fig. 5, the calibrated parameters are in Table 1. The UQI values of reconstructed images, the computation cost and accumulated cost at angle 60° are in Fig. 6, after different numbers of calibration iterations. The average angular errors and other calibrated geometric parameters are in Fig. 7, for different numbers of calibration iterations. The projection angular error distributions are in Fig. 8 with respect to the number of calibration iterations.

3.2 Geometric calibration with real projection data

Next, we applied our calibration algorithm in a case of real projection data with geometric parameter errors. The nominal SOD and ODD values were both 38 cm, the number of detector elements was 1024, and the length of detector was 13.0048 cm. The number of projection angles was 900 in the angular range [0°, 360°]. An object was on the rotation stage while the rotation center was not on the line between the X-ray source and the detector center, as shown in Fig. 9 where the red spot was the rotation center not on the system origin. Therefore, we used our algorithm to calibrate the rotation center offset. The sampling range was [–1 cm 1 cm] with a sampling step 0.05 cm. The number of nearest neighbors was still 2. The UQI stopping threshold was again 0.9 for reconstruction and 0.6 for calibration. After 5 calibration iterations, the reconstructed images with different iterations are in Fig. 10, the computation cost, accumulated cost at angle 60° and calibrated rotation center offset are in Fig. 11, after different numbers of calibration iterations. The calibrated rotation center offset is -0.3847 cm.

3.3 Calibration with rigid patient motion

As aforementioned, the rigid patient motion could be represented as object x coordinate offset error, object y coordinate offset error and projection angle error. Therefore, in this part we try to simulate with rigid patient motion calibration with abdomen phantom. The number of projection angle and other parameters are same with part 1. The object x coordinate offset error and object y coordinate offset error are all randomly selected in [–1 cm 1 cm], and projection angle error are randomly chose in [–1° 1°].

The calibration was in the following sequence: object x coordinate offset, object y coordinate offset and projection angle. The sampling rate for each parameter are 0.05 cm and 0.1°, the sampling ranges are [–1 cm 1 cm] and [–1° 1°], the number of nearest neighbors K = 2. The UQI stop value for reconstruction is 0.9 and for calibration is 0.8. After 20 calibration iterations, the reconstructed images with different iterations as shown in Fig. 12, the RMSE, UQI values of reconstructed images, average angel error, average object x coordinate offset error, the computation cost and accumulated cost at angle 60 as shown in Fig. 13, for different numbers of calibration iterations. Object x coordinate offset error distributions, object x coordinate offset error distributions and projection angular error distributions are in Figs. 14, 15 and 16 respectively, for different numbers of calibration iterations.

4 Discussions and conclusion

Our main contribution in this paper is a novel approach for 2D CT system calibration of view-wise random geometric parameters. The recently established dimensionality reduction theory LLE is used as the most important step to estimate geometric parameters subject to the inherent low-dimensional consistency, which has been demonstrated successful in numerical simulation and real experiments. LLE can find the K nearest re-projected projection vectors with the corresponding sampled parameters and update geometric parameters based on the linear combination of sampled parameters with weighting coefficients calculated via LLE.

Geometric parameters have different sensitivities on reconstruction quality. In the geometric calibration, the detector offset affects reconstruction very much. The detector tilt and projection angular error have stronger effects than the source-object distance error and the object-detector distance errors on reconstruction quality. This is also why the calibration results for detector offset, detector tilt and projection angle are better than that for source-object and object-detector distances. In the patient motion calibration, the object x and y coordinate offsets seemed more effective on a reconstructed image quality than the projection angular error.

A closely-related issue is the parametric sampling rate for re-projection. The greater the sampling rate, the more accurate the calibration results will be, but the higher computational cost will be involved. We need to choose a proper sampling rate to balance between calibration results and computational overhead.

An important issue in using TV is its weight. If it is too small, TV will not be able to reduce artifacts and noise. If it is too large, TV will over-smoothen CT images. The TV parameter depends on levels of artifacts and noise. In this study, we empirically set the TV parameter to 0.1 after geometric calibration. We could use a smaller TV parameter in a case of weaker noise.

In conclusion, we have proposed a unified image reconstruction and parameter estimation scheme in which both the projection matrix and an underlying image are updated iteratively by re-projection and LLE. This method can be applied to CT system calibration and rigid patient motion compensation. Realistic simulation studies and real experiments have been performed to verify the proposed method, with an encouraging performance. There are several issues worth further investigation, such as extending the results into the cases of cone-beam geometry and non-rigid patient motion.

Footnotes

Acknowledgments

This work was supported by NIH under Grants R01 EB016977 and U01 EB017140, and partly by the National Natural Science Foundation of China (No. 61201346, No. 61401049), the Fundamental Research Funds for the Central Universities (No. CDJZR14125501) and Graduate Innovative Research Projects of Chongqing (No. CYB14024).

References

Wang

, Yu

and Man

B.D.

, An outlook on x-ray CT research and development, Medical physics 35(3) (2008), 1051–1064.

Wang

, Zhao

, Yu

et al., Design, analysis and simulation for development of the first clinical micro-CT scanner, Academic Radiol 12(4) (2005), 511–525.

Rock

, Linsenmaier

, Brandl

et al., Introduction of a new mobile C-arm/CT combination equipment (ISO-C-3D). Initial results of 3-D sectional imaging, Der Unfallchirurg 104(9) (2001), 827–833.

Rajiv

, Grasruck

, Suess

et al., Ultra-high resolution flat-panel volume CT: Fundamental principles, design architecture, and system characterization, European Radiology 16(6) (2006), 1191–1205.

Noo

, Clackdoyle

, Mennessier

et al., Analytic method based on identification of ellipse parameters for scanner calibration in cone-beam tomography, Phys Med Biol 45 (2000), 3489–3508.

Panetta

, Belcari

, Guerra

A.D.

et al., An optimization-based method for geometrical calibration in cone-beam CT without dedicated phantoms, Med Biol 53 (2008), 3841–3861.

Wein

, Ladikos

and Baumgartner

. Self-calibration of geometric and radiometric parameters for cone-beam computed tomography. in: Fully 3D, Potsdam, (2011), 327–330.

Patel

, Chityala

R.N.

, Hoffmann

K.R.

et al., Self-calibration of a cone-beam micro-CT system, Med Phys 36 (2009), 48–58.

Man

B.D.

, Basu

, Bequé

et al., Multi-source inverse geometry CT: A new system concept for X-ray computed tomography. Medical Imaging. International Society for Optics and Photonics (2007).

10.

Liu

, Liu

, Wang

et al., Half-scan cone-beam CT fluoroscopy with multiple sources, Med Phys 28(7) (2001), 1466–1471.

11.

Hsieh

, Computed Tomography, Principles, Design, Artifacts, and Recent Advances. Bellingham. WA: SPIE, 2003.

12.

Crawford

C.R.

and King

K.F.

, Computed-tomography scanning with simultaneous patient translation, Med Phys 17(6) (1990), 967–982.

13.

Wang

and Vannier

M.W.

, Preliminary study on helical CT algorithms for patient motion estimation and compensation, IEEE Trans Med Imag 14(2) (1995), 205–211.

14.

W.G.

and Mackie

T.R.

, Tomographic motion detection and correction directly in sinogram space, Phys, Med, Biol 47(8) (2002), 1267–1284.

15.

and Wang

, Data consistency based rigid motion artifact reduction in fan-beam CT, Medical Imaging, IEEE Transactions on 26(2) (2007), 249–260.

16.

, Wei

, Hsieh

et al., Data consistency based translational motion artifact reduction in fan-beam CT, Medical Imaging, IEEE Transactions on 25(6) (2006), 792–803.

17.

Leng

, Nett

and Speidel

et al., Motion artifact reduction in fan-beam and cone-beam computed tomography via the fan-beam data consistency condition (FDCC). Medical Imaging. International Society for Optics and Photonics, 2007.

18.

Roweis

S.T.

and Lawrence

S.K.

, Nonlinear dimensionality reduction by locally linear embedding, Science 290(5500) (2000), 2323–2326.

19.

Roweis

S.T.

and Lawrence

S.K.

, An introduction to locally linear embedding. unpublished. Available at: http://www.cs.toronto.edu/~roweis/lle/publications.html (2000).

20.

Horn

R.A.

and Johnson

C.R.

, Matrix Analysis. Cambridge University Press, Cambridge, 1990.

21.

Siddon

R.L.

, Fast calculation of the exact radiological path for a three-dimensional CT array, Medical Physics 12(2) (1985), 252–255.

22.

Wang

, Bovik

A.C.

, A universal image quality index, signal processing letters, IEEE 9(3) (2002), 81–84.

23.

Wang

and Jiang

, Ordered-subset simultaneous algebraic reconstruction techniques (OS-SART), Journal of X-ray Science and Technology 12(3) (2004), 169–177.

24.

Sidky

Y.E.

, Kao

C.M.

and Pan

, Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT, Journal of X-ray Science and Technology 14(2) (2006), 119–139.