Monitoring tumor motion by real time 2D/3D registration during radiotherapy.

BACKGROUND AND
PURPOSE: In this paper, we investigate the possibility to use X-ray based real time 2D/3D registration for non-invasive tumor motion monitoring during radiotherapy.
MATERIALS AND METHODS: The 2D/3D registration scheme is implemented using general purpose computation on graphics hardware (GPGPU) programming techniques and several algorithmic refinements in the registration process. Validation is conducted off-line using a phantom and five clinical patient data sets. The registration is performed on a region of interest (ROI) centered around the planned target volume (PTV).
RESULTS: The phantom motion is measured with an rms error of 2.56 mm. For the patient data sets, a sinusoidal movement that clearly correlates to the breathing cycle is shown. Videos show a good match between X-ray and digitally reconstructed radiographs (DRR) displacement. Mean registration time is 0.5 s.
CONCLUSIONS: We have demonstrated that real-time organ motion monitoring using image based markerless registration is feasible.

Inter- and intrafractional organ motion in conformal radiation therapy is a major source of uncertainties in dose application. To eliminate interfractional motion, well developed methods for patient setup, including robotic patient tables and setup control using stereoscopic X-ray imaging [1] and cone-beam computed tomography (CBCT) [2,3] are used. Therefore, a number of methods for maintaining the registration of the target volume and the treatment beam were proposed [4]. While precise patient setup is an inevitable procedure to ensure treatment success, residual periodic and aperiodic intrafractional motion remain a problem. Several promising solutions were presented; these include (a) the use of X-ray imaging and implanted fiducial markers [5,6], (b) the tracking of passive electromagnetic transponders implanted close to the tumor tissue [7], (c) the correlation of surface landmarks with internal motion [8,9], and (d) the correlation of surface motion with lung motion models [10-12]. While these efforts document the need to compensate for intrafractional organ motion, most of these approaches are stricken with challenges such as the implantation of markers or transponders, which pose an additional burden in clinical routine and may hamper follow-up magnetic resonance imaging [13]. Respiratory models, on the other hand, may suffer from drift problems, an insufficient correlation of organ movement patterns and external surrogate markers, and aperiodic motion. For this reason, the possibility to utilize 2D/3D registration [14], a well-known method for static patient setup, of intrinsic features on X-ray images to artificially generated digitally rendered radiographs (DRR) was proposed for motion compensation [15,16]. Such a procedure, providing an estimate for all six degrees-of-freedom in rigid body motion, can be used to derive a motion model [11], or it may replace tracking methods based on markers. In [17,18], such setups are presented. The challenge, however, lies in an implementation of a 2D/3D registration routine that provides reliable robust position estimates at a sufficient update rate. In [11], an update rate of tens of seconds is reported, while the motion between those updates is predicted using a motion model. The authors in [15] report a registration speed in the range of 5–10 s.

Given the possibility of hypofractionated radiotherapy and the availability of modern linear accelerators with additional X-ray imaging systems, real-time 2D/3D registration with update rates up to 5 Hz may provide a possibility for non-invasive motion monitoring. In this paper, we present a proof of concept study which achieves registration update rates in the range of 2 Hz using standard equipment. The feasibility of the method is presented in a phantom study and on the retrospective evaluation of five clinical cases.

Materials and methods

Our approach for real-time tumor motion monitoring using 2D/3D registration works as follows: for each data set, fluoroscopic intrafractional X-ray images acquired at a frame rate of 5.4 Hz were available. The computed tomography (CT) or CBCT volumes were aligned to the first X-ray image and the translational and rotational parameters were recorded. These parameters were subsequently used as a starting point for the registration of the next X-ray image in an iterative manner. The registrations were performed using a region of interest (ROI) focused on the planned target volume (PTV) as provided from the planning data stored in DICOM-RT format. The translational and rotational parameters were recorded after each registration. From the rotational and translational parameters, the virtual displacement along cranial–caudal (CC), left–right (LR) and anterior–posterior (AP) directions of the centroid of the phantom or tumor is reconstructed.

