Evaluation of deformable image registration (DIR) methods for dose accumulation in nasopharyngeal cancer patients during radiotherapy

The study population consisted of five NPC patients treated using a helical tomotherapy unit (Tomo-Therapy, Inc., Madison, Wisconsin). All the patients underwent intensity modulated radiotherapy (IMRT) with a planned dose of 70 Gy delivered to the gross disease at 2.12 Gy/fraction for a total of 33 fractions with a simultaneous integrated boost technique (SIB) according to the RTOG 022512, keeping the mean parotid dose as low as could be possibly achieved and respecting the tissue tolerance of other normal structures. Patient positioning was ensured by appropriate headrest and a personalized HN and shoulder mask. This study received ethics approval, granted by the institutional research committee.

The planned kVCT images were acquired on a computerized tomography unit (Somatom, SIEMENS, Germany) by using a matrix of 512 × 512 with voxel dimension of 0.976 × 0.976 × 3mm³ for the treatment planning process. The 1^st day MVCT images were also acquired on the helical tomotherapy unit as source images for deformable investigation on the same day of acquisition of the planned kVCT images.

When the radiotherapy treatment was started, the daily MVCT images were acquired on the helical tomotherapy unit prior to each treatment fraction used for patient alignment by using a matrix of 512 x512 with voxel dimension of 0.763 × 0.763 × 4 mm³. Typically, the MVCT scan range covers the entirety of the gross tumour volume (GTV), the clinical target volume (CTV), and the parotid glands bilaterally. The weekly MVCT images on the 1^st, 6^th, 11^th, 16^th, 21^st, 26^th, and 31^st fractions were used as the target images to assess the dose accumulation in this study.

The regions of interest (ROIs) including the target and the organ at risks (OARs) were defined by the radiation oncologist on planned kVCT images for the treatment planning processes. The ROIs on the planned kVCT images were transferred to the first-day MVCT images as the source images for each image set. The same oncologist who localized the target and the OAR for the HT treatment planning process also contoured the GTV, CTV, the bilateral parotid glands, and the spinal cord on the weekly MVCT images as the reference images. These contours were compared to the automatic deformed structure generated by the deformable image registration software.

The popular software suite for DIR and ART, DIRART version 1a developed by Yang (2009) 11, was used to create automatic deformed contours and dose accumulation from the MVCT images of NPC patients. DIRART works complementarily with computational environment for radiotherapy research (CERR) to offer more functions9 and provide the capability of selecting various deformation algorithms, transformation frameworks, and mapping directions for providing the deformation vector field (DVF) in deformable registration procedures.

In this study, the algorithms used were Horn and Schunck (HS) optical flow and demons (DM) combined with the asymmetric (Asy) and symmetric (Sym) transformation framework. Therefore, the study was carried out with four DIR methods, including the asymmetric transformation with the Horn and Schunck optical flow (AsyHS), the asymmetric transformation with the demon algorithm (AsyDM), the symmetric transformation with the Horn and Schunck optical flow (SymHS), and the symmetric transformation with the demon algorithm (SymDM).

A multi-resolution technique was used in this study. For the optimum DIR performance for each algorithm, various parameters were systematically adjusted: four multigrids were used (n = 1, 2, 3, and 4) with 10n to 40n iterations per pass.9 The number of passes for the optical flow algorithm was 6 and the number of passes for the demon algorithm was between 2 and 6. Coarser stages were typically run with a greater number of passes to improve the agreement with the target image prior to resampling at finer resolutions.9

The objective of the validation technique was to evaluate the accuracy of the automatic deformed contour by four DIR methods on weekly MVCT images, including the terms of the volume-based criterion and the deformation field analysis.

Regarding the volume-based criterion, the most common overlap metric is the Dice similarity coefficient (DSC).13 DSC is the metric that computes the number of pixels that overlaps between two volumes. If the images have no overlap, then the DSC is 0, and as the contours become identical, the DSC approaches the value of 1.13 Zimring et al.14 suggested that the satisfactory volume matching should be 70% (DSC of 0.7) or more for adaptive radiotherapy application.

Inverse consistency error (ICE) was used to ensure that the transformations were physically invertible for the deformation field analysis. Optimal transformation is found when the ICE minimizes the distance error.15

The dose accumulation process relied on six steps, as illustrated in Figure 1. Firstly, the ROIs from the planned kVCT images were transferred to the 1^st day MVCT as source images for registration. The four DIR methods were performed between the 1^st day MVCT and the weekly MVCTs (step 2). The deformation vector field was applied for creating the automatic deformed structure of the ROIs and for propagation to the weekly MVCT images (step 3), and the weekly dose distribution was deformed (step 4). The weekly dose deformation values were summed to the accumulated dose (step 5) and compared to the initial planned dose distribution (step 6).

Figure 1

To evaluate the effect of dose accumulation from DIR errors, the reference accumulated dose on ROIs was computed by summing the weekly doses corresponding to the weekly MVCTs defined by the radiation oncologist. Moreover, to ensure that the weekly dose summation from the DIRART software was accurate, dose accumulation from independent software, Planned Adaptive software (TomoTherapy Inc., Madison, WI) was used to compare.

