A Novel Markerless Lung Tumor-Tracking Method Using Treatment MV Beam Imaging

: A novel method was developed to track lung tumor motion in real time during radiation therapy with the purpose to allow target radiation dose escalation while simultaneously reducing the dose to sensitive structures, thereby increasing local control without increasing toxicity. This method analyzes beam’s eye view radiation therapy treatment megavoltage (MV) images with simulated digitally reconstructed radiographs (DRRs) as references. Instead of comparing global DRRs with projection images, this method incorporates a technique that divides the global composite DRR and the corresponding MV projection into sub-images called tiles. Registration is performed independently on tile pairs in order to reduce the effects of global discrepancies due to scattering or imaging modality differences. This algorithm was evaluated by phantom studies while simulated tumors were controlled to move with various patterns in a complex humanoid torso. Approximately 15,000 phantom MV images were acquired at nine gantry angles, with different tumors moving within ranges between 10 and 20 mm. Tumors were successfully identiﬁed on every projection with a total maximum/average error of 1.84/0.98 mm. This algorithm was also applied to over 5000 frames of MV projections acquired during radiation therapy of ﬁve lung cancer patients. This tumor-tracking methodology is capable of accurately locating lung tumors during treatment without implanting any internal ﬁducial markers nor delivering extra imaging radiation doses.


Introduction
The fundamental goal of radiotherapy is to deliver targeted tumoricidal radiation doses while minimizing doses to critical structures and healthy tissues. Large planning margins may be required to expand the target volume since respiratory motion can lead to lung tumor displacements of up to 2-3 cm [1,2]. The margins must be large enough to encompass the entire tumor motion range; if the margin is incorrectly determined, the lung tumor may move out of the planning target volume (PTV) during treatment. However, large margins produce larger PTVs, limiting the ability to deliver high doses to targets and subsequently compromising dose sparing to critical structures and normal tissues. Although motion ranges are normally quantified by four-dimensional (4D) computed tomography (CT) scans, 4D CT results represent the regular motion patterns and lung tumor placements only at the time of CT scan planning. Actual tumor motions and locations may change significantly during daily treatment [3,4]. Monitoring lung tumor motion during radiation therapy treatment will verify if tumors have received the prescribed radiation doses and if critical structures have been spared. It will also help to reduce the planning margins for incoming radiation therapy to spare healthy tissues.
Tumor motion management techniques, such as respiratory-gating and breath hold, along with lung tumor-tracking techniques, such as internal fiducial markers, external markers, and surrogates, have been used to monitor tumors [5]. However, several difficulties exist when using motion management techniques, including patient discomfort, patient-physician learning curve, and reproducibility. These ultimately determine the effectiveness and feasibility of their usage [5]. Fiducial markers may be implanted in and around the tumor site in order to track tumor motion. However, the nature of this technique poses several issues [3,4]. First, the patient faces a significant risk of pneumothorax and infection due to the invasive nature of the procedure. Second, over the course of the treatment, the markers may migrate so that the accuracy of the correlation between the tumor and the marker could be significantly reduced. Third, the implanted markers provide motion information of the center of mass, with no indication of the extent of tumor deformation and migration in response to radiation therapy over time. For this purpose, a non-invasive real-time lung tumor motion tracking algorithm is crucial for overcoming these issues. Attempts have been made to localize lung tumors on megavoltage (MV) projection images. Multiple research groups have conducted mask or template-based tumor tracking in the last decade [6][7][8]. Berbeco and his colleagues reported a series of studies on tumor tracking from a landmark-based algorithm that matches simulated images with MV projections [9][10][11][12][13][14][15][16]. Zhang et al. reported a kernel-based markerless approach for detecting lung tumor boundaries [17]. However, these methods may not be reliable if the tumor signal is week in comparison to the background. Poor soft tissue resolution may be due to the global discrepancies that exist between different imaging modalities, radiation scatter, and strong signals superimposed on the tumor by neighboring anatomical structures, among others. It is essential to accurately detect tumor motions when tumor projections and background projections overlap with each other. We report a novel but simple algorithm to generate tumor and background projections separately and determine their independent motions by comparing the estimated final projections with the acquired MV images. In order to reduce the effects of global discrepancies, a tiling technique is applied that accurately estimates image similarities [18].

