Abstract
This study aims to employ 4-dimensional computed tomography to quantify intrafractional tumor motion for patients with lung cancer to improve target localization in radiation therapy. A multistage regional deformable registration was implemented to calculate the excursion of gross tumor volume (GTV) during a breathing cycle. GTV was initially delineated on 0% phase of 4-dimensional computed tomography manually, and a subregion with 20 mm margin supplemented to GTV was generated with Eclipse treatment planning system (Varian Medical Systems, Palo Alto, California). The structures, together with the 4-dimensional computed tomography set, were exported into an in-house software, with which a 3-stage B-spline deformable registration was carried out to map the subregion and warp GTV contour to other breathing phases. The center of mass of the GTV was computed using the contours, and the tumor motion was appraised as the excursion of the center of mass between 0% phase and other phases. Application of the algorithm to the 10 patients showed that clinically satisfactory outcomes were achievable with a spatial accuracy around 2 mm for GTV contour propagation between adjacent phases and 3 mm between opposite phases. The tumor excursion was determined in the vast range of 1 mm through 1.6 cm, depending on the tumor location and tumor size. Compared to the traditional whole image-based registration, the regional method was found computationally a factor of 5 more efficient. The proposed technique has demonstrated its capability in extracting thoracic tumor motion and should find its application in 4-dimensional radiation therapy in the future to maximally utilize the available spatial–temporal information.
Introduction
For radiotherapeutic management of thoracic cancer, one of the major challenges is the precise localization of the mobile tumor. Treatment outcome could be seriously compromised if the motion is not properly taken into consideration in planning and delivery. There have been a number of techniques proposed to address the issues. Some of them are (1) respiratory gating technique, 1,2 (2) incorporating respiratory motion into treatment volume, 3 -5 (3) breath holding technique, 6,7 and (4) tumor tracking strategy. 8 -14 Although there are pros and cons with each technique, they were developed to better control breathing motion. Four-dimensional computed tomography (4DCT) dividing a patient’s breathing cycle into multiple phases to better represent anatomy has become available. It provides spatial–temporal information of patient’s anatomy that can be utilized for better breathing control. As opposed to conventional treatment techniques where treatment planning and delivery were conducted in a static manner, 4D radiation therapy 15 -17 would be adopted as a standard strategy in the future. Under this scenario, knowing precisely intrafractional anatomy change becomes important. For instance, the definite motion pattern is required to bridge different breathing phases to develop a 4D plan. It is equitably essential that the automatic delineation of gross tumor volume (GTV) and organs at risk of 4DCT requires accurate information of the motion.
To estimate tumor motion based on the 4DCT sets, especially in the thoracic region, deformable image registration is a common choice. There are many deformable registration techniques. 18 However, most of them are whole image-based models. 19 -22 In this study, we are particularly interested in the tumor target motion; thus, there is no need to conduct the calculation on the entire image region. In addition, the registration accuracy may be compromised because the registration can be influenced unnecessarily by the image content distant from the region of interest, which would otherwise be irrelevant to the process. The tumor motion can be well represented as the location change in the GTV or the GTV excursion relative to a preselected breathing phase, for example, 0% phase of 4DCT. This can be readily computed with GTV contours. In order to attain GTV contours on all breathing phases, a straightforward solution is to delineate the tumor on a selected phase and then propagate the GTV contour onto the remaining phases through a deformable registration. In this study, we implemented a multistage regional B-spline deformable registration to conduct the 4DCT registration, aiming at extracting tumor motion with high efficiency and accuracy.
Methods and Materials
Four-dimensional computed tomography of the 10 patients with lung cancer was acquired and sorted into 10 respiratory phase bins. Figure 1 depicts the flowchart of the calculation algorithm. The GTV was manually delineated on 0% phase termed as template phase of 4DCT, and a subregion with 20 mm margin supplemented to the GTV as shown in Figure 2 was created on the Eclipse treatment planning system (Varian Medical Systems, Palo Alto, California ). Radiation therapy structures (GTV and the subregion), together with 4DCT, were exported into an in-house research software. A 3-stage B-spline deformable image registration was implemented to map the subregion to other breathing phases (coined as target phases). The 3-stage registration was introduced to further increase the calculation accuracy. The registration results from stage 1 were input to stage 2 and subsequently those out of stage 2 were fed to stage 3. At stage 1, the initial image set was downsampled to a quarter of the original size on the transverse plane and by 2 along the superior–inferior (SI) direction. At stage 2, images on transverse plane were downsampled by 2 but no downsampling along the SI direction. At the last stage, the original image set was utilized in the calculation. This coarse-to-fine scenario was more efficient in calculation and better in controlling local minimum compared to the single-stage model.

