LOCALIZED BLADDER DOSE ACCUMULATION IN MULTI-FRACTION CERVICAL CANCER BRACHYTHERAPY by Roja Zakariaee Kouchaksaraee B.Sc., Sharif University of Technology, 2004 M.Sc., Simon Fraser University, 2009 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Physics) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) December 2016 © Roja Zakariaee Kouchaksaraee, 2016 ii Abstract Radiation therapy in the definitive treatment of locally advanced carcinoma of the cervix consists of external beam radiotherapy (EBRT) combined with image-guided high-dose-rate (IG HDR) intracavitary brachytherapy (ICBT). IG HDR-ICBT is a relatively new, advanced form of brachytherapy treatment planning and delivery that still relies largely on dose criteria based either on clinical experience using older techniques, or on limited data. Calculation of cumulative dose received over multiple treatment fractions currently utilizes dose-volume-histograms (DVH), which do not provide information about the spatial distribution of dose within the structures of interest. Since most of the organs at risk (OAR) for this site, like the bladder and rectum, are highly deformable, DVH data summed over multiple treatment fractions do not provide accurate estimates of the cumulative dose to specific regions of the organ, and therefore may not be the most appropriate metric to use in treatment planning and dose assessment. The primary goal of this thesis was, therefore, to develop methods to more accurately quantify “locally accumulated” dose to the bladder-wall in multi-fraction IG HDR-ICBT for cervical cancer using deformable registration, and to apply these to a study of locally-accumulated dosimetric parameters as predictors of late urinary toxicity. To this end, different deformable image and point-set registration methods were evaluated, using phantom data, for their ability to register bladder-wall contours of various sizes. The best of these was retrospectively applied to point-sets representing the bladder-wall in multiple HDR treatment fractions for 60 cervical cancer patients treated at the BC Cancer Agency. The transformation maps obtained from the registrations were used to calculate cumulative dose parameters for the bladder-wall and the urethral opening to the bladder (bladder neck). Patients were divided into Case and Control groups based on urinary iii toxicity scores for different symptoms, and the ability of the cumulative parameters to predict toxicity was evaluated. It was concluded that some locally-cumulative parameters estimated using our method have improved potential for predicting urinary toxicity, as compared to traditional DVH-based parameters, in our data set. In addition, dose to small volumes around the bladder neck was found to be a predictor of incontinence. iv Preface A version of Chapter 3 has been published (Zakariaee, R., Hamarneh, G., Brown, C. J., & Spadinger, I. Validation of non-rigid point-set registration methods using a porcine bladder pelvic phantom. Phys. Med. Biol. 2016; 61(2):825-54. doi:10.1088/0031-9155/61/2/825). The porcine bladder was acquired from a pig that was used for teaching purposes for traumatic injuries at the Jack Bell Research Centre, approved by UBC Animal Care Committee for the Basic Surgical Techniques (BST) project with Protocol # A16-0029. I. Spadinger provided advice and guidance on the study design. C. Brown helped with a dry run of the phantom assembly and scan, performed identification of the landmarks as the second independent observer, and assisted in selecting and implementing some of the point-set registration algorithms. G. Hamarneh audited the selected point-set registration methods, suggested some of the methods and also provided technical insight into the subject of image and point-set registration. I conducted phantom and landmark preparation and assembly, carried out the experiments and measurements, implemented the registration methods, performed the registrations and completed the analyses. I wrote most of the manuscript but all the authors actively participated in editing and revising the manuscript for publication. A version of Chapter 4 has been accepted for publication (Zakariaee, R., Hamarneh, G., Brown, C. J., Gaudet, M., Aquino-Parsons, C. & Spadinger, I. Bladder Accumulated Dose in Image-Guided High-Dose-Rate Brachytherapy for Locally Advanced Cervical Cancer and its Relation to Urinary Toxicity. Phys. Med. Biol. 2016; 61(24): 8395–8411. doi:10.1088/0031-9155/61/24/8408). For the study in Chapter 4, I conducted all the bladder delineations, landmark identifications, updating the LENT-SOMA toxicity scores, completing registrations, and calculating and analyzing all the data and results. The MATLAB codes for implementing the v deformable registration, dose accumulation, and calculation of the dose volumes were constructed by C. Brown and myself. I. Spadinger provided advice and guidance on the study design and interpretation of the results. G. Hamarneh oversaw the technical aspects of the registration process. The initial LENT-SOMA urinary toxicity data were gathered by M. Gaudet and C. Aquino-Parsons, who also provided clinical guidance. In Chapter 5, the analysis of the individual urinary endpoint sub-scores and their relation with dose-volume parameters was devised by me. Analytical guidance was provided by I. Spadinger, while clinical guidance was provided by C. Aquino-Parsons and P. Lim at the BC Cancer Agency. The work reported in Chapters 4 and 5 in this dissertation was approved by the University of British Columbia Research Ethics Boards, certificate number H1202320. vi Table of Contents Abstract .......................................................................................................................................... ii Preface ........................................................................................................................................... iv Table of Contents ......................................................................................................................... vi List of Tables .............................................................................................................................. xiii List of Figures ...............................................................................................................................xv List of Abbreviations ............................................................................................................... xxiv Acknowledgements ................................................................................................................. xxvii Dedication ................................................................................................................................. xxix Chapter 1: Introduction ................................................................................................................1 1.1 Cervical Cancer ............................................................................................................... 1 1.2 Radiation Therapy ........................................................................................................... 2 1.2.1 Overview and Historical Perspective .......................................................................... 2 1.2.2 Physical Basis of Radiation Therapy .......................................................................... 4 1.2.2.1 Ionizing Radiation ............................................................................................... 4 1.2.2.2 Electron Interactions ........................................................................................... 5 1.2.2.3 Photon Interactions ............................................................................................. 7 vii 1.2.2.4 Radiation Intensity and Geometric Effects ....................................................... 12 1.2.2.5 Radiation Dose .................................................................................................. 13 1.2.2.6 Sources of Therapeutic Photons and Electrons ................................................. 14 1.2.2.6.1 Teletherapy (External Beam) Photons and Electrons ................................. 14 1.2.2.6.2 External Beam Field Shaping and Intensity Modulation ............................ 15 1.2.2.6.3 Brachytherapy Sources ................................................................................ 15 1.2.2.6.4 Brachytherapy Delivery Systems and Methods .......................................... 19 1.2.2.7 Patient Dose Calculations ................................................................................. 21 1.2.2.7.1 External Beam Dose Calculations ............................................................... 23 1.2.2.7.2 Brachytherapy Dose Calculations ............................................................... 24 1.2.3 The Radiobiology of Cancer Treatment ................................................................... 26 1.2.3.1 Mechanism of Proliferative Cell Death ............................................................ 27 1.2.3.2 Radiation Damage to Normal Tissue ................................................................ 27 1.2.3.3 The Linear Quadratic Cell Survival Model ...................................................... 28 1.2.3.4 Fractionation and Therapeutic Ratio ................................................................. 29 1.2.3.5 Modeling Biological Effect............................................................................... 31 1.2.4 Quantifying the Side Effects of Radiation Treatment ............................................... 32 1.2.5 Commonly Used Terms and Parameters in Radiotherapy Treatment Planning ....... 34 viii 1.3 Image-guided Intracavitary Brachytherapy for Locally Advanced Carcinoma of the Cervix ...................................................................................................................................... 37 Chapter 2: Medical Image Registration and Shape Correspondence.....................................43 2.1 Image Registration ........................................................................................................ 43 2.1.1 Transformation Models ............................................................................................. 45 2.1.1.1 Rigid-body Transformation .............................................................................. 45 2.1.1.2 Affine Transformation ...................................................................................... 46 2.1.1.3 Deformable Transformation.............................................................................. 47 2.1.1.3.1 Linear Elastic Model ................................................................................... 48 2.1.1.3.2 Fluid Models ............................................................................................... 49 2.1.1.3.3 Radial Basis Functions ................................................................................ 49 2.1.2 Similarity Metrics ..................................................................................................... 52 2.1.2.1 Landmark-Based Registration .......................................................................... 52 2.1.2.2 Surface-Based Registration ............................................................................... 53 2.1.2.3 Intensity-Based Registration ............................................................................. 54 2.1.3 Optimization Process ................................................................................................ 56 2.2 Shape Correspondence .................................................................................................. 57 2.3 Validation of Registration Outcomes............................................................................ 59 2.3.1 Visual Inspection ...................................................................................................... 59 ix 2.3.2 Numeric Validation Methods .................................................................................... 60 2.3.2.1 Target Registration Error .................................................................................. 61 2.3.2.2 Dice Similarity Coefficient ............................................................................... 62 2.3.2.3 Hausdorff Distance ........................................................................................... 63 2.3.2.4 Root-Mean-Square Error .................................................................................. 63 2.4 Project Overview .......................................................................................................... 63 Chapter 3: Validation of Non-Rigid Point-Set Registration Methods ....................................66 3.1 Introduction ................................................................................................................... 66 3.2 Point-Set Registration Methods .................................................................................... 71 3.2.1 Iterative Closest Point ............................................................................................... 72 3.2.2 Thin-Plate-Spline Robust Point Matching ................................................................ 73 3.2.3 Shape Context ........................................................................................................... 74 3.2.4 Coherent Point Drift .................................................................................................. 75 3.2.5 GMM Registration .................................................................................................... 76 3.3 Materials and Methods .................................................................................................. 76 3.3.1 Phantom Preparation ................................................................................................. 77 3.3.2 CT Scan ..................................................................................................................... 77 3.3.3 Point-Set Registration ............................................................................................... 81 x 3.3.4 MIM DIR .................................................................................................................. 86 3.3.5 Registration Accuracy Measures .............................................................................. 87 3.3.5.1 TREavg and NE .................................................................................................. 87 3.3.5.2 TREtangential and TREperpendicular .......................................................................... 88 3.3.5.3 RMS and HD..................................................................................................... 89 3.4 Results ........................................................................................................................... 90 3.4.1 Landmark Identification............................................................................................ 91 3.4.2 Robustness and Runtime ........................................................................................... 93 3.4.3 Registration Accuracy: TRE ..................................................................................... 94 3.4.4 Registration Accuracy: Effect of Volume Difference .............................................. 99 3.4.5 Registration Accuracy: RMS and HD ..................................................................... 101 3.4.6 Subsampling Effect ................................................................................................. 103 3.4.7 MIM DIR Outcomes ............................................................................................... 105 3.5 Discussion ................................................................................................................... 105 3.6 Conclusions ................................................................................................................. 114 Chapter 4: Bladder Accumulated Dose and Urinary Toxicity ..............................................117 4.1 Introduction ................................................................................................................. 117 4.1.1 Limitations of DVH Parameters ............................................................................. 118 xi 4.1.2 Practice of Varying Bladder Filling Status at BCCA ............................................. 118 4.1.3 Genitourinary Toxicity and Radiation Dose ........................................................... 119 4.1.4 Chapter Overview ................................................................................................... 121 4.1.5 Previous Work ........................................................................................................ 122 4.2 Materials and Methods ................................................................................................ 123 4.2.1 Patient Materials and Treatment Characteristics .................................................... 123 4.2.2 Urinary Toxicity Follow-up and Scoring ................................................................ 125 4.2.3 Contour Processing ................................................................................................. 126 4.2.4 CPD Deformable Point-Set Registration ................................................................ 128 4.2.5 Dose Accumulation ................................................................................................. 129 4.2.6 Registered Dosimetric and Volumetric Parameters ................................................ 130 4.2.7 Validation Tests ...................................................................................................... 133 4.2.8 Volume Variation Effect ......................................................................................... 135 4.3 Results ......................................................................................................................... 136 4.3.1 Registration and Dose Accumulation ..................................................................... 136 4.3.2 Contiguous and Non-Contiguous Volume Calculation .......................................... 137 4.3.3 Urinary Toxicity, Age, and EBRT dose.................................................................. 139 4.3.4 Registered Dose Parameters ................................................................................... 140 xii 4.3.5 Volume Variation.................................................................................................... 142 4.4 Discussion ................................................................................................................... 143 4.5 Conclusion .................................................................................................................. 152 Chapter 5: Dose Response Analysis of Specific Urinary Symptoms .....................................153 5.1 Introduction ................................................................................................................. 153 5.2 Materials and Methods ................................................................................................ 155 5.3 Results ......................................................................................................................... 159 5.4 Discussion ................................................................................................................... 165 5.5 Conclusion .................................................................................................................. 168 Chapter 6: Conclusion and Future Work ................................................................................170 6.1 Conclusions ................................................................................................................. 170 6.2 Future Work ................................................................................................................ 177 Bibliography ...............................................................................................................................179 Appendices ..................................................................................................................................207 Appendix A : The LENT-SOMA questionnaire used at BC Cancer Agency for collecting late GU and GI toxicity data for cervical cancer patients. ............................................................. 207 Appendix B : Boxplots of frequency, dysuria, nocturia, and stream, illustrating the descriptive statistics for the registered dosimetric, volumetric, and bladder-neck-dose parameters for the Case and Control in G2/G3-cutoff Groupings. ....................................................................... 213 xiii List of Tables Table 1.1 Subjective, objective, Management, and Analytic (SOMA) criteria for bladder and urethra late effects, based on [28]. ................................................................................................ 34 Table 1.2 Generic LENT-SOMA scoring system, based on [28]. ................................................ 34 Table 1.3 List of the treatment planning parameters and constraints, used for HDR ICBT for cervical cancer at the BC Cancer Agency. .................................................................................... 42 Table 3.1 Summary of the components of each point-set registration technique evaluated in this study. ............................................................................................................................................. 72 Table 3.2 The software source and available output type/s for each registration method. ‘Correspondence’ indicates the correspondence vector output (C), and ‘deformed points’ indicates the transformed point-set matrix (M) output. ............................................................................... 83 Table 3.3 The registration parameters, trial and optimal values for the MATLAB implementation of the point-set registration algorithms evaluated in the study. .................................................... 84 Table 3.4 The centroid coordinates of the contour point-sets for different bladder volumes. ...... 86 Table 3.5 The TRE values and the inter-observer and intra-observer fiducial localization errors (FLE) for all landmarks, in millimeters. ....................................................................................... 93 Table 3.6 Average TRE values for the fiducial landmarks categorized based on their anatomic positions (smallest bladder used as target). ................................................................................. 100 Table 4.1 Bladder wall dosimetric and volumetric parameters and their definitions listed at the top xiv and landmark-based parameters at the bottom. ........................................................................... 134 Table 4.2 The differences between MATLAB-based and TPS-based DVH values (in Gy and as a percentage of the prescription dose) for 1cc and 2cc volumes, across all fractions of all subjects...................................................................................................................................................... 136 Table 4.3 Descriptive statistics for contiguous rD2cc values based on three different regularization regimens used to calculate contiguous volumes. All units are in Gy. ........................................ 139 Table 4.4 The Case (CS) and Control (CN) statistics for different dosimetric and volumetric parameters based on registered dose, and the bladder volume variation in percentage. The p-values less than 0.1 are highlighted in bold text. ................................................................................... 145 Table 4.5 The Case and Control statistics for the bladder neck and left and right ureters’ orifices...................................................................................................................................................... 145 Table 5.1 The p-values and statistical power (1-β) values based on Student t-test or Mann-Whitney U-test for the differences in the means of Case and Control groups for incontinence, urgency, and hematuria, based on G2 and G3/G1-cutoff groupings, for all the registered parameters. .......... 166 xv List of Figures Figure 1.1 Illustration of the Bremsstrahlung process (a) and its prevalence (b). .......................... 7 Figure 1.2 The mass attenuation coefficient as a function of photon energy, for coherent scattering, incoherent (Compton) scattering, photoelectric effect, pair production, for water and lead. The total mass attenuation coefficient (sum of all coefficients) is also shown. Data obtained from [8]. ...... 9 Figure 1.3 Schematic diagrams of photon interactions with atoms, illustrating (a) coherent scattering, (b) photoelectric effect, (c) incoherent (Compton) scattering, and (d) pair production........................................................................................................................................................ 11 Figure 1.4 Illustration of an isotropic point source emitting photons. .......................................... 13 Figure 1.5 The dose distribution of a 4-field box treatment plan for cancer of cervix, shown as a colour wash. The bladder (yellow), rectum (green) and clinical target volume (red) contours can also be seen. .................................................................................................................................. 16 Figure 1.6 The GammaMed high dose rate Iridium-192 source. .................................................. 19 Figure 1.7 Diagram of the source drive mechanism of the Nucletron microSelectron HDR remote-control afterloading unit. The radiation source is permanently fixed to a cable and is driven to a number of different “dwell” positions within the applicator using a precision stepping motor. A second cable with a dummy source is used to check the integrity of the applicator tube before attempting to drive the active source through the tube. Each guide tube in the intracavitary brachytherapy applicator connects to one of the 18 channels....................................................... 22 xvi Figure 1.8 Top: Examples of tandem-ring (Left) and tandem-ovoid (Right) intracavitary applicator used for cervical cancer HDR brachytherapy. The radiation sources are guided through the tandem to irradiate the paracervical region and separately through the ring/ovoids to irradiate the cervix. The rectal paddle can be used optionally to push the rectum away from the radiation sources. Bottom: The tandem-ring applicator in-situ. The tandem lies inside the uterus and the ring sits at the cervix. ...................................................................................................................................... 23 Figure 1.9 The geometry of a brachytherapy cylindrical source, as described in the TG-43 formalism. ..................................................................................................................................... 25 Figure 1.10 Representative (not quantitative) survival curves for tumour and early-responding tissues (larger 𝜶/𝜷 ratio) and that of late-responding tissues (smaller 𝜶/𝜷 ratio), demonstrating higher therapeutic ratio at smaller doses compared to larger doses, in a single-hit delivery. ...... 29 Figure 1.11 Surviving fraction curves for tumour (solid blue) and normal tissue (solid black) cell lines, splitting a total dose of 12 Gy into (a) two dose fractions of 6 Gy each and (b) six fractions of 2 Gy each. The solid lines depict the unfractionated linear-quadratic (LQ) survival curves for the corresponding cells. The surviving fractions for each dose fraction are shown as dots connected by the dotted curves, which represent the initial portion of the LQ curve applied to each fraction, for each cell line. ........................................................................................................................... 31 Figure 1.12 Isodose distribution within cervical HDR brachytherapy treatment volume (CTV contour shown in orange).............................................................................................................. 37 Figure 1.13 Cumulative DVH plot for EBRT treatment of cervical cancer, showing the dose distributions for CTV, bladder, rectum, right and left kidneys, and spinal canal ROIs. The volumes xvii are expressed in terms of percentages of the total structure volumes and the total structure volumes are shown in the legend. The dose values are shown in absolute dose values in the top x-axis and in percentages of the prescribed dose (45 Gy) in the bottom x-axis. ............................................ 38 Figure 1.14 Sample HDR planning CT with labelled OAR contoured structures, showing bladder, sigmoid and rectum. The Isodose curves are also displayed. ....................................................... 39 Figure 1.15 Left: The ICRU point A and point B locations, defined based on the applicator alignment. Right: The bladder and rectum reference points. ........................................................ 40 Figure 2.1 Schematic of rigid, affine and non-rigid (deformable) transformations in 2D. ........... 46 Figure 2.2 A 2D illustration of different “hat and head” solutions for a feature-based registration, yielding similar RMSD distances, while only one is a perfect match. ......................................... 54 Figure 2.3 Example of possible non-anatomical mismatch between two brachytherapy planning CT images, showing air pocket in bladder (yellow contours) and rectum (orange contours) present only in one image. The brachytherapy applicator in the left (tandem-ring - green and magenta markings) is also different from the one in the right image (tandem-ovoids - green markings). In addition, a high degree of deformation in bladder is present close to the non-moving bony anatomy (cross section of the pelvic bones as shown). ............................................................................... 59 Figure 2.4 Sample of 2D deformable point-set registration, matching the blue fish to the red fish using the coherent point drift (CPD) algorithm (adapted from [94]). Both point-sets have missing points and spurious outlier points. ................................................................................................ 60 Figure 3.1 The study design flowchart showing the main data acquisition and analysis steps. ... 78 Figure 3.2 The pig bladder with urethra and ureters attached. The yellow rectangular fiducials are xviii rubber bars and the rest are plastic beads. Different groupings of each fiducial type were used. 79 Figure 3.3 Pictures of the pelvic phantom model. The hand-made plasticine uterus model is visible in the bottom image. The Foley-catheter and the CT-safe plastic clamp were used to fill the bladder with the contrast-water mixture and secure the liquid inside. All visible fiducials are circled in red in both pictures. ............................................................................................................................. 80 Figure 3.4 (a) The axial (top left), coronal (bottom left), and sagittal (bottom right) planes of the CT scan of the pelvis phantom, with the 180 mL bladder volume filling. The unlabeled structures are the water-filled balloons and bags. (b) The appearance of a plastic bead fiducial and two rubber fiducials on the CT image. While the two fiducial types are completely distinguishable, the rubber bars are more distinct because of their shape and high HU number. ............................................ 82 Figure 3.5 Sample modified CT images shown against the raw images. The modified and raw landmark areas are indicated by red arrows. The HU value of the uterus model region (denoted by green arrows) was also altered to a value closer to that of soft tissue. ......................................... 89 Figure 3.6 A frame of reference is created on the tangent plane to the surface of the bladder at landmark P on the reference bladder. Using the corresponding point in the other deformed point-set, P, one can find the tangential (orange) and perpendicular (pink) distances, i.e. TREtangential and TREperpendicular, by using the normal vector of the tangent plane. .................................................. 90 Figure 3.7 The fiducial landmarks shown on the surface contour point-cloud of the 360 mL bladder. The Anterior, Posterior, Superior, and Inferior directions are shown as A, P, S, and I, respectively. The landmark points facing the front of the view (anatomic Right) are shown in dark green and the ones on the back surface (anatomic Left), in light green. Note that pl-RPS is behind xix a fold on the surface that is part of a dent in the Superior–Posterior zone from the uterus indentation..................................................................................................................................... 92 Figure 3.8 The change in the registration runtime with sample size. The plot is in log scale for better representation of the variations. TPS-RPM was not tested with sample sizes larger than 4 K due to its very long runtime. ......................................................................................................... 95 Figure 3.9 (a) The point-clouds of the target (90 mL, in this example) point-set, in blue, and moving (360 mL), in red, before registration. ((b)–(e)) The target and deformed point-clouds after registration with (b) CPD, (c) GMM-reg, (d) TPS-RPM, and (e) ICP-finite. Note that the scale of the plot in (a) is different from the rest of the plots, in order to better see the details in the shape of the point-clouds. (f) The target and moving point-clouds matched by SC, showing the correspondences for 1/tenth of the points, as green vectors. ........................................................ 96 Figure 3.10 Registration failure of ICP-finite method, showing the deformed 180 mL point-set (blue) and the 90 mL reference structure (red). ............................................................................ 97 Figure 3.11 Comparison of the registered landmark locations across different registration methods. The reference landmarks of the 90 mL volume are depicted by red stars and are numbered in the axial, coronal, and sagittal views of the reference point cloud. The corresponding average location, across the three moving structures, computed using five registration methods are shown in different colours. The projection of the TREavg vectors for each method are also displayed in each view to facilitate comparison. ....................................................................................................... 98 Figure 3.12 The vector diagram representing TREavg (arrows), TREtangential (horizontal dashed lines), and TREperpendicular (vertical dashed lines) for all five point-set registration programs. The xx circular traces are shown for aiding the comparison of the TREavg among different methods. .. 101 Figure 3.13 (a) The bar graph of TREavg values for different registration methods, separated by the moving structure volumes (with the 90 mL structure used as the reference). The error bars in both plots represent the SD over different landmarks and do not necessarily represent the uncertainty in the TRE values. (b) The bar graph of normalized error (NE) values for different registration methods, separated by the reference structure volumes. The mean values over all reference volumes are shown above the bars for each method................................................................... 102 Figure 3.14 The bar graphs of the RMS error and HD for different registration methods, except SC, separated by the reference structure. Mean and SD values are shown above the bars for each method......................................................................................................................................... 104 Figure 3.15 The mean and SD (error bars) of TREavg for different sample sizes. TPS-RPM was not included this study due to its very slow processing. The ICP-finite data was excluded from this graph because it showed almost no variation in both mean and SD values and also it was out of the range of other data................................................................................................................. 106 Figure 3.16 Examples of the MIM DIR outcome, for registering the 480 mL raw CT image to that of 90 mL, showing the axial (left), sagittal (middle) and coronal (right) views of the deformed image. MIM’s ‘RegReveal’ tool is used, in three different positions on the bladder surface, to reveal the achieved alignment of the moving and target images inside the smaller inner box. .. 107 Figure 4.1 Cross section of a bladder-wall cumulative dosage mask, showing the distribution of voxel doses as different shadings in the mask, with the brighter regions having a larger registered dose than the darker ones. ........................................................................................................... 130 xxi Figure 4.2 Example of the registration outcomes for registering the reference bladder-wall contour point-set (FR4, shown in red) to two larger fractions (FR2 and FR3, shown in blue): ‘Before’ and ‘After’ scenarios.......................................................................................................................... 138 Figure 4.3 The per-fraction and accumulated point-dose distributions, on the same geometrical scale, for the representative case in Fig. 4.2. The fraction doses are locally accumulated on the reference fraction (FR4), showing the variations in the locations of hotspots, among the fractions and also the accumulated result. ................................................................................................. 140 Figure 4.4 Sample hotspot regions on the bladder wall for (a) tightly regularized contiguous 2-cc volume, (b) loosely regularized contiguous 2-cc volume, (c) non-regularized contiguous 2-cc volume, (d) non-contiguous 2-cc volume, and (e) non-contiguous 5-cc volume. (f) Corresponding dose distribution. ......................................................................................................................... 141 Figure 4.5 The scatter plot of toxicity score versus volume variation in terms of VolVar=(Max-Min/Min) of fraction volumes for the sixty subjects in the study, also showing the linear regression fit. Note that VolVar is in absolute values and not in percentage.............................................. 144 Figure 4.6 Schematic showing a case with 3 hot spots (A, B, and C): region A represents a contiguous 2cc region enclosed by the 28 Gy isodose line, whereas regions B and C enclosed by the 28 Gy isodose line are slightly smaller than A (therefore the contiguous D2cc for these two regions would be slightly lower than 28 Gy and would not be selected as the representative contiguous D2cc). The non-contiguous D5cc has contributions from all three hot spots, and is represented by the dashed red 29 Gy line. .................................................................................. 151 Figure 5.1 Bladder neck location determined using the location of the Foley catheter on three xxii orthogonal views and the bladder wall contour. Note that since the orthogonal views are rotated from the original planes, the bladder wall contour looks jagged. ............................................... 158 Figure 5.2 The toxicity grade distributions, separated by endpoints (top) and grade levels (bottom)...................................................................................................................................................... 160 Figure 5.3 Boxplot illustrating values of the different registered dosimetric and volumetric parameters (in Gy and cc, respectively) for Case and Control groups defined on the basis of Incontinence toxicity grades using G2 and G3-cutoffs. Red lines denote the median for each distribution, the top and bottom edges of the boxes indicate the 25th (𝒒𝟏) and 75th (𝒒𝟑) percentiles, the whiskers represent the minimum and maximum values excluding outliers, and the outliers (red crosses) are values > 𝒒𝟑 + 𝟏. 𝟓(𝒒𝟑 − 𝒒𝟏) or < 𝒒𝟏 − 𝟏. 𝟓(𝒒𝟑 − 𝒒𝟏). The integral dose values (I_dose) are divided by 100 in order to fit within the scale of the plot. .................. 162 Figure 5.4 The boxplot for different registered dosimetric and volumetric parameters (in Gy and cc, respectively), separated by Case and Control groups based on the Urgency toxicity grades using G2 and G3-cutoffs. The meaning of the symbols is the same as in Fig. 5.2. The integral dose values (I_dose) are divided by 100 in order to fit within the scale of the plot. ...................................... 163 Figure 5.5 The boxplot for different registered dosimetric and volumetric parameters (in Gy and cc, respectively), separated by Case and Control groups based on the Hematuria toxicity grades using G2 and G1-cutoffs. The meaning of the symbols is the same as in Fig. 5.2. The integral dose values (I_dose) are divided by 100 in order to fit within the scale of the plot. .......................... 164 Figure 6.1 A mock bladder dose distribution created in VariSeed, showing two separate high-dose regions. The larger maroon region represents the contiguous D2cc, and all doses in the maroon xxiii regions are > D2cc. ....................................................................................................................... 174 Figure 6.2 Comparison of the DVH curves for the larger and smaller contiguous volumes with the traditional non-contiguous DVH curve. The bottom graph is a close-up of the top graph. ........ 175 xxiv List of Abbreviations (including page numbers where the terms are defined in the body of the thesis) 2D two-dimensional 3D three-dimensional BCCA British Columbia Cancer Agency (page 39) BED biologically equivalent dose (page 30) BT brachytherapy (page 2) cc cubic centimeters (commonly used in lieu of cm3 or milliliters in the brachytherapy literature) CC Correlation coefficient (page 52) CPD coherent point drift (pages 56, 67, 73) CT computed tomography (page 33) CTCAE common Terminology Criteria for Adverse Events (page 32) CTV clinical target volume (page 34) DICOM Digital Imaging and Communications in Medicine DIR deformable image registration (page 46) DNA deoxyribonucleic acid DSB double strand break (page 27) DSC Dice similarity coefficient (pages 60, 106) DVH dose volume histogram (page 35) EBRT external beam radiotherapy (page 3) EORTC European Organization for Research and Treatment of Cancer (page 32) EQD2 biologically equivalent dose (page 30) ESTRO European Society for Radiotherapy & Oncology (page 36) eV electron volt FDM finite difference method (page 71) xxv FEM finite-element-method (page 47) FFD free form deformations (page 50) FLE fiducial localization error (page 89) FWER family wise error rate (page 134) GEC Groupe Européen de Curiethérapie (page 36) GMM Gaussian mixture model (page 73) GTV gross tumour volume (page 34) GU genitourinary (page 116) Gy Gray (page 12) GYN Gynaecological HD Hausdorff distance (pages 61, 87) HDR high-dose-rate (page 17) HPV human papillomavirus (page 1) HU Hounsfield unit (page 42) ICBT intracavitary brachytherapy (page 36) ICE inverse consistency error (page 59) ICP iterative closest point (pages 52, 56, 70) ICRU International Commission on Radiation Units and Measurements (page 34) kV kilovoltage LENT-SOMA Late Effects of Normal Tissue-Subjective, Objective, Management, and Analytic (page 32) LDR low-dose-rate (page 17) Linac linear accelerator (page 13) LQ linear quadratic (pages 27, 30) mL milliliters MLC multi-leaf collimators (page 14) MLE maximum likelihood estimation (page 73) MRI magnetic resonance imaging (page 33) xxvi NE normalized error (page 86) NRE normalized registration error (page 122) OAR organ at risk (page 31) PD prescribed dose (page 35) PDR pulsed-dose-rate (page 19) PET positron emission tomography PSR point-set registration (page 118) PTV planning target volume (page 34) RBF radial basis functions (page 49) RDE residual distance error (page 105) RMS root mean square (page 87) RMSD root-mean-square distance (page 51) RMSE root mean square error (page 62) ROI region of interest (page 33) RPM robust point matching (pages 67, 71) RSD root squared distance (page 51) RTOG Radiation Therapy Oncology Group (page 32) SC shape context (page 72) SDE surface distance error (page 105) SSD sum of squared differences (page 52) SSM statistical shape model (page 67) STPSRPM Symmetric Thin Plate Spline Robust Point Matching (page 68) TPS thin-plate-spline (pages 49, 72) TPS Treatment planning software (pages 22, 128) TRE target registration error (pages 60, 68) SF surviving fraction (page 27) SOR successive over-relaxation (page 48) xxvii Acknowledgements First and foremost, I would like to extend my deep and humble gratitude to my supervisor, Dr. Ingrid Spadinger. Without your enduring support and comprehensive and profound knowledge of the subject, I could have not completed this quest. Your constant guidance and insight, elaborate and sleek reviews, pleasant and positive attitude made this challenging journey a lot more enjoyable. I could have not asked for a better supervisor. My sincere appreciation also goes to Dr. Ghassan Hamarneh for playing such a valuable and crucial role in guiding the computing aspect of this study and always providing the most useful advice and insight. I would also like to use this opportunity to thank Dr. Christina Parsons for providing the ever so helpful clinical insight and feedback. I am grateful to my academic supervisor, Dr. Alex Mackay, and Dr. Carl Michal in my supervisory committee for commitment of their time and expertise to this project. My special thanks go to Colin Brown (PhD candidate in Dr. Hamarneh’s medical image processing lab at Simon Fraser University) because of his expertise in the field and helping me to learn different aspects of the field of image computing and registration, which was a completely new field for me to explore. Also, I would like to thank Alborz Amir Khalili at UBC for helping with one of the algorithms used in the study. I would also like to thank Dr. Marc Gaudet for sharing all his research data. I also thank staff members in various fields at the BC Cancer Agency, who taught me a lot during my time there. In particular, I want to genuinely thank two people: Rustom Dubash whose expertise in the field of HDR brachytherapy and Eclipse treatment planning software came to my xxviii rescue in a lot of occasions, and Dr. Monty Martin in Radiology, who taught me how to see things in a CT that I couldn’t see before. I gratefully acknowledge that this research has been funded by the BC Cancer Agency and for that I have to thank Dr. Cheryl Duzenli for her continuous support throughout this project. Lastly, I would like to express my heartfelt gratitude to my parents and my husband for their endless love and unwavering support during my studies. I could not have fulfilled my intellectual ventures without their enduring encouragement and desire to see me thrive. I am ever so grateful for that. xxix Dedication To the very special people in my life: My loving husband, Damon and My dear parents, Pari and Mohammad 1 Chapter 1: Introduction 1.1 Cervical Cancer Cancer is one of the most common life threatening diseases known to humans, claiming many lives globally every day. It is the leading cause of death in Canada and many other countries in the world. Cancer occurs when the genetic material of the cells in the body is damaged in such a way that the cells begin to multiply abnormally and create tumours. In advanced stages of the disease, abnormal cells may spread and create secondary tumours in other organs (metastasis). There are many different types of cancer, depending on the organ or system affected. Gynecological cancers, including carcinoma of the uterus, cervix, vagina, ovaries and vulva, are the third leading cause of cancer death in women in the Unites States and the fourth in Canada after lung, breast, and colorectal cancer [1, 2]. Carcinoma of the cervix (aka cervical cancer) is the fourth most common cancer in women worldwide, and it has the fourth highest mortality rate among cancers in women [3]. Each year, nearly half a million women worldwide are diagnosed with cervical cancer [4]. While exposure to human papillomavirus (HPV), mainly through sexual interaction, is the prevailing cause of cervical cancer, smoking, immune system disorder, and birth control pills are some other important risk factors [4, 5]. Regular screening, usually in the form of a Pap-smear test, is an effective method for detecting early-stage cellular changes, which can, if necessary, be treated and thereby prevent progression to cancer. If cancer is suspected, a cervical biopsy of a suspicious lesion can be used for diagnosis and staging. [4] Treatment options for cervical cancer vary, and include surgery, megavoltage external beam 2 radiation therapy, internal radiotherapy (i.e. brachytherapy), and chemotherapy, often in combination, depending on the stage of the disease (i.e. tumour size, tumour extension into adjacent tissues, and lymph node involvement) [4, 6]. For several decades before the mid-1950's, kilovoltage x-ray therapy together with brachytherapy was the standard care for treating cervical cancer, with average 5-year survival rate of 45% [7]. The current average 5-year survival rate for cervical cancer in Canada is 74%, ranging from 15% for the most advanced stage to 93% for the earliest stage [2]. 1.2 Radiation Therapy 1.2.1 Overview and Historical Perspective Radiation therapy specifically involves exposure of the cancerous tissue to ionizing radiation, either in the form of charged (i.e. electrons, protons) or uncharged particles (neutrons), or high energy electromagnetic radiation (x-rays or gamma rays). Treatment of cancer with radiation dates back to shortly after the discoveries of x-rays by Wilhelm Roentgen in 1895, which gave rise to both diagnostic imaging with x-rays and external beam radiation therapy, and radium by Marie and Pierre Curie in 1898, which led to the development of brachytherapy techniques. The first successful treatment for cancer by placing radium near tumours was in 1903, where two patients with basal cell carcinoma of the face were treated in St. Petersburg [8, 9]. This method of treatment was originally called “Radium Therapy”, but in 1931, Forsell proposed the term brachytherapy (BT), i.e. “short distance” therapy, as opposed to “long distance” therapy or teletherapy (called external beam radiotherapy nowadays) [9]. Brachytherapy treatment quickly evolved from surface applications to also include insertion of radioactive sources into body 3 cavities (intracavitary), or implanting them directly into the tumour site (interstitial brachytherapy). Surface and intracavitary treatments always end in removal of the radioactive material after the treatment is completed (temporary brachytherapy), whereas interstitial treatments can, but don’t always, involve permanent implantation. Permanent implantation requires use of a rapidly decaying radioactive source. Early external beam radiotherapy (EBRT) was limited to kilovoltage (kV) beams generated in x-ray tubes. X-rays at these energies are not sufficiently penetrating for optimal treatment of deep-seated tumours, such as those in the pelvis, without causing excessive damage to the overlying normal tissues. Pelvic sites that are accessible through intracavitary routes, such as the cervix and uterus, were therefore targeted with radium brachytherapy early on. Cervix cancer was treated with notable success in this way. Technological advances in the form of radioisotope teletherapy units (particularly cobalt-60) and megavoltage x-ray generators (e.g. betatrons, linear accelerators) during the 1940-1960’s brought high energy photon beam teletherapy to the clinic, making external beam treatment of deep seated tumours possible. EBRT soon became the radiotherapy treatment of choice for many tumour sites. However, for cervical cancer, brachytherapy has remained an important treatment modality, although it is now often used in combination with external beam therapy. Compared to conventional EBRT, brachytherapy allows for delivery of a higher total dose to a smaller volume over a shorter amount of time [10]. Due to the rapid fall-off of radiation intensity with inverse square of distance from a radiation source, a very high dose of radiation is delivered to the cancerous tissue close to the radioactive source, while surrounding normal tissues are largely spared. For bulky tumours and/or tumours that have extended beyond the cervix, EBRT is used 4 to treat the lymphatics and shrink the tumour, and brachytherapy is used to boost the area of gross tumour. 1.2.2 Physical Basis of Radiation Therapy 1.2.2.1 Ionizing Radiation Ionizing radiation is a beam of particles or electromagnetic radiation with sufficient energy to eject electrons from atoms in a medium, leaving them ionized. The energy threshold for ionization depends upon the material with which the radiation is interacting, but is generally considered to be on the order of 10-33 eV for biological systems. For electromagnetic radiation (photons), ionizing radiations above 124 keV are classified as either x-rays or gamma rays (γ-rays), with energies above approximately 20 keV being useful for radiotherapy purposes. Since the energy ranges of x- and -rays overlap to some degree, the most widely applied distinction between these two classifications is the physical process that gives rise to the emission. X-rays are defined as being produced by electron interactions, either through energy level transitions of orbital electrons (characteristic x-rays) or by interaction of incident high energy electrons with the electric field of atoms in the absorbing material (bremsstrahlung x-rays). Gamma rays are emitted during energy level transitions within atomic nuclei that are left in an excited state after events such as radioactive decay or nuclear fission. They may also result from subatomic particle decay (e.g. annihilation events). While gamma rays typically have higher energies than characteristic x-rays, modern technology has enabled the production of extremely high energy bremsstrahlung x-rays. In medical applications, x-ray beams with peak energies up to 25 MeV are not uncommon, whereas gamma 5 rays typically fall within the 0.1-2 MeV range. With some exceptions, x-rays are the primary modality in external beam radiotherapy and radiographic imaging, while gamma rays predominate in brachytherapy and nuclear medicine imaging. Charged particles with sufficient energy produce multiple excitations and ionizations of atoms in their path as they slow down. These are the events that give rise to biological damage in exposed tissues. In contrast, uncharged particles, such as photons (electromagnetic radiation interact sparsely with the medium, but many of these interactions result in a large energy transfer to a charged particle, which then produces multiple ionizations and excitations and, consequently, biological damage. The predominant ionizing radiations used in radiation therapy are photons, which set electrons in motion, and, less commonly, direct electrons. The interactions of photons and electrons most relevant to radiation therapy are described in the following sections. 1.2.2.2 Electron Interactions High-energy electrons create multiple excitations and ionizations in the medium through which they are traveling. Excitation occurs when the electron passes near an orbital electron and transfers a portion of its energy, raising the electron to a higher energy state but not ejecting it completely. Ionization occurs when the incident electron transfers enough energy to the orbital electron to eject it from the atom, leaving it in an excited state. If the acquired energy of the expelled electron is sufficient to cause secondary ionizations, it is called a 𝛿 ray. Each ionization or excitation event is termed a collisional event and results in an incremental loss of kinetic energy to the medium. The energy that is absorbed by the medium in this process is the radiation dose. A vacancy created in an inner atomic orbital by an ionization or excitation event is promptly 6 filled by an outer-shell electron. The energy lost in this transition is converted to a characteristic photon emitted with energy E=h, where E is equal to the difference in the energy levels of the involved electron orbitals, is the frequency and h is Planck’s constant. If the transition is of sufficient energy, a characteristic x-ray results. As mentioned previously, high-speed electrons can also lose their energy through radiative collisions with atomic nuclei. In this process, known as bremsstrahlung, an electron passing close to a nucleus is deflected by the strong electric field and loses some or all of its energy, which is converted to photon energy (Figure 1.1a). These photons can have energies up to the incident kinetic energy of the high-speed electron. Photons of sufficient energy are called bremsstrahlung x-rays. The probability of radiative collisions increases substantially with both electron energy and atomic number (Z) of the absorbing medium (Figure 1.1b). Because of the Z-dependence, radiative collisions are relatively uncommon in tissues at typical diagnostic and therapeutic energies. 0.010.101.0010.00100.000.01 0.10 1.00 10.00 100.00% of interactions that are radiativeelectron energy (MeV)waterlead(b) ℎ𝜈 𝑒− 𝐸 𝐸 − ℎ𝜐 (a) 7 Figure 1.1 Illustration of the Bremsstrahlung process (a) and its prevalence (b). 1.2.2.3 Photon Interactions Ionizing photons interact with probability µ per unit distance travelled in the absorbing medium. The parameter µ, known as the linear attenuation coefficient, is dependent on photon energy and the composition of the absorbing medium, as shown in Figure 1.2. Since this attenuation of the incident photon beam is highly dependent on the mass density (𝜌) of the medium, it is often expressed in the form 𝜇 𝜌⁄ (the mass attenuation coefficient). For water-equivalent absorbers, such as tissue, photon beams become more penetrating with increasing energy (i.e. 𝜇 𝜌⁄ decreases) up to about 30 MeV. Thus, high energy photons between approximately 1-30 MeV can effectively treat deep-seated tumours in the body. If each interaction is considered to remove a photon from the incident (primary) beam, attenuation of a monoenergetic beam of photons through thickness x of an absorber is governed by an exponential relationship, 𝑁 = 𝑁0𝑒−𝜇𝑥, (1.1) where N0 is the number of photons incident on the absorber and N is the number of transmitted photons. Several types of photon interaction processes contribute to µ. These include coherent scattering, Compton (incoherent) scattering, photoelectric effect, and pair production, with the latter three resulting in significant energy transfer to a charged particle. Each of these interaction types can be assigned its own attenuation coefficient, such that the total linear or mass attenuation coefficient of an absorbing medium is the sum of the individual coefficients for each interaction, 8 i.e.: 𝜇 𝜌 = 𝜎𝑐𝑜ℎ 𝜌⁄ + 𝜏 𝜌⁄ + 𝜎 𝜌⁄ + 𝜅 𝜌⁄⁄ , (1.2) where 𝜎𝑐𝑜ℎ 𝜌⁄ , 𝜏 𝜌⁄ , 𝜎 𝜌⁄ , and 𝜅 𝜌⁄ are the mass attenuation coefficients for coherent scattering, photoelectric effect, Compton scattering, and pair production, respectively. The energy dependence of these attenuation coefficients is illustrated in Figure 1.2 for water and lead. In this figure, water can be considered representative of soft tissue, which is composed primarily of water and other low atomic number constituents. 9 Figure 1.2 The mass attenuation coefficient as a function of photon energy, for coherent scattering, incoherent (Compton) scattering, photoelectric effect, pair production, for water and lead. The total mass attenuation coefficient (sum of all coefficients) is also shown. Data obtained from [8]. Figure 1.3 shows schematics of each type of photon interaction process. In coherent (Rayleigh) scattering, photons are scattered by passing near the electrons in the atomic cloud with negligible exchange of energy (Figure 1.3a). The atomic electrons oscillate due to the electromagnetic field of the radiation, but the atom does not reach an excited state (i.e. the electrons remain in their orbitals). Coherent scattering is prevalent at low energies, extending to somewhat 10 higher energies as the atomic number increases. In the photoelectric effect (Figure 1.3), the photon collides with an orbital electron of an atom, transferring all its energy (ℎ𝜈) to the bound electron. The ejected electron, called the photo-electron, leaves the atom with kinetic energy equal to ℎ𝜈 − 𝐸𝑏, where 𝐸𝑏 is the binding energy of the shell from which the electron is ejected. Photoelectric interactions are more likely to occur in the inner shells, provided ℎ𝜈 𝐸𝑏. As described previously in the context of electron interactions, characteristic photons are emitted when the resulting electron orbital vacancies are filled. The photoelectric effect is highly dependent on the atomic number (Z) of the absorber, with the probability of interaction increasing in proportion to Z3. This important property is utilized in radiographic imaging, which employs x-ray energies in the ~10-150 keV range, where the photoelectric effect is sufficiently dominant to provide excellent bone contrast, as bone has a higher effective Z value than soft tissue. In radiation therapy, where higher energies are typically used, enhanced absorption in bone is usually a disadvantage. Compton, or incoherent, scattering is a process in which the photon transfers some of its energy (ℎ𝜈) to an orbital atomic electron (mostly from outer orbitals with negligible binding energy) and scatters at an angle (𝜃) from its original direction with reduced energy ℎ𝜈′ (Figure 1.3c). The Compton effect predominates in the therapeutic energy range (0.1-25 MeV) in soft 11 Figure 1.3 Schematic diagrams of photon interactions with atoms, illustrating (a) coherent scattering, (b) photoelectric effect, (c) incoherent (Compton) scattering, and (d) pair production. tissue and is not dependent on atomic number. If the energy of the incident photon is greater than 1.02 MeV (i.e. twice the rest mass of an electron or positron) and it passes near the nucleus of an atom, pair production can occur. In this process, the photon interacts with the strong electromagnetic field of the nucleus and vanishes, giving up all its energy and creating an electron-positron pair (Figure 1.3d). The difference between the energy of the incident photon and the total rest mass of the electron and positron (i.e. 1.02 MeV) is shared between the kinetic energies imparted to the electron and the positron. The positron loses energy through ionization and excitation interactions with absorbing medium, and 12 is eventually annihilated in a final interaction with an electron that converts the rest mass of both particles into two gamma rays moving in opposite directions. Some pair production occurs in soft tissues exposed to very high energy therapeutic x-rays, but the total contribution of pair production events is small relative to Compton events. 1.2.2.4 Radiation Intensity and Geometric Effects The radiation intensity is the amount of energy passing per unit time per unit area perpendicular to the direction of the energy flow. For a point source (either a small focal spot on an x-ray target or a radioisotope point-source) that emits ?̇? photons of energy ℎ𝜈 per second, the intensity of radiation is defined as 𝐼 =?̇?ℎ𝜈𝐴, (1.3) where A is the surface area through which the photons pass. For an isotropic point source (Figure 1.4), in the absence of attenuation, the intensity decreases with the square of the distance from the source (𝑟): 𝐼 ∝1𝑟2 (1.4) This is known as the “inverse square law”, which is also valid for planar cross-sectional areas perpendicular to the beam direction in external beam radiotherapy. Because they are extreme at short distances, inverse square effects strongly govern the dose distribution around brachytherapy sources. 13 Figure 1.4 Illustration of an isotropic point source emitting photons. 1.2.2.5 Radiation Dose As described previously, energy is deposited by electrons that are either directly incident on the medium or have been set into motion by photon interactions. Energy deposited by ionization and excitation events along the electron tracks is quantified by the absorbed dose, which can be correlated with the biologically significant effects produced by the ionizing radiation, and is defined as: Absorbed Dose = 𝑑?̅?𝑑𝑚 , (1.5) where 𝑑𝜖 ̅ is the mean energy deposited by the ionizing radiation in mass 𝑑𝑚 of the absorbing medium [11]. The SI unit for the absorbed dose is gray (Gy), defined as 1 gray = 1 J/kg. An older unit for absorbed dose is rad, where 1 rad = 10-2 Gy = 1 cGy (i.e. centigray). In the case of photons, the initial photon interaction events that transfer energy to electrons are largely responsible for the pattern of dose deposition, and are described by a quantity called kerma, i.e. kinetic energy released in the medium. Kerma is defined as K = 𝑑𝐸𝑡𝑟𝑑𝑚, (1.6) 14 where 𝑑𝐸𝑡𝑟 is the total initial kinetic energy of all the electrons released in a medium of mass 𝑑𝑚 [11]. While the SI unit for kerma is the same as the absorbed dose, the two quantities represent distinct processes that have different, but closely related, spatial distributions. The distribution of kerma is governed by the inverse square and photon attenuation relationships, as applied to the primary and scattered photons, while the distribution of dose is further influenced by the less predictable, but spatially limited, electron collisional events that follow each kerma event. 1.2.2.6 Sources of Therapeutic Photons and Electrons Therapeutic photons and electrons are generated either in electron accelerators, or through radioactive decay processes. Both methods are or have been applied in the delivery of the two broad categories of radiation therapy described previously: external beam therapy (“teletherapy”), where the radiation is generated externally, collimated to various cross-sectional field sizes and shapes, and delivered to the tumour site from outside the patient, and brachytherapy, where an uncollimated radiation source (or sources) is deposited inside or in very close proximity to the tumour. 1.2.2.6.1 Teletherapy (External Beam) Photons and Electrons Teletherapy beams are usually produced by electron accelerators, although large radioactive sources are sometimes also used to generate photon beams (such as in Cobalt-60 teletherapy units). In modern radiotherapy, the most frequently used devices are linear accelerators (i.e. Linac), which produce both x-ray beams and direct monoenergetic electron beams, and x-ray tubes, which are only used for generating relatively low energy (<500 kV) x-rays. X-ray tubes are also heavily utilized in radiographic imaging. In both devices, x-rays are produced by bombarding a high Z 15 target (typically made of tungsten) with high-speed electrons. Although characteristic x-rays are also generated, bremsstrahlung x-rays predominate the emissions, giving rise to a beam composed of a spectrum of photon energies ranging from zero to the kinetic energy of the accelerated electrons and a mean energy of approximately one third of the maximum. X-ray beams are nominally described in terms of the maximum energy, expressed as the accelerating potential through which the electrons have traveled (e.g. 15 MV). 1.2.2.6.2 External Beam Field Shaping and Intensity Modulation In external beam therapy, beams are custom collimated and are usually delivered from multiple directions aimed at the centre of the tumour volume. Customized collimation is achieved with movable jaws integrated in the treatment unit head, which define rectangular fields. On a finer scale, individually controlled multi-leaf collimators (MLC) or fixed shielding can further conform the shape of the aperture to the target volume. The pattern of dose delivery can also be tailored by spatially modulating the beam intensity or intensity rate, by modulating the beam dose rate, or through the introduction of attenuating devices into the path of the beam, or by movement of the MLC's to predefined locations during beam delivery. Figure 1.5 shows an example of a 4-field EBRT treatment plan. This plan has conformal apertures but no intensity modulation, and represents a typical EBRT plan for treatment of cervix cancer. 1.2.2.6.3 Brachytherapy Sources With the exception of relatively recently developed electronic brachytherapy sources in the form of “mini” x-ray tubes, brachytherapy sources are made from radionuclides. The original radioactive elements used for brachytherapy were naturally occurring radium and its daughter pro- 16 Figure 1.5 The dose distribution of a 4-field box treatment plan for cancer of cervix, shown as a colour wash. The bladder (yellow), rectum (green) and clinical target volume (red) contours can also be seen. duct, radon. Due to several disadvantages of radium and radon, they were eventually replaced by selected artificially produced radionuclides after the development of nuclear reactors made these available. Commonly used radioisotopes in modern brachytherapy include Iridium-192, Cesium-137, Iodine-125, Palladium-103, and Cobalt-60. Radioactive decay occurs in atomic nuclei with an unfavourable balance of protons and neutrons. Amongst the possible types of decay, modern brachytherapy sources almost all decay either by – emission or electron capture. In – emission, a neutron within the nucleus is converted to a proton through emission of a negatively charged particle, the – particle, which is actually an 17 electron: ~011110 p n (1.7) Also emitted is an anti-neutrino ( ~ ), which interacts minimally and is therefore inconsequential to the aims of radiotherapy treatment. Often, the nucleus is left in an excited energy state after the decay, and emits one or more -rays in order to reach the ground state. It is these -rays that are the most useful for brachytherapy, although in limited applications (i.e. treatment of very shallow, usually superficial lesions) the -particle itself is preferred. In those instances, a radionuclide that is a pure emitter (i.e. it emits no -rays) is used. If -rays are preferred, the – particle can be stopped by encasing the radioactive material in a metallic capsule, which is also useful for other purposes, such as containing the radioactive material. Some modern brachytherapy sources, such as Iodine-125, undergo decay by electron capture, where an inner orbital electron is captured by the nucleus, converting a proton to a neutron: nep100111 (1.8) where is a neutrino. Depending on the radionuclide, electron capture may be accompanied by -ray emission. Characteristic x-ray emission also occurs when the vacancy in the inner electron orbital is filled. While the energies of the photon emissions determine, to some degree, the penetrating ability of the emitted radiation, this is important primarily for radiation safety purposes, as the dose distribution close to brachytherapy sources is predominated by inverse square fall-off. Many radionuclides undergo several different decay processes with differing probabilities. For example, Iridium-192, which is the most commonly used radionuclide for high-dose-rate 18 brachytherapy sources, has several – decay modes as well as 3 electron capture modes. Other important characteristics of a radionuclide include the source strength and the rate of decay. Source strength affects the dose rate delivered by the brachytherapy source, and is measured in terms of activity (number of decays per unit time) or air kerma strength (SK): 𝑆𝐾 = ?̇?𝛿(𝑑)𝑑2, (1.9) where 𝐾?̇? is the air kerma rate at distance d from the source. Source strength diminishes over time as the source decays, where the rate of decay is described by 𝑆𝐾 = 𝑆𝐾0𝑒−𝜆𝑡. (1.10) is the decay constant, representing the rate of decay per unit time, and 𝑆𝐾0 is the initial strength. An alternative measure of the rate of decay is the half-life (𝑇1/2), which is the time required for a sample to decay to half of its original strength. From Eq. 1.10, the half-life 𝑇1/2 = 𝑙𝑛2/𝜆. Traditionally, brachytherapy sources had sufficiently low source strengths that they could be handled by the practitioner without acute effects, but a tumouricidal dose required protracted treatments with multiple sources over the course of hours, days, or even as permanent implants. This is referred to as Low Dose Rate (LDR) brachytherapy, and is still effectively practiced today for a number of tumour sites. However, with the advent of artificially produced radionuclides in the middle of the 20th century, it became possible to produce small but very high strength radioactive sources. These cannot be handled manually, but additional technological advances allowed them to be partnered with computerized, remotely controlled delivery mechanisms. This is called High Dose Rate (HDR) brachytherapy, and can be delivered in minutes. Figure 1.6 shows an example of an HDR 19 source, which has an iridium-192 radioactive core encapsulated in stainless steel. Figure 1.6 The GammaMed high dose rate Iridium-192 source1. 1.2.2.6.4 Brachytherapy Delivery Systems and Methods There are many ways in which brachytherapy can be delivered, depending on the size and location of the target area, source design, source strength and half-life, and duration of treatment. In the early days of brachytherapy, radioactive sources, usually in the form of radium tubes or needles, were manually positioned by an operator, left in for a period of hours or days, and then removed. Permanent interstitial implants using small seeds containing radon gas, which has a short half-life, were also performed. Due to the risk of radiation exposure to the medical staff, manual positioning of radiotherapy sources is now limited to low dose rate applications using low energy radionuclides [12]. Manual afterloading was introduced in the 1950s, and was widely used afterwards because it reduced exposure to the caregivers and patients [13]. In this method, hollow needles, catheters, or 1 Figure adapted from: Perez-Calatayud J, Ballester F, Das RK, Dewerd LA, Ibbott GS, Meigooni AS, Ouhib Z, Rivard MJ, Sloboda RS, Williamson JF. Dose calculation for photon-emitting brachytherapy sources with average energy higher than 50 keV: report of the AAPM and ESTRO. Med Phys. 2012;39(5):2904-29.doi: 10.1118/1.3703892. 0.9 mm 0.7 mm 4.52 mm STEEL CABLE ACTIVE Ir-192 CORE LASER WELD 3.5 mm 0.3 mm STAINLESS STEEL CAPSULE 20 larger applicators containing a narrow channel are positioned in the treatment site. The location of the applicator can then be confirmed radiographically if required, after which the radioactive sources are introduced, still manually but much more rapidly than in manual loading, thereby reducing operator exposure. The concept of afterloading was soon applied to the development of remotely controlled afterloading devices, where the sources are introduced into the applicators using a remotely driven, shielded mechanical system. In addition to almost eliminating radiation exposure to the operator, these devices offer greater flexibility in tailoring the treatment plan, which, with the help of radiographic imaging, can be done on the basis of the actual applicator positioning in situ. Early systems, developed in the 1960's and 1970's, used multiple low dose rate sources that could be introduced into the applicators in various configurations. One of the most widely used systems, the Selectron LDR, used a pneumatic drive system combined with a magnetic sorting mechanism to drive combinations of active Cesium-137 pellets and inactive spacer pellets into the applicators under microprocessor control. Systems such as this one were used widely for intracavitary treatment of cervical cancer over the span of 30-40 years. The development of remote control afterloading technology ultimately led to the development of high dose rate brachytherapy delivery systems in the 1980's. In these systems, a single high intensity, cable-driven “stepping” source "dwells" for a pre-set amount of time in a sequential set of defined positions within the applicators (Figure 1.7). Treatment times are short and patient recovery is often very quick, therefore it can be done on an outpatient basis for many sites [14]. In addition, because dwell times are highly selectable, treatment plans are more customizable than is practically achievable using multi-source LDR techniques. 21 Despite early concerns regarding the radiobiological consequences of such a radical change in dose rate (discussed in more detail in Section 1.2.3), HDR is very popular currently, and has in recent years largely replaced LDR as the method of choice for brachytherapy of cervical cancer. A third option is pulsed-dose rate (PDR) brachytherapy, which uses short pulses of HDR radiation, usually about once per hour, in order to simulate the overall dose rate of LDR brachytherapy. Figure 1.8 shows an example of an intracavitary tandem and ovoid applicator system for the treatment of cervical cancer. This basic configuration, where the tandem passes through the cervix into the uterine canal, and the ovoids are nestled in the vaginal vault next to the cervix, was developed in the days of manual radium brachytherapy and has been gradually adapted for use with modern afterloading devices. A modification of this basic design, the tandem and ring (also shown in Figure 1.8), is used only with HDR units, which can accommodate the tight radius of curvature of the ring. Specialized applicators also exist for a number of other tumour sites. 1.2.2.7 Patient Dose Calculations Radiation therapy is delivered according to a prescribed dose and preferred dose distribution, the specifications for which can vary considerably according to the tumour type, stage, site, and usage of additional therapies (i.e. surgery and chemotherapy), as well as the properties and limitations of the radiation delivery method and modality. Additional factors such as patient size, body composition, tumour size and location, and proximity of sensitive normal tissues can also affect the design of the treatment plan. The parameters specified in an external beam therapy treatment plan include beam energy, dimensions, direction of incidence, and duration of exposure (amongst others). The parameters in brachytherapy treatment planning consist of the type, number and distribution of sources (or source 22 Figure 1.7 Diagram of the source drive mechanism of the Nucletron microSelectron HDR remote-control afterloading unit. The radiation source is permanently fixed to a cable and is driven to a number of different “dwell” positions within the applicator using a precision stepping motor. A second cable with a dummy source is used to check the integrity of the applicator tube before attempting to drive the active source through the tube. Each guide tube in the intracavitary brachytherapy applicator connects to one of the 18 channels. delivery positions in the case of HDR), source strength, and duration of exposure. Parameters are decided upon using a combination of prior experience and patient-specific dose calculations based on the chosen parameter values. These calculations provide an accurate estimate of the expected dose distribution. In both external beam therapy and brachytherapy, dose calculations are based on measurement of the radiation dose distributions produced by the radiation source, combined with consideration of the basic physical processes that underlie the pattern of dose deposition. In the case of photons, for example, these include attenuation, scatter, and inverse square dose fall-off. Of these, attenuation and inverse square effects are relatively easy to model mathematically. Photon scatter and electron interactions are more difficult to model. 23 Figure 1.8 Top: Examples of tandem-ring (Left) and tandem-ovoid (Right) intracavitary applicator used for cervical cancer HDR brachytherapy1. The radiation sources are guided through the tandem to irradiate the paracervical region and separately through the ring/ovoids to irradiate the cervix. The rectal paddle can be used optionally to push the rectum away from the radiation sources. Bottom: The tandem-ring applicator in-situ. The tandem lies inside the uterus and the ring sits at the cervix. 1.2.2.7.1 External Beam Dose Calculations The dosimetric data used for EBRT treatment planning can be acquired by the end user by taking measurements in a flat, homogeneous tissue-equivalent (usually water) phantom. For each modality and beam energy, data are collected for a series of standard field sizes. The measurements normally include machine output calibrations and relative dose measures along and away from the central beam axis (i.e. beam profiles). Application of these data to patient dose calculations 1 Open source: http://www.aboutcancer.com/intracavitary_radiation_treatments.htm 24 requires correction for varying field apertures, tissue composition and body surface contours, which is carried out by algorithms integrated into specialized treatment planning software (TPS) available from several vendors. 1.2.2.7.2 Brachytherapy Dose Calculations Measuring dose distributions around brachytherapy sources is a very challenging task because it needs to be done close to the sources, where dose gradients are very steep due to inverse square effects. Therefore, brachytherapy sources are characterized generically by one or two investigators through a carefully designed process that involves both measurement and computation of expected dose distributions in water using numerical (Monte Carlo) modeling of the photon interactions. These dose distributions are then parameterized using a defined set of factors as described by the American Association of Physicists in Medicine Task Group protocol AAPM TG-43 [15, 16]. Since the first publication in 1995 [15], the TG-43 protocol has been the accepted dose calculation formalism for brachytherapy sources in North America. The TG-43 framework accounts for anisotropy in the dose distributions caused by the source material and design, but does not correct for heterogeneity in tissue composition and density. However, this is usually not critical for brachytherapy treatments because the source(s) are in or around water-equivalent soft tissues. There are brachytherapy algorithms available that account for heterogeneities [17, 18], but they are not yet in common use. The cylindrical geometry of a brachytherapy source as used in the TG-43 protocol in shown in Figure 1.9. The dose is normalized with respect to a reference point (𝑟0, 𝜃0), where 𝑟0 = 1𝑐𝑚 and 𝜃0 = 90°. For a line source as shown in Figure 1.9, the dose rate ?̇?(𝑟, 𝜃) at the point 𝑄(𝑟, 𝜃) 25 Figure 1.9 The geometry of a brachytherapy cylindrical source, as described in the TG-43 formalism. in water is given by ?̇?(𝑟, 𝜃) = 𝑆𝐾Λ𝐺𝐿(𝑟,𝜃)𝐺𝐿(𝑟0,𝜃0)𝑔𝐿(𝑟)𝐹(𝑟, 𝜃). (1.11) where Λ is the dose rate constant (dose rate in water per unit air kerma strength at the reference distance 𝑟0), which is specific to the source model. 𝐺𝐿(𝑟, 𝜃) is a geometry factor for a line source, which accounts for the inverse square law effects on the dose distribution. It uses an approximate model of the spatial distribution of the radioactivity within the source and is calculated as 𝐺𝐿(𝑟, 𝜃) = {𝛽𝐿 𝑟 𝑠𝑖𝑛 𝜃 𝑖𝑓 𝜃 ≠ 00(𝑟2 − 𝐿2 4⁄ )−1 𝑖𝑓 𝜃 = 00 (1.12) Finally, 𝑔𝐿(𝑟) and 𝐹(𝑟, 𝜃) in Eq. 1.11, are the radial dose and the anisotropy functions and are both specific to the source model. The radial dose function, 𝑔𝐿(𝑟), accounts for the variation of dose with distance along the transverse axis (i.e. 𝜃 = 900) of the linear source due to the absorption and scatter in the medium. The anisotropy of the dose distribution, excluding the effects accounted for by 𝐺𝐿(𝑟, 𝜃), is accounted for by the anisotropy function 𝐹, where 𝐹(𝑟0, 𝜃0) = 1. active length, L long axis of source encapsulation r Q(r, ) transverse axis of source 26 The total dose delivered by a source is obtained by integrating the dose rate from Eq. 1.11, over the period of time during which the source is in place. If the source half-life is short relative to the treatment duration, the decay of the source strength over time must also be taken into account. Moreover, a typical brachytherapy treatment involves multiple sources or source dwell positions. The total dose delivered at each point in the region of interest must therefore be calculated by superposition of the dose contributions from all the source positions. 1.2.3 The Radiobiology of Cancer Treatment Ionizing radiation absorbed by living tissue can create a range of biological effects, including (but not limited to) gene mutation, carcinogenesis, loss of function, and proliferative cell death. The primary cause of these biological effects is damage to the cellular DNA [19]. Damage can result from either direct interaction of charged particles with the atoms in the DNA molecules, or by indirect damage mechanisms, where the charged particles interact with the water content of the tissues to produce potent free radicals that can damage DNA. Cells are capable of repairing some types of lesions (i.e. ruptures) in their DNA. These are referred to as sub-lethal lesions, as there is potential for them to become more serious if, for example, they are not repaired before a subsequent critical event (e.g. cell division or another lesion induced nearby). Other types of lesions are more severe, and if incorrectly repaired or left unrepaired, can lead to proliferative cell death. These are said to be lethal lesions. Of the biological effects caused by radiation exposure, proliferative cell death underlies the success of radiotherapy as a cancer treatment, and has therefore been extensively studied and modeled in order to gain a better understanding of the underlying mechanisms and thereby design more effective treatment regimens. 27 1.2.3.1 Mechanism of Proliferative Cell Death The primary mechanism of proliferative cell death is widely considered to be breaks in both strands of the double-strand DNA molecule in close proximity to each other [19], which can lead to physical dissociation of the two ends of the molecule. These are termed double-strand breaks, and, while they are occasionally repaired successfully by the cell, they often lead to irreparable damage that prevents cell division. Single strand breaks, on the other hand, can in most circumstances be repaired faithfully using the intact complementary section of the DNA molecule as a template for the repair enzymes. 1.2.3.2 Radiation Damage to Normal Tissue Cell killing in the normal tissues exposed during radiotherapy gives rise to different consequences depending on the affected organs and body systems. Some radiation side effects start shortly after onset of radiation treatment and are called acute or early effects. Acute effects are usually present when rapidly proliferating tissues have been exposed, and are often temporary if the damage is not too severe. Examples include skin reactions such as erythema (skin reddening), epilation (loss of hair), and gastrointestinal symptoms (nausea, diarrhea, and rectal bleeding). Some radiation effects do not develop until months or even years after irradiation, and usually last longer. These are called late or chronic effects and are attributed to accumulated dead or mutated cells in slowly proliferating tissues. While therapeutic radiation doses are limited so that late effects do not typically become life threatening, they do affect the quality of life and can cause serious health issues. Radiation can also cause genetic mutations that increase the risk of developing cancer in the exposed tissues, which may not occur until several years after radiation. Large clinical studies have shown that radiotherapy is associated with a statistically significant, 28 but very small chance of developing secondary cancers [19]. 1.2.3.3 The Linear Quadratic Cell Survival Model The double strand break as the lethal lesion is the basis of the most widely used model of proliferative cell death: the linear-quadratic (LQ), or - model. The LQ model provides a good fit to survival data obtained from mammalian cells irradiated in vitro by gamma or x-rays. It has also proven to be useful as a basis for the analysis and interpretation of far more complex in vivo dose response data. Assuming there are on average ?̅? lethal lesions per irradiated cell, the probability of a cell having zero lethal lesions can be described by the Poisson distribution. Therefore, the surviving fraction (SF) of the cells in an irradiated population decreases exponentially as the lethal lesions increase, i.e. 𝑆𝐹 = 𝑒−?̅? (1.13) The quantity ?̅? depends directly on the radiation dose. Two mechanisms of producing double strand DNA breaks (DSB) are considered: one directly proportional to dose (𝛼) and the other to the square of dose (𝛽). The 𝛼 coefficient represents the probability of DSB caused by a single ionization track, therefore it depends directly on the first power of dose [20, 21]. The 𝛽 coefficient represents the probability of two lesions produced by two separate ionization tracks, but in close enough proximity to cause a DSB. This parameter, therefore, depends on (dose)2 [22]. Thus, 𝑆𝐹 = 𝑒−(𝛼𝐷+𝛽𝐷2) [23]. (1.14) The ratio between 𝛼 and 𝛽 (known as the 𝛼/𝛽 ratio) characterizes the shape of the survival 29 Figure 1.10 Representative (not quantitative) survival curves for tumour and early-responding tissues (larger 𝜶/𝜷 ratio) and that of late-responding tissues (smaller 𝜶/𝜷 ratio), demonstrating higher therapeutic ratio at smaller doses compared to larger doses, in a single-hit delivery. curve (Figure 1.10), and is a reflection of the intrinsic radio-sensitivity of the affected cells. Although other factors also play a role, the radio-sensitivity of different cell types is strongly influenced by their rate of proliferation and degree of differentiation. Cells that proliferate faster and are less differentiated, such as bone marrow cells and tumour cells, are more radiosensitive. 1.2.3.4 Fractionation and Therapeutic Ratio The ultimate goal of radiation treatment is delivering enough radiation to the tumour in order to eradicate it without irradiating the normal tissue too much, causing serious complications in the affected organs. Therefore, there is always a trade-off between tumour control and normal tissue complication rate, which is recognized in the concept of therapeutic ratio. The therapeutic ratio or index is the ratio of the tumour control probability to a specified type and probability of normal tissue damage. 0.00010.0010.010.11Surviving FractionDoseLate Responding TissuesTumor & Early Responding Tissues30 A standard method that has been used to improve the therapeutic ratio in radiation therapy is fractionation, which exploits the difference in dose response between tumour and normal cells. Referring back to Figure 1.10, the initial shoulder of the survival curve is wider and less curved for more rapidly dividing tissues (larger 𝛼/𝛽 ratio) and is more curved for slowly proliferating late responding normal tissues (smaller 𝛼/𝛽 ratio). If it is assumed that sub-lethal damage can be fully repaired given a sufficient amount of time, dividing the radiation dose delivery up into fractions separated by the recovery period leads to a significant separation of the survival curves of the tumour and late-responding normal tissue, as shown in Figure 1.11, given that each fractional dose is less than the value at which the single fraction survival curves (such as those in Figure 1.10) cross. This difference in the dose response of the tumour and most normal tissues can be exploited to improve the therapeutic ratio of the radiation treatment. It can be seen in that an optimal therapeutic ratio is achieved when using a higher number of treatment fractions (Figure 1.11b), compared to fewer fractions with higher dose per fraction (Figure 1.11a). It can also be observed that the fractionated survival curves in Figure 1.11b become almost linear, suggesting that fractionation beyond the 2 Gy per fraction used in the figure will not entail a significant benefit. Therefore, fractionation of the prescribed dose into 2-Gy fractions is commonly used when radiotherapy involves significant dose to large volumes of normal tissue, such as in external beam radiotherapy. For practical reasons (i.e. due to often invasive nature of the procedure), larger doses per fraction are generally used in HDR brachytherapy, where normal tissue doses are lower. In contrast, continuous low dose rate brachytherapy is more like a series of very small fractions. 31 Figure 1.11 Surviving fraction curves for tumour (solid blue) and normal tissue (solid black) cell lines, splitting a total dose of 12 Gy into (a) two dose fractions of 6 Gy each and (b) six fractions of 2 Gy each. The solid lines depict the unfractionated linear-quadratic (LQ) survival curves for the corresponding cells. The surviving fractions for each dose fraction are shown as dots connected by the dotted curves, which represent the initial portion of the LQ curve applied to each fraction, for each cell line. 1.2.3.5 Modeling Biological Effect Because of the various ways in which radiation therapy can be delivered (e.g. spanning brachytherapy versus EBRT, total dose delivered, dose rate, and fractionation), having a means to compare the potential biological effects of different regimens is useful. Amongst other approaches, the LQ model has been expanded to provide a measure of biological effect. From Eq. 1.14, the biological effect (E) for a single acute dose can be defined as 𝐸 = − ln(𝑆𝐹) = 𝛼𝐷 + 𝛽𝐷2. (1.15) When n well-separated fractions with dose d per fraction are used, the biological effect is given by 𝐸 = 𝑛(𝛼𝑑 + 𝛽𝑑2). (1.16) This quantity is used to define two parameters that help quantify and compare the overall effect of fractionated treatment regimens. These are biologically effective dose and biologically equivalent 32 dose, which are only slightly different. Biologically effective dose (BED) is defined as 𝐵𝐸𝐷 = 𝐸𝛼 = (𝑛𝑑)×(1 + 𝑑𝛼/𝛽), (1.17) which is essentially the total dose (nd=D) times the “relative effectiveness”. An important advantage of BED is that it is cumulative when different radiation treatment regimens are used (such as external beam radiotherapy and fractionated HDR brachytherapy), if the proliferation during treatment is ignored [24]. Since 2-Gy-per-fraction is a standard fractionation regimen, it is advantageous to express the biological effect for other fractionation regimens in terms of d = 2 Gy per fraction. This quantity is called the biologically equivalent dose or EQD2, 𝐸𝑄𝐷2 = 𝐵𝐸𝐷(1+2𝛼 𝛽⁄) . (1.18) The 𝛼/𝛽 values used routinely in practice for tumour and normal cells are 10 and 3, respectively [23, 25]. These concepts can also be used to model the effect of protracted low dose rate brachytherapy, but, because the irradiation is continuous over long periods of time, this usually requires the inclusion of additional parameters to model repair rates and sometimes also repopulation. 1.2.4 Quantifying the Side Effects of Radiation Treatment The type and site of late effects after radiotherapy depends mostly on the healthy organ/s that must unavoidably be included in the irradiation volume, often referred to as the organs at risk (OAR). For example, if the radiation site is in the pelvic region, different pelvic organs such as bladder, bowel, or rectum may be affected. This can result in urinary disorders in bladder, such as 33 incontinence, blood in the urine, abnormal frequency of urination, etc., or gastrointestinal complications. In an effort to prevent or reduce these OAR complications in radiotherapy, organ-specific dose thresholds have been advised by different health authorities and research groups [26, 27]. However, many of these recommendations have relied on inconsistently reported historical data, and are not always applicable to modern techniques. Thus, in order to describe and assess the late effects in different anatomical sites more systematically, the European Organization for Research and Treatment of Cancer (EORTC) and the Radiation Therapy Oncology Group (RTOG) in North America, built a universal Subjective, Objective, Management, and Analytic (SOMA) classification system for scoring Late Effects of Normal Tissue (LENT) in radiation treatment [28, 29, 30]. This system, known as LENT-SOMA, was introduced in 1995 [28] and includes guidelines for subjective (by patient) and objective (by physician) assessment of symptoms, methods for managing the symptoms, and analytic laboratory and imaging tools for further inspection of the injuries. Various anatomic sites were included in the RTOG/EORTC report, such as brain, lung, liver, rectum, bladder, eye, etc. An example for the bladder/urethra SOMA is shown in Table 1.1. For each symptom of each organ, there are four possible scores (Grade 1-4) based on the severity of the complication, as shown in Table 1.2. Another important and widely used reference for scoring any unfavorable and unintended symptom associated with the use of a medical treatment has been published by the US Department of Health and Human Services, called Common Terminology Criteria for Adverse Events (CTCAE) [31, 32]. The Christie Hospital in UK has prepared comprehensive Subjective and 34 Objective questionnaires and forms based on both CTCAE and LENT-SOMA guidelines for head and neck, cervix, prostate, breast and bladder cancers. These questionnaires are available through their website [33], and were the ones used for cervix cancer patients treated with radiotherapy at the BC Cancer Agency (Appendix A). Table 1.1 Subjective, objective, Management, and Analytic (SOMA) criteria for bladder and urethra late effects, based on [28]. Subjective Objective Management Analytic Dysuria (painful urination) Frequency Hematuria (blood in urine) Incontinence Decreased stream Hemoglobin levels Endoscopy Maximum volume Residual volume Narcotic Alkalization / Cystectomy Iron therapy / Transfusion / Surgical intervention Pads, Catheter Cystography Volumetric analysis Contrast radiography Ultrasound Electromyography Table 1.2 Generic LENT-SOMA scoring system, based on [28]. Grade 1 Grade 2 Grade 3 Grade 4 Subjective Occasional & minimal Intermittent & tolerable Persistent & intense Refractory & excruciating Objective Microscopic (barely detectable) Intermittent macroscopic (easily detectable) Persistent macroscopic Refractory Management Occasional nonnarcotic Regular nonnarcotic Regular narcotic Surgical intervention 1.2.5 Commonly Used Terms and Parameters in Radiotherapy Treatment Planning There are a few important terms and parameters that are commonly used in the preparation of radiotherapy treatment plans. Because modern radiation therapy treatment plans are almost always based on a three-dimensional image (usually Computed Tomography (CT), but sometimes Magnetic Resonance Imaging (MRI)) of the treatment site, various regions of interest (ROI) can 35 be identified by user-defined overlays (“contours”) on the axial image slices. Specialized treatment planning software, which also overlays calculated dose distributions, is used for this purpose. To encourage consistency within the radiation therapy community, the International Commission on Radiation Units and Measurements (ICRU) has published two reports (ICRU Reports 50 and 62 [34, 35]) that rigorously define the volumes to be used. These include: Gross Tumour Volume (GTV): The gross palpable, visible, and demonstrable extent and location of the malignant growth. Clinical Target Volume (CTV): A tissue volume that contains GTV, plus a margin to account for the microscopic extent of the disease. This is the main anatomical-clinical target of the treatment and is usually the volume to which the radiation dose is prescribed. There might be more than one CTV defined for the same patient, treated separately (normally one is a “boost” to the other). Planning Target Volume (PTV): A volume that contains CTV with a margin that accounts for the effects of organ motion, patient setup, and other uncertainties in EBRT (usually not relevant for brachytherapy). The organ motion could be due to bowel movement, respiration, variation in the degree of filling of bladder and rectum, etc. Other sources of uncertainties include equipment and dosimetric uncertainties, human factors, etc. Additional ROIs including the OAR structures in the vicinity of the treatment site. For example, in HDR brachytherapy treatment for cervical cancer, typical OARs are bladder, bowel, sigmoid, and rectum. The dose distribution within the treatment site is evaluated using dose calculations to a matrix of points, usually at 1-5 mm spacing, that spans the irradiated volume. Typically, isodose overlays, 36 which are the points in 2D or 3D with the same dose, are computed from this matrix and assessed qualitatively (Figure 1.12). In addition, dose volume histograms (DVH) can be generated for each ROI. Usually, these are presented in cumulative (i.e integral) form (Figure 1.13), and certain dose-volume parameters are computed from them. Volumetric parameters, V%PD and VD, are the volumes of the corresponding structure receiving a minimum dose, expressed in terms of a percentage of the prescribed dose (PD) or an absolute dose threshold, respectively. The volumetric DVH parameters may be expressed in absolute or percentage volumes. The dosimetric parameters, D%V and DV, are the minimum dose to the hottest %V or absolute V of the structure volume, respectively. For example, for the spinal canal in Figure 1.13, 𝑉9𝐺𝑦 ≅ 37% ≅ 30 𝑐𝑚3 (82 𝑐𝑚3×0.37) and 𝐷50% = 2.4%𝑃𝐷 ≅1.1 𝐺𝑦. While providing quantitative data, DVH parameters have limitations, primarily because they provide no information about the spatial distribution of dose. For example, they cannot be used to determine if there are several small high-dose areas or a single large area. Also, in the case of fractionated radiotherapy, it is not possible to determine from DVH parameters whether the same or different sub-volumes of the OARs are receiving high doses across different treatment fractions. Another commonly used parameter in the dose assessment of image-based radiotherapy is integral dose, which is the total energy absorbed in the entire or sub-volumes of the treatment area. Integral dose is defined as 𝐼𝑛𝑡𝑒𝑔𝑟𝑎𝑙 𝐷𝑜𝑠𝑒 = ∑ 𝐷𝑖𝑚𝑖 (𝑓𝑜𝑟 𝑛 𝑡𝑖𝑠𝑠𝑢𝑒 𝑣𝑜𝑥𝑒𝑙𝑠)𝑛𝑖=1 , (1.19) where 𝐷𝑖 and 𝑚𝑖 are, respectively, the dose and mass of each tissue voxel in the volume of interest. Integral dose can be used for target volumes or OAR structures. 37 Figure 1.12 Isodose distribution within cervical HDR brachytherapy treatment volume (CTV contour shown in orange). 1.3 Image-guided Intracavitary Brachytherapy for Locally Advanced Carcinoma of the Cervix Current evidence continues to support the use of intracavitary brachytherapy (ICBT) in combination with EBRT in the management of locally advanced cervical cancer [36, 37, 38]. Moreover, HDR-ICBT has demonstrated several advantages over traditional LDR brachytherapy [39, 40]. Therefore, in 2000, the American Brachytherapy Society recommended concomitant multi-fraction HDR-ICBT and EBRT for definitive radiation treatment of cervical cancer [41]. Around the same time, the Groupe Européen de Curiethérapie (GEC) and the European Society for Radiotherapy & Oncology (ESTRO) promoted the use of 3D imaging based treatment planning for ICBT for cervical cancer treatment. They formed a working group (Gynaecological (GYN) GEC-ESTRO) of physicists and physicians that produced two comprehensive articles describing basic treatment planning concepts, terminology, and protocols for brachytherapy of cervical carcinoma [42, 43]. A L 38 Figure 1.13 Cumulative DVH plot for EBRT treatment of cervical cancer, showing the dose distributions for CTV, bladder, rectum, right and left kidneys, and spinal canal ROIs. The volumes are expressed in terms of percentages of the total structure volumes and the total structure volumes are shown in the legend. The dose values are shown in absolute dose values in the top x-axis and in percentages of the prescribed dose (45 Gy) in the bottom x-axis. HDR-ICBT for cervical cancer consists of four stages, Applicator Insertion, Imaging, Treatment Planning, and Delivery. The ICBT applicator is either a tandem-ring or tandem-ovoid, as shown previously in Figure 1.8. At the beginning of each brachytherapy treatment fraction, the applicator is inserted in the operating room. With the applicator in-situ, the patient is moved to the imaging suite, where a planning image set (either CT or MRI) is acquired. Treatment planning is then performed on the acquired images using the treatment planning software (Figure 1.14). The applicator is identified on the image and the potential source dwell positions along the applicator tubes are determined. The GTV, CTV and different OAR structures (bowel, rectum, sigmoid, and bladder) are also contoured (Figure 1.12 and Figure 1.14). The plan specifying the dwell positions and times required to produce the desired dose distribution is then Structure Volume CTV 522 cm3 Bladder 295 cm3 Rectum 94 cm3 Rt. Kidney 170 cm3 Lt. Kidney 178 cm3 Spinal Canal 82 cm3 2.4 37 50 39 Figure 1.14 Sample HDR planning CT with labelled OAR contoured structures, showing bladder, sigmoid and rectum. The Isodose curves are also displayed. created based on the prescription and OAR dose constraint criteria. Once the plan has been approved by the physician, the patient is moved to the HDR afterloader suite. Guide tubes attached to the applicator set are connected to the appropriate channels of the afterloader (as previously shown in Figure 1.7), and treatment commences. Historically, prior to the use of 3D imaging to guide treatment planning, gynecological brachytherapy planning was based on point-doses identified on planar radiographs, as proposed in the ICRU Report 38 [44], as well as older dosimetry systems, such as the Manchester system [45]. Two points (point A and B) were defined on each treatment plan based on the positioning of the applicator (Figure 1.15a), where point A was the prescription point and B represented dose to the pelvic wall. Two reference points intended to represent maximum dose to the principal OARs, i.e. 40 Figure 1.15 Left: The ICRU point A and point B locations, defined based on the applicator alignment. Right: The bladder and rectum reference points. bladder and rectum, were also identified (Figure 1.15b). This set of points is still commonly used as a link to the previous practice, and is sometimes used as a planning guideline. However, the correlation of these reference point doses with late effects, especially for bladder, is inferior compared to that of dose-volume parameters [46, 43]. This is mainly because point doses do not represent the dose distribution in the corresponding regions, especially with the high dose gradients present in brachytherapy. Moreover, dose-volume parameters allow easier summation when combining EBRT with BT. Therefore, using DVH-based dosimetric parameters for OARs, such as D1cc and D2cc, has been recommended by GYN GEC-ESTRO. In addition, it is recommended that the dose be prescribed to the CTV volume rather than to point A. Based on the above guidelines and recommendations, radiotherapy treatment of locally-advanced carcinoma of cervix at the British Columbia Cancer Agency (BCCA) currently consists of concomitant fractionated EBRT and multi-fraction image-guided HDR-ICBT. Typically, there 41 are 25 EBRT fractions with 180 cGy (total of 45 Gy) prescribed to the ICRU Report 62 reference point [35]. The EBRT planning for all fractions is performed on a single CT scan. The ICBT treatment typically consists of 5 fractions with 600 cGy peripheral dose to CTV. Occasionally, the ICBT prescription dose is modified (usually to 400, or 500, or 550 cGy) based on the suspected individual radio-sensitivity or anatomical constraints of the patient. For each ICBT treatment a planning CT is acquired with the applicator in-situ. The OAR structures, including bowel, rectum, sigmoid and bladder (Figure 1.14) are then manually contoured using tools available in the Eclipse™ planning system. Treatment planning is guided by a set of DVH dose-volume constraints, shown in Table 1.3. For each treatment fraction, the bladder is filled with water-contrast mixture using a Foley-catheter. The Foley balloon, visible inside the bladder in Figure 1.14, is used to determine the location of the bladder neck. The amount of bladder filling is controllable using the Foley-catheter system. Although consistent bladder filling has been recommended by GYN GEC-ESTRO for valid and reliable data collection [43], it has been hypothesized by Gynecological radiation oncologists at BCCA that varying the bladder volume can potentially change the location of the hotspots on the bladder wall from fraction to fraction, and thereby possibly reduce the observed urinary toxicity. Therefore, the practice of varying the degree of bladder filling at consecutive HDR fractions has been applied over the past few years. However, there have not been enough investigations to test this theory and filling this void was one of the motivations for designing this study. In order to record and monitor late morbidity due to radiation treatment, LENT-SOMA questionnaires similar to those designed by The Christie Hospital [33] are used at BCCA 42 (Appendix A). After completion of radiation treatment, the cervical cancer patients are asked to fill these questionnaires, which includes sections addressing urinary, gastrointestinal and reproductive symptoms. The bladder section encompasses all possible urinary complications: hematuria (blood in urine), nocturia (waking up during the night for voiding), dysuria (difficult/painful urination), incontinence, urgency, frequency (of urination), and change in stream. In spite of complying with the recommended dose-volume constraints for bladder (Table 1.3), undesired high-grade bladder toxicity (Grade 3-4) is observed in some patients, which cannot be easily explained by the ICRU and DVH-based parameters. This was another motivation for the study in this dissertation. Table 1.3 List of the treatment planning parameters and constraints, used for HDR ICBT for cervical cancer at the BC Cancer Agency. Plan parameter Constraint per fraction EQD2 Total ICRU Point A dose - - ICRU Rectum < 410 cGy - ICRU Bladder < 460 cGy - Bladder D2cc < 550 cGy < 90 Gy Rectum D2cc < 430 cGy < 75 Gy Sigmoid D2cc < 430 cGy < 75 cGy Bowel D2cc < 200 cGy - CTV V100% > 90% - CTV D90% > 625 cGy > 87 Gy 43 Chapter 2: Medical Image Registration and Shape Correspondence 2.1 Image Registration Image registration is the process of aligning two or more 2D or 3D images in a common coordinate system. One of the most important applications of image registration is in matching medical images, including those used in radiotherapy treatment planning and delivery. In medical image registration, images with different modalities (CT, MRI, PET, etc.) or taken at different points of time during the course of treatment are matched. These processes are called intra-patient image registration, which in radiation therapy can be used for treatment planning, dosimetric analysis, etc. Inter-patient registrations are also sometimes useful, such as in atlas-based applications where images from a population of patients are used to generate a model (atlas). In radiation therapy, inter-patient registration can be used for facilitating OAR contouring in treatment planning, and has been found especially helpful for head and neck radiotherapy. Most of the medical images used in radiotherapy applications are 3D volumes with grayscale intensity values, such as CT and MRI. In CT images, the intensity of the image voxels is described in terms of the Hounsfield unit (HU), which is a linear scale, defined as 𝐻𝑈 = 1000×𝜇−𝜇𝑤𝑎𝑡𝑒𝑟𝜇𝑤𝑎𝑡𝑒𝑟−𝜇𝑎𝑖𝑟, (2.1) where 𝜇, 𝜇𝑤𝑎𝑡𝑒𝑟, and 𝜇𝑎𝑖𝑟are the linear attenuation coefficients of the image voxel of interest, water and air, respectively. PET images are also important, but these are typically acquired in a PET-CT scanner, where the CT images are used for registration. The 3D image is a discrete map of intensity values to pixels, also known as voxels, in the image domain Ω: 𝐼 ∶ 𝒑(𝑥, 𝑦, 𝑧) ∈ Ω ⟼ 𝐼(𝒑). (2.2) 44 The image domain (or field of view) is usually a subset of the real patient space. The discrete sampling of the image domain can vary in the 𝑥, 𝑦, 𝑧 directions, i.e. anisotropic voxel resolution. For example, CT images have a higher resolution in the transverse (x-y) planes than in the longitudinal (z) direction, where the slice thickness defines the resolution. Image registration involves using a transformation model that describes the geometric mapping from the coordinate system of a moving (primary or source) image to a fixed (secondary or target) image. If we consider the moving image, 𝐼𝑚(𝒑), and fixed image, 𝐼𝑓(𝒒), with the domain Ω𝑚 and Ω𝑓, respectively, the transformation model (𝑻) is defined so that 𝑻: 𝒑(𝑥, 𝑦, 𝑧) ⟼ 𝒒(𝑥′, 𝑦′, 𝑧′) ⟺ 𝑻(𝒑) = 𝒒, (2.3) where (𝒑, 𝒒) ∈ Ω𝑚,𝑓 = Ω𝑚 ∩ Ω𝑓. (2.4) Since images have discrete domains that might even have different resolutions, the transformation 𝑻 might map a point in 𝐼𝑚 to a point that is not defined in the normal grid of 𝐼𝑓. Therefore, interpolation is used to estimate the values of the defined grid points in 𝐼𝑓. The most common interpolation methods are nearest neighbour, linear, bilinear, trilinear, and cubic interpolation. The image registration algorithm iteratively finds the optimal geometric transformation that maximizes the correspondence across the images. The degree of correspondence between the two images is measured in terms of a similarity metric. A cost function is defined based on this similarity measure and is optimized by varying the parameters of the transformation model in an iterative fashion. Any registration algorithm, therefore, involves three main components, transformation model, similarity metric, and optimization method. Also, validation of the registration outcomes must be performed to assess the accuracy and robustness of the results. 45 2.1.1 Transformation Models The transformation model can be rigid, affine, or deformable, depending on the type of transformations involved (Figure 2.1). In rigid and affine registrations, parallel lines are preserved, and a single transformation model can be defined by a single matrix that is applied to all the points in the moving image. Deformable registration, on the other hand, does not preserve parallel lines, and is therefore much more complex, requiring an anisotropic non-linear voxel-wise transformation. 2.1.1.1 Rigid-body Transformation Rigid-body registration includes translation and rotation transformations (Figure 2.1). It can be written as 𝑻𝑟𝑖𝑔𝑖𝑑(𝒑) = 𝑹. 𝒑 + 𝒕, (2.5) where the rotation matrix in 3D is built from the rotation angle (𝛼, 𝛽, 𝛾) and the translation is defined as 𝒕 = (𝑡𝑥, 𝑡𝑦, 𝑡𝑧). Therefore, rigid registration involves six degrees of freedom in 3D. The transformation can be represented as a 4×4 matrix: 𝑻𝑟𝑖𝑔𝑖𝑑(𝒑) = (𝑐𝑜𝑠𝛽 𝑐𝑜𝑠𝛾−𝑐𝑜𝑠𝛽 𝑠𝑖𝑛𝛾𝑠𝑖𝑛𝛽0𝑐𝑜𝑠𝛼 𝑠𝑖𝑛𝛾 + 𝑠𝑖𝑛𝛼 𝑠𝑖𝑛𝛽 𝑐𝑜𝑠𝛾𝑐𝑜𝑠𝛼 𝑐𝑜𝑠𝛾 − 𝑠𝑖𝑛𝛼 𝑠𝑖𝑛𝛽 𝑠𝑖𝑛𝛾−𝑠𝑖𝑛𝛼 𝑐𝑜𝑠𝛽0𝑠𝑖𝑛𝛼 𝑠𝑖𝑛𝛾 − 𝑐𝑜𝑠𝛼 𝑠𝑖𝑛𝛽 𝑐𝑜𝑠𝛾𝑠𝑖𝑛𝛼 𝑐𝑜𝑠𝛾 + 𝑐𝑜𝑠𝛼 𝑠𝑖𝑛𝛽 𝑠𝑖𝑛𝛾𝑐𝑜𝑠𝛼 𝑐𝑜𝑠𝛽0𝑡𝑥𝑡𝑦𝑡𝑧1) (𝑥𝑦𝑧1) (2.6) Rigid image registration is usually used for anatomical sites of the same patient with minor or no variations in patient position or internal anatomy, such as intra-patient registrations of intra- or inter-modality brain images. For images from different patients or other more deformable sites of the body, rigid registration fails to correctly align the images. However, it can be used to "pre-46 align" the image sets on the basis of features such as bony anatomy, or external or internal fiducial marker systems, if available. 2.1.1.2 Affine Transformation Affine registration has two additional transformation types compared to rigid-body. These are scaling and shearing (also referred to as skewing) (Figure 2.1). Therefore, affine registration has a total of twelve degrees of freedom in 3D. Scaling can be described in terms of three parameters, as follows: (𝑥′𝑦′𝑧′1) = (𝑠𝑥0000𝑠𝑦0000𝑠𝑧00001) (𝑥𝑦𝑧1). (2.7) Shearing is slightly more complicated and can be defined in different directions. For example, a shearing of the x coordinate along the y-axis can be described as: Figure 2.1 Schematic of rigid, affine and non-rigid (deformable) transformations in 2D. 47 (𝑥′𝑦′𝑧′1) = (1000𝑠ℎ𝑥𝑦10000100001) (𝑥𝑦𝑧1). (2.8) The affine transformation matrix can still be defined in a single 4×4 matrix that can be applied to all the voxels of the moving image 𝑻𝑎𝑓𝑓𝑖𝑛𝑒(𝒑) = (𝑎11 𝑎12 𝑎13 𝑎14𝑎21 𝑎22 𝑎23 𝑎24𝑎31 𝑎32 𝑎33 𝑎340 0 0 1) (𝑥𝑦𝑧1), (2.9) where the 𝑎𝑖𝑗 elements are complicated combinations of all the transformation parameters. Affine registration can be used in the presence of geometric image distortions. In CT images, for example, scaling and skewing distortions can occur due to mis-calibration of voxel resolutions, gantry tilt with respect to the axis of bed, etc. [47] Affine registration can also be useful for rough alignment of inter-subject brain images. However, accurate and robust inter-subject registrations and also intra-subject registrations of areas of the body with deformable organs, such as abdominal and pelvic organs, require more degrees of freedom. 2.1.1.3 Deformable Transformation Deformable transformation, or deformable image registration (DIR) [48, 49, 50], has been applied to registration of deformable structures in the body. Examples include lung, abdominal (e.g. liver) or pelvic organs (e.g. bladder and rectum) [51, 52, 53], creation of human brain and lung atlases [54, 55], tracking soft tissue deformation during surgery [56], and numerous other medical applications. Deformable or non-rigid transformation involves many more degrees of freedom compared to rigid and affine transformations. The transformation is such that it does not preserve parallel 48 lines. Some non-rigid transformation models can be described using a single matrix. For example, a quadratic transformation model can be defined by second-order polynomials: (𝑥′𝑦′𝑧′1) = (𝑎00 𝑎01 … 𝑎09𝑎10 𝑎11 … 𝑎19𝑎20 𝑎21 … 𝑎290 0 … 1) (𝑥2𝑦2⋮1), (2.10) with 30 degrees of freedom. However, polynomial transformations can only model global deformations. For local anatomical variations, more degrees of freedom and local control is required and more complex transformation models are used. Deformable image registration is a very complex and vast topic in itself and therefore, a complete and in-depth analysis of the DIR methods and algorithms is out of scope of this dissertation, and this section will be limited to a brief discussion of selected, commonly used DIR methods, including linear elastic, viscous fluid and radial basis functions models. Other important methods not included are optical flow and finite-element-method (FEM) using mechanical models. For more comprehensive and technical coverage of the subject of DIR, the reader is referred to a large body of literature on this topic, including review articles, book chapters and books, such as Hill et al [47], Rueckert [50], Sotiras et al [57], Brown [58], Modersitzki [59, 60], and Shackleford et al [61]. 2.1.1.3.1 Linear Elastic Model Registration using an elastic model was introduced by Bajcsy et al [62] for aligning human brain CT images. The basis of linear elastic registration is to model the transformation of the moving image to the fixed image as the deformation of an elastic body, such as rubber. The behaviour of the elastic object is governed by the Navier linear elastic equation: 49 𝜇𝜵2𝒖(𝑥, 𝑦, 𝑧) + (𝜆 + 𝜇)𝜵(𝜵. 𝒖(𝑥, 𝑦, 𝑧)) + 𝒇(𝑥, 𝑦, 𝑧) = 0, (2.11) where 𝛁 and 𝛁2 are the gradient and Laplacian operators, 𝐮(x, y, z) is the displacement field, and 𝐟(x, y, z) is the external force on the elastic body. The parameters λ and 𝜇 are elasticity constants describing the material properties of the elastic object and are defined in terms of the Young modulus and Poisson ratio. The external force is what drives the registration and there are different choices for the external force in different models, such as the gradient of local correlation of intensities [62]. 2.1.1.3.2 Fluid Models The viscous-fluid model was proposed by Christensen et al [63] for large and highly localized deformations while preserving continuity and smoothness. In the Euclidian reference frame, the deformation of viscous fluid is described by the Navier-Stokes equation: 𝜇𝜵2𝒗(𝑥, 𝑦, 𝑧) + (𝜆 + 𝜇)𝜵(𝜵. 𝒗(𝑥, 𝑦, 𝑧)) + 𝒇(𝑥, 𝑦, 𝑧) = 0, (2.12) similar to Eq. 2.11 except that the differentiation is based on the velocity field, 𝐯(x, y, z), and λ and 𝜇 are viscosity constants. Equation 2.12 was originally solved by Christensen et al using the successive over-relaxation (SOR) method [64], which is very slow and computationally heavy for 3D image volumes. A faster method was suggested by Bro-Nielsen [65] by using linear elastic deformation of the velocity field of the fluid and applying a convolution filter in a scale-space framework. 2.1.1.3.3 Radial Basis Functions A commonly used representation of deformable transformation is through linear combination of basis functions Θ𝑖(𝑥, 𝑦, 𝑧): 50 (𝑥′𝑦′𝑧′1) = (𝑎00 … 𝑎0𝑛𝑎10 … 𝑎1𝑛𝑎20 … 𝑎2𝑛0 … 1) (𝛩1(𝑥, 𝑦, 𝑧)𝛩2(𝑥, 𝑦, 𝑧)⋮1) (2.13) One group of basis functions is polynomials (as shown in Eq. 2.10, which are not commonly used due to their limitation for soft tissue deformation applications and because they can introduce artifacts such as oscillations [66]. Fourier (trigonometric) and wavelet basis functions, which are both orthonormal, are a common choice for deformable transformations [67, 68]. Another commonly used group is radial basis functions (RBF), in which the basis functions Φ𝑖(𝑟) are defined based on the distance from the origin. The family of splines is particularly useful for this purpose. The term “spline” is based on an analogy with long thin strips of metal or wood, which are flexible but have some resistance toward bending. In spline-based registration, a set of corresponding “control points” in the moving and fixed images is defined. The spline-based transformation finds the displacement between the moving and target control points to define a smooth displacement field that also describes the correspondence at locations away from these points, using interpolation. The control points can be anatomical, fiducial, or geometrical landmarks. If no landmarks are available, the control points can be defined using a uniform mesh across the images [69], and are called quasi- or pseudo- landmarks since they are not actual corresponding landmarks. The most common types of spline-based models are thin-plate-splines and B-splines. In the thin-plate-spline (TPS) method [70], the deformation of an infinite thin metal plate is modelled. In the original TPS formulation, a flat plate in the X-Y plane is exposed to known 51 vertical (along Z) displacements at particular control points in the plane causing the plate to deform. The vertical displacement 𝑡 at any location (x,y) in the plane, given the displacements at the control points i = 1 … n, is interpolated as a linear combination of radial basis functions: 𝑡(𝑥, 𝑦) = 𝑎1𝑥 + 𝑎2𝑦 + ∑ 𝑏𝑖𝑟𝑖2𝑙𝑛𝑟𝑖2𝑛𝑖 (2.14) where 𝑟𝑖 is the distance from location (x,y) to the control point i. These displacements are constrained so that the plate deforms in a way that minimizes its total bending energy, ∬ ((𝜕2𝑡(𝑥,𝑦)𝜕𝑥2)2+ 2 (𝜕2𝑡(𝑥,𝑦)𝜕𝑥𝜕𝑦)2+ (𝜕2𝑡(𝑥,𝑦)𝜕𝑦2)2) 𝑑𝑥 𝑑𝑦𝑅2. (2.15) In the context of TPS for deforming points in 3D space (as opposed to a flat plane), when 3D displacements are applied to control points in 3D (as opposed to vertical displacements), each of the three components (𝑑𝑥, 𝑑𝑦, 𝑑𝑧) of the displacement vector at any point (x,y,z) is interpolated in a manner as in equation (2.14), where 𝑟𝑖 becomes the distance in 3D between (x,y,z) and the i-th control point. In the TPS model, each basis function (i.e. each control point) has a global effect on the transformation so that varying the position of one control point also varies the positions of all the other points in the moving image. Increasing the number of control points significantly increases the computational cost associated with moving one control point. Therefore, TPS method might not be suitable for modeling very complex and localized deformations. The B-spline method, on the other hand, uses locally controlled functions. This type of deformation is also called free form deformations (FFD), which is widely used for medical DIR tasks [71, 72, 73]. In spline-based FFD, a regular mesh of control points with uniform spacing is 52 used, rather than the arbitrary arrangement of control points in RBF splines. The B-splines are only defined in the vicinity of each control point and varying the position of one control point only changes the transformation in the neighbourhood of that point. Therefore, the B-spline method is known to have “local support” and it is very popular for modeling deformations in the chest, heart, liver, breast, and brain, amongst other sites. 2.1.2 Similarity Metrics The degree of similarity between the two images is measured in terms of a similarity metric. Intensity-based registration uses a similarity metric defined in terms of the image voxel intensity values. The similarity measure can also be defined based on the geometric features in the image. These feature-based registrations include functionally important points (i.e. landmarks) or surfaces in the image and are, therefore, known as landmark-based and surface-based registration methods, respectively. 2.1.2.1 Landmark-Based Registration Landmark-based registration [74, 75, 76] uses a set of corresponding points, i.e. landmarks, in the moving and fixed images. In this method, the optimal spatial transformation (T) is found that minimizes the root squared distance (RSD) between n corresponding landmarks, calculated as 𝑅𝑆𝐷({𝒑𝑖}, {𝒒𝑖}) = √∑ [𝒒𝒊 − 𝑻(𝒑𝑖)]2𝑛𝑖=1 , (2.16) where 𝒑𝑖 and 𝒒𝑖 are the corresponding landmark points in the moving and fixed images, respectively. These corresponding points may be fiducial markers or anatomical landmarks identified in each image by the user. Most of the landmark-based registrations are either rigid or affine, but some DIR techniques also incorporate landmark information [77]. 53 2.1.2.2 Surface-Based Registration In surface-based registration [78, 79, 80], geometrical boundary features extracted from images in the form of edges, lines, or surfaces are used for the alignment. The feature extraction can be done through automatic or semi-automatic image segmentation methods [81], such as region growing [82] or thresholding [83], and is itself subject to error. The transformation is optimized for some measure of the distance between the corresponding lines/surfaces [84, 85]. The original and widely used surface-based method is called the “hat and head” algorithm [78], which was developed to align images of the brain acquired with different modalities. In this method, surface contours from one of the images form the “head” and a set of points that correspond to the same surface in the other image form the “hat”. The algorithm then tries to “fit a hat to head” through minimization of the root-mean-square distance (RMSD) between the points in the hat and the head surfaces, through progressive refinements (i.e. iterations). The “hat and head” method is prone to error. There are many ways to actually “fit a hat to head”, so there are often multiple alignments that produce similar RMSD values, as shown in a simple 2D illustration in Figure 2.2. Other feature-based methods are also currently used. Most are based on the iterative closest point (ICP) algorithm [86], which will be discussed in more detail in Chapter 3. Surface-based methods can incorporate either linear or non-linear transformation models. They can also combine landmark information in addition to the surfaces [80]. 54 Figure 2.2 A 2D illustration of different “hat and head” solutions for a feature-based registration, yielding similar RMSD distances, while only one is a perfect match. 2.1.2.3 Intensity-Based Registration The most commonly used similarity measure in image registration, especially for more complicated problems, is intensity of image voxels. Intensity-based registration is straightforward for intra-modality registrations, but can also be applied to inter-modality images if an appropriate similarity metric is defined. Various similarity metrics have been used for intensity-based registration, the most common of which are sum of squared differences, correlation coefficient, and mutual information. For N number of voxels in the image overlap region (i.e. Ω𝑚,𝑓 in Eq. 2.4), the sum of squared differences (SSD) is defined as 𝑆𝑆𝐷(𝐼𝑚, 𝐼𝑓) =1𝑁∑ [𝐼𝑚(𝒑𝑖) − 𝐼𝑓(𝑻(𝒑𝑖))]2𝑁𝑖=1 [87]. (2.17) A perfect alignment is achieved when this measure is zero (in the absence of noise). This definition of SSD is based on the assumption that the images have the same intensity characteristics and therefore, it is not suitable for inter-modality registration problems. Another measure for intra-modality registration is correlation coefficient [88], which assumes 55 a linear relationship between the image intensities. Correlation coefficient (CC) is defined by where 𝜇𝑚 and 𝜇𝑓 are the average voxel intensities in the moving and fixed images, respectively. The correlation coefficient metric is more flexible than SSD, but its applications are still widely limited to intra-modality registrations. The most commonly employed inter-modality intensity-based measure of alignment is mutual information, introduced by Viola and Wells [89] and also by Maes [90]. This method makes use of the information content or entropy of the images. In information theory, the information content of images 𝐼𝑚 and 𝐼𝑓 is measured by the Shannon-Wiener entropy [85], 𝐻(𝐼𝑚) = − ∑ 𝑃(𝐴)𝑙𝑜𝑔𝑃(𝐴)𝐴 (2.19) and 𝐻(𝐼𝑓) = − ∑ 𝑃(𝐵)𝑙𝑜𝑔𝑃(𝐵)𝐵 , (2.20) where 𝑃(𝐴) is the probability that a voxel in image 𝐼𝑚 has intensity A and 𝑃(𝐵) is the probability that a voxel in image 𝐼𝑓 has intensity B. The joint entropy of 𝐼𝑚 and 𝐼𝑓 is then given by 𝐻(𝐼𝑚, 𝐼𝑓) = − ∑ ∑ 𝑃(𝐴, 𝐵)𝑙𝑜𝑔𝑃(𝐴, 𝐵)𝐵𝐴 , (2.21) where 𝑃(𝐴, 𝐵) is the joint probability that a voxel in the overlapping region Ω𝑚,𝑓 has the intensity of A and B in images 𝐼𝑚 and 𝐼𝑓, respectively. The image alignment can then be determined using this joint entropy 𝐻(𝐼𝑚, 𝐼𝑓), by defining a mutual information (MI) similarity measure 𝐶𝐶(𝐼𝑚, 𝐼𝑓) = ∑(𝐼𝑚(𝒑𝑖) − 𝜇𝑚)(𝐼𝑓(𝑻(𝒑𝑖)) − 𝜇𝑓)𝑖 , (2.18) √(∑ 𝐼𝑚(𝒑𝑖) − 𝜇𝑚𝑖)2(∑ 𝐼𝑓(𝑻(𝒑𝑖)) − 𝜇𝑓𝑖)2 56 𝑀𝐼(𝐼𝑚, 𝐼𝑓) = 𝐻(𝐼𝑚 ) + 𝐻(𝐼𝑓 ) − 𝐻(𝐼𝑚, 𝐼𝑓 ), (2.22) which must be maximized in the registration process. 2.1.3 Optimization Process A cost function (aka energy function, E) is defined based on the similarity measure. In the case of affine and deformable transformations, a regularization term is typically included in the cost function, in order to assign a high cost to undesirable transformations, such as high local stretching or bending, based on a priori knowledge about the expected deformations [73, 67]. The general form of a cost or energy function can be written as 𝐸 = −𝐸𝑠𝑖𝑚𝑖𝑙𝑎𝑟𝑖𝑡𝑦 + 𝐸𝑖𝑟𝑟𝑒𝑔𝑢𝑙𝑎𝑟𝑖𝑡𝑦. (2.23) This cost function is minimized during the optimization process by searching over the parameters space defined by the transformation model, while constraining the search using the regularization criteria. Regularization is particularly important in deformable registration with numerous degrees of freedom available to the transformation model. In most deformable registration algorithms, the degree of regularization can be controlled by the user using a regularization parameter (often called regularization lambda, λ). The optimization is usually an iterative process, where an initial estimate (starting point) is gradually refined through multiple iterations, until a certain threshold is met. In each iteration, an estimate of the similarity measure (or the cost function) is calculated based on the current transformation. Then, the optimization algorithm finds another (ideally better) estimation of the transform and the corresponding similarity measure is determined again. This process is continued toward convergence, which is the global minimum of the cost function. 57 Finding this global minimum is especially hard for deformable registration tasks. Therefore, a finite threshold for the similarity metric (or cost function) or number of iterations serves as the “stopping criterion” for the iteration process of optimization. A common challenge for complicated non-rigid registration problems is “local minima”, where the optimizer gets trapped in a local minimum of the cost function and cannot converge toward the globally optimum solution. Different studies have attempted to address the “local minima” problem [91, 92, 85] and possible solutions include multi-start (i.e. multiple starting points) and multiscale or multi-resolution (i.e. using multiple sub-samplings of the original images) optimization methods, which provide at least multiple local minima if not the global minimum. 2.2 Shape Correspondence While image registration is the method of choice for many clinical and research applications, it has some limitations. Because image registration relies directly on the image data, its applicability and accuracy suffer when image quality is low, or in the presence of significant image artefacts. Also, there are scenarios where “non-corresponding” structures, which can compromise the registration, are present in the pair of images. For example, in HDR Gynaecological brachytherapy, there might be variation in the anatomical position of the HDR applicators or the Foley balloon inside the bladder; one structure in an image might be missing in the other image, such as air pockets in the bladder or rectum, or, in extreme cases, the size and/or type of HDR applicator might be different (Figure 2.3). The registration algorithm tries to match these “non- corresponding” structures, and yields inaccurate outcomes as a result. Furthermore, DIR usually cannot handle large local deformations in structures close to bony anatomy, such as the bladder (Figure 2.3). 58 Considering that the primary goal is to establish correspondence in the anatomical structures in the matching images, a well-suited alternative is Shape Correspondence [93]. Shape correspondence, or structure-based registration, is the process of establishing a point-correspondence between two defined 2D or 3D shapes. In imaging applications, the shapes are defined by outlines, or "contours", of the anatomical structures of interest in the image data sets. These contours are usually generated manually, although high contrast structures are amenable to the application of auto-segmentation algorithms. The structures can be represented as surfaces, volumes, or point sets. The process of aligning two sets of points is called point-set registration. Point-set registration can be used when seeking correspondence between two 2D or 3D point clouds. It can be either rigid, affine or deformable, similar to image registration. An example of a 2D point-set registration is illustrated in Figure 2.4. Many different algorithms and methods have been developed for point-set registration, such as iterative closest point (ICP), thin-plate-spline (TPS), coherent point drift (CPD), etc. These techniques are discussed in detail in Chapter 3. Point-set registration can easily be substituted for image registration when the correspondence is only required for certain ROI in the image volume. For example, in image-guided radiation therapy, it is often required to match contoured OAR in different sets of images of a patient, mostly for dose assessment purposes. While deformable image registration can be used for this task, for the reasons explained above, it may not render accurate results. In such situations, deformable point-set registration might be used with superior results. 59 Figure 2.3 Example of possible non-anatomical mismatch between two brachytherapy planning CT images, showing air pocket in bladder (yellow contours) and rectum (orange contours) present only in one image. The brachytherapy applicator in the left (tandem-ring - green and magenta markings) is also different from the one in the right image (tandem-ovoids - green markings). In addition, a high degree of deformation in bladder is present close to the non-moving bony anatomy (cross section of the pelvic bones as shown). 2.3 Validation of Registration Outcomes Registration methods must be assessed by the quality and accuracy of their output. Validating registration outcomes is crucial in medical applications, yet sometimes difficult due to lack or scarcity of reliable ground truth data. Nevertheless, a number of methods to assist in evaluating the quality or accuracy of the registration outcomes have been proposed. Most of these methods are applicable to both image and shape registrations and some are only applicable to image registration. These will be discussed in each context accordingly. 2.3.1 Visual Inspection A qualitative approach to assess registration results is visual inspection, where the user visually compares the registered (i.e. deformed) image with the target (fixed) image in 2D or 3D. 60 This validation method can be used by itself for simple registration tasks (such as for brain images) or as an initial check in the validation process. For more complicated image alignment tasks, some methods are available to facilitate the visual evaluation process, such as overlay (using image transparency), checkerboard (i.e. alternating boxes of deformed and fixed images), and image subtraction. While visual inspection is more feasible and meaningful for image registration, it can still be useful for shape correspondence given enough available features in the shapes, such as in Figure 2.4. For shape correspondence, other than feature matching, the distance between the surfaces of the target and deformed structures can be visually assessed as well. Figure 2.4 Sample of 2D deformable point-set registration, matching the blue fish to the red fish using the coherent point drift (CPD) algorithm (adapted from [94]). Both point-sets have missing points and spurious outlier points. 2.3.2 Numeric Validation Methods For complex registration problems, especially when using deformable registration, more analytic and robust validation is required. However, validating complex registrations is also more difficult. Therefore, in a lot of cases, multiple methods and metrics are applied. If some ground truth data are available, the registration accuracy can be assessed using metrics such as target 61 registration error and Dice similarity coefficient. In the absence of ground truth correspondence, quantitative validation of the registration performance can also be done with what are termed robustness and consistency tests. Robustness and precision in the registration outcomes is established by testing the bias and sensitivity of the results after adding random noise or choosing different priori or parameter settings. Consistency checks evaluate the performance of a registration algorithm in circular transformations. For example, the deformation vector field in DIR can be tested using an inverse consistency test between image A and B, where A is deformed to match B and B is separately matched to A. An inverse-consistent algorithm will produce a true inverse deformation field when the moving and fixed images are switched. An inverse consistency error (ICE) between the forward and inverse registration is defined as a measure of this quality [95, 96]. Other approaches to validating registration outcomes include gold standard transformation [76], biomechanical simulations [97], segmentation overlap [98], etc. There are also voxel/point-wise metrics, such as Hausdorff distance and root-mean-square error, that are mainly used for evaluating registration of volumes, surfaces, or point clouds. 2.3.2.1 Target Registration Error One of the most popular and systematic validation methods is based on a quantity called target registration error (TRE) [99]. Target registration error is only efficient when a) there is a set of relevant matching landmarks (anatomical or fiducial) accurately localized in both moving and fixed images (or structures), and b) these landmarks were not used in the registration process itself (e.g. in a landmark-based registration, they are not the same landmarks that were used for registration). 62 For a pair of matching landmarks 𝒑𝑙 and 𝒒𝑙 in the moving and fixed images (or structures), respectively, TRE is calculated by 𝑇𝑅𝐸 = ‖𝒒𝑙 − 𝑻(𝒑𝑙)‖, (2.24) i.e. the spatial distance between the landmark 𝒒𝑙 and the registered matching landmark position. Although TRE has the advantage of utilizing “ground truth”, it evaluates the registration accuracy only at certain points and not on a point-by-point basis. Therefore, without a uniform and dense distribution of the landmarks in the image (or structure) volume, it might misrepresent the accuracy of the registration results for some areas of the matched volumes where there are no landmarks available. Moreover, reliable landmarks might be limited or not available at all for some applications. 2.3.2.2 Dice Similarity Coefficient Dice similarity coefficient (DSC) was originally introduced by Dice [100] and can be used as a measure of similarity or match between two volumes, which could be binary segmentations, ROI volumes, etc. It is a measure of spatial overlap, where the overlap volume of the corresponding objects is determined in the transformed and fixed data sets and then normalized by the average of the volumes in the two data sets. If we define volumes 𝑉𝑚 and 𝑉𝑓 as a set of pixels in the moving and fixed data sets, the Dice coefficient can be defined as 𝐷𝑆𝐶 = (𝑉𝑚′∩𝑉𝑓)(𝑉𝑚′+𝑉𝑓)/2, (2.25) where 𝑉𝑚′ is the set of pixels in the deformed moving volume. The Dice coefficient may take any values from 0, indicating no overlap, to 1, being a complete overlap between the images, which are both unlikely to occur. In image data sets, DSC can be 63 applied to evaluation of auto-segmentation algorithms, and can also be used for evaluating image registration methods if matching objects are identified in each image either through auto-segmentation or user-defined segmentation. In such applications, the reliability of the result will depend on the accuracy of the segmentations. 2.3.2.3 Hausdorff Distance Hausdorff distance (HD) [101] is determined by calculating the distances between each point in the registered structure and the closest point in the target, and taking the maximum of these. This metric is mostly used for validation of shape correspondence and auto-segmentation methods, but similar to DSC, it can be applied to image registration using ROI contours or segmentations. The HD metric provides a measure of the worst-case-scenario for the registration uncertainty, which could be over-estimated in the presence of outliers. 2.3.2.4 Root-Mean-Square Error The root-mean-square error (RMSE) is a measure of the overall agreement between the deformed and target structures. It is defined by the mean distance between the transformed points to the corresponding (i.e. closest) points in the target structure. This measure is also mostly applicable to shape correspondence rather than image registration. 2.4 Project Overview The main hypotheses of this study were as follows: 1. Dose accumulation in bladder using deformably registered data sets gives substantially different results than DVH-based dose accumulation methods. 2. Changing the bladder volume across HDR ICBT treatment fractions may vary the location and 64 physical separation of hotspots on the bladder surface from fraction to fraction, and thereby potentially reduce the rate of associated urinary toxicity in cervical cancer patients treated with multi-fraction HDR ICBT. 3. Cumulative doses calculated in deformably registered data sets are superior to DVH-based dose accumulation methods in predicting urinary toxicity in cervical cancer patients treated with multi-fraction HDR ICBT. In order to test these hypotheses, it is crucial to identify a suitable and accurate deformable registration technique to track every point on the bladder between treatment fractions. This is a complicated task due to the high inter-fraction variability of bladder size and shape. In Chapter 3, a selection of deformable registration methods, including both DIR and point-set registration algorithms and software, are thoroughly evaluated and tested using a porcine bladder phantom. Then, the locally accumulated dose to the bladder (and bladder neck) over the course of treatment fractions must be determined. In Chapter 4, the superior point-set registration technique (CPD), based on the outcomes of Chapter 3, is used retrospectively on the treatment planning data of a cohort of 60 cervical cancer patients treated at BCCA. The bladder wall was carefully re-segmented and the ureteral and urethral points of attachment to the bladder were identified and used as an additional quality control measure of the registration outcomes. The bladder-wall contour point clouds were registered using CPD and the accumulated bladder-wall dose was determined. A set of cumulative dose parameters was then defined for the bladder wall and bladder neck. The association of the bladder volume variation with the movement of the bladder hotspots and urinary toxicity was also investigated. Urinary toxicity data collected using LENT-SOMA subjective questionnaires were analyzed 65 to define both aggregated and endpoint-based toxicity grades. These grades are studied against the cumulative parameters, in Chapters 4 and 5, in order to investigate the association of urinary complications with localized dose to bladder wall and bladder neck. Also, the correlation between bladder neck and bladder wall dose parameters was studied in Chapter 5. The conclusions and future directions based on this study are presented in Chapter 6. 66 Chapter 3: Validation of Non-Rigid Point-Set Registration Methods 3.1 Introduction Since the introduction of three-dimensional (3D) medical imaging in radiation therapy more than three decades ago [102, 103, 104] modern radiotherapy treatment methods, whether external beam radiotherapy (EBRT) or brachytherapy, rely increasingly on 3D imaging for treatment planning, delivery, and quality control. For most radiotherapy applications, 3D patient images must first be annotated by identifying and contouring regions of interest (ROI) such as the clinical target volume (CTV) and different organs at risk (OAR). The dose received by each structure can then be modeled using sophisticated dose calculation algorithms. However, while the precision of these algorithms continues to improve, dose calculations are typically limited to a single, pre-treatment image set used for planning purposes. These planned doses are not completely representative of the dose delivered by the treatments that follow, due to variations in patient set-up and anatomy during the course of treatment. At some institutions, a new set of planning images is obtained at each fraction of certain radiotherapy treatments, such as cervical cancer brachytherapy. This allows the accumulated dose to different ROI to be estimated by summing the dose values calculated from the different image sets. Nevertheless, the summed dose remains a poor representation of the actual accumulated dose to individual tissues, especially for deformable pelvic OAR, such as bladder, rectum, and bowel, which can have a significantly different size, shape and orientation in treatment fractions separated by one or more days. A popular solution to the problem of dose accumulation for deformable structures in radiation therapy is to use deformable image registration (DIR) when multiple intra-patient image sets are 67 available. DIR matches the coordinates of a source (deforming, secondary or moving) image to a target (reference, primary or fixed) image, providing a spatial transformation map. For the purpose of dose accumulation in EBRT for example, multiple treatment verification image sets in the form of cone beam CT’s, with accompanying calculated dose maps, can be registered to a single target image set (such as the planning CT) using DIR. The resulting spatial transformation maps can then be used to map the cumulative dose onto the target image, and provide a more accurate estimate of the actual delivered dose. Many DIR techniques have been developed over the years with different algorithms and applications [57, 105]. However, the accuracy of DIR is directly dependent on the image quality and is compromised when there are image artifacts. Moreover, in cases such as gynecology brachytherapy treatment for cervical cancer, there are non-matching structures, such as the brachytherapy applicators, Foley balloon, air bubbles, etc, In this milieu, DIR is considerably challenged, as we verified during the course of testing a range of DIR methods, including DIR software packages (MIM Maestro (MIM Software Inc. Cleveland, OH), Elastix1, MIRT2, Plastimatch3), algorithms (Demons, B-spline, etc), and similarity measures (Mutual Information, Sum of Squared Differences, etc), for the task of registering the bladder wall in intra-patient gynecology brachytherapy planning CT scans taken for different treatment fractions. We found this registration task to be very complex because of (i) the high degree of local deformation in the pelvic organs, (ii) the presence of the surrounding non-matching structures, as mentioned, and (iii) occasional presence of low image quality and image artifacts. Given this complexity, even with substantial tweaking of DIR parameters and settings, 1 http://elastix.isi.uu.nl 2 https://sites.google.com/site/myronenko/research/mirt 3 http://www.slicer.org/slicerWiki/index.php/Documentation/4.3/Extensions/Plastimatch 68 none of these DIR methods produced satisfactory results, as assessed by visual inspection. In an attempt to improve the DIR results for matching bladder, some of the surrounding structures were masked out of the images to negate their complicating effects on the DIR’s similarity metric calculation. This required extra data preparation effort and the final results still did not pass visual inspection. Therefore, we considered the alternative approach of aligning the ROI contours, which may be more suitable for highly elastic ROI, such as bladder, since the dose data are only required for the ROI structures, rather than the whole image. In this approach, the correspondence is sought between two 3D shapes, which can be represented by 3D point-clouds. The applicability of point-set correspondence is similar to DIR, with the difference being that point-set correspondence provides a transformation map between two point clouds that represent ROI contours, rather than the mapping between image voxels. Furthermore, these point clouds can, to a good approximation, be converted to image voxels to allow voxel-based dose accumulation. Shape correspondence [93] is particularly useful when there are image artifacts and/or surrounding non-matching structures that can confuse a DIR algorithm. Although manual intervention, and with it some degree of subjectivity, is generally required for creation of the ROI contours, the registration results are completely independent of the image data, so they are not affected by image artifacts or other features of the image itself. When using a deformable image or shape registration method, validating the outcome is a challenging yet crucial step, especially when (i) the underlying actual (‘ground truth’) anatomical correspondence is unknown and/or (ii) the objects lack distinctive anatomical landmarks (or such landmarks are crowded into a small localized region). Most OAR in the pelvic region are elastic structures. The bladder is one example. It has a rather uniform bean-like shape with only three 69 roughly identifiable landmarks: the points of attachment of the two ureters, and the urethra. Therefore, it is both important and difficult to validate the registration outcome for such a structure. The purpose of this study was to evaluate the performance of five different non-rigid point-set registration methods for registering bladder surface contours, with regard to the algorithms’ robustness, computation speed, and most importantly, the registration accuracy. The methods tested represent a range of non-rigid point-matching paradigms, and were selected based on their applicability to the given problem. One of the methods is based on affine registration, while the rest are deformable registration methods, which do not preserve parallel lines. All the selected algorithms were available as open-source implementations. In order to benchmark a currently available commercial package against methods presented here, we also used the DIR tool in the MIM Maestro software. A pelvic phantom enclosing a pig bladder with fiducial markers attached to it was created for validation purposes. CT scans of the phantom were acquired with the bladder at various filling volumes, in order to simulate the possible range of deformations present in multi-fraction radiotherapy in the pelvic region. 3D point-sets were extracted from the contours of the bladder surface in each scan. The pig bladder was selected as a surrogate of the human bladder specifically because it was relatively easy to acquire, and because its biomechanical and physiological properties are documented to be similar to those of the human bladder [106, 107, 108]. To the best of our knowledge, there have not been any studies that compare multiple point-set registration techniques by validating the results on a deformable structure lacking discriminative features. A relevant study has been reported by Hu et al [109], who compared the performance of their patient-specific statistical shape model (SSM) to two elastic registration methods (Coherent Point Drift (CPD) and Robust Point Matching (RPM), both of which were also 70 tested in our study), for surface-based deformable registration of prostate in trans-rectal ultrasounds. The accuracy of SSM was compared to the other two methods based on anatomical landmarks that were identified on the ultrasound images, using metrics such as surface distance and target registration errors. While it was found that a statistical motion model slightly outperformed the alternative deformable registration methods, CPD provided better registration accuracy compared to RPM. In our study, our aim was to evaluate methods for non-rigid registration of the bladder, which has a much higher degree of deformation and size variation than the prostate, rendering its registration more difficult. Therefore, we used fiducial landmarks, which, in contrast to anatomical landmarks, can be made more numerous, and can be placed in carefully selected positions on the object of interest. They are also typically much less subject to false identification, and can be localized with a high degree of certainty (i.e. low fiducial localization error). Therefore, they permit very accurate measurements of the target registration error (TRE) [99], which is a measure of the spatial accuracy of a registration. In a closely related study using fiducial markers on a porcine bladder phantom, Wognum et al [110] compared a variety of DIR techniques and one point-set registration method: Symmetric Thin Plate Spline Robust Point Matching (STPSRPM) (Bondar et al [111]), which is very similar to one of the methods evaluated in our study. Their study showed better performance of STPSRPM compared to the DIR methods. Our study is complementary to theirs and further demonstrates that the flexibility of deformable point-set matching algorithms makes them a suitable candidate for registering bladder contours and similar elastic structures. A slightly different approach, employed by Kirby et al [112], used a carefully constructed deformable phantom that closely resembled a single axial slice of pelvic anatomy to test the 71 performance of eleven different DIR methods (including MIM, which was also used in our study). This phantom was used to simulate 3D deformations, which, due to its 2-dimesional geometry, could be far from the actual deformations observed in human anatomy. The ability of this phantom, therefore, to mimic the large localized 3D deformation observed in the pelvic organs, such as bladder and rectum, was limited. Moreover, the amount of variation in the size, shape and position of the pelvic organs present in the two modes of the deformable phantom (i.e. empty and full bladder modes) was minimal compared to the actual scenario. Another related study was reported by Varadhan et al [113], where the accuracy of B-spline and diffeomorphic demons DIR algorithms was tested on geometrically deformed image data, as rendered by ImSimQA software. Although simulated phantoms have their benefits, the end results could be further from reality compared to a physical phantom that uses structures similar to the human organs (such as the porcine bladder pelvic phantom used in our study). 3.2 Point-Set Registration Methods In this section, the point-set registration methods evaluated in our study are described. Table 3.1 gives a high-level summary of the important attributes of each method. The chosen techniques were selected based on their acclaimed robustness, popularity, and applicability to 3D data. Furthermore, the selected methods represent a wide range of methodological paradigms (i.e. two methods were not both chosen if they were too similar). Although this study is not assessing the feasibility of these software packages for clinical use, the accessibility and usability of each software package is important for prospective studies like this one, as well as potential future application to clinical data sets. 72 Table 3.1 Summary of the components of each point-set registration technique evaluated in this study. Method Shape Representation Correspon-dence Model Cost Function Deformation Model / Regularization Optimization ICP-finite original ICP by Besl and McKay [86] 3D points Binary Euclidean distance Affine transform Iterative TPS-RPM, Chui and Rangarajan [91] 3D points Fuzzy Euclidean distance Thin-plate splines Deterministic annealing SC, Belongie et al. [114] Log-spherical histogram (shape context) Binary 𝜒2 histogram distance B-spline freeform deformation Hungarian algorithm CPD, Myronenko and Song [115] 3D points and Gaussian Mixture Model (GMM) Implicit Expectation of the log likelihood function Motion coherence Expectation-maximization (EM) GMM-reg, Jian and Vemuri [116] GMM Implicit L2 distance between distributions Thin-plate splines Iterative scale reduction with gradient descent sub-step 3.2.1 Iterative Closest Point Iterative closest point matching (ICP) was originally introduced by Besl and McKay [86] as a standard technique for establishing correspondence between two point-clouds. Since then, the ICP method has been used and/or cited in numerous studies in different fields, particularly medical image analysis applications [117, 118] and range images [119, 120] (i.e. 2D images with pixel values describing a 3D scene in terms of the distance from a fixed point, such as a camera or sensor device). The algorithm works by iteratively seeking a rigid transform to minimize the mean squared distance between pairs of closest points in the two point-sets. In each iteration, a correspondence between closest points in the two point-sets is found, and the distance between corresponding points is minimized. 73 The transformation in the original ICP is limited to translation and rotation only. The algorithm does not usually converge exactly to the global optimum, so the correspondence found is not always one-to-one. The algorithm is also ineffective in the presence of gross noise and outliers [86, 121]. Despite these drawbacks, ICP continues to be a frequently cited approach for seeking correspondence between 3D shapes, mainly because the algorithm is simple to understand and straightforward to implement, and also it can be modified to fit a task and/or be used in conjunction with other algorithms. Over the years, multiple modifications of the original method have been put forward to improve its performance [122, 123, 120, 124]. There are also studies that review and compare these ICP-variant algorithms [121, 125]. Due to the popularity of ICP for point-set and free-form curve and surface matching, it was included in this study even though it is not a deformable registration method. However, because of the large amount of deformation present in the shapes to be registered and anticipated poor performance of standard ICP, we used the ICP-finite method by Kroon [126], which is an improved version of the original rigid ICP. ICP-finite is an affine transformation since it incorporates scaling and shearing in the transformation model by using the finite difference method (FDM), proposed by Grossmann et al [127]. 3.2.2 Thin-Plate-Spline Robust Point Matching Robust point matching (RPM), introduced by Gold et al [128], is a deformable point-set registration method based on fuzzy point correspondence. It relies on optimizing a cost function that simultaneously discovers (i) a correspondence between points in the two point-clouds and (ii) a deformation of the points in the moving cloud that minimizes the Euclidean distance between corresponding point-pairs. Correspondence is represented as a value between 0 and 1 for every possible pair of corresponding points, with 0 indicating no possible match and 1 a perfect match. 74 An entropy term encourages points to partially correspond to multiple partners. A weight parameter, T, (a.k.a temperature parameter) on this entropy term controls the overall softness or crispness of the correspondence and, through an iterative deterministic annealing process [128, 129], first encourages and then discourages the fuzzy correspondence between points. Annealing proceeds as follows: T starts at a pre-defined Tinit value and is gradually reduced at a linear rate, i.e. Tnew = Told.r, where r is the annealing rate. At each iteration, the cost function is minimized for a given temperature. This process repeats until a preset final temperature, Tfinal, is reached, at which the correspondence values converge to either 0 or 1, i.e. binary assignment. TPS-RPM by Chui and Rangarajan [91] is an RPM algorithm based on thin-plate spline (TPS) [130, 131] parameterization of the deformable transformation. In TPS-RPM, the registration regularization parameter is also driven by the parameter, T, so that in the initial iterations, where the temperature is high, the regularization is also very high and the transformation is mainly rigid. As the temperature drops, the regularization also decreases and the registration becomes more and more deformable. TPS-RPM is widely used (more than 500 citations) and has been applied to various registration applications [132, 133, 134]. Its performance has been previously tested and compared to different techniques, such as those discussed in sections 3.2.4 and 3.2.5. 3.2.3 Shape Context Shape context (SC) was first introduced by Belongie et al [114] as a technique to find one-to-one correspondences between two shapes represented as point clouds. In SC, each point is associated with a point descriptor, which is represented by a histogram of the position of its neighbouring points in log-spherical coordinates (i.e. spherical coordinates but with log of the radial component). To achieve scale invariance, the radial components are first normalized by the mean distance between every pair of points in each shape. 75 Each point descriptor characterizes the local shape appearance of the point-cloud at that point. To judge the similarity between two descriptors, the 𝜒2 distance is used, which compares differences between each corresponding pair of bins in the two histograms. A matching between two point-clouds is achieved by finding a one-to-one correspondence with the minimum summed -𝜒2 distances. This matching is found using the Hungarian algorithm [135], a polynomial time algorithm that finds the optimal assignment on a bipartite graph. Kroon [126] has used this shape context descriptor in combination with a b-spline freeform deformation grid [136] to register two matching shapes with regularization over the transformation. They validated their SC approach on a set of mandible (lower jawbone) surfaces, which has a rather complex structure, for ten different subjects. A freely available implementation of their method is available online. 3.2.4 Coherent Point Drift Coherent point drift (CPD), introduced by Myronenko et al [115], is a point-cloud matching method that, contrary to ICP, TPS-RPM and SC, uses a probabilistic model of point positions and searches for the best registration via maximum likelihood estimation (MLE) [137, 138]. CPD first initializes a Gaussian mixture model (GMM) [139, 140] with isotropic Gaussians centered at the points in the moving structure point cloud. The GMM parameters are then updated to maximize the likelihood of the target points given the moving model. The likelihood function is regularized by a term that encourages points close to one-another to move coherently. It works by penalizing high frequencies in the displacement vector field and has the effect of preserving the topological structure of the moving point-set. This kind of regularization is based on ‘motion coherence theory’ [141, 142]. The deformable registration 76 capabilities of CPD have been tested on various synthetic and real models [115], including 3D left ventricle contours segmented from ultrasound images. The registration performance of CPD has also been compared to TPS-RPM on point-sets of simple 2D models with missing points and/or different levels of deformation, noise, and outliers-to-data ratios [115]. The registration errors for these tests demonstrated superior performance of CPD in the presence of high levels of noise and outliers. 3.2.5 GMM Registration Jian and Vemuri [116] presented a GMM-based point-set registration framework, which, in contrast to CPD, represents both point-clouds as GMMs. Alignment is done by searching for a thin-plate spline deformation of the moving point-cloud that minimizes the L2 distance for distributions between the two GMMs. The scale of the spherical Gaussian fit to each point is first set to be a large value in order to align the point-clouds coarsely, ignoring small discrepancies and noise. Similar to the temperature parameter in the annealing process used for RPM, the scale is incrementally decreased and the points are realigned at each step. In Jian and Vemuri [116], the performance of this GMM-based registration was compared against CPD and TPS-RPM on 2D fish and 3D face datasets. It was found that all the methods produced close to perfect alignment. Therefore, they concluded that these algorithms perform efficiently for most deformable point-set registration tasks. 3.3 Materials and Methods This section provides a detailed description of the materials, procedures and software used to perform this study. In section 3.3.1, we describe the pelvic phantom components and setup. Section 3.3.2 includes details about CT-scanning of the phantom with the pig bladder at different filling 77 volumes. In section 3.3.3, the procedures for point-set registrations using the software libraries are discussed. Figure 3.1 shows an overall flowchart of the study design including the important data acquisition and analysis steps for point-set registrations. Section 3.3.4 reviews the MIM DIR method. The methods for point-set registration accuracy evaluation are discussed in section 3.3.5. 3.3.1 Phantom Preparation The bladder of a female pig was freshly acquired, with the urethra and the two ureters attached, and was used within a few hours. To create landmarks that were clearly identifiable and distinguishable, two types of fiducial markers (made of rubber and plastic) were glued onto the bladder surface at different locations, each with different groupings (Figure 3.2). The fiducial markers were distributed in a reasonably uniform pattern on the bladder surface. The total number of fiducial landmarks used was a trade-off between sufficient coverage of the surface area and the need to maintain, as closely as possible, free deformation of the bladder as it would occur in situ. A pelvic phantom was built to house and support the bladder, and to model typical human female pelvic anatomy. The phantom consisted of a plastic container into which were placed an orthopaedic pelvis bone model (Model 4064 Pelvis fx, SYNBONE AG, Malans, Switzerland) made of solid foam, a hand-made uterus model made of plasticine, and a few water-filled plastic bags and balloons. The plastic bags and balloons were used to simulate surrounding soft tissue (e.g. muscle, bowel) and also to provide support and stability for the other parts. The marked pig bladder was then placed in the phantom as shown in Figure 3.3. 3.3.2 CT Scan A Foley-catheter was used to fill the bladder with a contrast-water mixture, consisting of 3 mL of ISOVUE contrast agent (Bracco Diagnostic, Inc. Monroe Township, NJ) per 60 mL of sterile 78 Figure 3.1 The study design flowchart showing the main data acquisition and analysis steps. water, in accordance with the gynecological brachytherapy treatment planning protocols in our center. The bladder was sequentially filled with 90 mL, 180 mL, 360 mL and 480 mL of contrast-water and the phantom was CT-scanned for each of the four levels of filling. As an example, Figure 3.4(a) shows three orthogonal cross-sectional planes of the CT image taken from the phantom containing the 180 mL-filled bladder, as displayed by the Eclipse™ treatment planning software (Varian Medical Systems, Inc. Palo Alto, CA). The Hounsfield unit (HU) of the bladder, along with its geometry and deformation due to filling and drooping around the uterus, appear similar to those of real patient images. The characteristic bean-like shape can be seen in the axial view. Also, the Foley balloon was the same as those used in the patients, so its appearance in the images matches real data exactly. The sizes and shapes of the pelvis and the 79 Figure 3.2 The pig bladder with urethra and ureters attached. The yellow rectangular fiducials are rubber bars and the rest are plastic beads. Different groupings of each fiducial type were used. uterus model roughly matched what we expected to see in a CT of real female anatomy. However, the HU value of the foam pelvis bone model was much lower than that of actual bone and the HU of the uterus model was very high because of the high density of plasticine. The bladder surface was contoured in each of the four CT scans using the MIM Maestro software (MIM Software Inc. Cleveland, OH). The CT images were then carefully studied in MIM in order to identify the fiducial markers and track them across different images. Figure 3.4(b) shows an axial slice containing the two fiducial marker types as they appear on CT. Fiducial marker tracking was based on their type, grouping, and rough anatomic position (Anterior-Posterior, Superior-Inferior, Left-Right). The coordinates of the grouped fiducial markers in each scan were averaged to find their centroid. Landmarks that are very close to each other would be expected to move coherently, considering the homogeneity of the bladder material, smooth surface gradients, and uniform external and/or internal forces on the bladder wall. Therefore, a centroid 80 Figure 3.3 Pictures of the pelvic phantom model. The hand-made plasticine uterus model is visible in the bottom image. The Foley-catheter and the CT-safe plastic clamp were used to fill the bladder with the contrast-water mixture and secure the liquid inside. All visible fiducials are circled in red in both pictures. 81 position that could be uniquely identified and localized in all four scans was considered a valid landmark. The bladder neck (point of attachment of urethra to the bladder) was also identified as an anatomical landmark in all data sets. 3.3.3 Point-Set Registration All the registrations were performed using an Intel® Core™ i5-4200U CPU @ 1.60GHz 2.29 GHz machine, with a Windows 8.1 64-bit Operating System. The bladder surface contours were exported into MATLAB (version R2013b, The MathWorks, Inc.: Natwick, MA, USA) as DICOM RT structure files. The DICOM files contain the 3D contour point-sets, which are N×3 matrices where each row contains the (x, y, z) coordinates of a point in the point-set of size N. The number of points in each point-set was different based on the size of the contoured structure, ranging from 19,506 for the smallest bladder volume to 40,568 for the largest. In order to reduce the computation time, the point-sets were down-sized by randomly sampling the points to a certain sample size, based on an empirically determined threshold of 4000, below which the accuracy of the models, and consequently the registration outcomes, were found to become compromised. The valid landmark points for each structure were added to the matrices after the subsampling, to be used for evaluation of the registration accuracy. The elapsed computation time for each method was measured as a performance metric. The registration methods used in this study were available as MATLAB libraries (Table 3.2), with each requiring different numbers and types of input parameters. Some of these parameters are registration parameters, which are the ones that directly affect the registration performance. In order to search for the optimal combination of all registration parameters (i.e. the combination that minimizes the average target registration error, TREavg) for each method, first all registration parameters were set to their default/recommended values. 82 (a) (b) Figure 3.4 (a) The axial (top left), coronal (bottom left), and sagittal (bottom right) planes of the CT scan of the pelvis phantom, with the 180 mL bladder volume filling. The unlabeled structures are the water-filled balloons and bags. (b) The appearance of a plastic bead fiducial and two rubber fiducials on the CT image. While the two fiducial types are completely distinguishable, the rubber bars are more distinct because of their shape and high HU number. 83 Then, using the smallest bladder volume as the reference structure, the value of one parameter was changed, according to the allowed ranges and recommended values, if available, and/or in the direction of minimizing TREavg, with increments established through observing the sensitivity of the outcomes. The parameter value providing the lowest TREavg was then selected as the ‘optimal value’, and the next parameter was adjusted and so on. Table 3.3 shows the registration parameters, their trial values, and optimal values for different algorithms. For the parameters that were assumed to be interacting, a finite set of all combinations of the selected values was tested. The parameter space search was done in increments, because a comprehensive search would require an enormous amount of runtime. However, the increment sizes were tuned based on the observed sensitivity of the outcome so that finer increments were used for larger sensitivities and vice versa. A local minimum was found for all parameters, which was selected as the optimal value. Table 3.2 The software source and available output type/s for each registration method. ‘Correspondence’ indicates the correspondence vector output (C), and ‘deformed points’ indicates the transformed point-set matrix (M) output. Method Source Correspondence Deformed points ICP-finite www.mathworks.com/matlabcentral/fileexchange/24301-finite-iterative-closest-point/ ✓ ✓ TPS-RPM www.cise.ufl.edu/~anand/students/chui/tps-rpm.html ✓ SC www.mathworks.com/matlabcentral/fileexchange/30845-shape-context-basedcorresponding-point-models ✓ CPD https://sites.google.com/site/myronenko/research/cpd ✓ ✓ GMM-reg https://code.google.com/p/gmmreg/ ✓ 84 Table 3.3 The registration parameters, trial and optimal values for the MATLAB implementation of the point-set registration algorithms evaluated in the study. Method Parameter Definition (allowed range, recommended and/or default values) Trial values Optimal value TPS-RPM T_init The initial temperature value (default: 0.5) 0.1, 0.2, 0.3, 1, 5 0.3 T_finalfac The ratio of the T_init to T_final, which is calculated by T_final= T_init/T_finalfac; within the algorithm (default:500) 100, 500, 1000, 5000 500 ICP-finite registration The registration transformation model (options: Rigid, Size (rigid+scaling), Affne – default: Rigid) Rigid, Size, Affine Affine optimizer The optimization method (options: fminlbfgs (finds a local minimum of a function of several variables), fminsearch (MATLAB's multidimensional unconstrained nonlinear minimization function), lsqnonlin (MATLAB's non-linear least-squares-problem solver function), default: fminlbfgs) fminlbfgs, fminsearch, lsqnonlin fminlbfgs TolP Allowed tolerance on distance error (allowed range: (0,1), default: 0.001) 0.99, 0.5, 0.1, 0.01,0.001, 0.001 0.99 TolX Registration position tolerance (default: max(maxP-minP)/1000; where maxP=max(target_points); and minP=min(target_points);) 100, 1, 0.1 (default), 0.01, 0.001 any GMM-reg alpha recommended: 1 0, 0.3, 0.5, 0.8, 1, 3, 5 5 beta The weight of the bending energy in the TPS nonlinear transformation (recommended: 0) 0, 2, 3, 5 0 interval The number of control points in x, y, z direction (recommended: 5) 3, 4, 5, 6, 7, 8, 9 5 scale A scalar parameter representing the covariance matrix of all the Gaussian components in the Gaussian mixtures, with all the components weighted equally (recommended: scale = power(det(model.model/n),1/(2^d));, where n is the number of points in the moving point-set (model) and d is the dimension (3)) 6, 10, 13.37 (from recommended), 15, 20, 30 20 max-iter Maximum number of iterations when finding the local minimum of the cost function (recommended: 100) 100, 200, 300, 400 100 85 Method Parameter Definition (allowed range, recommended and/or default values) Trial values Optimal value CPD Numit Maximum number of iterations allowed towards convergence 300, 400, 500 400 Lambda Regularization weight, i.e. the weight of a prior distribution over the displacements of the points (allowed range: > 0, default: 3) 1, 3, 5, 7 3 Beta Weight of Gaussian kernel, i.e. smoothing filter size that forces rigidity (allowed range: > 0, default: 2) 1, 2, 3, 4 3 Outliers The weight of noise and outliers (allowed range: [0,1], recommended: 0.7) 0.7, 0.8, 0.9 0.7 tol Tolerance stopping criterion, (allowed range: > 0, default: 1e-5) 1e-5, 1e-10 1e-10 fgt Fast Gauss transform (FGT) parameter (allowed values: 0, 1 or 2, Default: 0), 0 - Do not use FGT. 1 - Use FGT to speed up some matrix-vector products. 2 - Use FGT and fine tune (at the end) using truncated kernel approximations. 0, 1, 2 2 eigfgt Whether to use FGT to find the largest eigenvectors of the Gaussian affinity matrix for large data problems (allowed values: 0 or 1, default: 1) 0, 1 1 numeig The number of the largest eigenvectors to use (default: 30, recommended: N1/2, where N is the number of points in T or M) 10, 30, 60, 90 30 SC r_bins Number of log distance bins in the 3D log-spherical matching histogram, as described in Eq. (8.2) in [126] (default:15) 7, 15, 30, 60 30 a_bins Number of angle bins in the matching histogram (default: 15) 1, 5, 10, 15, 30 5 rotate Rotation normalization of the points using eigen vector analysis (0: No correction (default), 1: Minimal align rotation, 2: Rotation with Heaviside flip) 0, 1, 2 1 method The matching method used (0: Munkres (Hungarian [135]) one-to-one matching method, 1: Minimum one to multiple matching method (default)) 0, 1 1 maxdist The maximum matching distance between normalized points. The point-sets are normalized to have an average radius of one to their mean value. (default: 5) 1, 5, 10, 50 5 86 In the course of the CT-scan procedure, a change in the positioning of the bladder inside the phantom and the phantom inside the CT-scanner occurred for the 180 mL fill volume, resulting in a difference in the position and deformation shape of the 180 mL point-set compared to the other volumes. This caused a noticeable shift in the centroid of the 180 mL point-set along the y-axis relative to the other structures, as observed in Table 3.4. No initial alignment was performed prior to the registrations to account for this shift in order to evaluate the performance of the algorithms in the presence of such changes in the initial conditions of the structures to be aligned. Registration requires one of the datasets to be chosen as the reference (target) space that the other three (moving) structures are registered to. The registration algorithms take the (x,y,z) coordinates of the target and moving point-sets (𝑻 and 𝑴, respectively) as input. The output consists of the deformed point-set, 𝑴′, and/or the correspondence vector, 𝑪, which matches the points in 𝑴 to those in 𝑻. Table 3.2 shows the available output types for the implementation used for each method. Table 3.4 The centroid coordinates of the contour point-sets for different bladder volumes. Bladder Volume Centroid (x,y,z) mm 90 mL (4.0, 16.0, −32.9) 180 mL (7.2, 53.9, −39.0) 360 mL (0.6, 1.3, −36.0) 480 mL (3.3, −1.9, −41.7) 3.3.4 MIM DIR Since the phantom used in this study was primarily built for the purpose of point-set registration validation, it is not a suitable phantom for DIR validation, because: (i) to aid fiducial 87 identification, the bladder in the phantom was surrounded by air (rather than water, which would be a closer match to surrounding soft tissue in this body site); (ii) the HU value of the uterus model is very high, compared to real patient image values; and (iii) the presence of the fiducial landmarks (some with high HU values) compromises the credibility of the DIR registration outcomes by providing additional information that is not available in real patient data. Therefore, we modified the phantom CT images so that they better resembled an actual scenario of applying DIR on pelvic images, making a less biased comparison between DIR and point-set registration possible. In order to achieve that, MIM software was used to change the HU values of the image voxels of the fiducials, uterus model and surrounding air. Examples of raw and modified CT images are shown in Figure 3.5. Prior to performing DIR, a landmark point was generated for each fiducial and anatomical landmark in MIM. Then, the deformable intensity-based registration tool was used to register the 180 mL, 360 mL and 480 mL images to the 90 mL image. The landmark points were also transferred in the registered image, and the TRE values for each landmark was calculated by finding the distance between the transferred landmark and the known position of the corresponding landmark in the fixed image. This was done for both raw and modified image sets. 3.3.5 Registration Accuracy Measures 3.3.5.1 TREavg and NE For each landmark, the TRE was determined by calculating the distance between the known coordinates of the landmark in the target structure, and the coordinates of the corresponding landmark either in 𝑴′, if available, or the coordinates of the matched point in 𝑻, using 𝑪. The TRE values were then averaged over all landmarks for each method, to determine TREavg, which was 88 used to represent the spatial accuracy of the registration outcome. Since the choice of the reference structure can potentially affect the registration outcomes, registration accuracies were compared for all possible choices of the reference structure. Because the physical distance between any two corresponding pairs of points was naturally less for smaller reference structures, which would bias their TRE values to be smaller, we defined a normalized error (NE): 𝑁𝐸 = 𝑇𝑅𝐸𝑎𝑣𝑔𝑑𝑖𝑎𝑚𝑒𝑡𝑒𝑟 𝑜𝑓 𝑒𝑞𝑢𝑖𝑣𝑎𝑙𝑒𝑛𝑡 𝑠𝑝ℎ𝑒𝑟𝑒 , (3.1) where the equivalent sphere is a sphere with the same volume as the reference bladder. While bladder was the only ROI structure that was contoured and registered, the registration accuracy was also evaluated for different sub-ROI regions of the bladder based on anatomic sectioning (i.e. posterior sub-ROI, anterior sub-ROI, etc.) in order to compare the registration quality in different anatomic zones. 3.3.5.2 TREtangential and TREperpendicular The quality of registration can also be evaluated in terms of the tangential and expansion/shrinkage components of the registration error. Therefore, we divided the TRE errors for each landmark into components perpendicular and tangential to the bladder surfaces. To achieve this goal, we first found the tangent plane to the surface of the reference bladder at each landmark. Then, we used the frame of reference of the tangent plane to compute the ‘vertical’ (perpendicular) and ‘horizontal’ (tangential) distances to the corresponding point in the deformed point-set (see Figure 3.6). For SC, where no deformed point-set was available, the correspondence vector 𝑪 was used to find the point 𝑷′ in the reference point-set matched to the landmark 𝑷. 89 Figure 3.5 Sample modified CT images shown against the raw images. The modified and raw landmark areas are indicated by red arrows. The HU value of the uterus model region (denoted by green arrows) was also altered to a value closer to that of soft tissue. 3.3.5.3 RMS and HD Additional metrics for registration accuracy were determined using only the difference between point-clouds, without use of the landmarks, including the root-mean-square (RMS) error and the Hausdorff distance (HD). The RMS error was calculated as the root mean square of the distances between the closest-point pairs in the target and deformed point-sets. HD can be defined for two point-sets P and Q as the distance between the point p in point-set P that is farthest from all the points in Q, and its closest neighbour in Q. RMS and HD are commonly used to represent the accuracy in matching the structure shapes as a whole. Since the SC implementation used here only provides the point-correspondence and not the deformed points, it was not possible to calculate RMS and HD for this method. 90 Figure 3.6 A frame of reference is created on the tangent plane to the surface of the bladder at landmark P on the reference bladder. Using the corresponding point in the other deformed point-set, P, one can find the tangential (orange) and perpendicular (pink) distances, i.e. TREtangential and TREperpendicular, by using the normal vector of the tangent plane. 3.4 Results This section presents the outcomes of the different processes, data analyses, and the results of different registration accuracy measures. In section 3.4.1, we present the outcomes and accuracy of the landmark identification process. Robustness and runtime of different methods are discussed in section 3.4.2. The registration accuracy of different methods is evaluated through visual inspection of outcomes and TRE in section 3.4.3, with the smallest volume of bladder chosen as the reference (target) structure. In section 3.4.4, the effect of volume difference on registration accuracy is studied by evaluating TREavg for registrations with different volumes of bladder used as the reference. Section 3.4.5 presents the RMS and HD errors for different reference structures. The effect of random subsampling of the original contour point-sets on the registration accuracy is reported in section 3.4.6. Finally, the registration results of MIM’s DIR method are presented in section 3.4.7 91 3.4.1 Landmark Identification In total, eleven valid fiducial landmarks were successfully identified on all scans (Figure 3.7). The naming convention for the landmarks was based on the material (rub for rubber and pl for plastic fiducials markers) and the anatomic zones (Right/Left, Anterior/Posterior, and Superior/Inferior, or Medial for near the corresponding mid-plane). For example, rub-RPI indicates a rubber fiducial marker in the Right-Posterior-Inferior region of the bladder and pl-LAM indicates a plastic fiducial marker in its Left-Anterior-Medial zone. The landmarks were well distributed over the bladder surface as intended. It was also possible to identify the bladder neck as an anatomical landmark in all scans. In order to evaluate the accuracy and reproducibility of landmark identification in the CT images, two types of fiducial marker localization errors (FLE) were measured. For assessing inter-observer variation, all fiducial landmarks were localized by two different observers, and for intra-observer variation, two times by the same observer. Columns 2 and 3 of Table 3.5 list the average errors for each landmark, over the four scans. The average intra- and inter-observer FLE values over all fiducial landmarks were 0.43 mm and 0.54 mm, respectively, which demonstrates very small discrepancies in locating the landmarks. All the rubber-bar fiducial markers used were successfully identified by both users, since they were clearly visible on all scans due to their shape and high HU value. The plastic beads, however, demonstrated lower visibility and some of the fiducial marker locations were discarded by both users because they were not identifiable and/or traceable in all scans. 92 Figure 3.7 The fiducial landmarks shown on the surface contour point-cloud of the 360 mL bladder. The Anterior, Posterior, Superior, and Inferior directions are shown as A, P, S, and I, respectively. The landmark points facing the front of the view (anatomic Right) are shown in dark green and the ones on the back surface (anatomic Left), in light green. Note that pl-RPS is behind a fold on the surface that is part of a dent in the Superior–Posterior zone from the uterus indentation. Each user was also tasked to track an anatomical landmark in the bladder neck, to (i) compare its accuracy against the fiducial landmarks and (ii) explore the feasibility of using this landmark in human subjects where no fiducial markers are present. The bladder neck landmark demonstrated the highest inter-observer FLE errors, which can be explained by high subjectivity in its localization due to the orientation of the urethra in the phantom setup and the fact that the margin between the bladder and the opening of urethra was not well-defined. Therefore, this landmark was excluded from all TREavg values reported. 93 Table 3.5 The TRE values and the inter-observer and intra-observer fiducial localization errors (FLE) for all landmarks, in millimeters. Landmark Intra FLE Inter FLE ICP-finite TPS-RPM SC CPD GMM-reg Deformable methods rub-LAI 0.44 0.58 15.6 6.7 6.8 7.5 8.6 7.4±0.9 rub-RAI 0.27 0.33 20.1 14.7† 9.7† 15.3† 18.2† 14.5±3.5 rub-LAS 0.44 0.60 22.2 10.4†† 9.0 7.6†† 10.2†† 9.3±1.3 rub-RPI 0.44 0.53 9.0 6.7 7.7 6.7 8.7 7.5±0.9 rub-RPS 0.50 0.65 11.7 3.9 5.6 3.0 3.6 4.0±1.1 rub-LPS 0.57 0.47 9.0 5.8 6.6 5.1 6.9 6.1±0.8 pl-LAM 0.29 0.45 14.5 5.7 8.3 4.2 6.6 6.2±1.7 pl-RAS 0.33 0.28 16.2 6.8 6.9 7.6†† 8.0 7.3±0.6 pl-LPS 0.42 0.67 12.7 3.6 5.8 4.1 9.1 5.7±2.5 pl-RPS 0.59 0.95 8.1 3.1 9.5†† 2.0 2.4 4.3±3.5 pl-RMM 0.45 0.43 12.5 6.0 3.9 7.0 6.2 5.8±1.3 Blad-Neck 0.29 1.6 21.8 7.9 12.4 5.4 8.3 8.5±2.9 Fiducial landmarks 0.43±0.10 0.54±0.18 13.8±4.5 6.7±3.3 7.3±1.8 6.4±3.5 8.0±4.1 7.1±2.9 Note: Columns 4–8 list the TRE errors for the five methods assessed, averaged over all three registrations using 90 mL structure as the reference, and the last column lists average TRE values for each landmark across all methods except ICP-finite. The last row shows the average FLE and TRE values (i.e. TREavg) across all fiducial landmarks for each method. The lowest TRE values for each landmark and lowest TREavg is shown in bold. The largest and the second largest fiducial landmark TRE’s for each deformable method are marked by † and †† symbols, respectively. 3.4.2 Robustness and Runtime CPD, GMM-reg and SC ran successfully on all the registrations performed on the data, over a wide range of parameter settings. However, TPS-RPM and ICP-finite both demonstrated convergence issues. In some cases, TPS-RPM registration collapsed with small-increment changes of one of the registration parameters (T_init). Overall, the TPS-RPM software, was found to be highly sensitive to the initial parameter settings, to the point of total failure with slight changes in the parameters. ICP-finite, also, showed abnormal outcomes in registering some of the cases. These anomalies are discussed in more detail in section 3.4.3. Although the optimum sample size (i.e. 4 K) was used for the evaluation of the registration 94 accuracies, the algorithms’ runtimes were measured for different dataset sizes. Figure 3.8 shows the runtime for a single registration process for all methods versus the sample size (in log-scale in order to optimally illustrate the wide range of runtime values). We did not test TPS-RPM above 4K sample size, since it became infeasible to run. While ICP-finite was the fastest method, due likely to its affine transformation model, CPD was the fastest deformable method. The runtime increased with the number of points in the datasets for all methods, as anticipated. 3.4.3 Registration Accuracy: TRE Optimal settings for the registration parameters of each algorithm, as shown in the last column of Table 3.3, were used with the smallest bladder structure set as the reference. Visual inspection of all registration outputs, based on the overlapping of the target and deformed point-clouds (Figure 3.9), demonstrated an overall good match between the shapes with minor discrepancies for all deformable methods. On the other hand, ICP-finite showed a high degree of mismatch in registering 360 mL (Figure 3.9 (e)) and 480 mL structures to the 90 mL structure, and failed in matching the 180 mL volume to 90 mL (Figure 3.10). Among the three deformable methods that offered the transformed points as outputs (i.e. TPS-RPM, ICP-finite, and GMM-reg, as shown in Table 3.2), TPS-RPM and GMM-reg exhibited the least and most amount of discrepancy, respectively. The matching between target and moving point-sets by SC is illustrated by the green vectors connecting the corresponding points (for one-tenth of the points only for visibility). It can be seen that the matching vector field is rather well-regularized, so that the vectors disperse in a uniform radial pattern, with no large cross-matching observed except on the left side, where the point-sets overlap. ICP-finite performs poorly compared to the deformable methods, as expected. The last row in Table 3.5 shows the TREavg values achieved for each method. ICP-finite 95 demonstrates substantially higher TREavg values compared to the deformable methods, which is expected due to the fact that it is limited to affine transformations and cannot account for the non-linear deformations in the registered structures. While the TREavg values for the deformable methods were within one standard deviation of one another, CPD gave the smallest TREavg and also the lowest TRE values for the majority of the landmarks (7 out of 12). Figure 3.11 shows the registration results for the tested algorithms, for all the fiducial landmarks in three orthogonal views. Figure 3.8 The change in the registration runtime with sample size. The plot is in log scale for better representation of the variations. TPS-RPM was not tested with sample sizes larger than 4 K due to its very long runtime. 0.11101001000100000 2000 4000 6000 8000 10000 12000 14000Runtime (s)Sample size (#points)CPD GMM SCICP TPS-RPM96 (a) (b) (c) (d) (e) (f) Figure 3.9 (a) The point-clouds of the target (90 mL, in this example) point-set, in blue, and moving (360 mL), in red, before registration. ((b)–(e)) The target and deformed point-clouds after registration with (b) CPD, (c) GMM-reg, (d) TPS-RPM, and (e) ICP-finite. Note that the scale of the plot in (a) is different from the rest of the plots, in order to better see the details in the shape of the point-clouds. (f) The target and moving point-clouds matched by SC, showing the correspondences for 1/tenth of the points, as green vectors. 97 Figure 3.10 Registration failure of ICP-finite method, showing the deformed 180 mL point-set (blue) and the 90 mL reference structure (red). Columns 4–8 of Table 3.5 show the TRE values for individual landmarks in each row for all methods. The TRE in the ‘Bladder Neck’ anatomical landmark was comparable to the fiducial TREavg values for all methods but SC, where Blad-Neck had one of the highest landmark errors. Since for most landmarks, the TRE values for ICP-finite varied considerably from the deformable methods, it was excluded from the average landmark errors in the last column. Rub-RAI had the highest average TRE and was also the most erroneous landmark in each deformable method. The second highest error belonged to rub-LAS for three out of four deformable methods. On the other hand, pl-RPS and rub-RPS had the two lowest average TREs and were collectively the least erroneous landmarks for each method. Table 3.6 shows the TRE values averaged by left/right, anterior/posterior and superior/inferior halves to better show the spatial dependence of registration accuracy. Considering the significant difference in the TREs of the anterior and posterior landmarks (𝑝 < 0.05), and the observations in the previous paragraph, it was concluded that the anterior part of the bladder is much harder to re- 98 Figure 3.11 Comparison of the registered landmark locations across different registration methods. The reference landmarks of the 90 mL volume are depicted by red stars and are numbered in the axial, coronal, and sagittal views of the reference point cloud. The corresponding average location, across the three moving structures, computed using five registration methods are shown in different colours. The projection of the TREavg vectors for each method are also displayed in each view to facilitate comparison. gister than the posterior part. This could be due to (i) free expansion of the anterior section compared to the posterior section, and/or (ii) the posterior indentation of bladder surface where it drapes over the uterus, which might serve as a “guiding” feature that helps improve the registration 99 quality in that area. The inferior landmarks showed higher TREs than the superior ones; however, the difference was not quite significant (𝑝 ≅ 0.05). Finally, the TREs for the landmarks in the right and left sides were comparable, which is in agreement with the inherent left/right symmetry of the bladder. The tangential and perpendicular components of the TRE vectors for each landmark were also calculated for each method, as described in section 3.3.5.2. These TRE components were then averaged over the three registrations and eleven landmarks to calculate a mean TREtangential and mean TREperpendicular for each method. The results for mean TREtangential and TREperpendicular with respect to the TREavg for all methods are shown in Figure 3.12. Despite ICP-finite producing the largest TREavg magnitude (the blue circle with the largest radius in Figure 3.12), the relative contributions to this TRE from the tangential and perpendicular components were similar. On the other hand, for the four other methods (CPD, GMM-reg, SC, and TPS-RPM), the tangential components had clearly larger contributions to TREavg than the perpendicular ones. Thus, the deformable methods had a greater tendency towards tangential deformation misalignment, with TPS-RPM marginally showing the highest and GMM-reg, the least. 3.4.4 Registration Accuracy: Effect of Volume Difference The effect of volume difference for the 90 mL reference fraction is shown in Figure 3.13(a), which plots TREavg for different moving structures. It can be observed that the error for the 360 mL structure was lower than that of 480 mL for all the methods and also, for SC and GMM-reg, TREavg increased with the volume of the moving structure. However, for TPS-RPM, ICP-finite, and CPD, the registration errors for the 180 mL moving structure were higher than those of 360 mL. The 100 large spike observed in the 180 mL value for ICP-finite was due to the failure of the algorithm to properly register this case (as illustrated in Figure 3.10), which persisted for all registration parameter settings. All the results presented so far belonged to the case where the smallest bladder was used as the reference. Here, and in the next section, registration error results for the cases where different volumes of bladder were selected as the reference structure are presented, in order to also study the effect of reference volume choice. Figure 3.13(b) shows a bar graph of NE, in percentage, for different reference structures. TPS-RPM, SC, and CPD gave very similar mean errors, lower than that of GMM-reg. All methods showed a very similar pattern, with the 360 mL reference volume consistently giving the lowest NE value for all methods and the 180 mL volume the highest, except for SC. The recurring failed registration outcomes by ICP-finite, similar to Figure 3.10, were partly responsible for the high NE value for the 180 mL volume. Table 3.6 Average TRE values for the fiducial landmarks categorized based on their anatomic positions (smallest bladder used as target). Anatomic zone TRE (mm) Number Posterior 5.5±1.4 5 Anterior 8.9±3.3 5 Superior 6.1±2.0 6 Inferior 9.8±4.1 3 Right 7.2±3.8 6 Left 6.9±1.5 5 101 Figure 3.12 The vector diagram representing TREavg (arrows), TREtangential (horizontal dashed lines), and TREperpendicular (vertical dashed lines) for all five point-set registration programs. The circular traces are shown for aiding the comparison of the TREavg among different methods. 3.4.5 Registration Accuracy: RMS and HD The RMS error and HD values for all registration methods but SC are shown in Figure 3.14. Using the 90 mL reference structure, RMS and HD values were compared to values achieved using the other structures as reference. As expected, the values for ICP-finite were much larger than the deformable methods. The RMS values for the deformable models were quite small, demonstrating high quality shape-matching by all these methods. TPS-RPM gave the smallest mean RMS and 102 (a) (b) Figure 3.13 (a) The bar graph of TREavg values for different registration methods, separated by the moving structure volumes (with the 90 mL structure used as the reference). The error bars in both plots represent the SD over different landmarks and do not necessarily represent the uncertainty in the TRE values. (b) The bar graph of normalized error (NE) values for different registration methods, separated by the reference structure volumes. The mean values over all reference volumes are shown above the bars for each method. 103 HD values, as well as the smallest individual RMS and HD, except for one case. This was in agreement with observations based on visual inspection of the registration outputs in Figure 3.9. 3.4.6 Subsampling Effect In order to study the effect of the random subsampling of the original contour point-sets on the registration accuracy, an experiment using all methods except TPS-RPM, because of its very long runtimes (as discussed in section 3.4.2), was performed. The four point-sets were randomly sampled 10 times, to 2 K, 4 K, 8 K and 12 K points. The sub-sampled point-sets were then registered using each method (with the smallest bladder filling volume as target) and the mean TREavg values over the ten sets of samples were calculated together with the standard deviation (SD), as shown in Figure 3.15. The ICP-finite results showed a small but consistent increase in TREavg values with increasing sample size (13.88±0.06, 13.89±0.07, 13.90±0.02, and 13.94±0.03 mm, in the order of increasing sample size); however, the SD did not show any particular pattern. Since ICP-finite results were out of range of the other methods it was not included in the plot. The variations in TREavg values across different sample sizes were also relatively small for the deformable methods. The 4 K sample size showed the lowest TREavg values for CPD and GMM, while for SC, the 12 K dataset achieved the lowest error. The 2 K sample size was consistently the most erroneous of all, suggesting that the low density of the points at this sample size has a direct effect on the quality and accuracy of the registration. For samples larger than 2 K, the variation in TREavg did not follow an identifiable pattern and it is unclear what sample size is representing the most accurate registration results. 104 (a) (b) Figure 3.14 The bar graphs of the RMS error and HD for different registration methods, except SC, separated by the reference structure. Mean and SD values are shown above the bars for each method. 105 For all three deformable methods, the standard deviation of TREavg almost consistently decreased with increasing sample size (except the 8 K for SC). This is expected because the larger the sample size, the closer the sample is to the actual data, and there is less variation in the sampling of the points. Therefore, using a larger sample size would provide more regularity in the results. Relatively low TRE values and small variation in TREavg and SD values across all methods, combined with reasonably low runtimes, all factored into the decision to use a sample size of 4 K for comparison between methods. 3.4.7 MIM DIR Outcomes The average TRE values across all landmarks (TREavg) for registering the images corresponding to the 180 mL, 360 mL and 480 mL volumes to the 90 mL bladder volume image were 8.2±5.6 mm, 6.2±3.0 mm, and 8.7±3.2 mm, respectively, for the modified CT images, with an average of 7.7±1.3 mm over all three registrations. The corresponding values for the raw images were 7.2±4.3 mm, 6.0±3.3 mm, and 8.1±3.8 mm, respectively, with an average of 7.1±1.1 mm. Figure 3.16 shows examples of the MIM DIR outcome for registering the 480 mL raw CT image to that of 90 mL in three different positions on the bladder surface. Good alignment was achieved for the position in the images in the bottom row and acceptable alignment in the middle row; however, the alignments in the top row images were far from accurate. 3.5 Discussion In this study, eleven fiducial landmarks were identified on CT images of the phantom, and these landmark positions were used to estimate the TRE values for evaluation of five registration methods. The TREavg values for the four deformable techniques were in the range of 6.4–8.1 mm, 106 when the 90 mL bladder structure was used as the reference dataset. These TREavg values are, to a good approximation, representative of the registration accuracies across the bladder surface because: (i) if the surface area of the largest bladder size is estimated by the surface area of a sphere with radius ≃4.86 cm (A = 295.5 cm2), there is one landmark location per ∼27 cm2 area on the surface of the bladder, which would be equivalent to a square of side ∼5 cm (∼3 cm for the smallest bladder), (ii) the landmark locations were well distributed on the bladder surface, and (iii) due to the type of the existing external forces on the bladder as an elastic object, abrupt changes in its shape and surface are not expected over small distances on the order of centimeters. Therefore, points between the landmarks were not expected to have registration accuracies drastically outside the range of the observed landmark TRE values. Figure 3.15 The mean and SD (error bars) of TREavg for different sample sizes. TPS-RPM was not included this study due to its very slow processing. The ICP-finite data was excluded from this graph because it showed almost no variation in both mean and SD values and also it was out of the range of other data. 107 Figure 3.16 Examples of the MIM DIR outcome, for registering the 480 mL raw CT image to that of 90 mL, showing the axial (left), sagittal (middle) and coronal (right) views of the deformed image. MIM’s ‘RegReveal’ tool is used, in three different positions on the bladder surface, to reveal the achieved alignment of the moving and target images inside the smaller inner box. Nevertheless, we could not guarantee the TRE bounds to apply in all areas of the bladder (e.g. those most distant from any of the landmarks), but TRE remains the most commonly adopted and accepted measure for evaluating registration. One option for future work is to explore the possibility of having a much larger number of fiducial markers while ensuring minimal effect on 108 tissue deformation. That said, providing a precise and global (over the whole field of view) validation of registration remains an open challenge. In a study similar to this one, Wognum et al [110] evaluated one structure-based registration algorithm (STPSRPM) and seven different intensity-based DIR algorithms. They used 30–40 radiopaque fiducial markers on the exterior surface of a porcine bladder in a water phantom. CT scans for a range of different bladder filling volumes between 100 cc and 400 cc were used for registering with the DIR methods. The contours of the bladder surfaces were used for STPSRPM. Similar to our TRE calculations, they used the Euclidean distances between the known marker positions on the smallest volume as the reference and the transformed positions on the larger bladders, to calculate a residual distance error (RDE) to quantify the registration spatial accuracies. They also reported other error metrics, including HD, surface distance error (SDE, which is similar to RMS, except that instead of the root mean square it is the mean of the distances), and Dice similarity coefficient (DSC, usually used for DIR evaluation). The runtimes for different methods were not reported. In Wognum’s study, STPSRPM showed superior SDE and HD values (<0.5 mm, and <5 mm, respectively) compared to the DIR methods. The RDE values consistently increased with the volume difference between the reference (100 mL) and moving structures, for all methods. The best RDE results (∼3–8 mm for different moving structures) were obtained with STPSRPM. Although the highest accuracies achieved by STPSRPM in Wognum’s paper are better than those achieved in our study, our results are better than most of the DIR methods evaluated. Moreover, it must be noted that the shape of the bladder in the Wognum’s setup is simpler than in this study, and they also used a rigid pre-alignment using three reference markers, which can potentially aid 109 and improve the deformable registration outcomes. The small marker errors achieved by most of the DIR methods in Kirby et al [112] might poorly represent how these methods perform for actual patient data, because of larger deformations (e.g. varying bladder sizes up to three-fold in some practices), change in position of organs, and existence of surrounding structures, such as applicators. In Varadhan et al [113], the best achieved Hausdorff distance error for the bladder contour structure (7.7 mm, by B-spline algorithm) is higher than the average values achieved by two of the point-set registration methods evaluated in our study (6.8 mm for TPS-RPM and 7.3 mm for GMM-reg). The lowest root mean square error achieved for the prostate (with much less degrees of deformation compared to the bladder) is also more than seven-fold larger than the lowest RMS values in our study. In our study, CPD, TPS-RPM, and GMM-reg demonstrated good registration outcomes based on visual inspection of the target and deformed point-clouds. In addition, the RMS error values were quite low for these methods, suggesting a good overall match of the shapes. While HD values achieved by deformable methods were larger than the result achieved by STPSRPM in Wognum’s study, TPS-RPM showed reasonable results with HD errors in the range of 6.4–7.3 mm. CPD and GMM-reg demonstrated larger errors, suggesting that we have errors of up to 1 cm and more in matching some parts of the bladder surfaces. The ICP-finite error values were significantly larger than the deformable methods, especially for HD. It must be noted that although RMS and HD are commonly used to represent the accuracy in matching the structure shapes as a whole, small RMS and HD values do not necessarily mean an accurate registration, due to susceptibility to error in capturing the shear and rotational motions. It 110 was observed that RMS falls short in correctly representing the shape-matching accuracy, in the cases where ICP-finite failed to register properly, such as matching the 180 mL to the 90 mL structure (Figure 3.10). While the TRE and HD values are high for these registrations (Figure 3.13(a) and Figure 3.14(a)), as expected, RMS renders deceptively low values (Figure 3.14(b)). Therefore, it was concluded that RMS is only useful when the registration is successful and the shapes are at least roughly matched. The mean TREavg value achieved by MIM DIR for registering the modified CT images (where the larger bladder volumes were all registered to the smallest), was higher than three out of four deformable PSR methods. Also, the individual TREavg values for 180 mL, 360 mL and 480 mL volumes registered to the 90 mL image were 14%, 15%, and 32% larger than the best corresponding results achieved by CPD. Even the raw images yielded up to 23% larger TREavg values compared to the CPD results. Moreover, visual inspection of the registration results suggested unsatisfactory registration in some of the regions of the bladder surface. Therefore, the results of this study confirm our initial findings of poor performance of various DIR methods for this complicated task, particularly when it is considered that registering actual patient images is anticipated to be a much more complicated task for DIR due to potentially significant variations in the position, size, orientation, and/or shape of high contrast regions such as the bladder filling, brachytherapy applicators, Foley balloon and air bubbles, in the moving and target images. Since the surface tangential motion is generally hard to estimate using DIR and shape matching (given the typically homogenous appearance along the surface), the performance of the methods in handling tangential motion was also evaluated. An interesting observation was made from these results, shown in Figure 3.12. Despite its overall worse performance in TRE, the ICP-111 finite method did not show a stronger tendency to produce tangential error compared to inward-outward (perpendicular) error; contrarily, a clearly stronger tangential component was observed with the other techniques. The reason for this difference in performance may be attributed to the fact that ICP-finite is a linear transformation encoding a (linear) tangential component (along with rotation and scaling), whereas the deformable methods deform the point clouds non-linearly and do not have a similar pure tangential component. In current radiotherapy practice, dose parameters for the bladder are usually either calculated for the whole volume (whether including the filling or only the wall) or the hottest volumes of a certain absolute or percentage volumes [44, 143, 144, 145, 146]. The smallest high-dose volume that is usually considered clinically relevant is 2 mL, as reflected in the commonly reported D2cc dose-volume parameter (where cc is conventionally used to denote cm3, or mL), which is equivalent to a circular area of 11.3 mm radius on the bladder surface, with a typical wall thickness of 5 mm. The TREavg accuracies of 6.4–13.8 mm obtained by the non-rigid methods in this study were mostly smaller than the dimensions above. The maximum landmark TRE errors across all four deformable methods ranged from 9.7 to 18.2 mm. However, considering that rub-RAI consistently showed the highest TRE values for the deformable methods and appeared to be an outlier, the maximum landmark TRE errors excluding rub-RAI values ranged between 7.6–10.4 mm, which is smaller than the size of the high-dose area of interest. This suggests that the use of deformable registration would allow cumulative dose volume estimates to be made for the bladder wall with significantly better accuracy than the current practice, where parameters such as D2cc are routinely summed across fractions without any consideration of their location and whether or not the high 112 dose regions overlap across the fractions. Furthermore, by accumulating the fraction dose using the optimum registration results found in this study, we are guaranteed to identify localized cumulative high-dose regions on the bladder wall. It should also be possible to determine when high-dose regions are not in close enough proximity to overlap; thereby correcting the assumption of complete overlap of the high-dose areas that underlies the current practice of summing DVH-based dose-volume parameters over all fractions. Based on the observations from Figure 3.13(a) for SC and GMM-reg, and also the consistent trend between 360 mL and 480 mL volumes, we may conclude that an increase in the registration error occurs with increase in the volume difference between the moving and reference structures. Therefore, some of the observations for TPS-RPM, CPD and ICP-finite might have been attributable to the experimental setup for the 180 mL bladder, which required significant translation of the structure during registration, and may have caused distortion of the bladder shape in the process. While all registration techniques tested here do allow for translation, they may not be invariant to the initial conditions. Therefore, the observed high TRE values for the 180 mL structure for these three methods could be due to its different initial condition. Based on this theory, we can conclude that the registration outcomes of SC and GMM are less prone to become compromised by initial conditions of the aligning structures. More data would be required to test this theory, and to evaluate alternative means of correcting for translation prior to deformable registration of the point clouds. Trying different structures as the reference dataset showed that the 360 mL-volume provided 113 the lowest NE values for all deformable models. While the exact reason for this outcome is unclear, we could speculate that choice of a reference structure with a volume close to the mean of all structures minimizes the normalized error. While this is contradicted by high errors for the 180 mL reference volume, compared to 90 mL and/or 480 mL for most cases; this may have been due to its initial condition, as discussed previously. Looking at the distribution of the landmark TRE values based on their anatomic positions, it was found that the registration methods perform better in matching the posterior part of the bladder (TREavg = 5.6 mm) compared to the anterior, with statistical significance. This could be due to anatomical factors or by shape features, such as the indentation of the posterior bladder due, presumably, to pressure from posterior structures such as the uterus, that might assist the registration. However, in real patients, the uterus and other organs posterior to the bladder are mobile/deformable and may not cause similar deformations seen here. Therefore, the significant difference we observed in the TRE values of the posterior and anterior landmarks might be less pronounced in actual patient cases. The TRE values obtained for the bladder neck as an anatomical landmark were mostly near the TREavg for each method, except for SC. This suggests that while localization of this anatomical landmark was subject to some inaccuracy – manifested in its higher FLE value – it still provided a good measure of the overall performance of each method and could be used as a stand-alone measure for evaluating the registration accuracy. In the phantom images, it was not possible to identify the opening of the ureters to the bladder and use them as additional anatomical landmarks. However, in actual patient data, given good image quality and a sufficiently experienced observer, it is possible to identify these landmarks. Therefore, it could be feasible to use TRE values of these 114 three anatomical landmarks with the patient data as a measure of the registration accuracy. The computation time varied considerably for different software, with the ICP-finite being the fastest of all. CPD showed the shortest runtime among deformable registration software options, with a very large margin compared to the method with the longest runtime, i.e. TPS-RPM. The runtimes increased with the sample size for all the methods, as expected. While the sample size and random sub-sampling did not have a considerable effect on the TREavg values, the regularity increased with increasing sample size. However, for most applications, the degree of this improvement did not outweigh the increased computation time for very large sample sizes. Therefore, the choice of 4 K or 8 K sample size for representing the bladder contour data seems appropriate. 3.6 Conclusions In this work, five different point-set registration (PSR) methods and the DIR tool of a commercially available software (MIM Maestro) were evaluated for registering bladder surface contours of different sizes and shapes. Four of the PSR methods, CPD, GMM-reg, TPS-RPM, and SC are based on deformable registration algorithms. An affine version of the commonly used method ICP, ICP-finite, was also included in the study. Three different common error metrics were used to assess different aspects of the registration outcomes, i.e. registration spatial accuracy by TRE and shape-matching accuracy by RMS and HD. It was shown that a non-deformable registration algorithm is clearly not an appropriate choice for registering bladder surface point-sets. All deformable methods provided relatively good registration outcomes; therefore, we suggest that the choice of the most suitable deformable 115 method would depend on the application and the user’s criteria. CPD seems to be the most well-rounded method, considering its fast processing and low TRE values for the smallest reference structure. TPS-RPM provided better overall structure matching among all methods, but at the cost of substantially longer runtime and diminished robustness. The MIM DIR tool provided inferior results compared to most of the deformable PSR methods, both from TRE values and visual inspection of the outcomes. The ICP-finite affine transformation was found to be less prone to dominating perpendicular transformation and showed more balanced tangential error components compared to the deformable methods, but otherwise had poor performance overall. Therefore, addressing ways to reduce tangential error remains an area of future research for general shape/image registration. While the phantom used in this study does not provide a fair comparison of DIR performance compared to the deformable point-set registration methods, the DIR registration errors were inferior to the best TREavg results achieved by deformable point-set registration in this study. Further investigation using phantoms that model true tissue image densities more closely would be required to confirm this finding. The best TREavg registration accuracy found in this study was 6.4 mm, which is reasonable considering the size and elasticity of the bladder. This achievement improves on the current practice of summing volumetric dose parameters for bladder at each fraction in brachytherapy treatment planning, where high-dose regions can be well separated between some or all fractions. We therefore conclude that, while our current implementations of these point-matching methods are not yet fully developed for routine clinical use, incorporation of a suitable deformable point-116 set registration method into existing clinical workflows has the potential to improve the accuracy of cumulative dose estimates, for bladder at a minimum, in fractionated radiotherapy. While testing point-set registration on other surrounding tissues and organs was outside the scope of this study, an interested reader might test these methods on other structures through building phantoms similar to ours. It must also be recognized, however, that not all structures of interest in the irradiated volume may be suited to these methods, and that a combination of methods, including DIR, might be required for a complete analysis. 117 Chapter 4: Bladder Accumulated Dose and Urinary Toxicity 4.1 Introduction Although the use of 3D imaging in gynaecological brachytherapy treatment planning was first reported over 25 years ago [147] it has only become routine clinical practice in the past ten to fifteen years [148, 149]. During the same period, fractionated high dose rate (HDR) has slowly replaced low dose rate (LDR) as the brachytherapy component of standard regimens for treatment of locally advanced carcinoma of the cervix. HDR lends itself well to CT or MRI-based planning, as the brachytherapy applicators for this type of delivery can be constructed from materials compatible with these imaging modalities. Thus, at the BC Cancer Agency and many other centres, definitive treatment of locally advanced cervical carcinoma now consists of fractionated EBRT in combination with multi-fraction HDR image-guided intracavitary brachytherapy (IG-ICBT). With the advantages of image-based treatment planning, it is possible to calculate the dose at every voxel of the image and to delineate soft tissue structures. Therefore, within the limits of uncertainty in structure delineation, potential movement of the brachytherapy applicators between imaging and treatment, and the accuracy of the dose calculation algorithm, dose delivered to OAR such as the bladder, rectum, uterus, vagina, and paracervical tissues can be determined with reasonable precision for each treatment fraction. These capabilities have been recognized in the standard recommendations for image-based dose recording, and a number of DVH-based parameters for both target and OAR volumes have been recommended for dose reporting [150, 151, 42, 43]. 118 4.1.1 Limitations of DVH Parameters The ability to quantify the accumulated dose received by defined regions of tumour or normal structure volumes over multiple radiotherapy fractions remains limited due to significant variations in the size and relative positions of most, if not all, of these structures from one fraction to the next. These variations are caused by factors such as variable bladder, bowel, and rectal filling, tumour regression, and displacement of structures by the intracavitary brachytherapy apparatus. Many of these factors are not easily controlled, with the exception of bladder filling. Thus, the GYN GEC-ESTRO task group [43] has advised that consistent bladder filling volumes be used to aid in valid and reliable data collection. Nevertheless, even with consistent bladder filling, the DVH-summed dosimetric parameters may still not be accurate surrogates for cumulative dose because of possible organ deformations caused by the other factors mentioned above. An additional limitation of current methods of dose accumulation is that DVH parameters are usually calculated for the whole volume of the hollow OAR, including the contents, whereas only the wall tissue is of radiobiological interest. For these reasons, in order to more accurately assess tolerance doses and then apply them in practice, it is clear that improved methods are needed to replace summed DVH parameters for estimating accumulated dose within the hollow organ walls in situations of repeated fractions with changes in organ shape or size. 4.1.2 Practice of Varying Bladder Filling Status at BCCA It has been shown that bladder distension can significantly reduce the dose received by the bowel, with no significant difference in the dose distribution delivered to the target volume in cervical HDR ICBT [152, 153, 154, 155]. Bladder distension, however, might slightly increase the 119 bladder dose itself [97] and present a tradeoff between benefits for different OAR’s. Therefore, varying the amount of bladder filling in different fractions might provide a balanced sparing of the OAR’s. Moreover, it can be hypothesized that varying the bladder filling status might effectively change the location of dosimetric hotspots on the bladder in consecutive HDR-ICBT treatment fractions, thereby limiting the magnitude of cumulative dose maxima over the course of treatment. Therefore, in an effort to reduce the overall dose maxima to the bladder, the practice of varying the amount of bladder filling in consecutive fractions of high-dose-rate (HDR) brachytherapy for cervical cancer patients has been implemented at the British Columbia Cancer Agency (BCCA) since 2010. In doing so, it was recognized that the resulting bladder size variation is expected to cause even larger discrepancies between the actual cumulative dose parameters and the plan-summed DVH parameters used in practice [156, 157, 43]. This has motivated a more rigorous assessment of the accumulated bladder dose due to this practice, which is presented in this chapter. The effectiveness of varying bladder filling volume in reducing urinary toxicity is also investigated here. 4.1.3 Genitourinary Toxicity and Radiation Dose One of the known late effects of radiotherapy for locally advanced carcinoma of the cervix is genitourinary (GU) sequelae. In general, urinary morbidity in pelvic radiation therapy can be attributed to damage to the lower urinary tract, including bladder, urethra, and urethral opening (i.e. bladder neck) [158]. Ureteral injury and injury to the ureteral anti-reflux sphincter mechanism at the ureteral opening to bladder can also occur, although this is rare following radiotherapy. Radiation induced urinary sequelae include different symptoms, such as frequency (frequent 120 urination), urgency (sudden urge to urinate), dysuria (painful urination), nocturia (need to wake up during night to void), hematuria (bloody urine), and incontinence (involuntary voiding). In the LENT-SOMA subjective scoring system, these symptoms are given individual subscores (0, 1, 2, 3 or 4). There are slightly different approaches for combining urinary subscores into an overall “toxicity score/grade” and some studies do not even specify the approach taken [159, 160]. While some groups use the average of the subscores as a final toxicity score [161] as recommended by the working group that produced the LENT-SOMA system, others take the highest subscore as the “toxicity grade” (i.e. G0-G4) [46, 162]. The subtleties and controversies on this subject will be discussed in Chapter 5. GU symptoms might be due to “global” injuries, such as radiation-induced cystitis (bladder infection and inflammation), or “focal” injuries in more severe conditions, such as fistula (abnormal opening in bladder wall), necrosis, ulceration, etc. Some GU sequelae might in fact be the result of radiation damage to the urethra and bladder neck [158, 146]. It has been reported that dose to the urethra and bladder neck is an important predictor for late urinary toxicity in prostate brachytherapy [159] and interstitial gynecological brachytherapy [160]. However, only a very small body of literature exists on this topic and investigations using deformable dose accumulation in fractionated treatments are not available. Several studies have examined the relationship between dose and urinary toxicity in EBRT treatment or brachytherapy for prostate cancer [163, 164, 165, 145, 166, 167, 168], but there have been only a few dosimetric studies of urinary toxicity for BT in cervical cancer patients [46, 169, 170, 171, 162]. The most comprehensive of these [46, 169], analyzed DVH-summed D2cc, D1cc, 121 and D0.1cc parameters (in terms of equivalent dose in 2 Gy fractions) as predictors of late urinary toxicity in a cohort of 141 patients treated with EBRT and IG-ICBT. They found a significant dose effect for all summed DVH parameters when comparing groups with minor toxicity (Grade 0-1) to those with major toxicity (Grade 2-4). Mazeron et al. [162] found a statistically significant association between urinary morbidity and DVH-summed D2cc and D0.1cc in pulsed-dose rate brachytherapy. However, the sample sizes analyzed in these studies were small, and they did not examine the relationship between the true locally accumulated bladder dose and bladder complications. 4.1.4 Chapter Overview The main purpose of the study described in this chapter was to develop a method to more accurately accumulate the fractionated dose to the bladder wall in HDR IG-ICBT for carcinoma of the cervix. It is hypothesized that the locally-accumulated dose parameters are significantly different from the equivalent DVH-summed parameters, particularly when the amount of bladder filling is varied between fractions. In addition, the locally-accumulated dose parameters were retrospectively studied as predictors of the observed urinary toxicity. The effect of alternating the bladder filling status in consecutive HDR ICBT fractions on urinary toxicity was also investigated. Furthermore, the localized dose parameters around the bladder neck and ureteral openings and their association with urinary toxicity score were assessed. For dose accumulation, the correspondence between the bladder structures in different fractions was estimated using a deformable point-set registration (PSR) method, called coherent point drift (CPD). A software implementation of the CPD algorithm, proposed by Myronenko et 122 al. (2010), xis available as a toolbo for MATLAB (The MathWorks Inc., Natick, MA, version R2013b), with adjustable registration parameters. In Chapter 3, the performance of CPD was thoroughly evaluated for deformable bladder registration against four alternative PSR algorithms as well as a commercial DIR method (MIM). CPD showed the lowest average target registration error (TRE) of 6.4 mm amongst all methods. It was therefore selected as the method of choice for bladder dose accumulation in this study. 4.1.5 Previous Work Similar attempts to accumulate radiotherapy fractionated dose in highly deformable organs such as bladder and rectum have been reported in the literature [172, 173, 174, 156, 175, 176, 46, 177, 157, 178]. Three of the most relevant studies were by Andersen et al. (2013), Jamema et al. (2015), and Teo et al. (2015). Andersen et al. used a biomechanical DIR algorithm to locally accumulate the BT radiation dose on bladder in 47 cervical cancer patients receiving fractionated EBRT and two sessions of IGBT. They compared the DIR-based and DVH-based D2cc and D0.1cc parameters and found the DVH-based parameters to be (1.5 ± 1.8) % and (5.2 ± 4.2) % larger than the DIR-based parameters, respectively. Jamema et al. used an intensity-based and contour-based DIR for deformable dose accumulation for bladder and rectum in two-fraction IGBT treatment of cervical cancer. They compared the DIR-based D2cc values to the DVH-summed D2cc. Teo et al. performed a similar study on the data from 20 patients treated with EBRT and 4-5 fractions of IG-ICBT, using MIM’s DIR. They found a larger difference of (7.2 ± 6.3) % between DVH-based and DIR-based D2cc. Otherwise, there have been no studies reported in the literature that collectively: (i) look at 123 the accumulated dose in the ‘bladder wall’ rather than the ‘bladder surface’ or whole volume, (ii) apply deformable ‘point-set’ registration using ‘contour’ data rather than ‘image’ data, and (iii) study the relationship between the bladder ‘toxicity’ and its geometrically-localized accumulated dose. The present study aims to address all of the above. 4.2 Materials and Methods 4.2.1 Patient Materials and Treatment Characteristics The treatment planning data for 60 cervical cancer patients treated with fractionated EBRT plus five HDR CT-based IG-ICBT fractions (using a tandem/ring applicator for most cases and tandem/ovoid for a few) between August 2008 and June 2014, were studied retrospectively. The work reported herein was approved by the University of British Columbia BC Cancer Agency research ethics board, certificate number H1202320. Two groups of patients were included in this study: a. All cases with G2+ urinary toxicity observed (total of 17 subjects) b. A selection of the subjects treated during the same time period with G0-1 toxicity summing to 43 cases to make a sample size of 60, using two exclusion criteria: (i) less than 5 completed brachytherapy fractions and (ii) low-quality planning CT images. The mean age for this cohort was 48.9±15.6 years at the time of diagnosis. Based on the FIGO classification, 12 patients had stage IB1, 14 IB2, 5 IIA, 22 IIB, 3 IIIA, and 4 stage IIIB cervical tumours. Twenty-seven patients were treated for para-aortic nodal involvement using either a node removal surgery or nodal EBRT boost. While the bladder was mostly not included in the nodal EBRT boost volume, in cases that it was, the nodal boost dose was included in the mean bladder 124 EBRT dose. The prescribed minimum peripheral dose to the clinical target volume (CTV) was 600 cGy per HDR fraction for the vast majority of cases, with a few cases receiving 300-550 cGy in one, two or three out of five fractions due to OAR constrains. All patients underwent 25 EBRT fractions, with a total dose of 45 Gy prescribed to the ICRU reference point (ICRU Report 62, 1999). EBRT fractions were delivered on consecutive weekdays (except on BT treatment days) over five consecutive weeks. BT treatments were approximately 1 week apart and typically started 1.5-2.5 weeks after the first EBRT fraction. The EBRT planning was performed on a single planning CT image set and the patients were not imaged prior to each fraction. However, the patients were instructed to fill their bladder prior to treatment (by drinking 0.5 liter of fluid an hour before treatment and no voiding) and, therefore, roughly the same bladder filling regimen is assumed for all EBRT treatment fractions. The mean EBRT dose to bladder (including occasional dose from boost), based on the initial planning data, was on average 43.6 Gy. Due to the small inter-patient variations in the bladder EBRT dose (~5% based on the planning CT’s), the EBRT component was not included in the dosimetric analysis in this study. A Foley-Catheter was used to control bladder filling during HDR planning and delivery. Other than during a procedure (i.e. acquisition of the planning CT or treatment delivery), the bladder was allowed to drain freely into the urine bag. Immediately prior to each procedure, the drained bladder was filled with the desired volume of fluid using a syringe, after which the catheter tube was clamped in order to retain the fluid. A water-contrast mixture [3 cm3 of ISOVUE contrast agent (Bracco Diagnostic, Inc. Monroe Township, NJ] per 60 cm3 of sterile water) was used for 125 the planning CT, and water alone was used for treatment. The aim was to vary the bladder filling volume from fraction to fraction by alternating the amount of fluid injected between 60 cm3 and 300 cm3. In practice, this was not always achievable due to factors such as residual urine in bladder and intolerance of some patients to larger volumes. Also, in some patients, the bladder was purposely kept full to avoid higher doses to small bowel. The fluid volume was recorded for each fraction. 4.2.2 Urinary Toxicity Follow-up and Scoring Post-treatment late urinary effects were assessed based on the Late Effects in Normal Tissues–Subjective, Objective, Management and Analytic (LENT-SOMA) scale [179], collected from the patients over the course of their post treatment follow-up. The follow-up protocol consisted of an appointment at 6 weeks after completion of treatment, followed by visits every 2-3 months for 1 year, then every 3-4 months for year 2, and visits every 6 months for years 3-4, plus one more visit at year 5, but was not always fully complied with. This was mostly either because patients with symptoms were seen earlier than their scheduled follow-up appointment, or because patients who developed local, regional or metastatic recurrence were no longer given the LENT-SOMA questionnaire to complete. The last review of the LENT-SOMA questionnaires for this cohort was performed in November 2015. According to RTOG, late morbidity is defined as any event occurring or lasting 90 days after treatment initiation. The median follow-up for the subjects in this study was 23.6 months (2.3-66.2). Six control subjects had follow-up slightly below the minimum 3 months. Therefore, the analysis was also performed with exclusion of these cases. No baseline (pre-treatment) urinary 126 complication data had been collected for this cohort. The overall urinary toxicity scores were recorded as the average of the LENT-SOMA grades of all the side effects [161] and, therefore, ranged from 0.00 to 4.00. The subjects with overall toxicity scores >1.5 (i.e. Grade 2-4) were placed into the Case group and those with scores <1.5 (i.e. Grade 0-1) into the Control group. This choice of cutoff point for a clinically significant late toxicity was based on previous studies showing significant dose effect for complication grades ≥2 [46, 169]. The descriptive statistics for the dosimetric and volumetric parameters were recorded for the Case and Control groups separately. The age distribution of the Case and Control patients was also recorded. The toxicity analysis was then performed for all assessed parameters using Student t-test with the null hypothesis of equal means for the Case and Control groups. In order to check whether the requirement for normality of the Case and Control groups’ distributions was met for all parameters, the Shapiro-Wilk normality test was performed using SPSS Statistics software (IBM Corp. Armonk, NY, version 23). 4.2.3 Contour Processing For the purposes of this study, the bladder inner and outer surfaces were retrospectively re-contoured by a single observer on the planning CT images of all brachytherapy fractions using the MIM Maestro software. This was done to ensure the highest possible accuracy and consistency of contouring (as the contours used for treatment planning had been generated with tight time constraints by various observers), and to obtain both inner and outer bladder wall contours. As a quality assurance indicator, the volume of the bladder-wall tissue was monitored during the contouring process, with the assumption that it should remain reasonably constant for each patient 127 regardless of degree of distension – i.e. distension causes the wall to become thinner, but total wall volume stays constant. The points of attachment of the urethra, as well the two ureters for the cases with CT images of adequate quality, were identified at each scan, and used to monitor the accuracy of the deformable registration. These landmarks were well-suited to this purpose because they are located in the posterior part of the bladder, where the highest dose regions also tend to occur. As an additional quality assurance check, a random selection of the bladder-wall contours and landmark points, equivalent to 10% of the cases, were reviewed and approved by a radiologist. All dose calculations were performed using the Eclipse™ and BrachyVision™ treatment planning software (Varian Medical Systems Inc., Palo Alto, CA). For each patient, the dose matrix associated with each HDR fraction was imported into MATLAB with a minimum resolution (i.e. voxel size) of 2.5 mm by 2.5 mm by 1.25 mm. The bladder inner and outer surfaces were also imported into MATLAB, but as 3D point-sets in DICOM RT files. Since the resolution of the treatment plan CT images was typically ~0.5mm by 0.5 mm by 1.25 mm, the exported contour points were much denser within each slice (i.e. in the x-y plane) compared to between slices (i.e. along the z-axis). Therefore, in order to create a more homogenous distribution of the points in all directions and also reduce the computation time, the points in each inner and outer contour point-set were down-sampled to 4000 (4K) points (from ~20K-40K). The bladder-wall structure was then built by linearly interpolating between neighbouring inner- and outer-surface points so as to form four additional shells in between. The points were then down-sampled once more to 6000 (6K) points to form the final bladder-wall model. The 4K 128 and 6K sample sizes were selected empirically based on a trade-off between computation speed and accuracy in representing the structure. 4.2.4 CPD Deformable Point-Set Registration The detailed description, attributes and optimal parameter settings of the CPD toolbox for registering bladder contour point-sets were reported in Chapter 3. The same procedure and optimal settings of the registration parameters were used in this study, with one main difference. In this study, the reference point-set was the moving (i.e. deforming) structure, in contrast to the previous work, where the reference point-set was the target (i.e. fixed) structure. This was because in order to locally accumulate (i.e. register) the fraction doses on the contour points of a single (reference) structure, a correspondence must be found for all the points in that point-set. Therefore, CPD was used to register the reference point-set to the other four target point-sets, finding a match for every reference point to another point in each of the target point-sets. The TRE value for each anatomical landmark was calculated, for every registration, as the distance between the matched point in the target point-set and the known position of the corresponding landmark in reference. In order to remove the effect of distension of the chosen reference structure on the scaling of the TRE values, a normalized registration error (NRE) was defined using the equation, 𝑁𝑅𝐸 = 𝑇𝑅𝐸×𝐸𝑆 𝑟𝑎𝑑𝑖𝑢𝑠 𝑜𝑓 𝑠𝑚𝑎𝑙𝑙𝑒𝑠𝑡 𝑣𝑜𝑙𝑢𝑚𝑒𝐸𝑆 𝑟𝑎𝑑𝑖𝑢𝑠 𝑜𝑓 𝑟𝑒𝑓𝑒𝑟𝑒𝑛𝑐𝑒 𝑣𝑜𝑙𝑢𝑚𝑒 , (4.1) where ES is the equivalent sphere of the corresponding volume. The NRE values for all four registrations and across three landmarks were averaged to derive a NREavg. For each subject, the bladder point-set derived from each of the five fractions was tested as the moving structure and 129 NREavg was calculated for each scenario. The fraction yielding the smallest NREavg value was chosen as the reference fraction for the dose registration, in order to minimize the uncertainty in the registered dosimetric parameters propagated from the registration uncertainty. 4.2.5 Dose Accumulation After establishing the point correspondence, the point-dose values were calculated for each brachytherapy fraction by linear interpolation of the voxel-dose values, assuming locally-linear decay of dose in regions away from the steep dose gradients closer to the brachytherapy sources. This was confirmed to be a reasonable approximation based on the comparison of point dose values in MATLAB with doses calculated by BrachyVision for a series of test points. The registered point-dose values for all the points in the reference point-set were calculated using 𝑟𝐷𝑝 = ∑ 𝑑𝑝,𝑖5𝑖=1 , (4.2) where rDp is the registered dose at point p and dp,i is the dose at point p in fraction i. The registered doses were also calculated as biologically weighted, 2 Gy fraction equivalents (EQD2), derived using the linear quadratic model (α/β = 3 Gy) [23, 25], 𝑟𝐸𝑄𝐷2𝑝 = 𝑟𝐵𝐸𝐷𝑝/(1 +2𝛼/𝛽), (4.3) where 𝑟𝐵𝐸𝐷𝑝 = ∑ [𝑑𝑝,𝑖(1 +𝑑𝑝,𝑖𝛼 𝛽⁄)]5𝑖=1 [180]. (4.4) In order to create a ‘solid’ bladder-wall and accumulate the dose at each ‘voxel’ of this bladder-wall, the reference bladder-wall contour was rasterized into an image volume (binary mask of the bladder-wall) with the same size and resolution as the planning CT images. The accumulated voxel-dose values were then approximated by linear interpolation of the rDp values 130 of five nearest neighbouring points at the position of each bladder-wall mask voxel, creating a cumulative dosage mask (Figure 4.1). The same method was applied to find the EQD2 dosage mask. Note that, as discussed above, the bladder tissue volume is consistent and equivalent across different fractions, so the choice of reference fraction should not affect the absolute values of the computed dose-volume parameters. Figure 4.1 Cross section of a bladder-wall cumulative dosage mask, showing the distribution of voxel doses as different shadings in the mask, with the brighter regions having a larger registered dose than the darker ones. 4.2.6 Registered Dosimetric and Volumetric Parameters Because the registered cumulative dose distributions often exhibited separated, non-contiguous high dose regions, the search for critical high-dose volumes on the bladder wall that might explain the observed severe bladder toxicity in some patients included dosimetric parameters calculated for contiguous high dose regions. However, defining a biologically meaningful contiguous high dose volume is not always a simple task, particularly in the presence of "noise" caused by structure contours that are not smooth, and in situations where two (or more) separate high-dose regions are beginning to merge (as shown later in Figure 4.4). 131 Thus, a novel region-growing algorithm was developed and evaluated for the purpose of finding contiguous dosimetric parameters. This method uses a software implementing a random walker algorithm, called Fast Random Walker (FRW)1, to enforce smoothness on the calculated regions around the high-dose voxels, using the following algorithm to find convex-shaped contiguous high-dose volumes of desired sizes: Step 1. Find the voxels, d_h and d_l, containing the highest and lowest doses, respectively. Step 2: Apply FRW to assign probabilities to each voxel in the bladder wall mask, with d_h regarded as a “source” and assigned probability 1, and d_l a “sink” worth probability zero. The FRW algorithm ensures a smooth transition in probabilities. Step 3. Compute a path from d_h to d_l by descending down (i.e. along the negative gradient) the probability map. This path is effectively a regularized sorted list of the voxels by dosage from highest to lowest. Step 4. Compute the desired Dncc for this iteration by taking the top N voxels (corresponding to a n-cc volume) in the sorted list and picking the lowest value. Step 5. Mask out the top N0 voxels that amount to >= 1cc volume temporarily and find the next highest dosage point d_h2. Step 6. Unmask the masked values and repeat steps 2-4 for d_h2. Step 7. Compare the new Dncc from d_h2 to the previous Dncc volume. Step 8. If the Dncc is higher than the previous value, then the 1cc volume of all computed d_h points are masked and steps 4-6 are repeated until otherwise or if the preset number of trials is exhausted. 1 Fast Random Walker by Shawn Andrews et al. (2010), available at http://fastrw.cs.sfu.ca/. 132 Step 9. The highest Dncc is returned. The spatial regularization is controlled by a regularization weight parameter (beta), which can be adjusted, where greater regularization weights lead to a smoother, more circular surface area of the calculated volume. Three different levels of regularization were applied and compared visually to the corresponding non-contiguous volumes. A preferred regularization level was chosen and used to calculate the minimum registered dose to the contiguous high-dose n-cc of the bladder wall, rDncc, for 1 and 2 cc volumes. In addition, conventional DVH-based rDncc were calculated for volumes recommended by GYN GEC-ESTRO guidelines for OARs (n = 0.1, 1, and 2 cc; optional 5 and 10 cc) [43], with the exclusion of the 0.1 cc volume because its dimensions were comparable to the registration uncertainties. These volumes were calculated by sorting the bladder wall dosage mask and picking the top N voxels, making up the desired volume. Therefore, this method is a non-contiguous method. Dose-volume analysis was limited to the non-contiguous method for the larger volumes because the results obtained by FRW were more sensitive to the effects of regularization at larger volumes, and because the conventional DVH-based rDncc represented a contiguous volume for most patients at 5cc and all patients at 10cc. Also calculated were the non-contiguous volume of the bladder wall receiving a minimum registered dose of n Gy, rVnGy (evaluated for n = 3, 5, and 10 Gy), rEQD2ncc parameters for 1and 2 cc (contiguous), and 5 and 10 cc (DVH-based) volumes in the EQD2 dosage mask, and the mean and integral dose to the bladder wall. Conventional plan-summed parameters, sDncc, were calculated using the DVH data from the treatment planning software (TPS) for each fraction and 133 summing across all fractions. Because they were based solely on DVH data, sDncc were not spatially defined. The minimum cumulative dose to small registered volumes around the bladder neck (BNDncc), left and right ureter orifices (LUDncc and RUDncc, respectively) were also determined. Table 4.1 shows a list of all the dosimetric and volumetric parameters calculated in this study, with their definitions. 4.2.7 Validation Tests In order to validate different steps of the complex calculations performed in this study, multiple validation tests were completed. Most of the validations were in form of benchmarking the computation outcomes in MATLAB against the corresponding quantities calculated by the MIM Software. Bladder Wall Mask The volume of the calculated bladder wall mask, for the reference fraction of each subject, was calculated in MATLAB and compared to the volume of the corresponding Bladder_wall contour structure created in MIM Software, by using the Boolean operator tool to subtract the Inner_wall contour from the Outer_wall contour. While some discrepancies were anticipated due to differences in the way MIM creates the wall contour and calculates its volume and the method that is used in MATLAB to rasterize the wall mask from the calculated wall point-set and sum the mask voxel volumes, a close match was achieved between the volumes. The difference between the two volumes across all subjects ranged between 0.04-1.9 cc, with an average of 0.6±0.3 cc. 134 Table 4.1 Bladder wall dosimetric and volumetric parameters and their definitions listed at the top and landmark-based parameters at the bottom. Parameter Definition rDncc Minimum registered dose to a high dose n-cc volume of bladder wall (n=1*, 2*, 5†, 10†) rEQD2ncc Minimum registered EQD2 dose to the high dose n-cc volume of bladder wall (n=1*, 2*, 5†, 10†) rVnGy Volume of bladder wall receiving a minimum registered dose of n Gy (n=3†, 5†, 10†) Mdose Mean registered dose value in bladder wall Idose Integral registered dose (∑ 𝑑𝑖𝑚𝑖𝑖 , where mi is the mass of each voxel and di is its registered dose) sDncc Summed Dncc from DVH data (∑ 𝐷𝑛𝑐𝑐,𝑖5𝑖 =1 , where i=fraction number) sEQD2ncc Summed EQD2ncc from DVH data (∑ 𝐸𝑄𝐷2𝑛𝑐𝑐,𝑖5𝑖 =1 , where i=fraction number) BNDncc Minimum dose to registered n-cc volume of bladder wall centered at bladder neck (n=0.1*, 0.5*, 2*) LUDncc Minimum dose to registered n-cc volume of bladder wall centered at left ureter orifice (n= 0.5*, 2*) RUDncc Minimum dose to registered n-cc volume of bladder wall centered at right ureter orifice (n=0.5*, 2*) * Volume calculated using contiguous method † Volume calculated using non-contiguous method, based on sorting voxel doses Dose Interpolation and DVH Calculation Two linear interpolation steps were performed during the MATLAB-based dose calculation procedure, which introduce some approximation into the calculated dosage mask. The first interpolation was performed to re-sample the dose matrix to match the resolution of the image space, since the voxel size of the default dose matrix exported from BrachyVision in the x-y direction (typically 2.5 mm × 2.5 mm) is much larger than the image volume voxels (typically 0.4 mm × 0.4 mm or 0.5 mm × 0.5 mm). This interpolation was validated by comparing different reference point doses calculated in MATLAB to those in BrachyVision, for all subjects. It was ensured that the MATLAB-based values were within 3% of the TPS-based values. For the cases were this limit was exceeded, a finer dose matrix was exported to provide more accurate interpolation. Finer dose matrices are large files that require a lot of memory to store and process 135 and therefore, were not practical to use for all the cases. The second interpolation was performed after establishing the point correspondence from registration, where the bladder wall contour point-dose values were calculated for each fraction by linear interpolation of the voxel-dose values. This interpolation was validated together with the non-contiguous DVH calculation method (described in Section 4.2.6) by comparing the per-fraction and fraction-summed DVH values for 1cc and 2cc volumes (Table 4.2). It can be seen in Table 4.2 that the maximum fraction-summed deviation observed was 0.61 Gy, i.e. 2% of the total prescription dose. Contiguous Volume Calculation Different aspects of the region growing algorithm were evaluated. The calculated regions were inspected visually, as discussed in Section 4.3.2, to ensure the regions encompass a meaningful region of the bladder wall. Also, the volumes of the calculated hotspot regions were calculated, by accumulating the voxel sizes in the calculated mask, and cross-checked with the nominal values. The contiguous registered dose parameters were compared against both the corresponding registered non-contiguous and unregistered DVH-summed values to ensure that they are always smaller, which was the predominant outcome. The results for the differences between contiguous and both non-contiguous and DVH-summed parameters are presented in Section 4.3.4. 4.2.8 Volume Variation Effect The intra-patient variation of the bladder volume (i.e. wall tissue plus lumen) over the five fractions, VolVar, was defined as (maximum-minimum)/minimum of the filling volumes. Also, in 136 Table 4.2 The differences between MATLAB-based and TPS-based DVH values (in Gy and as a percentage of the prescription dose) for 1cc and 2cc volumes, across all fractions of all subjects. Per fraction Fraction summed 1cc 2cc 1cc 2cc Mean 0.05 (0.8%) 0.05 (0.8%) 0.16 (0.5%) 0.18 (0.6%) Median 0.05 (0.8%) 0.04 (0.7%) 0.14 (0.5%) 0.15 (0.5%) Stdev 0.03 (0.5%) 0.03 (0.5%) 0.13 (0.4%) 0.13 (0.4%) Min 0.00 (0.0%) 0.00 (0%) 0.00 (0.0%) 0.01 (0.0%) Max 0.48 (8.0%) 0.33 (5.5%) 0.61 (2.0%) 0.46 (1.5%) Note: For the per-fraction data, the min and max values are across all individual fractions of all subjects, whereas the mean values are across the average of five fractions for each case. order to assess the anatomical movement of the hotspots on the bladder, the location of the hottest point on each fraction was mapped to the reference bladder using the registration outputs. Then, the mean and maximum of all the distances between different pairs of the hottest points for the five fractions were calculated. To evaluate the effectiveness of varying the bladder filling in moving the hotspot regions on the bladder and reducing urinary toxicity, a linear regression analysis was performed between (i) the mean and maximum hotspot distances and VolVar, and (ii) urinary toxicity score (on a 0.00-4.00 scale) and VolVar. For both analyses, the p-values for the non-zero slope of the linear regression were calculated. 4.3 Results 4.3.1 Registration and Dose Accumulation All registrations were visually inspected and found to be reasonable (i.e. no extreme outliers could be identified). For most of the cases, the overall match was quite satisfactory with small apparent discrepancies (see Figure 4.2 for an example). TRE data was not available for one out of 137 the 60 subjects since no landmarks were identifiable for that subject. For 9 other subjects, the ureter landmarks were not identifiable and therefore, only bladder neck was used for TRE calculations. Quantitatively, the NREavg values (for the optimal reference fraction, as described previously) ranged from 2.4 mm to 9.2 mm across all 59 subjects with available data, with a mean of 4.8±1.5 mm. The corresponding mean TREavg value was 5.5±2.0 mm. The per-fraction and registered dose distributions on the bladder-wall contour point clouds, for the same case as in Figure 4.2, are shown in Figure 4.3. It can be seen that the locations of the hotspots varied significantly for most of the fractions and, therefore, the location of the summed hotspot was quite different from the per-fraction hotspot locations. This confirms the conjecture that simply adding the unregistered dosimetric parameters would not provide an accurate estimate for the magnitude or spatial distribution of the cumulative dose over all fractions. 4.3.2 Contiguous and Non-Contiguous Volume Calculation Table 4.3 shows the descriptive statistics for contiguous rD2cc values across all cases, estimated using three different regularization schemes. Figure 4.4(a-c) and (d) show the surface doses of contiguous and non-contiguous 2 cc high-dose regions, respectively, for a sample case that represented one of the more difficult scenarios for defining a contiguous high-dose volume. The contiguous 2cc volumes shown in the figure were defined using the three regularization schemes and the non-contiguous 2-cc region was defined using the DVH-based method. It can be seen that there are two separate hotspots at small-volume thresholds for this case. As a result, the non-regularized scheme generates a highly irregular shape in the form of two islands connected by a thin strip, which presumably does not represent a biologically meaningful region. 138 Figure 4.2 Example of the registration outcomes for registering the reference bladder-wall contour point-set (FR4, shown in red) to two larger fractions (FR2 and FR3, shown in blue): ‘Before’ and ‘After’ scenarios. On the other hand, the tightly regularized regimen outputs a well-rounded region. However, the strict regularity constraint on this region prevents it from closely following the spread of the dose distribution, producing underestimated rDncc values, as shown in Table 4.3. The loose regularization, on the other hand, provided a moderate degree of smoothness and reasonable rD2cc values. This regularization scheme was therefore used for the calculation of the contiguous rDncc values presented in Table 4.4. It can be seen in Figure 4.4(e) that, even in the presence of two separate hotspots the 5cc has become contiguous. Moreover, the effects of 139 regularization were found to be more pronounced for the larger volumes. Therefore, non-contiguous rDncc values (equivalent to the values for non-regularized contiguous volume) were used for 5 and 10 cc volumes. Table 4.3 Descriptive statistics for contiguous rD2cc values based on three different regularization regimens used to calculate contiguous volumes. All units are in Gy. Tight Regularization Loose Regularization No Regularization Mean 17.1 19.3 22.9 Median 17.3 19.5 23.3 Stdev 4.0 4.3 5.1 Min 8.2 8.8 10.3 Max 27.1 29.1 37.7 4.3.3 Urinary Toxicity, Age, and EBRT dose The observed bladder LENT-SOMA toxicity scores across all subjects ranged from 0.00 to 3.57 (i.e. Grade 0-4). There were 17 subjects in the Case group and 43 Controls. The toxicity scores in the Case group ranged from 1.57 to 3.57, with an average of 2.02±0.61, and the Controls ranged from 0.00 to 1.43, with an average of 0.81±0.43. The difference between the toxicity scores of Case and Control groups was highly significant at p<0.000001. The mean EBRT dose (mEBRT) received by the bladder was determined based on the initial planning data. The coefficient of variation of the mEBRT values was only 5% across all cases. The mean mEBRT values for the Control and Case groups were 43.3±9% Gy and 44.2±3% Gy, respectively. The mean age for Case and Control groups was 52.1±17.5 (26-83) yrs and 47.6±14.9 (19-82) yrs, respectively. The difference between the Case and Control means for age was not significant (p=0.13). 140 Figure 4.3 The per-fraction and accumulated point-dose distributions, on the same geometrical scale, for the representative case in Fig. 4.2. The fraction doses are locally accumulated on the reference fraction (FR4), showing the variations in the locations of hotspots, among the fractions and also the accumulated result. 4.3.4 Registered Dose Parameters The normality requirement was met for all parameters in both groups, except for the Case group BND0.5cc and BND2cc. The mean registered non-contiguous rD1cc and rD2cc values were 24.9±7.3 Gy and 23.1±6.9 Gy for Case, and 25.2±4.9 Gy and 23.1±4.5 Gy for Control, respectively. The registered non-contiguous rD1cc and rD2cc values were, respectively, (5±3) % and (3±2) % lower 141 Figure 4.4 Sample hotspot regions on the bladder wall for (a) tightly regularized contiguous 2-cc volume, (b) loosely regularized contiguous 2-cc volume, (c) non-regularized contiguous 2-cc volume, (d) non-contiguous 2-cc volume, and (e) non-contiguous 5-cc volume. (f) Corresponding dose distribution. than the corresponding unregistered DVH-based sDncc values, across all subjects. The contiguous rD1cc and rD2cc values were (17±5) % and (18±6) % smaller than the equivalent sDncc parameters. The contiguous rDncc values for 1 and 2 cc were, on average, (13±5) % and (15±5) %, smaller than the corresponding non-contiguous parameters, suggesting that multiple hotspots with volumes of at least 2cc were relatively common in the registered dose distributions. The contiguous rEQD2ncc values were, on average, (26±10) % and (27±11) %, (a) (b) (c) (d) (e) (f) 142 smaller than the equivalent DVH-based sEQD2ncc parameters for 1 and 2cc volumes, respectively. Table 4.4 shows the descriptive statistics for the registered dosimetric and volumetric parameters for the Case and Control groups. For all parameters, except one, the mean value for the Case group was higher than the Control group, as expected. V3Gy demonstrated a significant difference between Case and Control means based on a two-sample t-test at a 5% level of significance (p=0.014). rV5Gy and Idose showed a significant difference at the 10% level. The t-test analysis was also performed for the sDncc and sEQD2ncc parameters, with the smallest p-value being observed for sEQD210cc (p=0.097). The toxicity analysis was repeated excluding the six Control subjects with <3 months follow-up and very similar results were obtained (p=0.01 for rV3Gy and non-significant slightly increased p-values for the rest,). The increase in the p-values was expected because of the smaller sample size. Table 4.5 shows the landmark-based dosimetric parameters for the bladder neck, right and left ureter orifices. None of these parameters demonstrated a significant difference in the Case and Control means based on a t-test analysis. In fact, the average of the dosimetric parameters for smaller volumes around the ureteral orifices in the Case group is smaller than the Control group. As mentioned previously, no landmarks were identifiable in one subject (in Case group) and the ureters landmarks were not identifiable in 9 other subjects (2 Case and 7 Control). 4.3.5 Volume Variation The per-patient, inter-fraction variation of bladder-wall tissue volume was low with an average of (5.4±1.6) % and a maximum of 8%, which suggested that the internal and external bladder contours were consistent from fraction to fraction. The average bladder-wall thickness 143 over all fractions of all subjects was found to be 3±1 mm. The maximum intra-patient volume variation (maximum volume – minimum volume) ranged between 4.5-427 cc, with a mean of 152±112 cc. Volume variation, as measured by VolVar, ranged between 3.7% and 353.3% across all cases, with an average of (107±87) %. Figure 4.5 shows urinary toxicity score versus VolVar in a scatter plot. The linear regression fit is also shown, where the slope of the regression line is −0.23 (p=0.03). The ranges of the mean and maximum hotspot separation distances across all cases were 3-38 mm and 5-67 mm, respectively. It must be noted that these separation distances are Euclidean distances and not a geodesic shortest path on the surface on the bladder. Therefore, these distances may be an underestimation of the actual distances between the hotspots along the bladder surface, particularly for cases where the hotspots occur in different lateral horns of the bladder. The slopes of the regression lines for the mean and maximum hotspot separations as a function of VolVar were, respectively, 3.8 mm (p=0.03) and 7.1 mm (p=0.002). 4.4 Discussion In this study, the deformable PSR utility of the CPD toolbox was applied to bladder-wall contour point sets of five treatment fractions for 60 cervical HDR brachytherapy patients. The quality of the registration outcomes was confirmed both visually and through the landmark registration errors calculated for three anatomical points on the bladder surface. The PSR results were used to accumulate the dose from all fractions onto the bladder wall of a reference fraction, and registered dosimetric and volumetric parameters were computed. The unregistered, DVH-summed dosimetric parameters commonly used in clinical practice were also calculated and com- 144 Figure 4.5 The scatter plot of toxicity score versus volume variation in terms of VolVar=(Max-Min/Min) of fraction volumes for the sixty subjects in the study, also showing the linear regression fit. Note that VolVar is in absolute values and not in percentage. pared to the registered-dose parameters. The dosimetric and volumetric parameters were then evaluated as predictors of urinary toxicity grades. Multiple validation tests for different steps of the complex dose accumulation procedure were carried out to ensure the accuracy and quality of the outcomes. All tests produced satisfactory results with reasonable uncertainties. The existing small discrepancies were expected due to factors such as the effects of rounding, down-sampling, point-dose interpolation, slightly different approaches in calculating DVH, and minor discrepancies between bladder wall volumes calculated by TPS and MATLAB (up to ~2cc). The registration uncertainty achieved in this study (4.8±1.5 mm) is amongst the best reported values for bladder deformable registration [181, 110, 182]. It was shown previously (Chapter 3 and [182]) that CPD has superior registration accuracy compared to five other PSR techniques, y = -0.2356x + 1.40590.000.501.001.502.002.503.003.504.000.00 0.50 1.00 1.50 2.00 2.50 3.00 3.50 4.00Toxicity ScoreVolVar145 Table 4.4 The Case (CS) and Control (CN) statistics for different dosimetric and volumetric parameters based on registered dose, and the bladder volume variation in percentage. The p-values less than 0.1 are highlighted in bold text. rD1cc* (Gy) rD2cc* (Gy) rD5cc (Gy) rD10cc (Gy) rEQD21cc* (Gy) rEQD22cc* (Gy) rEQD25cc (Gy) rEQD210cc (Gy) rV3Gy (cc) rV5Gy (cc) rV10Gy (cc) Idose (Gy) Mdose (Gy) VolVar (%) CS CN CS CN CS CN CS CN CS CN CS CN CS CN CS CN CS CN CS CN CS CN CS CN CS CN CS CN Mean 21.7 21.8 19.8 19.1 19.0 18.5 14.1 13.0 34.7 33.1 29.5 27.0 27.9 26.0 17.7 15.3 50.9 42.7 35.5 30.5 16.0 13.9 484 414 9.0 8.9 79.5 117.3 Median 21.9 21.5 20.3 19.4 19.3 18.9 13.5 12.6 34.0 33.6 28.3 25.9 27.4 27.3 16.3 15.4 50.2 40.8 32.6 31.6 15.1 14.1 431 419 9.1 8.7 62.1 87.8 Stdev 6.1 4.1 5.7 3.7 5.8 4.2 4.6 3.8 15.0 9.6 12.5 8.6 12.8 9.6 8.3 7.0 13.1 8.9 13.0 9.2 7.4 5.0 173 122 2.4 2.0 69.4 91.3 Min 9.9 12.8 8.8 10.9 8.5 10.9 6.6 6.9 9.7 13.5 8.4 11.0 8.2 8.6 5.7 4.4 27.2 23.9 17.2 15.1 2.4 5.7 233 195 4.9 5.2 8.5 3.7 Max 32.5 30.1 29.1 28.4 29.7 29.5 21.4 24.5 68.7 51.6 55.0 50.9 55.4 53.8 31.2 40.0 75.9 62.5 61.7 55.6 30.1 29.7 778 829 12.7 13.3 254.4 353.3 p-value 0.482 0.317 0.375 0.188 0.343 0.223 0.294 0.156 0.014 0.080 0.154 0.071 0.493 0.046 * parameter is calculated for a contiguous dose volume Table 4.5 The Case and Control statistics for the bladder neck and left and right ureters’ orifices. BND0.1cc (Gy) BND0.5cc (Gy) BND2cc (Gy) LUD0.5cc (Gy) LUD2cc (Gy) RUD0.5cc (Gy) LUD2cc (Gy) CS CN CS CN CS CN CS CN CS CN CS CN CS CN Mean 7.3 6.1 6.4 5.4 5.3 4.7 12.3 13.0 8.5 7.4 12.9 13.7 8.2 7.9 Median 6.2 5.3 5.4 4.7 4.6 4.2 11.4 12.6 7.6 6.8 13.8 14.5 9.0 7.6 Stdev 3.6 2.5 3.4 2.3 2.6 2.1 3.3 4.6 2.9 3.2 3.9 4.3 2.2 2.7 Min 3.6 3.0 3.2 2.7 2.8 2.3 6.9 5.1 4.8 3.2 4.9 5.6 3.9 3.9 Max 15.8 13.2 14.9 12.6 11.5 11.5 17.6 22.3 12.8 17.2 18.1 20.5 11.1 13.6 p-value 0.106 0.129 0.188 - 0.128 - 0.322 Note: For parameters where the mean for Control group was higher than the Case, no p-value was calculated.146 including a widely-used method called thin-plate spline robust point matching (TPS-RPM). In a recent study, Chen et al. (2016) developed an improved TPS-RPM method with local topology preservation, for the purpose of bladder dose accumulation in cervical cancer HDR brachytherapy. They tested their method on three different sets of data: synthetic data, porcine bladder phantom (using fiducials), and patient data (using bladder neck landmark). While they reported residual distance errors (RDE) of 3.7±2 mm on average, for both porcine bladder and patient data, the amount of variation in bladder size in that study (80-300 cc total volume) was much less than the variations tested in the current study (i.e. 71-562 cc). Given that registration error has been shown to become larger with the variation in structure dimensions [181, 110, 182], a slightly larger registration error is expected in our work. According to the phantom-based PSR validation study in Chapter 3, the registration uncertainty for CPD is more tangential than transverse. Therefore, the average NRE error reported in our current clinical study can be compared to the planar dimensions of the measured hotspot volumes. Considering the average bladder-wall thickness of 3 mm, the circular diameter of the smallest high-dose volume considered (i.e. 1cc) is about 20.3 mm. This size is more than 4 times the average registration error of 4.8 mm, suggesting that although this registration uncertainty might cause some smearing and blurring of the edges of the hotspot, the hotspot regions from each fraction would still be reasonably well aligned. Moreover, since the dose fall-off is more gradual in the bladder than it is in regions closer to the radioactive source delivery, there are no drastic variations in dose values across a 5-mm distance on the bladder surface. Furthermore, as the hotspot volume increases to 2, 5 and 10 cc, 147 the circular diameters expand to 29, 46, and 65 mm, diminishing the impact of the registration uncertainty even further. Revisiting Figure 4.3, it is clear that accumulation of dose, after performing registration, is far more accurate than assuming total overlap of the hotspots and simply summing their DVH dose parameters, despite minor uncertainties in registration. Based on the results of this study, the unregistered, summed D2cc overestimates the cumulative contiguous hotspot dose by (18±5)%. This difference was much lower for the registered non-contiguous dose, i.e. (3±2)%, which is comparable to the results previously reported by Andersen et al. (2013), Jamema et al. (2015), and Teo et al. (2015), i.e. (1.5±1.8)%, (5.2±5.1)%, and (7.2±6.5)%, respectively. As previously discussed in Andersen et al. (2013), some overestimation of dosimetric parameters from DVH summation is expected because adding, for example, the per-fraction D2cc values assumes that the same 2 cc on the bladder wall is receiving the highest dose over all treatment fractions. This assumption is referred to as the “worst case assumption” in the GYN GEC-ESTRO recommendations [183, 184]. Application of OAR dose thresholds based on unregistered DVH-summed parameters has reduced the rate of observed urinary toxicity in some institutions [46, 169], but it failed to explain the observed toxicity in the cervical cancer patients treated with HDR brachytherapy at our center, perhaps because hotspots were moved around by the practice of varying the bladder filling between treatment fractions. Being able to correctly accumulate the dose across treatment fractions in this situation may provide a more meaningful dosimetric analysis. An analysis was therefore performed for 13 different registered dose parameters by separating the subjects into Case and Control groups based on their overall bladder toxicity grades. The 148 normality requirement for Student t-test was met for all bladder-wall parameters. Only one of the parameters, rV3Gy, yielded a significant rejection of the null hypothesis (p=0.01). A possible underlying reason for this finding might be that the urinary endpoints that more effectively influence the average toxicity score for our data set (i.e. frequency and urgency) are more strongly associated with larger volumes of the bladder receiving accumulated doses beyond a certain threshold. Further investigation of this effect would benefit from inclusion of baseline urinary complication data in the analysis in order to account for the effect of higher baseline values that might exist for some patients. Overall, however, it was found that 17 out of 20 dosimetric parameters evaluated had higher means for the Cases compared to the Controls. Based on binomial probability test, the probability of finding the mean in the Case group higher than the Control, for 17 out of 20 parameters, by random chance is 0.001 (equivalent to p=0.004, using corresponding z-score). This is a significant finding, but it must be noted that these parameters are not all independent (e.g. there is a correlation between rD1cc and rD2cc). The statistical power (1-β, i.e. the probability that the test correctly rejects the null hypothesis) based on the Case and Control sample sizes for α=0.05 was 0.95 for rV3Gy, which is reasonably strong considering the accepted standard minimum of 0.80 for avoiding a type II error (i.e. false negative), according to Cohen [185]. However, given the total number of parameters tested, we recognize that our findings could be subject to the multiple comparisons problem [186], causing a type I error (i.e. false positive). Therefore, the Holm's step-down approach [187] for controlling the family wise error rate (FWER) [188] was applied, and it determined that no true rejected null hypotheses can be established from the results of this study. However, this FWER calculation 149 would represent the worst case scenario of considering all 13 dosimetric parameters to be independent, which is not the case as stated above. Nevertheless, the findings at this stage should be regarded as hypothesis-generating only. Analysis of a larger data set would be required in order to reach more definitive conclusions. Since six Control subjects had follow-up slightly shorter than 3 months, which is the minimum time required for assessing late complications according to RTOG, the toxicity association analysis was repeated excluding these cases and similar results were achieved. However, for more robust and reliable toxicity data, a cohort including cases with longer follow-up is desirable. The results reported in Table 4.4 show a few apparent anomalies, which appear mostly as rD2cc ≤ rD5cc (rather than rD2cc > rD5cc, as would be expected from conventional DVH analysis) for some of the median, maximum and minimum values. Because these apparent anomalies occur for statistics that are determined from single data points, one contributing factor is re-assortment of the data order for different parameters. For example, the Min(rD2cc) and Min(rD5cc) in the Control group do not come from the same case. Statistics computed from all data values (i.e. the mean), follow the expected pattern across all parameters for both Case and Control groups. Additional contributing factors to these anomalies include (i) uncertainties in the computed values caused by noise and irregularity in the contours and dose distribution, (ii), the effect of regularization in the contiguous volume calculation, and (iii) high dose gradients present in the accumulated dose distribution, especially for cases where the hotspots of different fractions mostly overlap. In the presence of these high gradients, a slight displacement in the edge of the calculated region can result in rather drastic changes in the calculated Dncc value. Finally, because rD2cc is 150 based on a contiguous volume and rD5cc is not, it is possible for rD2cc to be less than rD5cc. For example, if there are three or more separate high-dose regions, the contiguous rD2cc may be based on one of these, while the non-contiguous rD5cc is defined by a slightly higher isodose line, as shown in Figure 4.6. A separate analysis was also performed, looking at the minimum dose received at small volumes around the anatomically important points on the bladder, i.e. bladder neck and ureteral openings. Although no significant difference in the means of Case and Control was observed for bladder neck dose parameters, the mean of the Case group was higher than that of the Controls. The p-value for the 0.1cc volume around the bladder neck was the smallest amongst all the landmark-based dose parameters. The fact that an inverse relationship was observed between Case and Control means for LUD0.5cc and RUD0.5cc could be partly due to the fact that the ureteral landmarks were not available for 10 out of the 60 patients, where 3 were Cases, making the total available Case subjects as low as 14. The normality requirement for t-test analysis was not met for the Case group for bladder neck 0.5cc and 2cc parameters, but this was not of importance here since no significant findings were achieved for these parameters. In this study, the effectiveness of varying the bladder filling volume across different HDR BT treatment fractions in moving the location of the hotspots and reducing urinary toxicity was also evaluated. It was observed that varying the bladder filling volume correlated significantly (p=0.03) with both hotspot separation and urinary toxicity. For the urinary toxicity, the slope of the linear fit was -0.23 with VolVar expressed in absolute values. This means that if the maximum filling volume of bladder is increased three-fold compared to the minimum, the urinary toxicity score is only reduced by 0.73 (0.23×3). Therefore, the result of this study suggests that varying the bladder 151 Figure 4.6 Schematic showing a case with 3 hot spots (A, B, and C): region A represents a contiguous 2cc region enclosed by the 28 Gy isodose line, whereas regions B and C enclosed by the 28 Gy isodose line are slightly smaller than A (therefore the contiguous D2cc for these two regions would be slightly lower than 28 Gy and would not be selected as the representative contiguous D2cc). The non-contiguous D5cc has contributions from all three hot spots, and is represented by the dashed red 29 Gy line. filling regimen has the potential to improve urinary toxicity by moving the location of the maximally irradiated regions on the bladder across different brachytherapy treatment sessions. However, this improvement is rather minimal. Although there was no difference observed in the bladder mean EBRT dose from the planning data between the Case and Control groups, the EBRT dose can play an important role in the observed urinary late effects. Since there was only a single planning image available for each patient and the patients were not imaged prior to each EBRT fraction, the intra-fraction variations could not be accounted for. Therefore, it was not possible to register the EBRT dose across treatment fractions. However, the mean EBRT dose was more than 95% of the prescription dose in most cases, implying that the entire bladder was covered in most of the cases, rendering the effect of organ motion across EBRT fractions less important. Moreover, considering that a reasonably consistent bladder filling regimen can be assumed for EBRT treatment fractions, the A B C 28 Gy 29 Gy 152 inter-fraction variations are believed to be small compared to the variations seen for the brachytherapy treatments. That being said, no quality assurance or auditing of the bladder filling during EBRT was performed. A rigorous study of the dose-effect relationship of the urinary late effects would require data with repeated EBRT imaging in order to find the total accumulated dose distributions on the bladder from both EBRT and brachytherapy. Moreover, as mentioned above, baseline urinary complication data and sufficient follow up is required for a more accurate evaluation of post-treatment toxicity. Such information could potentially lead to a more definitive association between urinary toxicity and registered dose parameters. 4.5 Conclusion The results of this study suggest that that deformable point-set registration appears to be a viable tool for more accurate dose accumulation in fractionated image-guided brachytherapy, and that DVH-summed parameters calculated from unregistered data sets overestimate the bladder dose in intracavitary cervical HDR brachytherapy by at least 3%. This overestimation increases when contiguous high dose volumes are considered and when large variations in bladder filling volume are present. While more research is required in order to fully characterize the dose effect relationship for registered dose accumulations, some evidence suggesting that volumetric registered-dose parameters with lower dose thresholds may have the potential to explain some of the severe urinary toxicity observed in certain patients was found. 153 Chapter 5: Analysis of Specific Urinary Symptoms 5.1 Introduction Although urinary dysfunction has long been recognized as one of the side effects following pelvic radiotherapy, the pathophysiology of radiation-induced damage to the bladder and the rest of the lower urinary tract has still not been fully understood [158]. However, based on several animal studies using rat and mouse bladder, it has been established that the cell renewal rate in bladder transitional epithelium (a.k.a urothelium) is low and therefore, stimulated proliferation does not start until months after irradiation [189, 190, 191]. As a result, late radiation effects are commonly observed in the bladder and urethra. Radiation damage may affect the membrane lining at the surface of the bladder that is resistant to urine, causing irritation and damage to the underlying tissue layers, which may result in urinary sequelae such as infection, pain, and hematuria. The smooth muscle fibers are also prone to radiation damage, causing edema and cellular destruction. This may change the bladder’s capacity, creating symptoms such as frequency and nocturia. Damage to the bladder vasculature may cause vascular occlusion and ischemia, which can lead to late bladder fibrosis and reduction of bladder capacity. Radiation to the urethra and bladder neck may also be responsible for some of the observed late urinary effects [158, 159, 160, 146]. In particular, incontinence might more significantly correlate with dose to bladder neck, compared to other endpoints, such as frequency, urgency, and hematuria, which are most likely caused by radiation-induced cystitis [158, 146]. Hematuria can also be speculated to correlate with high dose to small volumes of the bladder wall. However, there 154 has not been enough investigation to strongly support this theory. When recording the urinary sequelae in practice, based on the LENT-SOMA system, the toxicity is summarized in one grade, which might be the average or maximum of individual symptom scores [161, 46, 162]. However, considering that different urinary dysfunction endpoints might have different underlying causes and mechanisms, combining individual toxicity scores might inhibit our understanding of the radiation dose effects in urinary morbidity. In fact, the focal and global injuries in bladder are believed to have different dose-volume relationships [158]. Moreover, symptoms that are grouped together as one toxicity grade, which can be dichotomized at a certain cut-point, have different relative importance and effect on the quality of life of the patient. For example, a Grade-3 hematuria (frequently having blood in the urine with the presence of clots) has a considerably different clinical impact than Grade-3 frequency (need to void every 1-2 hour). While even a Grade-1 hematuria needs direct medical attention, Grade-3 frequency might be considered “normal” for a lot of women. The drawback of “aggregating” large volumes of morbidity information into a single “statistic” has also been described by Rosewall et al. [164] both as “an imprecise method of describing the actual type of dysfunction experienced by a patient” and “detrimental when attempting to find a link between dysfunction and radiotherapy dose”. The same argument has also been promoted by Bentzen et al. [192], where the idea of averaging the individual LENT-SOMA organ scores is strongly rejected. Furthermore, it has been shown that different urinary endpoints have different temporal variations [193, 194, 195], which would affect the consistency of aggregated grades with different follow-up intervals. 155 It is therefore worthwhile to study each symptom individually in relation to radiation dose. While there have been a few studies exploring the dose relationship for individual urinary symptoms, such as frequency and nocturia, following prostate radiotherapy, [193, 194], there have not been any studies investigating dose effects of individual urinary endpoints in HDR BT for cervical cancer, to the best of author’s knowledge. There have also been controversies and different approaches toward dichotomizing aggregated urinary toxicity grades using a cutoff point [164]. While most studies use G2+ urinary toxicity as a meaningful complication grade cutoff, some use both G2+ and G1+ [46] and some either G3+ [197] or G1+ [145]. When the toxicity grades from different endpoints are averaged to get an overall toxicity grade, the ratio of subjects with Grade 3+ and even 2+ is usually very low compared to the control group. However, the distribution for each symptom varies, with some symptoms such as frequency showing many Grade 3+ subjects while hematuria shows very few. As discussed above, the relative importance of these symptoms is also quite different. This must also be taken into consideration when performing a dichotomized toxicity grade study. This study aimed to look into the toxicity scores of each urinary symptom in relation to the dosimetric and volumetric parameters determined in Chapter 3, as well as the bladder neck dose parameters. For all endpoints, except hematuria, two dichotomization regimens were used, one with G2+ and the other with G3+ cutoff. For hematuria only, G1+ and G2+ cutoffs were used, due to the nature of the complication and the observed distribution of toxicity grades. 5.2 Materials and Methods The LENT-SOMA subjective urinary toxicity data collected for a cohort of 60 cervical cancer 156 patients (same cohort as in Chapter 3) were used retrospectively. Patient selection was based on the average toxicity grades; all patients treated over the past few years with aggregated toxicity grades 2 and higher were included (total of 17) and the rest of the subjects (Grade 0-1) were randomly selected from the rest of the cases treated with HDR-BT at the BC Cancer Agency, Vancouver Centre, to make a sample size of 60. The toxicity grades in each LENT-SOMA questionnaire were assessed on a 0-4 discrete scale for seven urinary sequelae through the relevant measure for each symptom, including dysuria (severity), urgency (prevalence), hematuria (prevalence), frequency (time separation), nocturia (number of times), incontinence (prevalence), and reduced urine stream (prevalence). For the majority of the patients, there was more than one questionnaire collected over the course of treatment follow-up. There were six Control subjects with <3 months follow-up. Since the toxicity analysis results with and without these cases produced very similar results in the previous Chapter, these subjects were included in this analysis. The assigned toxicity grade for each of the seven urinary dysfunction endpoints was the “maximum” reported grade across all collected questionnaires. The distribution of the toxicity grades for every endpoint across all subjects was studied. For all the endpoints the patients were dichotomized in two different cutoff points based on the grade level. The cutoff points were selected differently for hematuria due to the nature and scarcity of the complication. For all endpoints one traditional G2-cutoff with Case: Grade 2+ and Control: Grade 0-1 groups and an additional grouping was used. For hematuria, a G1-cutoff grouping was used, i.e. Case: Grade 1+ and Control: Grade 0, and for all other endpoints, G3-157 cutoff grouping, i.e. Case: Grade 3+ and Control: Grade 0-2. The bladder neck location was identified in the planning CT for each fraction of each patient. This was performed by first locating the Foley balloon and the urethral catheter, which is clearly visible inside the bladder when it is filled with contrast material (Figure 5.1). Since the Foley balloon is positioned in the bladder so that it sits at the bladder neck, the axial planes were scrolled through to find the slice at the bottom of the Foley balloon. Then the sagittal and coronal orthogonal planes were rotated so that the Foley catheter was cut perpendicularly. The cursor was then positioned at the middle of the catheter when all three views were centered on the edge of the bladder outer contour. This position was recorded as the bladder neck, as shown in Figure 5.1. The bladder neck positions were audited and approved by a Radiologist at our center, for 10% of the cases. 158 Figure 5.1 Bladder neck location determined using the location of the Foley catheter on three orthogonal views and the bladder wall contour. Note that since the orthogonal views are rotated from the original planes, the bladder wall contour looks jagged. After the CPD point-set registration was used to accumulate the HDR BT fraction doses in a reference fraction, the bladder neck point in the reference fraction was used to calculate the contiguous volume of high-dose centered around this point. Three different contiguous volumes (0.1, 0.5 and 2cc) were evaluated. Bladder-wall dosimetric and volumetric parameters assessed in Chapter 4 were also included in this study. An analysis was performed to determine whether any of the bladder-wall and bladder-neck dose parameters were predictive of individual urinary sequelae. The normality of the Case and Control groups’ distributions were tested for both groupings and all parameters, using the Shapiro-Wilk normality test, which was performed using SPSS Statistics software (IBM Corp. Armonk, NY, version 23). For the parameters/endpoints that passed the normality test, the parametric Student t-test was used, with the null hypothesis that the means of Case and Control groups are equal. For those that did not pass the normality test, a non-parametric Mann-Whitney U-test was performed. The p-values calculated by these methods were used as a measure of the significance in the difference between Case and Control means. The statistical power (1-β) was also measured using G*Power [199] (Version 3.0.10), with the corresponding sample sizes, means, and standard deviations in the Case and Control groups applied for each endpoint and each parameter. The ureteral dose parameters were not studied as complication predictors both because the ureteral data were only available for 50 out of 60 cases and also because there was no promising trend observed in the aggregated toxicity grade analysis 159 for these parameters in Chapter 4. 5.3 Results The toxicity grade distribution for all seven urinary symptoms across the cohort is shown in Figure 5.2. The maximum number of Grade 4 and Grade 3 incidents was reported for the urgency and stream endpoints, respectively. The minimum number of Grade 4 and Grade 3 and the maximum number of Grade 0 toxicity levels were all observed for hematuria. On the other hand, only 5% of the subjects had Grade 0 frequency, whereas >50% reported no signs of incontinence, dysuria, hematuria and stream complications. 160 Figure 5.2 The toxicity grade distributions, separated by endpoints (top) and grade levels (bottom). Since the bladder neck was not identifiable in one of the subjects, the total number of subjects included in the bladder neck dose analysis was 59 rather than 60. Figures 5.3-5.5 illustrate the boxplot representation of the descriptive statistics for the registered dosimetric, volumetric, and bladder-neck-dose parameters for the Case and Control in G2-cutoff and G3/G1-cutoff groupings of incontinence, urgency, and hematuria endpoints. The results for the rest of the endpoints, which demonstrated non-significant outcomes, are included in Appendix A. Table 5.1 lists the p-value and s50tatistical power (1-β) results for the difference in the means 0%10%20%30%40%50%60%70%80%90%100%Percentage of CountsGrade 0 Grade 1 Grade 2 Grade 3 Grade 40%20%40%60%80%100%0 1 2 3 4Percentage of CountsToxicity GradeIncontinence UrgencyFrequency NocturiaDysuria HematuriaStream161 of the Case and Control groups, based on the G2 and G3/G1-cutoff groupings of incontinence, urgency, and hematuria endpoints, which showed the most [statistically] significant results, with some of the parameters in each grouping. Among all the urinary endpoints, incontinence showed the lowest p-values. The significant findings (p<0.05) for incontinence were observed for all the bladder-neck dose parameters in both groupings. The most significant result was achieved for BND0.1cc and BND0.5cc in the G3-cutoff grouping, with p<0.01 and the highest statistical powers achieved among all the parameters, endpoints, and groupings. The G2 grouping for incontinence showed p<0.1 for two of the bladder-wall dosimetric parameters and p<0.05 for V3Gy. The G3 grouping revealed p<0.1 for other four wall parameters. For urgency, the G2 grouping yielded the most significant results (p<0.05) for V3Gy, EQD22cc, and Idose. However, the statistical power achieved for these parameters was smaller than 0.8. The G1 grouping for hematuria did not yield any significant findings; however, G2 grouping for hematuria showed p<0.1 for four of the bladder-wall parameters. 162 Figure 5.3 Boxplot illustrating values of the different registered dosimetric and volumetric parameters (in Gy and cc, respectively) for Case and Control groups defined on the basis of Incontinence toxicity grades using G2 and G3-cutoffs. Red lines denote the median for each distribution, the top and bottom edges of the boxes indicate the 25th (𝒒𝟏) and 75th (𝒒𝟑) percentiles, the whiskers represent the minimum and maximum values excluding outliers, and the outliers (red crosses) are values > 𝒒𝟑 + 𝟏. 𝟓(𝒒𝟑 − 𝒒𝟏) or < 𝒒𝟏 − 𝟏. 𝟓(𝒒𝟑 − 𝒒𝟏). The integral dose values (I_dose) are divided by 100 in order to fit within the scale of the plot. 163 Figure 5.4 The boxplot for different registered dosimetric and volumetric parameters (in Gy and cc, respectively), separated by Case and Control groups based on the Urgency toxicity grades using G2 and G3-cutoffs. The meaning of the symbols is the same as in Fig. 5.2. The integral dose values (I_dose) are divided by 100 in order to fit within the scale of the plot. 164 Figure 5.5 The boxplot for different registered dosimetric and volumetric parameters (in Gy and cc, respectively), separated by Case and Control groups based on the Hematuria toxicity grades using G2 and G1-cutoffs. The meaning of the symbols is the same as in Fig. 5.2. The integral dose values (I_dose) are divided by 100 in order to fit within the scale of the plot.165 5.4 Discussion Seven urinary endpoints are defined according to the LENT-SOMA system. The subjective toxicity levels are collected and graded for each cervical cancer patient treated with multi-fraction HDR-ICBT at BCCA. For this study, all the high-toxicity cases reported over recent years and a selection of low-toxicity cases were used to form a cohort of 60 patients. The toxicity grade distributions varied significantly among different urinary endpoints. Urgency and stream showed the highest Grade 4 and Grade 3 complications, respectively. Incontinence, dysuria, hematuria, and stream had more Grade 0 cases compared to urgency, frequency, and nocturia. These observations suggest that some complications, particularly frequency, urgency and nocturia, might be more common in women due to other factors or just naturally due to, for example, a high level of liquid consumption. Therefore, it would have been helpful to have access to baseline (i.e. pre-treatment) data, which was lacking in this study, so that the pre-treatment urinary function could be accounted for in the analysis. Urgency and incontinence were the only symptoms that appeared to be significantly associated (p<0.05) with certain registered dose parameters. Hematuria showed some association based on p<0.1 significance level. No sign of association was observed for frequency, nocturia, dysuria, and stream. This observation might be explained, for frequency and nocturia, by limited follow-up data, lack of baseline data, and the likelihood that the baseline values for these complications might be more variable than they are for the other parameters. On the other hand, the lack of significant results for dysuria and stream might simply be due to lack of a significant dose effect at the levels to which this patient group was exposed. Larger sample sizes 166 Table 5.1 The p-values and statistical power (1-β) values based on Student t-test or Mann-Whitney U-test for the differences in the means of Case and Control groups for incontinence, urgency, and hematuria, based on G2 and G3/G1-cutoff groupings, for all the registered parameters. Parameter Urgency Incontinence Hematuria G2 G3 G2 G3 G2 G1 p 1-β p 1-β p 1-β p 1-β p 1-β p 1-β V3Gy 0.042 0.49 0.226 - 0.049* 0.52 0.290 - (0.188*) - (0.175) - V5Gy 0.087 0.35 0.204 - 0.067 0.50 0.109 - (0.468) - (0.247) - V10Gy 0.089 0.48 0.209 - 0.136 - 0.132 - 0.422* - (0.391) - D1cc 0.119 - 0.285 - 0.222 - 0.183 - 0.058* 0.46 0.432 - D2cc 0.068 0.39 0.153 - 0.150 - 0.119 - 0.079* 0.39 0.449 - D5cc 0.097 0.38 0.227 - 0.140 - 0.164 - 0.348 - (0.254) - D10cc 0.151 - 0.291 - 0.126 0.117* - 0.422* - (0.231) - EQD21cc 0.069 0.38 0.158 - 0.144 - 0.131 - 0.092* 0.40 (0.497) - EQD22cc 0.034 0.51 0.085 0.38 0.114 - 0.098 0.44 0.077 0.42 0.466 - EQD25cc 0.070 0.44 0.185 - 0.111 - 0.075* - 0.421 - (0.191) - EQD210cc 0.137 - 0.260 - 0.120* - 0.098* - 0.422* - (0.214) - Idose 0.026 0.57 0.122 - 0.057 0.55 0.125 - 0.471* - (0.356) - Mdose 0.253 - 0.285 - 0.139 - 0.057 0.57 0.487 - (0.168) - BND0.1cc 0.205 - 0.122* 0.47 0.031* 0.73 0.001* 0.98 0.343 - 0.482 - BND0.5cc 0.430* - 0.196* 0.43 0.015* 0.79 0.001* 0.98 0.338 - (0.497) - BND2cc 0.313* - 0.330* - 0.036* 0.69 0.003* 0.94 0.399 - (0.322) - Notes: The p-values smaller than 0.1 are highlighted in bold and those smaller than 0.05 are underlined. The asterisks mark the p-values based on the non-parametric U-test (performed when at least one of the Case and Control distributions was not normal) and the rest are based on t-test. For parameters where the mean of the Case group is larger than that of the Control group, the p-values are shown in parentheses. The statistical power values are only shown for results with p<0.1. The power values larger than 0.8 are shown in bold. and baseline information are required to validate these theories. Overall, the G2 cutoff demonstrated the most p<0.1 and p<0.05 results (19 and 7, respectively) compared to G3 cutoff (8 and 3, respectively). This observation is in agreement with 167 the choice of G2 cutoff used in most other urinary toxicity studies and also the cutoff showing most positive findings for urinary disorders [46, 164]. However, the most significant results of all (p<0.01 and 1-β>0.9) were observed in the G3 grouping for bladder neck dose parameters as predictors of incontinence. Therefore, based on the results of this study it can be recommended to use both G2 and G3 cutoffs for similar dose-effect studies, but if a single threshold is used, the traditional G2 cutoff seems to be a more suitable choice. All parameters except D10cc, showed at least one p<0.1 result in one of the endpoint groupings. The only bladder-wall parameters showing significant association with a urinary endpoint were V3Gy, EQD22cc, and Idose, in the G2 grouping. The most significant results across all parameters were observed for incontinence for the G3 cut-off with both BND0.1cc and BND0.5cc as predictors (p=0.001). The statistical power was also 0.98 for both parameters. Significant results (p<0.05) for incontinence were also obtained for all the bladder neck dose parameters in the G2 grouping. This finding strongly suggests that dose to small volumes in the region of the bladder neck predicts for incontinence, which supports our initial hypothesis. However, BND2cc also appeared to be a strong predictor for incontinence. Considering that the volume of the bladder centered around the bladder neck would represent part of the “trigone” structure, these findings suggest that dose received by trigone in bladder may be a predictor for urinary incontinence. Not finding any significant relationship between dose to bladder neck and the rest of urinary endpoints was also as anticipated. There were some signs of an association between hematuria and accumulated dose to small volumes (i.e. 1-2 cc) as hypothesized, but the corresponding p-values were not significant and further investigation for this effect is required, preferably with larger sample size, longer follow-up, and with baseline (pre-treatment) scores available. The G1 grouping for hematuria did not 168 seem to uncover any important predictors and, therefore, is not recommended based on this study. Urgency showed the most significant association with Idose, which might suggest that the dose to large volumes of bladder wall could be a better predictor for the urgency endpoint, as might be expected. However, on the same basis, one would expect the frequency endpoint to provide the same result. In either case, lack of baseline information and small sample size precludes drawing any conclusions at this point. Moreover, there was only a small group of high-toxicity cases, which affects the presence, significance and statistical power of the observations. Therefore, a larger sample size with more high-toxicity cases would help clarify these findings. It must also be acknowledged that the number of dose-volume parameters included in this study is large for the size of the data set, so the findings may be subject to the multiple comparisons problem [186], causing a type I error (i.e. false positive for null hypothesis rejection). However, some of the parameters are highly correlated and therefore cannot be considered completely independent – e.g. R = 0.99 for BND0.5cc versus BND0.1cc, 0.92 for D2cc versus D1cc, and 0.95 for EQD22cc versus EQD22cc, where R is the correlation coefficient). If considered independent, application of the Holm's step-down approach [187] for controlling the family wise error rate (FWER) [25] determined that true rejection of null hypothesis for bladder neck parameters and incontinence was established based on the G3 grouping. 5.5 Conclusion This study looked into the association between locally accumulated dose parameters and urinary endpoints. It was found that dose to small volumes of 0.1-2cc around the bladder neck significantly predict for incontinence. On the other hand, there was no significant association 169 between bladder neck dose and other endpoints, which was not unexpected. Urgency was the most associated with integral dose in the bladder wall, as was hematuria with doses to small wall volumes; however, the statistical power was smaller than the accepted threshold of 0.8. Therefore, while these findings for urgency and hematuria are of interest and promising, more investigation and larger sample sizes are required before drawing conclusions. Lack of baseline data and insufficient follow-up for a few subjects in this study may have prevented detection of a dose response for certain endpoints such as frequency, which could have had high baseline values in this population. Also, lack of sufficient high-toxicity cases for some of the endpoints, such as nocturia, dysuria, stream and hematuria, while clinically desirable, makes the dose-response of these endpoints difficult to study. Validation of the results of this study with a larger cohort, ideally with baseline toxicity grades and longer follow-up available, is therefore recommended. 170 Chapter 6: Conclusion and Future Work 6.1 Conclusion The main objectives of this dissertation were to (i) develop a method to track bladder wall movement between multi-fraction HDR ICBT treatments for cervical cancer patients, (ii) accumulate the HDR fraction doses in the bladder tissue, and (iii) examine associations between patient-reported LENT-SOMA subjective urinary toxicity grades and the locally-accumulated bladder-wall dose. Both deformable image registration and point-set registration were evaluated for objective (i) and different methods and algorithms were examined. As a definitive step in this process, a porcine bladder phantom with external fiducial markers was used to test the performance of a selection of deformable registration methods applied to the specific task of registering the bladder contours with varying sizes due to the amount of liquid filling. In the study described in Chapter 3, the CPD deformable point-set registration toolbox showed superior registration accuracy results for bladder contour registration, and was therefore chosen for the analysis of patient data. The average registration uncertainty for CPD was found to be 6.4 mm. While there is room for improvement, this registration uncertainty was felt to be acceptable for dose-volume regions of interest greater than 1 cm3, given the small thickness of the bladder wall (which leads to the volume of interest having a large corresponding surface area), and the fact that the uncertainty is primarily in a direction tangential to the surface. The feasibility of using anatomical landmark TRE values for the patient data as a measure of the registration accuracy was established through comparison of the bladder neck TRE with the average fiducial TRE values in the porcine phantom, providing a means of assessing the 171 performance of the algorithm in real patient data. In the study in Chapter 4 with the patient data, the average landmark TRE values over 59 cases was 5.5±2.0 mm, which was very close to the bladder neck landmark TRE value achieved by CPD for the phantom data (5.4 mm) in this chapter. Also, this value was smaller than CPD’s average fiducial TRE value on phantom data (6.4±3.5). Therefore, contrary to expectations, CPD’s registration accuracy for the patient data appeared to be as good as the phantom results, at least in the posterior portion of the bladder where the anatomical landmarks lie. The planning CT scans of the HDR ICBT treatment fractions for 60 cervical cancer patients treated at BCCA were used, together with the collected LENT-SOMA urinary toxicity data for these patients. The bladder wall contours were created manually and then registered using CPD. The fraction doses were accumulated on the bladder wall and different registered dosimetric and volumetric parameters were calculated for the bladder wall, in addition to dosimetric parameters for small volumes around the ureteral and urethral opening points. The results of comparing DVH-summed and registered dose parameters supported the original hypothesis that the deformable dose accumulation is expected to be significantly different from DVH-based dose accumulation. A novel in-house region-growing algorithm was developed in the study in Chapter 4 to calculate contiguous volumes of high dose. Defining biologically meaningful volumes is a challenging task, especially due to the ambiguity associated with defining contiguous high-dose volumes when multiple hot-spot volumes begin to merge, as was shown in Figure 4.4. Based on the results of this study, DVH-summed and locally-accumulated dose parameters were found to be different and the differences were even larger when doses to contiguous hot spots were considered. The differences between the contiguous registered dose and the DVH-summed 172 parameters appear to be largely attributable to the inherent difference in these two parameters when hotspots are distributed over the volume of interest, since the contiguous volume analysis captures some aspects of the spatial distribution while conventional DVH-based analysis does not. The traditional DVH curves are based on non-contiguous volumes and can therefore, misrepresent a biologically meaningful volume of tissue. This may also be of relevance in dosimetric analysis of single treatment fractions, since it is not uncommon for the bladder to have two separate hotspots due to the location of the HDR applicator (and its corresponding dose distribution) with respect to the so-called ‘lateral horns’ [184] of the bladder. However, routine adoption of contiguous volume analysis would have to be undertaken with some caution due to the ambiguity associated with defining contiguous high-dose volumes, as mentioned above. Consensus on an appropriate level of spatial regularization would have to be achieved if contiguous volumes results were to be compared across multiple institutions. The difference between contiguous and non-contiguous DVH curves in the presence of multiple high-dose regions can be further illustrated through an example of a synthetic case mimicking a typical summed bladder dose distribution for cervical cancer brachytherapy with two separate hotspots (Figure 6.1). This figure was generated using VariSeed™ (Varian’s low-dose-rate prostate brachytherapy treatment planning system), which performs contiguous dose-volume analysis. Unfortunately, this software is designed for planning seed implants only, so it can only be used to “simulate” cervical brachytherapy dose distributions, but it offers an independent means of displaying and calculating contiguous dose volumes. Figure 6.2 shows the dose-volume histogram results generated by VariSeed for the 173 distribution of Figure 6.1. It can be seen in Figure 6.2 that at the dose threshold where the two hotspots connect, the “contiguous” DVH curves join the traditional “non-contiguous” DVH in a “step function”. At all dose values above this threshold, the corresponding non-contiguous volume is the sum of the two contiguous ones, and the contiguous Dncc values are clearly much lower than the corresponding non-contiguous one. For example, the D1cc value for the smaller and larger contiguous volumes are about 20% and 12%, respectively, smaller than the corresponding non-contiguous D1cc value. These values agree with our overall findings in the registered datasets. From these curves, it is also obvious that a volume such as 3cc would require contributions from both of the “hot spot” regions, and that there is more than one way to combine them. Thus, the definition of a contiguous volume in this scenario is open to a degree of subjectivity. Ideally, one would endeavor to choose a volume definition that is the most biologically meaningful. Refining and standardizing the definition of the contiguous volume is a potential direction for future studies following the work reported in this dissertation. The association between the assessed parameters and the urinary toxicity was studied through two separate analyses. In Chapter 4, the average toxicity scores were used to dichotomize the cases into Control and Case groups. For the locally-accumulated dose parameters, the only significant (p<0.05) difference between Case and Control groups was observed for V3Gy. In Chapter 5, the dose parameters were studied as predictors of individual urinary endpoints, using two different dichotomization regimens. A highly significant association was found between dose to small volumes around the bladder neck and incontinence, although no relationship was found between bladder neck dose and total toxicity score. Also, it was found that accumulated integral dose to the bladder wall was the most predictive of urgency, while equivalent dose to smaller volumes of blad- 174 Figure 6.1 A mock bladder dose distribution created in VariSeed, showing two separate high-dose regions. The larger maroon region represents the contiguous D2cc, and all doses in the maroon regions are > D2cc. der wall was predictive of hematuria based on 0.1 significance level. However, because of the multiple comparison issue, most of these findings are hypothesis-generating only. While the drawback of aggregating different morbidity scores has been argued by Rosewall et al. [164] and Bentzen et al. [192], no comprehensive studies to investigate and validate this argument have been available in the literature to date. The current study provided a thorough investigation, supporting the fact that aggregating individual morbidity information into a single statistic conceals the possible links between normal tissue complication and radiation dose. Lack of baseline data for urinary sequelae was an important limitation of the toxicity-dose relation analysis in this study. Some of the symptoms, such as frequency, may have had relatively large background values, which might have confounded the dose-response analysis for these ques- 175 Figure 6.2 Comparison of the DVH curves for the larger and smaller contiguous volumes with the traditional non-contiguous DVH curve. The bottom graph is a close-up of the top graph. tionnaire available. Typically, patients were lost to follow-up either because they did not come for their follow-up appointments, or because completion of the LENT SOMA questionnaires was aborted at disease progression. While using length of follow-up as an exclusion criterion in 0204060801001201401600 20 40 60 80 100 120Volume (cc)Dose (%)larger contiguous volumeSmaller contiguous volumenon-contiguous DVH01234567891080 85 90 95 100 105 110 115 120 125 130Volume (cc)Dose (%)larger contiguous volumesmaller contiguous volumenon-contiguous DVH176 selecting the sample population was an option, this could have introduced bias as well. For example, patients who are less well, but not so poorly that they would be classified into Case group, may be more likely to return for follow-up appointments than patients who are doing very well and not experiencing significant side-effects from their treatment. In this study, the effectiveness of varying the bladder filling volume across different HDR BT treatment fractions for the purpose of moving the location of the hotspots and reducing urinary toxicity was also evaluated. It was observed that varying the bladder filling volume correlated significantly (p=0.03) with both hotspot separation and urinary toxicity. This result suggests that varying the bladder filling regimen has the potential to reduce the urinary toxicity by moving the location of the maximally irradiated regions on the bladder across different brachytherapy treatment sessions. In conclusion, this study presents a novel method to locally accumulate fractionated dose using data from multiple CT image sets acquired at or near the time of dose delivery. A comprehensive analysis of urinary toxicity and dose relationship in this study suggests that studying the individual urinary symptom levels provides a better understanding of the dose relationship for urinary dysfunction after HDR ICBT. Also, as hypothesized, cumulative dose calculated using deformably registered data sets appeared to be superior to DVH-based summed dose in predicting urinary toxicity in our data set. However, further investigation is required to arrive at a definitive conclusion. While some promising and interesting results were achieved here, there were aspects of this research that could benefit from future work, as discussed below. 177 6.2 Future Work The study presented here demonstrated a novel method to locally accumulate fractionated dose in a highly deformable OAR. This method was used to study morbidity in relation to accumulated dose in the bladder, and can potentially also be used for other highly deformable OARs in the multi-fraction HDR ICBT treatment, such as the rectum, which is also subject to radiation side effects [169], although additional validation studies would be required. Also, the dose accumulation method carried out here can be used for other fractionated treatment regimens with multiple images, such as interstitial brachytherapy [160] or EBRT treatments with repeated imaging. Furthermore, the capabilities of deformable point-set registration can potentially be expanded to other radiotherapy applications, such as multi-modality registrations of ROIs. An investigation aimed at refining and standardizing the definition of a biologically meaningful contiguous volume is a potential direction for future work following this study. These investigations may benefit from multi-centre studies to achieve a consensus on an appropriate level of spatial regularization by comparing contiguous volumes results across multiple institutions. Continuation of the current research could include incorporation of additional potential co-morbidity factors, such as demographic background, body mass index (BMI), and baseline functional scores into the analysis. Also, the timeline of ICBT with respect to the concomitant EBRT treatment may be studied as an additional factor. A multi-variate analysis between these factors and the accumulated dose parameters could then be performed. Further improvement to the data set could also be sought through longer and more consistent follow-up of the patients, which would improve the reliability of the toxicity information. 178 However, contacting and inviting the patients to check in and fill questionnaires is challenging and not always ethically allowed. A larger sample size is also advised in order to improve the statistical strength of succeeding studies. 179 Bibliography [2] Jiaquan Xu, S. L. Murphy, K. D. Kochanek and B. A. Bastian, "Deaths: Final Data for 2013," National vital statistics reports, vol. 64, no. 2, 2010. [2] "Statistics, Canadian Cancer Society’s Advisory Committee on Cancer. Canadian Cancer Statistics 2015," Canadian Cancer Society, Toronto, 2015. [3] J. Ferlay, I. Soerjomataram, M. Ervik and e. al., Cancer Incidence and Mortality Worldwide, Lyon, France: International Agency for Research on Cancer, 2013. [4] S. E. Waggoner, "Cervical cancer," LANCET, vol. 361, no. 9376, pp. 2217-2225, 2003. [5] Canadian Cancer Society, "Risk factors for cervical cancer," 2016. [Online]. Available: http://www.cancer.ca/en/cancer-information/cancer-type/cervical/risks/?region=bc. [Accessed 24 April 2016]. [6] M. E. Powell, "Modern Radiotherapy and Cervical Cancer," Int J Gynecol Cancer, vol. 20, no. S2, pp. S49-S51, 2010. [7] R. Phillips, Super Voltage x-ray Therapy, London, UK: HK Lewis, 1944. [8] S. Goldberg and E. London, "Zur Frage der Beziehungen zwischen Becquerelstrahlen und Hautaffectionen," Dermatologische Zeit, vol. 10, p. 457, 1903. [9] B. Hilaris, D. Mastoras and L. Shih, "History of Brachytherapy: The Years After the 180 Discovery of Radium and Radioactivity," in Principles and Practice of Brachytherapy, N. Subir, Ed., Armonk, NY, Futura, 1997. [10] P. Finger, "Radiation Therapy for Orbital Tumors: Concepts, Current Use, and Ophthalmic Radiation Side Effects," Survey of Opthalmology, vol. 54, no. 5, pp. 545-568, 2009. [11] ICRU Report 33, Radiation quantities and units, Washington, DC: International Comission on Radiation Units ad Measurements, 1980. [12] A. Flynn, "Isotopes and delivery systems for brachytherapy," in Radiotherapy in practice: brachytherapy, C. C. Hoskin P., Ed., New York, Oxford University Press, 2005. [13] V. Gupta, "Brachytherapy – past, present and future," Journal of Medical Physics,, vol. 20, pp. 31-38, 1995. [14] S. Nag, "High dose rate brachytherapy: its clinical applications and treatment guidelines," Technology in Cancer Research and Treatment, vol. 3, no. 3, p. 269–87, 2004. [15] R. Nath, L. L. Anderson, G. Luxton, K. A. Weaver and J. F. Williamson, "Dosimetry of Interstitial Brachytherapy Sources: Recommendations of the AAPM Radiation Therapy Committee Task Group no. 43. American Association of Physicists in Medicine," Med Phys, vol. 22, pp. 209-34, 1995. [16] M. J. Rivard, B. M. Coursey, L. A. DeWerd, W. F. Hanson, M. S. Huq and G. S. Ibbott, "Update of AAPM Task Group no. 43 Report: A Revised AAPM Protocol for Brachytherapy Dose Calculations," Med Phys, vol. 31, pp. 633-74, 2004. 181 [17] J. F. Williamson, H. Perera, Z. Li and W. R. Lutz, "Comparison of calculated and measured heterogeneity correction factors for I-125, Cs-137, and Ir-192 brachytherapy sources near localized heterogeneities," Med. Phys., vol. 20, p. 209–222, 1993. [18] J. F. Williamson, R. Baker and Z. Li, "A convolution algorithm for brachytherapy dose computations in heterogeneous geometries," Med. Phys., vol. 18, p. 1256–1265, 1991. [19] E. Hall and A. Giaccia, Radiobiologty for radiologist, sixth ed., Philadelphia, PA: Lippincott Williams & Wilkins, 2006. [20] K. Chadwick and H. Leenhouts, "A molecular theory of cell survival," Phys Med Biol, vol. 18, p. 78, 1973. [21] J. Chapman, C. Gillespie, A. Reuvers and D. DL, "The inactivation of chinese hamster cells by x-rays: the effects of chemical modifiers on single- and double-events," Radiat Res, vol. 64, p. 365, 1975. [22] J. Chapman, "Biological models of mammalian cell inactivation by radiation," in Radiation Biology in Cancer Research, New York, Raven Press, 1980, p. 21. [23] J. Fowler, "The linear-quadratic formula and progress in fractionated radiotherapy," Br J Radiol, vol. 62, no. 740, pp. 679-694, 1989. [24] C. Beskow, A. Ågren-Cronqvist, R. Lewensohn and I. Toma-Dasu, "Biological effective dose evaluation and assessment of rectal and bladder complications for cervical cancer treated with radiotherapy and surgery," J Contemp Brachytherapy, vol. 4, no. 4, p. 205–182 212, 2012. [25] R. Dale and B. Jones, "The clinical radiobiology of brachytherapy," Br. J. Radiol., vol. 71, no. 845, pp. 465-483, 1998. [26] B. Emami, J. Lyman and e. a. Brown A, "Tolerance of normal tissue to therapeutic irradiation," International Journal of Radiation Oncology, Biology, Physics, vol. 21, pp. 109-122, 1991. [27] S. Bentzen, L. Constine and J. Deasy, "Quantitative Analyses of Normal Tissue Effects in the Clinic (QUANTEC): An Introduction to the Scientific Issues," International Journal of Radiation Oncology, Biology, Physics, vol. 76, pp. S3-S9, 2010. [28] "Lent soma scales for all anatomic sites," International Journal of Radiation Oncology, Biology, Physics., vol. 31, pp. 1049-1091, 1995. [29] P. Rubin, L. Constine, L. Fajardo, T. Phillips and T. Wasserman, "Overview: Late effects of normal tissues (LENT) scoring system," International Journal of Radiation Oncology, Biology, Physics, vol. 31, pp. 1041-1042, 1995. [30] J. Pavy, J. Denekamp and J. Letschert, "Late effects toxicity scoring: the soma scale," International Journal of Radiation Oncology, Biology, Physics, vol. 31, pp. 1043-1047, 1995. [31] U.S. DEPARTMENT OF HEALTH AND HUMAN SERVICES, "Common Terminology Criteria for Adverse Events (CTCAE)," 14 June 2010. [Online]. Available: 183 http://ctep.cancer.gov/protocolDevelopment/electronic_applications/ctc.htm#ctc_40. [Accessed 20 April 2016]. [32] S. Davidson, M. Burns and J. Routledge, "Short Report: A Morbidity Scoring System for Clinical Oncology Practice: Questionnaires produced from the LENT SOMA scoring system," Clinical Oncology, vol. 14, pp. 68-69, 2002. [33] The Christie Hospital, "CTCAE/LENT SOMA questionnaires," [Online]. Available: http://www.christie.nhs.uk/about-us/our-standards/clinical-outcomes-unit/clinical-oncology-scoring-treatment-effects/ctcaelent-soma-questionnaires/. [Accessed 21 April 2016]. [34] ICRU Report 50, Prescribing, recording and reporting photon beam therapy, Bethesda, MD: International Commission on Radiation Units and Measurements, 1993. [35] ICRU Report 62, Prescribing, recording and reporting photon beam therapy (supplement to ICRU report 50), Bethesda, MD: International Commission on Radiation Units and Measurements, 1999. [36] G. Montana, W. Fowler, M. Varia, L. Walton, Y. Mack and L. Shemanski, "Carcinoma of the cervix, stage III. Results of radiation therapy," Cancer, vol. 57, pp. 148-154, 1986. [37] R. Lanciano, M. Won, L. Coia and G. Hanks, "Pretreatment and treatment factors associated with improved outcome in squamous cell carcinoma of the uterine cervix: A final report of the 1973 and 1978 patterns of care studies," Int J of Radiat Oncol, Biol Phys, vol. 20, pp. 184 667-676, 1991. [38] C. Perez, S. Breaux and H. Madoc-Jones, "Radiation therapy alone in the treatment of carcinoma of the uterine cervix. I. Analysis of tumor recurrence," Cancer, vol. 51, p. 1393–1402, 1983. [39] D. O’Connell, N. Howard and C. Joslin, "A new remotely controlled unit for the treatment of uterine carcinoma," Lancet, vol. 2, p. 570 –571, 1965. [40] M. Wakabayashi, G. Irie and T. Sugawara, "The trial production of remote after-loading system unit. Japan.," Jpn J Clin Radiol, vol. 11, p. 678–648., 1966. [41] S. Nag, B. Erickson, B. Thomadsen, C. Orton, J. Demanes and D. Petereit, "The American Brachytherapy Society recommendations for high-dose-rate brachytherapy for carcinoma of the cervix," Int J Radiat Oncol Biol Phys, vol. 48, pp. 201-211, 2000. [42] C. Haie-Meder, R. Pötter, E. Van Limbergen and e. al., "Recommendations from Gynaecological (GYN) GEC-ESTRO Working Group ☆ (I): concepts and terms in 3D image based 3D treatment planning in cervix cancer brachytherapy with emphasis on MRI assessment of GTV and CTV," Radiotherapy and Oncology, vol. 74, pp. 235-245, 2005. [43] R. Pötter, C. Haie-Meder and E. Limbergen, "Recommendations from gynaecological (GYN) GEC ESTRO working group (II): Concepts and terms in 3D image-based treatment planning in cervix cancer brachytherapy—3D dose volume parameters and aspects of 3D image-based anatomy, radiation physics, radiobiology," Radiotherapy and Oncology., vol. 185 78, pp. 67-77, 2006. [44] ICRU Report 38, Dose and volume specification for reporting intracavitary, Bethesda, MD: International Commission on Radiation Units and Measurements, 1985. [45] R. Pötter, E. Van Limbergen, N. Gerstner and A. Wambersie, "Survey of the use of the ICRU 38 in recording and reporting cervical cancer brachytherapy," Radiotherapy and Oncology, vol. 58, pp. 11-18, 2001. [46] P. Georg, S. Lang, J. Dimopoulos and W. Dörr, "Dose–Volume Histogram Parameters and Late Side Effects in Magnetic Resonance Image–Guided Adaptive Cervical Cancer Brachytherapy," Int. J. Radiat. Oncol. Biol., vol. 79, no. 2, pp. 356-362, 2011. [47] D. L. G. Hill, P. G. Batchelor, M. Holden and D. J. Hawkes, "Medical image registration," Phys. Med. Biol., vol. 46 , p. R1–R45, 2001. [48] D. Rueckert, "Nonrigid Registration: Concepts, Algorithms, and Applications," in Medical Image Registration, Boca Raton, FL, CRC Press LLC, 2001, pp. 281-301. [49] W. R. Crum, T. Hartkens and D. L. G. Hill, "Non-rigid image registration: theory and practice," The British Journal of Radiology, vol. 77, p. S140–S153, 2004. [50] D. Rueckert and P. Aljabar, "Nonrigid Registration of Medical Images: Theory, Methods, and Applications," IEEE Signal Processing Magazine, vol. 27, pp. 113-119, 2010. [51] J. Gee, T. Sundaram, I. Hasegawa, H. Uematsu and H. Hatabu, "Characterization of regional pulmonary mechanics from serial magnetic resonance imaging data," Acad Radiol, 186 vol. 10, p. 1147–52, 2003. [52] G. Christensen, B. Carlson, K. Chao, P. Yin and P. Grigsby, "Image-based dose planning of intracavitary brachytherapy: registration of serial imaging studies using deformable anatomic templates," Int J Radiat Oncol Biol Phys, vol. 51, p. 227–43, 2001. [53] T. Rohlfing, C. Maurer, W. O’Dell and J. Zhong, "Modeling liver motion and deformation during the respiratory cycle using intensity-based nonrigid registration of gated MR images," Med Phys, vol. 31, p. 427–32, 2004. [54] P. Thompson, R. Woods, M. Mega and A. Toga, "Mathematical/computational challenges in creating deformable and probabilistic atlases of the human brain," Human Brain Mapping, vol. 9, p. 81–92, 2000. [55] B. Li, G. Christensen, E. Hoffman and G. McLennan, "Establishing a normative atlas of the human lung: intersubject warping and registration of volumetric CT images," Acad Radiol, vol. 10, p. 255–65, 2003. [56] M. Ferrant, A. Nabavi, B. Macq, P. Black and F. Jolesz, "Serial registration of intraoperative MR images of the brain," Med Image Anal, vol. 6, p. 337–59, 2002. [57] A. Sotiras, C. Davatzikos and N. Paragios, "Deformable Medical Image Registration: A Survey," IEEE Transactions on Medical Imaging, vol. 32, pp. 1153-1190, 2013. [58] B. LG, "A survey of image registration techniques," Computing Surveys, vol. 24, p. 325–76, 1992. 187 [59] J. Modersitzki, Numerical Methods for Image Registration, New York: Oxford Univ. Press, 2004. [60] J. Modersitzki, FAIR: Flexible Algorithms for Image Registration, Philadelphia, PA: SIAM, 2009. [61] J. Shackleford, N. Kandasamy and G. Sharp, High-Performance Deformable Image Registration Algorithms for Manycore Processors, Waltham, MA: Elsevier/MK, 2013. [62] R. Bajcsy and S. Kovacˇicˇ, "Multiresolution elastic matching," Comp. Vision, Graphics and Image Process., vol. 46, p. 1–21, 1989. [63] G. Christensen, R. Rabbitt and M. Miller, "Deformable templates using large deformation kinematics," IEEE Trans. Image Process., vol. 5, no. 10, p. 1435–1447, 1996. [64] W. Press, B. Flannery, S. Teukolsky and W. T. Vetterling, Numerical Recipes in C, 2nd ed. ed., Cambridge University Press, 1989. [65] M. Bro-Nielsen and C. Cramkow, "Fast Fluid Registration of Medical Images," Visualization in Biomedical Computing, Lecture Notes in Computer Science, vol. 1131, p. 267–76, 1996. [66] R. Szeliski and S. Lavallée, "Matching 3-D anatomical surfaces with non-rigid deformations using octree-splines," IEEE Workshop on Biomedical Image Analysis, p. 144–153, 1994. [67] J. Ashburner and K. J. Friston, "Nonlinear spatial normalization using basis functions," 188 Human Brain Mapping, vol. 7, p. 254–66, 1999. [68] Y. Amit, "A nonlinear variational problem for image matching," SIAM J. Scien Comp., vol. 15, no. 1, p. 207–224, 1994. [69] M. Davis, A. Khotanzad, D. Flaming and S. Harms, "A physics-based coordinate transformation for 3-D image matching," IEEE Trans. Med. Imaging, vol. 16, no. 3, p. 317–328, 1997. [70] B. FL, "Principal warps – thin-plate splines and the decomposition of deformations," IEEE Trans Pattern Anal, vol. 11, p. 567–85, 1989. [71] D. Mattes, D. Haynor, H. Vesselle and T. Lewellen, "PET-CT image registration in the chest using free-form deformations," IEEE Trans Med Imaging, vol. 22, p. 120–8, 2003. [72] T. Rohlfing, C. J. Maurer, D. Bluemke and M. Jacobs, "Volume-preserving nonrigid registration of MR breast images using free-form deformation with an incompressibility constraint," IEEE Trans Med Imaging, vol. 22, p. 730–41, 2003. [73] D. Rueckert, L. Sonoda, C. Hayes, D. Hill and M. Leach, "Non-rigid registration using free-form deformations: application to breast MR images," IEEE Trans. Med. Imaging, vol. 18, p. 712–21, 1999 . [74] C. Maurer, J. Fitzpatrick, M. Wang, R. Galloway and R. Maciunas, "Registration of head volume images using implantable fiducial markers," IEEE, Transactions on Medical Imaging, vol. 16, pp. 447-462, 1997. 189 [75] L. Schad, R. Boesecke, W. Schlegel, G. Hartmann, G. Sturm and L. Strauss, "Three dimensional image correlation of CT, MR, and PET studies in radiotherapy treatment of brain tumors," Journal of Computer Assisted Tomography, vol. 11, pp. 948-954, 1987. [76] J. B. West, J. M. Fitzpatrick, M. Y. Wang, B. M. Dawant, C. R. Maurer and R. M. Kessler, "Comparison and evaluation of retrospective intermodality image registration techniques," J. Comput. Assist Tomogr., vol. 21, no. 4, p. 554–566, 1997. [77] T. Lange, N. Papenberg and S. Heldmann, "3D ultrasound-CT registration of the liver using combined landmark-intensity information," International Journal of Computer Assisted Radiology and Surgery, vol. 4, pp. 79-88, 2009. [78] C. Pelizzari, G. Chen, D. Spelbring, R. Weichselbaum and C. Chen, "Accurate three-dimensional registration of CT, PET, and or MR images of the," Computer assisted tomography, vol. 13, pp. 20-26, 1989. [79] J. Herring, B. Dawant, C. J. Maurer, D. Muratore and R. Galloway, "Surface-based registration of CT images to physical space for image-guided surgery of the spine: a sensitivity study," IEEE, Transactions on Medical Imaging, vol. 17, p. 743 –752, 1998. [80] C. J. Maurer, R. Maciunas and J. Fitzpatrick, "Registration of head CT images to physical space using a weighted combination of points and surfaces," IEEE, Transactions on Medical Imaging, vol. 17, p. 753 –761, 1998. [81] D. Pham, C. Xu and J. Prince, "Current methods in medical image segmentation," Annual 190 review of biomedical engineering, vol. 2, pp. 315-37, 2000. [82] R. Adams, "Seeded region growing," IEEE transactions on pattern analysis and machine intelligence, vol. 16, pp. 641-647, 1994 . [83] M. Sezgin and B. Sankur, " Survey over image thresholding techniques and quantitative performance evaluation," Journal of Electronic Imaging, vol. 13, pp. 146-168, 2004. [84] G. Borgefors, "Distance transformations in digital images," Comput. Vision, Graphics Image Process, vol. 34, p. 344–371, 1986. [85] H. Jiang, R. Robb and K. Holton, "New approach to 3-D registration of multimodality medical images by surface matching," SPIE, vol. 1808, p. 196–213, 1992. [86] P. Besl and N. McKay, "A method for registration of 3D shapes.," IEEE Trans. Pattern Anal. Machine Vision, vol. 14, p. 239–256, 1992. [87] J. Hajnal, N. Saeed, E. Soar, A. Oatridge, I. Young and G. Bydder, "A Registration and Interpolation Procedure for Subvoxel Matching of Serially Acquired MR Images," Journal of Computer Assisted Tomography, vol. 19, pp. 289-296,, 1995. [88] C. Studholme, D. Hill and D. Hawkes, "Automated 3D MR and PET brain image registration," Journal of Computer Assisted Radiology, vol. 1, pp. 248-253, 1995. [89] W. M. Wells and P. Viola, "Alignment by maximization of mutual information," Int. J. Comput. Vis., vol. 24, no. 2, pp. 137-154, 1997. 191 [90] F. Maes, A. Collignon, D. Vandermeulen, G.Marechal and R. Suetens, "Multimodality image registration by maximization of mutual information," IEEE Trans. Med. Imag., vol. 16, no. 2, pp. 187-198, 1997. [91] H. Chui and R. A., "A new point matching algorithm for non-rigid registration," Comput Vision Image Understanding, vol. 89, p. 114–41, 2003. [92] A. Fitzgibbon, "Robust registration of 2D and 3D point sets," Image Vision Comput, vol. 21, p. 1145–53, 2003. [93] O. van Kaick, H. R. Zhang, G. Hamarneh and D. Cohen-Or, "A survey on shape correspondence," Comput. Graph. Forum., vol. 30, p. 1681–707, 2011. [94] "CPD_fish_nonrigid.gif," By Dllu (Own work) [CC BY-SA 3.0 (http://creativecommons.org/licenses/by-sa/3.0)], via Wikimedia Commons, [Online]. Available: https://upload.wikimedia.org/wikipedia/commons/b/b2/Cpd_fish_nonrigid.gif. [95] E. Bender and W. Tome, "The utilization of consistency metrics for error analysis in deformable image registration," Phys Med Biol, vol. 54, no. 18, p. 5561–67, 2009. [96] D. Yang, H. Li, D. Low, J. Deasy and I. Naqa, "A fast inverse consistent deformable image registration method based on symmetric optical flow computation," Phys Med Biol, vol. 53, no. 21, p. 6143–65, 2008. [97] J. A. Schnabel, C. Tanner, A. D. Castellano-Smith, A. Degenhard and M. O. Leach, "Validation of non-rigid image registration using finite element methods Application to 192 breast MR images," IEEE Trans. Med Imag., vol. 22, no. 2, p. 238–247, 2003. [98] A. Klein, J. Andersson, B. A. Ardekani, J. Ashburner and B. Avants, "Evaluation of 14 nonlinear deformation algorithms applied to human brain MRI registration," Neuroimage, vol. 46, no. 3, p. 786–802, 2009. [99] M. J. Fitzpatrick, "Detecting failure, assessing success," in Medical Image Registration, J. V. Hajnal and D. L. Hill, Eds., Boca Raton, FL, CRC, 2001, p. 117–36. [100] L. Dice, "Measures of the amount of ecologic association between species," Ecology, pp. 297-302, 1945. [101] D. Huttenlocher, G. Klanderman and W. Rucklidge, "Comparing images using the Hausdorff distance," IEEE Trans Pattern Anal Mach Intell, vol. 15, no. 9, p. 850–63, 1993. [102] B. Bauer-Kirpes, W. Schlegel and R. Boesecke, "Display of organs and isodoses as 3d objects for 3d therapy planning," Int. J. Radiat. Oncol. Biol. Phys., vol. 13, p. 135–40, 1987. [103] R. Mohan, G. Barest, L. J. Brewster, C. S. Chui and G. J. Kutcher, "A comprehensive three-dimensional radiation treatment planning system," Int. J. Radiat. Oncol. Biol. Phys., vol. 15, p. 481–95, 1988. [104] G. E. Hanks, T. E. Schultheiss, A. L. Hanlon, M. Hunt and W. R. Lee, "Optimization of conformal radiation treatment of prostate cancer: report of a dose escalation study," Int. J. Radiat. Oncol. Biol. Phys., vol. 37, p. 543–50, 1997. [105] L. Tang and G. Hamarneh, "Medical image registration: a review," in Medical Imaging: 193 Technology and Applications, Boca Raton, FL, CRC, 2013, p. 619–60. [106] G. Sibley, "A comparison of spontaneous and nerve-mediated activity in bladder muscle from man, pig and rabbit," J. Physiol., vol. 354, p. 431–43, 1984. [107] R. Crowe and G. Burnstock, "A histochemical and immunohistochemical study of the autonomic innervation of the lower urinary tract of the female pig. Is the pig a good model for the human bladder and urethra?," J. Urol., vol. 141, p. 414–22, 1989. [108] S. Dahms, H. Piechota, R. Dahiya and T. Lue, "Composition and biomechanical properties of the bladder acellular matrix graft: comparative analysis in rat, pig and human," Br. J. Urol., vol. 82, pp. 411-9, 1998. [109] Y. Hu, T. J. Carter, H. U. Ahmed, M. Emberton and C. Allen, "Modelling prostate motion for data fusion during image-guided interventions," IEEE Trans. Med. Imaging, vol. 30, p. 1887–900, 2011. [110] S. Wognum, S. Heethuis, T. Rosario and M. Hoogeman, "Validation of deformable image registration algorithms on CT images of ex vivo porcine bladders with fiducial markers," Med. Phys., vol. 41, no. 7, p. 071916, 2014. [111] L. Bondar, M. S. Hoogeman, E. M. V. Osorio and B. J. Heijmen, "A symmetric nonrigid registration method to handle large organ deformations in cervical cancer patients," Med. Phys., vol. 37, p. 3760–72, 2010. [112] N. Kirby, C. Chuang, U. Ueda and J. Pouliot, "The need for application-based adaptation 194 of deformable image registration," Med. Phys., vol. 40, p. 011702, 2013. [113] R. Varadhan, G. Karangelis, K. Krishnan and S. Hui, "A framework for deformable image registration validation in radiotherapy clinical applications," J. Appl. Clin. Med. Phys., vol. 14, p. 4066, 2013. [114] S. Belongie, J. Malik and J. Puzicha, "Shape matching and object recognition using shape contexts," IEEE Trans. Pattern Anal. Mach. Intell., vol. 24, p. 509–22, 2002 . [115] A. Myronenko and X. Song, "Point-Set Registration: Coherent Point Drift," IEEE Trans. Pattern Anal., vol. 32, no. 12, pp. 2262-2275, 2010. [116] B. Jian and B. C. Vemuri, " Robust point set registration using gaussian mixture models," IEEE Trans. Pattern Anal. Mach. Intell., vol. 33, p. 1633–45, 2011. [117] V. Chalana and Y. Kim, "A methodology for evaluation of boundary detection algorithms on medical images," IEEE Trans. Med. Imaging, vol. 16, p. 642–52, 1997. [118] A. F. Frangi, W. J. Niessen and M. A. Viergever, "Three-dimensional modeling for functional analysis of cardiac images, a review," IEEE Trans. Med. Imaging, vol. 20, p. 2–5, 2001. [119] R. T. Whitaker, "A level-set approach to 3D reconstruction from range data," Int. J. Comput. Vis., vol. 29, p. 203–31, 1998. [120] G. C. Sharp, S. W. Lee and D. K. Wehe, "ICP registration using invariant features," IEEE Trans. Pattern Anal. Mach. Intell., vol. 24, p. 90–102, 2002. 195 [121] J. Salvi, C. Matabosch, D. Fofi and J. Forest, "A review of recent range image registration methods with accuracy evaluation," Image Vis. Comput., vol. 25, p. 578–96, 2007. [122] M. Greenspan and G. Godin, "A nearest neighbor method for efficient ICP," IEEE Proc. Third Int. Conf. on 3D Digital Imaging and Modeling, p. 161–8, 2001. [123] T. Jost and H. Hugli, "A multi-resolution scheme icp algorithm for fast shape registration," IEEE Proc. First Int. Symp. on 3D Data Processing Visualization and Transmission, p. 540–3, 2002. [124] T. Zinβer, J. Schmidt and H. Niemann, "A refined ICP algorithm for robust 3D correspondence estimation," IEEE Proc. 2003 Int. Conf. on Image Processing, vol. 2, p. II–695, 2003. [125] S. Rusinkiewicz and M. Levoy, "Efficient variants of the ICP algorithm on 3D Digital Imaging and Modeling," IEEE Proc. Third Int. Conf., p. 145–52, 2001. [126] D. J. Kroon, "Segmentation of the mandibular canal in cone-beam CT data," 2011. [Online]. Available: http://eprints.eemcs.utwente.nl/21056/01/thesis_D_Kroon.pdf. [127] C. Grossmann, H. G. Roos and M. Stynes, Numerical Treatment of Partial Differential Equations, Berlin: Springer, 2007. [128] S. Gold, R. A, C. P. Lu, S. Pappu and E. Mjolsness, "New algorithms for 2d and 3d point matching: pose estimation and correspondence," Pattern Recognit., vol. 31, p. 1019–31, 1998. 196 [129] A. L. Yuille, "Generalized deformable models, statistical physics, and matching problems," Neural Comput, vol. 2, p. 1–24, 1990 . [130] F. L. Bookstein, "Principal warps: thin-plate splines and the decomposition of deformations," IEEE Trans. Pattern Anal. Mach. Intell., vol. 11, p. 567–85, 1989 . [131] G. Wahba, Spline Models for Observational Data, vol. 59, Philadelphia, PA: Society for Industrial and Applied Mathematics, 1990. [132] P. Claes, D. Vandermeulen, S. De Greef and G. Willems, "Craniofacial reconstruction using a combined statistical model of face shape and soft tissue depths: methodology and validation," Forensic Sci. Int., vol. 159, p. S147–58, 2006. [133] H. Chen and B. Bhanu, "Human ear recognition in 3d," IEEE Trans. Pattern Anal. Mach. Intell., vol. 29, p. 718–37, 2007. [134] E. M. V. Osorio, M. S. Hoogeman, A. Al-Mamgani and D. N. Teguh, "Local anatomic changes in parotid and submandibular glands during radiotherapy for oropharynx cancer and correlation with dose, studied in detail with nonrigid registration," Int. J. Radiat. Oncol. Biol. Phys., vol. 70, p. 875–82, 2008. [135] J. Munkres, "Algorithms for the assignment and transportation problems," J. Soc. Ind. Appl. Math., vol. 5, p. 32–8, 1957. [136] S. Lee, G. Wolberg and S. Y. Shin, "Scattered data interpolation with multilevel b-splines," IEEE Trans. Vis. Comput. Graph., vol. 3, p. 228–44, 1997. 197 [137] S. R. Eliason, Maximum Likelihood Estimation: Logic and Practice, vol. 96, Newbury Park, CA: Sage, 1993. [138] F. M. Dekking, C. Kraaikamp, H. P. Lopuha and L. E. Meester, Maximum Likelihood, Berlin: Springer, 2005, p. 313–27. [139] B. S. Everitt and D. J. Hand, Finite Mixture Distributions, London: Chapman and Hall, 1981. [140] B. G. Lindsay, Mixed Models: Theory, Geometry, and Applications, Hayward, CA: Institute of Mathematical Statistics, 1995. [141] A. L. Yuille and N. M. Grzywacz, "A mathematical analysis of the motion coherence theory," Int. J. Comput. Vis., vol. 3, p. 155–75, 1989. [142] A. L. Yuille and N. M. Grzywacz, "The motion coherence theory," Second Int. Conf. on Computer Vision, p. 344–53, 1989. [143] T. P. Hellebust, E. Dale, A. Skjønsberg and D. R. Olsen, "Inter fraction variations in rectum and bladder volumes and dose distributions during high dose rate brachytherapy treatment of the uterine cervix investigated by repetitive ct-examinations," Radiother. Oncol., vol. 60, p. 273–80, 2001. [144] C. E. Pelloski, M. Palmer, G. M. Chronowski, A. Jhingran and J. Horton, "Comparison between CT-based volumetric calculations and ICRU reference-point estimates of radiation doses delivered to bladder and rectum during intracavitary radiotherapy for cervical 198 cancer," Int. J. Radiat. Oncol. Biol. Phys. , vol. 62, p. 131–7, 2005. [145] M. Cheung, S. Tucker and L. Dong, "Investigation of bladder dose and volume factors influencing late urinary toxicity after external beam radiotherapy for prostate cancer," Int. J. Radiat. Oncol. Biol., vol. 67, no. 4, pp. 1059-1065, 2007. [146] A. Viswanathan, E. Yorke, L. Marks, P. Eifel and W. Shipley, "Radiation Dose–Volume Effects of the Urinary Bladder," Int J Radiat Oncol Biol Phys, vol. 76, pp. S116-S122, 2010. [147] C. Ling, M. C. Schell, K. R. Working, K. Jentzsch, L. Harisiadis, S. Carabell and C. C. Rogers, "CT-assisted assessment of bladder and rectum dose in gynecological implants," Int. J. Radiat. Oncol. Biol. Phys., vol. 13, no. 10, pp. 1577-1582, 1987. [148] K. Tanderup, D. Georg, R. Pötter, C. Kirisits, C. Grau and J. Lindegaard, "Adaptive management of cervical cancer radiotherapy," Semin. Radiat. Oncol., vol. 20, no. 2, pp. 121-129, 2010. [149] S. Pavamani, D. D'Souza, L. Portelance, P. S. Craighead, A. G. Pearce, L. L. Traptow and C. M. Doll, "Image-guided brachytherapy for cervical cancer: a Canadian Brachytherapy Group survey," Brachytherapy, vol. 10, no. 5, pp. 345-351, 2011. [150] S. Nag, H. Cardenes, S. Chang, I. Das and B. Erickson, "Proposed guidelines for image-based intracavitary brachytherapy for cervical carcinoma: report from Image-Guided Brachytherapy Working Group," Int. J. Radiat. Oncol. Biol. Phys., vol. 60, no. 4, pp. 1160-72, 2004. 199 [151] S. Nag, "Controversies and new developments in gynecologic brachytherapy: Image-based intracavitary brachytherapy for cervical carcinoma," Sem. Rad. Onc, vol. 16, no. 3, pp. 164-167, 2006. [152] H. Yamashita, K. Nakagawa, K. Okuma and A. Sakumi, "Correlation Between Bladder Volume and Irradiated Dose of Small Bowel in CT-based Planning of Intracavitary Brachytherapy for Cervical Cancer," Jpn. J. Clin. Oncol., vol. 42, no. 4, pp. 302-308, 2012. [153] N. Patra, K. Manir, S. Basu and J. Goswami, "Effect of bladder distension on dosimetry of organs at risk in computer tomography based planning of high-dose-rate intracavitary brachytherapy for cervical cancer," J. Contemp. Brachytherapy, vol. 5, no. 1, pp. 3-9, 2013. [154] G. Harmon, B. Chinsky, M. Surucu and M. Harkenrider, "Bladder distension improves the dosimetry of organs at risk during intracavitary cervical high-dose-rate brachytherapy," Brachytherapy, vol. 15, no. 1, pp. 30-34, 2016. [155] A. Whitley, A. Dragovic, S. Shen and S. Wu, "Effects of Full Versus Empty Bladder on Total Equivalent Dose to Organs at Risk in 3D Image-Based Planning of High-Dose-Rate Intracavitary Brachytherapy for Cervical Cancer," Int. J. Radiat. Oncol. Biol. Phys., vol. 87, no. 2, p. Supplement 1 S404, 2013. [156] E. Andersen, K. Noe and T. Sørensen, "Simple DVH parameter addition as compared to deformable registration for bladder dose accumulation in cervix cancer brachytherapy," Radiother. Oncol., vol. 107, no. 1, pp. 52-57, 2013. 200 [157] B. Teo, L. Millar, X. Ding and L. Lin, "Assessment of cumulative external beam and intracavitary brachytherapy organ doses in gynecologic cancers using deformable dose summation," Radiother. Oncol., vol. 115, no. 2, pp. 195-202, 2015. [158] L. Marks, P. Carroll, T. Dugan and M. Anscher, "The response of the urinary bladder, urethra, and ureter to radiation and chemotherapy," Int J Radiat Oncol Biol Phys, vol. 31, pp. 1257-1280, 1995. [159] L. Hathout, M. Folkert, M. Kollmeier, Y. Yamada, G. Cohen and M. Zelefsky, "Dose to the Bladder Neck Is the Most Important Predictor for Acute and Late Toxicity After Low-Dose-Rate Prostate Brachytherapy: Implications for Establishing New Dose Constraints for Treatment Planning," Int J Radiat Onc Phys, vol. 90, pp. 312-319, 2014. [160] M. Amsbaugh, N. Bhatt, T. Hunter and e. al., "Computed tomography-planned interstitial brachytherapy for locally advanced gynecologic cancer: Outcomes and dosimetric predictors of urinary toxicity," Brachytherapy, vol. 15, pp. 49-56, 2016. [161] J. Routledge, M. Burns, R. Swindell and V. Khoo, "Evaluation of the LENT-SOMA scales for the prospective assessment of treatment morbidity in cervical carcinoma," Int. J. Radiat. Oncol., vol. 56, no. 2, pp. 502-510, 2003. [162] R. Mazeron, I. Dumas, E. Rivin and e. al., "D2cm³/DICRU ratio as a surrogate of bladder hotspots localizations during image-guided adaptive brachytherapy for cervical cancer: assessment and implications in late urinary morbidity analysis," Brachytherapy, vol. 14, pp. 300-307, 2015. 201 [163] S. Nilsson, B. Norlén and A. Widmark, "A systematic overview of radiation therapy effects in prostate cancer," Acta. Oncol., vol. 43, no. 4, pp. 316-381, 2004. [164] T. Rosewall, C. Catton, G. Currie and A. Bayley, "The relationship between external beam radiotherapy dose and chronic urinary dysfunction–a methodological critique," Radiother. Oncol., vol. 97, no. 1, pp. 40-47, 2010. [165] L. Boersma, M. V. d. Brink and A. Bruce, "Estimation of the incidence of late bladder and rectum complications after high-dose (70-78 Gy) conformal radiotherapy for prostate cancer, using dose-volume histograms," Int. J. Radiat. Oncol. Biol., vol. 41, no. 1, pp. 83-92., 1998. [166] A. Harsolia, C. Vargas, D. Yan and D. Brabbins, "Predictors for chronic urinary toxicity after the treatment of prostate cancer with adaptive three-dimensional conformal radiotherapy: Dose-volume analysis of a phase II dose-escalation study," Int. J. Radiat. Oncol. Biol., vol. 69, no. 4, pp. 1100-1109, 2007. [167] C. Cozzarini, C. Fiorino, L. D. Pozzo and F. Alongi, "Clinical factors predicting late severe urinary toxicity after postoperative radiotherapy for prostate carcinoma: a single-institute analysis of 742 patients," Int. J. Radiat. Oncol. Biol., vol. 82, no. 1, pp. 191-199, 2012. [168] K. Konishi, Y. Yoshioka, F. Isohashi and I. Sumida, "Correlation between dosimetric parameters and late rectal and urinary toxicities in patients treated with high-dose-rate brachytherapy used as monotherapy for prostate cancer," Int. J. Radiat. Oncol. Biol., vol. 75, no. 4, pp. 1003-1007, 2009. 202 [169] P. Georg, R. Pötter, D. Georg and S. Lang, "Dose Effect Relationship for Late Side Effects of the Rectum and Urinary Bladder in Magnetic Resonance Image-Guided Adaptive Cervix Cancer Brachytherapy," J. Rad. Onc. Biol. Phys., vol. 82, no. 2, p. 653–57, 2012. [170] E. M. Wiebe, K. Surry and D. P. D’Souza, "Rectal and Bladder Toxicity in Patients Receiving CT-guided HDR Interstitial Brachytherapy for Gynecologic Malignancies," Int. J. Radiat. Oncol. Biol., vol. 78, no. 3, pp. Supplement, 1 S422, 2010. [171] S.-W. Chen, J.-A. Liang, Y.-C. Hung, L.-S. Yeh and W.-C. Chang, "Geometrical Sparing Factors for the Rectum and Bladder in the Prediction of Grade 2 and Higher Complications After High-Dose-Rate Brachytherapy for Cervical Cancer," Int. J. Radiat. Oncol. Biol., vol. 75, no. 5, p. 1335–1343, 2009. [172] D. Yan, D. Jaffray and J. Wong, "A model to accumulate fractionated dose in a deforming organ," Int. J. Radiat. Oncol. Biol., vol. 44, no. 3, pp. 665-675, 1999. [173] K. Brock, D. McShan and R. T. Haken, "Inclusion of organ deformation in dose calculations," Med. Phys., vol. 30, no. 3, pp. 290-295, 2003. [174] E. Andersen, L. Muren and T. Sørensen, "Bladder dose accumulation based on a biomechanical deformable image registration algorithm in volumetric modulated arc therapy for prostate cancer," Phys. Med. Biol., vol. 57, no. 21, p. 7089, 2012. [175] K. Noe, K. Tanderup and T. Sorensen, "Surface membrane based bladder registration for evaluation of accumulated dose during brachytherapy in cervical cancer," IEEE Int. Sym. 203 on Biomedical Imaging: From Nano to Macro, pp. 1253-1256, 2011. [176] L. Xiong, A. Viswanathan, A. Stewart and S. Haker, "Deformable structure registration of bladder through surface mapping," Med. Phys., vol. 33, no. 6, pp. 1848-1856, 2006. [177] S. Sabater, I. Andres, M. Sevillano and R. Berenguer, "Dose accumulation during vaginal cuff brachytherapy based on rigid/deformable registration vs. single plan addition," Brachytherapy, vol. 13, no. 4, pp. 343-351, 2014. [178] S. Jamema, U. Mahantshetty, E. Andersen, K. Noe and T. Sørensen, "Uncertainties of deformable image registration for dose accumulation of high-dose regions in bladder and rectum in locally advanced cervical cancer," Brachytherapy, vol. 14, pp. 953-962, 2015. [179] P. Rubin, L. S. Constine, L. F. Fajardo and T. L. Phillips, "Overview of late effects normal tissues (LENT) scoring system," Radiother. Oncol., vol. 35, pp. 9-10, 1995. [180] S. Lang, C. Kirisits, J. Dimopoulos and D. Georg, "Treatment Planning for MRI Assisted Brachytherapy of Gynecologic Malignancies Based on Total Dose Constraints," Int. J. Radiat. Oncol., vol. 69, no. 2, pp. 619-627, 2007. [181] H. Chen, Z. Zhong, Y. Liao and A. Pompoš, "A non-rigid point matching method with local topology preservation for accurate bladder dose summation in high dose rate cervical brachytherapy," Phys. Med. Biol. , vol. 61, no. 3, p. 1217, 2016. [182] R. Zakariaee, G. Hamarneh, C. Brown and I. Spadinger, "Validation of non-rigid point-set registration methods using a porcine bladder pelvic phantom," Phys. Med. Biol., vol. 61, 204 no. 2, pp. 825-54, 2016. [183] C. Haie-Meder, R. Pötter and E. V. Limbergen, "Recommendations from Gynaecological (GYN) GEC-ESTRO Working Group (I): concepts and terms in 3D image based 3D treatment planning in cervix cancer brachytherapy with emphasis on MRI assessment of GTV and CTV," Radiother. Oncol., vol. 74, no. 3, pp. 235-245, 2005. [184] R. Pötter, C. Haie-Meder and E. V. Limbergen, "Recommendations from Gynaecological (GYN) GEC ESTRO Working Group (II): concepts and terms in 3D image based treatment planning in cervix cancer brachytherapy: aspects of 3D imaging, radiation physics, radiobiology, and 3D dose volume parameters," Radiother. Oncol., vol. 78, no. 1, pp. 67-77, 2006. [185] J. Cohen, Statistical power analysis for the behavioral sciences, Hillsdale, New Jersey: Lawrence Erlbaum Associates, 1988. [186] R. Miller, Simultaneous Statistical Inference, 2nd ed., New York: Springer Verlag, 1981. [187] S. Holm, "A simple sequentially rejective multiple test procedure," Scandinavian Journal of Statistics, vol. 6, no. 2, pp. 65-70, 1979. [188] A. Tamhane, "Multiple comparisons," Handbook of Statistics , vol. 13, pp. 587-630, 1996. [189] T. Farsund, "Cell kinetics of mouse urinary bladder epithelium," Arch. B. Cell Pathol, vol. 18, pp. 35-49, 1975. [190] F. A. Stewart, J. Denekamp and D. G. Hirst, "Proliferation kinetics of the mouse bladder 205 after irradiation," Cell Tissue Kinet., vol. 13, pp. 75-89, 1980. [191] F. A. Stewart, "The proliferative and functional response of mouse bladder to treatment with radiation and cyclophosphamide," Radiother. Oncol., vol. 4, pp. 353-362, 1985. [192] S. Bentzen, W. Dorr, M. Anscher and e. al, "Normal tissue effects: reporting and analysis," Semin Radiat Oncol, vol. 13, p. 189–202, 2003. [193] D. Dearnaley, E. Hall, D. Lawrence and e. al., "Phase III pilot study of dose escalation using conformal radiotherapy in prostate cancer: PSA control and side effects," Br J Cancer, vol. 92, p. 488–98, 2005. [194] S. Peeters, W. Heemsbergen, W. van Putten and e. al, "Acute and late complications after radiotherapy for prostate cancer: results of a multicenter randomized trial comparing 68 Gy to 78 Gy," Int J Radiat Oncol Biol Phys, vol. 61, p. 1019–34, 2005. [195] V. Fonteyne, G. Villeirs, N. Lumen and e. al., "Urinary toxicity after high dose intensity modulated radiotherapy as primary therapy for prostate cancer," Radiother Oncol, vol. 92, p. 42–7, 2009. [196] J. Cox, J. Stetz and P. TF, "Toxicity criteria of the Radiation Therapy Oncology Group (RTOG) and the European Organization for Research and Treatment of Cancer (EORTC)," Int J Radiat Oncol Biol Phys, vol. 31, p. 1341–6, 1995. [197] S. Mangar, K. Foo, A. Norman and e. al, "Evaluating the effect of reducing the highdose volume on the toxicity of radiotherapy in the treatment of bladder cancer," Clin Oncol , vol. 206 18, p. 466–73, 2006. [198] K. Nkiwane, R. Pötter, L. Fokdal and P. Hoskin, "se of bladder dose points for assessment of the spatial dose distribution in the posterior bladder wall in cervical cancer brachytherapy and the impact of applicator position," Brachytherapy, vol. 14, no. 2, pp. 252-259, 2015. [199] F. Faul, E. Erdfelder, A. G. Lang and A. Buchner, " G*Power 3: A flexible statistical power analysis program for the social, behavioral, and biomedical sciences," Behavior Research Methods, vol. 39, pp. 175-191, 2007. 207 Appendices Appendix A : The LENT-SOMA Questionnaire Used at BC Cancer Agency for Collecting Late GU and GI Toxicity Data for Cervical Cancer Patients. 208 209 210 211 212 213 Appendix B : Boxplots of Frequency, Dysuria, Nocturia, and Stream, Illustrating the Descriptive Statistics for the Registered Dosimetric, Volumetric, and Bladder-Neck-Dose Parameters for the Case and Control in G2/G3-Cutoff Groupings. 214 215 216 217 218 219 220
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Localized bladder dose accumulation in multi-fraction...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Localized bladder dose accumulation in multi-fraction cervical cancer brachytherapy Zakariaee Kouchaksaraee, Roja 2016
pdf
Page Metadata
Item Metadata
Title | Localized bladder dose accumulation in multi-fraction cervical cancer brachytherapy |
Creator |
Zakariaee Kouchaksaraee, Roja |
Publisher | University of British Columbia |
Date Issued | 2016 |
Description | Radiation therapy in the definitive treatment of locally advanced carcinoma of the cervix consists of external beam radiotherapy (EBRT) combined with image-guided high-dose-rate (IG HDR) intracavitary brachytherapy (ICBT). IG HDR-ICBT is a relatively new, advanced form of brachytherapy treatment planning and delivery that still relies largely on dose criteria based either on clinical experience using older techniques, or on limited data. Calculation of cumulative dose received over multiple treatment fractions currently utilizes dose-volume-histograms (DVH), which do not provide information about the spatial distribution of dose within the structures of interest. Since most of the organs at risk (OAR) for this site, like the bladder and rectum, are highly deformable, DVH data summed over multiple treatment fractions do not provide accurate estimates of the cumulative dose to specific regions of the organ, and therefore may not be the most appropriate metric to use in treatment planning and dose assessment. The primary goal of this thesis was, therefore, to develop methods to more accurately quantify “locally accumulated” dose to the bladder-wall in multi-fraction IG HDR-ICBT for cervical cancer using deformable registration, and to apply these to a study of locally-accumulated dosimetric parameters as predictors of late urinary toxicity. To this end, different deformable image and point-set registration methods were evaluated, using phantom data, for their ability to register bladder-wall contours of various sizes. The best of these was retrospectively applied to point-sets representing the bladder-wall in multiple HDR treatment fractions for 60 cervical cancer patients treated at the BC Cancer Agency. The transformation maps obtained from the registrations were used to calculate cumulative dose parameters for the bladder-wall and the urethral opening to the bladder (bladder neck). Patients were divided into Case and Control groups based on urinary toxicity scores for different symptoms, and the ability of the cumulative parameters to predict toxicity was evaluated. It was concluded that some locally-cumulative parameters estimated using our method have improved potential for predicting urinary toxicity, as compared to traditional DVH-based parameters, in our data set. In addition, dose to small volumes around the bladder neck was found to be a predictor of incontinence. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2017-01-21 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0340529 |
URI | http://hdl.handle.net/2429/60173 |
Degree |
Doctor of Philosophy - PhD |
Program |
Physics |
Affiliation |
Science, Faculty of |
Degree Grantor | University of British Columbia |
GraduationDate | 2017-02 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2017_february_zakariaeekouchaksaraee_roja.pdf [ 5.94MB ]
- Metadata
- JSON: 24-1.0340529.json
- JSON-LD: 24-1.0340529-ld.json
- RDF/XML (Pretty): 24-1.0340529-rdf.xml
- RDF/JSON: 24-1.0340529-rdf.json
- Turtle: 24-1.0340529-turtle.txt
- N-Triples: 24-1.0340529-rdf-ntriples.txt
- Original Record: 24-1.0340529-source.json
- Full Text
- 24-1.0340529-fulltext.txt
- Citation
- 24-1.0340529.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0340529/manifest