IMAGE FOR LUNG CANCER RADIOTHERAPY

CROSS REFERENCE TO RELATED APPLICATION

[0001] This application claims priority to U.S. Provisional Patent Application No.

61/325,700, filed April 19, 2011, entitled "Real-Time Volumetric MRI Image Reconstruction and 3D Tumor Localization Based on a Single X-Ray Projection Image for Lung Cancer Radiotherapy" . The entire content of the before-mentioned patent application is incorporated by reference as part of the disclosure of this application.

BACKGROUND

[0002] This application relates to devices and techniques that use radiotherapy techniques.

[0003] Several categories of lung tumor localization methods are available to identify lung tumor location during radiotherapy treatment including direct tumor localization, localization via breathing surrogates, and localization of implanted markers. Radiotherapy treatment machines equipped with an on-board imaging system, which includes one kV x-ray source and one flat panel imager can be used to track lung tumors directly without implanted markers in

fluoroscopic images acquired by the on-board imaging system. Direct localization techniques are mainly for images acquired from the anterior-posterior (AP) direction and the tumors with well-defined boundaries. Also, direct localization techniques involve acquisition and labeling of the training fluoroscopic data prior to treatment. In direct tumor localization, the third dimension of tumor motion perpendicular to the x-ray imager cannot be resolved. For localization based on implanted markers, invasive procedure of marker implantation has an associated risk of pneumothorax. Surrogate-based localization sometimes involves no additional radiation dose to the patient (e.g., optical localization) and comes with a varying relationship between tumor motion and surrogates both intra- and inter- fractionally.

SUMMARY

[0004] Techniques and systems and apparatus are disclosed for implementing real-time volumetric image reconstruction and 3D tumor localization based on a single x-ray projection image for lung cancer radiotherapy. [0005] In one aspect, a method for implementing real-time volumetric image reconstruction and 3D tumor localization includes performing deformable image registration between a reference phase and remaining N-l phases for a set of volumetric images of a patient at N breathing phases to obtain a set of N-l deformation vector fields (DVFs). The set of DVFs are represented using eigenvectors and coefficients obtained from principle component analysis (PC A). The PCA coefficients are varied to generate new DVFs. The generated new DVFs are applied on a reference image to generate new volumetric images.

[0006] Implementations can optionally include one or more of the following features. The method can include generating a reconstructed volumetric image from a single projection image by optimizing the PCA coefficients such that a computed projection image matches a measured one of the new volumetric images. The method can include deriving a 3D location of a tumor by applying an inverted DVF on a position of the tumor in the reference image. The method can include using x-ray projections with rotational geometry. The method can include using patient's 3D surface images measured with a range camera. The method can include using one or multiple markers placed on the patient's chest or abdomen. The method can include using one or multiple fiducial markers implanted inside the patient. The method can include using patient's anatomic features such as diaphragm. The method can include using pressure measured with a belt. The method can include using lung volume measured with a spirometer. The method can include using lung air flow measured with a spirometer. The method can include incorporating respiratory motion. The method can include using fixed-angle geometry, such as fluoroscopy.

[0007] In another aspect, a method for implement real-time volumetric image reconstruction and 3D tumor localization based on a single x-ray projection image for lung cancer radiotherapy includes generating a parameterized PCA model. Generating the parameterized PCA model includes approximating a DVF relative to a reference image as a function of space and time by a linear combination of the sample mean vector and a few eigenvectors corresponding to the largest eigenvalues.

[0008] Implementations can optionally include one or more of the following features. The method can include generating a set of optimal PCA coefficients such that the projection of the corresponding volumetric image matches with the measured x-ray projection. The method can include identifying a tumor position, which can include using a push forward DVF, which is defined on the coordinate of the reference image and tells how each voxel in the reference image moves by adopting an efficient fixed-point algorithm for deformation inversion and calculating tumor position.

[0009] In another aspect, a graphics processing unit (GPU) is configured to perform operations including performing deformable image registration between a reference phase and remaining N-1 phases for a set of volumetric images of a patient at N breathing phases to obtain a set of N-1 deformation vector fields (DVFs). The GPU is configured to represent the set of DVFs using eigenvectors and coefficients obtained from principle component analysis (PC A). The GPU is configured to vary the PC A coefficients to generate new DVFs. The GPU is configured to apply the generated new DVFs on a reference image to generate new volumetric images.

[0010] Implementations can optionally include one or more of the following features. The GPU can be configured to generate a reconstructed volumetric image from a single projection image by optimizing the PCA coefficients such that a computed projection image matches a measured one of the new volumetric images. The GPU can be configured to derive a 3D location of a tumor by applying an inverted DVF on a position of the tumor in the reference image.

[0011] In another aspect, a method of manufacturing a GPU includes using an algorithm to perform the following operations using a programming environment selected from a group comprising compute unified device architecture (CUD A), Brook+, OpenCL and DirectCompute. The operations include performing deformable image registration between a reference phase and remaining N-1 phases for a set of volumetric images of a patient at N breathing phases to obtain a set of N-1 deformation vector fields (DVFs); representing the set of DVFs using eigenvectors and coefficients obtained from principle component analysis (PCA); varying the PCA

coefficients to generate new DVFs; and applying the generated new DVFs on a reference image to generate new volumetric images.

[0012] Implementations can optionally include one or more of the following features. The manufactured GPU can perform operations including generating a reconstructed volumetric image from a single projection image by optimizing the PCA coefficients such that a computed projection image matches a measured one of the new volumetric images. The manufactured GPU can perform operations comprising deriving a 3D location of a tumor by applying an inverted DVF on a position of the tumor in the reference image. [0013] The subject matter described in this specification potentially can provide one or more of the following advantages. For example, the described techniques, systems and apparatus can use a single projection image for real-time volumetric image reconstruction and 3D tumor localization.

BRIEF DESCRIPTION OF THE DRAWINGS

[0014] FIGS. 1A and IB are flow charts illustrating a principal component analysis (PCA) lung motion model based process 100 for real-time volumetric image reconstruction and 3D tumor localization based on a single x-ray projection image for lung cancer radiotherapy.

[0015] FIGS. 2A, 2B and 2C are flow charts that summarize an exemplary respiratory motion prediction incorporating process for volumetric image reconstruction and tumor localization based on a single x-ray projection.

[0016] FIG. 3A shows a "measured" projection of a test image at a right posterior oblique (RPO) angle.

[0017] FIG. 3B shows an objective value (top) and a relative image reconstruction error (bottom) at each iteration given the image in FIG 1A.

[0018] FIGS. 4A, 4B and 4C show a sagittal view of an absolute difference image between: A) ground truth test and reference images, B) ground truth test image and image reconstructed using a single projection, and C) ground truth test image and the deformed reference image using demons. Tumor is a round object near the center of the slice.

[0019] FIGS. 4D, 4E and 4F show the same as FIGS. 4A, 4B and 4C, except for coronal view. Tumor is a round object in the right lung.

[0020] FIG. 5 shows on top row: relative image error between the ground truth test image and: reference image; image reconstructed using the proposed algorithm as a function of cone beam projection angle; and on bottom row: same as top row, except for 3D localization error.

[0021] FIGS. 6A, 6B, 6C, 6D, 6E, 6F, 6G, and 6H show the SI tumor position estimated by the respiratory motion prediction incorporating algorithm as well as the ground truth for the digital respiratory phantom for Cases 1 through 8.

[0022] FIGS. 7A and 7B show the estimated and ground truth tumor position in the SI direction for both regular and irregular breathing.

[0023] Results for the five lung cancer patients are shown in FIGS. 8A, 8B, 8C, 8D, 8E, 8F, 8G, 8H, 81 and 8J, where the tumor positions on the axial and tangential directions are shown separately.

[0024] FIGS. 9A, 9B, 9C, 9D, 9E, and 9F show the image reconstruction and tumor localization results for patient 1 at an EOE phase and an EOI phase (indicated in FIG. 8, Case 11), as well as the corresponding cone beam x-ray projections.