Materials and Method
Continuous MV projection images were acquired in beam's eye view (BEV) by an electronic portal imaging device (EPID) during radiation therapy treatment. Figure 1 illustrates the geometry of the imaging system. A 2-dimensional (2D) axial CT slice represents a patient. Typically, a patient is positioned head first and facing up while lying on a treatment couch to locate the treatment target at the isocenter of a linear accelerator. Both an X-ray source and an EPID imager are mounted on a gantry which rotates around a horizontal axis corresponding to the patient's superior-inferior (SI) or head-foot direction. At any gantry angle (GA), a central radiation beam passes through the isocenter. The imager has two dimensions. As shown in Figure 1, at any gantry angle, one axis of the imager always coincides with the gantry rotation axis, the SI direction. The other orthogonal (Orth) component changes in space depending on the gantry angle. The Orth direction coincides with the patient's anterior-posterior (AP) direction at GA = 270 • or lateral direction at GA = 0 • . Normally, the most significant component of lung tumor motions occurs in the SI direction. However, depending on the location of the tumor within the thorax, the AP or lateral motion contributions may be more pronounced than the SI direction.
An overview of the markerless tiling technique flow chart is illustrated in Figure 2. Planning kilovoltage (kV) CT images was converted to MV CT image by introducing a kV-MVCT corresponding curve in order to match the digitally reconstructed radiographs (DRR) with MV treatment EPID images. The fitting corresponding curve was obtained by scanning the same CT electron density phantom with a kV CT scanner and an MV CT scanner (Tomotherapy, Accuray Inc., Sunnyvale, CA, USA).  During a conventional treatment planning process, tumors and surrounding critical anatomical structures were contoured by physicians on the basis of the maximum intensity projection (MIP) images from the 4D CT scan data. MIPs provide an overall moving envelope of tumors and other regions of interest (ROIs) but do not give an exact snapshot at a particular timestamp. Lung tumors were removed from the anatomy, and the cavities were filled with tissues with averaged lung tissue CT numbers. The complete anatomy was divided into two separate parts, i.e., the tumor and the other structures. Two sets of DRRs were generated from the converted MV CT images using a ray-tracing algorithm, which simply integrates the CT numbers along the X-ray trace from the X-ray source point to every imager pixel. For every MV projection image acquired, a pair of DRRs (background DRR consisting of the surrounding ROI with no tumor and tumor DRR consisting of only the tumor, as illustrated in Figure 3a,b) were generated. The same geometric configurations as those for acquiring MV projections after a pre-process were used to detect and remove the boundary edges (as illustrated in Figure 3c). A composite DRR was formed by a superimposition of the tumor DRR on a similar sized window of the background DRR with individual positions, as illustrated in Figure 3d. The tumor location on each MV projection frame was determined by the relative position of the tumor DRR fused on the background DRR that yielded the highest similarity between the composite DRR and the MV projection frame. (e) MV projection with 2 × 2 titling (two rows and two columns); (f) composite DRR with 2 × 2 tiling (two rows and two columns). All four tiles were normalized individually so that their average intensities were zero. The widow level of the upper left tile changed significantly, much closer to the corresponding tile in (e), and better similarity was observed for the upper left tiles using tiling-based matching than using global matching between (c,d).
Additionally, an image enhancement filter was used to improve the contrast ratio between the ROIs and the noise. A variation of the contrast limited adaptive histogram equalization (CLAHE) [19] filter was used to enhance poor soft tissue resolution and introduce high-gradient curves between the ROIs while suppressing the noise. The noise suppression was achieved using a contrast-limiting factor that allows the user to control the strength of the contrast by specifying the capacity of each bin. In general, CLAHE splits the global image into small contextual regions, performs histogram equalization on each region individually, and then stitches them back together. CLAHE performs well when the tumor signal is weak and when there are strong signals superimposed on it by neighboring anatomical structures.
A tiling strategy was introduced where the global MV projection image and the corresponding composite DRR were divided into sub-images called tiles [18]. The tiles of the MV projection image and the composite DRR were evaluated to find similarities between corresponding tiles. The tile pairs of the MV projection and composite DRR were shifted globally and evaluated independently of each other to maximize the local similarities. As shown in Figure 3e,f after 2 × 2 tiling (two rows and two columns), four tiles were normalized individually so that their average intensities were zero. The widow level of the upper left tile of DRR changed significantly, much closer to the corresponding MV image tile in Figure 3e. As a result, the similarity achieved was greater than the similarity between Figure 3c,d (without tiling). Normalized cross-correlation (NCC) was used to account for image similarity that exists between a pair of DRR and projection tiles: Projections are represented by P, and DRRs are represented by D. P and D were represented as matrices of tiles U × V (i.e., U columns and V rows tiles). The starting pixel (p u , p v ) of a projection tile P (u, v)(u,v) was determined by the tiling configuration. On the other hand, the starting pixel was not only determined by the tile configuration but also affected by the motions. P u,v and D u,v are the average intensities of the MV projection and DRR tiles (u, v), respectively. In order to maximize the effectiveness of the algorithm, the tile configuration (matrix of small tiles) needs to be carefully designed. Factors such as tumor location, size, and MV projection dimensions are considered when splitting the global image into tiles.
The global NCC between the MV projection and the composite DRR is the cumulative NCC computed between the individual tile pairs: The global NCC was maximized by adjusting the relative position of the tumor DRR with the background DRR window. The coordinates of the tumor DRR, fused onto a background that yields the highest similarity match between the global projection and the DRR, represents the tumor location on the projection.
Different tiling matrix have been tested, and it was found that the global NCC could be increased when more tiles are used. As an example, the global NCC could be increased from 0.66 to 0.79 when a 2 × 2 tiling is applied for identical geometric matching. However, when the tile size is too small, a mismatch of the tumor might occur. The smallest tile size should be larger than the tumor size. This study used a 2 × 2 tiling technique.
The proposed lung tumor-tracking algorithm was evaluated on a sophisticated respiratory torso phantom to determine its accuracy. A chest cavity phantom constructed with materials equivalent to natural bones and soft tissues was used. Tumors were simulated by two cylindrical bulks consisting of pliable bolus. Their volumes were 10.6 cm 3 and 6.9 cm 3 , simulating a large and a small lung tumor, respectively. Two phantoms were configured with the large tumor or the small tumor at different locations. A tungsten sphere with a diameter of 2 mm was attached to the tumors. The tumors were fixed at one end of a Styrofoam rod and positioned on a computer-controlled moveable platform, such that the tumor and Styrofoam rod were independently movable inside the chest cavity of the phantom. Five tumor motion patterns were programmed to control a three-dimensional (3D) motion platform. They were all 3D sinusoidal respiratory motions with maximum motion ranges between 10 mm and 20 mm. All phantom studies were performed on a Truebeam linear accelerator (Varian Medical Systems, Palo Alto, CA, USA), equipped with a MV imager (EPID) with pixel sizes of 0.4 mm. MV treatment images were acquired continuously by an EPID in cine mode. Varian iTool (Varian Medical Systems, Palo Alto, CA, USA) was used to receive streaming MV BEV images instantly.
Two 3D conformal radiation therapy plans were generated for corresponding phantoms/tumors. Each plan had nine treatment fields at different gantry angles. Each plan was delivered five times, while the tumor was driven by five different motion patterns. On average, there were about 170 images acquired for every treatment beam. Extra pre-processing was performed for the phantom study. This marker was detected by a pattern-matching algorithm in both CT images and MV images as ground truth of motions. The marker was removed from both images, and the holes were filled by expanding gradually from the edges of the holes. The markerless tumor tracking was then started on the modified images.
BEV images of five lung cancer patients were acquired when patients were treated on a SynergyS (Elekta, Stockholm, Sweden) by stereotactic body radiation therapy (SBRT). The tumor average position and the cumulated probability within various ranges were studied to quantify tumor motion. The cumulated probability within a range, in a statistical manner, is the percentage of time when the tumor motion is equal to, or smaller than, a given range [20]. For an institutional protocol, an internal target volume (ITV) was contoured to include the motion of the tumor based on the 4D CT image set. The ITV is defined by the gross tumor volume (GTV) plus the tumor motion range, which was 5 mm according to the 4D CT results. An additional margin of 5 mm was added in all directions over the ITV resulting in the PTV. Hence, the contour of the PTV was similar to that of the GTV with an additional 10 mm margin. The cumulated probability within a range of 10 mm represents the ratio of time when the tumor is within the planning margin. An optimal planning margin will be recommended in future studies.

