A Comparative Evaluation of 3 Different Free-Form Deformable Image Registration and Contour Propagation Methods for Head and Neck MRI: The Case of Parotid Changes During Radiotherapy.

PURPOSE: To validate and compare the deformable image registration and parotid contour propagation process for head and neck magnetic resonance imaging in patients treated with radiotherapy using 3 different approaches-the commercial MIM, the open-source Elastix software, and an optimized version of it.
MATERIALS AND METHODS: Twelve patients with head and neck cancer previously treated with radiotherapy were considered. Deformable image registration and parotid contour propagation were evaluated by considering the magnetic resonance images acquired before and after the end of the treatment. Deformable image registration, based on free-form deformation method, and contour propagation available on MIM were compared to Elastix. Two different contour propagation approaches were implemented for Elastix software, a conventional one (DIR_Trx) and an optimized homemade version, based on mesh deformation (DIR_Mesh). The accuracy of these 3 approaches was estimated by comparing propagated to manual contours in terms of average symmetric distance, maximum symmetric distance, Dice similarity coefficient, sensitivity, and inclusiveness.
RESULTS: A good agreement was generally found between the manual contours and the propagated ones, without differences among the 3 methods; in few critical cases with complex deformations, DIR_Mesh proved to be more accurate, having the lowest values of average symmetric distance and maximum symmetric distance and the highest value of Dice similarity coefficient, although nonsignificant. The average propagation errors with respect to the reference contours are lower than the voxel diagonal (2 mm), and Dice similarity coefficient is around 0.8 for all 3 methods.
CONCLUSION: The 3 free-form deformation approaches were not significantly different in terms of deformable image registration accuracy and can be safely adopted for the registration and parotid contour propagation during radiotherapy on magnetic resonance imaging. More optimized approaches (as DIR_Mesh) could be preferable for critical deformations.

Introduction

Large anatomical variations can be observed during head and neck cancer intensity-modulated radiation therapy, such as weight loss, primary tumor shrinkage, parotid gland (PG) volume reduction, and decrease of the neck diameter.[1-3] Consequently, the actual (accumulated) delivered dose does not correspond to the planned one, with a risk of overdosing organs at risk, in particular PGs.[4,5]

The application of deformable image registration (DIR) has become increasingly present in modern image-guided radiation therapy (IGRT[6-8])—for routine treatment planning, it is very common to register magnetic resonance imaging (MRI) or positron emission tomography (PET) images with planning computed tomography (CT) images for better tumor delineation. In stereotactic body radiotherapy for lung or liver lesions, DIR is used to propagate contours in 4-dimensional CT images to assess tumor motion and create internal target volumes.[9] Daily in-room images, such as cone-beam CT (CBCT) or mega voltage CT (MVCT) images, are usually fused with the planning CT image to verify and correct the patient setup during IGRT. Planning contours may be propagated to the CBCT or second planning CT image for adaptive treatment planning. In case of retreatment, both contours and radiation dose can be mapped from the previous treatment to the new set of planning images for optimizing the retreatment plan.

Besides CT images, MRI provides additional information related to anatomical, morphological, and functional properties of the considered structures that can be extracted using different image acquisition, such as T1- and T2-weighted (T1w- and T2w-) MRI, diffusion and perfusion MRI, and tractography. Magnetic resonance imaging may also be considered as a noninvasive marker for functional changes of PGs during and after RT treatments. In this context, although many efforts have been focused on developing and evaluating the accuracy of deformable algorithms by considering CT and/or CBCT images, few papers[10-12] dealt with DIR algorithms with MRI in the head and neck (HN) district. The latter may become increasingly important with the recent introduction of in-room MRI guidance techniques in radiotherapy.

Several DIR algorithms have been proposed in this context, including free-form deformation (FFD),[13] thin-plate spline,[14] Thirion demon,[15] and viscous fluid.[16] Among them, FFD is one of the most used approaches, due to its ability in estimating a smoothed and regular vector field and its good performances.[17] However, many different implementations of this method are available, and results in terms of accuracy may vary depending on the clinical scenario, imaging modality, and also on the chosen DIR parameters.