2D/3D registration

Intensity DRR-based methods were the first and the most intuitive solution proposed to solve the 2D/3D registration problem [19]. In the first step, a simulated X-ray image – the DRR – is derived from a CT volume by simulating the attenuation of virtual X-rays. The DRR image is then compared to the X-ray image by means of a merit function. Rigid body motion of the tumor is then estimated by an iterative optimization which yields the transformation of the volume T.

The generation of DRRs is typically the most time consuming step. In this evaluation, a dedicated ray-casting algorithm implemented on the graphical processing unit (GPU) was used [20]. The bounding structures of the objects in the volume were refined using the cuberille algorithm [21].

Another important choice for the registration is the merit function. We used the Normalized Mutual Information metric, where a statistical dependence between the distributions of intensities is used to achieve registration [22]. This approach proved successful and robust in our registration problem in comparison with other metrics [23,24].

In our case, the z direction coincided with the imaging beam axis. Our registration scheme uses 2D projection geometry with one projection, therefore the movement along the z-axis cannot be resolved [25]. The translation along the z direction has been fixed and the registration scheme was performed using five degrees of freedom (translations along x and y directions and rotations).

Image data

Experiments were carried out with a phantom and with five clinical patient data sets. The data sets consist of X-ray images acquired during tumor irradiation for visual monitoring purposes, CBCT volume data acquired prior to each fraction and CT images taken for the purpose of treatment planning. No motion monitoring was performed during patient data acquisition. All planning CT and CBCT volumes were interpolated to an isotropic voxel size of 1.0 mm3; X-ray images were also interpolated to a pixel size of 1.0 mm2. Regions of interest (ROI) were defined on the X-ray images for two reasons: (a) to reduce the computation time needed to generate the DRR since only pixels of the ROI were actually rendered and (b) to focus the registration process on the PTV so that rigid motion is a valid assumption.

The kV X-ray as well as the CBCT images were obtained with an Elekta Synergy linear accelerator (LINAC) equipped with an electronic portal imaging device (EPID) allowing fluoroscopic X-ray acquisition and a CBCT system (XVI). The CBCT X-ray unit uses a Perkin Elmer XRD amorphous silicon detector with an active surface of 410 × 410 mm2 and 1024 × 1024 pixels resolution. The detector panel for the CBCT and kV X-ray imaging is located 536 mm from the axis of rotation. The source is located at a distance of 1536 mm from the imaging panel. In the case of CBCT imaging, the images were captured at a fixed frame rate of 2.7 Hz. During the 360° rotation the system acquired approximately 650 planar images which were used to make a full 3D image. The fluoroscopic X-ray images were captured during regular treatment with a frame rate of 5.4 Hz always in the AP direction; for each data set more than 100 X-ray images were available. The images are available on hard-disk for registration immediately after acquisition. The CT images were obtained by a Siemens Somatom Plus 4 Volume Zoom (Siemens AG, Erlangen, Germany) at 120 kVp and 156 mAs with intraslice resolution of 0.97 × 0.97 mm2 and 4 mm slice thickness.

Phantom data sets

The respiratory gating system AZ-733V (ANZAI Medical, Tokyo, Japan) was used in the initial part of the experiments. This phantom constitutes three spheres of same diameter and different materials embedded in a circular cylinder. The cylinder is fixed to an axle and translates in one direction. The translation movement mimics the respiratory signal. The phantom was positioned on the couch so that the isocenter of the gantry system was located in the center of the cylinder. The axis of translation was positioned along the axis of rotation of the LINAC system. The phantom was set to the respiration mode at 15 cycles/min for acquisition of the X-ray images. The movement of the phantom was stopped to acquire the CBCT image. The phantom motion was also recorded in six degrees-of-freedom by means of a tracking system (Polaris, Northern Digital Inc., Waterloo, Canada) for validation purposes. The experiment protocol is essentially identical to the one used for the patients described in the next section.

Patient data sets