Regarding the evaluation of accumulated dose, for the target volume, the median absorbed dose (D_50%), the near-minimum (D_98%) absorbed dose, and the near-maximum (D_2%) absorbed dose values from each DIR method were assessed, the mean absorbed dose (D_mean) of the bilateral parotid glands and D_2% of the spinal cord were compared to the original planned dose for the OARs investigation. The one-way analysis of variance (ANOVA) test and the paired sample t-test were carried out on each set of comparison metrics to determine the statistical significance, with a threshold of p < 0.05; SPSS statistical software version 17 was used to compare and assess the impact of each of the DIR methods.

Regarding volume variations during the radiotherapy, the percent ratio to the volume at the initial treatment planning of five NPC patients is illustrated in Figure 2. The averages of the NPC patients for the volume variation in the initial plan were significantly different from the averages after the treatment in 3 weeks for GTV, with p-value = 0.025, as demonstrated in Figure 2A, and for CTV, with p-value = 0.020, as demonstrated in Figure 2B. The volume was observed to have decreased by an average of 29.8% (GTV) and 21.0% (CTV) at the end of the treatment course.

Figure 2

As regards the OARs, the right and the left parotid volume variations were significantly different from those of the initial plan after 5 weeks and 4 weeks of treatment, with p-values of 0.017 and 0.026, respectively. The average volume decreased by 40.3% (right) and 43.6% (left) at the end of the treatment.

The results of DIR accuracy were consistent between the volume-based criterion, DSC, and the deformation field analysis, ICE. Figure 3A demonstrates the histogram of the DSC values for all of the ROIs by four DIR methods. The SymDM methods showed significant difference from other methods by the one-way ANOVA analysis, with p-value = 0.00, with the worst performance in terms of volume-based criterion by mean values of DSC = 0.50 ± 0.30, 0.56 ± 0.34, 0.67 ± 0.19, 0.65 ± 0.28, and 0.69 ± 0.19 for GTV, CTV, right parotid, left parotid, and spinal cord, respectively. The average of DSC value was less than 0.7 for all the ROIs, that represented unsatisfactory volume matching for adaptive radiotherapy application.14 Regarding the ICE analysis, the results were consistent with DSC value. Figure 3B illustrates the histogram of the ICE values for all of the ROIs by the four DIR methods. The SymDM method also showed the maximum error in the deformation field analysis in terms of ICE for all of the ROIs throughout the treatment. The one-way ANOVA analysis showed that the SymDM method was significantly different from other methods, with p-value = 0.00. Therefore, the SymDM failed to adequately register on weekly MVCT for the dose accumulation application in this study. This paper focuses on the three highest performing algorithms in the AsyHS, AsyDM, and SymHS methods for dose accumulation.

Figure 3

As regards target dose variation, the median GTV and CTV doses received at the end of treatment were slightly different from those in the initial plan. They were 0.11% (range: 0–0.29%) lower than the initial planned dose. The median dose variations of the GTV and CTV were significantly different from the initial planned dose after 6 weeks of treatment, with p-value = 0.016. Regarding the minimum and the maximum doses, they are represented by near-minimum dose (D_98%) and near-maximum dose (D_2%), respectively. As for the D_98%, they received slightly higher doses than the initial plan deals, with an average variation less than 0.5% (range: 0.29–1.60%). However, the dose at the end of treatment received slightly decreased doses of 0.45% (GTV) and 0.28% (CTV) from the initial doses planned for the D₂.

Regarding organ dose variation, the dose differences tended to increase as the treatment progressed. For the bilateral parotid gland, the discrepancy between the delivered and the planned mean doses was found to have increased by 6.8% (range: 2.2 to 10.9%) for the right parotid and by 15.2% (range: -1.7 to 36.3%) for the left parotid. The average mean parotid dose increased in the ranges of 2.24 ± 0.97 Gy (right) and 5.70 ± 4.12 Gy (left) at the end of the treatment. The mean parotid dose variations were significantly different from the initial plan after 6 weeks (right) and 5 weeks (left) of the treatment, with p-value = 0.049 (right) and p-value = 0.010 (left). The spinal cord dose received increased by 6.4% (range: -1.6 to 13.2%) from the initial plan, with the average near-maximum dose increasing in the range of 1.83 ± 1.5 Gy at the end of the treatment.

For each patient, the running cumulative doses were calculated using the CERR software through the three deformable image registration methods carried out by the DIRART software. Figure 4 demonstrates the 1^st day MVCT image with original bilateral parotid gland (A) and the MVCT image at the 31^st fraction with the automatic deformed contour obtained using the AsyDM method (B). The initial planned dose distribution on the 1^st day MVCT image, as illustrated in Figure 4C, was used for comparison with the accumulated dose distribution at the end of the treatment, as demonstrated in Figure 4D. The variations in the cumulative doses between the delivered dose and the initial planned dose are illustrated in Figure 5 and Figure 6.

Figure 4

Figure 5

Figure 6