Nowadays, these methods are commercially available and implemented in clinical software, and thus, they are widely used in clinical practice. This kind of software provides the advantage of being easily implemented in the clinical settings as well as a reduced time required for all processes, which makes this approach very feasible for use in daily routine. Moreover, they are built to be robustly effective in a large spectrum of different applications. At the same time, there is the potential risk that the desired reduction in computational time and the not specific setup could lead to a diminished DIR accuracy. In fact, it has been shown that, in critical situations (ie, large deformation, discrepancy between images, etc), commercial software can lose some DIR accuracy.[18] Among the available commercial software solutions, MIM (v6.3; MIM Software, Cleveland, Ohio) is increasing its relevance in radiotherapy applications, and many validation studies and clinical applications for CT/CT monomodality DIR have already been reported[19,20]; however, the registration between MRI acquired before and after the treatment has not been completely evaluated yet.

On the contrary, with an open-source software, it is possible to modify the algorithm in order to achieve an acceptable DIR accuracy, even in more critical situations. Moreover, these open-source packages can be even more optimized for a specific application, by their inclusion in homemade tools. However, this type of software is not so practical for use in daily routine due to a generally higher requirement of computational time, a less user-friendly interface, and a more difficult integration with the clinical workflows. There is an important balance between time and DIR accuracy that varies for the different approaches.

The aim of this work was thus to evaluate and compare the FFD approach and parotid contour propagation process for head neck T1w-MRI of patients treated with radiotherapy with 3 different tools—the commercial MIM, the open-source Elastix software[21] in its conventional form, and an optimized version of it. In particular, this optimization concerns the contour propagation method, which was improved in a homemade version based on mesh deformation.

Materials and Methods

Patients and Imaging Data Set

Twelve patients (median age: 55 years; range: 33-79) previously treated with IGRT at our institute were selected for this study; at the time of treatment, all patients signed an informed consent form including the use of imaging information. All patients were treated with a simultaneous integrated boost approach, delivering 54 Gy to the elective lymph nodes, 61.5/66 Gy to the tumor bed/tumor, and 69 Gy to the PET-positive volume.[22] Seven oropharynx and 5 nasopharynx cases were considered; for 7 patients, the tumor was bilateral, whereas for 5, it was unilateral. Table 1 summarizes the main characteristics of the considered patients.

Table 1.

Main Patient Characteristics.

Characteristics	N
Age, years
Median (range)	55 (33-79)
Sex (M/F)	7/5
Surgery (y/n)	1/11
Chemo
Adjuvant	7 (58.3%)
Concomitant	12 (100%)
Tumor
Oropharynx	7 (58.3%)
Nasopharynx	5 (41.7%)
Stage
I-II	2 (16.7%)
III	10 (83.3%)

Abbreviations: F, female; M, male; n, no; y, yes.

Deformable registration and contour propagation were evaluated by considering the MRI acquired before treatment (MRI_1) and a MRI scan acquired 2 to 3 months after the end of the treatment (MRI_2; days after end of RT: mean = 70, range = 42-99). Magnetic resonance imaging was performed using a 1.5-T scanner (Achieva Nova; Philips Medical Systems, Best, the Netherlands) equipped with HN combined coil. The MRI protocol included a fat-suppressed 3-dimensional (3D) T1w fast field echo sequences (Repetition Time [TR] and Echo Time [TE] shortest; field of view = 200 mm × 200 mm; matrix 224 × 221; slice thickness: 2 mm; Number of Signal Averages [NSA]: 1) acquired after the intravenous injection of gadobutrol (Gadovist; Bayer Schering, Germany) at a dose of 0.1 mmol/kg body weight. On both MRI scans, an expert radiation oncologist contoured the PGs.

MIM Software

Version 6.3 of MIM was used in this work. First, the optimal setting and choice of parameters for the monomodal MRI rigid registration (RR) and DIR were investigated and evaluated in a subset of 10 patients (see Supplemental Material for a detailed description).

Rigid registration uses an intensity-based algorithm with the maximization of the mutual information index, as its optimization routine. The DIR consists of an intensity-based FFD registration algorithm,[13] in which the similarity and the smoothness criteria were combined into an energy function that is minimized. Several different parameter settings were evaluated for RR and DIR in terms of accuracy. The best results were obtained by the combination of an RR focused on a user-defined region of interest, including PGs, and a DIR on the entire image, with a deformable smoothness factor (DFS) equal to 0.5 (default value). However, no significant differences were found among the evaluated settings. In 2 patients with very high deformations, the DFS was set equal to 2. In fact, as reported in the Supplementary Material, an increase in DFS effectively improves the DIR quality in patients who present relevant nonrigid transformations. The generated deformation field was then used for contour propagation.