[0025] FIG. 10 shows the average localization error on tangential and axial directions as well as computational time versus the downsampling factor.

[0026] FIGS. 11A, 1 IB, 11C, 1 ID, 1 IE and 1 IF show scatter plots between the image intensities of the preprocessed cone beam projection and of the computed projection of the reconstructed image at two different projection angles for patient 1, corresponding to the two projections shown in FIGS. 9A, 9B, 9C, 9D, 9E and 9F.

[0027] FIG. 12 is a block diagram of an exemplary system 1200 for implementing real-time volumetric image reconstruction and 3D tumor localization for lung cancer radiotherapy.

[0028] Like reference symbols and designations in the various drawings indicate like elements.

DETAILED DESCRIPTION

[0029] The techniques, structures and apparatus described in this application can be used to implement real-time volumetric image reconstruction and 3D tumor localization based on a single x-ray projection image for lung cancer radiotherapy. For example, an algorithm can be provided to reconstruct volumetric image and to locate 3D tumor position from a single x-ray projection image. Additionally, implementations can include application of the algorithm on a graphics processing unit (GPU) to obtain real-time (e.g. sub-second) performance. Also, in addition to using the projection image, the described algorithm and respiratory motion prediction incorporating algorithm can also use surface image, surface marker, implanted marker, diaphragm, lung volume, and air flow to reconstruct a volumetric image and to derive the 3D tumor position.

[0030] In one aspect, FIGS. 1A and IB are flow charts illustrating a principal component analysis (PCA) lung motion model based process 100 for real-time volumetric image

reconstruction and 3D tumor localization based on a single x-ray projection image for lung cancer radiotherapy. The process 100 can be implemented using a Graphics Processing Unit (GPU) as described below. Given a set of volumetric images of a patient at N breathing phases as the training data, deformable image registration can be performed between a reference phase and the other N-1 phases, resulting in N-1 deformation vector fields (DVFs) (110). A parametric PCA lung motion model can be obtained (120). These DVFs can be represented efficiently by a few eigenvectors and coefficients obtained from the PCA lung motion model (130). The GPU can perform image reconstruction using the obtained parametric PCA lung motion model (140). The GPU can identify a set of optimal PCA coefficients such that the projection of the corresponding volumetric image matches with the measured x-ray projection (142). In identifying the set of optimal PCA coefficients, the PCA coefficients are initialized to those obtained for the last projection (144). Then the PCA coefficients are varied until the projection of the corresponding volumetric image matches with the measured x-ray projection (146). By varying the PCA coefficients, new DVFs can be generated, which, when applied on the reference image, can lead to new volumetric images. A volumetric image can be construed from a single projection image by optimizing the PCA coefficients such that its computed projection matches the measured one. The 3D location of the tumor can be derived by applying the inverted DVF on its position in the reference image (150). The described algorithm or process can be implemented on graphics processing units (GPUs) to achieve real-time efficiency. The training data using a realistic and dynamic mathematical phantom with 10 breathing phases are described in this document. The testing data were 360 cone beam projections corresponding to one gantry rotation, simulated using the same phantom with a 50% increase in breathing amplitude.

[0031] As described below, the average relative image intensity error of the

reconstructed volumetric images is 6.9% ± 2.4% . The average 3D tumor

localization error is 0.8 mm ± 0.5 mm. On an NVIDIA Tesla C 1060 GPU card, the average computation time for reconstructing a volumetric image from each projection is 0.24 seconds (range: 0.17 and 0.35 seconds). Accordingly, volumetric images can be

reconstructed and tumor positions can be localized in 3D in near real-time from a single x-ray image.

[0032] PCA LUNG MOTION MODEL BASED LUNG CANCER RADIOTHERAPY

[0033] An algorithm based on the PCA lung motion model can be used to reconstruct volumetric images and extract 3D tumor motion information in real time from a single x-ray projection in a non-invasive (no marker implantation required), accurate, and efficient way. The described techniques for volumetric image reconstruction and 3D tumor localization can use 4DCT from treatment simulation as a priori knowledge. Given a set of volumetric images of a patient at N breathing phases as the training data, deformable image registration can be performed between a reference phase and the other N-l phases, resulting in a set of N-l deformation vector fields (DVFs). This set of DVFs can be represented efficiently by a few eigenvectors and coefficients obtained from PCA. By changing the PCA coefficients, a new DVF can be generated. A new volumetric image can be obtained when the DVF is applied on the reference image. A volumetric image can be reconstructed from a single projection image by optimizing the PCA coefficients such that its computed projection image matches the measured one of the new volumetric image. The 3D location of the tumor can be derived by applying the inverted DVF on its position in the reference image.

[0034] In the PCA lung motion model, the DVF relative to a reference image as a function of space and time can be approximated by a linear combination of the sample mean vector and a few eigenvectors correspondin to the largest eigenvalues, i.e.,

where x(t) is the parameterized DVF as a function a space and time, x is the mean DVF with respect to time; u,t are the eigenvectors obtained from PCA and are functions of space; the scalars w _k (f ) are PCA coefficients and are functions of time. The eigenvectors are fixed after

PCA and the temporal evolution of the PCA coefficients can drive the new volumetric image and the dynamic lung motion in real time.

[0035] In fact, the PCA motion model is so flexible that it is capable of represent tumor motion which is beyond that of the training set. For instance, if the PCA coefficients are allowed to change arbitrarily, then given two eigenvectors, the tumor can move anywhere in a plane defined by the corresponding entries in the eigenvectors, and thus hysteresis motion can be handled; given three eigenvectors, the tumor can move anywhere in the space. Each voxel is associated with three entries (corresponding to three canonical directions in space) in each and every eigenvector. Therefore, each eigenvector defines a vector (or a direction) in space along which each voxel moves; two eigenvectors span a plane; and three eigenvectors span the entire 3D space where the tumor can move. This important feature can allow the PCA motion model to capture a wide range of tumor motion trajectories. [0036] In one example, the first three PCA coefficients associated with the largest eigenvalues in the PCA lung motion model can be retained. Because the PCA model attempts to efficiently encode the complex spatiotemporal patterns of the entire lung motion, the selection of three PCA coefficients is based on both temporal and spatial considerations. On the temporal side, by studying the motion trajectory of a single point (tumor), the first three PCA components can be interpreted in the temporal domain as three physical quantities: the mean, velocity, and acceleration of the tumor motion. Although the motion of the entire lung is more complicated due to the high spatial dimensionality, because of the high degree of redundancy and correlation, the entire lung motion may still be well approximated by the three physical components. For a hypothetical lung motion with sinusoidal motion and arbitrary phase for each voxel, the entire lung motion can be completely represented by as few as two PCA components, and for more realistic and complicated lung motion, the PCA spectrum is still dominated by only two or three modes for most breathing patterns. On the spatial side, three is the least number of PCA coefficients that are capable of capturing a full 3D tumor motion trajectory. On the other hand, using more PCA coefficients, although more flexible, could lead to the overfitting problem and may not be necessary.

[0037] After obtaining a parameterized PCA lung motion model, a set of optimal PCA coefficients are identified such that the projection of the reconstructed volumetric image corresponding to the new DVF matches with the measured x-ray projection. Denote f ₀ as the reference image, f as the reconstructed image, y as the measured projection image, and P as a projection matrix which computes the projection image of f . The cost function can be represented using Equation (2):

m ^' m J(w, a, b) = \\P †{x,† ₀ ) - a - y - b U (2)

s. t. x = x + U w

where, U is a matrix whose columns are the PCA eigenvectors and w is a vector comprised of the PCA coefficients to be optimized, x is the parameterized DVF, and p denotes the standard vector / -norm. Because the computed and measured projection images may have different pixel intensity levels, it can be assumed that there exists a linear relationship between them and introduce two auxiliary parameters a and b to account for the differences. The optimization of the cost function with respect to w, a,b is described below. [0038] To optimize the cost function in Eq. (2) with respect to w, a, b , the optimization algorithm alternates between the following 2 steps:

s ^te P A ₊₁ ) ^T = Y ^T Y pf _B (3)

, , _T r , d/ dx df dJ _itT df „ _T „ , _{^ \}

where, Y = y,l , and— = = 2- lf P ^r (P f -a - y -b- 1) .

dw dw dx df dx

[0039] The above alternating algorithm consisting of steps 1 and 2 is guaranteed to converge.

In step 1, the update for a, b is the unique minimizer of the cost function with fixed w . Step 2 is a gradient descent method with variable w and fixed a, b, where the step size μ _η is found by

Armijo's rule for line search. Therefore, the cost function always decreases at each step. The cost function is lower bounded by zero. The above alternating algorithm can converge for all practical purposes. The parameters selected for line search are: initial step size μ ₀ = 2 and a = 0.1,/? = 0.5 as defined by Armijo's rule. These values generally lead to fast convergence and were used for all experiments. The algorithm stops whenever the step size for current iteration is sufficiently small (e.g., fixed at 0.1 for one aspect). The reason is that each PCA coefficients has been preconditioned so that their standard deviation is 10 in the training set; thus a step size of 0.1 can be considered sufficiently small and do not alter the results significantly. The algorithm is also terminated when the maximum number of iterations is reached (e.g., fixed at 10 for one aspect). In the case of largest deformation from the reference image, it can take around 8 or 9 iterations for the algorithm to 'converge,' in the sense that the cost function does not change much after that. Therefore, a maximum of 10 iterations should be sufficient to get convergent results.

[0040] At each iteration, given the updated PCA coefficients in Eq. (4), the DVF is updated according to Eq. (2) and the reconstructed image f _n+1 is obtained through trilinear interpolation.

Accordingly, df/dx has to be consistent with the interpolation process in order to get the correct gradient. Next, the expression for df/dx can be derived. Let f ₀ (i, j, k) and f(i, j, k) be the reference image and reconstructed image at iteration n indexed by the integer set

{ (i , j, k) 1 1 < < N ₁ , 1 < j≤ N ₂ ,l≤k≤N ₃ } within the bounding box of its volume, and X _j (z, j,k),x ₂ (i, j,k), x ₃ (z, j,k) be the deformation vector fields corresponding to the three canonical directions, namely, lateral, AP, and SI, respectively. They can be equivalently reorganized into a vector whose length is the total number of voxels in the image. Then di/dx = [di/dx ₁ ;di/dx ₂ ;di/dx ₃ ] .

f (z, j, k) = f ₀ (i + x^z, j, k), j + x ₂ (i, j,k), k + x ₃ (z, j, k)) (5)

Then:

This means that the Jacobian matrix df/dx _l is a diagonal matrix.

[0042] If the continuous patient geometry can be approximated by trilinear interpolation of the discrete volumetric image, then we have:

f(i,j,k) =f ₀ (l,m,n)(l-z ₁ )(l-z ₂ )(l-Z ₃ ) + f ₀ (l+l,m,n)z ₁ (l-z ₂ )(l-z ₃ )

+f ₀ (l,m + l,n)(l-z ₁ )z ₂ (l-z ₃ )+f ₀ (l,m,n + l)(l-z ₁ )(l-z ₂ )z " ₃ 3

+f ₀ (/ + l,m+l,n)z ₁ z ₂ (l-z ₃ )-l-f ₀ (I + l,m,n + l) z^l- z ₂ )z - ₃ 3 +f ₀ (/, m + l,n + l) (!-¾)¾¾ +f ₀ (l + l,m + l,n + l)z _l z ₂ z ₃

(7) where,

/ = |_z + Χ ₁ (Ϊ, j,k)\,m = \_ j + x ₂ (i, j,k)\,n = \_k + x ₃ (i, j,k)j are the integer parts of the DVF and z = {x _q (i, j,k)} ,<? = 1,2,3 are the fractional parts of the DVF.

[0043] Notice that:

dl/dx ₁ = dm/dx ₁ = dn/dx ₁ = 0 ^ and ^dz i/ ^dx i ⁼¹ ' 5 ^ζ ₂ ^χ ι ₍₉₎

where, for simplicity, the index (z, j, k) for x ₁ has been omitted.

From Eq. (7)-(9), the following can be obtained:

df (ζ ^' , j,k)/dx ₁ (ζ ^' , j, k) = [f ₀ (/ + 1, m, n ) - f ₀ (/, m, n )] (1 - z ₂ )(l - z ₃ )

+ | ^~ f ₀ (/ + l,m + l,«)-f ₀ (/,m + l,«)lz ₂ (l-z ₃ )

Γ / / Ί ⁽¹⁰⁾

+ f ₀ (/ + l,m,« + l)-f ₀ (/,m,« + l)J(l-z ₂ )z ₃

+ [f ₀ (/ + l,m + l,« + l)-f ₀ (/,m + l, « + !)]¾¾ Similar expressions can be derived for df/dx ₂ and df/dx ₃ .

[0044] Equation (10) is the exact derivative given the continuous patient geometry obtained from trilinear interpolation in Eq. (7). Of course, any reasonable kind of interpolation can be used here, e.g. , tricubic interpolation. Then the Jacobian matrix df/dx _l does not have a nice diagonal form anymore, but may be block diagonal. It is clear from Eq. (10) that di/dx is a linear combination of the spatial gradients of the reference image evaluated at the neighboring eight grid points, weighted by the appropriate fractional parts of the DVF. At this stage, dt/dx contains only local information at the level of individual voxels. It is the eigenvectors that effectively combine all the local information into global information and lead to correct update for the PCA coefficients.

[0045] For trilinear interpolation, f is continuous but not differentiable at the boundaries of the voxel grids, so that Eqs. (8)-(9) do not hold in this situation. However, this happens with a probability of 0. If it does happen, Eq. (10) simply computes the right-side derivative and will not create numerical instability. What's more, even if it happens that the DVFs are integers at a few voxels, the pooling effects of the eigenvectors will almost eliminates its influence by a vast majority of other voxels.

[0046] Selecting good initial conditions for the parameters may help speed up the convergence of the algorithm. In one aspect, for the PCA lung motion model based process 100 described with respect to FIGS. 1A and IB, the PCA coefficients are initialized to those obtained for the last projection. The results of PCA coefficients initialized to those obtained for the last projection are shown in FIGS. 3A, 3B, 4A, 4B, 4C, 4D, 4E, 4F and 5.

[0047] In another aspect, an algorithm that incorporates respiratory motion prediction (see FIGS. 2A and 2B) considers the last few projections and initializes the PCA coefficients to their predicted values, taking into account the respiratory dynamics. This might become important when the imaging frequency is low. There is a vast amount of literature on respiratory motion prediction models that can be used. In this document, as an example, a simple linear prediction model is used, which predicts the current sample using a linear combination of several previous samples, i.e. ,

where, w _k0 (t) is the initial guess for the k-th PCA coefficients for the current projection; w _k ^* (t - 1) is the optimized k-th PCA coefficients for previous /-th projections. Note that each

PC A coefficient has a prediction model associated with it since they have different dynamics. The respiratory motion prediction incorporating process and results are shown with respect to FIGS. 2A, 2B, 2C and 6A-11F.

[0048] The training set (i.e. , 4DCT) can be used to build the prediction model (alternatively, the first few breathing cycles during the CBCT scan may be used to build the prediction model). Because the sampling rate for 4DCT and CBCT scan is different, the PCA coefficients can first be interpolated from the training set and re- sampled to have the same frequency as in CBCT scan. For example, for a patient with a 4-sec breathing period, the sampling rate for 4DCT is 2.5 Hz if 10 phases were reconstructed. On the other hand, the sampling rate for CBCT projections is faster, usually on the order of 10 Hz for a typical 1-min scan with around 600 projections. Then the model fitting process can be performed separately for each of the PCA coefficients since they are assumed to be independent in the PCA motion model. The prediction model can be used to obtain a "good" initial guess for the PCA coefficients, so that the efficiency of the algorithm may be improved. For this purpose, the linear model used here is sufficient.