The five patients suffering from non-small cell lung tumor or lung metastases were treated in a routine procedure at the department of radiotherapy. In this study, the planning CT, CBCT and intrafractional anterior–posterior X-ray images of each patient were used off-line.

CBCT images were acquired during the treatment of the patient prior to irradiation. The alignment of the CBCT to the X-rays is usually easy since these are generated with the same device and with the same patient position. On the other hand, the planning CT offers better image quality compared to the CBCT images, and it contains also the important structures from the treatment planning. For these reasons we used the planning CT volume for the registration with the patient X-ray data. In order to map the CT coordinate system into the coordinate system of the CBCT (thus having easy alignment between CT and X-rays) we used the CBCTs performing a 3D/3D registration using the normalized mutual information metric as implemented in AnalyzeAVW 10.0 (Biomedical Imaging Resource, Mayo Clinic, Rochester, MN). The initial matrix transformation Tinit that matched the DRR to the first X-ray was subsequently generated using the projection geometry of the LINAC system.

We decided to use the available structures from the CT planning in order to define the X-ray and DRR ROI. The PTV appears to be a reasonable ROI choice since it defines the area actually irradiated and is focused on the tumor position. Fig. 1 shows how the ROI was extracted from the contours. Fig. 2 shows one representative CT planning slice for each patient with the lung, CTV and PTV contours. Below each CT slice, the simplified projection of the contours to the X-ray image is shown. The coordinate system used during the registration is also shown in Fig. 2.

Fig. 1

Illustration of the method used to define the ROI for the clinical patient data sets. The contours of the planning CT from DICOM-RT data (a) were used to derive a 3D surface model of the organs, PTV and CTV (b, only the 3D model of the right lung is shown); the points that define the vertices of the 3D surface model were projected into the 2D X-ray imaging plane using the first transformation matrix Tinit (c). The convex hull that enclosed all points of the organs is shown in (d). The convex hull of the PTV was used as the ROI for the 2D/3D registration.

Fig. 2

The patient data sets used for off-line validation of the method. One representative slice of the CT planning data for each patient with the contours of the right lung (green), left lung (magenta), CTV (cyan) and PTV (red) is shown in the first row. The second row shows the initial X-rays with the contours projected by the help of the initial transformation matrix Tinit, the projected PTV is the ROI used for the registration of the X-ray to the CT planning.

Soft- and hardware

The program used for intensity-based 2D/3D registration was developed in C++ using gcc 4.3 under Ubuntu 9.10. The minimizer chosen for the optimization was the NEUWOA software proposed recently by Powell [26]. The GPU parts were programmed through OpenGL interface and Cg language. The program was run on a standard personal computer equipped with an Intel Core 2 Duo CPU of 3 GHz each and an NVIDIA GeForce 8800 GTS series graphics card with 640 MB RAM.

Results

Phantom results

Fig. 3 shows the rotational and translational parameters extracted from the registration of the successive X-ray images (Fig. 3a) and the reconstructed CC, LR and AP displacement of the centroid of the phantom cylinder (Fig. 3b). From Fig. 3a, it is clear that the pure translation along x axis followed the tracked phantom signal but with latency in the steepest part of the phantom signal leading to a root mean square error of 2.7 mm. The amplitude of the displacement of 20 mm was well reproduced. There were residual displacements along the y axis with an rms of 0.1 mm. Rotational movements with a rms of 0.5°, 0.3° and 0.2° for the rotation around x axis, y axis and z axis, respectively, were also recorded. The influence of the rotations is enhanced in the reconstructed centroid motion (Fig. 3b) where the reconstructed movement along CC direction showed more discrepancy to the phantom signal than the only translational parameter along x axis. The root mean square error between the true phantom displacement and the reconstructed displacement by registration was 2.6 mm. The motion along the LR and AP directions was quite low with root mean square values of 0.6 and 0.4 mm, respectively. The mean registration time for one X-ray image was 0.74 s.

Fig. 3