Phantom Study Results
Phantom studies simulated two lung cancer patients treated with a five-fraction SBRT technique. About 15,000 images were acquired and analyzed in the phantom study. Markers were successfully detected and removed. Tumors were identified on every single image even in the presence of strong signals caused by the surrounding anatomical structures.

Patient Data Study Results
This algorithm was also retrospectively applied to the BEV images of five lung cancer patients receiving SBRT treatment. Table 1 lists the characteristics of the patients. There were 5040 frames of BEV images that were acquired and analyzed. Figure 5 shows the tumor trajectories in the SI direction for the same patient (#5) in the same treatment field (GA = 270 • ) during three fractions/days of treatment. It is noted that both the motion ranges and the average positions of the tumor vary significantly. More details are listed in Table A3 in the Appendix A.

Discussion and Conclusions
This novel lung tumor-tracking algorithm is capable of accurately detecting lung tumor motion during treatment delivery. All tumor tracking was based on the simulation/planning CT images, which typically have a voxel size of 1.1 by 1.1 by 3.0 mm. As a result, the uncertainties of DRRs could be estimated at about 3.4 mm. At the same time, MV images had a pixel size of 0.25 mm projecting back to the isocenter. Our phantom studies demonstrate that the maximum 2D errors were 1.8 mm, which is half of the sum of the image uncertainties. This means that the reported method reached the maximum detection accuracy.
This motion management technique does not require implanted fiducial markers, adds no additional imaging dose, and provides a daily assessment of the tumor motion pattern as a means to ensure that the tumor motion is within prior estimates and is covered by the treatment planning margin. This method is particularly useful for lung SBRT with abdominal compression as a means to guarantee reproducibility of the compression application. This motion-tracking algorithm offers significant potential for future real-time adaptive radiotherapy strategies. It will enable clinicians to reduce margins and safely deliver higher radiation doses by accounting for the effects of respiratory lung tumor motion during treatment.
Author Contributions: T.R. developed the software and analyzed the data. T.D.C. carried out the experiments and collected the data. M.C., X.J., W.L., and S.B. helped to develop the algorithm and analyze the data. W.M. initiated and supervised this study.
Funding: This project has been partially supported by a Varian research grant, an Elekta research grant, and NSF award CCF-1718994.