[0049] Tumor Localization via Deformation Inversion

[0050] To identify the current tumor position, it is important to distinguish between two different kinds of DVFs: push-forward DVF and pull-back DVF. The DVF found by Eq. (2) is a pull-back DVF, which is defined on the coordinate of the new image and describes how each voxel in the new image moves. The DVF cannot be used directly to calculate the tumor position in the new image unless the DVF is rigid everywhere. To identify the correct tumor position, the inverse, i.e. , the push-forward DVF should be obtained. The inverse can be defined on the coordinate of the reference image and can describe how each voxel in the reference image moves. In the respiratory motion prediction incorporating algorithm described in this document, an efficient fixed-point algorithm can be adapted for deformation inversion, which can be about 10 times faster and yet 10 times more accurate than other gradient-based iterative algorithms. For tumor motion, it is not necessary to apply the deformation inversion procedure on the entire image. Instead, the deformation inversion procedure can be applied on only one voxel, which corresponds to the tumor center of volume in the reference image.

[0051] In another aspect, the above described PCA lung motion model based algorithm can be modified to incorporate respiratory motion prediction. The accuracy and efficiency of the respiratory motion prediction incorporating algorithm for 3D tumor localization were evaluated on 1) a digital respiratory phantom, 2) a physical respiratory phantom, and 3) five lung cancer patients. Described in this document are these evaluation cases that include both regular and irregular breathing patterns that are different from the training dataset.

[0052] For the digital respiratory phantom with regular and irregular breathing, the average 3D tumor localization error is less than 1 mm which does not seem to be affected by amplitude change, period change, or baseline shift. On an NVIDIA Tesla C1060 GPU card, the average computation time for 3D tumor localization from each projection ranges between 0.19 and 0.26 seconds, for both regular and irregular breathing, which is about a 10% improvement over first algorithm described above. For the physical respiratory phantom, an average tumor localization error below 1 mm was achieved with an average computation time of 0.13 and 0.16 seconds on the same GPU card, for regular and irregular breathing, respectively. For the five lung cancer patients, the average tumor localization error is below 2 mm in both the axial and tangential directions. The average computation time on the same GPU card ranges between 0.26 and 0.34 seconds.

[0053] Through a comprehensive evaluation of the respiratory motion prediction

incorporating algorithm, the accuracy of this algorithm in 3D tumor localization is described to be on the order of 1 mm on average and 2 mm at 95 percentile for both digital and physical phantoms, and within 2 mm on average and 4 mm at 95 percentile for lung cancer patients. The results also indicate that the accuracy is not affected by the breathing pattern, be it regular or irregular. High computational efficiency can be achieved on GPU, requiring 0.1 - 0.3 s for each x-ray projection.

[0054] FIGS. 2A, 2B and 2C are flow charts that summarize an exemplary respiratory motion prediction incorporating process 200 for volumetric image reconstruction and tumor localization based on a single x-ray projection. A GPU can be implemented to perform the process 200 as described below. The process 200 is similar to the process 100 except for incorporating respiratory motion prediction. The process 200 can include using 4DCT from treatment simulation as a priori knowledge, deformable image registration can be performed between a reference phase and other phases, resulting in a set of deformation vector fields (DVFs) (210). A parametric PCA lung motion model can be obtained (220). The set of DVFs can be represented efficiently by a few eigenvectors and coefficients obtained from the PCA lung motion model (230). The GPU can perform image reconstruction using the obtained parametric PC A lung motion model (240). The GPU can identify a set of optimal PC A coefficients such that the projection of the corresponding volumetric image matches with the measured x-ray projection (242). In identifying the set of optimal PCA coefficients, the PCA coefficients are initialized to predicted values of last few projections and taking into account the respiratory dynamics (244). Then the PCA coefficients are varied until the computed projection image matches measured one of the new volumetric images (246). By varying the PCA coefficients, new DVFs can be generated, which, when applied on the reference image, can lead to new volumetric images. A volumetric image can be construed from a single projection image by optimizing the PCA coefficients such that its computed projection matches the measured one. The 3D location of the tumor can be derived by applying the inverted DVF on its position in the reference image (250).

[0055] FIG. 2C is another flow chart illustrating addition processes for performing image reconstruction using the obtained parametric PCA lung motion model that incorporates respiratory motion prediction. A GPU can load image data including the projection matrix (P) which computes the projection image of reference image f, measured projection image (y), and reference image (fo) (260). The GPU can set a step size μ=μ ₀ where the step size is found by Armijo's rule for line search (261). The GPU can initialize the variable w according to equation 11 (262). This takes account of the respiratory motion. The GPU can update a and b according to equation (3). The GPU can compute gradient J/w (264). The GPU can apply Armijo's rule to find μ (265). The GPU can determine whether the found value for μ is less than 0.1 (266). This is to determine whether the step size for the current iteration is sufficiently small (e.g., 0.1). Each PCA coefficient is preconditioned so that their standard deviation is 10 in the training set, and thus a step size of 0.1 can be considered sufficiently small and do not alter the results significantly. When determined that the step size is sufficiently small (266- Yes), the process terminates and the GPU can compute the inverse DVF (270). The inverse DVF can be used to generate the reconstructed image (271), and the 3D location of tumor can be derived by applying the inverted DVF on its position in the reference image (250).

[0056] When the GPU determines that the step size is not sufficiently small (266-No), the GPU can upgrade w according to equation 4. Then the GPU can update x and f (268). The GPU can also terminate the process when the maximum number of iteration is reached (e.g., 10) (269). When the maximum number of iteration has been reached (269- Yes), the GPU can compute the inverse DVF (270) and output the reconstructed image (271). The 3D location of tumor can be derived by applying the inverted DVF on its position in the reference image (250).

[0057] When the maximum number of iteration has not been reached (269-No) and the optimum coefficients have not been identified, the process can go back to update a and b according to equation (3) (263) and continue the process as shown in FIG. 2C.

[0058] GPU Implementation

[0059] To achieve real-time efficiency, the described algorithms can be implemented on graphic processing units (GPUs). Various programming environments and Graphics Cards can be used. For the examples provided in this document, compute unified device architecture (CUD A) was used as the programming environment and an NVIDIA Tesla CI 060 GPU card was used as the Graphics Card. Other examples of programming environment can include GPU programming languages such as OpenCL and other similar programming languages.

[0060] To evaluate the PCA lung motion model based algorithm and the respiratory motion prediction incorporating algorithm, the end of exhale (EOE) phase in 4DCT were used as the reference image and deformable image registration (DIR) was performed between the EOE phase and all other phases. The DIR algorithm used here is a fast demons algorithm

implemented on GPU. Then PCA was performed on the nine DVFs from DIR and three PCA coefficients and eigenvectors were kept in the PCA lung motion model (except for the physical phantom, where the motion is along the longitudinal axis so that only one PCA coefficient is necessary). For the prediction model in Eq. (11), L = 2 was selected in all the cases.

[0061] Digital Respiratory Phantom-PCA Lung Motion Model Based Algorithm