Flowchart of the calculation algorithm.

Schematic of creation of the subregion that is formed on gross tumor volume (GTV) with a 20 mm margin. Top: axial view; bottom: sagittal view.
The Mattes mutual information (MMI)
23
was used as the metric function, which was optimized by a limited memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) algorithm optimizer with respect to the displacement parameters of the nodes. The central concept of mutual information (MI) is the calculation of entropy. For an image A, the entropy is defined as
where
where
The MI measures the level of information that a random variable (eg, gray value of image set a:
The L-BFGS algorithm
25
was used to optimize the MMI metric function with respect to the displacement parameters of the nodes
26
and to find the transformation matrix Determine the descent direction line search with a step size update compute
At each iteration, a backtracking line search is used in L-BFGS to determine the step size of movement to reach the minimum of f along the ray
or a preset maximum number of iterations is reached. In this study, we set
The registration was carried out between the template phase and the remaining 9 target phases. Upon completion of the successful registration, 9 displacement vector fields were obtained and employed to warp GTV contours to the target phases. Using the resultant contours, we proceed to calculate the tumor excursion.
Figure 3 depicts the schematic of calculating GTV excursion. In digital image processing, contours are often represented as polygons on individual CT slices. As shown in Figure 3, polygon ABCDE represents the manual contour on the template phase and polygon A′B′C′D′E′ is the warped contour on one of the target phases. The centers of mass (COMs) O and O′ can be computed as

Schematic describing the calculation of the mean and maximum excursions between the template phase and target phase. GTV indicates gross tumor volume.
Ten patients with lung cancer treated with either stereotactic body radiation therapy or intensity-modulated radiation therapy, named as patient 1 to 10, respectively, were enrolled in the study. Table 1 shows the tumor target information including locations and sizes, together with the patients’ demographic data. Four-dimensional computed tomography images were acquired with a Philips Brilliance 16-slice CT scanner (Philips Medical System, Andover, Massachusetts). The collected data were sorted into 10 phase bins. The GTV on the template phase was manually segmented with the Eclipse treatment planning system. The GTV contour is shown in blue curve together with the expanded subregion (uncolored area) on the axial and sagittal views in Figure 2 top and bottom, respectively.
List of patients enrolled in the study and the tumor target information.
Abbreviations: yo, years old.
The proposed registration algorithm was implemented using the Insight Toolkit 27 and the Visualization Toolkit, 28 which are open source cross-platform C++ software toolkits sponsored by the National Library of Medicine. The image and contour display was based on the open source software VV, the 4D Slicer (http://www.creatis.insa-lyon.fr/rio/vv/).
Results
The proposed 3-stage regional deformable registration was performed between the template phase and 9 target phases of 4DCT sets for the enrolled patients. An example of the warped GTV contour overlaid on 50% phase CT from patient 2 is demonstrated in Figure 4. It is displayed in axial (A), sagittal (B), and coronal (C) views. Visual inspection of the contour shows clearly the excellent delineation of the GTV, indicating that the regional registration yielded reasonable accuracy. Going through warped contours on all patients, we observed that normally around 2 mm accuracy could be achieved for GTV contour propagation between adjacent phases and 3 mm between opposite phases. Manual delineation of GTV on the 9 target phases for the 10 patients was also completed as a reference contour to more quantitatively assess the accuracy of the COM excursion. The results are shown in Figure 5. The residual is defined as the COM difference between the warped GTV contour and the manual contour. For all patients, the largest residual is 1.89 mm, from phase 5 of patient 4, whereas the smallest is 0.59 mm, from phase 7 of patient 3, showing that our calculation accuracy is 2 mm.

Superimposition of the warped gross tumor volumes (GTV) contour on the CT from 50% phase for patient 2. They are displayed in (A) axial view, (B) sagittal view, and (C) coronal view.

Residual, calculated as the difference between the centers of mass (COMs) from manual contouring and calculation using our method, for all 9 target phases for all 10 patients.
More quantitative evaluation of the tumor excursion during the breathing cycle is demonstrated for all patients in Figure 6. Of them, the bulkiest tumor and the largest displacement take place in patient 1. The mean excursion spans from 2.4 to 5.7 mm, whereas the maximum range is between 7.7 and 15.1 mm. Both maximum and average target motions follow a similar pattern over the breathing period. The scatterplots also indicate that the tumor target moves in a relatively symmetric fashion. Although patient 1 represents an extreme case, patients 3, 5, and 10 experienced minimum movements, with mean motion less than 3.5 mm and maximum smaller than 4.0 mm, reflecting the fact that they are located on the upper lobes where breathing movement tends to be smaller. Motion patterns for patients 4, 6, and 9 exemplify relatively larger motion across various respiratory phases, as indicated by their attachment or next to the chest wall and diaphragm. It is of clinical interest to observe the largest excursion the tumor target could experience over the entire breathing cycle. Among them, the largest excursion is as great as 1.5 cm from patient 1 as discussed earlier. Closer inspection reveals that the motion was largely along the superior–inferior direction as demonstrated in the literature. 29,30 Sitting on the apex of the lung and with a relatively smaller size (diameter ∼1.6 cm), the least motion was 3.5 mm from patient 3.

The mean and maximum tumor excursion relative to 0% phase for the 10 enrolled patients. On the left hand side is the mean excursion and the right hand side is the maximum excursion for the target phases. The horizontal axis represents 9 breathing phases from 10% to 90% corresponding to numbers 2 to 10, respectively.
The margin size was varied for sensitivity analysis. It is found that if the size is less than 5 mm, then the calculation started to give large errors owing to insufficient image contents or simply the image registration failed. On the other hand, if the margin is chosen to be larger than 20 mm, no further improvement was observed in terms of calculation accuracy, but more computing resources were needed. Therefore, we chose to reduce the margin size down to 15 and 10 mm, respectively. It is found that if the tumor size is smaller than 2 cm, the registration had hardship in converging because insufficient surrounding contents were included even with 15 mm margin, for instance, patients 3 and 10. The analysis outcomes are demonstrated in Figure 7 for patients 1, 4, and 5 and Figure 8 for patients 6, 7, and 9. Figure 7A and B represent patient 1, Figure 7C and D represent patient 4, and Figure 7E and F represent patient 5. The left hand column plots are the mean excursion and the right hand column are the maximum excursion for 3 situations: 20, 15, and 10 mm margins, respectively. No significant difference was observed for patient 1 in terms of the margin choices; however, for patients 4 and 5, selecting different margins had significant impact on the motion calculation. For patient 5, using the 10 mm margin even led to wrong estimation of the motion, that is, the large digressions as shown in 80% phase and 90% phase in Figure 7F. Similar study outcomes are shown in Figure 8, where Figure 8A and B represent patient 6, Figure 8C and D represent patient 7, and Figure 8E and F represent patient 9. Significant difference was observed in patient 9. The large inconsistencies among different margins might be due to the relatively smaller tumor size (diameter ∼2.3 cm).

The study of varying the supplementary margins to the gross tumor volume (GTV) and comparison to the standard 20 mm margin for patients 1, 4, and 5. The legend shown in (A) is valid for (B) to (F) as well. The horizontal axis represents breathing phases from 10% to 90% corresponding to numbers 2 to 10, respectively.

The study of varying the supplementary margins to the gross tumor volume (GTV) and comparison to the standard 20 mm margin for patients 6, 7, and 9. The legend shown in (A) is valid for (B) to (F) as well. The horizontal axis represents breathing phases from 10% to 90% corresponding to numbers 2 to 10, respectively.
In terms of computing efficiency, a timer was set to count the time consumption for the regional technique as opposed to the whole image registration. For all patients, the average time for the whole image calculation required about 5 minutes to accomplish the 3-stage registration process, whereas the subregional technique only needed less than a minute to achieve similar accuracy, clearly showing the high efficiency.
Discussion
Tumor target localization is critical in administrating radiation dosage to patients with cancer. For patients with thoracic cancer, this entails precise respiratory motion handling. 31 -34 With spatiotemporal information revealed by 4DCT, one could extract tumor motion with various techniques. 12 -14,16,29,30 In this study, we calculated tumor excursion using a multiple stage regional deformable registration method. This subregional approach relies on the assumption that a subregion surrounding the manually segmented GTV contour can capture sufficient information to drive the finding of its counterparts on other phases of the 4DCT while unnecessary distant image contents have diminutive undesirable contribution to the calculation. This assumption is valid only when the region is sufficiently large so that a great number of voxels are involved in the registration calculation. As illustrated in Figures 6 to 8, the registration becomes reliable when the margin added to GTV is chosen no smaller than 15 mm, obviously depending on the tumor locations and sizes. When the tumor size is less than 2 cm and if there is not enough anatomy surrounding it as in patients 3 and 10, ample margin larger than 15 mm is necessary for successful registration. However, in the cases where tumor is small but it is immediately next to or attached to anatomy with good image contrast, the calculation could still be accomplished with a relatively smaller margin as shown in patients 4, 5, and 9. Computationally, the proposed approach is more economical and efficient. A comparison study was conducted with a timer using the whole image sets and revealed that the calculation time increased by a factor of 5 compared to our regional method, which makes it more attractive in eliminating the need for a global registration of the input images.
As illustrated in this work and others elsewhere, 12,13,16,29 the estimation of tumor motion using 4DCT is solely dependent on the accuracy of the deformable registration. A well-known problem in image registration and contour segmentation studies is the lack of quantitative validation. In the studies of Lu et al, 21,23,35 for examples, the accuracy of a deformable model-based contour mapping technique was evaluated purely based on visual inspection. Although it is a convenient way for rapid assessment of a registration calculation, especially in a case where the “ground truth” does not exist, the method falls short in quantization. Another way of assessing the registration accuracy is to perform a synthetic image experiment. 26,29 The idea is to warp one of the input image sets with known deformation vector fields that serve as the “ground truth.” Deformable image registration is then carried out between the synthetic image sets, and the calculation outcome is compared to the known “ground truth” so that a quantitative evaluation can be achieved. Another method of objectively evaluating registration accuracy was to utilize a physical deformable phantom 22,36 -38 where the deformation can be preprogrammed. Regardless of these efforts, it is still a challenging task to precisely assess the accuracy of a deformable image registration model in reality owing to so many complications and issues. Active research is still ongoing to tackle the problem.
It is worthwhile to point out that tumor motion extraction is based wholly on the GTV contours on 4DCT and it only gives the motion of the tumor as a single point, neglecting its volumetric changes. Although this is a simplification of the realistic question, it has useful application in 4D radiation therapy. For example, in the work of Suh et al, 39 a 4D plan was created by copying the multileaf collimator (MLC) from a template phase to target phases using the single point of tumor target motion. This planning strategy came up with a comparable dosimetric accuracy compared to the full 4D planning on each phase while making the MLC delivery possible. However, in reality, the tumor target in patients with lung cancer is a 4D object that has both positional change and shape variations over time. Therefore, a more realistic motion modeling has to take that into consideration, although it is quite challenging. 40 Another practical strategy is to implant fiducial marker into tumor target and using fluoroscopic images to track the tumor motion. 41,42 This is an invasive procedure and may not be applicable to every patient. Fluoroscopic images are often acquired before or during treatment, and they normally provide 1-dimensional information of the tumor target motion. At any rate, a more precise description of the tumor motion would require a 3-dimensional volumetric tracking. This would be a research topic for the future.
Conclusions
We have developed a multistage regional deformable registration method to calculate the tumor motion using 4DCT. The main idea is that a subregion encompassing GTV carries sufficient neighborhood information to establish a reliable association between 2 phase-specific image sets. A 2 mm accuracy can be achieved with the proposed strategy. Compared to the whole image-based deformable registration, a great reduction in computational burden and efficiency enhancement are observed, making it a useful tool for better control of respiratory motion in lung cancer radiotherapy.
Footnotes
Declaration of Conflicting Interests
The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The authors received no financial support for the research, authorship, and/or publication of this article.