Elastix Software

The deformable transformation between MRI_1 and MRI_2 was estimated using the DIR method implemented in the open-source software Elastix. The chosen registration method is similar to the one implemented in MIM and consists of a rigid realignment followed by the classical FFD based on B-splines, with the parameters optimization available in the Elastix implementation. In particular, the following setup was optimized for a monomodal MRI registration—normalized MI was chosen as the similarity metric, the adaptive stochastic gradient descent as the adopted optimization algorithm, and the stopping criterion was the maximum number of iterations, which were chosen to be sufficiently high to guarantee registration convergence. The complete set of parameters is reported in Table 2.

Table 2.

Parameters Chosen for Image Registration Algorithm Implemented in Elastix.

	Rigid Registration	Elastic Registration
Similarity metric	NMI	NMI + bending energy penalty
Number of histogram bin for NMI calculation	32	32
Transformation	Euler transform (6 parameters)	Free-form deformation based on cubic B-splines
Final grid dimension	–	10 mm
Optimization algorithm	Adaptive stochastic gradient descent	Adaptive stochastic gradient descent
Maximum number of iterations	1000	5000
Number of multiresolution levels	4	5

Abbreviations: NMI, normalized mutual information

Once the image registration step was performed, contours delineated on MRI_1 were automatically propagated on MRI_2 using 2 different approaches, both starting with the deformation field estimated by Elastix.

In the first one (named DIR_Trx), the contour propagation method—as proposed in the Elastix toolbox—was adopted; in particular, the Transformix function was used to apply the deformation field estimated using MRI_2 as fixed image and MRI_1 as moving image to the binary segmentation of PGs made on MRI_1.

The second approach (named DIR_Mesh) was an optimized homemade version of the Transformix function, based on mesh deformation. This method was previously proposed and evaluated,[23,24] and it applied the deformation field (estimated using MRI_1 as fixed image and MRI_2 as moving image) to the vertices of a 3D mesh generated from the binary segmentation of PGs on MRI_1. These meshes were constructed using a wavelet-based surface reconstruction method, where, starting from manual cross-sectional contours, image voxel segmentation processing was combined with implicit surface streaming methods using wavelets. This method was proved to generate smoothed and regular surfaces, with a mesh quality comparable or even better than other state-of-art methods.[23] Once the mesh was deformed, it was cut on the corresponding planes of MRI_2 to obtain the propagated binary masks on each slice (see Figure 1). The method was implemented based on the open-source software packages Insight Segmentation and Registration Toolkit (www.itk.org)[25] and Visualization Toolkit (www.vtk.org).[26]

Figure 1.

Schematic description of the DIR_Mesh contour propagation algorithm. From the delineated contour on MRI_1, a mesh was generated (A); the estimated vector field was applied to it (B) in order to deform the vertices of the mesh (C); the obtained deformed mesh (D) was finally cut on the correspondent slice on MRI_2.

Evaluation of the Contour Propagation Methods

In order to compare the accuracy of different DIR methods, classical indices of contours’ distance and volume overlap were used. In particular, for each PG, the manually delineated contours on MRI_2 (A) were compared to the automatically deformed binary masks (B); for this evaluation, the following indices were calculated:

Average symmetric distance (ASD; mm): the average distance between 2 contours A and B, calculated as:

where C(A) and C(B) indicate the set of voxels belonging to contours A and B, respectively; d(v, C(A)) and d(v, C(B)) indicate the shortest distances between an arbitrary voxel and the contours A and B, respectively, calculated using a 3D Euclidean distance transform. Average symmetric distance is defined as symmetric, since it is first calculated from A to B, then from B to A, and finally averaged over both.[27]

Maximum symmetric distance (MSD; mm): the maximum distance (or Hausdorff distance) between 2 contours A and B, calculated as:

Dice similarity coefficient (DSC): the most widely used index for measuring the agreement between 2 volumes[28], calculated as:

Dice similarity coefficient ranges from 0 to 1, where 1 indicated a perfect overlap.

Sensitivity index (Sens): it reflects the probability that the automatically propagated contours match the manual contour[29]; it is defined as:

Inclusiveness index (Incl): a surrogate of the specificity and reflects the probability that a voxel of the propagated mask is really a voxel of the manual mask[29]; it is defined as:

From these 2 last indices, a modified version of the receiver operating characteristic (ROC) analysis was done, as proposed in La Macchia et al,[29] by plotting Sens versus 1-Incl.