Recorded displacements of the phantom extracted by 2D/3D registration with 5 dofs. The first row (a) depicts the translations and rotations parameters obtained after each X-ray registration. (b) The reconstructed motion of the centroid of the cylinder along CC (blue line), LR (green line), AP (red line) directions. The black dotted line represents the movement of the phantom recorded by a tracking system. (c) A screenshot of video1.mpg (available online) showing checkerboard images of X-ray and corresponding registered DRR images acquired during phantom motion. Please note that the frame rate of the reconstructed videos is 5 Hz (about the acquisition rate of the X-ray imager), and not the actual speed of the registration (which is about 2 Hz).

The video available on-line shows reconstructed checkerboard images of the X-ray and corresponding registered DRR images and reveals the good match between the two image pairs all along the phantom motion.

Patient results

Fig. 4 presents the reconstructed motion of the tumor centroid along CC (blue line), LR (green line) and AP (red line) directions of patients 1, 2, 3 and 5. The diaphragm motion extracted via X-ray image by the help of edge detection algorithm was added to each plot (black dotted line). Table 1 gives the rms and maximum amplitude of the displacements. For the patient 4, the tumor was situated in the apex of the lung (see Fig. 2), where the motion is low, moreover the ROI depicted a low contrast. Therefore, the detected displacements by registration were below 1 mm in any direction and therefore, the tumor motion for this patient is not shown in Fig. 4 and Table 1.

Fig. 4

Reconstructed motion of the centroid of the tumor along CC (blue line), LR (green line), AP (red line) directions for patients 1, 2, 3 and 5. The diaphragm motion of each patient is also shown (except for patient 4) as a black dotted line. Note that the displacement scale is different for each plot.

Table 1

rms and maximum amplitude of the tumor motion along CC, LR and AP directions for all the patients. If the extracted displacement did not feature a sinusoidal like signal, the statistics are not relevant and N.A. (non-applicable) is written in the table.

	RMS			Amplitude (min–max)			Mean reg. time (s)
	CC (mm)	LR (mm)	AP (mm)	CC (mm)	LR (mm)	AP (mm)
Patient 1	2.1	0.7	0.4	9.2	3.2	1.8	0.4
Patient 2	1	0.5	0.3	4.0	2.2	1.1	0.6
Patient 3	0.4	0.7	0.2	1.9	2.8	1.0	0.5
Patient 5	7.6	N.A.	N.A.	33.5	N.A.	N.A.	0.4

From Fig. 4 it is clear that the reconstructed CC displacements were in phase with the diaphragm motion revealing a consistent extraction of the lung tumor movement by 2D/3D image registration. Except for patient 3 for which CC and LR movements were of the same amplitude (1.9 and 2.8 mm, respectively, Table 1), the motion was generally larger in the CC directions than LR or AP directions. The large CC displacements of patient 1 (9.2 mm) and patient 5 (33.5 mm) are linked to the position of the tumors being located near the diaphragm. The CC displacement of patient 1 and patient 5 exhibits secondary peaks with higher frequency. A Fourier analysis revealed a frequency of 1.35 Hz for patient 1 and 1 Hz for patient 5 which corresponds to the frequency of the beating heart. For these patients, the ROI encompassed the heart wall (Fig. 2), whose movement was therefore detected during the registration. A drift in the displacements along LR and AP direction is visible for the patient 5 for the first 40 X-ray images. This is not likely a real displacement but rather a misregistration.

Finally, the DRR and X-ray images were having a good match as shown in the video of the patient 1 PTV available on-line. The mean registration time for the registration of one X-ray was between 0.39 and 0.55 s and depends on the ROI size (Table 1, last column).

Discussion

Radiation therapy is a common and cost-effective technology for the treatment of malignant tumors in lung, either as a palliative or curative approach. The level of sophistication in dose escalation planning is already considerable [27-30], and it is expected that the broad introduction of multimodal imaging and deformable registration techniques [31,32] will push the frontier of dose planning further. Still, all efforts are in vain if one cannot monitor the motion of the target volume.