[0062] The PCA lung motion model based algorithm was tested using a non-uniform rational B-spline (NURBS) based cardiac-torso (NCAT) phantom. This mathematical phantom is based on data from the Visible Human Project, and is very flexible, maintaining a high level of anatomical realism (e.g., a beating heart, detailed bronchial trees, etc.). The respiratory motion was developed based on basic knowledge of respiratory mechanics. The NCAT phantom was used instead of real clinical data because with the NCAT phantom, the ground truth volumetric image and the 3D tumor location are available at any given time of the breathing. This is important for an objective assessment of the accuracy of our algorithm. A dynamic NCAT phantom composed of 10 breathing phases was generated as the training data, with a 3D tumor motion magnitude of 1.6 cm and a breathing period of 4 seconds. The dimension of each volumetric image is: 256x256x 120 (voxel size: 2x2x2.5 mm ). The end of exhale (EOE) phase was used as the reference image and deformable image registration (DIR) between the EOE phase and all other phases was performed. The DIR algorithm used here is a fast demons algorithm implemented on GPU. Then PCA was performed on the nine DVFs from DIR and three PCA coefficients and eigenvectors were kept in the PCA lung motion model. For testing purposes, a volumetric image was reconstructed and derived the 3D tumor location from each of 360 cone beam projections, which are generated from the NCAT phantom and uniformly distributed over one full gantry rotation. The phantom has a 4 s breathing period and a 50% increase in breathing amplitude (3D tumor motion magnitude: 2.4 cm) relative to the training data. The gantry rotation lasted one minute, resulting in 15 breathing cycles and 24 projections per cycle. Corresponding to 24 projections, 24 volumetric images at 24 breathing phases were obtained. These 24 volumetric images were used as ground truth test images to evaluate the accuracy of the reconstructed images. The imager has a physical size of 40x30 cm . For efficiency considerations, every measured projection image was down- sampled to a resolution of 200x 150 (pixel size: 2x2 mm ).