The quantitative comparative analysis was performed using a homemade software implemented in Matlab (Release 8.5.0; The Mathworks Inc, Natick, Massachusetts).

A nonparametric Kruskal-Wallis test was performed for each calculated comparative index to determine whether significant differences were observed between the 3 considered methods, and a Dunn correction was eventually used to individuate which groups were significantly different.

Results

Qualitative results of contour propagation using MIM, DIR_Trx, and DIR_Mesh are presented in Figure 2, showing a good agreement between the manual reference contours and the propagated ones, without evident differences among the 3 methods. Quantitative results of contours evaluation are reported in Table 3. Average symmetric distance, MSD, and DSC indices calculated on both PGs for all patients and for all 3 considered DIR methods are plotted in Figure 3. DIR_Mesh presented the lowest average values of ASD and MSD and the highest average value of DSC over the population among the 3 methods, although nonsignificant (P values of the nonparametric Kruskal-Wallis test = 0.57, 0.40, and 0.27 for ASD, MSD, and DSC, respectively). The average propagation errors with respect to the reference contours are lower than the voxel diagonal (2 mm), and DSC is around 0.8, generally considered as the threshold for good agreement.

Figure 2.

Examples of contour propagation using MIM (left), DIR_Trx (center), and DIR_Mesh (right) methods. Manual reference contours are also reported in red.

Table 3.

Results of Contour Evaluation for the 3 Methods.a

	ASD, mm	MSD, mm	DSC	Sens	Incl
MIM	1.83 (1.03)	12.40 (5.37)	0.76 (0.10)	0.87 (0.07)	0.70 (0.14)
DIR_Trx	1.52 (0.56)	11.02 (4.25)	0.79 (0.06)	0.88 (0.06)	0.73 (0.09)
DIR_Mesh	1.46 (0.50)	10.00 (3.11)	0.81 (0.05)	0.85 (0.07)	0.78 (0.09)

Abbreviations: ASD, average symmetric distance; DSC, Dice similarity coefficient; Incl, inclusiveness index; MSD, maximum symmetric distance; Sens, sensitivity index.

aValues are reported as mean (standard deviation) within the considered population.

Figure 3.

Average symmetric distance (ASD), maximum symmetric distance (MSD), and Dice similarity coefficient (DSC) indices calculated for each patient and each parotid with the 3 different DIR methods (MIM, DIR_Trx, DIR_Mesh). LP indicates left parotid; RP, right parotid.

In Figure 4, the modified ROC curve is shown; as already reported in Table 3, DIR_Mesh performed slightly better than the other 2 methods, having an area under the curve higher than the others, although nonsignificant (P value = .26).

Figure 4.

Modified receiver operating characteristic (ROC) curve for the 3 methods.

In few cases, contour propagation failed, in particular when MIM was considered; this happened when big deformations and changes in parotid shape occurred (see Figure 5A and B, where a clear difference in head and parotid shape is shown). In this case, MIM had a greater difficulty in recovering parotid contours (green line in Figure 5C), with respect to the other 2 methods (light blue and purple line in Figure 5C).

Figure 5.

Example of a failed contour propagation, especially for MIM (green line). A, Magnetic resonance imaging (MRI) acquired before radiation therapy (RT) with manual parotid contours. B, Magnetic resonance imaging acquired after RT with manual parotid contours. C, Magnetic resonance imaging acquired after RT with manual and deformed contours.

Discussion

Three different implementations of the FFD method for the MRI–MRI registration were qualitatively and quantitatively evaluated in the HN region, in terms of PGs contour propagation. The 3 methods resulted to be nearly comparable to each other, with a trend of the optimized homemade approach in better recovering strong deformations.

In general, the estimated DIR registration accuracy was in agreement with the results already published in the literature in the HN district. However, quantitative comparisons to other works are quite difficult due to the several DIR algorithms used and the different metrics and structures considered. Few papers[10,30,31] evaluated multimodal and multitemporal DIR accuracy by considering MRI for HN region—a Dice coefficient between 0.6 and 0.9 is reported on PGs and on other structures by using open-source and homemade DIR approaches. Regarding FFD approach, similar accuracy was reported by considering monomodality (CT–CT)[32] or quasi monomodality (CT–CBCT) studies,[33] by using both open-source DIR approaches and commercial software. In these studies, a mean Dice value around 0.8 and a distance between contours of 2 mm are generally reported by considering soft tissues (as PGs), which are values similar to those found in our work.