In this paper we propose a methodology for image-based tracking of internal organs in IGRT, and we prove that the concept is feasible and exhibits limitations. Our first experiment deals with a simple unidirectional movement produced by a respiratory phantom. The results show that using an adequate ROI our 2D/3D registration scheme is able to follow the motion of the phantom with an rms error of 2.6 mm. Regarding the patient study, we first propose an efficient way to select an appropriate ROI for organ motion monitoring. The projection of the PTV seems to be an appropriate choice since it focuses on the area where the tumor is present and which is actually irradiated. The PTV projection is retrieved from the contour of the CT planning provided that the appropriate projection matrix is available.

Without a known ground truth, it is difficult to assess the accuracy of the registration on the real patient data set, however tumor motion that clearly correlates to the breathing of the patient is extracted. Moreover, our results are consistent with the results reported by Seppenwoolde et al. [33] and Suh et al. [34]: most of the patients have a predominant movement in the cranio-caudal direction and the motion amplitude depends on the tumor site. Virtually no motion is detected for a tumor in the apex of lung in patient 4. Large motion is detected for tumors in the lower lobes (patient 1 and patient 5).

Stereotactic body radiotherapy is done with 6–8 different gantry angles, however the registration will not perform uniformly in all directions as accuracy directly depends on soft tissue contrast in the X-ray images of the lung. Therefore, in this study all fluoroscopic X-rays were acquired in a gross AP direction, during lateral treatment beam placement, where best tissue differentiation was achieved. Since the kV imaging system in our setup has always a fixed 90° angle in relation to the treatment beam, this work is, for the time being, valid only for this beam position. A lower tissue differentiation with imaging from other directions can be expected but should nevertheless be sufficient for robust registration. In the future we aim at evaluating the system in all treatment directions as well as in other non-optimal gantry angles.

One of the major limitations of our approach is the insufficient information on PTV motion in the direction of the normal beam of the imaging device leading to misregistration in the AP direction. Also, rotations around LR and CC axes induce translations in the AP direction as the center of rotation in not located at the centroid of the tumor. Moreover, and with the small ROI, rotations around x-axis result in apparent translation in the y-axis and vice versa. This was the likely cause for misregistration in some cases (e.g. the drifts observed in patient 5). Another limitation is the inability to distinguish between movements due to respiration or heartbeat. However, since the major component of lung motion is the cranio-caudal direction, it is clear that in the treatment of lung tumors, we will not face the problem of missing the most important part of the tumor trajectory. Another limitation is the additional dose to the patient from permanent X-ray imaging. This problem can only be solved if further clinical studies will show that a significant reduction of the PTV can be achieved by permanent motion monitoring.

With our implementation the registration speed was as fast as 2 Hz with off-line images. For an on-line application, communication latency between the software and the Elekta system is inevitable and must be measured and added to the time necessary for the registrations themselves. The on-line real-time 2D/3D registration scheme, however, should be performed at a speed of 5.4 Hz, which is the acquisition rate of the X-ray imagers.

Our procedure can be improved by different ways. First, non-rigid registration algorithms can be used for the initial mapping of the CT to the CBCT volume for daily basis adaptation. Then, in order to overcome the problem of unresolved displacement along the beam axis, the volume could be registered to several X-ray images taken with different angles. However this would result in an increase of the computation time [18]. At the moment, only the rendering is performed using GPGPU programming techniques, the speed of the whole procedure can therefore be further improved by implementing the whole registration procedure using dedicated GPGPU systems. Finally, to fully validate our approach, a larger clinical study should be undertaken. The ground truth necessary in such a validation, could be recorded by following the movement of implanted markers and comparing it to the motion retrieved by 2D/3D registrations.

Conclusions

In this paper, we presented a feasibility study for the monitoring of real-time organ motion during radiotherapy based on image content using 2D/3D registration. The motion of a respiratory phantom and the tumor motion of patient were retrieved at a speed rate of 2 Hz. Even if a more detailed clinical study is needed in order to further validate our approach, we believe this first study lay the foundations for non-invasive real time organ motion tracking.

Conflict of interest statement

Any actual or potential conflict of interest does not exist.