Figure 5A illustrates the weekly GTV dose variation from the initial plan with three DIR methods. The average of the median dose difference for all methods at the end of the treatment was lower than that in the initial plan, with 0.34 Gy (0.5%), 0.04 Gy (0.1%), and 0.30 Gy (0.4%) for the AsyHS, AsyDm, and SymHS DIR methods, respectively. However, the reference dose of GTV was found to have decreased by 0.11%, with the accumulated GTV dose at 70.12 Gy (range: 69.9–70.4 Gy), at the end of treatment. The median dose variations of the GTV was significantly different from the initial planned dose after 6 weeks of treatment, with p-value = 0.016.

Regarding the near-minimum dose and the near-maximum dose, the D_98% of GTV in three DIR methods were found to be lower than that in the initial plan, as illustrated in Figure 5B; the average discrepancy of the three DIR methods between the planned dose and the delivered dose was 0.33 Gy (0.5%) at the end of the treatment with the high differences of 1.0% after 5 week for AsyHS method. However, the reference near-minimum dose was found to have increased by 0.3% of the initial planned dose, with 69.2 Gy (range: 68.7–70 Gy). As for the maximum GTV dose consideration, the three DIR methods of D_2% are presented in Figure 5C with gradually decrease in dose difference from initial plan for all methods when the time increase. The average D_2% from the three methods was lower than the initial D_2%, with 0.76 Gy (1.1%), at the end of the treatment. The reference D_2% was found to have decreased by 0.45% of the initial planned dose, with 71.5 Gy (range: 70.9–72.5 Gy).

As regards the median CTV dose, the dose variations tended to be similar to the dose variations of GTV. Figure 5D illustrates the very small median CTV differences of various deformable registration methods from the initial planned dose. The discrepancy at the end of the treatment was lower than that in the initial planned dose, by 0.34 Gy (0.4%), 0.02 Gy (0%), and 0.26 Gy (0.4%) for AsyHS, AsyDm, and SymHS DIR methods, respectively. However, the reference median dose was found to have decreased by 0.11%, with 70.12 Gy (range: 69.9–70.4 Gy) at the end of the treatment. The median dose variations of the CTV was significantly different from the initial planned dose after 6 weeks of treatment, with p-value = 0.016.

Regarding the D₉₈ dose analysis of CTV, the discrepancy between the initial and the delivered dose of the three DIR methods are shown in Figure 5E. The average of the D_98% variations from the initial D_98% in the three DIR methods was 0.89 Gy (1.2%). Figure 5F shows the D_2% of CTV; the average variation from the initial D_2% of the three DIR methods was 0.76 Gy (1.1%).

As regards organ dose accumulation in the three DIR methods, Figure 6 illustrates the weekly dose difference from the initial plan in the three DIR methods for the bilateral parotid gland and the spinal cord. Overall, the dose differences tended to increase as the treatment progressed. Figure 6A shows the mean right parotid dose (D_mean) to be higher than the initial planned dose, by 5.38 Gy (16.0%), 3.38 Gy (10.1%), and 4.84 Gy (14.4%) for the AsyHS, AsyDm, and SymHS methods, respectively. However, the reference mean dose was found to have increased by 2.24 Gy (range: 0.8–3.7 Gy), at 6.82%, at the end of the treatment.

For the left parotid mean dose, as illustrated in Figure 6B, these variations were higher after treatment than those for the initial planned dose, and the discrepancy was by 6.88 Gy (18.3%), 4.12 Gy (11.0%), and 6.82 Gy (18.1%) for the AsyHS, AsyDm, and SymHS DIR methods, respectively. However, the reference mean dose was found to have increased by 5.7 Gy (range: -0.6 to12.4 Gy), at 15.2% at the end of the treatment. The mean parotid dose variations were significantly different from the initial plan after 6 weeks (right parotid) and 5 weeks (left parotid) of the treatment, with p-value of 0.049 and 0.010, respectively.

As regards spinal cord weekly dose accumulation, the variations tended to increase in all the three DIR methods, by 2.33 Gy (7.9%), 1.46 Gy (4.9%) and 0.60 Gy (5.4%) for AsyHS, AsyDm, and SymHS, respectively. However, the reference cord dose variation was found to have increased by 1.83 Gy (range: -0.5 to 4.0 Gy), at 6.37% at the end of the treatment.

To ensure dose summation on the DIRART software, the independent planned adaptive software on helical tomotherapy treatment planning was used to compare the dose accumulation values. The same data set of all the reference ROIs on weekly MVCTs defined by the radiation oncologist was transferred to the planned adaptive software. The comparison of the weekly accumulated doses between DIRART and planned adaptive is illustrated in Figure 7. The variations in the accumulated median parotid doses of DIRART were not significantly different according to the planned adaptive software, with p-value = 0.972 for the right parotid gland, as shown in Figure 7A, and p-value = 0.958 for the left parotid gland, as shown in Figure 7B. The consistency in the dose variations between the two independent types of software demonstrates that the dose accumulation of the DIRART software can be applied for use in dose accumulation studies.

Figure 7