Good results for DIR registration were also reported in studies where the performance of DIR algorithms was investigated in physical or virtual head–neck phantoms with known deformations,[20,34,35] especially by considering monomodality (CT–CT, MVCT–MVCT) or quasi monomodality (CT–CBCT) image acquisition. In addition to Dice analysis, the difference between the applied deformation vector field and the ground truth was evaluated voxel-by-voxel within the image set. A mean difference within 4 mm was reported for monomodality image registration; larger differences (up to 40 mm) were reported for different imaging modalities (as CT–MVCT) or in case of relevant neck deformations. In our study, we have found that values of mean distance are similar to the voxel diagonal, which is a result even better than those reported in these mentioned studies.

Looking at the differences among the 3 considered FFD implementations, when small deformations were present, they showed similar accuracy. On the contrary, MIM presented a worsening in DIR accuracy for patients with large deformations and discrepancy between images. In fact, as already mentioned, commercial packages are generally built to reach acceptable performance very fast and in many different situations with an easy and user-friendly parameter optimization. However, the setup, although possible, is limited and the complete implementation is unknown, as they are “black-box” tools; therefore, the accuracy in a specific application, as the MRI–MRI parotid registration, may be reduced in critical cases. For these specific patients, the optimized Elastix algorithm showed a slightly improved performance with better results in terms of parotids propagation, thanks to the integration of a homemade implementation in an open-source code. This suggests that the contour propagation method has a strong impact on the final accuracy, together with the image registration algorithm. In fact, the choice of using a mesh-based deformation, combined with a mesh generation algorithm, which warrants regular and smoothed surfaces, has highlighted its ability in correctly deforming parotid contours, even in challenging cases.

In this work, we have found that the evaluated DIR methods were accurate for the specific registration of T1w-MRI acquired before and after RT. In this context, the use of MRI in the clinical practice seems to be very promising. In fact, there is a growing interest in further introducing multimodal and functional MRI into the patient’s treatment process. The additional and complementary information can be potentially used to improve and personalize the treatment strategies, by considering the patient’s biological characteristics. For example, target delineation, dose-painting applications, treatment strategies decision based on biological patient’s characteristics, outcome/toxicity predictive model assessment can take advantage by MRI. Moreover, MRI could provide a huge amount of information that can be extracted and used to correlate with the received physical dose and with the biological properties of the pretherapy tissue, in order to characterize and develop predictive models of treatment outcome.[36] Deformable image registration and contour propagation between MRI acquired before treatment start, during RT, and at follow-up can be widely used to evaluate the patient-specific treatment response, by considering morphological and functional properties evaluated on these images. For example, for salivary gland toxicity, Dirix et al and Zhang et al[37,38] showed that a variation of the apparent diffusion coefficient (ADC) value in diffusion-weighted MRI, measured before and after RT, is significantly correlated with the incidence of xerostomia. The ADC changes resulted also correlated with changes in saliva stimulation measurements, gustatory scores, and scintigraphy examination.

The major limitation of this study was the small number of patients, which limited the statistical power of the analysis. For the DIR validation, we have not performed any evaluation about vector field consistency and feasibility; however, we are confident that the adopted approach, based on FFD, B-splines, and a regularization term, can ensure a smoothed and regular deformation field without folding, at least in the volume of interest. Finally, we have not evaluated the intra- and interobserver variability in contouring PGs, due to the very long time needed on these high-resolution images. However, based on previous works reporting a DSC of 0.8,[39] a conformity index of 0.66[40] for the interobserver contouring of parotids in CT images, and an interobserver parotid volume variation of about 5% in MRI,[41] the accuracy found here can be considered acceptable.

Conclusion

Based on the results reported in this study, we may state that a commercial clinical software such as MIM can reach adequate MRI–MRI registration accuracy. In fact, it combines the advantage of time saving (due to the fast DIR registration and contour propagation algorithm) with an easy, though reduced, setup and a user-friendly interface, which, most of the time, is already integrated in the clinical routine. For critical situations (ie, large deformations, discrepancy between images, significant noise), a customized software, such as the presented modified version of Elastix software, should be preferred. With this kind of approach, the increase in the required time and the major difficulty to find the optimal configuration and to introduce it in the clinical setting are counterbalanced by a gain in DIR and contour propagation accuracy.