[0063] To quantify the accuracy of the reconstruction algorithm, the relative image error was defined as: e = . ( t _t - f- ) ² / y\ t _t ² , where/is the reconstructed volumetric image, /* is the ground truth image from the NCAT phantom, and the summation is taken over all the voxel indices. The localization accuracy can be quantified by the 3D root-mean-square (RMS) error.

[0064] The accuracy of the described algorithm can be limited by the accuracy of training DVFs derived with the demons algorithm. The potential accuracy that can be achieved is the accuracy of the DIR between the ground truth test image and the reference image using demons. Accordingly, the reference image is deformed to the ground truth test image using the same demons algorithm and the relative image error of the deformed image is computed against the ground truth test image. This error, termed as deformation error, in contrast to the reconstruction error of the reconstructed image, is used as the benchmark to evaluate our algorithm. Similarly, the 3D RMS deformation error can be computed.

[0065] PCA Lung Motion Model Based Algorithm-Data [0066] FIGS. 3 A, 3B, 4A, 4B, 4C, 4D, 4E, 4F, and 5 show the data for the PCA lung motion model based process. FIGS. 3A and 3B show the "measured" projection image, simulated from the NCAT phantom at the end of the inhale phase, so that it is maximally different from the reference image. Specifically, FIG. 3A shows a "measured" projection of the test image at a right posterior oblique (RPO) angle (300). FIG. 3B shows an objective value (310) and the relative image reconstruction error (320) at each iteration given the image in FIG. 3 A. Ten iterations were performed and further iterations were found to have little influence on the results. The relative image reconstruction error is initially 35% and approaches to the deformation error (10.9% compared with 8.3%) after 9 iterations. A relative image error below 10% is usually visually acceptable.

[0067] FIGS. 4A, 4B, 4C, 4D, 4E and 4F show the sagittal and coronal views of the absolute difference images between the reference image and three other images, namely the ground truth test image, the image reconstructed by the described algorithm, and the deformed reference image using demons. The 3D RMS tumor localization error is about 0.9 mm.

(anterior-posterior: 0.2 mm; lateral: 0.8 mm; superior-inferior: 0.4 mm). Specifically, FIGS. 4A, 4B and 4C show: sagittal view of the absolute difference image between: ground truth test and reference images (400), ground truth test image and image reconstructed using a single projection, (410), and ground truth test image and the deformed reference image using demons (420). Tumor is a round object near the center of the slice. FIGS. 4D, 4E and 4F show the coronal view (430, 440 and 450 respectively) of the same absolute difference image as in FIGS. 4A, 4Band 4C. Tumor is a round object in the right lung.

[0068] Overall, for the 360 cone beam projections, the average relative image

reconstruction error is 6.9% ± 2.4%. The result is comparable to that (8.6%) obtained using 30 cone beam projections without prior knowledge. The average 3D tumor localization error is 0.8 mm ± 0.5 mm and is not affected by projection angles (see FIG. 5). In comparison, the average relative image deformation error is 5.4% ± 2.2% and 3D RMS tumor localization error from deformation is 0.2 mm ± 0.1 mm. The periodic pattern in the relative image reconstruction error of the described algorithm is similar to the breathing pattern. This may be due to the accuracy of the reconstructed image being dependent on the difference between the current patient geometry and the reference one. So if the difference is large, then the image error will also be large, and vice versa. There does not appear to be such a pattern in the 3D localization error.

[0069] FIG. 5 shows on the top row: the relative image error between the ground truth test image and reference image (500) compared to the relative image error between an image reconstructed using the described algorithm and reference image (502) as a function of cone beam projection angle. On the bottom row, FIG. 5 shows the 3D localization error between the ground truth test image and reference image (510) compared to the 3D localization error between an image reconstructed using the described algorithm and reference image (512) as a function of cone beam projection angle.

[0070] The PCA coefficients are initialized to those at its previous frame considering the high-frequency image acquisition (6 Hz) relative to breathing. The image reconstruction and tumor localization for each projection can be achieved in less than 0.4 seconds using an

NVIDIA Tesla C1060 GPU card. Particularly, depending on the difference between the test image and reference image, the algorithm converges between 0.17 seconds and 0.35 seconds (average 0.24 seconds). Compared with an implementation on MATLAB 7.7 running on a PC with a Quad 2.67 GHz CPU, which takes around 15 minutes to converge, the GPU version can achieve a speedup factor of more than 2500.

[0071] Various implementations have been described to illustrate a GPU implemented algorithm to reconstruct volumetric images and localize 3D tumor positions from a single x-ray projection in a non-invasive, accurate, and efficient way. The relative image reconstruction error is only marginally larger than its lower bound defined by the deformation error. The average 3D RMS tumor localization error is below 1 mm. By utilizing the massive computing power of GPUs, a volumetric image can be reconstructed and 3D tumor locations can be derived from one projection within 0.4 seconds.

[0072] To investigate whether the projection image resolution of 200x 150 is sufficient, the algorithm was tested on a higher resolution of 1000x750 (pixel size: 0.4x0.4 mm ). The improvement of using the fine-resolution projection images is negligible (both image error and localization error are the same as above). In fact, further reduction in imager resolution does not have noticeable effects on the results unless it drops below 40x30, for the phantom case tested in this document.

[0073] As described herein, the dynamic NCAT phantom generated for testing purposes has a regular breathing pattern. Additionally, a dynamic NCAT phantom can be generated with irregular breathing patterns for testing purposes. For example, the accuracy of the algorithm can be tested on clinical data. Factors to consider for the algorithm can include: 1) the x-ray energies used to generate the training volumetric images (such as 4DCT) and cone beam projection images may be different; and 2) there will be scattering and other physical effects that may degrade the quality of the cone beam projection images in patient data. In both cases, a linear relationship between the image intensity of the computed and measured projection images may not be accurate. Some preprocessing (e.g. , a nonlinear transform) may be needed. The application of more general similarity measures such as mutual information can be implemented. While a few embodiments have been described using x-ray projections with rotational geometry, the same principle can be easily applied to those with fixed-angle geometry, such as fluoroscopy. The only difference is that the projection matrix will be constant for fixed-angle geometry instead of varying at each angle for the rotational geometry considered in this work. To further speed up the computation, the predicted PCA coefficients from previous histories can be used to serve as the initial guess for optimization. This procedure might be important when the imaging frequency is low.

[0074] Digital Respiratory Phantom - Respiratory Motion Prediction Incorporating Algorithm

[0075] In another aspect, the respiratory motion prediction incorporating algorithm was first tested using a non-uniform rational B-spline (NURBS) based cardiac-torso (NCAT) phantom. Three fundamental breathing parameters were considered: namely amplitude, period, and baseline shift relative to the training dataset. In the case of regular breathing, the same breathing parameters were first used as those in the training set and generated testing volumetric images. Then the breathing pattern was systematically changed such that each time, only one of the three parameters was changed while the other two were kept the same as those in the training set. In doing so, the effects of each of the three breathing parameters can be evaluated. The different breathing parameters considered include: amplitude (diaphragm motion: 1 cm and 3 cm), period (3 s and 5 s), and baseline shift (1 cm toward the EOE phase and EOI phase). In the case of irregular breathing, a scaled version of the breathing signal of a real patient acquired by the Realtime Position Management (RPM) system (Varian Medical Systems, Palo Alto, CA, USA) was used. The dynamic NCAT phantom motion can be generated by specifying the diaphragm SI motion. In this case, the scaled RPM signal was then used as a substitute for the diaphragm SI motion of the NCAT phantom. The breathing signal from RPM has a varying breathing amplitude (between 1.0 and 2.5 cm), period (between 3.5 and 6.2 s), and a baseline shift of about 1 cm.

[0076] After obtaining the dynamic NCAT phantom, 360 cone beam projections were simulated using the Siddon's algorithm from angles that are uniformly distributed over one full gantry rotation (i.e., the angular spacing between consecutive projections is 1°). For instance, the cone beam projection at angle 180° corresponds to the dynamic phantom at 30 sec, and the projection at angle 360° corresponds to the dynamic phantom at 60 sec, and so on. Since the gantry rotation lasts 60 seconds, the sampling rate of the cone beam projections is 6 Hz.

Quantum noise was not added in the simulated x-ray projections for the digital phantom cases. However, intrinsic quantum noise is present in all physical phantom and patient projections and no noise reduction was attempted on those images. The imager has a physical size of 40x30 cm . For efficiency considerations, every measured projection image was down-sampled to a resolution of 200x150 (pixel size: 2x2 mm ).

[0077] Physical Respiratory Phantom

[0078] The respiratory motion prediction incorporating algorithm was also tested on a simple physical respiratory phantom. The phantom consisted of a cork block resting on a platform that can be programmed to undergo translational motion in one dimension. This platform was used to simulate the SI respiratory motion of a lung tumor. Inside the cork block were embedded several tissue-like objects including a 2.5 cm water balloon which was used as the target for localization.

[0079] 4DCT of the physical phantom was acquired using a GE four-slice LightSpeed CT scanner (GE Medical Systems, Milwaukee, WI, USA) and the RPM system. During the 4DCT scan, the physical phantom moved according to a sinusoidal function along the SI direction, with a peak-to-peak amplitude of 1.0 cm. In order to calculate the correct projection matrix and image for the algorithm, the CT image corresponding to the EOE phase in 4DCT was re- sampled in space to have the same spatial resolution as CBCT and then rigidly registered to the CBCT. Since only rigid registration was performed at this stage, in reality, paired kV radiographs may also be used to achieve this task during patient setup. Cone-beam projections were acquired using Varian On-Board Imager 1.4 in full-fan mode with 100 kVp, 80 mA and 25 ms exposure time. The x-ray imager has a physical size of 40x30 cm and each projection image has a resolution of 1024x768. [0080] For testing purposes, two CBCT scans were performed for the physical phantom: one with regular breathing, and the other with irregular breathing. For the case of regular breathing, the phantom moved according to a sinusoidal function along the SI direction, with an increased peak-to-peak amplitude of 1.5 cm. For the case of irregular breathing, the phantom moved according to the same RPM signal used for the digital phantom as described before, with the peak-to-peak amplitude scaled to 1.5 cm. In both cases, 359 bone beam x-ray projections were acquired over an arc of around 200° with a frequency of about 10.7 Hz. The CBCT scans for the physical phantom lasted for about 36 s, so only the first 36 s of the RPM signal was used during the scan. Since the motion of the physical phantom during the CBCT scan is known, it was used as the ground truth to evaluate the accuracy of target localization.

[0081] The following preprocessing was performed on each of the cone beam x-ray projections. First, due to the presence of the full-fan bow-tie filter near the x-ray source, a systematic bias will be introduced in the linear relations between the image intensities of the simulated and measured projections. So an in-air x-ray projection with exactly the same imaging protocol was acquired in order to correct for the systematic bias. This was done by dividing each pixel intensity value of the projection images acquired for the physical phantom by those in the in-air projection image. Each corrected projection image was then down-sampled by a factor of 4, resulting in a resolution of 256x192 (pixel size: 1.56x1.56 mm ).

[0082] Patient data

[0083] The respiratory motion prediction incorporating algorithm was evaluated with five patient data sets. Both 4DCT and CBCT were acquired using the same imaging systems as for the physical phantom. The cone beam projections were taken with the imaging system in half- fan mode with 110 kVp, 20 mA and 20 ms exposure time. Similarly to the physical phantom, the CT at the EOE phase was re- sampled in space and then rigidly registered to the CBCT according to the bones. The same preprocessing on each of the raw cone beam projections of the patient was used as for the physical phantom, namely, cutoff, downsampling, and a logarithm transform, except for the bow-tie filter correction, which was based on the half-fan mode. For patient 2, however, since the computed projections of the 4DCT did not have sufficient longitudinal coverage compared with the cone beam projections, an additional 3 cm was cut off from the raw cone beam projections in both superior and inferior directions.

[0084] For all patient data there were no implanted fiducial makers and real-time 3D location of the tumor was not available to evaluate the respiratory motion prediction incorporating algorithm. Instead, the estimated 3D tumor location was projected onto the 2D imager and compared with that manually defined by a clinician. The patients chosen had tumors that were visible to the clinician in some of the cone beam projections. From the clinician-defined contour, the tumor centroid position was calculated for each projection and used as the ground truth to evaluate the respiratory motion prediction incorporating algorithm. The tumor localization error was calculated along the axial and tangential directions, both scaled back to the mean tumor position. The axial direction is defined to be along the axis of rotation on the imager and the tangential direction is perpendicular to the axis of rotation on the imager and tangential to the gantry rotation. For a tumor located near the isocenter, the tumor motion along the axial direction on the imager is roughly a scaled version of its SI motion in the patient coordinate system. On the other hand, the tumor motion along the tangential direction is a mixture of its AP and LR motions depending on the imaging angle. In some special cases, the tumor motion along the tangential direction can be roughly a scaled version of either the AP motion (for imaging angles near LR or RL) or LR motion (for imaging angles near AP or PA). Since the tumor is hard to see in some of the projections, the uncertainty in the definition of ground truth can be large, and may depend on the angle at which a projection is taken. A rough estimate of the uncertainty in the ground truth based on comparing contours drawn from two observers suggests that the error in centroid position is about 1-2 mm on average.

[0085] For each patient, approximately 650 projections were acquired over a full gantry rotation with a frequency of about 10.7 Hz. However, since the cone beam scans were performed in the half-fan mode (with the imager shifted about 14.8 cm laterally), and isocenter of the scan was not placed at the tumor, the tumor is only visible in a subset of these projections. For each patient, the tumor was marked by the clinician in the largest continuous set of projections in which the tumor was visible. For the five patients, the number of cone beam projections used for this study ranged between 70 and 281, corresponding to angles of 39° and 155°, respectively.

[0086] Summary of testing cases

[0087] Table 1 summarizes all the testing cases for the respiratory motion prediction incorporating algorithm, including 8 digital phantom cases, 2 physical phantom cases, and 5 lung cancer patient cases. The table lists the relevant breathing parameters that were used during the CBCT scan for tumor localization purposes. The table also includes information on the tumor size, location, as well as motion characteristics as measured in 4DCT for the 5 lung cancer patients. For the digital phantom, the breathing parameters in the training dataset are: 2 cm amplitude, 4 s period, and 0 cm baseline shift; for cases 2 through 7, all breathing parameters are the same as in training dataset except the one indicated in the table. All parameters for lung cancer patients were measured in 4DCT.

[0088] Table 2 summarizes the localization accuracy and computational time for all testing cases. For the description of each case, refer to Table 1. For Cases 1 through 10 (digital and physical phantoms), the localization errors were computed in 3D; for Cases 11 through 15 (lung cancer patients), the localization errors were computed on the axial and tangential directions with respect to the imager.

[0089] Results for the Digital Respiratory Phantom

[0090] FIGS. 6A, 6B, 6C, 6D, 6E, 6F, 6G, and 6H show the SI tumor position estimated by the respiratory motion prediction incorporating algorithm as well as the ground truth for the digital respiratory phantom for Cases 1 through 8 (600, 610, 620, 630, 640, 650, 660, and 670). The numerical results for localization accuracy and computational time are summarized in Table 2. The average 3D tumor localization error for regular breathing with different breathing parameters is less than 1 mm, which does not seem to be affected by amplitude change, period change, or baseline shift. The same holds true for irregular breathing.

[0091] From the re-sampled PCA coefficients in the training set, the prediction models were estimated as: w _lfl (0 = 1.90 (t - 1) - 0.98 (t - 2)

w _2f) (t) = 1.5 bv ₂ ^* (t - 1) - 0.84¼> ₂ ^* (i - 2) (12) w _3i0 (i ) = 0.97 w ₃ ^* (i - 1) - 0.84M¾* (t - 2)

[0092] By initializing the PCA coefficients to their predicted values using the above equations, the average computation time for both regular and irregular breathing ranges between 0.19 and 0.26 seconds using an NVIDIA Tesla C1060 GPU card. Particularly, the average computation time is 0.22 seconds in the case of regular breathing with a breathing amplitude of 3 cm. This is about 10% less than that using the PCA coefficients from the previous projection. In the case of irregular breathing, the average computation time is 0.26 + 0.06 seconds. The slight increase in computation time for irregular breathing was mainly due to a decreased accuracy of the prediction model used for parameter initialization.

[0093] Results for the physical respiratory phantom

[0094] FIGS. 7A and 7B show the estimated and ground truth tumor position in the SI direction for both regular (700) and irregular (710) breathing. The numerical results for localization accuracy and computational time are summarized in Table 2 (Cases 9 and 10). The average 3D tumor localization error for both regular and irregular breathing is less than 1 mm.

[0095] Since only one PCA coefficient was used for the physical phantom, the prediction model was estimated as:

w ₀ (t) = 1.96V (t - 1) - 0.99V (t - 2) (13)

[0096] When PCA coefficients were initialized to their predicted values, the average computation time on GPU is 0.13 seconds and 0.16 seconds, for regular and irregular breathing, respectively. The reduced computation time compared with the NCAT phantom is mainly due to the less number of PCA coefficients used for the physical phantom.

[0097] Results for the Patient Data

[0098] Results for the five lung cancer patients are shown in FIGS. 8A, 8B, 8C, 8D, 8E, 8F, 8G, 8H, 81 and 8J, where the tumor positions on the axial (800, 820, 840, 860 and 880) and tangential directions (810, 830, 850, 870 and 890) are shown separately. The solid circles represent the tumor positions estimated by the respiratory motion prediction incorporating algorithm and the ground truth is shown as the solid lines. The two arrows (a and b) in 800 indicate the two angles where the image reconstruction results are shown in FIGS. 9A, 9B, 9C, 9D, 9E, 9F. All five patients had somewhat irregular breathing during the CBCT scans especially for patient 2 (Case 12) and the last breathing cycle of patient 4 (Case 14). Again, the numerical results for localization accuracy and computational time are summarized in Table 2 (Cases 11 through 15). For all five patients, the average 3D tumor localization error is less than 2 mm and the absolute error at 95 percentile is less than 4 mm on both axial and tangential directions. When PC A coefficients were initialized to their predicted values, the average computation time on GPU is between 0.26 seconds and 0.34 seconds for the five patients.

[0099] FIGS. 9A, 9B, 9C, 9D, 9E, and 9F show the image reconstruction and tumor localization results for patient 1 at an EOE phase and an EOI phase (indicated in FIG. 8, Case 11), as well as the corresponding cone beam x-ray projections. Although it is hard to see the tumor on the projections, it is clearly visible in the coronal and sagittal slices of the reconstructed images. FIG. 9A shows the raw cone beam projection (900), and FIGS. 9C and 9E show the corresponding coronal (920) and sagittal (940) views of the reconstructed image at angle - 146.5° at an EOE phase (arrow "a" in FIG. 8). FIG. 9B shows the raw cone beam projection (910), and FIGS. 9D and 9F show the corresponding coronal (930) and sagittal (950) views of the reconstructed image at angle - 74.5° at an EOI phase (arrow "b" in FIG. 8). The arrows indicate where the tumor was. Solid lines represent the estimated tumor SI position at the current phase, and the dashes lines represent the estimated tumor SI position at the other phase.

[00100] There is a tradeoff between the localization accuracy and computational time, which is dependent on the resolution of the projection image. Intuitively, as the image resolution is increased, the localization error would decrease; at the same time, the resulting computational time will naturally increase. To investigate the effect of the projection image resolution, we took the original cone beam x-ray projection images with a resolution of 1024x768 for patient 1 (Case 11) and down-sample them by a factor of 1, 2, 4, 8, 16, 32 on both dimensions of the projection images. The resulting images have the resolution of 1024x768, 512x384, 256x192, 128x96, 64x48 and 32x24, respectively. FIG. 10 shows the average localization error on tangential and axial directions as well as computational time versus the downsamphng factor (1000). As expected, the average localization error increases monotonically as the downsamphng factor increases. However, downsamphng the projection image by a factor of 4 does not increase the localization error noticeably. On the other hand, the computation time decreases monotonically as the downsamphng factor increases, but flattens out once it reaches 8. This supports the choice of a resolution of 256x192 (downsamphng factor: 4) for the projection image used in this study, because it achieves a reasonable balance between the localization accuracy and computational time.

[00101] Applications

[00102] Various implementations have been described for the respiratory motion prediction incorporating algorithm to localize 3D tumor positions in near real-time from a single x-ray image. The respiratory motion prediction incorporating algorithm was systematically tested on digital and physical respiratory phantoms as well as lung cancer patients. For the digital phantom, the average 3D tumor localization error is below 1 mm in the case of digital and physical phantoms. For the five lung cancer patients, the average tumor localization error is below 2 mm in both the axial and tangential directions. By utilizing the massive computing power of GPUs, the 3D tumor locations can be derived from one projection around 0.3 seconds on average. Successful applications of the respiratory motion prediction incorporating algorithm in the clinic could lead to accurate 3D tumor tracking from a single x-ray imager. Furthermore, volumetric images can be reconstructed from a single x-ray projection image. This is potentially useful for many clinical applications such as dose verification in lung cancer radiotherapy.

[00103] The tumor localization error can arise from two different processes: 1), the PCA lung motion model; and 2), the 2D to 3D deformable registration between the projection and the reference image. In order to investigate how much each process contributes to the overall tumor localization error, the ground truth of the 3D lung motion were used in the digital NCAT phantom and solved for the optimal PCA coefficients by directly minimizing the mean square error between the left and right sides in Eq. (1). This procedure bypassed the 2D to 3D registration process, and thus gives the tumor localization error solely due to the PCA lung motion model. In Case 1, this error is about 0.28 mm over all the projections. Assuming that the above two processes are statistically independent, the error due to the 2D to 3D registration process can be obtained, which is about 0.75 mm. Therefore, the overall tumor localization error of about 0.8 mm may be mainly attributed to the 2D to 3D registration process. In comparison, the error introduced by the PCA lung motion model is insignificant.

[00104] For the respiratory motion prediction incorporating process, a linear relationship between the image intensities of the cone beam projection and of the computed projection of the reconstructed image is assumed, in order to account for the differences between the two image "modalities". There are several physical processes in reality (e.g., scattering) which are not simulated in the described computation of simulated x-ray projection. FIGS. 11 A, 11B, 11C, 11D, HE and 11F show scatter plots between the image intensities of the preprocessed cone beam projection and of the computed projection of the reconstructed image at two different projection angles for patient 1, corresponding to the two projections shown in FIGS. 9A, 9B, 9C, 9D, 9E and 9F. Specifically, FIGS. 11A and 11B show measured cone bean projections after preprocessing of the image intensities (logarithm transform followed by reversion of the sign) (1100 and 1110). FIGS. 11C and 11D show computed projections of the reconstructed image (1120 and 1130). FIGS. HE and 11F show scatter plots between the intensities of the images shown in FIGS. 11A, 11B, 11C and 1 ID (1140 and 1150). FIGS. 11A, 11C and HE correspond to the EOE phase in FIGS. 8A, 8B, 8C, 8D, 8E, 8F, 8G, 8H, 81 and 8J (arrow "a"). FIGS. 11B, 11D and 11F correspond to the EOI phase in FIGS. 8A, 8B, 8C, 8D, 8E, 8F, 8G, 8H, 81 and 8J (arrow "b"). A linear model between the two image intensities is able to explain 94% and 91% of the variance in the left and right subplot, respectively. Both linear models are significant (p<0.001).

[00105] A linear model is sufficient in both cases: more than 90% of the variance can be explained by a simple linear model. However, in the right subplot, because of more apparent interference from the treatment couch in the cone beam projection (which is absent in the computed projections), especially near the right edge, the linear model is less accurate than in the left subplot. This effect can be mitigated by pre- scanning the treatment couch and adding it to the reference image so that the couch will be also present in the computed projections.

[00106] While various implementations have bee described based on x-ray projections with rotational geometry. The same principle can be easily applied to those with fixed-angle geometry, such as fluoroscopy. The only difference is that the projection matrix will be constant for fixed-angle geometry instead of varying at each angle for the rotational geometry considered in this document.

[00107] It is possible that during the course of the treatment patients may undergo anatomical changes on an inter-fractional basis, which could make the results worse if the lung motion model is built based on 4DCT acquired during patient simulation. In this case, a PCA lung motion model built from 4DCBCT during patient setup may overcome this difficulty.

[00108] A respiratory motion prediction incorporating algorithm has been described for 3D tumor localization from a single x-ray projection image by incorporating respiratory motion prediction. Through a comprehensive evaluation of the respiratory motion prediction incorporating algorithm, the localization error can be identified to be less than 1 mm on average and around 2 mm at 95 percentile for both digital and physical phantoms, and within 2 mm on average and 4 mm at 95 percentile for five lung cancer patients on both axial and tangential directions. In addition, the localization accuracy is not affected by the breathing pattern, be it regular or irregular. The 3D localization can be achieved on a sub-second scale on GPU, requiring approximately 0.1 - 0.3 s for each x-ray projection.

[00109] FIG. 12 is a block diagram of an exemplary system 1200 for implementing real-time volumetric image reconstruction and 3D tumor localization for lung cancer radiotherapy. The system 1200 can be implemented as a computing system that includes at least a central processing unit (CPU) 1210, a graphics processing unit (GPU) 1220, an output unit 1230 and a memory/storage unit 1240. These system components can communicate with each other over an interconnect system such as a bus 1250. The GPU 1220 can perform the PCA lung motion model based process (see FIGS. 1A and IB) and the PCA lung motion model based process that incorporates respiratory motion prediction (see FIGS. 2A, 2B and 2C) for real-time volumetric image reconstruction and 3D tumor localization for lung cancer radiotherapy. In performing these processes for real-time volumetric image reconstruction and 3D tumor localization for lung cancer radiotherapy, the GPU 1220 can load the needed data from the memory/storage 1240 and output the results of the processes through the output unit 1230. Various programming environments and Graphics Cards can be used. For the examples provided in this document, compute unified device architecture (CUD A) was used as the programming environment and an NVIDIA Tesla C1060 GPU card was used as the Graphics Card. Other examples of

programming environment can include GPU programming languages such as OpenCL and other similar programming languages. The memory/storage unit 1240 can include various volatile and non-volatile memory devices including the various types of read only memories (ROMs), random access memories (RAMs), hard disks, Flash memories, etc. The output unit 1230 can include various output devices including liquid crystal displays, printers, storage devices, etc.

[00110] Implementations of the subject matter and the functional operations described in this specification can be implemented in various imaging and radiation therapy systems, digital electronic circuitry, or in computer software, firmware, or hardware, including the structures disclosed in this specification and their structural equivalents, or in combinations of one or more of them. Implementations of the subject matter described in this specification can be implemented as one or more computer program products, i.e., one or more modules of computer program instructions encoded on a tangible and non-transitory computer readable medium for execution by, or to control the operation of, data processing apparatus. The computer readable medium can be a machine-readable storage device, a machine-readable storage substrate, a memory device, a composition of matter effecting a machine-readable propagated signal, or a combination of one or more of them. The term "data processing apparatus" encompasses all apparatus, devices, and machines for processing data, including by way of example a

programmable processor, a computer, or multiple processors or computers. The apparatus can include, in addition to hardware, code that creates an execution environment for the computer program in question, e.g., code that constitutes processor firmware, a protocol stack, a database management system, an operating system, or a combination of one or more of them.

[00111] A computer program (also known as a program, software, software application, script, or code) can be written in any form of programming language, including compiled or interpreted languages, and it can be deployed in any form, including as a stand alone program or as a module, component, subroutine, or other unit suitable for use in a computing environment. A computer program does not necessarily correspond to a file in a file system. A program can be stored in a portion of a file that holds other programs or data (e.g., one or more scripts stored in a markup language document), in a single file dedicated to the program in question, or in multiple coordinated files (e.g., files that store one or more modules, sub programs, or portions of code). A computer program can be deployed to be executed on one computer or on multiple computers that are located at one site or distributed across multiple sites and interconnected by a communication network.

[00112] The processes and logic flows described in this specification can be performed by one or more programmable processors executing one or more computer programs to perform functions by operating on input data and generating output. The processes and logic flows can also be performed by, and apparatus can also be implemented as, special purpose logic circuitry, e.g., an FPGA (field programmable gate array) or an ASIC (application specific integrated circuit).

[00113] Processors suitable for the execution of a computer program include, by way of example, both general and special purpose microprocessors, and any one or more processors of any kind of digital computer. Generally, a processor will receive instructions and data from a read only memory or a random access memory or both. The essential elements of a computer are a processor for performing instructions and one or more memory devices for storing instructions and data. Generally, a computer will also include, or be operatively coupled to receive data from or transfer data to, or both, one or more mass storage devices for storing data, e.g., magnetic, magneto optical disks, or optical disks. However, a computer need not have such devices.

Computer readable media suitable for storing computer program instructions and data include all forms of non volatile memory, media and memory devices, including by way of example semiconductor memory devices, e.g., EPROM, EEPROM, and flash memory devices. The processor and the memory can be supplemented by, or incorporated in, special purpose logic circuitry.

[00114] While this specification contains many specifics, these should not be construed as limitations on the scope of any invention or of what may be claimed, but rather as descriptions of features that may be specific to particular embodiments of particular inventions. Certain features that are described in this specification in the context of separate embodiments can also be implemented in combination in a single embodiment. Conversely, various features that are described in the context of a single embodiment can also be implemented in multiple

embodiments separately or in any suitable subcombination. Moreover, although features may be described above as acting in certain combinations and even initially claimed as such, one or more features from a claimed combination can in some cases be excised from the combination, and the claimed combination may be directed to a subcombination or variation of a subcombination.

[00115] Similarly, while operations are depicted in the drawings in a particular order, this should not be understood as requiring that such operations be performed in the particular order shown or in sequential order, or that all illustrated operations be performed, to achieve desirable results. In certain circumstances, multitasking and parallel processing may be advantageous. Moreover, the separation of various system components in the embodiments described above should not be understood as requiring such separation in all embodiments.

[00116] Only a few implementations and examples are described and other implementations, enhancements and variations can be made based on what is described and illustrated in this application. For example, in addition to using the projection image, the described PCA lung motion model based algorithm and respiratory motion prediction incorporating algorithm can also use surface image, surface marker, implanted marker, diaphragm, lung volume, and air flow to reconstruct a volumetric image and to derive the 3D tumor position.