- Email: [email protected]
Regression model-based real-time markerless tumor tracking with fluoroscopic images for hepatocellular carcinoma
Physica Medica 70 (2020) 196–205
Contents lists available at ScienceDirect
Physica Medica journal homepage: www.elsevier.com/locate/ejmp
Original paper
Regression model-based real-time markerless tumor tracking with fluoroscopic images for hepatocellular carcinoma
T
Ryusuke Hiraia, , Yukinobu Sakataa, Akiyuki Tanizawaa, Shinichiro Morib ⁎
a b
Corporate Research and Development Center, Toshiba Corporation, Kanagawa 212 8582, Japan Research Center for Charged Particle Therapy, National Institute of Radiological Sciences, Chiba 263 8555, Japan
ARTICLE INFO
ABSTRACT
Keywords: Markerless tumor tracking Regression model Particle beam therapy Image guidance
Purpose: We have developed a new method to track tumor position using fluoroscopic images, and evaluated it using hepatocellular carcinoma case data. Methods: Our method consists of a training stage and a tracking stage. In the training stage, the model data for the positional relationship between the diaphragm and the tumor are calculated using four-dimensional computed tomography (4DCT) data. The diaphragm is detected along a straight line, which was chosen to avoid 4DCT artifact. In the tracking stage, the tumor position on the fluoroscopic images is calculated by applying the model to the diaphragm. Using data from seven liver cases, we evaluated four metrics: diaphragm edge detection error, modeling error, patient setup error, and tumor tracking error. We measured tumor tracking error for the 15 fluoroscopic sequences from the cases and recorded the computation time. Results: The mean positional error in diaphragm tracking was 0.57 ± 0.62 mm. The mean positional error in tumor tracking in three-dimensional (3D) space was 0.63 ± 0.30 mm by modeling error, and 0.81–2.37 mm with 1–2 mm setup error. The mean positional error in tumor tracking in the fluoroscopy sequences was 1.30 ± 0.54 mm and the mean computation time was 69.0 ± 4.6 ms and 23.2 ± 1.3 ms per frame for the training and tracking stages, respectively. Conclusions: Our markerless tracking method successfully estimated tumor positions. We believe our results will be useful in increasing treatment accuracy for liver cases.
1. Introduction One of the major characteristics of charged particle beam therapy is the improvement it offers in dose conformation. This is facilitated by the Bragg peak allowing for the beam to stop at the distal edge of a tumor, rather than proceeding through the full thickness of the patient, as occurs with photon beams. However, achieving this requires exquisitely accurate treatment planning [1,2]. For beam treatment of thoracic and abdominal tumors, such as liver and lung cancers, monitoring the tumor position during the treatment is very important for achieving high accuracy in dose conformation because the tumors are subject to respiratory motion. In the absence of tumor position monitoring an accurate dose is not guaranteed. Several treatment centers perform liver and lung cancer treatments controlled by respiratory gating, wherein the “beam-on” time is restricted to the time when the tumor is within the planning target volume (PTV). Tumor position during treatment can be captured using both internal and external methods. One example of external tracking is to monitor
⁎
body surface movement using infrared cameras to follow markers placed on the patient’s abdomen [3]. However, this type of system observes respiratory motion at a single point, which may be subject to inaccuracy [4], and accordingly requires that the artificial markers be positioned to ensure a good relationship between the marker and tumor motion [5]. The internal method tracks fiducial markers, which are implanted adjacent to or in the tumor and are visible on X-ray images. When the 3D coordinates of a marker are within the gating window, the treatment beam is turned on. A high accuracy marker tracking technique for liver cancers was proposed by Li et al [6]. Although this method allows the gated irradiation of mobile targets despite irregular respiratory motion, it carries the patient risks of surgery, including pneumothorax and infection [7,8]. Several methods that do not require fiducial markers have been developed and are already being used in clinical applications [9]. Cui et al developed a tumor tracking system as a function of the respiratory cycle using multi-template matching [10]. The tumor position can be estimated by comparing the incoming fluoroscopic images with the
Corresponding author. E-mail address: [email protected] (R. Hirai).
https://doi.org/10.1016/j.ejmp.2020.02.001 Received 17 May 2019; Received in revised form 27 December 2019; Accepted 1 February 2020 1120-1797/ © 2020 Associazione Italiana di Fisica Medica. Published by Elsevier Ltd. All rights reserved.
Physica Medica 70 (2020) 196–205
R. Hirai, et al.
Fig. 1. Tracking calculation procedures. (a) Step 1: Operator-set ROIs include the diaphragm boundary for a single respiratory cycle on DRR images and the training fluoroscopic images at the same position. (b) Step 2: The trajectories of the diaphragm on the respective lines (U_(1…l), V_(1…l)) are obtained along the SI direction ̃ as a function of a time. (c) Step 3: A suitable trajectory (U_l ,V_l )̃ is chosen from the lines. (d) Step 4: A linear regression model of the positional relationship between the detected diaphragm and the tumor on 4DCT designated for making a treatment plan is built. (e) Step 5: The tumor positions on the tracking X-ray sequences are calculated by applying the built model to the diaphragm on each frame.
template fluoroscopic images. However, this may require slightly more time to set the tumor position on the respective template images before treatment. In addition, it is difficult in some cases to manually input liver tumor position on the fluoroscopic images due to low contrast resolution. For the same reason, it is difficult to apply the method of tracking liver tumors with the image pattern as a clue, assuming that the tumor is clearly imaged such as lung cancers [11–13]. To solve these problems, Li et al proposed that the template images be made from planning four-dimensional computed tomography (4DCT) data [7]. However, a drawback of this method is the reduced tracking accuracy when interfractional changes, such as patient positional changes and bowel gas location changes, occur between the planning CT and treatment stages. To reduce the effect of interfractional changes, an indirect method of tracking the tumor using the boundary of the diaphragm as a surrogate has been proposed [14,15]. These authors found a close correlation between diaphragmatic movement and tumor position on fluoroscopic images, and showed that the positional relationship could be formulated by a linear regression model. They also studied the interfractional changes, and reported that their influence can be reduced by applying the linear regression model to each treatment fraction. Nevertheless, this method still required tumor positions on fluoroscopic images to build the model file. Several groups have proposed model-based methods to track the diaphragm [16–18] using extraction of a region with high intensity gradient on the image and then applying the model in which the shape of the diaphragm boundary is represented by a circle or a quadratic function. However, the accuracy of these methods is vulnerable to degradation by four-dimensional computed tomography (4DCT) breathing artifact. Respiratory-induced organ motion compromises CT geometrical accuracy because the CT reconstruction algorithm was developed for static objects. Tumor tracking method using a deep neural network (DNN) has been proposed [19]. This method is applicable to lung and liver tumors. However, it takes much time to reconstruct the DNN's model file, so it is
difficult to make a modified model immediately when tracking accuracy is poor in actual treatment. To solve this problem, we developed a new method to track tumor position using a linear regression model from 4DCT image data. We evaluated the accuracy using liver image data sets. 2. Materials and methods 2.1. Image acquisition A retrospective analysis was performed on a total of 21 fluoroscopic image data sets (=seven liver data sets on each of three different days), which were processed as part of our treatment workflow at National Institute of Radiological Science [9]. These images were acquired during treatment using a dynamic flat panel detector (DFPD) (PaxScan 3030®, Varian Medical Systems, Palo Alto CA, USA). 4DCT images were acquired under free-breathing conditions using a 320-detector row CT (Aquilion One Vision®, Canon Medical Systems, Otawara, Japan). 4DCT images were subdivided into 10 phases (T00: peak inhalation). CT section thickness was 1.0 mm. A board-certified oncologist manually contoured the gross tumor volume (GTV) on the CT data at T50, and these contours at other respiratory phases were then automatically calculated using B-Spline-based deformable image registration [22]. Reference tumor positions in the respective phases were defined as the center of mass of the GTV. Fluoroscopic image data were acquired during treatment delivery using the dynamic flat panel detectors installed on each side of the vertical irradiation port, with the corresponding x-ray tubes installed under the floor of our treatment room. The distance from the tube to the room isocenter and the detector were 170 cm and 240 cm, respectively. Each fluoroscopic image consisted of about 200 frames acquired at 15 fps. Fluoroscopic image matrix size was 768 × 768, with pixel sizes of 0.388 × 0.388 mm.
197
Physica Medica 70 (2020) 196–205
R. Hirai, et al.
2.2. Tracking algorithm
diaphragm using the whole frame.
We introduced five procedures to implement our tracking algorithm, as follows (Fig. 1).
2.2.3. Step 3: Selecting optimum calculation lines A single suitable trajectory was selected from respective SI calculation lines by calculating diaphragmatic positional correlation in between the trajectories of the ROIs and diaphragmatic position as calculated in Step 2 in DRR and the training fluoroscopic images. We can optimize line l using the following equation:
2.2.1. Step 1: Setting an identical region of interest in DRR and fluoroscopic images Digitally-reconstructed radiography (DRR) images were projected using the treatment planning 4DCT data along a virtual X-ray source beam path to the target isocenter. An identical region of interest (ROI) were set on both the 4D-DRR images derived from the 4DCT planning data sets and the training fluoroscopic images, taking care to include the diaphragm boundary on both images for a single respiratory cycle (Fig. 1a).
l = arg max R (Ul, U ) R (Vl , V ) l = 1,
where Ul = {u1, l, , uN , l and Vl = {v1, l, , vM , l are the diaphragmatic trajectories of l in the DRR and fluoroscopic images, respectively. N and M are the number of DRR and training fluoroscopic images. U = {u1, , uN } and V = {v1, , vM }T are the trajectories of the paired ROIs in the DRR and training fluoroscopy calculated in Step 1, respectively. R (x , y ) is the correlation coefficient of the two vectors x and y ,
F
, xF , l ) = arg
x 1, l ,
max
, xF , l Rl
t=1
R (x , y ) =
d (xt + 1, l , xt + 1, l ) t=1
Yi = Ax i , i = 1,
(1)
2
1, l
+ (x t
xt 1)
N
(yi
m y )2 N
(6)
(7)
,N
(8)
where YD = [Y1, , YN is anN × 3 matrix that represents the tumor trajectory, which was set on the planning 4DCT and XD = [1N , Ul1, Ul2] is anN × 3 matrix that represents the trajectory of the diaphragm on each optimum line, where 1N is theN × 1 vector, all of whose elements = 1.
]T
2.2.5. Step 5: Tracking the tumor position on fluoroscopic images To estimate the tumor position on the fluoroscopic images acquired during treatment, we calculate the position by obtaining the tumor position Yt in the 3D position at the time the t-th fluoroscopic image is acquired by the following equation:
(3)
where is a parameter adjusted to be inversely proportional to the frame rate of the image acquisition. x t , l is calculated by
x t,l = xt
m y) N N i=1
A = (XDT XD ) 1XDT YD
x t , l )2 2
mx
m x )(yi )2
where N is the number of 4DCT frames, Yi = [Xi , Yi , Zi is the 3D position of the tumor, and A is the 3 × 3 matrix, which represents a linear regression model. x i = [ui, l1, ui, l2 , 1]T , where ui, l1 and ui, l2 are the positions of the diaphragm on each DRR image. We can obtain the matrix A using the following equation:
(2)
(x t , l
(x i
(x i
]T
where dIt (xl ) dxl is the spatial derivative of image intensity of the tth frame with respect to the SI direction on the line l . That is, in the Eq. (2), the evaluation becomes higher as the intensity gradient is stronger. The second summation in Eq. (1) becomes higher in value as it closely follows the smooth diaphragmatic trajectory. The functiond (x t , l , xt , l ) to calculate the distance of x t , l and x t , l is defined as:
d (x t , l , xt , l ) = exp
N i=1
2.2.4. Step 4: Designing the linear regression model We obtained suitable lines l1 andl2 by Steps 1–3 in the two DRR images projected by the planning 4DCT along the two virtual X-ray source positions. A linear regression model was designed using the detected diaphragm positions along the optimum calculation line and the tumor position on the planning 4DCT.
where l is the index of a line within the ROI, Rl is the set of pixel positions on the line l within the ROI, x t , l represents the diaphragm position on the line l in the tth frame, F is the frame number, and is a parameter that adjusts the weights of the two evaluation functions. The function c (xt , l ) is defined as:
dIt (xl ) c (x t , l ) = dxl
N i=1
where N is the vector length of x and y , and x i and yi are the ith elements of the vectors x and y , and m x and m y is the mean values of x i and yi respectively.
F 1
c (x t , l ) +
}T
}T
2.2.2. Step 2: Detecting diaphragmatic motion along multiple lines The boundary of the diaphragm searches for a high intensity gradient as in the conventional method [14,15], but image noise may cause faulty detection. The accuracy of tracking therefore can be expected to improve the estimation of the trajectory of the boundary of the diaphragm where the high in correlation with the trajectory of the ROI including the diaphragm. The ROI trajectory ( X = {x1, , xF }T ) on both the 4D-DRR and fluoroscopic images were calculated in the superior-inferior (SI) direction using template matching based on the zero-mean normalized cross-correlation (Lewis et al 1995), where x i , i = 1, , F is the ROI position of the ith frame and F and T represent the number of 4D-DRR or fluoroscopic image frames and the permutation, respectively. Then, the diaphragm trajectory ( X = {x1, , xF }T ) on respective lines along the SI direction in the ROI can be obtained by maximizing the following objective function:
G (x1, l ,
(5)
,L
vt , l1 Yt = A vt , l2 1
(4)
(9)
where A is the linear regression model obtained in step 4, and vt , l1, vt , l2 is the position of the diaphragm on the optimum lines in the two fluoroscopic images. vt , l can be computed by:
where x t is the ROI position. That is, in Eq. (3), the evaluation becomes higher as the trajectory of the boundary of the diaphragm is higher in correlation with the trajectory of the ROI. To solve Eq. (1), we exploited dynamic programming [20,21], which allowed us to obtain the best tracking trajectory of the
vt , l = arg max c (vt , l ) + d (vt , l, vt vt , l R l
198
1, l )
(10)
Physica Medica 70 (2020) 196–205
R. Hirai, et al.
Table 1 Tracking accuracy of diaphragm edge position in pixels by our method for seven cases. (mean + standard deviation) [mm]. Patient number
Tracking error
1 2 3 4 5 6 7 Average
0.58 0.65 0.57 0.66 0.46 0.59 0.51 0.57
± ± ± ± ± ± ± ±
Table 3 Tumor tracking accuracy for 21 fluoroscopic sequences. (mean ± standard deviation) [mm].
0.52 0.55 0.90 0.95 0.65 0.38 0.40 0.62
Patient number
1st treatment
2nd treatment
3rd treatment
Total
1 2 3 4 5 6 7 Average
0.62 1.19 1.92 1.32 1.24 1.26 0.47 –
1.17 1.88 1.64 1.13 1.15 1.15 1.03 –
0.68 2.51 1.65 1.69 1.13 0.52 1.88 –
0.82 1.86 1.74 1.38 1.17 0.98 1.13 1.30
± ± ± ± ± ± ±
0.18 0.56 0.87 0.58 0.56 0.53 0.19
± ± ± ± ± ± ±
0.38 1.14 0.67 0.55 0.69 0.56 0.24
± ± ± ± ± ± ±
0.37 0.62 0.74 0.59 0.72 0.25 0.30
± ± ± ± ± ± ± ±
0.31 0.77 0.76 0.57 0.66 0.45 0.24 0.54
the tracking accuracy of the algorithm. To quantify this, we manually input reference diaphragm edge positions and compared them to the calculated diaphragm edge positions on the fluoroscopic images. We measured the distance between the reference and calculated diaphragm edge positions on the DRR image in pixels, which allowed us to calculate the mean ( µe ) and standard deviation ( e ) of the error. To evaluate the model and setup error including the diaphragmatic edge detection error, we added random values generated by normal distribution, N(µe , e ), to the diaphragmatic edge positions, and compared the recalculated tumor positions with each reference position.
where l is the index of a line within the ROI, and Rl is the set of pixel positions on the line l within the ROI, and the first function is the same as Eq. (2) using treatment fluoroscopic images, and the second function is the same as Eq. (3). 2.3. Evaluation According to clinical conditions, tracking error could be caused by any of four factors: diaphragm edge detection error, modeling error, patient setup error, and inter-fractional changes. The first three metrics assess inaccuracies in tracking. To evaluate the impact of inter-fractional changes, we measured the tumor tracking accuracy for 3 treatments. To emphasize real-time processing, we measured the computation time. All evaluations were performed with the following parameters: λ in Eq. (1) was set to 0.8, and σ in Eq. (4) was set to 5 and 2 for the DRR and fluoroscopic images, respectively, because the frame rate in the fluoroscopic image acquisitions was higher than that in the DRR images.
2.3.2. Modeling error Our tracking algorithm made a linear regression model by using the diaphragm and tumor positions on DRR images. The diaphragm positions were set manually for this test. The tumor positions were derived from projection of the tumor centroid, which was contoured in the 4DCT onto the DRR images and were regarded as the reference tumor positions. This model is, however, just approximated input data, and the model might therefore not always export the same tumor positions even though the same diaphragm positions are imported. To quantify this, we compared the tumor positions exported from the model file and the diaphragm position with and without the diaphragmatic edge detection error to the reference tumor positions on DRR images, in terms of their Euclidean distance.
2.3.1. Diaphragmatic edge detection error Since the quality of fluoroscopic images showed high spatial resolution but low contrast resolution compared to DRR images, diaphragmatic edge detection error on the fluoroscopic images could affect
Table 2 Tracking accuracy with the modeling and setup error. (mean ± standard deviation) [mm]. The first and second rows show diaphragm position without and with the tracking error. Pt.
Diaphragm position tracking error
Model error
no. 1 2 3 4 5 6 7 Total
Setup error in LR ± 1 mm
without with without with without with without with without with without with without with without with
0.64 0.67 0.89 0.93 0.81 0.88 0.48 0.67 0.62 0.79 0.55 0.99 0.41 0.44 0.63 0.77
± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±
0.23 0.28 0.45 0.49 0.45 0.47 0.20 0.36 0.24 0.39 0.37 0.47 0.15 0.21 0.30 0.38
1.37 1.42 1.52 1.54 1.44 1.49 1.23 1.29 1.52 1.56 1.40 1.58 1.28 1.23 1.39 1.45
± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±
Setup error in AP ± 2 mm
0.37 0.44 0.53 0.59 0.43 0.52 0.20 0.32 0.46 0.61 0.68 0.80 0.20 0.20 0.41 0.50
2.27 2.31 2.42 2.44 2.23 2.26 2.20 2.24 2.79 2.80 2.39 2.45 2.15 2.11 2.35 2.37
Abbreviations: LAT = lateral; AP = anterior-posterior; SI = superior-inferior.
199
± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±
± 1 mm 0.33 0.37 0.44 0.50 0.33 0.38 0.20 0.32 0.56 0.70 0.74 0.87 0.19 0.21 0.40 0.48
1.35 1.33 1.42 1.46 1.34 1.37 1.19 1.24 1.43 1.46 1.39 1.58 1.08 1.10 1.31 1.36
± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±
Setup error in SI ± 2 mm
0.52 0.58 0.53 0.62 0.43 0.52 0.26 0.34 0.40 0.46 0.61 0.64 0.24 0.26 0.43 0.49
2.27 2.30 2.22 2.24 2.25 2.29 2.15 2.18 2.41 2.42 2.15 2.25 2.06 2.07 2.21 2.25
± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±
± 1 mm 0.68 0.71 0.49 0.56 0.51 0.59 0.26 0.32 0.34 0.40 0.33 0.40 0.22 0.24 0.40 0.46
0.73 0.69 0.98 1.05 0.86 0.93 0.67 0.78 0.80 0.90 0.92 1.16 0.75 0.74 0.81 0.89
± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±
± 2 mm 0.23 0.30 0.58 0.62 0.48 0.51 0.23 0.40 0.36 0.45 0.58 0.61 0.20 0.21 0.38 0.44
0.78 0.74 1.16 1.22 0.97 1.02 1.01 1.05 0.99 1.09 1.30 1.53 1.22 1.24 1.22 1.13
± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±
0.28 0.39 0.62 0.67 0.48 0.53 0.19 0.42 0.39 0.57 0.90 0.90 0.18 0.19 0.18 0.52
Physica Medica 70 (2020) 196–205
R. Hirai, et al.
Fig. 2. Example of the best results for the tracking position along SI direction (Case 1) (left axis) and tracking positional error in Euclidian distance (right axis), (a) 1st treatment, (b) 2nd treatment, (c) 3rd treatment.
2.3.3. Setup error The diaphragmatic position could change due to patient setup error. To simulate the setup error quantitatively, the 4DCT images were shifted ± 1 mm and ± 2 mm along the lateral (LAT), SI, and anteriorposterior (AP) directions. Then the tumor positions were calculated
using the same model as that used to evaluate the modeling error and the diaphragmatic positions that were manually set on the DRR images projected by the shifted 4DCT. The reference tumor positions were acquired from the shifted 4DCT and were compared to those with and without the diaphragmatic edge detection error in terms of their
200
Physica Medica 70 (2020) 196–205
R. Hirai, et al.
Fig. 3. Example of the worst results for the tracking position along SI direction (Case 2) (left axis) and tracking positional error in Euclidian distance (right axis), (a) 1st treatment, (b) 2nd treatment, (c) 3rd treatment.
Euclidean distance.
medical physicist with comparing the image pattern around GTV on 4DDRR. We evaluated the tumor tracking positional error using the fluoroscopic images (200 frames × 3 treatments for respective cases). The model file was calculated using 4DCT images without any setup error. The tumor position was then derived from the model file by importing fluoroscopic images. The tumor tracking error was evaluated
2.3.4. Tumor tracking accuracy Reference tumor positions on each DFPD image were calculated by 2D-2D image registration using the 4D-DRR with GTV contours and DFPD images, and were modified manually by both the oncologist and
201
Physica Medica 70 (2020) 196–205
R. Hirai, et al.
3.3. Setup error For the setup error along the LAT and AP directions, the positional error between the reference and calculated tumor positions without the diaphragmatic edge detection error was 1.39 ± 0.41 mm and 1.45 ± 0.50 mm with ± 1 mm shift and 2.35 ± 0.40 mm and 2.37 ± 0.48 mm with ± 2 mm shift (Table 2). That with the diaphragmatic edge detection error was 1.31 ± 0.43 mm and 1.36 ± 0.49 mm with ± 1 mm shift and 2.21 ± 0.40 mm and 2.25 ± 0.46 mm with ± 2 mm shift (Table 2). For the setup error in the SI direction, the error was smaller without and with the diaphragmatic edge detection error (=0.81 ± 0.38 mm and 0.89 ± 0.44 mm with ± 1 mm shift and 1.22 ± 0.18 and 1.13 ± 0.52 mm with ± 2 mm shift) than those along the AP and LAT directions (Table 2). 3.4. Tumor tracking error The tumor tracking positional error averaged over all cases was 1.30 ± 0.54 mm (Table 3). The positional error averaged over all treatment fractions was 0.82–1.86 mm and that of the standard deviation was 0.31–0.77 mm. Figs. 2 and 3 show the best (Case 1) and worst (Case 2) cases of tumor tracking along the SI direction by our method. According to the best case (Fig. 2), the ranges of tumor movement (max − min reference position of reference SI direction) for three treatments were 9.8 mm, 10.9 mm, and 11.2 mm, for a mean range of 10.7 ± 0.60 mm. The values for the worst case (Fig. 3) were 28.7 mm, 13.7 mm, and 14.1 mm, for a mean of 18.8 ± 7.00 mm. Fig. 4 shows two examples of diaphragm tracking results on the DRR images with the line that our method selected in Step 3. These DRR images showed 4DCT artifact in the gray circle. It can be seen that the location of the selected lines by our algorithm was intended to avoid 4DCT artifact.
Fig. 4. (a) Examples of diaphragmatic detection error on DRR images: (left panel) with a 4DCT artifact at T10 (marked with white arrow), and (right panel) without a 4DCT artifact at T20 (Case 2). (b) Examples of diaphragmatic detection error on DRR images: (left panel) with a 4DCT artifact at T20 (marked as white arrow), and (right panel) without a 4DCT artifact at T30 (patient 5).
using the calculated and the reference tumor positions in Euclidean distance. 2.3.5. Computation time Computation time for the training (steps 1–4) and tracking (step 5) stages was averaged over all fluoroscopic images. Our tracking algorithm was programmed using the C++ program language (Microsoft Visual Studio 2012®, Microsoft, Redmond WA, USA). It works under a Windows 7 environment and is installed on a workstation (Intel Xeon® [email protected] GHz, 32 GB physical memory, Hewlett-Packard, Inc., Palo Alto CA, USA).
3.5. Computation time The mean computation time averaged over all cases was < 114 ms (69.0 ± 4.6 ms) and 25 ms (23.2 ± 1.3) per frame for training (Steps 2–4) and tracking stages (Step 5), respectively. 4. Discussion
3. Results
We evaluated the basic performance of our tumor tracking algorithm using seven liver cases. Overall, the positional error averaged over all data sets was 1.30 ± 0.54 mm.
3.1. Diaphragm edge detection error Table 1 shows the diaphragmatic tracking accuracy by our method for seven cases. The diaphragmatic edge detection on the fluoroscopic images in comparison with the reference and calculated tumor position was 0.57 ± 0.62 mm.
4.1. Diaphragmatic edge detection error In all cases, the diaphragmatic edge detection error was small; therefore the increase in the mean tracking error was small. Thus, our method is robust for any case and will track the diaphragm accurately. Our method is also robust for correction of 4DCT artifacts, since it avoids use of a diaphragm contour susceptible to artifact (Fig. 4) with selecting optimum calculation lines from step 3 in our algorithm. Using a larger ROI reduces the likelihood of artifacts and increases the accuracy of the results calculated by our algorithm, while the computation time would be longer. Adjusting the parameters of 4DCT reconstruction would also be effective in reducing this artifact.
3.2. Modeling error The mean and standard deviation of the positional error between the reference and calculated tumor positions without and with the diaphragmatic edge detection error due to the modeling error for all the cases was 0.63 ± 0.30 mm without diaphragmatic edge detection error and 0.77 ± 0.38 mm with diaphragmatic edge detection error, respectively (Table 2).
202
Physica Medica 70 (2020) 196–205
R. Hirai, et al.
4.2. Modelling error
Therefore it can be considered that the tracking error increased over the treatment because the patient's breathing was becoming more irregular, leading to increased error in the model. We refer to the comparison with the fiducial markers for liver cases in term of the tracking accuracy. Tang et al. (2007) [6] report the maximum marker tracking errors are < 1.0 mm. Our markerless tumor tracking method maximum error was 1.3 mm. We considered this tracking accuracy difference not to be critical. In addition, fiducial markers, are vulnerable to changes in position over the course of therapy (> 2 mm 95% confidence interval [CI] over a few weeks) [23], which can affect dose distribution [24]. Since markerless tumor tracking is not subject to changes in target position, we believe that this advantage overcomes any differences in tracking accuracy.
To reduce tracking error, a more complex model rather than a linear regression model could be used. Cervino et al. (2009) reports results using the linear combination model of eigenvectors obtained by principal component analysis (PCA) of 200 fluoroscopic image sequences as a regression model of the tumor and the diaphragm [11]. In our case, however, it is difficult to use such a model from just ten 4DCT samples. 4.3. Setup error The magnitude of the tumor tracking error with the setup error along the LAT and AP directions was larger than that along the SI direction. This is because the tumor moved in parallel with the diaphragm edge along the SI direction. Thus, the error along the SI direction is caused by model error only, since the model was built neglecting setup errors, while the cause of tumor tracking error along the LAT and AP directions is not only the model error but also the diaphragmatic edge positional variation. For example, when the calculation line was set at the lateral aspect of the hemidiaphragm on the image, the edge position on the line was changed greatly if the diaphragm in the DRR image moved in parallel with a horizontal line; this is equal to the setup error along the LAT or AP directions. It is therefore better that the calculation line is set near the vertex of the diaphragm in the image. Our group previously reported an actual setup accuracy of 1.1 ± 1.2 mm and a tracking error of 1.8 mm [9]. The tracking error estimate for all patients calculated by simulating setup error on the DRRs was 0.89 ± 0.44 mm–2.37 ± 0.48 mm (Table 2). The estimation of the tracking error due to the setup error by the simulation using DRRs is therefore considered to be a reasonable result even when compared with the result in the actual treatment.
4.5. Computation time To reduce the effect of inter-fractional changes, we should use the training (Steps 1–4) at each fraction. If 100 frames are used for training (Step 2–4) stage, which is sufficient to see one respiratory cycle, it will take approximately 6.9 (=69.0 ms × 100 frames) seconds. Adding other processes, such as loading the sequence, would allow our method to be completed within 10 s. It is therefore possible to modify the position of the ROI and correct the model immediately in actual treatment. In this respect, the proposed method is superior to the DNN-based method, which takes hours to modify the model file. We have measured the mean computation time for tracking; it is within 33 ms (30 fps). Thus, our method could be used for real-time markerless tumor tracking in scanned carbon-ion beam treatment using fluoroscopic images for hepatocellular carcinoma. 5. Study limitations
4.4. Tumor tracking error
We had only seven patient data sets to evaluate. Other studies have also reported on approximately same number of patients [10,15].
Cervino et al. (2010) reports a mean positional error of a tumor as 1.3 mm for five cases using the linear regression model on fluoroscopic image sequences [15], while the mean positional error of the tumor was 1.30 mm using our method. It can be said that the accuracy of tracking tumors of both methods is almost equal, although the error was expected to be larger than that by Cervino et al [15] due to the influence of patient setup error. We, therefore, consider the accuracy of tracking the diaphragm by our method to be superior to that of Cervino et al. [15]. Comparing the best and worst cases, we found that the range of tumor movement over the three treatments was significantly different. In the best case, the change in range was small over the three treatments. In other words, the patient's breathing was stable and the model error would be small, therefore we can confirm that Case 1 was the best case. In the worst case, the change of the range over three treatments was large. In terms of the tracking error, the first treatment had the smallest error, and the error became larger after the second treatment.
6. Conclusion We present a markerless tumor position estimation method by tracking the diaphragmatic position and learning the relationship between the diaphragm and tumor positions from treatment planning. The high accuracy and short computation time of the proposed method was shown in experiments involving seven liver cases. The proposed method will be useful for treatment based on the above results. Acknowledgments Mr. Hirai, Dr. Sakata and Mr. Tanizawa are employed by the Toshiba Corporation, Kawasaki, Japan. We thank Libby Cone, MD, MA, from DMC Corp. (www.dmed.co.jp) for editing drafts of this manuscript.
Appendix Fig. 5.
203
Physica Medica 70 (2020) 196–205
R. Hirai, et al.
Fig. 5. Experimental results for (a) Case 1, (b) Case 2, (c) Case 3, (d) Case 4, (e) Case 5, (f) Case 6 and (g) Case 7. The left vertical axis shows the tumor position along the SI direction of the reference as estimated by our method. The right vertical axis shows the positional error of the two in 3D space.
204
Physica Medica 70 (2020) 196–205
R. Hirai, et al.
References [13]
[1] Jakel O, Karger CP, Debus J. The future of heavy ion radiotherapy. Med Phys 2008;35:5653–63. [2] Riboldi M, Orecchia R, Baroni G. Real-time tumour tracking in particle therapy: technological developments and future perspectives. Lancet Oncol 2012;13:e383–91. [3] Kubo HD, Hill BC. Respiration gated radiotherapy treatment: a technical study. Phys Med Biol 1996;41:83–91. [4] Meschini G, Seregni M, Pella A, Ciocca M, Fossati P, Valvo F, et al. Evaluation of residual abdominal tumour motion in carbon ion gated treatments through respiratory motion modelling. Phys Med Biol 2017;34:28–37. [5] Shirato H, Harada T, Harabayashi T, Hida K, Endo H, Kitamura K, et al. Feasibility of insertion/implantation of 2.0-mm-diameter gold internal fiducial markers for precise setup and real-time tumor tracking in radiotherapy. Int J Radiat Oncol Biol Phys 2003;56:240–7. [6] Tang X, Sharp GC, Jiang SB. Fluoroscopic tracking of multiple implanted fiducial markers using multiple object tracking. Phys Med Biol 2007;52:4081–98. [7] Li R, Lewis JH, Cervino LI, Jiang SB. A feasibility study of markerless fluoroscopic gating for lung cancer radiotherapy using 4DCT templates. Phys Med Biol 2009;54:N489–500. [8] Lin T, Li R, Tang X, Dy JG, Jiang SB. Markerless gating for lung cancer radiotherapy based on machine learning techniques. Phys Med. Biol. 2009;54:1555–63. [9] Mori S, Karube M, Shirai T, Tajiri M, Takekoshi T, Miki K, et al. Carbon-ion pencil beam scanning treatment with gated markerless tumor tracking: an analysis of positional accuracy. Int J Radiat Oncol Biol Phys 2016;95:258–66. [10] Cui Y, Dy JG, Sharp GC, Alexander B, Jiang SB. Multiple template-based fluoroscopic tracking of lung tumor mass without implanted fiducial markers. Phys Med Biol 2007;52:6229–42. [11] Zhang X, Homma N, Ichiji K, Abe M, Sugita N, Takai Y, et al. A kernel-based method for markerless tumor tracking in kV fluoroscopic images. Phys Med Biol 2014;59:4897–911. [12] Shieh CC, Caillet V, Dunbar M, Keall PJ, Booth JT, Hardcastle N, et al. A Bayesian
[14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]
205
approach for three-dimensional markerless tumor tracking using kV imaging during lung radiotherapy. Phys Med Biol 2017;62:3065–80. Sakata Y, Hirai R, Taguchi Y, Mori S. Marker-less tumor tracking for lung cancer by tumor image pattern learning. Int J Radiat Oncol 2016;96. E651–E651. Cervino LI, Chao AK, Sandhu A, Jiang SB. The diaphragm as an anatomic surrogate for lung tumor motion. Phys Med Biol 2009;54:3529–41. Cervino LI, Jiang Y, Sandhu A, Jiang SB. Tumor motion prediction with the diaphragm as a surrogate: a feasibility study. Phys Med Biol 2010;55:N221–9. Petkov S, Romero A, Suarez XC, Radeva P, Gatta C. Robust and accurate diaphragm border detection in cardiac X-ray angiographies. Proceedings of the 3rd STACOM, part of MICCAI 2012 pp. 225–234. Petkov S, Carrillo X, Radeva P, Gatta C. Diaphragm border detection in coronary Xray angiographies: new method and applications. Comput Med Imaging Graph 2014;38:296–305. Bögel M, Hofmann HG, Hornegger J, Fahrig R, Britzen S, Maier A. Respiratory motion compensation using diaphragm tracking for CBCT a simulation and a phantom study. Int J Biomed Imaging 2013. Article #6, 10 pages. Hirai R, Sakata Y, Tanizawa Mori S. Real-time tumor tracking using fluoroscopic imaging with deep neural network analysis. Physica Med 2019;59:22–9. Barniv Y. Dynamic Programming solution for detecting dim moving target. IEEE Trans Aerosp Electron Syst 1985:144–56. Dreuw P, Deselaers T, Rybach D. Keysers D and Ney H Tracking using dynamic programming for appearance-based sign language recognition. Proc FGR 2006:293–8. Sharp GC, Kandasamy N, Singh H, Folkert M. GPU-based streaming architectures for fast cone-beam CT image reconstruction and demons deformable registration. Phys Med Biol 2007;52:5771–83. Imura M, Yamazaki K, Shirato H, Onimaru R, Fujino M, Shimizu S, et al. Insertion and fixation of fiducial markers for setup and tracking of lung tumors in radiotherapy. Int J Radiat Oncol Biol Phys 2005;63:1442–7. Newhauser WD, Koch NC, Fontenot JD, Rosenthal SJ, Gombos DS, Fitzek MM, et al. Dosimetric impact of tantalum markers used in the treatment of uveal melanoma with proton beam therapy. Phys Med Biol 2007;52:3979–90.
Contents lists available at ScienceDirect
Physica Medica journal homepage: www.elsevier.com/locate/ejmp
Original paper
Regression model-based real-time markerless tumor tracking with fluoroscopic images for hepatocellular carcinoma
T
Ryusuke Hiraia, , Yukinobu Sakataa, Akiyuki Tanizawaa, Shinichiro Morib ⁎
a b
Corporate Research and Development Center, Toshiba Corporation, Kanagawa 212 8582, Japan Research Center for Charged Particle Therapy, National Institute of Radiological Sciences, Chiba 263 8555, Japan
ARTICLE INFO
ABSTRACT
Keywords: Markerless tumor tracking Regression model Particle beam therapy Image guidance
Purpose: We have developed a new method to track tumor position using fluoroscopic images, and evaluated it using hepatocellular carcinoma case data. Methods: Our method consists of a training stage and a tracking stage. In the training stage, the model data for the positional relationship between the diaphragm and the tumor are calculated using four-dimensional computed tomography (4DCT) data. The diaphragm is detected along a straight line, which was chosen to avoid 4DCT artifact. In the tracking stage, the tumor position on the fluoroscopic images is calculated by applying the model to the diaphragm. Using data from seven liver cases, we evaluated four metrics: diaphragm edge detection error, modeling error, patient setup error, and tumor tracking error. We measured tumor tracking error for the 15 fluoroscopic sequences from the cases and recorded the computation time. Results: The mean positional error in diaphragm tracking was 0.57 ± 0.62 mm. The mean positional error in tumor tracking in three-dimensional (3D) space was 0.63 ± 0.30 mm by modeling error, and 0.81–2.37 mm with 1–2 mm setup error. The mean positional error in tumor tracking in the fluoroscopy sequences was 1.30 ± 0.54 mm and the mean computation time was 69.0 ± 4.6 ms and 23.2 ± 1.3 ms per frame for the training and tracking stages, respectively. Conclusions: Our markerless tracking method successfully estimated tumor positions. We believe our results will be useful in increasing treatment accuracy for liver cases.
1. Introduction One of the major characteristics of charged particle beam therapy is the improvement it offers in dose conformation. This is facilitated by the Bragg peak allowing for the beam to stop at the distal edge of a tumor, rather than proceeding through the full thickness of the patient, as occurs with photon beams. However, achieving this requires exquisitely accurate treatment planning [1,2]. For beam treatment of thoracic and abdominal tumors, such as liver and lung cancers, monitoring the tumor position during the treatment is very important for achieving high accuracy in dose conformation because the tumors are subject to respiratory motion. In the absence of tumor position monitoring an accurate dose is not guaranteed. Several treatment centers perform liver and lung cancer treatments controlled by respiratory gating, wherein the “beam-on” time is restricted to the time when the tumor is within the planning target volume (PTV). Tumor position during treatment can be captured using both internal and external methods. One example of external tracking is to monitor
⁎
body surface movement using infrared cameras to follow markers placed on the patient’s abdomen [3]. However, this type of system observes respiratory motion at a single point, which may be subject to inaccuracy [4], and accordingly requires that the artificial markers be positioned to ensure a good relationship between the marker and tumor motion [5]. The internal method tracks fiducial markers, which are implanted adjacent to or in the tumor and are visible on X-ray images. When the 3D coordinates of a marker are within the gating window, the treatment beam is turned on. A high accuracy marker tracking technique for liver cancers was proposed by Li et al [6]. Although this method allows the gated irradiation of mobile targets despite irregular respiratory motion, it carries the patient risks of surgery, including pneumothorax and infection [7,8]. Several methods that do not require fiducial markers have been developed and are already being used in clinical applications [9]. Cui et al developed a tumor tracking system as a function of the respiratory cycle using multi-template matching [10]. The tumor position can be estimated by comparing the incoming fluoroscopic images with the
Corresponding author. E-mail address: [email protected] (R. Hirai).
https://doi.org/10.1016/j.ejmp.2020.02.001 Received 17 May 2019; Received in revised form 27 December 2019; Accepted 1 February 2020 1120-1797/ © 2020 Associazione Italiana di Fisica Medica. Published by Elsevier Ltd. All rights reserved.
Physica Medica 70 (2020) 196–205
R. Hirai, et al.
Fig. 1. Tracking calculation procedures. (a) Step 1: Operator-set ROIs include the diaphragm boundary for a single respiratory cycle on DRR images and the training fluoroscopic images at the same position. (b) Step 2: The trajectories of the diaphragm on the respective lines (U_(1…l), V_(1…l)) are obtained along the SI direction ̃ as a function of a time. (c) Step 3: A suitable trajectory (U_l ,V_l )̃ is chosen from the lines. (d) Step 4: A linear regression model of the positional relationship between the detected diaphragm and the tumor on 4DCT designated for making a treatment plan is built. (e) Step 5: The tumor positions on the tracking X-ray sequences are calculated by applying the built model to the diaphragm on each frame.
template fluoroscopic images. However, this may require slightly more time to set the tumor position on the respective template images before treatment. In addition, it is difficult in some cases to manually input liver tumor position on the fluoroscopic images due to low contrast resolution. For the same reason, it is difficult to apply the method of tracking liver tumors with the image pattern as a clue, assuming that the tumor is clearly imaged such as lung cancers [11–13]. To solve these problems, Li et al proposed that the template images be made from planning four-dimensional computed tomography (4DCT) data [7]. However, a drawback of this method is the reduced tracking accuracy when interfractional changes, such as patient positional changes and bowel gas location changes, occur between the planning CT and treatment stages. To reduce the effect of interfractional changes, an indirect method of tracking the tumor using the boundary of the diaphragm as a surrogate has been proposed [14,15]. These authors found a close correlation between diaphragmatic movement and tumor position on fluoroscopic images, and showed that the positional relationship could be formulated by a linear regression model. They also studied the interfractional changes, and reported that their influence can be reduced by applying the linear regression model to each treatment fraction. Nevertheless, this method still required tumor positions on fluoroscopic images to build the model file. Several groups have proposed model-based methods to track the diaphragm [16–18] using extraction of a region with high intensity gradient on the image and then applying the model in which the shape of the diaphragm boundary is represented by a circle or a quadratic function. However, the accuracy of these methods is vulnerable to degradation by four-dimensional computed tomography (4DCT) breathing artifact. Respiratory-induced organ motion compromises CT geometrical accuracy because the CT reconstruction algorithm was developed for static objects. Tumor tracking method using a deep neural network (DNN) has been proposed [19]. This method is applicable to lung and liver tumors. However, it takes much time to reconstruct the DNN's model file, so it is
difficult to make a modified model immediately when tracking accuracy is poor in actual treatment. To solve this problem, we developed a new method to track tumor position using a linear regression model from 4DCT image data. We evaluated the accuracy using liver image data sets. 2. Materials and methods 2.1. Image acquisition A retrospective analysis was performed on a total of 21 fluoroscopic image data sets (=seven liver data sets on each of three different days), which were processed as part of our treatment workflow at National Institute of Radiological Science [9]. These images were acquired during treatment using a dynamic flat panel detector (DFPD) (PaxScan 3030®, Varian Medical Systems, Palo Alto CA, USA). 4DCT images were acquired under free-breathing conditions using a 320-detector row CT (Aquilion One Vision®, Canon Medical Systems, Otawara, Japan). 4DCT images were subdivided into 10 phases (T00: peak inhalation). CT section thickness was 1.0 mm. A board-certified oncologist manually contoured the gross tumor volume (GTV) on the CT data at T50, and these contours at other respiratory phases were then automatically calculated using B-Spline-based deformable image registration [22]. Reference tumor positions in the respective phases were defined as the center of mass of the GTV. Fluoroscopic image data were acquired during treatment delivery using the dynamic flat panel detectors installed on each side of the vertical irradiation port, with the corresponding x-ray tubes installed under the floor of our treatment room. The distance from the tube to the room isocenter and the detector were 170 cm and 240 cm, respectively. Each fluoroscopic image consisted of about 200 frames acquired at 15 fps. Fluoroscopic image matrix size was 768 × 768, with pixel sizes of 0.388 × 0.388 mm.
197
Physica Medica 70 (2020) 196–205
R. Hirai, et al.
2.2. Tracking algorithm
diaphragm using the whole frame.
We introduced five procedures to implement our tracking algorithm, as follows (Fig. 1).
2.2.3. Step 3: Selecting optimum calculation lines A single suitable trajectory was selected from respective SI calculation lines by calculating diaphragmatic positional correlation in between the trajectories of the ROIs and diaphragmatic position as calculated in Step 2 in DRR and the training fluoroscopic images. We can optimize line l using the following equation:
2.2.1. Step 1: Setting an identical region of interest in DRR and fluoroscopic images Digitally-reconstructed radiography (DRR) images were projected using the treatment planning 4DCT data along a virtual X-ray source beam path to the target isocenter. An identical region of interest (ROI) were set on both the 4D-DRR images derived from the 4DCT planning data sets and the training fluoroscopic images, taking care to include the diaphragm boundary on both images for a single respiratory cycle (Fig. 1a).
l = arg max R (Ul, U ) R (Vl , V ) l = 1,
where Ul = {u1, l, , uN , l and Vl = {v1, l, , vM , l are the diaphragmatic trajectories of l in the DRR and fluoroscopic images, respectively. N and M are the number of DRR and training fluoroscopic images. U = {u1, , uN } and V = {v1, , vM }T are the trajectories of the paired ROIs in the DRR and training fluoroscopy calculated in Step 1, respectively. R (x , y ) is the correlation coefficient of the two vectors x and y ,
F
, xF , l ) = arg
x 1, l ,
max
, xF , l Rl
t=1
R (x , y ) =
d (xt + 1, l , xt + 1, l ) t=1
Yi = Ax i , i = 1,
(1)
2
1, l
+ (x t
xt 1)
N
(yi
m y )2 N
(6)
(7)
,N
(8)
where YD = [Y1, , YN is anN × 3 matrix that represents the tumor trajectory, which was set on the planning 4DCT and XD = [1N , Ul1, Ul2] is anN × 3 matrix that represents the trajectory of the diaphragm on each optimum line, where 1N is theN × 1 vector, all of whose elements = 1.
]T
2.2.5. Step 5: Tracking the tumor position on fluoroscopic images To estimate the tumor position on the fluoroscopic images acquired during treatment, we calculate the position by obtaining the tumor position Yt in the 3D position at the time the t-th fluoroscopic image is acquired by the following equation:
(3)
where is a parameter adjusted to be inversely proportional to the frame rate of the image acquisition. x t , l is calculated by
x t,l = xt
m y) N N i=1
A = (XDT XD ) 1XDT YD
x t , l )2 2
mx
m x )(yi )2
where N is the number of 4DCT frames, Yi = [Xi , Yi , Zi is the 3D position of the tumor, and A is the 3 × 3 matrix, which represents a linear regression model. x i = [ui, l1, ui, l2 , 1]T , where ui, l1 and ui, l2 are the positions of the diaphragm on each DRR image. We can obtain the matrix A using the following equation:
(2)
(x t , l
(x i
(x i
]T
where dIt (xl ) dxl is the spatial derivative of image intensity of the tth frame with respect to the SI direction on the line l . That is, in the Eq. (2), the evaluation becomes higher as the intensity gradient is stronger. The second summation in Eq. (1) becomes higher in value as it closely follows the smooth diaphragmatic trajectory. The functiond (x t , l , xt , l ) to calculate the distance of x t , l and x t , l is defined as:
d (x t , l , xt , l ) = exp
N i=1
2.2.4. Step 4: Designing the linear regression model We obtained suitable lines l1 andl2 by Steps 1–3 in the two DRR images projected by the planning 4DCT along the two virtual X-ray source positions. A linear regression model was designed using the detected diaphragm positions along the optimum calculation line and the tumor position on the planning 4DCT.
where l is the index of a line within the ROI, Rl is the set of pixel positions on the line l within the ROI, x t , l represents the diaphragm position on the line l in the tth frame, F is the frame number, and is a parameter that adjusts the weights of the two evaluation functions. The function c (xt , l ) is defined as:
dIt (xl ) c (x t , l ) = dxl
N i=1
where N is the vector length of x and y , and x i and yi are the ith elements of the vectors x and y , and m x and m y is the mean values of x i and yi respectively.
F 1
c (x t , l ) +
}T
}T
2.2.2. Step 2: Detecting diaphragmatic motion along multiple lines The boundary of the diaphragm searches for a high intensity gradient as in the conventional method [14,15], but image noise may cause faulty detection. The accuracy of tracking therefore can be expected to improve the estimation of the trajectory of the boundary of the diaphragm where the high in correlation with the trajectory of the ROI including the diaphragm. The ROI trajectory ( X = {x1, , xF }T ) on both the 4D-DRR and fluoroscopic images were calculated in the superior-inferior (SI) direction using template matching based on the zero-mean normalized cross-correlation (Lewis et al 1995), where x i , i = 1, , F is the ROI position of the ith frame and F and T represent the number of 4D-DRR or fluoroscopic image frames and the permutation, respectively. Then, the diaphragm trajectory ( X = {x1, , xF }T ) on respective lines along the SI direction in the ROI can be obtained by maximizing the following objective function:
G (x1, l ,
(5)
,L
vt , l1 Yt = A vt , l2 1
(4)
(9)
where A is the linear regression model obtained in step 4, and vt , l1, vt , l2 is the position of the diaphragm on the optimum lines in the two fluoroscopic images. vt , l can be computed by:
where x t is the ROI position. That is, in Eq. (3), the evaluation becomes higher as the trajectory of the boundary of the diaphragm is higher in correlation with the trajectory of the ROI. To solve Eq. (1), we exploited dynamic programming [20,21], which allowed us to obtain the best tracking trajectory of the
vt , l = arg max c (vt , l ) + d (vt , l, vt vt , l R l
198
1, l )
(10)
Physica Medica 70 (2020) 196–205
R. Hirai, et al.
Table 1 Tracking accuracy of diaphragm edge position in pixels by our method for seven cases. (mean + standard deviation) [mm]. Patient number
Tracking error
1 2 3 4 5 6 7 Average
0.58 0.65 0.57 0.66 0.46 0.59 0.51 0.57
± ± ± ± ± ± ± ±
Table 3 Tumor tracking accuracy for 21 fluoroscopic sequences. (mean ± standard deviation) [mm].
0.52 0.55 0.90 0.95 0.65 0.38 0.40 0.62
Patient number
1st treatment
2nd treatment
3rd treatment
Total
1 2 3 4 5 6 7 Average
0.62 1.19 1.92 1.32 1.24 1.26 0.47 –
1.17 1.88 1.64 1.13 1.15 1.15 1.03 –
0.68 2.51 1.65 1.69 1.13 0.52 1.88 –
0.82 1.86 1.74 1.38 1.17 0.98 1.13 1.30
± ± ± ± ± ± ±
0.18 0.56 0.87 0.58 0.56 0.53 0.19
± ± ± ± ± ± ±
0.38 1.14 0.67 0.55 0.69 0.56 0.24
± ± ± ± ± ± ±
0.37 0.62 0.74 0.59 0.72 0.25 0.30
± ± ± ± ± ± ± ±
0.31 0.77 0.76 0.57 0.66 0.45 0.24 0.54
the tracking accuracy of the algorithm. To quantify this, we manually input reference diaphragm edge positions and compared them to the calculated diaphragm edge positions on the fluoroscopic images. We measured the distance between the reference and calculated diaphragm edge positions on the DRR image in pixels, which allowed us to calculate the mean ( µe ) and standard deviation ( e ) of the error. To evaluate the model and setup error including the diaphragmatic edge detection error, we added random values generated by normal distribution, N(µe , e ), to the diaphragmatic edge positions, and compared the recalculated tumor positions with each reference position.
where l is the index of a line within the ROI, and Rl is the set of pixel positions on the line l within the ROI, and the first function is the same as Eq. (2) using treatment fluoroscopic images, and the second function is the same as Eq. (3). 2.3. Evaluation According to clinical conditions, tracking error could be caused by any of four factors: diaphragm edge detection error, modeling error, patient setup error, and inter-fractional changes. The first three metrics assess inaccuracies in tracking. To evaluate the impact of inter-fractional changes, we measured the tumor tracking accuracy for 3 treatments. To emphasize real-time processing, we measured the computation time. All evaluations were performed with the following parameters: λ in Eq. (1) was set to 0.8, and σ in Eq. (4) was set to 5 and 2 for the DRR and fluoroscopic images, respectively, because the frame rate in the fluoroscopic image acquisitions was higher than that in the DRR images.
2.3.2. Modeling error Our tracking algorithm made a linear regression model by using the diaphragm and tumor positions on DRR images. The diaphragm positions were set manually for this test. The tumor positions were derived from projection of the tumor centroid, which was contoured in the 4DCT onto the DRR images and were regarded as the reference tumor positions. This model is, however, just approximated input data, and the model might therefore not always export the same tumor positions even though the same diaphragm positions are imported. To quantify this, we compared the tumor positions exported from the model file and the diaphragm position with and without the diaphragmatic edge detection error to the reference tumor positions on DRR images, in terms of their Euclidean distance.
2.3.1. Diaphragmatic edge detection error Since the quality of fluoroscopic images showed high spatial resolution but low contrast resolution compared to DRR images, diaphragmatic edge detection error on the fluoroscopic images could affect
Table 2 Tracking accuracy with the modeling and setup error. (mean ± standard deviation) [mm]. The first and second rows show diaphragm position without and with the tracking error. Pt.
Diaphragm position tracking error
Model error
no. 1 2 3 4 5 6 7 Total
Setup error in LR ± 1 mm
without with without with without with without with without with without with without with without with
0.64 0.67 0.89 0.93 0.81 0.88 0.48 0.67 0.62 0.79 0.55 0.99 0.41 0.44 0.63 0.77
± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±
0.23 0.28 0.45 0.49 0.45 0.47 0.20 0.36 0.24 0.39 0.37 0.47 0.15 0.21 0.30 0.38
1.37 1.42 1.52 1.54 1.44 1.49 1.23 1.29 1.52 1.56 1.40 1.58 1.28 1.23 1.39 1.45
± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±
Setup error in AP ± 2 mm
0.37 0.44 0.53 0.59 0.43 0.52 0.20 0.32 0.46 0.61 0.68 0.80 0.20 0.20 0.41 0.50
2.27 2.31 2.42 2.44 2.23 2.26 2.20 2.24 2.79 2.80 2.39 2.45 2.15 2.11 2.35 2.37
Abbreviations: LAT = lateral; AP = anterior-posterior; SI = superior-inferior.
199
± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±
± 1 mm 0.33 0.37 0.44 0.50 0.33 0.38 0.20 0.32 0.56 0.70 0.74 0.87 0.19 0.21 0.40 0.48
1.35 1.33 1.42 1.46 1.34 1.37 1.19 1.24 1.43 1.46 1.39 1.58 1.08 1.10 1.31 1.36
± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±
Setup error in SI ± 2 mm
0.52 0.58 0.53 0.62 0.43 0.52 0.26 0.34 0.40 0.46 0.61 0.64 0.24 0.26 0.43 0.49
2.27 2.30 2.22 2.24 2.25 2.29 2.15 2.18 2.41 2.42 2.15 2.25 2.06 2.07 2.21 2.25
± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±
± 1 mm 0.68 0.71 0.49 0.56 0.51 0.59 0.26 0.32 0.34 0.40 0.33 0.40 0.22 0.24 0.40 0.46
0.73 0.69 0.98 1.05 0.86 0.93 0.67 0.78 0.80 0.90 0.92 1.16 0.75 0.74 0.81 0.89
± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±
± 2 mm 0.23 0.30 0.58 0.62 0.48 0.51 0.23 0.40 0.36 0.45 0.58 0.61 0.20 0.21 0.38 0.44
0.78 0.74 1.16 1.22 0.97 1.02 1.01 1.05 0.99 1.09 1.30 1.53 1.22 1.24 1.22 1.13
± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ±
0.28 0.39 0.62 0.67 0.48 0.53 0.19 0.42 0.39 0.57 0.90 0.90 0.18 0.19 0.18 0.52
Physica Medica 70 (2020) 196–205
R. Hirai, et al.
Fig. 2. Example of the best results for the tracking position along SI direction (Case 1) (left axis) and tracking positional error in Euclidian distance (right axis), (a) 1st treatment, (b) 2nd treatment, (c) 3rd treatment.
2.3.3. Setup error The diaphragmatic position could change due to patient setup error. To simulate the setup error quantitatively, the 4DCT images were shifted ± 1 mm and ± 2 mm along the lateral (LAT), SI, and anteriorposterior (AP) directions. Then the tumor positions were calculated
using the same model as that used to evaluate the modeling error and the diaphragmatic positions that were manually set on the DRR images projected by the shifted 4DCT. The reference tumor positions were acquired from the shifted 4DCT and were compared to those with and without the diaphragmatic edge detection error in terms of their
200
Physica Medica 70 (2020) 196–205
R. Hirai, et al.
Fig. 3. Example of the worst results for the tracking position along SI direction (Case 2) (left axis) and tracking positional error in Euclidian distance (right axis), (a) 1st treatment, (b) 2nd treatment, (c) 3rd treatment.
Euclidean distance.
medical physicist with comparing the image pattern around GTV on 4DDRR. We evaluated the tumor tracking positional error using the fluoroscopic images (200 frames × 3 treatments for respective cases). The model file was calculated using 4DCT images without any setup error. The tumor position was then derived from the model file by importing fluoroscopic images. The tumor tracking error was evaluated
2.3.4. Tumor tracking accuracy Reference tumor positions on each DFPD image were calculated by 2D-2D image registration using the 4D-DRR with GTV contours and DFPD images, and were modified manually by both the oncologist and
201
Physica Medica 70 (2020) 196–205
R. Hirai, et al.
3.3. Setup error For the setup error along the LAT and AP directions, the positional error between the reference and calculated tumor positions without the diaphragmatic edge detection error was 1.39 ± 0.41 mm and 1.45 ± 0.50 mm with ± 1 mm shift and 2.35 ± 0.40 mm and 2.37 ± 0.48 mm with ± 2 mm shift (Table 2). That with the diaphragmatic edge detection error was 1.31 ± 0.43 mm and 1.36 ± 0.49 mm with ± 1 mm shift and 2.21 ± 0.40 mm and 2.25 ± 0.46 mm with ± 2 mm shift (Table 2). For the setup error in the SI direction, the error was smaller without and with the diaphragmatic edge detection error (=0.81 ± 0.38 mm and 0.89 ± 0.44 mm with ± 1 mm shift and 1.22 ± 0.18 and 1.13 ± 0.52 mm with ± 2 mm shift) than those along the AP and LAT directions (Table 2). 3.4. Tumor tracking error The tumor tracking positional error averaged over all cases was 1.30 ± 0.54 mm (Table 3). The positional error averaged over all treatment fractions was 0.82–1.86 mm and that of the standard deviation was 0.31–0.77 mm. Figs. 2 and 3 show the best (Case 1) and worst (Case 2) cases of tumor tracking along the SI direction by our method. According to the best case (Fig. 2), the ranges of tumor movement (max − min reference position of reference SI direction) for three treatments were 9.8 mm, 10.9 mm, and 11.2 mm, for a mean range of 10.7 ± 0.60 mm. The values for the worst case (Fig. 3) were 28.7 mm, 13.7 mm, and 14.1 mm, for a mean of 18.8 ± 7.00 mm. Fig. 4 shows two examples of diaphragm tracking results on the DRR images with the line that our method selected in Step 3. These DRR images showed 4DCT artifact in the gray circle. It can be seen that the location of the selected lines by our algorithm was intended to avoid 4DCT artifact.
Fig. 4. (a) Examples of diaphragmatic detection error on DRR images: (left panel) with a 4DCT artifact at T10 (marked with white arrow), and (right panel) without a 4DCT artifact at T20 (Case 2). (b) Examples of diaphragmatic detection error on DRR images: (left panel) with a 4DCT artifact at T20 (marked as white arrow), and (right panel) without a 4DCT artifact at T30 (patient 5).
using the calculated and the reference tumor positions in Euclidean distance. 2.3.5. Computation time Computation time for the training (steps 1–4) and tracking (step 5) stages was averaged over all fluoroscopic images. Our tracking algorithm was programmed using the C++ program language (Microsoft Visual Studio 2012®, Microsoft, Redmond WA, USA). It works under a Windows 7 environment and is installed on a workstation (Intel Xeon® [email protected] GHz, 32 GB physical memory, Hewlett-Packard, Inc., Palo Alto CA, USA).
3.5. Computation time The mean computation time averaged over all cases was < 114 ms (69.0 ± 4.6 ms) and 25 ms (23.2 ± 1.3) per frame for training (Steps 2–4) and tracking stages (Step 5), respectively. 4. Discussion
3. Results
We evaluated the basic performance of our tumor tracking algorithm using seven liver cases. Overall, the positional error averaged over all data sets was 1.30 ± 0.54 mm.
3.1. Diaphragm edge detection error Table 1 shows the diaphragmatic tracking accuracy by our method for seven cases. The diaphragmatic edge detection on the fluoroscopic images in comparison with the reference and calculated tumor position was 0.57 ± 0.62 mm.
4.1. Diaphragmatic edge detection error In all cases, the diaphragmatic edge detection error was small; therefore the increase in the mean tracking error was small. Thus, our method is robust for any case and will track the diaphragm accurately. Our method is also robust for correction of 4DCT artifacts, since it avoids use of a diaphragm contour susceptible to artifact (Fig. 4) with selecting optimum calculation lines from step 3 in our algorithm. Using a larger ROI reduces the likelihood of artifacts and increases the accuracy of the results calculated by our algorithm, while the computation time would be longer. Adjusting the parameters of 4DCT reconstruction would also be effective in reducing this artifact.
3.2. Modeling error The mean and standard deviation of the positional error between the reference and calculated tumor positions without and with the diaphragmatic edge detection error due to the modeling error for all the cases was 0.63 ± 0.30 mm without diaphragmatic edge detection error and 0.77 ± 0.38 mm with diaphragmatic edge detection error, respectively (Table 2).
202
Physica Medica 70 (2020) 196–205
R. Hirai, et al.
4.2. Modelling error
Therefore it can be considered that the tracking error increased over the treatment because the patient's breathing was becoming more irregular, leading to increased error in the model. We refer to the comparison with the fiducial markers for liver cases in term of the tracking accuracy. Tang et al. (2007) [6] report the maximum marker tracking errors are < 1.0 mm. Our markerless tumor tracking method maximum error was 1.3 mm. We considered this tracking accuracy difference not to be critical. In addition, fiducial markers, are vulnerable to changes in position over the course of therapy (> 2 mm 95% confidence interval [CI] over a few weeks) [23], which can affect dose distribution [24]. Since markerless tumor tracking is not subject to changes in target position, we believe that this advantage overcomes any differences in tracking accuracy.
To reduce tracking error, a more complex model rather than a linear regression model could be used. Cervino et al. (2009) reports results using the linear combination model of eigenvectors obtained by principal component analysis (PCA) of 200 fluoroscopic image sequences as a regression model of the tumor and the diaphragm [11]. In our case, however, it is difficult to use such a model from just ten 4DCT samples. 4.3. Setup error The magnitude of the tumor tracking error with the setup error along the LAT and AP directions was larger than that along the SI direction. This is because the tumor moved in parallel with the diaphragm edge along the SI direction. Thus, the error along the SI direction is caused by model error only, since the model was built neglecting setup errors, while the cause of tumor tracking error along the LAT and AP directions is not only the model error but also the diaphragmatic edge positional variation. For example, when the calculation line was set at the lateral aspect of the hemidiaphragm on the image, the edge position on the line was changed greatly if the diaphragm in the DRR image moved in parallel with a horizontal line; this is equal to the setup error along the LAT or AP directions. It is therefore better that the calculation line is set near the vertex of the diaphragm in the image. Our group previously reported an actual setup accuracy of 1.1 ± 1.2 mm and a tracking error of 1.8 mm [9]. The tracking error estimate for all patients calculated by simulating setup error on the DRRs was 0.89 ± 0.44 mm–2.37 ± 0.48 mm (Table 2). The estimation of the tracking error due to the setup error by the simulation using DRRs is therefore considered to be a reasonable result even when compared with the result in the actual treatment.
4.5. Computation time To reduce the effect of inter-fractional changes, we should use the training (Steps 1–4) at each fraction. If 100 frames are used for training (Step 2–4) stage, which is sufficient to see one respiratory cycle, it will take approximately 6.9 (=69.0 ms × 100 frames) seconds. Adding other processes, such as loading the sequence, would allow our method to be completed within 10 s. It is therefore possible to modify the position of the ROI and correct the model immediately in actual treatment. In this respect, the proposed method is superior to the DNN-based method, which takes hours to modify the model file. We have measured the mean computation time for tracking; it is within 33 ms (30 fps). Thus, our method could be used for real-time markerless tumor tracking in scanned carbon-ion beam treatment using fluoroscopic images for hepatocellular carcinoma. 5. Study limitations
4.4. Tumor tracking error
We had only seven patient data sets to evaluate. Other studies have also reported on approximately same number of patients [10,15].
Cervino et al. (2010) reports a mean positional error of a tumor as 1.3 mm for five cases using the linear regression model on fluoroscopic image sequences [15], while the mean positional error of the tumor was 1.30 mm using our method. It can be said that the accuracy of tracking tumors of both methods is almost equal, although the error was expected to be larger than that by Cervino et al [15] due to the influence of patient setup error. We, therefore, consider the accuracy of tracking the diaphragm by our method to be superior to that of Cervino et al. [15]. Comparing the best and worst cases, we found that the range of tumor movement over the three treatments was significantly different. In the best case, the change in range was small over the three treatments. In other words, the patient's breathing was stable and the model error would be small, therefore we can confirm that Case 1 was the best case. In the worst case, the change of the range over three treatments was large. In terms of the tracking error, the first treatment had the smallest error, and the error became larger after the second treatment.
6. Conclusion We present a markerless tumor position estimation method by tracking the diaphragmatic position and learning the relationship between the diaphragm and tumor positions from treatment planning. The high accuracy and short computation time of the proposed method was shown in experiments involving seven liver cases. The proposed method will be useful for treatment based on the above results. Acknowledgments Mr. Hirai, Dr. Sakata and Mr. Tanizawa are employed by the Toshiba Corporation, Kawasaki, Japan. We thank Libby Cone, MD, MA, from DMC Corp. (www.dmed.co.jp) for editing drafts of this manuscript.
Appendix Fig. 5.
203
Physica Medica 70 (2020) 196–205
R. Hirai, et al.
Fig. 5. Experimental results for (a) Case 1, (b) Case 2, (c) Case 3, (d) Case 4, (e) Case 5, (f) Case 6 and (g) Case 7. The left vertical axis shows the tumor position along the SI direction of the reference as estimated by our method. The right vertical axis shows the positional error of the two in 3D space.
204
Physica Medica 70 (2020) 196–205
R. Hirai, et al.
References [13]
[1] Jakel O, Karger CP, Debus J. The future of heavy ion radiotherapy. Med Phys 2008;35:5653–63. [2] Riboldi M, Orecchia R, Baroni G. Real-time tumour tracking in particle therapy: technological developments and future perspectives. Lancet Oncol 2012;13:e383–91. [3] Kubo HD, Hill BC. Respiration gated radiotherapy treatment: a technical study. Phys Med Biol 1996;41:83–91. [4] Meschini G, Seregni M, Pella A, Ciocca M, Fossati P, Valvo F, et al. Evaluation of residual abdominal tumour motion in carbon ion gated treatments through respiratory motion modelling. Phys Med Biol 2017;34:28–37. [5] Shirato H, Harada T, Harabayashi T, Hida K, Endo H, Kitamura K, et al. Feasibility of insertion/implantation of 2.0-mm-diameter gold internal fiducial markers for precise setup and real-time tumor tracking in radiotherapy. Int J Radiat Oncol Biol Phys 2003;56:240–7. [6] Tang X, Sharp GC, Jiang SB. Fluoroscopic tracking of multiple implanted fiducial markers using multiple object tracking. Phys Med Biol 2007;52:4081–98. [7] Li R, Lewis JH, Cervino LI, Jiang SB. A feasibility study of markerless fluoroscopic gating for lung cancer radiotherapy using 4DCT templates. Phys Med Biol 2009;54:N489–500. [8] Lin T, Li R, Tang X, Dy JG, Jiang SB. Markerless gating for lung cancer radiotherapy based on machine learning techniques. Phys Med. Biol. 2009;54:1555–63. [9] Mori S, Karube M, Shirai T, Tajiri M, Takekoshi T, Miki K, et al. Carbon-ion pencil beam scanning treatment with gated markerless tumor tracking: an analysis of positional accuracy. Int J Radiat Oncol Biol Phys 2016;95:258–66. [10] Cui Y, Dy JG, Sharp GC, Alexander B, Jiang SB. Multiple template-based fluoroscopic tracking of lung tumor mass without implanted fiducial markers. Phys Med Biol 2007;52:6229–42. [11] Zhang X, Homma N, Ichiji K, Abe M, Sugita N, Takai Y, et al. A kernel-based method for markerless tumor tracking in kV fluoroscopic images. Phys Med Biol 2014;59:4897–911. [12] Shieh CC, Caillet V, Dunbar M, Keall PJ, Booth JT, Hardcastle N, et al. A Bayesian
[14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]
205
approach for three-dimensional markerless tumor tracking using kV imaging during lung radiotherapy. Phys Med Biol 2017;62:3065–80. Sakata Y, Hirai R, Taguchi Y, Mori S. Marker-less tumor tracking for lung cancer by tumor image pattern learning. Int J Radiat Oncol 2016;96. E651–E651. Cervino LI, Chao AK, Sandhu A, Jiang SB. The diaphragm as an anatomic surrogate for lung tumor motion. Phys Med Biol 2009;54:3529–41. Cervino LI, Jiang Y, Sandhu A, Jiang SB. Tumor motion prediction with the diaphragm as a surrogate: a feasibility study. Phys Med Biol 2010;55:N221–9. Petkov S, Romero A, Suarez XC, Radeva P, Gatta C. Robust and accurate diaphragm border detection in cardiac X-ray angiographies. Proceedings of the 3rd STACOM, part of MICCAI 2012 pp. 225–234. Petkov S, Carrillo X, Radeva P, Gatta C. Diaphragm border detection in coronary Xray angiographies: new method and applications. Comput Med Imaging Graph 2014;38:296–305. Bögel M, Hofmann HG, Hornegger J, Fahrig R, Britzen S, Maier A. Respiratory motion compensation using diaphragm tracking for CBCT a simulation and a phantom study. Int J Biomed Imaging 2013. Article #6, 10 pages. Hirai R, Sakata Y, Tanizawa Mori S. Real-time tumor tracking using fluoroscopic imaging with deep neural network analysis. Physica Med 2019;59:22–9. Barniv Y. Dynamic Programming solution for detecting dim moving target. IEEE Trans Aerosp Electron Syst 1985:144–56. Dreuw P, Deselaers T, Rybach D. Keysers D and Ney H Tracking using dynamic programming for appearance-based sign language recognition. Proc FGR 2006:293–8. Sharp GC, Kandasamy N, Singh H, Folkert M. GPU-based streaming architectures for fast cone-beam CT image reconstruction and demons deformable registration. Phys Med Biol 2007;52:5771–83. Imura M, Yamazaki K, Shirato H, Onimaru R, Fujino M, Shimizu S, et al. Insertion and fixation of fiducial markers for setup and tracking of lung tumors in radiotherapy. Int J Radiat Oncol Biol Phys 2005;63:1442–7. Newhauser WD, Koch NC, Fontenot JD, Rosenthal SJ, Gombos DS, Fitzek MM, et al. Dosimetric impact of tantalum markers used in the treatment of uveal melanoma with proton beam therapy. Phys Med Biol 2007;52:3979–90.