A Robotic Needle Guide for Prostate Brachytherapy with Pre-operative to Intra-operative Prostate Volumes Registration by Thomas Diego Prananta B.A.Sc., The University of British Columbia, 2006 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in The Faculty of Graduate Studies (Electrical Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) August 2010 c Thomas Diego Prananta 2010 ° Abstract The conventional prostate brachytherapy approach is limited by needle positioning accuracy, needle trajectory option, and prostate motion and deformation between the pre-operative volume study and the seed implant procedure. These limitations increase the risks of post implant complications. In this thesis we develop a robotic needle guide to improve prostate brachytherapy needle placement accuracy and trajectory option as well as a pre-operative to intra-operative prostate volume registration algorithm to address the issue of prostate motion and deformation. Our four degrees of freedom robot provides X-Y axes translational accuracy of 0.12 and 0.1 mm compared to the 5 mm accuracy of the standard needle guide. The robot also provides yaw and pitch angulations with 0.050 accuracy which can be used to reach prostate regions blocked by pubic arch interference. The robot is adaptable to conventional brachytherapy apparatus without adding the clinical procedure time and can be used manually in the case of electronic control failure. The registration approach is based on fitting prostate surfaces into ellipsoids. Pre-operative and intra-operative sagittal view-based volume data are contoured using a novel semi automatic sagittal view-based segmentation algorithm. The resulting contours are fit into ellipsoids whose parameters centers, orientations, and radii - are used to calculate the registration matrix. The accuracy of the registration algorithm was compared with Optotrak measurement as the gold standard and with the Iterative Closest Point (ICP) algorithm. The result shows that the orientation of the ellipsoid fit is sensitive to user initialization points causing up to 5 mm translational errors and 5.50 angular error. The comparison with ICP shows that the ellipsoid fitting based algorithm is faster but less accurate. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Acronyms . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Prostate Brachytherapy . . . . . . . . . . . . . . . . . . . . . 1.2 Thesis Motivation and Objective . . . . . . . . . . . . . . . . 1.2.1 Motivation for a Robotic Needle Guide . . . . . . . . 1.2.2 Motivation for Pre-operative to Intra-operative Prostate Volume Registration . . . . . . . . . . . . . . . . . . . 1.2.3 Thesis Objectives . . . . . . . . . . . . . . . . . . . . 1.3 Literature Review and Thesis Approach . . . . . . . . . . . . 1.3.1 Literature Review . . . . . . . . . . . . . . . . . . . . 1.3.1.1 Literature Review on Robotic Needle Guides 1.3.1.2 Literature Review on Prostate Volumes Registration . . . . . . . . . . . . . . . . . . . . 1.3.2 Thesis Approach . . . . . . . . . . . . . . . . . . . . . 1.4 Thesis Organization . . . . . . . . . . . . . . . . . . . . . . . 1 1 6 6 8 10 10 10 10 12 13 15 2 Brachyguide: a Robotic Needle Guide for Prostate Brachytherapy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.1 Robot Design and Implementation . . . . . . . . . . . . . . . 16 2.1.1 Control and Graphical User Interface (GUI) . . . . . 21 2.2 Brachyguide Characteristics . . . . . . . . . . . . . . . . . . 23 iii Table of Contents 2.2.1 Devia. . . . . . . . . . . . . . . . . . . . . . . . . 25 26 26 26 34 3 Segmentation Based on Sagittal Images . . . . . . . . . . . . 3.1 The Transversal Segmentation Algorithm . . . . . . . . . . . 3.1.1 User Initialization . . . . . . . . . . . . . . . . . . . . 3.1.2 Image Warping . . . . . . . . . . . . . . . . . . . . . 3.1.3 Initial Ellipse Fit . . . . . . . . . . . . . . . . . . . . 3.1.4 IMMPDAF Edge Detection . . . . . . . . . . . . . . . 3.1.5 Symmetrization and Final Ellipse Fit . . . . . . . . . 3.1.6 Image Unwarping . . . . . . . . . . . . . . . . . . . . 3.2 Modifications on the Transversal Segmentation Algorithm . . 3.2.1 Initialization Points . . . . . . . . . . . . . . . . . . . 3.2.2 Warping . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.3 Selective IMMPDAF Edge Detector . . . . . . . . . . 3.2.4 Stable Ellipse Fitting Algorithm . . . . . . . . . . . . 3.3 The Sagittal View-based Segmentation Algorithm . . . . . . 3.4 Extending the Sagittal Segmentation Algorithm . . . . . . . 3.5 Segmentation Results on Phantom and Patient Volume Data 3.5.1 Segmentation on CIRS Phantom Volume Data . . . . 3.5.2 Segmentation on Patient Volume Data . . . . . . . . 37 37 38 38 40 40 40 40 40 41 41 45 45 46 49 51 51 56 2.3 2.4 Translational Accuracy and Straight Motion tion . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 Angular Accuracy . . . . . . . . . . . . . . . 2.2.3 Stiffness . . . . . . . . . . . . . . . . . . . . Seed Implant Study on a Phantom . . . . . . . . . . Additional Features and Updates on the Robot . . . 4 Surface-based Pre-operative to Intra-operative Registration of the Prostate . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Registration Using Ellipsoid Parameters . . . . . . . . . . . . 4.2 Registration Experiment: Simulation . . . . . . . . . . . . . 4.3 Registration Experiment: CIRS Phantom . . . . . . . . . . . 4.3.1 Experiment Setup . . . . . . . . . . . . . . . . . . . . 4.3.2 TRUS Probe Calibration . . . . . . . . . . . . . . . . 4.3.3 Registration Experiment Results . . . . . . . . . . . . 4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.1 Ellipsoid Orientation Fit Differences Cause Translation and Distance Registration Errors . . . . . . . . . 4.4.2 Sensitivity of Final Ellipsoid Parameters to Initialization Points . . . . . . . . . . . . . . . . . . . . . . . . 4.4.3 Role of the Kalman Filter Segmentation . . . . . . . 61 61 68 74 74 79 79 84 84 86 90 iv Table of Contents 4.5 Iterative Closest Point (ICP) . . . . . . . . . . . . . . . . . . 93 5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 5.1 A Robotic Needle Guide for Prostate Brachytherapy . . . . . 97 5.1.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . 97 5.1.2 Contributions . . . . . . . . . . . . . . . . . . . . . . 98 5.1.3 Future Work . . . . . . . . . . . . . . . . . . . . . . . 99 5.2 Pre-operative to Intra-operative Registration of Prostate Volumes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.2.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.2.2 Contributions . . . . . . . . . . . . . . . . . . . . . . 100 5.2.3 Future Work . . . . . . . . . . . . . . . . . . . . . . . 101 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 Appendices A Ellipse Fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . A.1 Direct Least Squares Fitting of Ellipses . . . . . . . . . . . A.2 Numerically Stable Least Squares Fitting . . . . . . . . . . A.3 Converting Algebraic Parameters to Geometric Parameters . . . . 112 112 114 116 B Ellipsoid Fitting . . . . . . . . . . . . . . . . . . . . . . . . . . B.1 Direct Least Square Fitting of Ellipsoids . . . . . . . . . . B.2 Converting Algebraic Parameters to Geometric Parameters B.3 Generating Ellipses from an Ellipsoid . . . . . . . . . . . . . . . . 120 120 121 124 C The IMMPDAF Edge Detector . . . . . . . . . . . . . . . C.1 State Estimate Initialization . . . . . . . . . . . . . . . . C.2 Discrete Time Kalman Filter . . . . . . . . . . . . . . . . C.2.1 State Space Model . . . . . . . . . . . . . . . . . . C.2.2 Kalman Filter Iteration . . . . . . . . . . . . . . . C.2.3 Obtaining the Measurements: the PDAF Concept C.3 Combining the Estimates From the Two Modes . . . . . . . . . . . . 127 127 128 128 130 131 132 . . . . . . . D Probe Calibration Using Stradwin and the Single Wall Technique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 D.1 The Framework of Single Wall Calibration Technique . . . . 134 D.2 US Image and Optotrak Data Collection . . . . . . . . . . . 136 D.3 Single Wall Calibration Technique Using Stradwin . . . . . . 140 v Table of Contents D.4 Automatic Line Detection and Optimization . . . . . . . . . 142 E Registration Experiment Results . . . . . . . . . . . . . . . . E.1 Registration Experiment Results: the Sparse Volume Registration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . E.2 Registration Experiment Results: the Full Volume Registration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . E.3 ICP Registration Results . . . . . . . . . . . . . . . . . . . . 144 144 147 150 F Ellipsoid Fitting Results . . . . . . . . . . . . . . . . . . . . . 153 F.1 Ellipsoid Fitting Results: the Sparse Volume Registration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 F.2 Ellipsoid Fitting Results: the Full Volume Registration . 155 G Determining the Motor Specifications G.1 The Rotation Stages . . . . . . . . . . G.2 The X-axis Linear Stage . . . . . . . G.3 The Y-axis Linear Stage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157 157 157 160 H Initial Ellipsoid Fit Sensitivity Analysis . . . . . . . . . . . . 163 I Ethics Certificates for Full Board Approval . . . . . . . . . 167 vi List of Tables 2.1 2.2 2.3 3.1 3.2 4.1 4.2 4.3 4.4 4.5 Brachyguide characteristics table. . . . . . . . . . . . . . . . . Seed placement errors based on post implant CT scan. . . . . Updated translational accuracy and repeatability measurement results. . . . . . . . . . . . . . . . . . . . . . . . . . . . Validation of segmentation results on CIRS phantom volume data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Validation of segmentation results on patient volume data. . . Ellipsoid fitting error under various simulation cases. . . . . . Registration errors under various simulation cases. . . . . . . Summary of registration errors: the sparse volume registration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Summary of registration errors: the full volume registration. Summary of registration errors: ICP. . . . . . . . . . . . . . . E.1 Registration tration. . . E.2 Registration tion. . . . . E.3 Registration 24 32 35 56 59 72 73 80 82 95 experiment results: the sparse volume regis. . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 experiment results: the full volume registra. . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 experiment results - ICP. . . . . . . . . . . . . . 151 F.1 Ellipsoid fit results for the sparse volume registration. . . 154 F.2 Ellipsoid fit results for the full volume registration. . . . . 156 vii List of Figures 1.1 1.2 1.3 1.4 Pre operative (pre-op) volume study. . . . . . . . . . . The conventional brachytherapy apparatus. . . . . . . Pre-op volume study images from the treatment plan. The conventional needle guide template. . . . . . . . . . . . . 3 4 5 7 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 2.10 2.11 2.12 2.13 2.14 2.15 2.16 Brachyguide kinematics design and implementation. . . . . . Translation stages. . . . . . . . . . . . . . . . . . . . . . . . . Robot wrist. . . . . . . . . . . . . . . . . . . . . . . . . . . . . Brachyguide unique features. . . . . . . . . . . . . . . . . . . Control block diagram. . . . . . . . . . . . . . . . . . . . . . . Brachyguide Graphical User Interface (GUI). . . . . . . . . . Translational accuracy measurement setup. . . . . . . . . . . Translational accuracy measurement results. . . . . . . . . . . Angular accuracy measurement results. . . . . . . . . . . . . Brachyguide stiffness test setup. . . . . . . . . . . . . . . . . . Brachyguide stiffness measurement results. . . . . . . . . . . . Seed implant study on phantom. . . . . . . . . . . . . . . . . Seed placement errors histograms. . . . . . . . . . . . . . . . Post implant seed detection results. . . . . . . . . . . . . . . . New digital linear calipers. . . . . . . . . . . . . . . . . . . . . Updated translational accuracy and repeatability measurements. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 18 19 20 22 23 25 27 28 29 30 31 33 33 34 Transversal segmentation algorithm steps. . . . . . . . . . . . Warping factor calculation. . . . . . . . . . . . . . . . . . . . Comparisons of initialization points. . . . . . . . . . . . . . . Modifications in image warping. . . . . . . . . . . . . . . . . . Warping on the sagittal view based images. . . . . . . . . . . Steps of the sagittal segmentation algorithm. . . . . . . . . . Initialization points for the extended sagittal segmentation algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 39 42 43 44 48 3.1 3.2 3.3 3.4 3.5 3.6 3.7 . . . . . . . . . . . . 36 50 viii List of Figures 3.8 Extended sagittal segmentation results on two phantom volume data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Distance calculation between contour points. . . . . . . . . . 3.10 Comparison of interpolated sagittal segmentation results and transversal segmentation results on phantom data. . . . . . . 3.11 Extended sagittal segmentation results on two phantom volume data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.12 Comparison of interpolated sagittal segmentation results and transversal segmentation results on patient data. . . . . . . . 4.1 4.2 4.3 4.4 4.5 4.6 4.7 Coordinate system convention. . . . . . . . . . . . . . . . . . Ellipsoid fitting based registration. . . . . . . . . . . . . . . . Ellipsoid parameters definition. . . . . . . . . . . . . . . . . . Simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . Experiment setup. . . . . . . . . . . . . . . . . . . . . . . . . The coordinate frames in the experiment setup. . . . . . . . . Ellipsoid parameters sensitivity to changes in initialization point location in X-axis. . . . . . . . . . . . . . . . . . . . . . 4.8 Ellipsoid parameters sensitivity to changes in initialization point location in Y-axis. . . . . . . . . . . . . . . . . . . . . . 4.9 Ellipsoid parameters sensitivity to changes in initialization point location in Z-axis. . . . . . . . . . . . . . . . . . . . . . 4.10 Ellipsoid parameters sensitivity comparison. . . . . . . . . . . 4.11 Comparison of registration results: the full volume registration and ICP. . . . . . . . . . . . . . . . . . . . . . . . . . 52 54 55 58 60 62 63 65 71 75 77 87 88 89 92 94 A.1 Ellipse fit coordinate frames. . . . . . . . . . . . . . . . . . . 117 B.1 Coordinate frames for ellipsoid parameter conversion. D.1 D.2 D.3 D.4 D.5 D.6 Framework of US calibration. . . . . . Calibration setup. . . . . . . . . . . . Calibration setup: water bath. . . . . Calibration probe poses. . . . . . . . . Setting the calibration image scales. . Automatic line detection in Stradwin. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 137 138 139 141 142 G.1 Estimating the friction rotation stage static friction torque. . 158 G.2 Estimating the static friction coefficient of the X-axis stage’s slider. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159 G.3 Estimating the static friction of the Y-axis stage’s slider. . . . 161 ix List of Figures H.1 Initial ellipsoid parameters sensitivity to changes in initialization point location in X-axis. . . . . . . . . . . . . . . . . . . 164 H.2 Initial ellipsoid parameters sensitivity to changes in initialization point location in Y-axis. . . . . . . . . . . . . . . . . . . 165 H.3 Initial ellipsoid parameters sensitivity to changes in initialization point location in Z-axis. . . . . . . . . . . . . . . . . . . 166 I.1 I.2 I.3 I.4 Ethics Approval Certificate for patients data use: 2007-2008 . Ethics Approval Certificate for patients data use: 2008-2009 . Page 1 of the Ethics Approval Certificate for patients data use: 2009-2010 . . . . . . . . . . . . . . . . . . . . . . . . . . Page 2 of the Ethics Approval Certificate for Patients Data use: 2009-2010 . . . . . . . . . . . . . . . . . . . . . . . . . . 168 169 170 171 x List of Acronyms 2D 3D BAT BCCA CT DC DOF EBRT GUI HDR ICP IMMPDAF LCD LDR LOG MRI OEM PAI PC PD PDAF PI PVC Two dimensional Three dimensional B-mode Acquisition and Targeting British Columbia Cancer Agency Computed Tomogragphy Direct Current Degrees Of Freedom External Beam Radiation Therapy Graphical User Interface High Dose Rate Iterative Closest Points Interacting Multiple Model Probabilistic Data Association Filter Liquid Crystal Display Low Dose Rate Laplacian Of Gaussian Magnetic Resonance Imaging Original Equipment Manufacturer Pubic Arch Interference Personal Computer Proportional Derivative Probabilistic Data Association Filter Proportional Integral Polyvinyl Chloride xi List of Acronyms RANSAC RAPID RCM TRUS UMI US Random Sample Consensus Robot-Assisted Platform for Intratumoral Delivery Remote Center of Motion Transrectal Ultrasound Ultrasound-guided Motionadaptive needle-insertion Instrument Ultrasound xii Acknowledgements First of all, I would like to thank God for his blessings. Without Him, I would never grew up to be the man I am today and to be able to present this work. I would also thank my parents who have supported me in every way they can. I dedicate my thesis work to make them proud. To my brother, thanks for being an inspiration for me to do well. You have studied abroad before I did and that inspired me to follow your steps. To Dr. Salcudean (Tim) and Dr. Rohling (Rob), it has been an honour to be able to work with the brightest minds in this field. The guidance and help that they have given me are invaluable. I would like to thank Dr. William J. Morris, Dr. Ingrid Spadinger, Dan Gelbart, Victor Wang, Honza Vohradsky, Xu Wen, Mehdi Moradi, Sara Mahdavi, Sara Badiei, and Nick Chng for their assistance in parts of this project. To my dragon boat team, I have learnt a lot on the meaning of team work, sacrifice, and leadership. Lastly, I would like to thank my closest friends: Denis, Andre, Budi, Joanna, Danny, Ali, and Yoenhee Ahn. They have taken care of me in my lowest periods, kept my Matlab skills sharp, shared pints of beer, and helped me in prayers and emotional support. I am lucky to be blessed with so many good friends. Diego xiii Chapter 1 Introduction 1.1 Prostate Brachytherapy Prostate cancer is the most commonly diagnosed cancer among men in Canada with 25500 new cases and 4400 related deaths reported in 2009 [1]. At its early stage, the cancer is localized within the prostate and several treatment options are available [2]: 1. Prostatectomy 2. Radiation therapy: External beam radiation therapy (EBRT) Brachytherapy 3. Cryotherapy 4. Hormonal therapy 5. Laser therapy 6. Watchful waiting (conservative management) There are two types of brachytherapy: HDR (High Dose Rate) and LDR (Low Dose Rate). The differences between the two are the number of fractions and the method of application. In HDR brachytherapy, radioactive dosage are delivered into the prostate with removable catheters, leaving no radioactive material inside the prostate. The treatment is done in several fractions. In LDR brachtyherapy, radioactive seeds are permanently implanted into the prostate, requiring no treatment fractions. In this thesis, the term brachytherapy refers to LDR brachytherapyy. Prostatectomy, EBRT, and brachytherapy are the standard treatment options for localized prostate cancer. Among the three, brachytherapy is the most convenient treatment option [3]. The procedure takes only about an hour and is usually done on an outpatient basis. The patient can resume daily activities afterwards. In comparison, prostatectomy takes 2-3 1 1.1. Prostate Brachytherapy hours and requires a hospital stay and patient activity restrictions, while EBRT requires regular patient visits over a period of 7-8 weeks. This convenience factor along, with aging demography and increase in early detection of prostate cancer [1], makes brachytherapy a popular treatment choice in Canada. In brachytherapy, radioactive seeds (125 I or 103 P d) are implanted through the perineum into the prostate using surgical needles under template guidance and transrectal ultrasound (TRUS) imaging. The method delivers localized dosage of radiation while minimizing damage to healthy tissues and organs. This dose profile results in lower complications rates compared to prostatectomy or EBRT. Such complications include impotence, urinary toxicity (incontinence, retention, and dysuria), and renal toxicity [4, 5]. Every brachytherapy patient must undergo a prostate volume study (the pre-operative volume study), typically 2-3 weeks before the seed implant procedure [6]. The study’s purpose is to obtain the volume and shape of the patient’s prostate. This information will be used to generate radioactive dose and seed distribution plans for the implant procedure. During the pre-operative volume study, the patient lies on the operating table in a dorsolitothomy position with his leg supported by the bed’s legs stand (illustrated in Figure 1.1). A radiation oncologist mounts a TRUS probe onto the brachytherapy stepper and mount, such as the EXII stepper R LP mount (CIVCO Medical, Kalona, Iowa) [8] and [7] and Micro-Touch° inserts the probe into the patient’s rectum. Both apparatus are shown in Figure 1.2. The probe is inserted until the oncologist is able to identify the base of the prostate (the part of the prostate closest to the bladder). Often, the prostate is viewed sagittally to confirm the base location. The TRUS probe is also adjusted so that the prostate is symmetric about the center line of the image. Then, the oncologist collects a set of 2D transverse B-mode US images starting from 5 mm superior to the base up to 5 mm inferior to the prostate’s apex. The apex is the part of the prostate closest to the perineum (also illustrated in Figure 1.1). The images are captured in 5 mm steps provided by the brachytherapy stepper. In total, 13 2D B-Mode US images are collected. The 2D images are loaded into a treatment planning system such as VariSeed (Varian Medical Systems, Chalottesville, VA) [9]. Such a system, enables the medical physicist to manually contour (segment) the prostate boundary on each of the 2D images. The contours are interpolated to define the 3D volume and shape of the patient’s prostate. Based on this information, the medical physicist creates the plan for the implant which consists of the radioactive dose distribution profile, the seed location, and the needle 2 1.1. Prostate Brachytherapy Figure 1.1: A Typical pre-operative volume study setup. The patient is put in litothomy position. The TRUS probe is inserted into the rectum and transverse images of the prostate are taken at 5 mm intervals. 3 1.1. Prostate Brachytherapy Figure 1.2: The conventional EXII Brachytherapy stepper with the MicroR Touch LP° Mount. 4 1.1. Prostate Brachytherapy Figure 1.3: Pre-operative volume study images of a patient taken from a VariSeed treatment plan. In the images we can see the prostate boundary (red line) from manual segmentation and the brachytherapy seeds location (light blue dots). Oncologists use these images to manually register the probe during the implant procedure. 5 1.2. Thesis Motivation and Objective insertion plan as shown in Figure 1.3. During the implant procedure, the patient is put under anesthesia, general or regional (spinal or epidural), and is positioned in the same dorsolithotomy position during the pre-operative volume study. The radiation oncologist sets up the TRUS probe on the brachytherapy stepper and mount, then inserts the probe into the patient’s rectum. To make sure that the stepperprobe-patient setup is the same as in the pre-operative volume study, the oncologist manually registers the prostate midgland (base+25 mm) B-Mode image with the midgland image on the printed treatment plan (illustrated by Figure 1.3). The prostate motion is recognized by changes in its shape and the location of its discernible boundary. The oncologist rotates and translates the TRUS probe via the stepper’s and mount’s adjustment dials to reproduce the pre-operative volume study setup. After the patient setup is confirmed, the radiation oncologist performs the needle insertions one by one, according to the plan, as read out by the medical physicist. Needles are inserted through the grid holes in the template guide. As shown in Figure 1.4, the holes are spaced 5 mm apart and constrain the needle trajectory to be parallel to the TRUS probe axis. The oncologist confirms the needle locations using real time B-mode US images. After deploying the seeds, the oncologist checks the seed locations using both B-mode US images and real time fluoroscopy imaging. 1.2 1.2.1 Thesis Motivation and Objective Motivation for a Robotic Needle Guide Improving needle placement accuracy should result in a more conformal radioactive dosage and, consequently, less risk of post brachytherapy complications. Needle placement error is one of the sources of seed placement deviation from the treatment plan [10]. The deviation can cause significant radiation overdose [11] to regions of the prostate causing increased risk of post implant complications [5]. Pubic arch interference (PAI) is one of the main concerns in prostate brachytherapy. It is caused by the pubic bone partially blocking the anterior portion of the prostate. The bone will obstruct parallel needle trajectories thus preventing seeds from being implanted in the blocked region. If PAI is detected in the pre-operative volume images, the patient may be prevented from following the brachytherapy treatment [12] or may be prescribed hormone therapy to reduce his prostate volume. Hormone therapy, however, comes with its own side effects. Figure 1.1 illustrates the occurrence of PAI. 6 1.2. Thesis Motivation and Objective Figure 1.4: The conventional needle guide template with 5 mm grid hole spacing. The holes enforce parallel needles trajectories. 7 1.2. Thesis Motivation and Objective The pubic bone (red) obstructs the needle from some of the anterior parts of the prostate. Unexpected PAI during the implant procedure is typically overcome by performing angled needle insertions. The oncologist adjusts the brachytherapy stepper to angle the template guide with respect to the insertion point to provide better needle access. This approach is limited by the stepper’s workspace and the template guide’s 5 mm hole spacing. Moreover, since the TRUS probe is fixed on the brachytherapy stepper, the already registered probe-patient setup changes. The change introduces seed misplacement errors on the remaining needle insertions. Replacing the current template guide with a robotic needle guide would help reduce the seed placement errors and solve the PAI problem. The robot can potentially provide less than 5 mm needle spacing resolution for higher needle placement precision and can add two more degrees of freedom (pitch and yaw) to angulate the needle. When PAI occurs, a new needle path incorporating pitch and yaw can be calculated. Since the robot provides needle angulation, treament plans can be optimized to include needle angulation in order to lower seed misplacement errors [13]. With such treatment plans, patients with high risk of PAI may be more likely to be admitted for brachytherapy with or without hormone therapy. The use of the needle guide robot does not change the operating room workflow and implant procedure time. The robot can easily be customized to fit the current brachytherapy apparatus (the brachytherapy stepper and stand). The robot helps the oncologist guide the needle but he or she will still have control over the needle insertion. Manual controls can be incorporated should the oncologist wish to rely on manual control without reverting back to the conventional template guide. 1.2.2 Motivation for Pre-operative to Intra-operative Prostate Volume Registration The issue of pre-operative to intra-operative registration occurs in both brachytherapy and EBRT. In EBRT, the patient’s prostate suffers from interfraction motion and deformation. The position, pose, and shape of the prostate changes between the planning phase and the radiation phase. Unlike EBRT, prostate brachytherapy treatment does not have any fractions. Yet, between the pre-operative volume study and the implant procedure there is a significant time period of time - typically weeks - in which prostate motion and deformation can occur. They contribute to seed placement errors and radioactive dose deviation since the treatment plan is based on the 8 1.2. Thesis Motivation and Objective pre-operative prostate volume. Errors and deviations affect the effectiveness of the treatment and the risk of post implant complications. Prostate motion and deformation are caused by external and internal factors. Patient setup error, caused by difficulties in reproducing the preoperative volume setup, is one of the external factors. Although the patient’s legs support and brachytherapy mount are fixed on the operating table, subtle variations (a few mm) in the patient’s dorsolitothomy position are not discernible by the naked eye. Another external factor is the TRUS probe. The pressure it exerts dislocates and deforms the prostate. Internal factors include the change in bladder filling, which causes bladder expansion and pushes the prostate in the anterior-posterior direction, and the relaxation of the pelvic muscles due to anesthesia. As mentioned before, to correct for prostate motion and deformation, oncologists manually register the real time midgland B-Mode image with its pre-operative one by matching the prostate’s discernible boundary. This method, although claimed to be accurate, suffers from both intra-user and inter-user variability. Poor US visualization of the prostate makes boundary recognition difficult. It is even more difficult to determine the changes in prostate pose, location, and shape based on just its boundaries. Experience and practice dictates the precision of manual registration. To the author’s knowledge, there have not been extensive publications regarding pre-operative to intra-operative registration of prostate volumes for brachytherapy. The author found related studies focusing on prostate motion due to brachytherapy needle insertion [14–16], interfraction motion compensation for EBRT treatment [17, 18] and US localization in prostate biopsy [19]. A dedicated study to quantify interfraction motion in prostate brachytherapy is warranted. Such study can asses the accuracy of manual registration and eventually provide solutions to correct for interfraction prostate motion. There are several specific benefits of developing an algorithm to register pre-operative and intra-operative prostate volumes. First, the algorithm can quantify prostate motion and deformation in brachytherapy treatment, providing a more consistent benchmark and understanding of the motion. Second, the transformation parameters from the algorithm can be used for training of manual registration. Third, applying the transformation parameters to a needle guide device or to the treatment plan may improve the quality of the implant. 9 1.3. Literature Review and Thesis Approach 1.2.3 Thesis Objectives To summarize, this thesis project has two objectives. First is to design, build, and characterize a robotic needle guide for prostate brachytherapy. Second is to develop an ultrasound B-mode image-based pre-operative to intra-operative prostate volumes registration algorithm. 1.3 1.3.1 1.3.1.1 Literature Review and Thesis Approach Literature Review Literature Review on Robotic Needle Guides Along with robotic approaches to needle guidance, several groups have explored the idea of modifying the needle guide template to provide angulation and to improve needle placement accuracy. Green et al. [20] presented a pivoting needle template that provides manual yaw angulation of needles with respect to the US probe. Hogendijk et al. also explored template guide modifications providing needle angulation and translation [21]. To the author’s knowledge, none of these ideas have been implemented for clinical use. Schneider et al. developed a three degrees of freedom (DOF) robotic needle guide for transrectal prostate brachytherapy [22]. The needle is guided by a motorized half cylinder TRUS probe sheath through one of its side channels. The needle tip location is calculated based on three parameters: the sheath rotation angle, the sheath insertion depth, and the needle insertion depth. The later two are measured manually since only the probe rotation is motorized. Although the proposed needle trajectory is novel, it is completely different compared to the conventional template guidance approach. Thus, direct use by physicians may be difficult. Bassan et al. developed a needle guide robot based on a Macro-Micro system design [23]. The seven DOF passive macro manipulator is used to grossly position the needle tip on the skin entry point. The active (motorized) micro manipulator has five DOF. Two for needle insertion axis orientation, one for needle rotation, one for needle insertion, and one for seed ejection. The micro manipulator is also equipped with a force sensor to measure the needle insertion force. Although the micro-manipulator is light weight (less than 5kg), the overall system is large and may obstruct the transperineal view and/or the fluoroscopy machine. RAPID (Robot-Assisted Platform for Intratumoral Delivery) [24, 24], developed by Yu et al., is a cart-mounted surgery module which consists 10 1.3. Literature Review and Thesis Approach of a two DOF TRUS probe driver, a three DOF gantry, a two DOF needle driver, and a one DOF seed dispenser. The gantry module enables motorized x-y translation and pitch angulation of the needle. The cart itself provides manual control of three axes of translations and three axes of rotations for the probe driver and the gantry. The module is equipped with two force sensors and one torque sensor to measure the needle insertion force and torque. RAPID completely replaces the conventional brachytherapy stepper and stand which may make clinical acceptance more difficult. Fenster et al. developed a complete 3D TRUS guided system for prostate brachytherapy which includes a 3D US imaging system and a needle guide robot [25, 26]. In their earlier work, the group used a six DOF commercial robot (Model A465 from Thermo-CRS, Burlington Ontario, Canada). A needle guide for manual insertion is attached to the robot’s arm. The robot positions the needle guide onto the skin insertion point and a physician performs the needle insertion. In their later publication, the group developed a custom made four DOF eleven revolute cylindrical closed-kinematic robot [27]. The robot is compact and equipped with mechanical brake clamps for manual positioning. The encoders’ angular positions define the position and orientation of the needle guide through the linkage’s coordinate transformation. Since the encoder readings do not directly translate to X-Y axes translation or the pitch and yaw angle of the needle, manual ease of use by physicians might be an issue. Fichtinger et al. developed several needle guidance robots. The “Acubot” [28] is a robot for percutaneous needle intervention using fluoroscopy and computed tomography (CT) - originally intended for applications such as biopsy and ablation. The design comprises of a three DOF gantry attached to the operating table, a seven DOF passive arm and a three DOF needle insertion device with a Remote Center of Motion (RCM). The same robot is later applied for prostate brachytherapy [29]. The AcuBot’s RCM design is novel but the three DOF gantry has a limited range of motion and a large structure which makes it difficult to integrate into the clinical setting. Hence, the group recently designed a lighter and more compact custom parallel robot for prostate brachytherapy [30]. The robot is compatible with a conventional brachytherapy stepper and provides four DOF needle movement (X-Y translations and pitch-yaw angulations). The robot is very light (1.3 kg) with decoupled kinematics which makes manual positioning easy. However, if the robot’s electronic control fails, the solution is to revert back to the manual template. The same group also developed a Magnetic Resonance Imaging (MRI) compatible robot for prostate brachytherapy [31]. Other robots were designed for needle guidance but not specifically for 11 1.3. Literature Review and Thesis Approach prostate brachytherapy application. Phee et al. developed a nine DOF passive robot for prostate biopsy [32]. The robot has a large structure which makes it not adaptable to the typical brachytherapy apparatus. The UMI (ultrasound-guided motion-adaptive needle-insertion instrument) robot is a two DOF needle driving manipulator attached to a five DOF passive arm developed by Hong et al. for percutaneous cholecystostomy [33]. Kettenbach et al. developed a seven DOF robot, “B-Rob I”, for general percutaneous biopsy [34]. None of these robots have been used for the brachytherapy application, possibly due to the use of different needle insertion approach that limits the operational space and dictates specific needle trajectories. 1.3.1.2 Literature Review on Prostate Volumes Registration Several research groups studied inter-fraction prostate motion and prostate localization for EBRT by using implanted markers. Wu et al. [35] and Crook et al. [36] used three implanted gold markers and portal film imaging to measure the prostate motion and patient setup error by identifying the shifts in markers’ positions. Shimizu et al. [37] developed a real time tracking radiotherapy system to localize prostate cancer by tracking a single marker using two X-ray images. The Calypso system by Willoughby et al. [38] uses three electro magnetic markers to measure the patient’s prostate motion during EBRT treatment. The group claimed that the system can also be used for real-time tracking of the prostate. Marker-based registration approaches, beside being invasive, suffer from the possibility of marker migration [39]. It is also not suitable for registering pre-operative prostate volume to the intra-operative one as the presence of the markers might interfere with the brachytherapy implant procedure. Other research groups studied prostate motion using Cine-MRI [40, 41] and CT [42–45] imaging modalities. Although both CT and Cine-MRI are not the main imaging modality used in prostate brachytherapy, it is interesting to note that these research groups, with the exception of Ghilezan et al. [40], follow similar approaches in studying prostate motion. They quantify the motion based on the changes of prostate contours which are obtained by manually contouring the prostate in their CT and Cine-MRI images. Most of these research groups reported the prostate motion in terms of translation in the anterior-posterior, superior-inferior, and left-right direction. Only vanHerk et al. [45] reported the translation and rotation of the prostate with respect to the pelvic bone. Several groups studied prostate motion and localization using ultrasound as the main imaging modality. The Fox Chase group developed a stereo12 1.3. Literature Review and Thesis Approach tactic guidance system for prostate localization during radiation therapy treatment, called B-mode Acquisition and Targeting (BAT) [43, 43, 46]. The system combines a transabdominal US probe with a position sensing arm and a touch screen. The prostate contours from the CT treatment planning system is displayed on the touch screen along with the transverse and sagittal suprapubic ultrasound images. The contours are manually shifted to overlay the prostate in the ultrasound images. The motion of the prostate is quantified by the shift. The system was compared with CT based localization by other groups [18, 47]. Xu et al. [19] register pre-operative US volume data with intra-operative US images in fused MR-TRUS image guided prostate biopsy. A few 2D intra-operative images are registered to the pre-operative volume with a stereotactic tracking system’s assistance. The registration is used in a closed-loop control to correct for prostate motion in the fused MR-TRUS images. Krupa et al. [48] developed a prostate tracking algorithm based on region tracking and speckle decorrelation for visual servoing purposes. Prostate brachytherapy is identified as one of the possible clinical applications of the work. 1.3.2 Thesis Approach The conventional needle guide template, illustrated in Figure 1.4, has 5 mm needle placement accuracy and limits the needle trajectory to be parallel. As mentioned, improving the needle placement accuracy (less than 5 mm) will help reduce the risk of post implant complications while adding the degrees of freedom of the needle guide may help patients with pubic arch interference. Most of the existing robotic needle guides are not directly adaptable to this clinical setting due to their weight, space requirements, ease of use, and adaptability to the conventional brachytherapy apparatus e.g. the EXII R stepper and Micro-Touch ° stand. The design by Fichtinger et al. [30] does not suffer from these limitations and has a light weight. However, if power or electronic control failure occurs, the oncologist has to revert back to the conventional needle guide template. Our design aims to avoid such limitations while still being adaptable to clinical apparatus and work flow. Our proposed approach to improve prostate brachytherapy needle placement is to design and implement a four DOF robotic needle guide. The four DOF consist of two translation axes, which mimic the XY (horizontal and vertical) template grid positioning method, and two rotation axes which provide pitch and yaw angulations to the needle. The robot is light and is designed to be compatible with the conventional brachytherapy stepper 13 1.3. Literature Review and Thesis Approach R and mount (the EXII Stepper and the Micro-Touch° stand). The robot’s design also allow manual positioning capabilities in the event of electronic control or power failures. To register the pre-operative prostate volume to the intra-operative volume, we propose a new surface-based registration technique using ellipsoid fitting. Surface-based registration is suitable for the prostate as it is a soft organ without any distinct landmarks that can be used as references. Our approach is to contour the pre-operative and intra-operative prostate volumes using a semi automatic algorithm. The contours are fitted into ellipsoids and their parameters are used to calculate the registration matrix. There are two reasons we are using an ellipsoid as the model of the prostate. First, it is common to use an ellipsoid to calculate the volume of the prostate [49]. Second, ellipsoid fitting can be solved with robust, fast, and computationally inexpensive algorithm such as [50]. In our registration approach, we opt to use the sagittal images instead of the transverse ones due to several reasons. First, our research group has developed an automatic US volume data acquisition system by motorizing the TRUS probe rotation. With this system, the prostate volume can be finely sampled using the sagittal images as the system allows rotation angle increments as low as 0.5o . This angle interval translates to the worst case scenario of 0.52 mm lateral resolution with 60 mm imaging depth. As a comparison, the conventional technique using the EXII stepper samples the prostate with 5 mm increments. Second, the rotation of the TRUS probe does not change its position with respect to the prostate. So, prostate motion due to probe angular rotation should be small. Third, the sagittal images provide direct depiction of pitch rotation motion [14], which has been identified as a significant motion in EBRT procedures [45]. For our image segmentation purposes, we opt to use the semi-automatic segmentation algorithm developed by Badiei et al. [51]. The algorithm ensures smooth and continuous contours, only requires the user to input five initialization points, and does not require training models like [52]. It only takes about 1 second to segment one US image. However, the algorithm by Badiei et al. is intended for prostate segmentation in its transverse images. Modifications have to be made to apply the algorithm for sagittal US images. A survey of available segmentation algorithms for US images can be found in [53] 14 1.4. Thesis Organization 1.4 Thesis Organization This thesis is divided into five chapters. Chapter 1 explains the background knowledge, motivation, and approach of the thesis. In the chapter we give an overview of the general prostate brachytherapy procedure. We also identify some of its key challenges that motivate the thesis project to develop a robotic needle guide and a prostate volume registration algorithm. Based on related publications, we explain the general approach of this thesis. Chapter 2 centers around the robot design and implementation. First, we explain the robot design: its mechanical components, its controller, and its Graphical User Interface (GUI). Then, we performed several experiments to quantify the characteristics of the robot. The chapter also explains a seed implant study conducted on a prostate phantom using the robot prototype. The chapter finishes with updates and changes to the robot since its first prototype. Chapter 3 explains the sagittal view based segmentation algorithm whose contour results are used for registration. The first three sections of the chapter explain the origins of the algorithm [51] and the modifications made to develop the algorithm. We extend the algorithm, in the fourth section of the chapter, to segment 3D-sagittal-view-based US data. In the last section, we validate the algorithm by comparing it with a recently published transverse view based segmentation algorithm [54]. Chapter 4 describes our registration work: simulation, experiment setup, and results. In the chapter’s first section, we explain the idea of using ellipsoid fitting on prostate contours to calculate the registration matrix. Then, we simulate the ellipsoid-fitting algorithm using noiseless but partial data (the second section of the chapter). To measure the registration accuracy, we designed and implemented an experimental setup using an accurate optical tracking as the gold standard - explained in the third section of the chapter. We analyzed and discussed the results in the fourth section of the chapter. At the last section, we explored the use of ICP as an alternative registration algorithm. Chapter 5 concludes the thesis with a summary, list of contributions, and possible future work. 15 Chapter 2 Brachyguide: a Robotic Needle Guide for Prostate Brachytherapy This chapter is divided into 4 sections. Section 2.1 explains the robotic needle guide design and implementation. Section 2.2 explains the characterization tests done on the robot and their results. Section 2.3 describes the first use of the robot during an implant study. 2.1 Robot Design and Implementation We designed and implemented a four DOF serial manipulator (Brachyguide) that addresses the needs of brachytherapy: 1. It facilitates angled needle insertion 2. It has a translational work space of at least 60 × 60 mm2 3. Its positioning resolution is less than 5 mm 4. It is light weight and easily mounted on the commonly used EXII stepper 5. It does not interfere with TRUS insertion or view of the surgical space 6. It provides manual override during power or electronic control failure The translational stage positions the robot wrist in the X-Y direction. The wrist carries the needle guide (see Figure 2.1). The robot wrist provide pitch and yaw angulation to the needle guide. The Brachyguide’s translational stages are made of two A15 UniSlide linear drives (Velmex Inc., Bloomfield, NY). Each drive consists of a low friction lead screw driven slider on a nickel-plated aluminum dovetail base. 16 2.1. Robot Design and Implementation Figure 2.1: Brachyguide kinematics design and implementation. Brachyguide kinematics design. Right: Brachyguide. Left: 17 2.1. Robot Design and Implementation The Y-Axis stage is mounted on the EXII brachytherapy stepper with its slider carrying the X-Axis stage (see Figure 2.2). Each UniSlide has a translation range of 150 mm giving a larger workspace than needed. The M6 lead screws gives 1 mm/ turn resolution to the sliders. They are actuated by 23 mm coreless DC Micromotors (Series 2342-024CR, Faulhaber, Germany) through spur gears with a 3.125 ratio. The motors use a 512 count magnetic encoder (series IE2-512) from Faulhaber. Due to quadrature signal processing, each encoder resolution is actually 2048 counts per turn. Theoretically, the motor - spur gears - lead screw configuration gives 0.65 µm X-Axis and Y-Axis translation resolution. The actual resolution of the translational axes are affected by other factors such as backlash. Figure 2.2: Brachyguide translation stages. The robot’s wrist - mounted on the X-Axis stage - consists of two precision rotation stages (Model 7R174-11, Standa, Vilnius, Lithuania). The stages are mounted orthogonally with the yaw rotation stage carrying the pitch rotation stage (see Figure 2.3). The latter rotates the needle guide through a Kevlar belt parallelogram linkage. A distance of 10 mm separates the center of rotation of the stage and the needle guide. The parallelogram 18 2.1. Robot Design and Implementation linkage limits the pitch angle range within −180 to +240 . While the yaw angle is limited by interference with the translation stages, it still reaches ±300 . The servomotors of the rotation stages are 17 mm DC micromotors (Series 1727-24CR) with IE2-512 encoders from Faulhaber. Each motor drives the rotation stage through a built in worm gear with a 3 degrees/ turn ratio, giving it a theoretical angular resolution of 0.00580 . Figure 2.3: Robot wrist. The needle guide is an acrylic hollow tube (see Figure 2.3). The guide bore is milled with sufficient clearance for an 18 G brachytherapy needle (1.27 mm diameter). The needle guide is “snap mounted” onto the needle holder flexure at the bottom of the parallelogram linkage. The unique feature of the Brachyguide is its manual positioning capability in the case of electronic control failure or power outage. This is achieved through manual cranks and quick release mechanisms on the translation stages, translational readings via digital calipers, and turn dials on the rotation stages (see Figure 2.4). Manual cranks are fixed on the lead screws of the translational stages. Each crank directly rotates the screw, giving fine manual adjustment of 1 mm per turn. However, if the robot’s operator wishes to save time by coarsely adjusting the stages, the quick release mechanism can be used. A spring loaded lever is fixed onto the brass nut of the stage slider. The nut 19 2.1. Robot Design and Implementation Figure 2.4: Brachyguide manual positioning features. Left: rotational stage dial (top) and translation stage manual crank (bottom). Middle: quick release mechanism. Right: digital calipers. 20 2.1. Robot Design and Implementation was custom milled with sufficient clearance at an angle with respect to the M6 lead screw. Thus, the nut can rotate around an axis orthogonal to the slider and lead screw. When the spring lever is pulled, the nut disengages the lead screw allowing the slider (stage) to move freely. Linear digital calipers (Model 303-9303, Shars Tools Inc., IL) gives the XY position of the needle guide. These calipers can be interfaced to the computers through a simple serial connection; and serve as redundant sensors to the linear stages’ motor encoders. In the absence of electronic control, the calipers help the operator position the needle guide (similar to the use of conventional grid template). The LCD displays on the calipers give 0.01 mm reading resolution while the engraved markings give 1 mm reading resolution. Turning dials are installed between the rotation stages and their motors. The dials provide 0.050 angular reading accuracy. We do not expect needle angulations of more than ±150 as it translates to an endpoint position change of tens of mms. One turn of the dial translates to three degrees rotation and with five turns of the dial, one can easily obtain ±150 needle angulation. Hence, quick release mechanisms for the rotation stages are not necessary. 2.1.1 Control and Graphical User Interface (GUI) Each motor is individually controlled using the Faulhaber’s MCDC 3006S microcontroller. Four MCDC 3006S units are daisy-chained together and connected to a PC through a single serial port (see Figure 2.5). Each of the controllers employs PI control system for speed control and PD control system for position control. Controller parameters were set according to the user manual using Faulhaber’s user interface. A medical grade power supply (SL Industries, CA, US) powers up the controllers. 21 2.1. Robot Design and Implementation Figure 2.5: . The block diagram diagram of Brachyguide controller (top). The block diagram for the MCDC 3006S (bottom),taken from its manual [55]. 22 2.2. Brachyguide Characteristics The operator can control the robot’s motion through a simple GUI (Figure 2.6). The interface provides individual control to each of the robot’s axis. The operator can type in the end point X-Y coordinates (in mm) and pitch-yaw angle values (in degrees) and command the stages to move accordingly. As a safety feature, the encoder readings are monitored through multiple threads and the motors are disabled when the desired coordinate is reached. Scaled motor encoder readings and digital caliper readings are displayed as feedback to the operator. The GUI also provides an option to load pre-determined needle plans from VariSeed. Needle plans from VariSeed can be exported into text files and then loaded into the GUI for the robot to follow. These plans do not incorporate angled needle insertions. Figure 2.6: Brachyguide GUI. The needle plan from VariSeed is shown in the center. Individual motor control is shown on the left (green) while calipers readings are shown on the top right (red). 2.2 Brachyguide Characteristics We performed several characteristics tests on the Brachyguide prototype. Table 2.1 summarizes the test results. 23 2.2. Brachyguide Characteristics Mass Translational workspace Pitch Range Yaw Range Translational Accuracy X-axis Max Mean Std. Dev Y-axis Max Mean Std. Dev. Max X-axis Mean Std. Dev Caliper accuracy Max Y-axis Mean Std. Dev Max Rotational accuracy Mean Std. Dev X Axis Stiffness Y Axis Stiffness 60 mm Point to point move time 150 mm 30 degree 2.6 kg 150 × 150 mm2 +240 to -180 +300 to -300 0.292 mm 0.120 mm 0.102 mm 0.89 mm 0.607 mm 0.299 mm 0.059 mm 0.022 mm 0.014 mm 0.066 mm 0.017 mm 0.015 mm 0.042o 0.018o 0.008o 9200 N/m 15500 N/m <2 s <4.6 s 0.7 s Table 2.1: Brachyguide characteristics table. 24 2.2. Brachyguide Characteristics 2.2.1 Translational Accuracy and Straight Motion Deviation X-axis and Y-axis translation accuracy was tested using a laser interferometer (ML10, Renishaw PLC, Gloucesteshire, UK) in a setup shown in Figure 2.7. Each axis’ stage travels unidirectionally while caliper, encoder, and interferometer readings were recorded at 5 mm intervals. Due to reflector interference, the Y-Axis range was limited from 0-110 mm. Figure 2.7: Translational accuracy measurement setup. As seen from Figure 2.8, maximum caliper errors were within 0.06 mm of the interferometer reading with an average error of about 0.02 mm (both directions). Encoder readings shows that the motor control always follow the commanded motion input. However, the maximum actual motion errors are within 0.3 mm (X-Axis) and 0.9 mm (Y-Axis). Motion errors are primarily caused by backlash in the lead screw nuts of the translational stages. The error is more severe on the Y-axis stage because its nut is loaded by the combined weight of the X-axis stage and the robot wrist. Different friction levels at different point of contacts between the slider and the UniSlide’s base caused error variations within a single motion. In this initial work, we controlled the motor without using the digital calipers as redundant sensors. Using the caliper readings in the control loop will help reduce motion errors due to backlash. 25 2.3. Seed Implant Study on a Phantom A digital dial gauge (Mitutoyo Co., UK) was mounted onto the needle holder to measure straight motion deviation. The dial gauge’s arm was rested against a granite table while the stage traverse its full motion range (0-150 mm) back and forth. A deviation of less than 0.1 mm was recorded. 2.2.2 Angular Accuracy Angular accuracy was measured using the interferometer which has a limited motion range within ±1.50 (see Figure 2.9). The maximum angular error was 0.0420 , the average error was 0.0180 and standard deviation was 0.0080 . The result confirms the 3 arcmin stage resolution specification (Standa, Vilnius, Lithuania). The graph showed several random noises (spikes) which we suspected was due to vibrations on the optical table. 2.2.3 Stiffness We determined the stiffness by pressing a load cell (MDB-2.5, Transducer techniques, Temecula, CA) against the Brachyguide’s wrist and measuring its deflection using the interferometer (see Figure 2.10 for illustration). The force data was plotted against the measured displacement and a straight line is fitted onto the data (see Figure 2.11). The slope of the line gives the measured stiffness. Both X and Y axes stiffnesses were measured. In a similar fashion, we measured the stiffness of the Brachystand (EXII stepper R mount). mounted on the MicroTouch° The Brachyguide’x X-Axis and Y-Axis measured stiffness are 9200 N/m and 15500 N/m while the Brachystand’s are 2600 N/m and 6500 N/m respectively. The Brachyguide frame is significantly stiffer than the Brachystand. Consequently the accuracy of the needle placement will mostly be affected by flexing of the Brachystand when the Brachyguide is mounted on it. This flexing comes mainly from the base of the stepper: the connection point between the EXII stepper and the MicroTouch stand. The problem can be reduced by adding a support mechanism to the stepper which will be used in conjunction with the Brachyguide. 2.3 Seed Implant Study on a Phantom We conducted a seed implant study [56] on a prostate phantom (Model 053, CIRS Inc., Norfolk, Virginia) using a setup consisting of: the Brachyguide R protoype, the EXII Stepper, the MicroTouch° stand, and an US machine (Sonix RP, Ultrasonix Medical Corp., Richmond, BC). The purpose of this 26 2.3. Seed Implant Study on a Phantom Figure 2.8: Translational accuracy measurement results. X+ direction is the direction away from the ML10 laser and Y+ direction is the upward direction. 27 2.3. Seed Implant Study on a Phantom Figure 2.9: Angular accuracy measurement measurement results. 28 2.3. Seed Implant Study on a Phantom Figure 2.10: Brachyguide stiffness measurement setup. 29 2.3. Seed Implant Study on a Phantom Figure 2.11: Brachyguide stiffness measurement results. The colored lines (red and blue) represent the linear fit to the data. 30 2.3. Seed Implant Study on a Phantom implant study was to test whether or not the robot use disrupts the clinical work flow or adds to implant procedure time. Several days before the implant, a medical physicist conducted a volume study on the phantom as typically done in the B.C. Cancer Agency (BCCA). Using VariSeed, a seed implant plan consisting 136 dummy seeds in 26 needles were generated. Both loose and stranded seeds (RAPIDStrand, Oncura, Plymouth Meeting, PA) were used. The implant plan was loaded into the Brachyguide’s motion controller Before the actual implant, the Brachyguide was registered to the US image by using a water bath. The implant grid was overlaid onto the BMode transversal images. The robot was used to guide a needle onto the center of the grid. Once the needle tip was located, the encoder and caliper readings were zero-ed. During the seed implant study, a medical physics student read out the needle position from the printed plan, an electrical engineering student operated the robot, and the radiation oncologist performed the needle insertion into the phantom (see Figure 2.12). The study was finished in 32 minutes: comparable to conventional procedure time (30-40 minutes). The study showed that Brachyguide use did not alter implant procedure time although it was the first time use of the prototype by the oncologist. Figure 2.12: Brachytherapy seed implant study on a phantom. Guided by the US image, the radiation oncologist is inserting the needle using the Brachyguide. A post implant CT scan (in-slice resolution 0.0293 mm/pixel, thickness 31 2.3. Seed Implant Study on a Phantom Maximum absolute error (mm) X-axis 3.3301 Y-axis 2.0464 Distance 3.4273 in XY plane Average error (mm) 0.108 0.1993 1.2296 Standard deviation of errors (mm) 1.1048 0.8689 0.7102 Table 2.2: Seed placement errors based on post implant CT scan. 1.5 mm) was performed on the phantom. Xu Wen, a colleague from the UBC Robotics and Control Lab, applied his seed detection algorithm [57] to quantify the seed placement errors (see Table 2.2) measured relative to the planned seed location. Figure 2.13 shows the seed placement error histograms. As illustrated in Figure 2.14, the implanted seeds are projected onto a plane orthogonal to the needles direction and then compared with the pre-plan seeds locations. Although the results were affected by registration error, backlash, and needle deflection and deformation, the errors are small, ranging from 2 - 3.5 mm. In this study, we did not perform needle insertion with oblique trajectories since they have not been incorporated within the general clinical procedures. 32 2.3. Seed Implant Study on a Phantom Figure 2.13: Seed placement errors histograms. Figure 2.14: Projection of the CT scan of the prostate phantom onto the plane orthogonal to the planned needle insertions. 33 2.4. Additional Features and Updates on the Robot 2.4 Additional Features and Updates on the Robot Since its first prototype trial, changes and additional features have been added to the robot with the help of Victor Wang and Honza Vrodazky. First, new linear digital calipers were installed into the robot. It was found that whenever the calipers from Shars Tools Inc. are interfaced with the computer, their measurement values are noisy (they deviate by 0.04 mm under stationary positions). The noise is caused by a design problem within the caliper serial interface unit. Thus, we replaced the calipers with OEM scales from Fowler (Fred V. Fowler Co. Inc., Newton, MA) as shown in Figure 2.15. Figure 2.15: The robots calipers were replaced by ones from Fowler Second, control loop incorporating the caliper readings was programmed into the robot to compensate for backlash on the translational stages. When each of the robot’s translational stages is actuated, the controller samples the encoder reading at 500 ms intervals. Each time, the controller calculates the distance difference between the target position and the current position. Once the difference is within 0.01 mm, the controller switches to sampling the caliper’s reading at 600 ms intervals. The stage is moved until the distance difference between the target position and the caliper’s reading are within 0.03 mm. Unlike the encoder readings, the caliper’s readings are affected by the stage’s backlash. Hence, the switch from encoder to caliper readings provides backlash adjustments. Last, Victor Wang added a path-finding module and updated the GUI to incorporate angled needle plans. The path-finding module helps the robot 34 2.4. Additional Features and Updates on the Robot avoid collision with the TRUS probe as it moves to its designated position. The details of these improvements are not described in this thesis. Victor Wang repeated the translational accuracy measurement using the same setup in Figure 2.7. Table 2.3 and Figure 2.16 show the updated X-axis and Y-axis accuracy measurement and their repeatability. The Renishaw Laser 10 linear measurement software (ML10, Renishaw PLC, Gloucesteshire, UK) was used to plot the measurement results and to calculate the repeatability values. Figure 2.16 shows four plots in each box. The blue-cyan and red-green pairs describe two sets of uni-directional accuracy error measurements. The error values are the differences between the linear encoder’s readings and the laser interferometer measurements. With the new calipers and an improved control system, we achieved less than 0.12 mm and less than 0.1 mm X-axis and Y-axis translational accuracies. Such accuracies cannot be achieved using the current template grid approach as it limits the X-Y spacing to 5 mm. Table 2.3: Updated translational accuracy and repeatability measurement results. The error values in the table show the differences between linear caliper readings and laser interferometer measurements. X Axis Y Axis Property Value Accuracy 113.6 µm Positive Direction Repeatability 13.1 µm Negative Direction Repeatability 11.8 µm Bi-directional Repeatability 70.2 µm Accuracy 93.1 µm Positive Direction Repeatability 8.8 µm Negative Direction Repeatability 8.1 µm Bi-directional Repeatability 39.2 µm 35 2.4. Additional Features and Updates on the Robot Figure 2.16: ments. Updated translational accuracy and repeatability measure- 36 Chapter 3 Segmentation Based on Sagittal Images The sagittal view-based segmentation algorithm (“sagittal segmentation algorithm”), was developed by modifying Badiei et al.’s algorithm (“transversal segmentation algorithm”) [51]. The transversal segmentation algorithm outlines the prostate boundaries in 2D transverse B-mode transrectal ultrasound (TRUS) images based on image warping, ellipse fitting, and the IMMPDAF (Interacting Multiple Model Probabilistic Data Association Filter) edge detector [58]. This chapter starts by introducing the transversal segmentation algorithm in Section 3.1. A brief explanation of the algorithm steps will help the readers follow the modifications explained in Section 3.2. These modifications convert the transversal segmentation algorithm to the sagittal one (Section 3.3). Section 3.4 explains the extended sagittal segmentation algorithm. Section 3.5 gives validation of the algorithm in phantom and patient volume data. 3.1 The Transversal Segmentation Algorithm The transversal segmentation algorithm consists of the following steps (see Figure 3.1 for illustrations of each steps) 1. User initialization 2. Image warping 3. Initial ellipse fit 4. IMMPDAF edge detection 5. Symmetrization and final ellipse fit 6. Image unwarping 37 3.1. The Transversal Segmentation Algorithm Figure 3.1: Six steps of the transversal segmentation algorithm by Badiei et al. The algorithm is applied on a transverse image of the CIRS phantom. 3.1.1 User Initialization The transversal segmentation algorithm asks the user to input five initialization points at specific locations on the prostate’s midgland transverse view as illustrated in Figure 3.3(a). The order of the initialization point is predetermined. The first initialization point is the TRUS probe center. The second one is the center of the prostate. The last three points are located at the prostate’s right most edge, mid-bottom, and top. These points are used to calculate the warping factor in step 2. Using 5 initialization points, 3 additional points along the prostate boundary are interpolated. It is assumed that the prostate boundary is symmetric along the prostate center (initialization point 2). 3.1.2 Image Warping A novel warping function transforms the shape of the prostate to resemble that of an ellipse. It converts the pixel coordinate system (u, v) into a polar coordinate system (R, θ) with the TRUS center as the coordinate center. A warping function that is Gaussian in the radial direction and sinusoidal in 38 3.1. The Transversal Segmentation Algorithm the angular direction is then applied: −R2 Rwarped = −R · [1 − sin θ · e 2σ2 ] (3.1) The standard deviation (σ) is calculated using the coordinates of the initialization points (the TRUS location, the seed location, the top initialization point location and, the bottom initialization point location) with formula 3.2 (see Figure 3.2). s −(vT − v3 )2 σ= (3.2) 3 − 2vs 2 · ln(1 − vT +vTv−v ) 3 Figure 3.2: Warping factor calculation. The transverse image was generated using the CIRS phantom. The warping function is applied to the image pixels, located radially as far as 356 pixels. The warping is also applied to the initialization points, but excludes the probe center and the prostate center. 39 3.2. Modifications on the Transversal Segmentation Algorithm 3.1.3 Initial Ellipse Fit An ellipse fitting algorithm [59] is applied to the warped initialization points. The detailed mathematical background of the algorithm is provided in Appendix B.1. Using the ellipse parameter results, 120 ellipse sample points are drawn on the warped image. The ellipse fit ensures smooth contours of the prostate. 3.1.4 IMMPDAF Edge Detection The ellipse sample points are used to guide the IMMPDAF edge detector. The full mathematical explanation of the detector is provided in Appendix C. 3.1.5 Symmetrization and Final Ellipse Fit The algorithm enforces symmetry on the IMMPDAF result by mirroring edge points about the prostate center line (a vertical line crossing through the prostate center). After the symmetrization process, a second ellipse fit is applied to the boundary points and 120 ellipse sample points are generated from the ellipse parameters. 3.1.6 Image Unwarping The last step of the transversal segmentation algorithm is the unwarping of final ellipse sample points. The unwarped ellipse sample points become the final segmentation result. All the ellipse sample points are expressed in the same polar coordinate system used in the warping step. The radius of each sample point corresponds to Rwarped and we wish to find the unwarped radius (R) - see equation (3.1). Since we found no closed form solution, we applied the Newton Rhapson method with R = 356 pixels as the initial guess. 3.2 Modifications on the Transversal Segmentation Algorithm The following subsections explain the modifications on the transversal segmentation algorithm that leads to the sagittal segmentation algorithm. 40 3.2. Modifications on the Transversal Segmentation Algorithm 3.2.1 Initialization Points We have modified the requirements of the initialization points. In the sagittal segmentation algorithm, we use the midsagittal US image of the prostate. We ask the user to input the first point at the anterior midgland of the prostate. The second point is at either the base or the apex of the prostate (depending on which part is more visible). The third point is at the posterior midgland of the prostate. We give the flexibility of putting the last two points along the boundaries of the prostate. In the transversal segmentation algorithm, the seed point is defined as the center of the prostate and is the second initialization point. This point is important for the calculation of the warping factor (vs in equation 3.2). We decided to remove the seed point from the user initialization for sagittal segmentation because sometimes it is hard to determine the prostate center. For example, this occurs when only part of the prostate is visible, or when the prostate is rotated. Therefore, the seed location is calculated using the first three initialization points: ( u1 +2 u3 , v1 + 2v42 + v3 ) where u refers to the X-axis image coordinate and v refers to the Y-axis image coordinate (as shown in Fig 3.2). The first three initialization points for the sagittal viewbased segmentation algorithm is as shown in Figure 3.3(b). We expect that the user will place the second initialization point around the same height as the prostate’s center. Hence, more weight is given on the coordinate v2 . By applying symmetry around the seed location, three additional points are interpolated, also illustrated in Figure 3.3(b). The TRUS probe center is not part of the user initialization points as it can always be calculated by using the probe’s radius and US pixel conversion ratio. The TRUS probe center is at the bottom left corner of the image, opposite the image origin (the top left corner), unlike in the transverse images. The probe that we use has a 9.87 mm radius [60]. In our setting, the imaging depth is set to 60 mm with 156 µm per pixel conversion ratio. If the image settings change, the location of the probe can always be recalculated. 3.2.2 Warping The sinusoidal component of the warping function in the sagittal case is dependent on the rotation of the probe around its center axis (as illustrated in Figure 3.4). Therefore, we removed the sinusoidal dependency from equation 3.2 and substituted the probe rotation angle as θ into equation 3.1. Following equation 3.2, the warping factor (σ) is calculated using the seed point location, the TRUS location, the top initialization point (first initialization 41 3.2. Modifications on the Transversal Segmentation Algorithm (a) Transverse midgland plane of the CIRS phantom. (b) Midsagittal plane of the CIRS phantom. Figure 3.3: Modifications in initialization points. Figure (a) and (b) show the initialization points of the transversal and sagittal segmentation algorithms. The red dots are the user initiated points. The green dots are the interpolated points. 42 3.2. Modifications on the Transversal Segmentation Algorithm point), and the bottom initialization point (third initialization point). The warping function is applied uniformly across the width of the image. It is only applied to the lower half of the image. This approach makes the warping step faster without sacrificing the accuracy because warping at the top half is negligible compared to the pixel size (see Figure 3.5 for illustration). Figure 3.4: Image warping on sagittal images. The images show the consideration for modifications. 43 3.2. Modifications on the Transversal Segmentation Algorithm Figure 3.5: Uniform warping (radial compression) across the bottom half of the image. 44 3.2. Modifications on the Transversal Segmentation Algorithm 3.2.3 Selective IMMPDAF Edge Detector On occasion, the IMMPDAF edge detector encounters a division by zero problem. This happens when the detector strays outside the image boundary due to the prostate being larger than the US array width as in our CIRS phantom images. Division by zero also happens when the detector encounters a completely black prostate region (image pixel magnitude of zero). The problem is caused by the PDAF (Probabilistic Data Association Filter), refer to Appendix C.2.3 for details. The PDAF collects 5 candidate edge points with the largest edge magnitudes. It assigns a Gaussian probability density function which depends on the edge magnitude to each of the candidates. When the image region is completely black, the five candidate edge points have edge magnitudes of zero, causing divisions by zero in the probability density calculation. To solve this problem the IMMPDAF detector has been made selective: whenever all the five candidate edge points have zero edge magnitudes, the detector skips to the next Kalman filter step. We also have removed the symmetrization step from the transversal segmentation algorithm as we do not expect the prostate to be symmetric in the sagittal view. 3.2.4 Stable Ellipse Fitting Algorithm A more stable ellipse fitting algorithm [61] is employed in the segmentation algorithm. The full mathematical formulation of this algorithm is provided in Appendix A.2. The ellipse fitting problem comes down to finding the eigenvalues and eigenvectors pair of the following mathematical formulation: → − → − SP = λ CP (3.3) where S is the Scatter Matrix, C is the Constraint Matrix, and P is the parameter vector (a six element vector containing the algebraic parameter of the ellipse). The Scatter Matrix and the Constraint Matrix are formulated as follows: x21 x22 x1 y1 x2 y2 2 y1 y22 S = x1 x2 y1 y2 1 1 ··· ··· ··· ··· ··· ··· x2n x21 x1 y1 y12 x1 y1 1 x n yn 2 x2 x2 y2 y22 x2 y2 1 yn2 .. .. .. .. .. .. xn . . . . . . yn x2 xn yn y 2 xn yn 1 n n 1 (3.4) 45 3.3. The Sagittal View-based Segmentation Algorithm Where (xn , yn ) are the coordinate sample points to be fit into an ellipse. 0 0 2 0 0 0 0 −1 0 0 0 0 2 0 0 0 0 0 C = (3.5) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 As noted by Halir and Flusser [61], the Constraint Matrix is singular and when the set of points are very close to an ellipse, the Scatter Matrix is close to singular: in occasions, leading to the wrong sets of eigenvalues and eigenvectors. They circumvent the problem by decomposing the Scatter matrix, and the Constraint matrix: · S1 S2T S2 S3 ¸ · P1 P2 ¸ · = λ C1 0 0 0 ¸ · P1 P2 ¸ (3.6) The solution to ellipsoid fitting boils down to solving the following equations and eigensystem: S1 P1 + S2 P2 S2T P1 C1−1 (S1 − = λC1 P1 (3.7) + S3 P2 = 0 S2 S3−1 S2T ) P2 = −S3−1 S2T P1 (3.8) P1 = λP1 (3.9) Since C1 and (S1 − S2 S3−1 S2T ) are no longer singular, the eigensystem is more numerically stable. 3.3 The Sagittal View-based Segmentation Algorithm The overall flow of the sagittal segmentation algorithm is illustrated in Figure 3.6. The algorithm steps are as follows: 1. Input 5 initialization points 2. User accepts or declines the extrapolated initialization points. If the user accepts, move to next step. Otherwise, pick 5 new initialization points 46 3.3. The Sagittal View-based Segmentation Algorithm 3. Warp the image and initialization points 4. Fit an ellipse to the warped initialization points 5. Apply the selective IMMPDAF edge detector guided by the initial ellipse fit 6. Fit a final ellipse fit to the detector results 7. Unwarp the image and final ellipse fit 47 3.3. The Sagittal View-based Segmentation Algorithm Figure 3.6: The steps of the sagittal segmentation algorithm. 48 3.4. Extending the Sagittal Segmentation Algorithm 3.4 Extending the Sagittal Segmentation Algorithm By controlling the TRUS probe rotation, we can collect a set of 2D B-mode sagittal images with finely sampled rotation angle intervals. The collected images make up a 3D US volume and can be visualized using software such as Stradwin (an open source Windows based research software developed from Stradx [62]). It offers some flexibility in viewing the prostate in 3D. We can use the sagittal segmentation algorithm to contour the prostate boundary from each individual 2D sagittal image in the acquired volume. However, it is very tedious and time consuming to segment the images one by one as the user must input initialization points at each image. For example, a finely sampled volume data (e.g. 0.5o rotation angle intervals) typically consists of 200 images. We developed a method to segment all the sagittal images automatically by using the ellipse slice-generation method described in Appendix B.3 and Stradwin. The result is an extended segmentation method based on sagittal image-based volume data. This method is hereafter called the “extended sagittal segmentation algorithm”. The algorithm steps are as follows: 1. The volume data are exported into Stradwin to be visualized in 3D. 2. Using Stradwin, 14 initialization points are selected (seven on the mid sagittal view and seven on one of the coronal views, see Figure 3.7 for illustration). 3. The initialization points coordinates are read by the extended sagittal segmentation algorithm. 4. The warping factor is calculated using the first three initialization points in the mid sagittal view. 5. Warping is applied to the initialization points and an ellipsoid is fit into the warped initialization points. 6. With the known probe rotation angles, ellipse slices are generated from the ellipsoid. 7. Each ellipse slice becomes the initial ellipse fit for the sagittal segmentation algorithm (step 4 in Section 3.3). 8. Each sagittal images in the volume is segmented using step 5 - 7 of the sagittal segmentation algorithm (see Section 3.3). 49 3.4. Extending the Sagittal Segmentation Algorithm Figure 3.7: In the extended sagittal segmentation algorithm, fourteen initialization points are selected from the mid-sagittal view and one of the coronal views. 50 3.5. Segmentation Results on Phantom and Patient Volume Data In the extended sagittal segmentation algorithm, Stradwin is used to visualize the volume data and to choose the fourteen initialization points which are required for warping and ellipsoid fitting. These fourteen points are selected based on the mid sagittal image and one of the coronal view images. The coronal view images are interpolated by Stradwin based on the US volume data. We choose the coronal view with the clearest prostate image. Instead of one of the coronal images, we can opt to use one of the transverse images. However, in our phantom images, the coronal view shows sharper prostate boundaries compared to the transverse ones. In summary, we extended our 2D sagittal segmentation algorithm to be able to segment sagittal image-based US volume data. The extension requires a larger number of initialization points (14 instead of 5) and the Stradwin 3D visualization software. However, it offers the convenience of not having to initialize each individual image. 3.5 3.5.1 Segmentation Results on Phantom and Patient Volume Data Segmentation on CIRS Phantom Volume Data We tested the extended sagittal segmentation algorithm on 25 sagittal viewbased volume data sets of the CIRS phantom. The volume data was collected under the following conditions: imaging depth of 60 mm, 100o sector angle, and 6.6 MHz center frequency. These settings translate to horizontal and vertical conversion ratios of 156 µm per pixel. We used the Ultrasonix US machine (Ultrasonix, Richmond, BC,CA) with the BIP 6.5/R10/2x128-472 TRUS probe (Vermon, Tours, France) [60]. Note, we segmented the same volume data set that is used for our registration experiment (explained in Section 4.3). Figure 3.8 shows the Stradwin visualization of the sagittal view-based segmentation results of two CIRS phantom volumes. 51 3.5. Segmentation Results on Phantom and Patient Volume Data Figure 3.8: Extended sagittal segmentation results on two CIRS phantom volume data as visualized using Stradwin. Left: mid sagittal view, center: mid transverse view, right: one of the coronal views of the. The red lines on the transverse and coronal view show the location of the mid sagittal plane. 52 3.5. Segmentation Results on Phantom and Patient Volume Data In the absence of a gold standard for the true segmentation, we compared the results with the transverse image-based segmentation results obtained using the 3D segmentation algorithm by Mahdavi et al. [54]. First, the sagittal view-based US volume data and the extended sagittal algorithm contour results were re-sampled in the transverse views. This step involved scan converting the volume data to produce sector transverse images views and linearly interpolating the sagittal contour points along the transverse axis direction (TRUS probe axis direction). Second, we manually determined the base of the phantom from the resampled volume data (the transverse view-based US volume data). From the transverse view-based volume data, we extracted several transverse images. The images were sampled using 5 mm step intervals along the transverse view axis (TRUS probe axis) starting from the base. Last, we segmented the extracted transverse images using the recently published segmentation algorithm by Mahdavi et al. [54]. The transversal segmentation results were then compared with the interpolated sagittal segmentation results. The comparison results were affected by the re-sampling of the transverse view-based volume from the original sagittal view-based volume. In each transverse image, we calculated the mean absolute distance (MAD) and maximum absolute distance (MAXD) between the transversal segmentation contour points and the interpolated sagittal segmentation contour points. The distance is defined based on the point to line distance calculation (see Figure 3.9 for illustration). The green curve in the figure illustrates the transversal segmentation results while the red one illustrates the interpolated sagittal segmentation results. For each transversal segmentation results, there are 120 contour points in total. In our comparison calculations, we only use three transverse images for each volume data: the images at base+20mm, base+25mm (midgland), and base+30mm. The reason is the phantom size is larger than the region of interest in the sagittal view. Hence, the apex of the phantom is not visible in the sagittal US images. The transversal segmentation algorithm requires full view of the prostate (from base to apex) because it fits the contours into a tapered ellipsoid. The partial view of the phantom moves the apex forward, giving an impression of a smaller prostate. Hence, the segmentation results are only reliable near the midgland. Table 3.1 and Figure 3.10 show the comparison between the segmentation results. The overall maximum absolute distance is 2.01 ± 0.79 mm while the mean average distance is 0.72 ± 0.37 mm. An interesting observation is the jagged edges on the interpolated sagittal segmentation contours are due 53 3.5. Segmentation Results on Phantom and Patient Volume Data Figure 3.9: Distance calculation between two segmentation results. The red curve illustrates the interpolated sagittal segmentation results while the green one illustrates the transversal segmentation results. to interpolation. The contours are sampled finely in the sagittal view but not in the transverse view. Even with these jagged edges, the mean average distance is still relatively small: 0.72 ± 0.37 mm. Excluding the user’s time for initialization, the extended segmentation algorithm can contour a volume data of 200 ± 0.5 US images in 73 ± 3 seconds. In other words, it takes about 0.37 ± 0.01 seconds to segment each individual US sagittal image. The algorithm was implemented in Matlab running on a 1.66GHZ computer with 1 GB RAM. 54 3.5. Segmentation Results on Phantom and Patient Volume Data Figure 3.10: Comparison of interpolated sagittal segmentation results (red) and transverse segmentation results (green). Notice that the green line is smooth while the red line has jagged edges. 55 3.5. Segmentation Results on Phantom and Patient Volume Data Table 3.1: Validation of segmentation results on CIRS phantom volume data. Base+20mm Base+25mm (Midgland) Base+30mm MAXD (mm) MAD (mm) MAXD (mm) MAD (mm) MAXD (mm) MAD (mm) Vol data 1 1.84 0.51 1.29 0.36 0.90 0.34 Vol data 2 2.21 0.65 2.14 0.49 1.62 0.50 Vol data 3 1.49 0.71 1.48 0.57 1.27 0.44 Vol data 4 2.10 0.95 1.95 0.76 1.88 0.70 Vol data 5 2.47 0.51 2.12 0.36 2.68 0.46 Vol data 6 1.87 0.47 1.74 0.47 1.71 0.61 Vol data 7 3.47 0.63 0.99 0.41 0.90 0.30 Vol data 8 1.46 0.39 1.76 0.47 1.56 0.63 Vol data 9 1.55 0.64 0.98 0.39 0.76 0.24 Vol data 10 2.44 0.97 2.44 1.07 2.49 1.26 Vol data 11 2.86 0.88 2.32 0.70 1.84 0.58 Vol data 12 3.89 1.69 3.53 1.45 3.04 1.34 Vol data 13 2.18 0.71 1.59 0.54 2.42 0.55 Vol data 14 5.19 1.64 3.15 1.28 2.20 0.87 Vol data 15 1.50 0.47 1.44 0.52 1.47 0.55 Vol data 16 1.59 0.54 1.65 0.54 1.60 0.50 Vol data 17 2.20 0.91 1.67 0.70 1.19 0.56 Vol data 18 1.26 0.43 1.21 0.44 1.11 0.39 Vol data 19 1.79 0.55 1.51 0.41 0.99 0.37 Vol data 20 3.97 1.69 3.07 1.43 2.45 1.23 Vol data 21 2.50 1.01 2.58 0.84 2.32 0.75 Vol data 22 2.79 1.68 2.56 1.35 2.14 0.96 Vol data 23 2.48 0.94 1.97 0.73 1.74 0.62 Vol data 24 1.27 0.51 1.86 0.52 2.40 0.98 Vol data 25 1.97 0.56 1.17 0.38 1.24 0.36 Average (mm) 2.33 0.83 1.93 0.69 1.76 0.64 Std. Dev (mm) 0.95 0.42 0.67 0.35 0.63 0.30 Overall MAXD MAD Average (mm) 0.72 2.01 Std. Dev (mm) 0.79 0.37 3.5.2 Segmentation on Patient Volume Data We also tested the extended sagittal segmentation algorithm on sagittal view-based volume data of 6 patients. The volume data were collected during an ongoing elastography study within our group by Mehdi Moradi and Xu Wen in the BCCA. They comprised of pre-processed B-mode sagittal images captured at 0.4o probe rotation angle increments. The volume data were collected under the following conditions: imaging depth of 50 mm, 100o sector angle, and 5.0 MHz center frequency. These settings translate to horizontal and vertical conversion ratios of 175 µm per pixel. We used the Ultrasonix US machine (Ultrasonix, Richmond, BC,CA) with the BIP 6.5/R10/2x128-472 TRUS probe (Vermon, Tours, France) [60]. The ethics approval certificates for the use of these patient data can be found in Appendix I. 56 3.5. Segmentation Results on Phantom and Patient Volume Data Prior to segmentation, low level image processing was applied to each pre-processed B-mode sagittal images. The processing consisted of image enlargement and interpolation. The enlargement was needed to maintain the uniform horizontal and vertical pixel ratio (square pixels). Figure 3.11 shows the Stradwin visualization of the sagittal view-based segmentation results of two patient volume data. We compared the segmentation results with transverse image-based segmentation results obtained using the 3D segmentation algorithm by Mahdavi et al. [54] - similar to the steps for the phantom data. The only difference is the method of extracting the transverse images. To extract the transverse images, first we manually chose the location of both the base and the apex from the resampled volume data (the transverse view-based US volume data). Once these two have been identified, we extracted 9 transverse images in between the base and the apex with equal spacing. Hence, the transversal segmentation algorithm was applied to a set of 11 transverse images for each patient. However, we chose to exclude the base and apex segmentation results when calculating the distances between the transversal segmentation contour points and the interpolated sagittal segmentation contour points. Due to poor visibility of the base and the apex of the prostate, we do not include their contours in our comparison. Table 3.2 and Figure 3.12 show the comparison between the segmentation results. The overall maximum absolute distance is 3.54 ± 1.17 mm while the mean average distance is 1.05 ± 0.31 mm. From the transverse images of patient 2 in Figure 3.12 and Figure 3.11, we can see the short coming of the extended sagittal segmentation algorithm. When the prostate is significantly “squished” by the probe, the contours at the bottom left and bottom right tailings of the prostate are slightly inaccurate. This is due to the probe “squishing” effect that lowers the rightmost and leftmost prostate edges. If the coronal view is chosen above these edges, the initial ellipsoid fit will give a smaller ellipse fit that does not cover the bottom right and bottom left tailings which affects the final contour. In comparison, the transverse view usually gives a better view of the probe “squishing” effect hence the transverse images-based segmentation tends to give better results around the bottom right or bottom left of the prostate when the “squishing” effect is significant (as illustrated in Figure 3.12). Excluding the user initialization time, the extended segmentation algorithm can contour a volume data of 234 ± 46 US images in 112 ± 24 seconds. In other words, it takes about 0.5 ± 0.02 seconds to segment each individual US sagittal image. The algorithm was implemented in Matlab running on a 1.66 GHz computer with 1 GB RAM. 57 3.5. Segmentation Results on Phantom and Patient Volume Data Figure 3.11: Extended sagittal segmentation results on two patient volumes (top row and bottom row) as visualized using Stradwin. Left images: mid sagittal view, center images: mid transverse view, right images: one of the coronal views. The red lines on the transverse and coronal view show the location of the mid sagittal plane. 58 Table 3.2: Validation of segmentation results on patient volume data. Base 5 mm 10 mm 15 mm 20 mm 25 mm 30 mm 35 mm 40 mm 45 mm Mean (mm) Std (mm) Overall Average Std. Dev Patient Data 1 MAXD MAD Patient Data 2 MAXD MAD Patient Data 3 MAXD MAD Patient Data 4 MAXD MAD Patient Data 5 MAXD MAD Patient Data 6 MAXD MAD (mm) 3.15 3.45 3.11 3.08 3.21 3.12 3.22 3.77 4.06 3.35 0.34 (mm) (mm) 2.09 1.01 2.28 0.86 2.23 0.69 2.85 0.65 3.23 0.77 2.93 0.96 3.31 1.22 3.69 1.66 5.70 1.94 3.15 1.08 1.10 0.45 MAXD (mm) 3.54 1.05 (mm) 4.89 3.33 3.09 6.12 6.58 5.89 3.72 3.85 4.06 4.62 1.30 (mm) 3.77 2.88 2.51 3.56 3.54 2.23 2.50 3.33 6.11 3.38 1.16 (mm) (mm) 2.33 1.07 3.14 1.04 3.02 1.07 2.30 1.00 2.43 1.05 2.80 1.12 3.29 1.23 3.92 1.27 4.71 1.34 3.10 1.13 0.80 0.12 MAD (mm) 1.17 0.31 (mm) 4.42 2.30 2.82 3.21 3.94 3.99 3.66 3.59 4.58 3.61 0.74 (mm) 1.32 1.25 1.27 1.16 1.08 0.98 1.02 1.12 1.62 1.20 0.19 (mm) 2.31 1.23 1.04 1.21 1.47 1.34 1.27 1.38 1.44 1.41 0.36 (mm) 1.37 1.20 1.08 0.94 0.94 0.76 0.65 0.72 1.34 1.00 0.27 (mm) 1.79 0.87 0.80 0.93 1.17 1.27 1.30 1.19 1.25 1.18 0.30 3.5. Segmentation Results on Phantom and Patient Volume Data Location from 59 3.5. Segmentation Results on Phantom and Patient Volume Data Figure 3.12: Comparison of interpolated sagittal segmentation results (red) and transversal segmentation results (green) on patient data. 60 Chapter 4 Surface-based Pre-operative to Intra-operative Registration of the Prostate 4.1 Registration Using Ellipsoid Parameters Our approach is to fit the pre-operative and intra-operative prostate boundary points to two different ellipsoids, using the direct least square ellipsoid fitting [50], and then use the ellipsoids’ parameters to calculate the registration matrix. The long-term goal is to use the contours from only a few pre-operative and intra-operative sagittal images in order to perform the registration. However, as an initial step we will use finely sampled preoperative and intra-operative US volume data. An US volume is a set of 2D sagittal B-mode US images which are acquired while rotating the TRUS probe around its center axis. The images are captured at 0.5 degree intervals and each of them are tagged with the probe’s angular position. Our group has developed an automated US volume data acquisition system. However, if such system is not available, the conventional EXII stepper from the commercial brachytherapy apparatus is equipped with an encoder that provides probe rotation angle readings. While rotating the TRUS probe manually, sagittal images can be acquired and tagged with the angle readings from the encoder. Hence, our registration algorithm can also be used with the conventional brachytherapy system. From the volume data, one selects multiple pre-operative and intraoperative sagittal images, and contours the prostate using the sagittal segmentation algorithm from Chapter 3. Using the known probe rotation angle and US machine’s pixel to mm conversion factor, the resulting contours can be reconstructed into pre-operative and intra-operative point clouds in the 3D US coordinate frame {US}. The Z-axis of this coordinate frame is the probe’s rotation axis and is pointing away from the TRUS distal end towards the TRUS proximal end. The Y-axis of the frame is perpendicular to 61 4.1. Registration Using Ellipsoid Parameters the sagittal US crystal array and is pointing away from it in the axial imaging direction. The two axes subtend the sagittal image plane. The frame’s origin is located along the probe’s rotation axis, one probe radius distance below the first element of the sagittal US crystal array. See Figure 4.1 for illustration. Figure 4.1: Coordinate system convention for the US frame. As illustrated in Figure 4.2, the goal is to find the rigid body transformation matrix, TReg , which transforms the pre-operative point cloud to the intra-operative one. By fitting ellipsoids to the two point clouds [50], we obtain the pre-operative ellipsoid “PreOpEll” and intra-operative ellipsoid “IntraOpEll” and use them as the equivalents of the point clouds. We express the registration as follows: US XIntraOpEll = TReg US XP reOpEll (4.1) The superscript U S indicates that the points in the point clouds are expressed in the US coordinate frame. Homogeneous coordinate convention is used to express the coordinate points (U S XP reOpEll and U S XIntraOpEll ) and the transformation matrix (TReg ). From the pre-operative and intra-operative ellipsoids we obtain two sets of parameters: (xcp , ycp , zcp , Axp , Ayp , Azp , rxp , ryp , rzp ) and (xci , yci , 62 4.1. Registration Using Ellipsoid Parameters Figure 4.2: Pre-operative to intra-operative prostate registration based on ellipsoid fitting. The two ellipsoids represent the pre-operative and intraoperative prostates. 63 4.1. Registration Using Ellipsoid Parameters zci , Axi , Ayi , Azi , rxi , ryi , rzi ). These parameters are derived by defining the pre-operative ellipsoid coordinate frame {PreOpEll} and intra-operative ellipsoid coordinate frame {IntraOpEll}, as illustrated in Figure 4.3. The origins of the coordinate frames are the ellipsoids’ centers and the coordinate axes are the ellipsoids’ principal axes. The parameters (xcp , ycp , zcp ) and (xci , yci , zci ) describe the location of the ellipsoids’ centers, in mm, with respect to the origin of {US}. The parameters (Axp , Ayp , Azp ) and (Axi , Ayi , Azi ) describe the orientations, in degrees, of the principal axes of the ellipsoids’ with respect to the axes of {US}. The orientations are expressed in Euler roll-pitch-yaw convention with Axp and Axi as the yaw angles, Ayp and Ayi as the pitch angles, and Azp and Azi as the roll angles [63]. The parameters (rxp , ryp , rzp ) and (rxi , ryi , rzi ) are the lengths of the preoperative and intra-operative ellipsoids principal radii in mm. Then, we define two homogeneous transformation matrices, U S TP reOpEll and U S TIntraOpEll . These matrices describe {PreOpEll} and {IntraOpEll} relative to {US} following equations 4.2 and 4.3. US US TP reOpEll TIntraOpEll cyp czp sxp syp czp − cxp szp cxp syp czp + sxp szp xcp cyp szp sxp syp szp + cxp czp cxp syp szp − sxp czp ycp = −syp sxp cyp cxp cyp zcp 0 0 0 1 (4.2) cyi czi sxi syi czi − cxi szi cxi syi czi + sxi szi xci cyi szi sxi syi szi + cxi czi cxi syi szi − sxi czi yci = −syi sxi cyi cxi cyi zci 0 0 0 1 (4.3) Where sxp is sin(Axp ), cxp is cos(Axp ), syp is sin(Ayp ), and so on. The first three columns of U S TP reOpEll and U S TIntraOpEll express the direction of the axes of the pre-operative and intra-operative ellipsoid coordinate frames relative to the US coordinate frame’s. Their values are calculated using the ellipsoids’ orientations - (Axp , Ayp , Azp ) and (Axi , Ayi , Azi ) using the Euler roll-pitch-yaw rotation matrix [63]. The last column of these matrices express the location of the ellipsoids origins’ relative to the origin of the US coordinate frame. Using U S TP reOpEll and U S TIntraOpEll , equation 4.1 can be rewritten as 64 4.1. Registration Using Ellipsoid Parameters Figure 4.3: Ellipsoid parameters definition. 65 4.1. Registration Using Ellipsoid Parameters follows: US TIntraOpEll IntraOpEll XIntraOpEll = TReg US TP reOpEll TintraOpEll = TReg US TP reOpEll P reOpEll XP reOpEll (4.4) We model the prostate motion as a rigid motion: the registration matrix is a rigid body transformation matrix. Ideally, the pre-operative and intraoperative ellipsoids represent the same prostate (only rotated and translated) and the point cloud distribution with respect to its own ellipsoid coordinate frame stays constant: P reOpEll XP reOpEll = IntraOpEll XIntraOpEll . Hence, equation 4.4 can be simplified into: US (4.5) And the registration matrix TReg can be calculated as follows: TReg = US TIntraOpEll P reOpEll TU S (4.6) The prostate, however, suffers from shape deformation and the segmentation algorithm will not always produce the same contour points for a given prostate. The two conditions affect our assumption and translate to errors in ellipsoid fitting. Consequently, they affect the accuracy of TReg . Note, we can also obtain TReg through the typical composition of homogeneous coordinate frame transformations. Equation 4.7 shows the composition of homogeneous transformations between the US coordinate frame, the pre-operative ellipsoid coordinate frame, and the intra-operative ellipsoid coordinate frame. US P reOpEll TIntraOpEll = US TIntraOpEll = P reOpEll TP reOpEll TU S P reOpEll US TIntraOpEll TIntraOpEll (4.7) (4.8) The transformation by P reOpEll TIntraOpEll in equation 4.7 rotates and translates the pre-operative ellipsoid frame to the intra-operative ellipsoid frame. As well, P reOpEll TIntraOpEll rotates and translates the pre-operative point cloud to the intra-operative one. However, P reOpEll TIntraOpEll is expressed relative to the pre-operative ellipsoid coordinate frame while the pre-operative and intra-operative contour points are expressed in the US coordinate frame. To find the registration matrix, TReg , we need to express P reOpEll TIntraOpEll in the US frame, i.e. we need to find the similar matrix of P reOpEll TIntraOpEll in the US frame. TReg can be calculated by applying a similarity transfor- 66 4.1. Registration Using Ellipsoid Parameters mation as follows: TReg = US TReg = US TReg = US TP reOpEll P reOpEll TP reOpEll P reOpEll TIntraOpEll TU S US P reOpEll TIntraOpEll TU S P reOpEll (4.9) TU S (4.10) TIntraOpEll P reOpEll TU S (4.11) Equation 4.11 gives the same final result as equation 4.6. As the ellipsoid represent the patient’s prostate, the transformation matrix TReg registers the pre-operative prostate to the intra-operative one. The oncologist can directly use the registration parameters to manipulate the US frame by actuating the TRUS probe through the EXII Stepper and R MicroTouch° LP mount since the registration is expressed in the US frame. To give a physical meaning to the transformation, we can express the rigid motion in terms of translation (tx , ty , tz ) and Euler roll-pitch-yaw angles (Ax , Ay , Az ), using the following formulae: tx TReg 14 ty = TReg 24 tz TReg 34 ´ ³ TReg 32 arctan2 Rx TReg 33 Ry = arcsin (− TReg 31 ) ´ ³ T Rz arctan2 Reg 21 ) TReg (4.12) (4.13) 11 Where arctan2 is the four quadrant inverse tangent. In summary, our registration approach is straight forward. Ellipsoids are fit to pre-operative and intra-operative prostate contours, obtained using our sagittal segmentation algorithm. Using their parameters, the registration matrix can be calculated through matrix inverse and multiplication. For homogeneous matrices, inverse operator just involves matrix transpose and multiplication operators. 67 4.2. Registration Experiment: Simulation 4.2 Registration Experiment: Simulation Before we test the registration algorithm on the CIRS phantom, we wish to know the feasibility and accuracy of the ellipsoid fitting algorithm. If we are given a set of ellipse sample points (“ellipse slices”) from an ellipsoid with known parameters, can we accurately retrieve those parameters and with a relatively small number of ellipse slices compared to the volume size. Then, if transformations (rotations and translations) are applied to these ellipse slices, will we be able to obtain the transformation values accurately. Note that the “ellipsoid” represents the prostate volume while the “ellipse” slices represent the prostate contours on each sagittal images which are obtained at different probe rotation angles. We simulated the ellipsoid fitting based registration algorithm with the following steps: 1. We specified the pre-operative ellipsoid’s parameters (xcp , ycp , zcp , Axp , Ayp , Azp , rxp , ryp , rzp ). 2. We calculated UST P reOpEll following equation 4.2. 3. We generated the desired rigid body translation (tx , ty , tz ) and rotation (Rx , Ry , Rz ) in mms and degrees respectively. 4. We obtained the transformation operator TReg by substituting (xcp , ycp , zcp , Axp , Ayp , Azp ) in equation 4.2 with (tx , ty , tz , Rx , Ry , Rz ). 5. We calculated UST IntraOpEll following equation 4.5. 6. We calculated the intra-operative ellipsoid’s center location (xci , yci , zci ) and principal axes orientation (Axi ,Ayi , Azi ) by substituting TReg , in equation 4.12 and 4.13, with U S TIntraOpEll . 7. We generated the pre-operative ellipsoid sample points U S XP reOpEll at specific rotation angles using the method described in Appendix B.3. 8. We applied the transformation operator TReg to each of the pre-operative ellipsoid sample points to obtain the intra-operative ellipsoid sample points U S XIntraOpEll , following equation 4.1 . 9. We fitted U S XIntraOpEll into an ellipsoid and labeled its parameters as the “fitted intra-operative ellipsoid” parameters (xif it , yif it , zif it , Axif it , Ayif it , Azif it , rxif it , ryif it , rzif it ). 68 4.2. Registration Experiment: Simulation 10. We calculated U S TIntraOpEllF it using equation 4.3 and the fitted intraoperative ellipsoid parameters. 11. We calculated the “fitted registration matrix”, TRegF it by substituting UST UST IntraOpEll with IntraOpEllF it in equation 4.6. 12. We extracted the fitted rigid body translation (txf it , tyf it , tzf it ) and rotation (Rxf it , Ryf it , Rzf it ) by substituting TReg , in equation 4.12 and 4.13, with TRegF it . 13. We calculated the ellipsoid fitting errors and registration errors by comparing the intra-operative ellipsoid parameters with the fitted ones and the desired rigid motion (translation and rotation) with the fitted one. As shown in Figure 4.4, the simulation was performed with 4 different cases: 1. Case 1: using 2 ellipse slices (at probe rotation angles 95o and 85o ), 2. Case 2: using 3 ellipse slices (at probe rotation angles 95o ,90o and 850 ), 3. Case 3: using 5 ellipse slices (probe rotation 100o , 95o , 90o , 85o , and 80o ), 4. Case 4: using all ellipse slices, typically 134 slices, with 0.50 rotation angle interval between slices. At each cases, we simulated 10 random translations on each axes (X,Y , and Z), 10 random rotations on each axes (yaw, pitch, and roll), and 20 instances of random transformations (combinations of translations and rotations), called “combination”. The translations and rotation numbers were randomly selected between [-10, 10] mm and [-10, 10] degrees, respectively. We assumed these values represent realistic prostate motion values [43, 45]. Table 4.1 shows the ellipsoid fitting errors in our simulation. We report the mean squared errors of the ellipsoid’s center (in mm), orientation (in degrees), and radii length (in mm). The errors are calculated as follows: 69 4.2. Registration Experiment: Simulation Ecenter = 1 N Eorientation = 1 N Eradii = 1 N ¯¯ ¯¯ N ¯¯xi − xif it ¯¯ X ¯¯ ¯¯ ¯¯ yi − yif it ¯¯ ¯¯ ¯¯ i=1 ¯¯ zi − zif it ¯¯ ¯¯ ¯¯ N ¯¯Axi − Axif it ¯¯ X ¯¯ ¯¯ ¯¯Ayi − Ayif it ¯¯ ¯¯ ¯¯ i=1 ¯¯ Azi − Azif it ¯¯ ¯¯ ¯¯ N ¯¯rxi − rxif it ¯¯ X ¯¯ ¯¯ ¯¯ryi − ryif it ¯¯ ¯¯ ¯¯ i=1 ¯¯ rzi − rzif it ¯¯ (4.14) (4.15) (4.16) Table 4.2 shows the registration errors in our simulation. We report the mean squared errors of the translation (in mm), orientation (in degrees). The errors are calculated as follows: Etrans = 1 N Erot = 1 N ¯¯ ¯¯ N ¯¯tx − txf it ¯¯ X ¯¯ ¯¯ ¯¯ty − tyf it ¯¯ ¯¯ ¯¯ i=1 ¯¯ tz − tzf it ¯¯ ¯¯ ¯¯ N ¯¯Rx − Rxf it ¯¯ X ¯¯ ¯¯ ¯¯Ry − Ryf it ¯¯ ¯¯ ¯¯ i=1 ¯¯ Rz − Rzf it ¯¯ (4.17) (4.18) 70 4.2. Registration Experiment: Simulation Figure 4.4: Four different cases of simulated pre-operative (blue) and intraoperative (red) ellipsoids: 2 ellipse slices (top left), 3 ellipse slices (top right), 5 ellipse slices (bottom left), and all possible ellipse slices (bottom right) before registration. 71 Table 4.1: Ellipsoid fitting error under various simulation cases. X translation 10 runs 8.22E-01 Ellipsoid Fitting Error Y translation Z translation Ax (yaw) Ay (pitch) Az (roll) Combination 10 runs 10 runs 10 runs 10 runs 10 runs 20 runs 6.92E-01 7.87E-01 7.88E-01 1.21E+00 1.02E+00 2.06E+00 Eorientation (o) 6.67E+00 5.19E+00 6.19E+00 6.25E+00 8.42E+00 4.99E+00 8.67E+00 Eradii (mm) 8.19E+00 6.82E+00 7.81E+00 Ecenter (mm) 2.42E-09 6.08E-10 3.38E-10 8.06E+00 9.90E+00 1.03E+01 2 non 3 non 9 non ellipsoid ellipsoid ellipsoid 4.73E-10 1.39E-10 1.75E-10 1.30E+01 9 non ellipsoid 2.05E-09 Eorientation (o) 3.39E-08 1.26E-08 1.55E-08 2.83E-08 1.80E-09 4.94E-09 5.05E-08 Eradii (mm) 2.03E-08 1.84E-09 9.52E-10 1.11E-09 1.26E-09 1.19E-09 1.46E-08 Ecenter (mm) 1.01E-11 3.87E-12 1.35E-11 6.78E-12 3.54E-11 1.96E-10 2.71E-11 Eorientation (o) 1.22E-10 3.91E-11 1.44E-10 8.05E-11 1.21E-09 3.92E-10 3.01E-10 Eradii (mm) 8.36E-11 3.22E-11 1.11E-10 5.23E-11 1.17E-10 5.01E-10 2.08E-10 All Slices Ecenter (mm) 2.99E-11 7.37E-12 6.59E-12 1.21E-11 8.26E-13 2.15E-12 1.94E-12 Eorientation (o) 8.08E-10 1.94E-11 5.57E-11 4.18E-11 3.47E-11 1.93E-11 5.39E-11 Eradii (mm) 2.66E-10 3.08E-11 2.87E-11 4.97E-11 5.96E-12 9.96E-12 1.45E-11 3 Slices 5 Slices 4.2. Registration Experiment: Simulation Motion Type Cases MSE Fitting 2 Slices Ecenter (mm) 72 Table 4.2: Registration errors under various simulation cases. Motion Type MSE Cases Registration ERot (o) 3 Slices 5 Slices All Slices Y Z Ax (yaw) Ay (pitch) Az (roll) Combination 10 runs 10 runs 10 runs 10 runs 10 runs 10 runs 20 runs 1.32E+00 9.41E-01 1.01E+00 2.01E+00 2.22E+00 1.36E-01 5.32E+00 1.89E+00 1.31E+00 1.66E+00 6.27E+00 9 non ellipsoid 2.06E-08 Etrans (mm) 1.77E-08 6.39E-09 3.93E-09 2.48E+00 4.15E+00 3.92E-01 2 non 3 non 9 non ellipsoid ellipsoid ellipsoid 8.47E-10 1.00E-08 9.12E-09 ERot (o) 4.88E-08 1.70E-08 1.13E-08 2.33E-09 2.78E-08 2.54E-08 5.71E-08 Etrans (mm) 3.92E-11 6.80E-11 2.75E-11 8.17E-11 5.14E-10 3.96E-10 1.01E-10 o ERot ( ) 1.06E-10 2.05E-10 7.57E-11 2.67E-10 1.29E-09 4.90E-10 3.39E-10 Etrans (mm) 6.67E-10 1.88E-11 4.14E-11 2.82E-11 4.86E-11 3.50E-11 5.76E-11 o 8.29E-10 1.97E-11 6.03E-11 3.13E-11 5.14E-11 4.06E-11 6.61E-11 ERot ( ) 4.2. Registration Experiment: Simulation 2 Ellipse Slices Etrans (mm) X 73 4.3. Registration Experiment: CIRS Phantom The simulation results show that a minimum of three ellipse slices are needed to recover the ellipsoid parameters reliably. As the number of ellipse slices increase, the ellipsoid fitting error decreases. When only two ellipse slices were used, the ellipsoid fitting algorithm may fail to produce an ellipsoid. In such occurrences, the ellipse slices make the ellipsoid look thin-long or compressed. The ellipsoid fitting algorithm can’t fit these types of ellipsoids. In the simulation, such problems did not occur when three or more ellipse slices were used. Furthermore, we assume that the prostate will not assume thin-long or compressed shapes. As illustrated in Figure 4.4, our ellipse slices are partial contours that are not closed. We purposely simulated this condition to mimic our CIRS phantom images, where the prostate phantom is larger than the sagittal view of the TRUS probe. Thus, in the images, the prostate seem to be “cut off”. However, our simulation shows that we still can accurately recover the ellipsoid parameters from partial data. The small errors in Table 4.2 show that if the ellipsoid parameters are accurately retrieved, the registration algorithm gives good results. 4.3 4.3.1 Registration Experiment: CIRS Phantom Experiment Setup We next tested the registration algorithm on the CIRS phantom. The phantom was mounted on a platform that provides three axes of translation, plus yaw rotation and pitch rotation. We could change the position and pose of the phantom by manipulating the platform. This platform was constructed using three precision translation stages (Newport Corp, CA, USA) and two precision rotation stages (Parker Positioning System, CA, US). Each of the translation stage has 2 µm resolution and each of the rotation stage has 1o resolution. Figure 4.5 shows the experiment setup. Measurements from the 3020 Optotrak (NDI, ON, Canada) system was used as the gold standard for registration because the platform was not calibrated to the US coordinate frame. Four Optotrak markers were attached to the phantom and six markers were attached to the TRUS probe. These markers established coordinate frames for the phantom and the probe relative to the Optotrak’s camera coordinate frame. From these two, the Optotrak can directly report the position and orientation of the phantom in the probe coordinate frame. A minimum of three markers must be visible for the Optotrak to establish the position and orientation of each of these coordinate frames. Since the TRUS probe’s rotation will occlude some of 74 4.3. Registration Experiment: CIRS Phantom Figure 4.5: Top: the CIRS phantom mounted on our custom made platform. Optotrak markers were attached to the phantom and probe. Bottom: Optotrak camera with its data acquisition unit. 75 4.3. Registration Experiment: CIRS Phantom its Optotrak markers from the camera, a larger number of markers were used for the probe. The Optotrak reported the probe and phantom frames in terms of three axes of translation and yaw-pitch-roll angles - similar to equation 4.12 and 4.13. Figure 4.6 illustrates the coordinate frames involved in the experimental setup. The phantom’s position and pose are changed using the 5 DOF platform. The labels “pre-operative” and “intra-operative” are used to differentiate the original pose and position of the phantom with the newer ones. Equations 4.19 - 4.21 express composition of homogeneous coordinate transformations in the experimental setup. US TP reOpEll = US TP robe P robe TP reOpP tm P reOpP tm TP reOpEll (4.19) US TIntraOpEll = US TP robe P robe TIntraOpP tm IntraOpP tm TIntraOpEll (4.20) IntraOpP tm TIntraOpEll = P reOpP tm TP reOpEll = P tm TEll (4.21) P tm T Ell expresses the pose and position of the of the synthetic prostate (inside the phantom) with respect to the phantom’s coordinate frame which is established by the four Optotrak markers. The value of this matrix is unknown but is assumed to be constant as the synthetic prostate is rigidly positioned within the phantom. Both P robe TP reOpP tm and P robe TIntraOpP tm were directly measured using the Optotrak. The transformation mapping the US coordinate frame to the probe coordinate frame, P robe TU S , was found by calibrating the probe with the Optotrak (briefly explained in the next sub section and in Appendix D). Similar to equation 4.1, we calculate the registration matrix TReg : TReg = US TReg = US TIntraOpEll TP robe P robe P reOpP tm TReg = US P reOpEll TU S TIntraOpP tm P tm TEll Ell TP tm P robe TP robe TU S P robe TP robe TIntraOpP tm P reOpP tm TP robe P robe TU S (4.22) TReg in equation 4.22 can be directly compared with the one obtained using the ellipsoid fitting based registration algorithm (equation 4.6) because both of them are expressed in the same reference frame (the US coordinate frame). The experimental data collection was carried out as follows: 76 4.3. Registration Experiment: CIRS Phantom Figure 4.6: The coordinate frames used in the experiment and their respective transformations. 77 4.3. Registration Experiment: CIRS Phantom 1. We started with a particular pose and position for the phantom and recorded the phantom’s and probe’s Optotrak readings. 2. We captured the US volume datum at its starting position and pose. 3. We adjusted the stage dials on the platform to change the phantom’s position and pose (while still being able to acquire US images of the prostate). 4. We recorded the new Optotrak readings of the phantom and the probe. 5. We captured a new US volume datum at the new position and pose. 6. We repeated step 3 to 5 for the next volume datum until 25 volume data for 25 phantom position and poses are collected. The 25 poses were divided into seven categories. From pose 1, the starting pose, to pose 3, no motion was applied to the phantom. In poses 4 - 7, the phantom was translated 5 mm back and forth in the X-axis direction twice. In poses 8 - 11 and poses 12 - 15 the phantom was translated in similar manner but in Y-axis and Z-axis direction. In poses 16 - 19 back and forth 5o yaw angulation was applied to the phantom. Similarly, in poses 20 - 23 pitch angulation was applied to the phantom. In the remaining two poses, rotation and translation combinations were applied to the phantom. From each volume, we contoured the prostate phantom with two approaches. The first approach was to individually segment ten to eleven sagittal images using the sagittal segmentation algorithm (explained in Section 3.3). The images were selected to be 5o apart and were symmetric around the 90o probe rotation angle, e.g. at rotation angles 90o , 85o , 95o , 100o , 80o and so on. The second approach was to segment the whole volume using our extended sagittal segmentation algorithm (explained in Section 3.4). As mentioned before, our ultimate long-term goal is to be able to do registration with ellipsoid fitting by just using several pre-operative and intra-operative images (similar to the first approach). However, if the goal is not feasible, we can use the whole volume to generate a complete set of contours for registration. Thus, in this experiment we compared the two approaches. Hereafter the first approach is called the sparse volume registration while the second one is called the full volume registration. Using the ellipsoid fitting-based registration, the contours were used to compute the transformation between adjacent poses (pose 1 and 2, pose 2 and 3, and so on). The results were compared with transformations calculated using Optotrak readings of the phantom. The transformation param78 4.3. Registration Experiment: CIRS Phantom eters were expressed in terms of three axes translations (tx ,ty ,tz ) and three axes rotations (Rx ,Ry ,Rz ). 4.3.2 TRUS Probe Calibration We opted to use the “single wall ultrasound calibration technique” [64] due to its convenience. The technique only required the use of a water bath with a flat floor and it is already incorporated within Stradwin. We performed the calibration off-line as Stradwin does not support real time probe calibration using the 3020 Optotrak. Optotrak data and US images had to be collected independently and then combined into the Stradwin’s standard file format. The file was then uploaded into Stradwin and the calibration steps followed its built in functions. The step by step calibration process is described in Appendix D. 4.3.3 Registration Experiment Results Tables 4.3 and Table 4.4 summarize the errors of the ellipsoid fitting registration based on the two approaches: the sparse volume registration and the full volume registration. In the tables, we report the average, standard deviation, and maximum of the registration errors for the five different categories of phantom motion (no motion, X-axis translation, Y-axis translation, and so on). We also report the bounding box registration error “bounding box error”) and the Hausdorff distance between the registered surfaces. The bounding box errors were calculated by defining two smallest rectangular boxes that cover the intra-operative and pre-operative contour points. The boxes’ vertices are the permutations of the X-axis, Y-axis, and Z-axis extrema of the intra-operative and pre-operative contour points. The later is registered using the registration matrix obtained from the ellipsoid fitting based registration algorithm. Then, the bounding box error is calculated by taking the root mean square distance between corresponding vertices. The complete registration experiment results can be found in Appendix E. 79 Table 4.3: Summary of registration errors: the sparse volume registration. Motion Type Transformation Title Absolute Tx (mm) value from optotrak Error |Ellipsoid Fit Reg - No Motion Std. Dev. Max X axis translation Average Std. Dev. Max Y axis translation Average Std. Dev. Max Z axis translation Average Std. Dev. Max 0.00 0.00 0.01 3.02 0.05 3.06 0.14 0.00 0.14 0.51 0.00 0.51 Ty (mm) 0.01 0.01 0.01 0.37 0.03 0.41 2.92 0.04 2.96 0.12 0.02 0.15 Tz (mm) 0.00 0.00 0.00 0.26 0.01 0.26 0.03 0.01 0.04 5.01 0.00 5.01 Rx ( o ) 0.00 0.00 0.00 0.01 0.00 0.02 0.36 0.00 0.36 0.02 0.01 0.03 Ry ( o ) 0.00 0.00 0.00 0.25 0.00 0.25 0.06 0.01 0.07 0.03 0.00 0.04 Rz ( o ) Tx (mm) 0.01 0.01 0.02 0.20 0.01 0.20 0.18 0.01 0.19 0.02 0.01 0.03 1.38 0.72 2.08 8.56 2.41 11.74 4.60 3.22 7.36 11.67 6.85 17.43 1.64 Ty (mm) 0.36 0.27 0.55 0.14 0.16 0.35 0.68 0.48 1.25 1.10 0.52 Tz (mm) 0.65 0.44 1.07 1.48 0.83 2.41 1.20 0.84 2.04 4.99 2.36 7.46 Rx ( o ) 0.68 0.48 1.03 0.49 0.30 0.73 0.99 0.96 2.24 1.59 1.34 3.31 Ry ( o ) 3.37 35.44 2.04 5.07 16.63 3.16 20.91 7.83 4.38 11.83 25.23 11.77 Rz ( o ) Optotrak| Bounding Box Error (mm) 0.81 0.67 1.24 0.82 0.80 1.90 1.07 0.85 2.12 2.09 1.70 4.13 1.00 0.36 1.32 3.16 0.85 4.41 2.30 0.11 2.37 4.65 2.25 7.52 Hausdorff Distance (mm) 4.32 0.28 4.57 12.49 1.11 14.01 9.71 1.61 11.06 15.38 3.78 18.16 4.3. Registration Experiment: CIRS Phantom motion Average 80 Table 4.3 continued from previous page Motion Type Title Transformation Absolute Tx (mm) motion optotrak Error |Ellipsoid Fit Reg - Y axis rotation (pitch) Average Std. Dev. Max Combination Average Std. Dev. Max 0.19 0.01 0.22 0.20 0.11 0.31 0.70 0.21 0.85 Ty (mm) 0.33 0.10 0.43 0.23 0.03 0.25 0.46 0.26 0.65 Tz (mm) 1.64 0.01 1.65 4.75 0.05 4.79 6.72 0.05 6.76 Rx ( o ) 4.95 0.03 4.99 0.35 0.01 0.36 3.38 0.01 3.39 Ry ( o ) 0.30 0.02 0.32 4.72 0.09 4.79 2.85 0.03 2.87 Rz ( o ) Tx (mm) 0.38 0.01 0.39 0.21 0.02 0.22 0.87 0.09 0.93 3.76 2.53 6.92 9.28 5.70 13.12 6.51 7.33 11.69 Ty (mm) 0.50 0.23 0.81 1.64 0.63 2.19 1.08 0.79 1.64 Tz (mm) 3.20 0.89 4.00 2.46 3.24 7.19 1.52 1.04 2.26 Rx ( o ) 1.40 0.58 1.76 2.42 2.22 5.66 1.93 0.46 2.25 Ry ( o ) 10.90 8.74 4.99 12.30 11.92 7.43 19.50 7.25 5.16 Rz ( o ) Optotrak| Bounding Box Error (mm) 0.99 0.61 1.82 3.09 2.11 5.26 6.03 2.28 7.64 4.22 0.12 4.31 3.69 2.20 6.42 4.26 2.24 5.84 Hausdorff Distance (mm) 11.55 0.49 12.19 14.69 2.19 16.00 8.51 2.07 9.97 4.3. Registration Experiment: CIRS Phantom value from X axis rotation (yaw) Average Std. Dev. Max 81 Table 4.4: Summary of registration errors: the full volume registration. Motion Type Transformation Title Absolute Tx (mm) value from optotrak Error |Ellipsoid Fit Reg - No Motion Std. Dev. Max X axis translation Average Std. Dev. Max Y axis translation Average Std. Dev. Max Z axis translation Average Std. Dev. Max 0.00 0.00 0.01 3.02 0.05 3.06 0.14 0.00 0.14 0.51 0.00 0.51 Ty (mm) 0.01 0.01 0.01 0.37 0.03 0.41 2.92 0.04 2.96 0.12 0.02 0.15 Tz (mm) 0.00 0.00 0.00 0.26 0.01 0.26 0.03 0.01 0.04 5.01 0.00 5.01 Rx ( o ) 0.00 0.00 0.00 0.01 0.00 0.02 0.36 0.00 0.36 0.02 0.01 0.03 Ry ( o ) 0.00 0.00 0.00 0.25 0.00 0.25 0.06 0.01 0.07 0.03 0.00 0.04 Rz ( o ) Tx (mm) 0.01 0.01 0.02 0.20 0.01 0.20 0.18 0.01 0.19 0.02 0.01 0.03 1.72 1.29 2.61 2.33 0.44 2.95 0.94 0.44 1.46 0.64 0.93 2.03 Ty (mm) 0.31 0.15 0.48 1.05 0.52 1.61 0.68 0.48 1.11 0.92 0.22 1.25 Tz (mm) 0.38 0.25 0.55 2.08 0.40 2.65 1.48 1.27 2.86 2.00 1.14 3.46 Rx ( o ) 0.41 0.21 0.64 2.02 0.65 2.82 1.18 1.42 2.91 1.88 0.67 2.78 Ry ( o ) 5.45 3.19 1.79 4.79 4.33 0.96 5.51 2.08 2.15 4.95 2.87 2.04 Rz ( o ) Optotrak| Bounding Box Error (mm) 0.84 0.36 1.26 0.95 0.25 1.30 0.72 0.38 1.23 1.35 0.71 2.33 0.76 0.23 1.01 1.03 0.08 1.12 1.05 0.15 1.25 3.23 0.92 3.87 Hausdorff Distance (mm) 2.76 0.21 2.97 3.97 0.25 4.22 3.41 0.35 3.70 7.34 0.45 7.93 4.3. Registration Experiment: CIRS Phantom motion Average 82 Table 4.4 continued from previous page Motion Type Title Transformation Absolute Tx (mm) motion optotrak Error |Ellipsoid Fit Reg - Y axis rotation (pitch) Average Std. Dev. Max Combination Average Std. Dev. Max 0.19 0.01 0.22 0.20 0.11 0.31 0.70 0.21 0.85 Ty (mm) 0.33 0.10 0.43 0.23 0.03 0.25 0.46 0.26 0.65 Tz (mm) 1.64 0.01 1.65 4.75 0.05 4.79 6.72 0.05 6.76 Rx ( o ) 4.95 0.03 4.99 0.35 0.01 0.36 3.38 0.01 3.39 Ry ( o ) 0.30 0.02 0.32 4.72 0.09 4.79 2.85 0.03 2.87 Rz ( o ) Tx (mm) 0.38 0.01 0.39 0.21 0.02 0.22 0.87 0.09 0.93 0.62 0.44 0.99 1.70 0.46 2.21 3.73 0.25 3.91 Ty (mm) 0.71 0.44 1.17 0.60 0.26 0.93 0.74 0.70 1.24 Tz (mm) 4.34 0.79 5.02 1.37 0.58 2.08 0.60 0.62 1.04 Rx ( o ) 2.11 0.69 2.72 1.51 0.76 2.39 0.71 0.32 0.93 Ry ( o ) 1.88 2.06 1.99 4.58 1.16 0.26 1.44 1.79 0.39 Rz ( o ) Optotrak| Bounding Box Error (mm) 1.01 0.71 1.67 1.54 0.59 2.34 3.05 0.25 3.23 4.73 0.34 5.01 3.09 0.95 4.15 2.61 1.20 3.46 Hausdorff Distance (mm) 12.93 0.18 13.12 8.34 0.94 9.37 7.36 0.06 7.40 4.3. Registration Experiment: CIRS Phantom value from X axis rotation (yaw) Average Std. Dev. Max 83 4.4. Discussion The results in Table 4.3 and 4.4 show that the sparse volume registration results in significant translational, angle, and bounding box errors. Noticeably, there is a 17.43 mm error in X-axis translation, 35.44o pitch (Y-axis) angle error, and 7.52 mm bounding box error when the phantom was only moved 5 mm in Z-axis direction. These error values are the largest amongst the other motion categories. The full volume registration produces smaller maximum error values: 5.02 mm Z-axis translation error, 5.51o pitch (Y-axis) angle error, and 5.01 mm bounding box error. We compared the bounding box error and Hausdorff distance between the sparse volume registration and the full volume registration using the paired t-test. The hypothesis is that the sparse volume registration has significantly larger bounding box error and significantly larger Hausdorff distance. The paired t-test on the bounding box error concludes that the sparse volume registration has a larger bounding box error with P-Value of 0.0015 (less than 0.05) and mean difference range of [0.3995 1.4838] within the 95% confidence interval. The paired t-test on the Hausdorff distance concludes that the sparse volume registration has a larger Hausdorff distance with P-Value of 5.2586 .10−6 (less than 0.05) and mean difference range of [3.0553 6.4093] within the 95% confidence interval. Hence, the sparse volume registration does not seem to be a feasible approach for ellipsoid fitting based registration. For registration purposes, we need to use the whole US volume. It is possible that the large difference in the number of sample points is the reason behind the difference in errors. On average, the sparse volume registration fits about 1100 contour sample points into an ellipsoid; compared to 13800 points fit by the full volume registration (see Appendix F). Using more sample points seems to produce a better ellipsoid fit on real US contour data as shown by the root mean square algebraic distance of the ellipsoid fit of the two approaches. The ellipsoid fit for the sparse volume registration has a root mean square algebraic distance of 0.011 ± 0.001 mm. While, the fit for the full volume registration has a root mean square algebraic distance of 0.008 ± 0.0005 mm. 4.4 4.4.1 Discussion Ellipsoid Orientation Fit Differences Cause Translation and Distance Registration Errors As shown in Table 4.4, we have as much as (2.61, 0.48, 0.55) mm maximum X-Y-Z axes translational errors for a (0.64o , 4.79o , 1.26o ) maximum angular 84 4.4. Discussion errors even when no motion was applied to the phantom (pose 1 to pose 3). These translation errors correspond to about 2.7 mm absolute distance error. On closer inspection of the ellipsoid fit results (Table F.2), we noticed that the ellipsoid centers in pose 1 to pose 3 are almost constant. The maximum absolute difference between the three ellipsoid centers are (0.37, 0,12, 0.07) mm in the X-Y-Z axes respectively - corresponding to a 0.4 mm absolute distance difference. However, the maximum orientation angle differences between the three ellipsoids are (0.60, 4.80, 1.3) degrees. We believe that the 2.7 mm absolute distance error in pose 1 to pose 3 are due to differences in the fit of the ellipsoids’ orientations. This is more apparent when equation 4.6 is expressed as: TReg = TReg = US · TIntraOpEll RReg tReg 0 1 P reOpEll ¸ TU S (4.23) (4.24) Where RReg = US tReg = US RIntraOpEll P reOpEll US tIntraOpEll − RU S (4.25) RIntraOpEll P reOpEll RU S US tP reOpEll (4.26) And US tIntraOpEll xci = yci zci US tP reOpEll xcp = ycp zcp (4.27) The rotation matrices U S RIntraOpEll and U S RP reOpEll are calculated using the intra-operative and pre-operative ellipsoids’ orientations respectively. Equation 4.26 shows that the translation part of the registration matrix is affected by both the orientation of the ellipsoids and the location of the ellipsoids’ centers. Although the ellipsoid centers are close together, the translational errors will be large if the orientation angles differences are large. The further away the ellipsoid’s center from the coordinate center of the US frame, the larger the translational errors will be due to lever arm effects. Hence, the key to getting an accurate registration results using ellipsoid fitting is to make the ellipsoid orientation parameter as accurate as possible. This is a very difficult task since the prostate is deformable and, being a soft tissue, it does not have any distinct landmarks. Moreover, the ellipsoid fit depends on the contours from the segmentation algorithm which includes some errors. 85 4.4. Discussion 4.4.2 Sensitivity of Final Ellipsoid Parameters to Initialization Points Large differences (up to 4.8o ) in the ellipsoid orientation parameters in volume 1-3 seem to be caused by sensitivity to initialization points. Volumes 1 to 3 are collected with the phantom being still. Hence, the US images in the volumes should be similar. To confirm, we cross correlated corresponding images between the volumes. Images correspond to each other when they are taken at the same probe rotation angle. The cross correlation coefficient is 0.9967 ± 0.0027. Thus, the only difference between the three volumes’ contours are their choices for initialization points. This affects the results of the segmentation and hence the ellipsoid fit. We performed a simple sensitivity test to confirm this hypothesis. We segmented volume 1 with the extended segmentation algorithm by varying the location of the first (and only the first) initialization point. This point is the one located at the top of the prostate in its mid sagittal view. From its original location, we independently moved the point in the X,Y,and Z axes directions and observed the resulting ellipsoid parameters. At most, the point is moved 1.56 mm (10 pixels) away from its original location. Figures 4.7 to 4.9 show the observation results. 86 4.4. Discussion Figure 4.7: Effects of changing one initialization point in X-axis direction to ellipsoid fit parameters. The plots show linear fits of the observation along with its R2 (the square of the correlation coefficient) and the line equation. The middle graph at the bottom row shows that the roll angle is linearly dependent to X-axis location of the initialization point with a slope of -6.28 o /mm. 87 4.4. Discussion Figure 4.8: Effects of changing one initialization point in Y-axis direction to ellipsoid fit parameters.The plots show linear fits of the observation along with its R2 (the square of the correlation coefficient) and the line equation. 88 4.4. Discussion Figure 4.9: Effects of changing one initialization point in Z-axis direction to ellipsoid fit parameters. The plots show linear fits of the observation along with its R2 (the square of the correlation coefficient) and the line equation. 89 4.4. Discussion The graphs show several strong linear dependencies (R2 >0.8) between several ellipsoid parameters and the initialization point location. For example, the slope of the linear fit between the ellipsoid roll angle (Z-axis angle) and the initialization point’s X-axis location is -6.28 (Figure 4.7). This means that a one mm change in the initialization point results in about -6 degrees change in ellipsoid roll angle. In our extended sagittal segmentation algorithm, the effects of the 14 initialization points locations toward the ellipsoid fit parameters are combined in a complex and non linear manner. Nevertheless, our simple sensitivity test shows that the ellipsoid fit parameters are, indeed, sensitive to the locations of the initialization points. This is something that we did not expect prior to implementing the ellipsoid fitting based registration algorithm. The user determines the choice of initialization points based on how well they can discern the prostate boundaries in the US images. This judgment depends on the US image quality which is affected by speckle noise and various imaging artifacts. For example, in our phantom images the top boundary of the prostate is unclear (see Figure 3.3) due to depth artifacts and shadowing which affects the choice of the first initialization point. In real prostate data, the patient’s physical conditions also affect the image quality, for example, gas, liquid, or calcification might be present in the rectum and prostate reducing the quality of the images. Such occurrences are mostly uncontrollable. Hence, variability in user initialization points is unavoidable. We can also improve the registration algorithm using other surface based registration approach that do not depend on specific model fitting, for example iterative closest point (ICP), as explored in Section 4.5 4.4.3 Role of the Kalman Filter Segmentation To explore the role of the Kalman filter based segmentation in the final ellipsoid fit, we fit the fourteen user initialization points into ellipsoids and observed the sensitivity of their parameters. Recall, in the extended segmentation algorithm, the user selected fourteen initialization points that guided the IMMPDAF edge detector. We used the same data set from Section 4.4.2. The complete results of this observation can be found in Appendix H. In Figure 4.10, we select and plot the ellipsoid fit parameters sensitivity on several cases. We compare the sensitivity of fitting the user initialization points (magenta line) and the sensitivity of fitting the final contour points after the Kalman filter segmentation (blue line). The graphs show that the Kalman filter segmentation helps reduce the sensitivity of the ellipsoid fit 90 4.4. Discussion parameters with respect to the user initialization points locations. In the ellipsoid’s roll axis (Az ) versus X-axis location of the initialization point plot, the slope of the linear fit decreases from −8o to −6o . In the ellipsoid’s yaw axis (Ax ) versus Y-axis location of the initialization point plot, the linear fit measure (R2 ) decreases from 0.97 to 0.57; meaning that there is no longer a strong linear dependency between the two. The less sensitive the ellipsoid fit results with respect to the initialization points locations, the more accurate is the registration results. 91 4.4. Discussion Figure 4.10: Comparison between the ellipsoid parameters sensitivity. The magenta line shows the results of fitting an ellipsoid to the user initialization points under different initialization point locations. The blue one shows the result of fitting the final contour (after the Kalman filter segmentation). As shown in the plots, the Kalman filter segmentation reduces the sensitivity of the ellipsoid parameters to user initialization points locations. 92 4.5. Iterative Closest Point (ICP) 4.5 Iterative Closest Point (ICP) We re-registered our prostate contours using the Iterative Closest Point algorithm [65] - a commonly used registration algorithm. We used the contours which were obtained using the extended sagittal segmentation algorithm. In our ICP registration, the pre-operative point cloud becomes the model for the registration while the intra-operative point cloud is the data to be registered. We under-sampled the data’s sample points to make the ICP algorithm run significantly faster. Since ICP does not require specific model fitting, we expect it to be less sensitive to the location of the initialization points and to give more accurate registration results. Table 4.5 summarizes the ICP registration errors, reported in a similar manner as Tables 4.3 and 4.4. The complete ICP registration results can be found in Appendix E.3. As expected, the maximum errors for the no-motion category, 1.42 mm X-axis translational error, 1.37o yaw (X-axis) angular error, and 0.76 mm bounding box error, are less than the full volume registration’s. The overall maximum errors, across all categories, are 3.28 mm Z-axis translational error, 4.43o pitch (Y-axis) angular error, and 4.19 mm bounding box error. They are also less than the full volume registration results. We also compared the bounding box error and Hausdorff distance between the the full volume registration and ICP registration using the paired t-test. The hypothesis is that the full volume registration has larger bounding box error and larger Hausdorff distance. The paired t-test on the bounding box error concludes that the full volume registration has a larger bounding box error with P-Value of 0.0073 and mean difference range of [0.1011 0.5825] within the 95% confidence interval. The paired ttest on the Hausdorff distance concludes that the full volume registration has a larger Hausdorff distance with P-Value of 0.0032 and mean difference range of [0.1740 0.7700] within the 95% confidence interval. Compared to the full volume registration, ICP registration proves to give better results, although not significantly (based on the 95% confidence interval range). Figure 4.11 illustrates that the registration results between ICP and the full volume registration are visually similar. The smaller errors in the no motion category suggest that ICP is less sensitive to initialization points. However, there is a trade off between registration error and execution time. ICP registration produces more accurate results as measured but it takes longer to execute. Overall, it takes ICP 76.49 ± 28.84 seconds to register the prostate contours even after under-sampling the data points. It only takes 0.015 ± 0.014 seconds to fit an ellipsoid to the 93 4.5. Iterative Closest Point (ICP) contours and 0.02 ± 0.025 seconds to register the prostate using the ellipsoid fitting registration algorithm. Both algorithm were implemented in Matlab running on a 1.66 GHz computer with 1 GB RAM. The ICP algorithm is sensitive to the presence of noise and outliers in the data [66–68]. In our study, we have not measured the robustness of the ellipsoid fitting registration towards outliers nor compared it with the ICP algorithm. Figure 4.11: Results of the full volume registration (left column) and ICP registration (right column) when X-axis rotation (top row) and X-axis translation (bottom row) are applied to the phantom. The blue, red, and green contours are the pre-operative, intra-operative, and registered preoperative volumes. 94 Table 4.5: Summary of registration errors: ICP. Motion Type Title Transformation Absolute Tx (mm) No Motion Average Std. Dev. Max X axis translation Average Std. Dev. Max Y axis translation Average Std. Dev. Max Z axis translation Average Std. Dev. Max 0.00 0.01 3.02 0.05 3.06 0.14 0.00 0.14 0.51 0.00 0.51 0.01 0.01 0.01 0.37 0.03 0.41 2.92 0.04 2.96 0.12 0.02 0.15 Tz (mm) 0.00 0.00 0.00 0.26 0.01 0.26 0.03 0.01 0.04 5.01 0.00 5.01 Rx (o) 0.00 0.00 0.00 0.01 0.00 0.02 0.36 0.00 0.36 0.02 0.01 0.03 Ry (o) 0.00 0.00 0.00 0.25 0.00 0.25 0.06 0.01 0.07 0.03 0.00 0.04 Rz (o) Tx (mm) 0.01 0.01 0.02 0.20 0.01 0.20 0.18 0.01 0.19 0.02 0.01 0.03 1.18 0.22 1.42 0.66 0.47 1.21 0.50 0.44 1.15 0.45 0.48 1.16 Ty (mm) 0.25 0.11 0.33 0.73 0.42 1.15 0.46 0.30 0.77 0.74 0.34 1.10 Tz (mm) 0.36 0.20 0.58 1.18 0.64 1.71 0.86 0.96 2.18 1.59 1.00 2.88 Rx (o) 0.50 0.22 0.65 1.18 0.78 1.75 0.86 1.05 2.40 1.23 0.65 2.08 Ry (o) 1.17 0.18 1.37 1.53 0.84 2.70 0.75 0.83 1.96 0.82 0.65 1.61 Rz (o) Optotrak| Bounding Box Error (mm) 0.42 0.32 0.75 0.72 0.67 1.45 0.11 0.13 0.26 0.71 0.87 2.00 0.55 0.25 0.76 0.67 0.11 0.79 0.88 0.10 1.03 3.15 0.70 3.81 Hausdorff Distance (mm) 2.69 0.17 2.87 4.13 0.18 4.33 3.35 0.32 3.71 6.78 0.75 7.37 motion value from optotrak Error |ICP Reg - 4.5. Iterative Closest Point (ICP) 0.00 Ty (mm) 95 Motion Type Title Transformation Absolute Tx (mm) motion value from Reg - 0.19 0.01 0.22 Y axis rotation (combination) Average Std. Dev. Max 0.20 Combination Average Std. Dev. Max 0.11 0.31 0.70 0.21 0.85 0.65 Ty (mm) 0.33 0.10 0.43 0.23 0.03 0.25 0.46 0.26 Tz (mm) 1.64 0.01 1.65 4.75 0.05 4.79 6.72 0.05 6.76 Rx (o) 4.95 0.03 4.99 0.35 0.01 0.36 3.38 0.01 3.39 Ry (o) 0.30 0.02 0.32 4.72 0.09 4.79 2.85 0.03 2.87 Rz (o) Tx (mm) 0.38 0.01 0.39 0.21 0.02 0.22 0.87 0.09 0.93 0.36 0.26 0.65 0.38 0.29 0.73 1.29 0.22 1.45 Ty (mm) 0.60 0.51 1.33 0.88 0.70 1.88 0.54 0.34 0.78 Tz (mm) 1.11 0.54 1.64 1.62 1.12 3.28 2.28 0.64 2.73 Rx (o) 0.66 0.83 1.90 1.70 1.40 3.57 1.74 0.39 2.02 Ry (o) 1.03 0.93 2.39 3.08 1.06 4.43 1.11 0.10 1.18 Rz (o) Optotrak| Bounding Box Error (mm) 0.56 0.32 0.81 1.15 0.42 1.46 1.90 0.32 2.13 3.25 0.16 3.41 3.30 0.87 4.19 2.45 0.75 2.98 Hausdorff Distance (mm) 12.20 0.11 12.32 7.47 0.71 8.51 5.67 0.26 5.86 4.5. Iterative Closest Point (ICP) optotrak Error |ICP X axis rotation (yaw) Average Std. Dev. Max 96 Chapter 5 Conclusions This thesis project was carried out with two main objectives. The first objective is to design, implement, and characterize a robotic needle guide for prostate brachytherapy. The second is to develop an ultrasound Bmode images-based pre-operative to intra-operative prostate registration algorithm. The summary, contributions, and future work of this thesis are divided in the following sections based on the two objectives. 5.1 5.1.1 A Robotic Needle Guide for Prostate Brachytherapy Summary The author has participated in the design and implementation of a novel robotic needle guide for prostate brachytherapy (Brachyguide). The four degree of freedom robot is made of two linear and rotation stages and driven by commercial off-the-shelf DC motors. The robot is designed with custom features - manual cranks, digital linear calipers, and quick release mechanisms - that enable it to be operated during control system or electrical failure. The robot is controlled with a new GUI design that provides needle plan extraction from the conventional treatment planning software (VariSeed) and digital caliper readings display. Upon completion of its first prototype, the robot’s characteristics were studied. Then a seed implant study on a commercial prostate phantom involving an oncologist and a medical physicist from the BCCA was conducted to assess the usability of the robot. The study showed that integration of the robot does not add to the implant procedure time even when it was the first use of the robot by the oncologist. To be specific, the author’s participation in the design and implementation of the robot can be divided as follows: 1. The determination and specification of suitable DC Motors for the robot. 97 5.1. A Robotic Needle Guide for Prostate Brachytherapy 2. The implementation of the robot’s electrical hardware: cables, connectors, and controller box. 3. The tuning of the robot’s controller’s. 4. The design and implementation of the robot’s graphical user interface. 5. The design and implementation of the robot’s characterization test. 6. Participation as the robot’s operator in the usability test with the oncologist from the BCCA. The mechanical design and implementation of the robot was contributed by the author’s supervisor, Dr. Tim Salcudean. 5.1.2 Contributions The robotic needle guide has resulted in the following contributions: • Better needle positioning accuracy for prostate brachytherapy: The robot has a needle positioning accuracy along the X-Y axes of less than 0.12 mm and 0.1 mm respectively . These values are better than the 5 mm needle positioning accuracy of the conventional template grid. • The capability of performing angled needle insertions: The robot provides yaw and pitch angulations (-300 to +300 and +240 to -180 respectively) to the needle with up to 0.05o accuracies. These added degrees of freedom are an improvement over the parallel trajectory of the template grid based approach which help with pubic arch interference in patients. Moreover, seed position and needle trajectory plan can now be optimized by including angled needle trajectories. • A convenient method of performing needle insertion: The GUI enables the needle insertion plan to be extracted from VariSeed - the conventional and most commonly used prostate brachytherapy planning software. Once imported into the robot’s GUI, the oncologist can step through the individual needle positions with a click of a button: eliminating the need to manually locate the grid location in the template grid approach. 98 5.2. Pre-operative to Intra-operative Registration of Prostate Volumes 5.1.3 Future Work The following are the possible future work and improvements involving the robotic needle guide: • Clinical integration and patient study: As proven in our phantom study, the robot’s usage does not add to the overall seed implant procedure time. This result paves way for the robot to be used in a clinical setting for in-vivo seed accuracy placement studies. • Future needle path planning studies: With the capability of angled needle trajectories, the robot can be used to aid with needle path planning studies such as suggested in [13]. • A lighter robot design: The current robot weighs about 2.6 kg. Possible design changes might reduce the weight of the robot to make it comparable to other designs such as the design found in [30]. • Integration of real time US visualization into the robot’s GUI : Adding real time US visualization to the robot’s GUI may be beneficial. On the real time US images, the target location of the needle, the preoperative prostate contours, and dosimetry plans can be displayed to give a holistic view of the needle insertion procedure. 5.2 5.2.1 Pre-operative to Intra-operative Registration of Prostate Volumes Summary We have developed a novel US image-based algorithm to register a preoperative prostate volume to an intra-operative volume. The algorithm involves a semi automatic segmentation of sagittal US images and registration of the resulting contour points based on ellipsoid fitting. We have developed a semi automatic algorithm to segment sagittal US images of the prostate. Note, transverse images-based segmentation is the most commonly used approach for prostate contouring in brachytherapy. Initially, we only used the segmentation algorithm to segment 2D sagittal images. Then, we extended the algorithm to be able to segment sagittal view-based 3D US volume data (each volume comprising a set of 2D US sagittal images). We validated the extended sagittal segmentation algorithm by comparing it with transverse image-based segmentation using the algorithm by Mahdavi et. al. [54]. 99 5.2. Pre-operative to Intra-operative Registration of Prostate Volumes We registered the pre-operative and intra-operative prostate contours using a novel ellipsoid fitting based algorithm. The contours are fit into two different ellipsoids and their parameters are used to calculate the registration matrix. The registration accuracy has been tested in an experimental setup using the CIRS prostate phantom and the Optotrak. Unfortunately, the results show that our ellipsoid fitting based algorithm suffers from sensitivity to user initialization during segmentation. The conventional ICP registration method produced more accurate results, but at greater computational cost. This means there is a current tradeoff between accuracy and speed. 5.2.2 Contributions The pre-operative to intra-operative registration work makes the following contributions: • A novel prostate segmentation algorithm and approach: We have developed a novel semi automatic segmentation algorithm based on sagittal US B-mode images. The algorithm is capable of segmenting sagittal view-based 3D US volume data. On average, it takes the algorithm 73 seconds to segment one volume which consists of 200 images. The algorithm was compared to a published transverse view-based segmentation algorithm using patient and phantom volume data. The mean average distance of the contours is 1.17 ± 0.31 mm. The algorithm can be used to aid prostate segmentation during brachytherapy preoperative volume studies as well as in intra-operative planning. Contouring the sagittal US volume data provides a better identification of the base and apex of the prostate which has always been difficult to do with the transverse images. The sagittal contours provide a more holistic view of the prostate since the sagittal US volume data can be more finely sampled than the transverse US volume data which are sampled at 5 mm depth intervals. • A novel ellipsoid fitting-based registration algorithm: We have developed a novel registration algorithm based on ellipsoid fitting. The algorithm execution time is fast: 0.02 ± 0.025 seconds. The registration accuracy however, is, sensitive to user initialization during the image segmentation. This leads to a 5.02 mm maximum translation error, 5.51 o maximum angular error, 5.01 mm maximum bounding box error, and 6.68 mm average Hausdorff distance between the registered surfaces. In our study, the main application of the registration 100 5.2. Pre-operative to Intra-operative Registration of Prostate Volumes algorithm is in the field of prostate brachytherapy. It is possible to apply the algorithm for other prostate cancer treatment options such as EBRT. • Possible registration solution using ICP: We have explored one possible solution to improve the registration accuracy by using ICP to register the prostate contours. ICP takes longer to execute (76 ± 29 seconds) but result in lower registration errors: 3.28 mm maximum translation error, 4.43 o maximum angular error, 4.19 mm maximum bounding box error, and 6.20 mm average Hausdorff distance between the registered surfaces. These lower error values suggest that registration without specific model fitting produces lower errors but at greater computational cost. • Minor contributions in the design of registration experiment setup and calibration setup using the Optotrak : To measure the accuracy of our registration algorithm we have designed an experimental setup using the Optotrak as the gold standard. This experimental setup can be re-used for future registration studies or for real time tracking experiments. Currently, Stradwin does not support probe calibration using the 3020 Optotrak. During our probe calibration, we developed the protocols for combining the Optotrak data into Stradwin standard files in order to perform off-line calibration. The protocols can be used for future Optotrak-based calibration or to develop online calibration using the Optotrak. 5.2.3 Future Work The following are the possible directions for future work and improvements involving the pre-operative to intra-operative prostate volume registration: • Improving the extended sagittal segmentation algorithm: Our extended sagittal segmentation algorithm requires the aid of Stradwin to visualize the volume data and to select 14 initialization points. Developing a custom made visualization module that can be integrated into the segmentation algorithm offers a smoother and more unified work flow. As mentioned in Section 3.5.2, one of the shortcoming of the algorithm is its inability to contour the prostate tailings at its bottom sides when significant prostate “squishing” is introduced by the probe. Placing extra initialization points by using the transverse views may improve the result. 101 5.2. Pre-operative to Intra-operative Registration of Prostate Volumes • Improving the registration accuracy: The major issue of our ellipsoid fitting-based registration is its sensitivity to user initialization. Possible solutions that can be studied include using a fully automatic segmentation algorithm such as [69] and [70] and using non-ellipsoid fitting models (e.g. super ellipsoids). A completely different approach is to use a model-less registration approach such as ICP which was briefly studied in this thesis. In addition, deformable registration can also be explored. • A more extensive study of the extended sagittal segmentation algorithm: We validated the extended sagittal segmentation algorithm on patient volume data by comparing it with a transverse view-based segmentation algorithm. However, the gold standard of segmentation is manual contouring by the oncologist. More proper validation can be explored by designing an extensive in-vivo comparison study where the oncologist manually contours the sagittal images - not the transverse ones - and compare the results to our extended sagittal segmentation algorithm. This study would also serve as a usability study of the algorithm. 102 Bibliography [1] C. C. Society, “Canadian cancer statistics 2009,” Statistics Canada, Tech. Rep., 2009. [2] R. Kirby, “Treatment options for early prostate cancer,” Urology, vol. 52, no. 6, pp. 948 – 962, 1998. [3] A. B. Jani and S. Hellman, “Early prostate cancer: clinical decisionmaking,” The Lancet, vol. 361, no. 9362, pp. 1045 – 1053, 2003. [4] R. M. Benoit, M. J. Naslund, and J. K. Cohen, “Complications after prostate brachytherapy in the medicare population,” Urology, vol. 55, no. 1, pp. 91 – 96, 2000. [5] P. Wust, D. von Borczyskowski, T. Henkel, C. Rosner, R. Graf, W. Tilly, V. Budach, R. Felix, and F. Kahmann, “Clinical and physical determinants for toxicity of 125-i seed prostate brachytherapy,” Radiotherapy and Oncology, vol. 73, no. 1, pp. 39 – 48, 2004. [6] Y. Yu, L. L. Anderson, Z. Li, D. E. Mellenberg, R. Nath, M. C. Schell, F. M. Waterman, A. Wu, and J. C. Blasko, “Permanent prostate seed implant brachytherapy: Report of the american association of physicists in medicine task group no. 64,” Medical Physics, vol. 26, no. 10, pp. 2054–2076, 1999. [7] C. M. Solutions, “Exii stepper,” Online Catalogue. R lp,” Online Catalogue. [8] ——, “Micro-touch ° [9] V. M. Systems, “Variseed,” Online Catalogue. [10] P. L. Roberson, V. Narayana, D. L. McShan, R. J. Winfield, and P. W. McLaughlin, “Source placement error for permanent implant of the prostate,” Medical Physics, vol. 24, no. 2, pp. 251–257, 1997. [11] J. Dawson, T. Wu, T. Roy, J. Gu, and H. Kim, “Dose effects of seeds placement deviations from pre-planned positions in ultrasound guided 103 Bibliography prostate implants,” Radiotherapy and Oncology, vol. 32, no. 3, pp. 268 – 270, 1994. [12] M. Glegg, D. Dodds, and G. M. Baxter, “Prostate brachytherapy: An option in the treatment of early stage prostate cancer,” Urooncology, vol. 4, no. 2, pp. 53 – 69, 2004. [13] E. Dehghan and S. Salcudean, “Needle insertion point and orientation optimization in non-linear tissue with application to brachytherapy,” in Robotics and Automation, 2007 IEEE International Conference on, April 2007, pp. 2267–2272. [14] V. Lagerburg, M. A. Moerland, J. J. Lagendijk, and J. J. Battermann, “Measurement of prostate rotation during insertion of needles for brachytherapy,” Radiotherapy and Oncology, vol. 77, no. 3, pp. 318 – 323, 2005. [15] J. Sherman, T. Podder, V. Misic, L. Fu, D. Fuller, B. Winey, E. Messing, D. Rubens, J. Strang, R. Brasacchio, and Y. Yu, “Efficacy of prostate stabilizing techniques during brachytherapy procedure,” in Engineering in Medicine and Biology Society, 2006. EMBS ’06. 28th Annual International Conference of the IEEE, 30 2006-Sept. 3 2006, pp. 563–566. [16] R. Taschereau, J. Pouliot, J. Roy, and D. Tremblay, “Seed misplacement and stabilizing needles in transperineal permanent prostate implants,” Radiotherapy and Oncology, vol. 55, no. 1, pp. 59 – 63, 2000. [17] K. M. Langen and D. T. L. Jones, “Organ motion and its management,” International Journal of Radiation Oncology*Biology*Physics, vol. 50, no. 1, pp. 265 – 278, 2001. [18] C. F. Serago, S. J. Chungbin, S. J. Buskirk, S. A. Vora, M. P. McLaughlin, and G. A. Ezzell, “Initial experience with ultrasound localization for positioning prostate cancer patients for external beam radiation therapy,” International Journal of Radiation Oncology*Biology*Physics, vol. 51, no. 3, Supplement 1, pp. 94 – 94, 2001. [19] X. Sheng, J. Kruecker, P. Guion, N. Glossop, Z. Neeman, P. Choyke, A. Singh, and B. Wood, “Closed-loop control in fused mr-trus imageguided prostate biopsy,” Medical Image Computing and ComputerAssisted Intervention MICCAI 2007, vol. 24, no. 2, pp. 128–135, 2007. 104 Bibliography [20] T. C. Green and M. J. Horzewski, “Pivoting needle template apparatus for brachytherapy treatment of prostate disease and methods of use,” U.S. Patent 6 398 711B1, Jun. 4, 2002. [21] M. Hogendijk and R. P. Boucher, “Positioning needle guide for brachytherapy treatment of prostate disease and methods of use,” U.S. Patent Request 10/856 074, May 28, 2004. [22] C. Schneider, A. Okamura, and G. Fichtinger, “A robotic system for transrectal needle insertion into the prostate with integrated ultrasound,” in Robotics and Automation, 2004. Proceedings. ICRA ’04. 2004 IEEE International Conference on, vol. 1, April-1 May 2004, pp. 365–370 Vol.1. [23] H. Bassan, T. Hayes, R. Patel, and M. Moallem, “A novel manipulator for 3d ultrasound guided percutaneous needle insertion,” in Robotics and Automation, 2007 IEEE International Conference on, April 2007, pp. 617–622. [24] Y. Yu, T. Podder, Y. Zhang, W.-S. Ng, V. Misic, J. Sherman, L. Fu, D. Fuller, E. Messing, D. Rubens, J. Strang, and R. Brasacchio, Medical Image Computing and Computer-Assisted Intervention MICCAI 2006, ser. Lecture Notes in Computer Science. Springer Berlin / Heidelberg, September 2006. [25] Z. Wei, M. Ding, D. Downey, and A. Fenster, “Dynamic intraoperative prostate brachytherapy using 3d trus guidance with robot assistance,” in Engineering in Medicine and Biology Society, 2005. IEEE-EMBS 2005. 27th Annual International Conference of the, Jan. 2005, pp. 7429– 7432. [26] Z. Wei, G. Wan, L. Gardi, G. Mills, D. Downey, and A. Fenster, “Robotassisted 3d-trus guided prostate brachytherapy: System integration and validation,” Medical Physics, vol. 31, no. 3, pp. 539–548, 2004. [27] J. Bax, L. Gardi, J. Montreuil, D. Smith, and A. Fenster, “A compact robotic apparatus and method for 3-d ultrasound guided prostate therapy,” K. R. Cleary and M. I. Miga, Eds., vol. 6509, no. 1. SPIE, 2007, p. 65092H. [28] D. Stoianovici, K. Cleary, A. Patriciu, D. Mazilu, A. Stanimir, N. Craciunoiu, V. Watson, and L. Kavoussi, “Acubot: a robot for radiologi105 Bibliography cal interventions,” Robotics and Automation, IEEE Transactions on, vol. 19, no. 5, pp. 927–930, Oct. 2003. [29] G. Fichtinger, E. C. Burdette, A. Tanacs, A. Patriciu, D. Mazilu, L. L. Whitcomb, and D. Stoianovici, “Robotically assisted prostate brachytherapy with transrectal ultrasound guidance–phantom experiments,,” Brachytherapy, vol. 5, no. 1, pp. 14 – 26, 2006. [30] G. Fichtinger, J. P. Fiene, C. W. Kennedy, G. Kronreif, I. Iordachita, D. Y. Song, E. C. Burdette, and P. Kazanzides, “Robotic assistance for ultrasound-guided prostate brachytherapy,” Medical Image Analysis, vol. 12, no. 5, pp. 535 – 545, 2008, special issue on the 10th international conference on medical imaging and computer assisted intervention MICCAI 2007. [31] S. DiMaio, G. Fischer, S. Maker, N. Hata, I. Iordachita, C. Tempany, R. Kikinis, and G. Fichtinger, “A system for mri-guided prostate interventions,” in Biomedical Robotics and Biomechatronics, 2006. BioRob 2006. The First IEEE/RAS-EMBS International Conference on, Feb. 2006, pp. 68–73. [32] L. Phee, D. Xiao, J. Yuen, C. F. Chan, H. Ho, C. H. Thng, C. Cheng, and W. S. Ng, “Ultrasound guided robotic system for transperineal biopsy of the prostate,” in Robotics and Automation, 2005. ICRA 2005. Proceedings of the 2005 IEEE International Conference on, April 2005, pp. 1315–1320. [33] J. Hong, T. Dohi, M. Hashizume, K. Konishi, and N. Hata, “An ultrasound-driven needle-insertion robot for percutaneous cholecystostomy,” Physics in Medicine and Biology, vol. 49, no. 3, pp. 441–455, 2004. [34] J. Kettenbach, G. Kronreif, M. Figl, M. Frst, W. Birkfellner, R. Hanel, and H. Bergmann, “Robot-assisted biopsy using ultrasound guidance: initial results from in vitro tests.” European Radiology, vol. 15, no. 4, pp. 765 – 771, 2005. [35] J. Wu, T. Haycocks, H. Alasti, G. Ottewell, N. Middlemiss, M. Abdolell, P. Warde, A. Toi, and C. Catton, “Positioning errors and prostate motion during conformal prostate radiotherapy using on-line isocentre set-up verification and implanted prostate markers,” Radiotherapy and Oncology, vol. 61, no. 2, pp. 127 – 133, 2001. 106 Bibliography [36] J. M. Crook, Y. Raymond, D. Salhani, H. Yang, and B. Esche, “Prostate motion during standard radiotherapy as assessed by fiducial markers,” Radiotherapy and Oncology, vol. 37, no. 1, pp. 35 – 42, 1995. [37] S. Shimizu, H. Shirato, K. Kitamura, N. Shinohara, T. Harabayashi, T. Tsukamoto, T. Koyanagi, and K. Miyasaka, “Use of an implanted marker and real-time tracking of the marker for the positioning of prostate and bladder cancers,” International Journal of Radiation Oncology*Biology*Physics, vol. 48, no. 5, pp. 1591 – 1597, 2000. [38] T. R. Willoughby, P. A. Kupelian, J. Pouliot, K. Shinohara, M. Aubin, M. R. III, L. L. Skrumeda, J. M. Balter, D. W. Litzenberg, S. W. Hadley, J. T. Wei, and H. M. Sandler, “Target localization and realtime tracking using the calypso 4d localization system in patients with localized prostate cancer,” International Journal of Radiation Oncology*Biology*Physics, vol. 65, no. 2, pp. 528 – 534, 2006. [39] K. Kitamura, H. Shirato, S. Shimizu, N. Shinohara, T. Harabayashi, T. Shimizu, Y. Kodama, H. Endo, R. Onimaru, S. Nishioka, H. Aoyama, K. Tsuchiya, and K. Miyasaka, “Registration accuracy and possible migration of internal fiducial gold marker implanted in prostate and liver treated with real-time tumor-tracking radiation therapy (rtrt),” Radiotherapy and Oncology, vol. 62, no. 3, pp. 275 – 281, 2002. [40] M. J. Ghilezan, D. A. Jaffray, J. H. Siewerdsen, M. V. Herk, A. Shetty, M. B. Sharpe, S. Z. Jafri, F. A. Vicini, R. C. Matter, D. S. Brabbins, and A. A. Martinez, “Prostate gland motion assessed with cine-magnetic resonance imaging (cine-mri),” International Journal of Radiation Oncology*Biology*Physics, vol. 62, no. 2, pp. 406 – 417, 2005. [41] A. R. Padhani, V. S. Khoo, J. Suckling, J. E. Husband, M. O. Leach, and D. P. Dearnaley, “Evaluating the effect of rectal distension and rectal movement on prostate gland position using cine mri,” International Journal of Radiation Oncology*Biology*Physics, vol. 44, no. 3, pp. 525 – 533, 1999. [42] C. J. Beard, P. Kijewski, M. Bussire, R. Gelman, D. Gladstone, K. Shaffer, M. Plunkett, P. Costello, and C. N. Coleman, “Analysis of prostate and seminal vesicle motion: Implications for treatment planning,” International Journal of Radiation Oncology*Biology*Physics, vol. 34, no. 2, pp. 451 – 458, 1996. 107 Bibliography [43] J. Lattanzi, S. McNeely, A. Hanlon, I. Das, T. E. Schultheiss, and G. E. Hanks, “Daily ct localization for correcting portal errors in the treatment of prostate cancer,” International Journal of Radiation Oncology*Biology*Physics, vol. 41, no. 5, pp. 1079 – 1086, 1998. [44] J. C. Roeske, J. D. Forman, C. F. Mesina, T. He, C. A. Pelizzari, E. Fontenla, S. Vijayakumar, and G. T. Y. Chen, “Evaluation of changes in the size and location of the prostate, seminal vesicles, bladder, and rectum during a course of external beam radiation therapy,” International Journal of Radiation Oncology*Biology*Physics, vol. 33, no. 5, pp. 1321 – 1329, 1995, implementation of Three Dimensional Conformal Radiotherapy. [45] M. van Herk, A. Bruce, A. P. G. Kroes, T. Shouman, A. Touw, and J. V. Lebesque, “Quantification of organ motion during conformal radiotherapy of the prostate by three dimensional image registration,” International Journal of Radiation Oncology*Biology*Physics, vol. 33, no. 5, pp. 1311 – 1320, 1995, implementation of Three Dimensional Conformal Radiotherapy. [46] J. Lattanzi, S. McNeeley, S. Donelly, E. Palacio, A. Hanlon, and T. E. S. amd Gerald E. Hanks, “Ultrasound-based stereotactic guidance in prostate cancerquantification of organ motion and set-up errors in external beam radiation,” Computer Aided Surgery, vol. 5, no. 4, pp. 289 – 295, 2000. [47] F. Trichter and R. D. Ennis, “Prostate localization using transabdominal ultrasound imaging,” International Journal of Radiation Oncology*Biology*Physics, vol. 56, no. 5, pp. 1225 – 1233, 2003. [48] A. Krupa, G. Fichtinger, and G. D. Hager, “Real-time motion stabilization with b-mode ultrasound using image speckle information and visual servoing,” Int. J. Rob. Res., vol. 28, no. 10, pp. 1334–1354, 2009. [49] E. R. Jr., D. Skarecky, N. Narula, and T. E. Ahlering, “Prostate volume estimation using the ellipsoid formula consistently underestimates actual gland size,” The Journal of Urology, vol. 179, no. 2, pp. 501 – 503, 2008. [50] Q. Li and J. Griffiths, “Least squares ellipsoid specific fitting,” in Geometric Modeling and Processing, 2004. Proceedings, 2004, pp. 335–340. 108 Bibliography [51] S. Badiei, S. E. Salcudean, J. Varah, and W. J. Morris, “Prostate segmentation in 2d ultrasound images using image warping and ellipse fitting,” in MICCAI (2), 2006, pp. 17–24. [52] L. Gong, S. Pathak, D. Haynor, P. Cho, and Y. Kim, “Parametric shape modeling using deformable superellipses for prostate segmentation,” Medical Imaging, IEEE Transactions on, vol. 23, no. 3, pp. 340 –349, march 2004. [53] J. A. Noble and D. Boukerroui, “Ultrasound image segmentation: a survey,” Medical Imaging, IEEE Transactions on, vol. 25, no. 8, pp. 987–1010, July 2006. [54] S. S. Mahdavi, W. J. Morris, I. Spadinger, N. Chng, O. Goksel, and S. E. Salcudean, “3d prostate segmentation in ultrasound images based on tapered and deformed ellipsoids,” in MICCAI ’09: Proceedings of the 12th International Conference on Medical Image Computing and Computer-Assisted Intervention. Berlin, Heidelberg: Springer-Verlag, 2009, pp. 960–967. [55] Faulhaber, “Mcdc 3006s user manual,” Online Manual. [56] S. Salcudean, T. Prananta, W. Morris, and I. Spadinger, “A robotic needle guide for prostate brachytherapy,” in Robotics and Automation, 2008. ICRA 2008. IEEE International Conference on, May 2008, pp. 2975–2981. [57] X. Wen, “Towards ultrasound-based intraoperative dosimetry for prostate brachytherapy,” Ph.D. dissertation, The University of British Columbia, 2010. [58] P. Abolmaesumi and M. R. Sirouspour, “An interacting multiple model probabilistic data association filter for cavity boundary extraction from ultrasound images,” IEEE Trans. Med. Imaging, vol. 23, no. 6, pp. 772–784, 2004. [59] A. Fitzgibbon, M. Pilu, and R. Fisher, “Direct least square fitting of ellipses,” Pattern Analysis and Machine Intelligence, IEEE Transactions on, vol. 21, no. 5, pp. 476–480, May 1999. [60] VERMON, “Bip 6.5/r10/2 x 128 - 472 probe specifications,” Specifications. 109 Bibliography [61] R. Halir and J. Flusser, “Numerically stable direct least squares fitting of ellipses,” 1998. [62] R. W. Prager, A. H. Gee, L. Berman, R. Prager, A. Gee, L. Berman, and C. C. Pz, “Stradx: Real-time acquisition and visualisation of freehand 3d ultrasound,” Medical Image Analysis, vol. 3, pp. 129–140, 1998. [63] J. J. Craig, Introduction to robotics: mechanics and control. Saddle River, NJ, USA: Pearson/Prentice-Hall, 2005. Upper [64] P. R.W., R. R.N., G. A.H., and B. L., “Rapid calibration for 3-d freehand ultrasound - technique, possibilities and limitations,” Ultrasound in Medicine and Biology, vol. 24, pp. 855–869(15), July 1998. [65] P. Besl and H. McKay, “A method for registration of 3-d shapes,” Pattern Analysis and Machine Intelligence, IEEE Transactions on, vol. 14, no. 2, pp. 239 –256, feb 1992. [66] D. Brujic and M. Ristic, “Analysis of free form surface registration,” vol. 1, sep. 1996, pp. 393 –396 vol.2. [67] J. M. Phillips, “Outlier robust icp for minimizing fractional rmsd,” ArXiv, vol. 606098, 2006. [68] E. Trucco, A. Fusiello, and V. Roberto, “Robust motion and correspondence of noisy 3-d point sets with missing data,” Pattern Recognition Letters, vol. 20, no. 9, pp. 889 – 898, 1999. http://www.sciencedirect.com/science/article/ [Online]. Available: B6V15-3XFTF4X-3/2/9d2386986ab0d48defc64f4ffaf9a941 [69] W. D. Richard and C. G. Keen, “Automated texture-based segmentation of ultrasound images of the prostate,” Computerized Medical Imaging and Graphics, vol. 20, no. 3, pp. 131 – 140, 1996. [70] A. Ghanei, H. Soltanian-Zadeh, A. Ratkewicz, and F.-F. Yin, “A threedimensional deformable model for segmentation of human prostate from ultrasound images,” Medical Physics, vol. 28, no. 10, pp. 2147–2153, 2001. [71] Z. Zhang, Z. Zhang, P. Robotique, and P. Robotvis, “Parameter estimation techniques: A tutorial with application to conic fitting,” Image and Vision Computing, vol. 15, pp. 59–76, 1997. 110 [72] Y. Bar-Shalom and T. E. Fortmann, Tracking and data association, ser. Mathematics in Science and Engineering. San Diego, CA, USA: Academic Press Professional, Inc., 1987, vol. 179. [73] PIC-Design, “Lead screws and nuts,” Online Catalogue. 111 Appendix A Ellipse Fitting − → − → − → The 2D algebraic distance F( X , P ) of a point X [x y]T to a conic with → − → − − → algebraic parameter P = [a b c d e f ]T is described by equation A.1. F( X , P ) = 0 when the point lies exactly on the conic’s boundary. The shape of the − → conic is an ellipse if P follows inequality (A.2) [71]. → − − → F ( X , P ) = ax2 + 2bxy + cy 2 + 2dx + 2ey + f 2 4ac − b > 0 (A.1) (A.2) − → By introducing an auxiliary vector L = [x2 xy y 2 x y 1], equation A.1 and constraint A.2 can be rewritten in matrix form as: − → − → − →→ − F (X , P ) = L P (A.3) 0 0 2 0 0 0 0 −1 0 0 0 0 → − →T 2 0 0 0 0 0 − P 0 0 0 0 0 0 P > 0 0 0 0 0 0 0 0 0 0 0 0 0 − →T → − P C P > 0 (A.4) The “Constraint Matrix” C is symmetric and singular. A.1 Direct Least Squares Fitting of Ellipses Fitzgibbon et.al. introduced two ideas to the ellipse fitting problem. First, − → − → fitting a set of n points X i=1..n into an ellipse is to optimize P values by minimizing the sum of squared algebraic distances between the points and the ellipse. Second, constraint (A.2) can be reformulated as an equation. Combining the two ideas, they developed a fast and robust ellipse fitting algorithm. 112 A.1. Direct Least Squares Fitting of Ellipses Referring to equation A.1, the sum of square algebraic distance is formulated as: n X → − → − − →2 Q( P ) = F ( L i, P ) i=1 →2 →2 → −T − →2 − → − → − − = (L1 P ) + (L2 P ) + ... + (Ln P ) ¯¯ − ¯¯ → ¯¯ L 1 → ¯¯ − ¯¯ L 2 ¯ ¯ = ¯¯ . ¯¯ .. ¯¯ → ¯¯ − Ln = ¯¯2 ¯¯ ¯¯ ¯¯ − →¯¯¯¯ → −T P ¯¯ = P ¯¯ ¯¯ ¯¯ − → L1 → − L2 . . . − → Ln T − → L1 → − L2 . . . − → Ln − → P → −T − → T P D | {zD} P (A.5) S D and S are referred as the “Design Matrix” and “Scatter Matrix” respectively. As shown below, the Scatter Matrix is symmetric - an important property for subsequent explanations. S= x21 x22 x1 y1 x2 y2 2 y1 y22 = x1 x2 y1 y2 1 1 = ··· ··· ··· ··· ··· ··· DT D x2n x21 x1 y1 y12 x1 y1 1 x n yn 2 x2 x2 y2 y22 x2 y2 1 yn2 .. .. .. .. .. .. xn . . . . . . yn x2 xn yn y 2 xn yn 1 n n 1 (A.6) Pn 4 n Pn i=1 x 3 i=1 xn yn Pn 2 2 i=1 xn yn P n 3 Pn i=1 xn 2 i=1 xn yn P n 2 i=1 xn | Pn 3 Pni=1 x2n yn2 xn yn Pi=1 n 3 Pni=1 x2n yn Pni=1 xn yn2 xn yn Pi=1 n i=1 xn yn Pn Pn Pn Pn 2 2 3 2 2 n i=1 xn yn P i=1 xn Pni=1 xn yn3 Pn i=1 x P n n 2 2 i=1 xn yn Pi=1 xn yn i=1 xn yn i=1 xn yn P P P n n n n 4 2 3 2 y x y y y i=1 n n P i=1 n Pni=1 n P Pn i=1 n 2 n n 2 xn x y x x y n n n n P i=1 n i=1 Pi=1 P Pni=1 n n n 2 3 y x y y x y n n n n n n i=1 i=1 i=1 i=1 P P Pn n n 2 n i=1 yn i=1 xn i=1 yn {z } Symmetric (A.7) 113 A.2. Numerically Stable Least Squares Fitting The ellipse constraint (inequality (A.2)) tells us that the left hand side must be a positive number: any positive number. Hence, we can set the ellipse constrain to be 4ac − b2 = k. By scaling down k to 1 and proper scaling of a,b, and c, the ellipse constraint is re-defined as an equality: 4ac − b2 = 1 − →T − → P C P =1 (A.8) The ellipse fitting problem boils down to minimizing the scalar function Q with constraint (A.8). Applying the Lagrange multiplier concept and taking advantage the symmetric property of C and S simplify the problem further: ∂Q − → ∂P − →T − → ∂( P S P ) → − ∂P − → 2S P − → SP − →T → − ∂( P C P − 1) = λ → − ∂P − → = λ 2C P − → = λ 2C P → − = λ CP (A.9) The solution amounts to finding the eigenvalues and eigenvectors pairs of the Scatter Matrix. Fitzgibbon et.al. proved that the only eigenvector corresponding to the positive eigenvalue gives the best ellipse fit. A.2 Numerically Stable Least Squares Fitting Fitzgibbon et.al.’s algorithm is, however, numerically unstable [61]. As noted by Halir and Flusser, the Constraint Matrix is singular and when the set of points are very close to an ellipse, the Scatter Matrix is close to singular: in occasions, leading to the wrong sets of eigenvalues and eigenvectors. To circumvent the problem Halir and Flusser developed a solution by − → decomposing P , the Design matrix, the Scatter matrix, and the Constraint matrix: 114 A.2. Numerically Stable Least Squares Fitting x21 x22 .. . D = x1 y1 x2 y2 .. . y12 y22 .. . x1 x2 .. . y1 y2 .. . 1 1 .. . x2n xn yn yn2 xn yn 1 £ ¤ D = D1(n×3) D2(n×3) (A.10) T S = D · TD¸ ¤ D1 £ D1 D2 S = T D2 ¸ · T D1 D1 D1T D2 S = D2T D1 D2T D2 · ¸ S1(3×3) S2(n×3) S = S2T S3(3×3) · C = C1 0 0 0 − → P = a b c d e f ¸ (A.11) 0 0 2 where C1 = 0 −1 0 2 0 0 (A.12) · ¸ = P1(3×1) P2(3×1) (A.13) The eigensystem (A.9) can then be reformulated as: · S1 S2T S2 S3 ¸ · P1 P2 ¸ · = λ C1 0 0 0 ¸ · P1 P2 ¸ (A.14) Which can be split into two equations: S1 P1 + S2 P2 S2T P1 = λC1 P1 (A.15) + S3 P2 = 0 P2 = −S3−1 S2T P1 (A.16) 115 A.3. Converting Algebraic Parameters to Geometric Parameters Substituting equation A.16 into (A.15) and taking the inverse of C1 , we have the following eigensystem: C1−1 (S1 − S2 S3−1 S2T ) P1 = λP1 (A.17) The eigenvector corresponding to the positive eigenvalue of equation A.17 yields the optimum value of P1 . The optimum solution of the ellipse fit is the combination of P1 and P2 - calculated using (A.16). Since C1 and (S1 − S2 S3−1 S2T ) are no longer singular, the eigensystem is more numerically stable. A.3 Converting Algebraic Parameters to Geometric Parameters − → The parameter P offers little insight to the geometric properties of the − → ellipse. An ellipse is more often described by its geometric parameters Q [xc yc ra rb θ]T . Where (xc , yc ) is the coordinate location of the ellipse center, ra is the length of the major axis, rb is the length of the minor axis, and θ is the tilt of the ellipse with respect to the X + axis. − → → − To convert P to Q , we define three coordinate frames {1,2,3} - see Figure A.1. Coordinate frame {1} is the stationary reference frame. Frame {2} is frame {1} rotated counter clockwise by angle θ. Frame {3} is frame {2} translated by (1 xc 1 yc ) and is attached to the ellipse. The following geometric transformation can be defined: ¸· ¸ ·2 ¸ · Cos(θ) Sin(θ) 1 x x 1y 2 y = −Sin(θ) Cos(θ) (A.18) ¸· ¸ · ¸· ¸ ·3 ¸ · Cos(θ) Sin(θ) 1 x Cos(θ) Sin(θ) 1 xc x 3 y = −Sin(θ) Cos(θ) 1 y − −Sin(θ) Cos(θ) 1y c (A.19) In frame {3}, the ellipse is unrotated and untranslated. Its equation can simply be written as: 3 x2 3y2 =1 (A.20) + ra2 rb2 The parametric equation of the ellipse in the reference frame (frame1) is: a 1 x 2 + b 1 x 1 y + c 1 y 2 + d 1 x + e1 y + f = 0 (A.21) 116 A.3. Converting Algebraic Parameters to Geometric Parameters Figure A.1: Different coordinate frames used to derive the ellipse parameter conversion 117 A.3. Converting Algebraic Parameters to Geometric Parameters Substituting equation A.19 into equation A.21 and grouping the terms according to the variables (x2 , xy, y 2 , x, y) we come up with the following set of equations: 2θ rb2 Cos2 θ ¡+ ra2 Sin ¢ a Sin2θ rb2 − ra2 b 2 2 2 2 r Sin θ + r Cos θ a b c 1 2 2 2 2 1 2 2 = −2 xc (ra Sin θ + rb Cos θ) + yc Sin2θ (ra − rb ) d 1 2 2 2 2 1 2 2 −2 y (r Sin θ + r Cos θ) + x Sin2θ (r − r ) c c a a b b e ¡ 2 ¢ 2 1 2 2 1 1 2 2 1 2 r Cos θ x + r Sin2θ x y + r Sin θ y c c c c b b ¡ 2 b 2 1 2 ¢ f 2 1 2 1 2 1 2 2 2 + ra Sin θ xc − ra Sin2θ xc yc + ra Cos θ yc − ra rb We can find the angle θ by from a, b, and c in equation A.22. a − c = (rb2 Cos2 θ + ra2 Sin2 θ) − (rb2 Sin2 θ + ra2 Cos2 θ) = Cos2θ(rb2 − ra2 ) ¡ ¢ b = Sin2θ rb2 − ra2 b θ = arctan( ) a−c (A.22) (A.23) (A.24) In frame {2}, the ellipse is translated but not rotated. The center of the ellipse is located at (2 xc 2 yc ). The equation of the ellipse can be described as: (2 x −2 xc )2 (2 y −2 yc )2 + =1 2 ra rb2 rb2 2 x2 + ra2 2 y 2 − 2rb2 2 xc 2 x − 2ra2 2 yc 2 y + (rb2 2 x2c + ra2 2 yc2 − ra2 rb2 ) = 0 (A.25) The variables ra and rb are the length of the major and minor axes of the ellipse. As rotation and translation do not affect geometric length, they are − → same as ra and rb in equation A.22. Let’s define P ’ as: 0 rb2 a b0 0 0 2 c r a = 2 2 d0 −2r x c b 2 2 e0 −2ra yc rb2 2 xc + ra2 2 yc − ra2 rb2 f0 (A.26) 118 A.3. Converting Algebraic Parameters to Geometric Parameters Substituting equation A.18 into equation A.21, we have the following equation: (a Cos2 θ + b SinθCosθ + c Sin2 θ) 2 x2 +(−2a SinθCosθ + b Cos2 θ − b Sin2 θ + 2c SinθCosθ) 2 x2 y +(a Sin2 θ − b SinθCosθ + c Cos2 θ) 2 y 2 +(d Cosθ + e Sinθ) 2 x + (−d Sinθ + e Cosθ) 3 y + f = 0 (A.27) We can solve for (2 xc 2 yc ) and (1 xc 1 yc ) from equations A.26 and A.27: ·2 ¸ · d0 ¸ xc − 2a0 = 2y e0 − 2c c 0 " # ·2 ¸ d Cosθ+e Sinθ − 2(a Cos2 θ+b 2 xc SinθCosθ+c Sin θ) = 2y Sinθ+e Cosθ − 2(a Sin2−d c θ−b SinθCosθ+c Cos2 θ) · ¸· ¸ ·1 ¸ Cos(θ) −Sin(θ) 2 xc xc = 2y 1y Sin(θ) Cos(θ) c c (A.28) The length of the ellipse’s radii can be found using the following equation: s −P + P402 + P502 6 4P 0 4P30 · ¸ 1 P10 ra s = P402 P502 rb −P6 + 4P + 0 4P 0 1 P10 3 (A.29) 119 Appendix B Ellipsoid Fitting → − − → → − The algebraic distance F( X , P ) of a point X [x, y, z]T to a 3D conic with − → parameter P = [a b c d e f g p q r s]T is defined as: → − − → F ( X , P ) = ax2 + by 2 + cz 2 + 2dxy + 2exz + 2f yz + 2px + 2qy + 2rz + s (B.1) − →− → Similar to the 2D conic equation, F( X , P ) = 0 if the point lies on the conic’s surface. The conic is an ellipsoid if J > 0 and I × K > 0 with I, J, and K defined as follows: I =a+b+c (B.2) − 2 (B.3) 2 J = ab + ac + ac − f e − d ¯ ¯ ¯a d e ¯ ¯ ¯ K = ¯¯d b f ¯¯ ¯e f c ¯ (B.4) − → Similar to the ellipse fitting problem, fitting a set of n points X i=1..n into an ellipsoid amounts to minimizing the sum of squared of the algebraic distance between the points to the ellipsoid while adhering to constraints B.2 - B.4. B.1 Direct Least Square Fitting of Ellipsoids Li and Griffiths [50] simplified ellipsoid the constraints into equation B.5. They proved that any conic that satisfies the equation is an ellipsoid and developed a robust ellipsoid fitting algorithm. 4J − I 2 = 0 (B.5) 120 B.2. Converting Algebraic Parameters to Geometric Parameters Equation B.5 can be rewritten in matrix form: − →T − → P CP = 0 · ¸ − →T → C1 0(6×4) − P P =0 0(4×6) 0(4×4) −1 1 1 0 0 0 1 −1 1 0 0 0 1 1 −1 0 0 0 C1 = 0 0 0 −4 0 0 0 0 0 0 −4 0 0 0 0 0 0 −4 (B.6) (B.7) (B.8) − → By redefining the auxillary vector L , redecomposing the Scatter Matrix, − → and redecomposing the vectot P , the solution of the ellipsoid fitting follows equation A.15) - A.17). £ ¤ − → L i=1..n = x2i yi2 zi2 xi yi xi zi yi zi xi yi zi 1 (B.9) # " T S1(6×6) S2(6×4) S = (B.10) T S2 S3(4×4) · ¸ − → P1(1×6) P = (B.11) P2(1×4)) The eigenvector corresponding to the positive eigenvalue of equation A.17 yields the optimum value of P1 . Although constraint B.5 guarantees an ellipsoid solution, it does not cover the whole set of the ellipsoids. In particular, Li and Griffiths cautions that the constraint cannot characterize ellipsoids that are “long-thin or compressed”. The ellipsoid fitting algorithm is used to model the prostate and we assume that normal prostate are not “long-thin or compressed”. B.2 Converting Algebraic Parameters to Geometric Parameters Similar to the ellipse fitting problem, we wish to express the ellipsoid using its geometric parameters: center location (xc , yc , zc ), radii length (rx , ry , rz ), orientation about the reference frame’s principal axis (Ax , Ay , Az ). We define the coordinate frames: {1} and {2} as illustrated in Figure B.1. Frame {1} is the reference frame. Frame {2} is frame {1} rotated, in 121 B.2. Converting Algebraic Parameters to Geometric Parameters order, by angle Ax , Ay , Az around its own x,y,and z-axis and then translated by (1 xc ,1 yc ,1 zc )T . Figure B.1: Coordinate frames for ellipsoid parameter conversion −→ For point 1 X on the conic’s surface, the parametric equation can be normalized into: a01 x2 + b01 y 2 + c01 z 2 + 2d01 x1 y + 2e01 x1 z + 2f 01 y 1 z +2p01 x + 2q 01 y + 2r01 z + 1 = 0 Equation B.12 0 T a − → 1 X d0 e0 can be expressed in matrix form as: d0 e0 −→ £ ¤ −→ b0 f 0 1 X + 2p0 2q 0 2r0 1 X + 1 = 0 f 0 c0 − →T −→ −→ 1 X A1 X + B 1 X + 1 = 0 (B.12) (B.13) The relationship between points in frame {1} and {2} can be described 122 B.2. Converting Algebraic Parameters to Geometric Parameters as: − → 1X = 1 R2 −→ −→ R2 2 X + 1 X c cy cz sx sy cz − cx sz cx sy cz + sx sz = cy sz sx sy sz + cx cz cx sy sz − sx cz −sy sx cy cx cy 1 (B.14) (B.15) R is the rotation matrix following the Euler yaw-pitch-roll convention. Where sx is sin(Ax ), cx is cos(Ax ), sy is sin(Ay ), and so on. Ax , Ay , Az are the yaw angle, pitch angle, and roll angle. They are the rotation angles about the x, y, and z axes of the reference frame1. We substitute equation B.14 into equation B.13 and normalize the result: −→ −−→T −→ (1 R2 2 X + 1 Xc ) A (1 R2 2 X + 1 Xc ) + −→ −−→ B (1 R2 2 X + 1 Xc ) + 1 −→ −→T −−→ −−→ −→ − →T 2 X 1 RT A 1 R 2 X + 2 X (1 RT A 1 X ) + (1 RT A 1 X )T 2 X + c c 2 2 2 2 T − → − − → − − → −−→ (B 1 R2 2 X) + (1 Xc A 1 Xc + B 1 Xc + 1) −→ −−→ −→ − →T 2 X 1 RT A 1 R 2 X + 2((1 RT A 1 X )T + B 1 R ) 2 X c 2 2 2 2 −−→T −−→ −−→ +(1 Xc A 1 Xc + B 1 Xc + 1) 1 RT A 1 R − →T − → 2 2 2X 2X + −−→T −−→ −−→ −(1 Xc A 1 Xc + B 1 Xc + 1) −−→ → (21 XcT A 1 R2 + B 1 R2 ) − 2X − 1 T − − → − − → − − → −(1 Xc A 1 Xc + B 1 Xc + 1) = 0 = 0 = 0 = 0 (B.16) In coordinate frame {2}, the ellipsoid is unrotated and untranslated. Its parametric equation is as follows: 2 x2 1 2 − → ra 2X 0 0 ra2 0 1 rb2 0 + 2y2 rb2 + 2z2 rc2 =1 0 −→ 0 2X − 1 = 0 (B.17) 1 rc2 Comparing equation B.16 and equation B.17, we can find 1 Xc (the center of the ellipsoid in the reference frame): 123 B.3. Generating Ellipses from an Ellipsoid −−→ 21 XcT A 1 R2 + B 1 R2 = 0 1 1 Xc = − A−1 B 2 (B.18) (B.19) The radii lengths of the ellipsoid and its orientation around the reference frame’s x, y, and z-axis can be solved using the eigen values and eigenvectors of the normalized A matrix: 1 0 0 ra2 A 1 1 1 T R2 R2 = 0 r2 0 T − − − → − → − − → b −(1 Xc A1 Xc + B 1 Xc + 1) 0 0 r12 c 1 0 0 ra2 1 1 T R2 D 1 R2 = 0 r2 0 (B.20) 0 b 0 1 rc2 Since the rotation matrix 1 R2 is orthogonal and the matrixA is symmetric (diagonalizable), the left hand side of equation B.20 becomes the diagonalization of the normalized matrix A. Matrix D is the result of this diagonalization. The inverse of square root of the eigenvalues of D are the radii lengths of the ellipsoid. The columns of 1 R2 are the eigenvectors of matrix D. The orientation angles of the ellipsoid can be solved from 1 R2 : ¶ µ1 R2 (3, 2) Ax = arctan2 1 (B.21) R2 (3, 3) ¡ 1 ¢ Ay = arcsin − R2 (3, 1) (B.22) µ1 ¶ R2 (2, 1) Az = arctan2 1 ) (B.23) R2 (1, 1) Note that arctan2 is the four quadrant inverse tangent. B.3 Generating Ellipses from an Ellipsoid For simulation purposes, we needed to generate ellipse sample points from a given ellipsoid. The ellipsoid represent the prostate and the ellipses represent the contours of the prostate as seen in different sagittal US image planes. 124 B.3. Generating Ellipses from an Ellipsoid We approached the problem as a plane (the “slicing plane”) intersecting − → an ellipsoid. We derived the algebraic parameter of the ellipse P ellipse , → − [a b c d e f ]T given the algebraic parameter of the ellipsoid P ellipsoid , [A B C D E F P Q R S]T , and the location and orientation of the “slicing plane” (x, y, z, yaw, pitch, roll). We start by expressing the 2D and 3D conic equations using the homogeneous coordinate convention. Using this convention, a point in 2D and − → −→ 3D can be expressed as X 0 [x y 1]T and X 00 [x y z 1]T . Consequently, the 2D (A.1) and 3D (B.1) conic equation can be written into matrices forms as: x ¤ a b d x y 1 b c e y d e f 1 − →0 T − → X Pellipse X 0 a d e p x £ ¤ d b f q y x y z 1 e f c r z p q r s 1 −→00 T −→00 X Pellipsoid X £ = 0 = 0 (B.24) = 0 = 0 (B.25) − → −→ Notice that we can obtain X 0 by pre-multiplying X 00 with a projection matrix: x x 1 0 0 0 y = 0 1 0 0 y z 1 0 0 0 1 1 − →0 −→00 (B.26) X = Pproj X The next step is to attach a coordinate frame {pl} to the desired slicing plane. The x-axis and y-axis of the frame are parallel to the orthogonal vectors subtending the plane and the z-axis is parallel to the plane’s normal vector. − → We can express points in the reference coordinate frame ref X 0 with respect to the slicing plane’s frame {pl} using the homogeneous transformation 125 B.3. Generating Ellipses from an Ellipsoid matrix ref T pl . →00 ref − X ref = ref = Tpl Tpl →00 pl − X ref t ref R pl 000 x ref t y ref t z (B.27) (B.28) 1 Where ref Rpl is the rotation matrix defined using the euler yaw-pitch-roll convention (equation B.15). The vector [ref tx ref ty ref tz ]T is the location of the center coordinate of {pl} in frame {ref }. We apply the coordinate transformation in equation B.28 to the ellipsoid equation: →00 T ref − ref ( Tpl →00 T pl − X Pellipsoid ref →00 ref − X →00 pl − = 0 X ) Pellipsoid ( Tpl X ) = 0 →T −→ pl − X 00 pl Pellipsoid pl X 00 = 0 (B.29) To obtain the ellipse parameters, we project the transformed 3D points into 2D following equation B.26: −→ −→ (Pproj pl X 00 )T pl Pellipsoid (Pproj pl X 00 ) = 0 →T −→ pl − T pl X 00 (Pproj Pellipsoid Pproj ) pl X 00 = 0 →T − → pl − X 0 Pellipse pl X 0 = 0 The parameter of the desired ellipse would be: Pellipse (1, 1) 2 Pellipse (1, 2) Pellipse (2, 2) − → P ellipse = 2 Pellipse (1, 3) 2 Pellipse (2, 3) Pellipse (3, 3) (B.30) (B.31) − → When a desired plane does not intersect a given ellipsoid, P ellipse will contain imaginary numbers. 126 Appendix C The IMMPDAF Edge Detector One of the components of the transversal segmentation algorithm is a the IMMPDAF edge detector. The detector is based on three concepts: discrete time Kalman filter, Probabilistic Data Association Filter (PDAF) [72], and multiple model interaction. The boundary of the prostate is modeled as the random trajectory of an object travelling around the seed point. The seed point is an imaginary point located at the prostate center. Two discrete time Kalman filter models are used to track the object’s trajectory. The measurement vector of the filters is obtained by combining multiple candidate edge points obtained from the US image’s pixel values. This idea is analogous to tracking a single target under a clutter of measurements (the PDAF). At each iteration step, the trajectory estimate can switch from one Kalman filter mode to the other based on some probability value. The final prostate boundary point is the combination of the estimate from the two modes. This combination process is the multiple model interaction. Note, the word “mode” is used to refer to each Kalman filter model. The algorithm can be divided into three main steps: state estimate initialization, discrete time Kalman filtering, and mode combination. Each step is described in the following subsections. C.1 State Estimate Initialization The algorithm starts by estimating the initial states of the Kalman filters by using previous step estimates X̂m (k − 1|k − 1). The index m represents the mode number (1 or 2). Each X̂m (k − 1|k − 1) is assigned a posterior transition probability µmn (k −1|k −1) (the probability of the filter switching from mode n to m at step k-1). It is calculated as follows: µmn (k − 1|k − 1) = pmn µm (k − 1) P2 i=1 pmi µi (k − 1) (C.1) 127 C.2. Discrete Time Kalman Filter The variable pmn is the prior transition probability that the filter switch from mode n to mode m. It is assumed that the prior probability of the system staying on the same mode is 75%: · ¸ · ¸ p11 p12 0.75 0.25 = (C.2) p21 p22 0.25 0.75 The variable µm (k − 1) is the mode occurence probability at step k-1. Its value is updated at the last step of the IMMPDAF algorithm. When k = 1, µ1 (0) = 1 and µ2 (0) = 0 The initial state estimate for mode m of the Kalman filter is the expected values of the distribution: X̂0m (k − 1|k − 1) = 2 X X̂n (k − 1|k − 1)µmn (k − 1|k − 1) (C.3) n=1 The initial state covariance is calculated as follows: P0m (k − 1|k − 1) = 2 X µmn (k − 1|k − 1) {Pn (k − 1|k − 1) + n=1 [X̂n (k − 1|k − 1) − X̂0m (k − 1|k − 1)] [X̂n (k − 1|k − 1) − X̂0m (k − 1|k − 1)]T } (C.4) Pn (k − 1|k − 1) is the state covariance at step k-1. When k = 1, P01 (k − 1|k − 1) = P02 (k − 1|k − 1) = I (the identity matrix). C.2 Discrete Time Kalman Filter The initial state estimate, X̂0m (k −1|k −1), and covariance, P0m (k −1|k −1), is passed to the Kalman filter iterations to obtain the current state estimate X̂m (k − 1|k − 1). C.2.1 State Space Model The general discrete time state space equation is given as: − → − → − → − → X (k) = A(k − 1) X (k − 1) + B U (k − 1) + V (k − 1) → − − → − → Z (k − 1) = H(k − 1) X (k − 1) + T (k − 1) (C.5) (C.6) − → − → Where X is the system state, A is the state matrix, B is the input matrix, U − → is the input vector, Z is the measurement vector, and H is the measurement 128 C.2. Discrete Time Kalman Filter − → → − matrix. V is the system noise vector with covariance matrix Q and T is the measurement noise vector with covariance matrix R. The prostate boundary is modeled as the trajectory of an object travelling, without extenal influence, at a random distance r(k∆θ) from the “seed point”. The variable k represents sampling of space (not of time). From the seed point, Nr radii to spread outwards, each separated by angle ∆θ = 360 Nr . In the algorithm Nr = 120 and the first radii (k=0) coincides with the x+ axis. As the object traverse from (k−1)∆θ to k∆θ the distance between the object and the seed point is expected to be constant: r(k∆θ) = r((k − 1)∆θ)). The model is linearized around r((k − 1)∆θ)) using a Taylor approximation: r(k∆θ) = r((k − 1)∆θ) + ṙ((k − 1)∆θ) [(k∆θ − (k − 1)∆θ] 1 + r̈((k − 1)∆θ) [k∆θ − (k − 1)∆θ]2 + .... 2 (C.7) For convenience, we drop the notation ∆θ from equation (C.7): r(k) = r(k − 1) + ṙ(k − 1) ∆θ + 1 r̈(k − 1) ∆θ2 + ... (C.8) 2 The first derivative of (C.8) is: ṙ(k) = ṙ(k − 1) + r̈(k − 1) ∆θ + 1 δ 3 r(k − 1) ∆θ2 + ... (C.9) 2 δk 3 The state space model of the object is constructed based on the first and second order derivatives (dropping third order derivatives and higher): · ¸ · ¸· ¸ · 1 2¸ r(k) 1 ∆θ r(k) ∆θ = + 2 r̈(k − 1) (C.10) ṙ(k) 0 1 ṙ(k) ∆θ The acceleration r̈(k − 1) is assumed to be white Gaussian noise with 2 . The index m represents the mode number. Each mode has variance of σm different variance number: σ12 = 10 pixel/rad2 and σ22 = 105 pixel/rad2 . The system noise covariance matrix Q(k) becomes: ¸ · 1 2¸ · 2 £1 2 ¤ σm 0 2 ∆θ Q(k) = 2 ∆θ ∆θ 2 0 σm ∆θ " 4 # 3 Q(k) = ∆ θ 4 ∆3 θ 2 ∆ θ 2 ∆2 θ 2 σm (C.11) 129 C.2. Discrete Time Kalman Filter At each space sample, we measure the distance of the object with respect to the seed point. The measurement is affected by white Gaussian noise with variance λ2 = 20 pixel/rad2 for both modes. Thus, the measurement vector model and measurement noise covariance are: £ ¤− − → → − → Z (k) = 1 0 X (k) + T (k) (C.12) R(k) = 20 (C.13) − → As apparent in subsequent explanations, the actual value of T (k) is not as important as its noise covariance R(k). To summarize, the following equations characterize the system’s dynamics: · ¸ r(k) X̂(k) = (C.14) ṙ(k) · ¸ 1 ∆θ A = (C.15) 0 1 £ ¤ H = 1 0 (C.16) " 4 # 3 Qm (k) = ∆ θ 4 ∆3 θ 2 ∆ θ 2 ∆2 θ 2 σm R(k) = 20 (C.17) (C.18) As a side note, Badiei et. al. compared different system models and found that the model described by equation (C.14) - (C.18) is the most suitable and stable for prostate image segmentation. The values used in the model were optimized through trial and error. C.2.2 Kalman Filter Iteration The iteration steps for the discrete Kalman filter are as follows: 1. Compute the prior estimate for the state at step k and its covariance using the initial state estimate and covariance: X̂m (k|k − 1) = A X̂0m (k − 1|k − 1) T Pm (k|k − 1) = A P0m (k|k − 1) A + Qm (k − 1) (C.19) (C.20) 2. Estimate the measurement at step k and its covariance using results from step 1: Ẑm (k|k − 1) = H X̂m (k|k − 1) (C.21) 0 Sm (k) = H Pm (k|k − 1) H + Rm (k − 1) (C.22) 130 C.2. Discrete Time Kalman Filter 3. Compute the posterior estimate for the state at step k and its posterior covariance using the filter gain W(k) and error vector E(k). −1 Wm (k) = Pm (k|k − 1) H 0 Sm (k|k − 1) (C.23) Em (k) = Z(k) − Ẑm (k|k − 2) (C.24) X̂m (k|k) = X̂m (k|k − 2) + Wm (k) Em (k) (C.25) Pm (k|k) = (I − Wm (k)H)Pm (k|k − 1) T + Wm (k) (y(k) − z 2 (k)) Wm (k) (C.26) The variables y(k) and z(k) will be explained in the section C.2.3. C.2.3 Obtaining the Measurements: the PDAF Concept The Kalman filter requires the tracked object’s measurement value: z(k). In this case, it is the distance of the edge point with respect to the “seed point”. The edge point is obtained by combining several candidate edge points using the PDAF - similar to tracking a single object under a clutter of measurements [72]. To extract the candidate edge points, the following 1 dimensional edge detection filter is applied along one radii: F (r(k), k∆θ) = 1 |{I(r(k) + 2, k∆θ) + I(r(k) + 1, k∆θ) + I(r(k), k∆θ) 3 −I(r(k) − 1, k∆θ) − I(r(k) − 2, k∆θ)}| (C.27) I(r(k),k) is the pixel value at distance r(k) from the seed points and at angle k ∆θ. As mentioned in the previous section, from the seed point a number of radii spreads out and each radii is separated by angle ∆θ. When k = 0, the radii coincides with the x+ axis. The value of r(k) is bounded by the ellipse fit of the user-initialized points. Like the radii, the ellipse is also sampled at k∆θ. The value of r(k) is limited to within 15 pixels within the ellipse sample point. This is how the ellipse fit guides the IMMPDAF edge detector. From the edge detection results, the 5 largest edge points are considered as candidate edge points. A Gaussian probability density function is assigned to the candidate edge points. The probability density function, βi , is proportional to the squared magnitude of the edge and centered around the prior estimate of the measurement Ẑ(k|k − 1): 131 C.3. Combining the Estimates From the Two Modes pi (k) = βi = Fi (ri (k), k∆θ) (ri (k) − Ẑ(k|k − 1))2 p exp− 2 Sm (k) 2πSm (k) pi P5 i=1 pi (C.28) (C.29) The true measurement z(k) is the expected value of the distribution: 5 X z(k) = ri (k)βi (k) (C.30) i=1 And y(k) from previous section is calculated as follows: y(k) = 5 X ri2 (k)βi (k) (C.31) i=1 C.3 Combining the Estimates From the Two Modes The last step of the IMMPDAF algorithm is combining the estimates from the two modes. Each mode is assigned a Gaussian probability density function centered around the measurement estimate and proportional to the their mode occurrence probability µm (k − 1) γm (k) = 1 1 −1 p (k) Em (k)} (C.32) exp{− Em (k)T Sm 2 2π|Sm (k)| µm (k) = P2 γm (k)µm (k − 1) n=1 γn (k)µm n(k − 1) (C.33) The final estimate for the object’s state is the expected value of the distribution: X̂(k|k) = 2 X X̂m (k|k) µm (k) (C.34) m=1 The boundary point location (segmentation result) the final measurement value Z(k|k): Ẑ(k + 1|k + 1) = H X̂(k|k) (C.35) 132 C.3. Combining the Estimates From the Two Modes X̂m (k|k) and µm (k) transfer over to the next step of the iteration as X̂m (k − 1|k − 1) and µm (k − 1). One can also calculate the final state and measurement covariance (P(k—k) and S(k—k)) but they are not needed for the next step of the iteration. 133 Appendix D Probe Calibration Using Stradwin and the Single Wall Technique We needed to determine the position and orientation of the US coordinate frame {U S} which was defined in Section 4.1 with respect to the position sensors (Optotrak markers attached on the probe), i.e. to calibrate the US probe. Calibration was performed using the single wall technique [64] implemented in Stradwin. D.1 The Framework of Single Wall Calibration Technique As illustrated in Figure D.1, there are 4 coordinate frames: the Optotrak camera’s coordinate frame {Opto}, the probe’s coordinate frame {P robe}, the 3D US coordinate frame {U S} and the US image coordinate frame {IM }. There is a 9.87 mm offset between {IM } and {U S} due to the TRUS probe’s radius as well as orientation difference. Equation D.1 expresses {IM } relative to {U S}. US TIM 0 0 = 1 0 0 −1 0 1 0 9.87 0 0 0 0 0 1 (D.1) US image calibration using the single wall technique gives the homogeneous transformation from {IM } to {P robe}: P robe TIM . The results are − → the eight values for the calibration parameter vector P [x y z Rx Ry Rz Sx Sy ] T . The parameters (x,y,z) - in mm - describe the location of the origin of {IM } with respect to the origin of {P robe}. The parameters (Rx , Ry , Rz ) - in degrees - describe the orientation of the axes of {IM } with 134 D.1. The Framework of Single Wall Calibration Technique respect to the axes of {P robe}. The parameter Sx and Sy are the scaling factor (mm/pixel). To obtain P robe TU S we post multiply P robe TIM with the inverse of U S TIM . Figure D.1: The framework of TRUS probe to Optotrak calibration. The single wall technique assigns an arbitrary coordinate frame {C} is attached to the floor of the water bath giving the following composition of homogoneous coordinate transformations: C IM x Sx x IM y C y C S = TOpto Opto TP robe P robe TIM y (D.2) 0 0 1 1 In equation (D.2), Opto T P robe and the pixel values (IM x, IM y) are the only 135 D.2. US Image and Optotrak Data Collection → − known values to solve for the calibration parameter P . Thus, only the zero component of (D.2) gives one equation in 11 unknowns (6 for P robe TImage and five for C TOpto ). In the US image, the water bath’s floor can be detected as a line. Since two pixel locations uniquely define the line, two equations can be extracted from each image frame. For m image frames, we have 2m equations: an over-determined system of equations. To solve for the unknowns, the single wall technique uses the Marquardt-Levenberg optimization algorithm. The − → algorithm gives the best values for P and calculated root mean squared error in mm. D.2 US Image and Optotrak Data Collection Six Optotrak markers were attached to the TRUS probe. The markers were glued onto the metallic casing that enclosed the probe. The casing were used to prevent direct damage to the probe due to the attachment of the markers. The markers are labeled from 1 - 6 and are connected to the Optotrak data acquisition unit via its strobers (see Figure D.2). We used the 6 markers configuration to define the probe as a rigid body using the NDI RigMaker software (NDI, ON, Canada). This software only needed a minimum of three markers to define an object as a rigid body. We used three extra markers in case some of them were out of the Optotrak camera’s field of view when the TRUS probe was translated or rotated. By defining the probe as a rigid body, a coordinate frame {P robe} is attached to the probe. The location of marker 1 defines the origin of this coordinate frame. The TRUS Probe was mounted on the EXII Stepper and Micro-Touch R mount and was partially submerged in a water bath. Due to the LP° stepper and mount’s workspace limitations, the probe was submerged in an upside down position. Hence, the wall of the water bath, not the floor, was used as the reflector for the single wall calibration technique. The temperature of the water bath was kept at 40o Celsius by frequent mixing with heated water. A thermometer was used to measure the temperature of the water bath at all time (see Figure D.3). The NDI Toolbench software (NDI, ON, Canada) collected the Optotrak readings. The data consisted of the location of each of 6 the markers (xi=1..6 , yi=1..6 , zi=1..6 ) and the position-orientation (x, y, z, Rx , Ry , Rz ) of {P robe} with respect to {Opto}. Position (x, y, z) were expressed in mm with respect to the origin of {Opto}. Orientation data (Rx , Ry , Rz ) were expressed in 136 D.2. US Image and Optotrak Data Collection Figure D.2: Calibration setup: US probe and Optotrak configuration. Top left: Optotrak markers attached on the probe. Top right: Optotrak strobers. Bottom: Optorak data acquisition unit and Optotrak camera. 137 D.2. US Image and Optotrak Data Collection Figure D.3: The water bath calibration setup. Left: the TRUS probe inside the water bath in upside down position. Right: the US machine, EXII stepper, and Micro-Touch mount. 138 D.2. US Image and Optotrak Data Collection degrees using the yaw-pitch-roll Euler angles convention. We used the Ultrasonix Ulterius research package to capture the post processed sagittal B-mode images. The US probe was set at 6.6 MHz center frequency with imaging depth of 60 mm. These parameters gave X-axis and Y-axis conversion ratio of 156 microns/pixel. The image width and height were 350 pixels and 384 pixels respectively. We collected 50 images of the water bath’s wall at different TRUS probe poses and positions. The probe was manipulated using the stepper and the mount to cover all the needed 6 degrees of freedom illustrated in Figure D.4. At each pose and location, the B-mode US image and Optotrak readings were captured simultaneously and independently as Stradwin doesn’t support the real time data capture using the Optotrak. Figure D.4: Types of poses and motions required for single wall calibration method. The Optotrak data and the post processed B-mode images were combined into Stradwin standard files format: a pair of files with “.sw” and “.sxi” extensions. The first is a text file containing information such as: the US imaging settings, position data from the location sensors, the calibration parameters, etc. The later is a binary image data file. A Matlab script was used to combine the B-mode images and Optotrak data into Stradwin standard files. First, header parameter tokens are written into the “.sw” file. The tokens’ name are self explanatory and the parameter values are taken from the actual US imaging parameters. In our calibration we used 51 image frames although we collected 50 images of the water bath. The first image, 139 D.3. Single Wall Calibration Technique Using Stradwin a phantom image, was an extra. It was only used to define the boundary of the image (will be explained later) but it is not used for calibration. The following are the parameter tokens written into the “.sw” file: RES RES RES RES RES RES RES RES RES BUF FRAMES 51 BUF WIDTH 350 BUF HEIGHT 384 POS REC true BUF RF false RF VECTORS 127 (irrelevant, used default value) RF SAMPLES 3128 (irrelevant, used default value) END HEADER BIN IM FILENAME Calibration.sxi Following the header parameter tokens, the Optotrak position data of the 50 images were written into the “.sw” file: IM 0 0 0 0 0 0 0 IM 1 -42.2469 0.40293 -189.9733 -4.324 -17.4499 -15.8196 IM 2 -42.2362 -1.356 -190.1175 -4.3067 -17.4411 -15.8755 .. . The phrase “IM” indicates that the information is tagged to each images inside the “.sxi” file. It is followed by the frame number and the Optotrak position reading at the instant image was captured in (x, y, z, Rz , Ry , Rx ) format. Note that Stradwin expects the image location in cm (instead of mm) and that there is a change in Rx and Rz order compared to the Optotrak data format. Since frame 0 (the first frame) is not used for calibration, its position and orientation is set to be zeros. D.3 Single Wall Calibration Technique Using Stradwin We loaded the “.sw” and “.sxi” files into the Stradwin research software and selected the probe calibration tab. First, we set the temperature to 40o C. The next step is to set the image scale (horizontal and vertical) using the image at frame 0. This image is an image of a PVC phantom with clear boundaries. It was recorded using the exact same imaging parameters as the calibration’s. The phantom image clearly showed the image frame 140 D.3. Single Wall Calibration Technique Using Stradwin top-bottom and left-right boundaries. Hence, we can set the horizontal and vertical scale more accurately. We set the horizontal image scale by clicking the “Mark ends of first line” button. Then we left clicked on the bottom left corner of the image (a red ‘x’ appeared) and right clicked on the bottom right corner of the image (a green ‘x’ appeared). Then, we clicked the “Finished locating ends of line” button and set the first line length to be 54.6 mm (350 pixels times 156 microns per pixels). Next, we set the vertical image scale by clicking the “Mark ends of second line” button. This step was followed by left clicking on the top right corner of the image and right clicking on the bottom right corner of the image. Red and green ‘x’s appeared once again. We finished by setting the length of the line to 60 mm (the imaging depth). Figure D.5: Setting the calibration image scales: marking the line (red and green x’s in the blue circles) and setting the line length (red circle) 141 D.4. Automatic Line Detection and Optimization D.4 Automatic Line Detection and Optimization The water bath’s wall shows up in the image as a thick bright line. The thickness of the line depends on the pose of the TRUS probe relative to the wall. Sometimes the image shows two bright lines due to reverberation. Stradwin applies an automatic line detection algorithm to find the wall of the water bath. The algorithm is a combination of Laplacian of Gaussian (LoG) edge detector and the Random Sample Consensus (RANSAC) algorithm. Figure D.6: Automatic line detection in Stradwin. The wall of the water bath appears as a thick bright line in the US image. Vertical red lines divides the image into bands. The line detection algorithm can be adjusted by the parameters in the blue circle The image is divided into several bands by several vertical red lines. The number of bands is adjustable through “Vertical analysis bands” parameter. 142 D.4. Automatic Line Detection and Optimization At each vertical line, the image pixel values are collected into a vector and a one dimensional LoG operator is applied to detect the strongest edge point. The user can change the variance of the Gaussian operator and the threshold of the gradient by adjusting the “Variance of Gaussian Kernel” and “Gradient threshold” parameters. Stradwin applies the RANSAC algorithm to the edge points to detect the best line fit that represent the wall. RANSAC selects all permutations of two edge points and fits lines to each of them. Each line is ranked based on its consensus score. An edge point contributes to the line consensus score if it is located within a certain distance from the line. This distance is adjustable through the “Pixel linearity threshold” parameter. The line with the highest consensus score is the best candidate for the wall. The consensus score is checked against the maximum possible consensus score (the number of vertical red lines). This number is not the same as the number of bands in the image. If the consensus score pass a certain fraction of the maximum possible score, the line will be used for calibration. This fractional number can be adjusted through the “Ransac accept threshold” parameter. For our calibration we used the default parameter values: 1. Variance of Gaussian kernel: 5.0 2. Vertical analysis band: 20 3. Gradient threshold 5.0 4. Pixel linearity threshold :3 5. Ransac accept threshold: 0.6 Fifty lines from fifty image frames were used to determine the calibration parameters. 143 Appendix E Registration Experiment Results E.1 Registration Experiment Results: the Sparse Volume Registration Table E.1 shows the registration experiment results using contours which were obtained by segmenting several sagittal US images from each volume data using the sagittal segmentation algorithm. In the table, we report relative transformation between adjacent poses. For example, the results in column “Vol 2” explain the transformation from pose 1 to pose 2, the results in column “Vol 3” explain the transformation from pose 2 to pose 3, and so on. The results in column “Vol 3-1” explain the relative transformation from pose 3 to pose 1 since the relative transformation from pose 1 to pose 1 is always zero. The transformation is expressed in three axes translation(Tx , Ty , Tz ) and three axes rotation (Rx , Ry , Rz ). The bounding box errors were calculated by defining two smallest rectangular boxes that cover the intra-operative and pre-operative contour points. The boxes’ vertices are the permutations of the X-axis, Y-axis, and Z-axis extrema of the intra-operative and pre-operative contour points. The later is registered using the registration matrix obtained from the ellipsoid fitting based registration algorithm. Then, the bounding box error is calculated by taking the root mean square distance between corresponding vertices. The overall average bounding box error is 3.34 mm while the overall average Hausdorff distance is 11.41 mm. The standard deviations are 1.71 mm and 3.89 mm respectively. We report the dial numbers on the motion stage, the transformation parameters calculated from the Optotrak measurements, and the transformation parameters calculated using our ellipsoid fitting registration algorithm. The registration error is calculated by taking the absolute difference between the later two. 144 Source Motion Platform Optotrak Ellipsoid Fitting Reg. Error |Ellipsoid Fit Reg. - No Motion Vol 3-1 Vol 2 Vol 3 Vol 4 X axis translation Vol 5 Vol 6 Vol 7 Vol 8 Y axis translation Z axis translation Vol 9 Vol 10 Vol 11 Vol 12 Vol 13 Vol 14 Vol 15 0.00 0.00 0.00 5.00 -5.00 5.00 -5.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 Ty (mm) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 5.00 -5.00 5.00 -5.00 0.00 0.00 0.00 0.00 0.00 Tz (mm) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 5.00 -5.00 5.00 -5.00 Rx ( o ) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 Ry ( o ) Tx (mm) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 -0.01 3.06 -3.06 2.98 -2.97 0.14 -0.14 0.13 -0.14 0.51 -0.50 0.51 -0.51 Ty (mm) 0.00 -0.01 0.01 0.41 -0.38 0.37 -0.33 2.86 -2.96 2.93 -2.95 -0.10 0.13 -0.12 0.15 Tz (mm) 0.00 0.00 0.00 -0.26 0.25 -0.26 0.26 0.03 -0.04 0.02 -0.03 5.01 -5.01 5.01 -5.01 Rx ( o ) 0.00 0.00 0.00 0.02 -0.01 0.01 0.00 -0.36 0.35 -0.36 0.36 -0.02 0.02 -0.02 0.03 Ry ( o ) 0.00 0.00 0.00 0.24 -0.25 0.25 -0.25 0.07 -0.06 0.05 -0.06 -0.04 0.03 -0.03 0.03 Rz ( o ) Tx (mm) 0.01 0.00 -0.02 0.20 -0.20 0.18 -0.20 0.19 -0.18 0.17 -0.18 -0.03 0.03 -0.02 0.01 -2.08 0.64 1.41 10.43 -12.02 14.71 -9.12 -7.21 4.80 0.13 -6.25 17.53 -17.94 -2.66 8.56 Ty (mm) 0.06 0.48 -0.53 0.40 -0.50 2.94 -3.73 3.54 -1.70 -0.61 -0.69 Tz (mm) -0.20 -0.69 1.07 0.13 2.66 1.24 1.87 0.53 2.01 -0.43 -1.84 11.56 2.45 1.95 -2.11 Rx ( o ) 0.14 0.87 -1.03 0.73 -0.40 0.74 0.10 -0.61 -0.89 -0.13 2.59 -1.32 -1.68 3.28 -0.03 Ry ( o ) 5.07 -3.94 -1.11 -13.63 16.79 -20.66 14.43 11.90 -10.92 2.41 6.21 -35.48 33.51 9.91 -22.03 Rz ( o ) Tx (mm) 0.05 -1.24 1.13 -0.72 -0.36 0.46 1.70 -1.15 -0.38 0.78 -2.30 0.53 0.86 2.81 -4.11 2.08 0.64 1.42 7.37 8.96 11.74 6.15 7.36 4.94 0.00 6.12 17.02 17.43 3.17 9.07 Ty (mm) 0.06 0.49 0.55 0.01 0.35 0.01 0.17 0.08 0.77 0.61 1.25 0.51 1.42 1.64 0.84 -0.03 0.38 -1.29 1.52 Tz (mm) 0.19 0.70 1.07 0.38 2.41 1.50 1.61 0.50 2.04 0.44 1.80 6.56 7.46 3.06 2.90 Rx ( o ) 0.14 0.88 1.03 0.71 0.39 0.73 0.11 0.25 1.24 0.23 2.24 1.31 1.70 3.31 0.05 Ry ( o ) 5.07 3.94 1.11 13.87 17.04 20.91 14.68 11.83 10.87 2.35 6.27 35.44 33.48 9.95 22.06 Rz ( o ) Optotrak| Bounding Box Error (mm) 0.04 1.24 1.15 0.92 0.16 0.28 1.90 1.35 0.20 0.60 2.12 0.56 0.83 2.84 4.13 1.32 1.08 0.61 2.87 2.85 4.41 2.51 2.37 2.35 2.13 2.35 7.52 4.29 4.73 2.04 Hausdorff Distance (mm) 4.37 4.57 4.01 12.50 12.05 14.01 11.40 10.55 11.06 7.43 9.78 18.16 18.09 10.12 15.17 145 E.1. Registration Experiment Results: the Sparse Volume Registration Table E.1: Registration experiment results: the sparse volume registration. Motion Type Transformation Tx (mm) Source Motion Platform Optotrak Ellipsoid Fitting Reg. Error |Ellipsoid Fit Reg. - X axis rotation (yaw) Y axis translation (pitch) Combination Vol 16 Vol 17 Vol 18 Vol 19 Vol 20 Vol 21 Vol 22 Vol 23 Vol 24 Vol 25 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 3.00 -3.00 Ty (mm) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 3.00 -3.00 Tz (mm) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 3.00 -3.00 Rx ( o ) 5.00 -5.00 5.00 -5.00 0.00 0.00 0.00 0.00 -3.00 3.00 Ry ( o ) Tx (mm) 0.00 0.00 0.00 0.00 5.00 -5.00 5.00 -5.00 3.00 -3.00 0.19 -0.22 0.19 -0.19 0.09 0.31 0.12 0.29 0.85 -0.55 Ty (mm) -0.20 0.38 -0.31 0.43 -0.24 0.25 -0.19 0.24 0.65 -0.28 Tz (mm) 1.63 -1.64 1.65 -1.63 4.67 -4.79 4.78 -4.74 6.69 -6.76 Rx ( o ) -4.94 4.99 -4.91 4.94 -0.35 0.36 -0.34 0.36 -3.37 3.39 Ry ( o ) -0.29 0.32 -0.27 0.30 4.59 -4.79 4.77 -4.74 2.87 -2.83 Rz ( o ) Tx (mm) 0.38 -0.39 0.38 -0.35 0.22 -0.22 0.18 -0.21 0.80 -0.93 -1.98 4.40 -1.13 -7.11 -0.89 10.41 -12.78 13.41 -10.84 0.78 Ty (mm) -0.45 -0.42 0.15 -0.07 -1.03 1.83 1.83 -1.95 0.13 1.36 Tz (mm) -0.32 2.36 -1.64 1.96 6.47 -5.60 4.75 2.45 8.94 -7.54 Rx ( o ) -4.40 3.22 -3.15 3.43 -1.25 1.39 1.76 -5.30 -5.63 4.99 Ry ( o ) 9.62 -11.03 1.13 12.60 6.56 -16.09 19.68 -24.24 13.78 0.77 Rz ( o ) Tx (mm) 2.20 -0.95 -0.66 -0.88 0.45 2.85 -5.09 3.58 -6.84 3.48 2.17 4.61 1.32 6.92 0.99 10.10 12.90 13.12 11.69 1.33 Ty (mm) 0.24 0.81 0.46 0.50 0.79 1.58 2.19 0.52 1.64 2.02 Tz (mm) 1.95 4.00 3.28 3.58 1.80 0.81 0.02 7.19 2.26 0.78 Rx ( o ) 0.54 1.76 1.76 1.51 0.90 1.03 2.10 5.66 2.25 1.60 Ry ( o ) 9.91 11.36 1.40 12.30 1.97 11.29 14.91 19.50 10.90 3.61 1.82 0.56 1.04 0.52 0.23 3.06 5.26 3.79 7.64 4.41 4.31 4.29 4.04 4.25 4.52 1.83 6.42 1.98 5.84 2.68 11.66 12.19 11.09 11.27 11.42 15.83 16.00 15.49 9.97 7.05 Rz ( o ) Optotrak| Bounding Box Error (mm) Hausdorff Distance (mm) 146 E.1. Registration Experiment Results: the Sparse Volume Registration Table E.1 continued from previous page Motion Type Transformation Tx (mm) E.2. Registration Experiment Results: the Full Volume Registration E.2 Registration Experiment Results: the Full Volume Registration Table E.2 shows the registration experiment results using contours which were obtained by segmenting the US volume data using the extended sagittal segmentation algorithm. In the table, we report relative transformation between adjacent poses. For example, the results in column “Vol 2” explain the transformation from pose 1 to pose 2, the results in column “Vol 3” explain the transformation from pose 2 to pose 3, and so on. The results in column “Vol 3-1” explain the relative transformation from pose 3 to pose 1 since the relative transformation from pose 1 to pose 1 is always zero. The transformation is expressed in three axes translation (Tx , Ty , Tz ) and three axes rotation (Rx , Ry , Rz ). The bounding box errors were calculated by defining two smallest rectangular boxes that cover the intra-operative and pre-operative contour points. The boxes’ vertices are the permutations of the X-axis, Y-axis, and Z-axis extrema of the intra-operative and pre-operative contour points. The later is registered using the registration matrix obtained from the ellipsoid fitting based registration algorithm. Then, the bounding box error is calculated by taking the root mean square distance between corresponding vertices. The overall average bounding box error is 2.40 mm while the overall average Hausdorff distance is 6.68 mm. The standard deviations are 1.53 mm and 3.49 mm respectively. We report the dial numbers on the motion stage, the transformation parameters calculated from the Optotrak measurements, and the transformation parameters calculated using our ellipsoid fitting registration algorithm. The registration error is calculated by taking the absolute difference between the later two. 147 Source Motion Platform Optotrak Ellipsoid Fitting Reg. Error |Ellipsoid Fit Reg. - No Motion X axis translation Y axis translation Z axis translation Vol 3-1 Vol 2 Vol 3 Vol 4 Vol 5 Vol 6 Vol 7 Vol 8 Vol 9 Vol 10 Vol 11 Vol 12 Vol 13 Vol 14 Vol 15 0.00 0.00 0.00 5.00 -5.00 5.00 -5.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 Ty (mm) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 5.00 -5.00 5.00 -5.00 0.00 0.00 0.00 0.00 0.00 Tz (mm) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 5.00 -5.00 5.00 -5.00 Rx ( o ) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 Ry ( o ) Tx (mm) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 -0.01 3.06 -3.06 2.98 -2.97 0.14 -0.14 0.13 -0.14 0.51 -0.50 0.51 -0.51 Ty (mm) 0.00 -0.01 0.01 0.41 -0.38 0.37 -0.33 2.86 -2.96 2.93 -2.95 -0.10 0.13 -0.12 0.15 Tz (mm) 0.00 0.00 0.00 -0.26 0.25 -0.26 0.26 0.03 -0.04 0.02 -0.03 5.01 -5.01 5.01 -5.01 Rx ( o ) 0.00 0.00 0.00 0.02 -0.01 0.01 0.00 -0.36 0.35 -0.36 0.36 -0.02 0.02 -0.02 0.03 Ry ( o ) 0.00 0.00 0.00 0.24 -0.25 0.25 -0.25 0.07 -0.06 0.05 -0.06 -0.04 0.03 -0.03 0.03 Rz ( o ) Tx (mm) 0.01 0.00 -0.02 0.20 -0.20 0.18 -0.20 0.19 -0.18 0.17 -0.18 -0.03 0.03 -0.02 0.01 0.24 2.32 -2.61 6.01 -5.00 5.11 -5.28 1.61 -0.83 -0.36 -1.25 2.54 -0.29 0.57 -0.79 Ty (mm) 0.48 -0.27 -0.17 1.45 -0.73 -0.81 1.28 1.81 -3.31 4.04 -3.13 0.66 -1.12 -0.96 0.97 Tz (mm) -0.50 0.55 0.09 -2.27 1.99 1.66 -2.39 2.89 0.30 -2.23 0.44 2.67 -1.55 5.99 -6.23 Rx ( o ) 0.64 -0.35 -0.22 2.25 -1.34 -1.68 2.81 -3.27 0.36 1.42 0.37 1.96 -2.76 -1.28 1.55 Ry ( o ) 1.25 -4.79 3.53 -2.98 3.76 -4.32 5.26 -4.88 0.00 2.44 0.85 -5.49 3.53 -1.68 -0.85 Rz ( o ) Tx (mm) 1.27 -0.65 -0.63 1.14 0.55 -0.62 1.10 -1.04 -0.55 0.97 -0.67 -0.66 2.35 -1.19 -1.25 0.24 2.32 2.61 2.95 1.93 2.13 2.31 1.46 0.68 0.49 1.12 2.03 0.21 0.06 0.28 Ty (mm) 0.48 0.26 0.19 1.05 0.35 1.18 1.61 1.06 0.35 1.11 0.18 0.76 1.25 0.84 0.82 Tz (mm) 0.49 0.55 0.09 2.01 1.74 1.92 2.65 2.86 0.33 2.24 0.47 2.33 3.46 0.98 1.22 Rx ( o ) 0.64 0.35 0.22 2.23 1.33 1.69 2.82 2.91 0.01 1.78 0.02 1.97 2.78 1.26 1.52 Ry ( o ) 1.26 4.79 3.52 3.22 4.01 4.57 5.51 4.95 0.05 2.39 0.91 5.45 3.50 1.65 0.88 Optotrak| Rz ( o ) Bounding Box Error (mm) 1.26 0.65 0.61 0.94 0.74 0.80 1.30 1.23 0.37 0.80 0.49 0.63 2.33 1.17 1.27 0.69 1.01 0.57 0.97 1.09 0.96 1.12 1.25 1.04 0.93 0.96 3.81 1.90 3.87 3.36 Hausdorff Distance (mm) 2.75 2.56 2.97 3.62 4.22 4.01 4.01 3.70 2.92 3.61 3.40 7.93 6.83 7.27 7.33 148 E.2. Registration Experiment Results: the Full Volume Registration Table E.2: Registration experiment results: the full volume registration. Motion Type Transformation Tx (mm) Source Motion Platform Optotrak Ellipsoid Fitting Reg. Error |Ellipsoid Fit Reg. - X axis rotation (yaw) Y axis translation (pitch) Combination Vol 16 Vol 17 Vol 18 Vol 19 Vol 20 Vol 21 Vol 22 Vol 23 Vol 24 Vol 25 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 3.00 -3.00 Ty (mm) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 3.00 -3.00 Tz (mm) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 3.00 -3.00 Rx ( o ) 5.00 -5.00 5.00 -5.00 0.00 0.00 0.00 0.00 -3.00 3.00 Ry ( o ) Tx (mm) 0.00 0.00 0.00 0.00 5.00 -5.00 5.00 -5.00 3.00 -3.00 0.19 -0.22 0.19 -0.19 0.09 0.31 0.12 0.29 0.85 -0.55 Ty (mm) -0.20 0.38 -0.31 0.43 -0.24 0.25 -0.19 0.24 0.65 -0.28 Tz (mm) 1.63 -1.64 1.65 -1.63 4.67 -4.79 4.78 -4.74 6.69 -6.76 Rx ( o ) -4.94 4.99 -4.91 4.94 -0.35 0.36 -0.34 0.36 -3.37 3.39 Ry ( o ) -0.29 0.32 -0.27 0.30 4.59 -4.79 4.77 -4.74 2.87 -2.83 Rz ( o ) Tx (mm) 0.38 -0.39 0.38 -0.35 0.22 -0.22 0.18 -0.21 0.80 -0.93 0.20 -1.20 0.75 -1.10 -1.12 2.27 -1.31 2.50 -2.71 3.36 Ty (mm) -0.06 -0.78 0.59 -0.21 0.69 -0.43 -0.59 0.65 0.40 0.96 -6.91 Tz (mm) -1.90 3.37 -3.38 2.19 2.59 -3.20 5.84 -5.51 5.65 Rx ( o ) -3.57 2.26 -2.25 3.27 2.04 -1.52 -1.32 1.14 -2.89 4.32 Ry ( o ) 4.30 -1.82 0.37 0.46 3.15 -3.49 3.89 -5.76 4.93 -4.34 Rz ( o ) Tx (mm) 2.05 -1.97 0.77 -0.77 -1.17 2.12 -1.31 0.72 -2.07 2.29 0.01 0.99 0.56 0.91 1.21 1.96 1.43 2.21 3.56 3.91 Ty (mm) 0.14 1.17 0.90 0.64 0.93 0.68 0.39 0.41 0.25 1.24 Tz (mm) 3.53 5.01 5.02 3.81 2.08 1.59 1.06 0.76 1.04 0.15 Rx ( o ) 1.37 2.72 2.66 1.67 2.39 1.88 0.98 0.78 0.48 0.93 Ry ( o ) 4.58 2.14 0.64 0.15 1.44 1.31 0.88 1.01 2.06 1.51 149 Optotrak| Rz ( o ) Bounding Box Error (mm) 1.67 1.58 0.39 0.41 1.40 2.34 1.49 0.92 2.87 3.23 4.29 5.01 4.96 4.65 3.65 2.39 4.15 2.19 3.46 1.76 Hausdorff Distance (mm) 12.83 13.12 13.03 12.73 7.24 7.95 8.82 9.37 7.32 7.40 E.2. Registration Experiment Results: the Full Volume Registration Table E.2 continued from previous page Motion Type Transformation Tx (mm) E.3. ICP Registration Results E.3 ICP Registration Results Table E.3 shows the registration experiment results using ICP on contours which were obtained by segmenting the US volume data using the extended sagittal segmentation algorithm. In the table, we report relative transformation between adjacent poses. For example, the results in column “Vol 2” explain the transformation from pose 1 to pose 2, the results in column “Vol 3” explain the transformation from pose 2 to pose 3, and so on. The results in column “Vol 3-1” explain the relative transformation from pose 3 to pose 1 since the relative transformation from pose 1 to pose 1 is always zero. The transformation is expressed in three axes translation (Tx , Ty , Tz ) and three axes rotation (Rx , Ry , Rz ). The bounding box errors were calculated by defining two smallest rectangular boxes that cover the intra-operative and pre-operative contour points. The boxes’ vertices are the permutations of the X-axis, Y-axis, and Z-axis extrema of the intra-operative and pre-operative contour points. The later is registered using the registration matrix obtained from the ellipsoid fitting based registration algorithm. Then, the bounding box error is calculated by taking the root mean square distance between corresponding vertices. The overall average bounding box error is 2.06 mm while the overall average Hausdorff distance is 6.20 mm. The standard deviations are 1.32 mm and 3.18 mm respectively. We also report the final root mean square closest points distance error “Final RMS Dist” as well as the ICP registration time. Overall, it takes 76.49 ± 28.84 seconds to register the prostate volumes even after under sampling the data points. The algorithm was implemented in Matlab running on a 1.66 GHz computer with 1 GB RAM. The final RMS distance error is 0.83 ± 0.39 mm. We report the dial numbers on the motion stage, the transformation parameters calculated from the Optotrak measurements, and the transformation parameters calculated using our ellipsoid fitting registration algorithm. The registration error is calculated by taking the absolute difference between the later two. 150 Table E.3: Registration experiment results - ICP. Source Motion Platform ICP Reg. Error |ICP Reg. - No Motion X axis translation Y axis translation Z axis translation Vol 3-1 Vol 2 Vol 3 Vol 4 Vol 5 Vol 6 Vol 7 Vol 8 Vol 9 Vol 10 Vol 11 Vol 12 Vol 13 Vol 14 Vol 15 0.00 0.00 0.00 5.00 -5.00 5.00 -5.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 Ty (mm) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 5.00 -5.00 5.00 -5.00 0.00 0.00 0.00 0.00 Tz (mm) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 5.00 -5.00 5.00 -5.00 Rx (o) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 Ry (o) Tx (mm) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 -0.01 3.06 -3.06 2.98 -2.97 0.14 -0.14 0.13 -0.14 0.51 -0.50 0.51 -0.51 Ty (mm) 0.00 -0.01 0.01 0.41 -0.38 0.37 -0.33 2.86 -2.96 2.93 -2.95 -0.10 0.13 -0.12 0.15 Tz (mm) 0.00 0.00 0.00 -0.26 0.25 -0.26 0.26 0.03 -0.04 0.02 -0.03 5.01 -5.01 5.01 -5.01 Rx (o) 0.00 0.00 0.00 0.02 -0.01 0.01 0.00 -0.36 0.35 -0.36 0.36 -0.02 0.02 -0.02 0.03 Ry (o) 0.00 0.00 0.00 0.24 -0.25 0.25 -0.25 0.07 -0.06 0.05 -0.06 -0.04 0.03 -0.03 0.03 Rz (o) Tx (mm) 0.01 0.00 -0.02 0.20 -0.20 0.18 -0.20 0.19 -0.18 0.17 -0.18 -0.03 0.03 -0.02 0.01 -0.97 1.14 -1.43 3.84 -3.13 4.19 -3.56 1.29 -0.50 -0.06 -0.43 0.76 -0.40 -0.65 -0.23 Ty (mm) 0.12 -0.31 -0.32 0.82 -0.05 -0.78 0.68 2.09 -3.34 3.52 -3.03 0.29 -0.97 -0.64 -0.80 Tz (mm) -0.19 0.58 -1.36 0.55 1.45 -1.35 2.21 0.23 -0.92 -0.09 4.09 -2.12 5.70 -3.15 0.33 Rx (o) 0.25 -0.59 -0.65 1.26 0.05 -1.74 1.66 -2.76 0.36 0.18 0.84 0.82 -2.07 -0.66 -1.34 Ry (o) 1.12 -1.37 1.02 -0.77 0.58 -1.34 2.45 -1.89 -0.10 0.51 0.50 -1.65 1.05 -0.61 -0.05 Rz (o) ICP Time (s) -0.39 0.11 -0.77 -0.04 0.93 0.12 1.25 0.19 -0.17 -0.01 0.07 -0.24 0.51 -2.02 0.18 110.7 89.9 83.8 59.1 119.2 132.9 83.2 24.7 40.9 42.5 64.8 48.8 115.6 32.1 Final RMS Dist (mm) Tx (mm) 0.58 0.48 0.47 0.55 0.57 0.56 0.61 0.55 0.58 0.54 0.57 0.71 0.99 0.60 1.17 0.97 1.14 1.42 0.78 0.07 1.21 0.58 1.15 0.36 0.19 0.29 0.25 0.10 1.16 0.28 94.8 Ty (mm) 0.12 0.30 0.33 0.42 0.33 1.15 1.02 0.77 0.39 0.59 0.08 0.39 1.10 0.52 0.94 Tz (mm) 0.18 0.58 0.33 1.11 0.30 1.71 1.60 2.18 0.26 0.94 0.06 0.91 2.88 0.69 1.86 Rx (o) 0.25 0.59 0.65 1.25 0.06 1.75 1.66 2.40 0.01 0.54 0.49 0.84 2.08 0.64 1.37 Ry (o) 1.12 1.37 1.02 1.02 0.83 1.59 2.70 1.96 0.05 0.46 0.56 1.61 1.02 0.58 0.08 0.40 0.11 0.75 0.24 1.13 0.06 1.45 0.00 0.01 0.19 0.26 0.21 0.49 2.00 0.16 0.76 2.65 0.61 2.87 0.27 2.54 0.79 4.21 0.62 4.33 0.54 4.03 0.72 3.93 0.83 3.31 1.03 2.94 0.80 3.45 0.88 3.71 3.63 6.61 2.32 5.80 3.81 7.35 2.83 7.37 Rz (o) Optotrak| Bounding Box Error (mm) Hausdorff Distance (mm) E.3. ICP Registration Results Optotrak Motion Type Transformation Tx (mm) 151 Table E.3 continued from previous page Source Motion Platform ICP Reg. Error |ICP Reg. - X axis rotation (yaw) Y axis translation (pitch) Combination Vol 16 Vol 17 Vol 18 Vol 19 Vol 20 Vol 21 Vol 22 Vol 23 Vol 24 Vol 25 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 3.00 -3.00 Ty (mm) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 3.00 -3.00 Tz (mm) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 3.00 -3.00 Rx (o) 5.00 -5.00 5.00 -5.00 0.00 0.00 0.00 0.00 -3.00 3.00 Ry (o) Tx (mm) 0.00 0.00 0.00 0.00 5.00 -5.00 5.00 -5.00 3.00 -3.00 0.19 -0.22 0.19 -0.19 0.09 0.31 0.12 0.29 0.85 -0.55 Ty (mm) -0.20 0.38 -0.31 0.43 -0.24 0.25 -0.19 0.24 0.65 -0.28 Tz (mm) 1.63 -1.64 1.65 -1.63 4.67 -4.79 4.78 -4.74 6.69 -6.76 Rx (o) -4.94 4.99 -4.91 4.94 -0.35 0.36 -0.34 0.36 -3.37 3.39 Ry (o) -0.29 0.32 -0.27 0.30 4.59 -4.79 4.77 -4.74 2.87 -2.83 Rz (o) Tx (mm) 0.38 -0.39 0.38 -0.35 0.22 -0.22 0.18 -0.21 0.80 -0.93 0.61 -0.87 -0.17 -0.17 0.04 -0.16 -0.61 0.02 -0.60 0.59 Ty (mm) -0.77 0.63 -0.56 1.77 0.62 -1.63 -0.57 -0.16 0.95 -1.06 Tz (mm) 0.67 -0.23 0.01 -2.04 3.48 -1.51 6.03 -3.98 4.86 -4.03 Rx (o) -5.25 4.80 -4.65 6.84 1.61 -3.21 -1.14 -0.12 -1.91 1.37 Ry (o) 2.10 -0.24 0.53 -0.04 2.14 -0.36 -1.36 1.69 -1.80 Rz (o) ICP Time (s) 1.19 -0.51 -0.39 0.16 -0.32 1.24 -1.20 1.01 -1.32 0.74 79.8 62.6 58.1 54.9 93.4 68.2 102.3 70.4 68.0 111.5 2.70 Final RMS Dist (mm) Tx (mm) 1.17 1.85 1.38 1.75 0.76 1.08 0.55 1.15 0.64 0.94 0.42 0.65 0.35 0.02 0.05 0.47 0.73 0.26 1.45 1.14 Ty (mm) 0.57 0.25 0.25 1.33 0.86 1.88 0.38 0.39 0.30 0.78 Tz (mm) 0.96 1.41 1.64 0.42 1.19 3.28 1.25 0.77 1.82 2.73 Rx (o) 0.31 0.18 0.26 1.90 1.95 3.57 0.79 0.48 1.46 2.02 Ry (o) 2.39 0.57 0.80 0.35 2.45 4.43 2.07 3.39 1.18 1.03 0.81 0.13 0.78 0.52 0.55 1.46 1.38 1.22 2.13 1.68 3.22 3.41 12.07 12.26 3.34 12.15 3.04 12.32 3.87 7.08 2.39 7.33 4.19 8.51 2.73 6.96 2.98 5.86 1.91 5.49 Rz (o) Optotrak| Bounding Box Error (mm) Hausdorff Distance (mm) E.3. ICP Registration Results Optotrak Motion Type Transformation Tx (mm) 152 Appendix F Ellipsoid Fitting Results F.1 Ellipsoid Fitting Results: the Sparse Volume Registration Table F.1 shows the ellipsoid fit results on the contours obtained by segmenting several sagittal US images from each volume data using the sagittal segmentation algorithm. In the table, the ellipsoid parameters (center location, orientation, radii lengths) are reported. We also report the sum of square of the algebraic distance (SSD) from the fit - calculated using equation A.5. The root mean squared of this algebraic distance (RMSD) is reported by dividing SSD with the number of the sample points in the contours and taking the square root of the division results. The ellipsoid fit takes 0.002 ± 0.004 seconds to execute and it produces a fit with root mean square distance of 0.011 ± 0.001 mm. The fit produces, on average, 1100 sample points. The algorithm was implemented in Matlab running on a 1.66 GHz computer with 1 GB RAM. 153 -0.74 -1.25 -1.30 2.50 -0.88 2.73 -0.43 -0.66 -0.96 -0.18 -1.12 -0.42 -0.81 0.04 -0.99 Yc (mm) 44.84 44.91 44.86 44.97 45.08 44.97 44.70 47.83 44.59 48.19 44.94 44.84 44.80 44.70 44.83 Zc (mm) 30.04 29.90 30.13 29.65 30.03 29.57 29.91 29.42 30.03 29.52 29.67 34.22 30.11 34.23 29.61 Ax (o) 1.75 2.81 1.86 3.26 1.84 3.94 3.20 2.27 1.91 1.65 3.99 5.83 1.71 4.58 4.43 Ay (o) 6.99 3.09 1.92 -11.66 5.08 -15.49 -1.08 10.80 -0.15 2.25 8.56 -26.90 6.41 16.46 -5.56 Az (o) Rx (mm) -2.77 -3.88 -2.80 -3.56 -3.77 -3.37 -1.57 -2.74 -3.24 -2.46 -4.69 -4.89 -2.77 0.33 -3.80 27.84 28.38 28.69 30.93 28.25 29.84 27.55 28.22 28.44 28.84 28.62 28.85 28.86 28.41 29.58 Ry (mm) 21.61 21.69 21.59 21.48 21.80 21.60 21.52 20.99 21.37 21.36 21.74 21.57 21.54 21.49 21.53 Rz (mm) Fit time (s) # of points SSD (mm2) RMSD (mm) 30.29 0.022 1150 0.148 29.81 0.001 1160 0.183 30.06 0.001 1147 0.165 29.80 0.001 1043 0.114 30.45 0.001 1146 0.173 29.96 0.001 1043 0.092 29.88 0.002 1175 0.156 30.03 0.001 1060 0.123 30.01 0.001 1149 0.168 30.14 0.001 1051 0.117 29.94 0.001 1167 0.107 29.49 0.001 1036 0.116 30.13 0.001 1141 0.143 29.39 0.001 1037 0.089 30.01 0.001 1161 0.191 0.011 0.013 0.012 0.010 0.012 0.009 0.012 0.011 0.012 0.011 0.010 0.011 0.011 0.009 0.013 Vol 1 No Motion Vol 2 Vol 3 Vol 4 X axis translation Vol 5 Vol 6 Vol 7 Vol 8 Y axis translation Z axis translation Vol 9 Vol 10 Vol 11 Vol 12 Vol 13 Vol 14 Vol 15 X axis rotation (yaw) Y axis rotation (pitch) Combination Motion Type Param. Vol 16 Vol 17 Vol 18 Vol 19 Vol 20 Vol 21 Vol 22 Vol 23 Vol 24 Vol 25 Xc (mm) -0.41 -0.64 -0.70 -0.89 1.17 -0.70 0.87 -1.05 0.39 -0.96 Yc (mm) 46.62 44.77 46.50 44.72 44.35 44.81 44.66 44.87 46.49 44.72 Zc (mm) 25.56 29.90 25.77 29.92 35.33 29.69 34.23 30.12 33.67 30.03 Ax (o) -0.55 2.86 -0.34 2.55 0.89 3.10 4.86 2.11 -2.76 2.12 154 Ay (o) 3.76 -7.21 -6.17 6.55 13.02 -2.97 16.71 -7.92 0.68 1.95 Az (o) Rx (mm) -1.15 -1.89 -2.14 -3.39 -3.15 0.09 -5.08 -2.86 -5.91 -2.35 28.59 28.64 27.96 28.34 28.67 28.51 27.44 28.56 27.34 28.06 Ry (mm) 21.56 21.61 21.45 21.57 21.30 21.72 21.67 21.74 21.53 21.79 Rz (mm) Fit time (s) # of points SSD (mm2) RMSD (mm) 30.62 0.001 1032 0.089 29.96 0.001 1158 0.215 30.64 0.001 1032 0.123 30.12 0.002 1156 0.171 29.82 0.002 1014 0.079 30.02 0.001 1163 0.219 29.59 0.001 1045 0.153 30.33 0.002 1145 0.124 30.02 0.002 955 0.147 30.30 0.001 1146 0.177 0.009 0.014 0.011 0.012 0.009 0.014 0.012 0.010 0.012 0.012 F.1. Ellipsoid Fitting Results: the Sparse Volume Registration Table F.1: Ellipsoid fit results for the sparse volume registration. Motion Type Param. Xc (mm) F.2. Ellipsoid Fitting Results: the Full Volume Registration F.2 Ellipsoid Fitting Results: the Full Volume Registration Table F.2 shows the ellipsoid fit results on the contours obtained by using the extended sagittal segmentation algorithm. In the table, the ellipsoid parameters (center location, orientation, radii lengths) are reported. We also report the sum of square of the algebraic distance (SSD) from the fit - calculated using equation A.5. The root mean squared of this algebraic distance (RMSD) is reported by dividing SSD with the number of the sample points in the contours and taking the square root of the division results. The ellipsoid fit takes 0.014 ± 0.013 seconds to execute and produces a fit with root mean square distance of 0.008 ± 0.0005 mm. The fit produces, on average, 13778 sample points. The algorithm was implemented in Matlab running on a 1.66 GHz computer with 1 GB RAM. 155 Vol 4 X axis translation Vol 5 Vol 6 Vol 7 Vol 8 Y axis translation Z axis translation Vol 9 Vol 10 Vol 11 Vol 12 Vol 13 Vol 14 Vol 15 2.11 -1.48 1.97 -1.19 -1.04 -1.40 -1.17 -1.43 -1.39 -1.61 -0.95 -1.32 44.75 44.72 44.81 44.65 48.12 44.63 47.90 44.58 44.25 44.69 44.43 44.46 Zc (mm) 29.81 29.86 29.79 29.14 29.87 30.02 29.47 29.56 30.16 29.06 29.82 33.72 30.03 34.96 29.90 Ax (o) 5.04 4.76 4.45 6.80 5.41 3.75 6.48 3.20 3.56 4.91 5.27 7.36 4.48 3.19 4.75 Ay (o) 6.55 1.75 5.27 2.37 6.11 1.78 7.08 2.20 2.20 4.69 5.54 0.09 3.53 1.86 1.03 Ax (o) Rx (mm) -0.78 -1.46 -2.10 -0.74 -0.26 -1.05 0.13 -1.31 -1.85 -0.83 -1.47 -1.93 0.41 -0.86 -2.06 25.83 25.31 25.15 25.80 25.75 25.69 25.39 25.92 25.11 25.23 25.55 25.39 25.06 24.80 25.56 Ry (mm) 21.61 21.64 21.52 21.53 21.49 21.53 21.57 21.41 21.65 21.34 21.54 21.28 21.46 21.46 21.43 Rz (mm) Fit time (s) # points SSD (mm2) RMSD (mm) 30.40 0.077 14650 1.002 30.48 0.020 14271 1.035 30.58 0.011 14033 1.019 30.11 0.012 14466 0.924 30.60 0.013 14340 0.871 31.05 0.012 14002 1.050 30.21 0.013 14437 1.022 30.65 0.011 13269 0.667 30.49 0.012 14098 1.011 30.11 0.011 13332 0.708 30.57 0.012 14351 0.928 28.79 0.011 13603 0.750 30.67 0.011 13947 0.864 30.13 0.010 12849 0.926 30.71 0.012 14211 1.003 0.008 0.009 0.009 0.008 0.008 0.009 0.008 0.007 0.008 0.007 0.008 0.007 0.008 0.008 0.008 X axis rotation (yaw) Y axis rotation (pitch) Combination Motion Type Param. Vol 16 Vol 17 Vol 18 Vol 19 Vol 20 Vol 21 Vol 22 Vol 23 Vol 24 Vol 25 Xc (mm) -0.74 -1.25 -0.92 -1.20 0.30 -1.10 0.58 -1.03 0.11 -1.05 Yc (mm) 46.17 44.38 46.11 44.39 43.97 44.34 44.40 44.31 45.90 44.07 Zc (mm) 25.20 30.33 25.20 29.99 34.16 29.74 34.56 29.80 33.46 29.82 Ax (o) 1.02 3.28 1.02 4.29 6.24 4.89 3.50 4.86 0.53 5.11 156 Ay (o) 5.19 3.38 3.67 4.19 7.40 3.84 7.71 2.00 5.83 1.75 Ax (o) Rx (mm) -0.08 -1.84 -1.21 -1.76 -2.79 -0.85 -2.26 -1.37 -3.66 -0.90 24.92 24.92 24.69 25.23 24.21 24.66 25.28 24.75 23.82 24.73 Ry (mm) 21.20 21.41 21.24 21.45 21.30 21.51 21.59 21.44 21.27 21.32 Rz (mm) Fit time # of points SSD (mm2) RMSD (mm) 31.84 0.010 12916 0.747 0.008 30.82 0.011 13896 0.974 0.008 31.52 0.010 13122 0.968 0.009 30.56 0.011 14173 0.954 0.008 29.36 0.010 12998 0.818 0.008 30.54 0.012 14050 1.052 0.009 29.63 0.011 13140 0.740 0.008 30.75 0.012 13947 1.032 0.009 30.30 0.011 12368 0.783 0.008 30.48 0.012 13971 0.934 0.008 F.2. Ellipsoid Fitting Results: the Full Volume Registration Table F.2: Ellipsoid fit results for the full volume registration. Motion Type No Motion Param. Vol 1 Vol 2 Vol 3 Xc (mm) -1.47 -1.10 -1.39 Yc (mm) 44.69 44.64 44.57 Appendix G Determining the Motor Specifications The Brachyguide robot is made of two linear stages and two rotation stages. Each stage is driven by a direct current (DC) motor with different torque requirements depending on the stage’s load. The motors were selected based on size constraints. First, their diameter must be less than or equal to 25 mm to keep the robot design compact and light. Second, the motors and gears configuration satisfy the load torque requirement. The following subsections explain how the motor torque requirements were estimated. G.1 The Rotation Stages The rotation stages’ static friction torques are the main load torques for the motors since the stages are already geared with high ratio bevel gears. Ten turns of the stages’ knobs (3600o ) rotates the stages’ plates by 30o . This high gear ratio makes any load torque on the stage’s plate negligible. To estimate the rotation stage static friction torque, an experiment - illustrated in Figure G.1 - was conducted. An object with a known mass was attached to the stage’s knob with a string. Different objects with increasing masses were placed until the stage reached “impending motion state” - when the knob was about to turn. The object’s weight multiplied by the knob’s radius gave the estimated value of the static friction torque: 3.4 mNm. The 1727024 DC micromotor from Faulhaber, satisfies this requirement and is used to drive the rotation stages. G.2 The X-axis Linear Stage For this stage, the DC motor has to generate enough torque to overcome the static friction of the stage’s slider, the friction in the M6 lead screw transmission, and the friction in the motor gear (if any gears are used). We assumed the latter is negligible. 157 G.2. The X-axis Linear Stage Figure G.1: An object with known mass is attached to the stage’s knob with a string to estimate the static friction torque. 158 G.2. The X-axis Linear Stage We conducted an experiment to estimate the slider’s static friction coefficient. An object with known mass (object A) was put on top of the stage’s slider. The stage’s slider generated a normal force equal to the combined weight of the slider and object A. Another object with known mass (object B) pulled the slider through a pulley system - see Figure G.2 for illustration. Different objects with increasing masses were substituted as object B. When the slider was at the stage of impending motion, the mass of object B was recorded and its weight was calculated. This weight divided by the normal force gave the estimated friction coefficient: 0.1637. The pulley easily turned when an object with 8 grams of mass ( 0.08 N weight) was attached to it. The maximum recorded static friction of the slider was about 3 N (almost forty time as large). Hence, in our estimation experiment, the static friction of the pulley is considered to be negligible. Figure G.2: Estimating the static friction coefficient of the stage’s slider using a pulley system. In the robot design, the X-axis slider carries the parallelogram linkage and the two rotation stages. The combined weight of these components are 2 kg at most, corresponding to a 19.6 N normal force exerted by the slider. 159 G.3. The Y-axis Linear Stage Multiplying this number with the estimated static friction coefficient, gives the estimated static friction force of the slider: 3.209 N. To calculate the friction torque of the M6 lead screw, we used the rotary to linear equation (equation G.1) from the PIC Design online catalogue [73]. The efficiency of the lead screw is estimated to be 44% based on the efficiency of a similar Acme lead screw in the catalogue. T orque = Load · lead 2 · π · ef f iciency (G.1) Where the load is 19.6 N, the lead is 1 mm, and the efficiency is 44 percent. The estimated lead-screw friction torque is 7.09 mNm. The DC motor torque requirement is calculated based on the conservation of energy principle (see equation G.2). Ff riction · l = (Tmotor − Tleadscrew ) · θmotor (G.2) Where Ff riction is 3.209 N, Tleadscrew is 7.09 mNm, l is 1 mm, and θmotor is 2π radians. The last two parameters are obtained from the lead value of the M6 lead screw. The required DC motor torque is 7.6 mNm. The 2342-024 DC micromotor from Faulhaber, satisfies this requirement as it has a stall torque of 85.4 mNm and a maximum operating torque of 16 mNm. The motor’s diameter is 23 mm (below our constraint value). As a safety factor, a 3.125 ratio gear is installed on the motor. G.3 The Y-axis Linear Stage The DC motor for the Y-axis stage must generate enough torque to overcome the static friction of the stage’s slider, the friction in the M6 lead screw transmission system, and the combined weight of the robot’s wrist and the X-axis stage. The Y-axis linear stage’s static friction was estimated using a digital balance. The stage’s slider was held manually in a vertical position and the stage was allowed to slide freely until the bottom side touched the balance, see Figure G.3 for illustration. The mass measured by the balance, minus the stage’s mass, multiplied by the gravitational acceleration gave the friction force. This crude estimation approach may over estimate the true value of the friction force. The estimated friction force of the Y-axis slider was 4.6021 N. 160 G.3. The Y-axis Linear Stage Figure G.3: Estimating the static friction forcr of the stage’s slider using a digital balance. 161 G.3. The Y-axis Linear Stage The Y-axis stage supports the combined weight of the parallelogram linkage, the two rotation axes, and the X-axis stage. The combined mass of these components are 2.5 kg at most (equivalent to 24.5 N weight). Substituting this value into equation G.1, we estimate the lead screw friction torque for the Y-axis stage to be 8.86 mNm. Then, the required DC motor torque for the Y-axis stage can be calculated as (Ff riction + Fweight ) · l = (Tmotor − TScrew−T hread ) · θmotor (G.3) where Ff riction is 4.6021 N, Fweight is 24.5 N, TScrew−T hread is 8.86 mNm, l is 1 mm, and θmotor is 2π radians. The required motor torque is 13.5 mNm. The 2342-024 DC micromotor from Faulhaber, satisfies this requirement as it has a stall torque of 85.4 mNm and a maximum operating torque of 16 mNm. The motor’s diameter is 23 mm (below our constraint value). As a safety factor, a 3.125 ratio gear is installed on the motor. 162 Appendix H Initial Ellipsoid Fit Sensitivity Analysis To explore the role of the Kalman filter based segmentation in the final ellipsoid fit, we fit the fourteen user initialization points into ellipsoids and observed the sensitivity of their parameters. Recall, in the extended segmentation algorithm, the user selected fourteen initialization points that guided the IMMPDAF edge detector. We used the same data set from Section 4.4.2. Figure H.1 to H.3 shows how the parameters of the ellipsoid fit to the user initialization point change as one of the initialization points is varied in X-Y-Z axes direction. 163 Appendix H. Initial Ellipsoid Fit Sensitivity Analysis Figure H.1: Effects of changing one initialization point in X-axis direction to the ellipsoid parameters when the user initialization points are fit into an ellipsoid . The plots show linear fits of the observation along with its R2 (the square of the correlation coefficient) and the line equation. 164 Appendix H. Initial Ellipsoid Fit Sensitivity Analysis Figure H.2: Effects of changing one initialization point in Y-axis direction to the ellipsoid parameters when the user initialization points are fit into an ellipsoid . The plots show linear fits of the observation along with its R2 (the square of the correlation coefficient) and the line equation. 165 Appendix H. Initial Ellipsoid Fit Sensitivity Analysis Figure H.3: Effects of changing one initialization point in Z-axis direction to the ellipsoid parameters when the user initialization points are fit into an ellipsoid . The plots show linear fits of the observation along with its R2 (the square of the correlation coefficient) and the line equation. 166 Appendix I Ethics Certificates for Full Board Approval This appendix contains the UBC Research Ethics Board certificates of Approval for the use of patients’ prostate volume data in this study. 167 Appendix I. Ethics Certificates for Full Board Approval Figure I.1: Ethics Approval Certificate for the use of patients’ prostate volume data for the year 2007 - 2008. 168 Appendix I. Ethics Certificates for Full Board Approval Figure I.2: Ethics Approval Certificate for the use of patients’ prostate volume data for the year 2008 - 2009. 169 Appendix I. Ethics Certificates for Full Board Approval Figure I.3: Page one of the Ethics Approval for the use of patients’ prostate volume data for the year 2009 - 2010. 170 Appendix I. Ethics Certificates for Full Board Approval Figure I.4: Page two of the Ethics Approval for the use of patients’ prostate volume data for the year 2009 - 2010. 171
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A robotic needle guide for prostate brachytherapy with...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A robotic needle guide for prostate brachytherapy with pre-operative to intra-operative prostate volumes… Prananta, Thomas Diego 2010
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | A robotic needle guide for prostate brachytherapy with pre-operative to intra-operative prostate volumes registration |
Creator |
Prananta, Thomas Diego |
Publisher | University of British Columbia |
Date Issued | 2010 |
Description | The conventional prostate brachytherapy approach is limited by needle positioning accuracy, needle trajectory option, and prostate motion and deformation between the pre-operative volume study and the seed implant procedure. These limitations increase the risks of post implant complications. In this thesis we develop a robotic needle guide to improve prostate brachytherapy needle placement accuracy and trajectory option as well as a pre-operative to intra-operative prostate volume registration algorithm to address the issue of prostate motion and deformation. Our four degrees of freedom robot provides X-Y axes translational accuracy of 0.12 and 0.1 mm compared to the 5 mm accuracy of the standard needle guide. The robot also provides yaw and pitch angulations with 0.05 degrees accuracy which can be used to reach prostate regions blocked by pubic arch interference. The robot is adaptable to conventional brachytherapy apparatus without adding the clinical procedure time and can be used manually in the case of electronic control failure. The registration approach is based on fitting prostate surfaces into ellipsoids. Pre-operative and intra-operative sagittal view-based volume data are contoured using a novel semi automatic sagittal view-based segmentation algorithm. The resulting contours are fit into ellipsoids whose parameters - centers, orientations, and radii - are used to calculate the registration matrix. The accuracy of the registration algorithm was compared with Optotrak measurement as the gold standard and with the Iterative Closest Point (ICP) algorithm. The result shows that the orientation of the ellipsoid fit is sensitive to user initialization points causing up to 5 mm translational errors and 5.5 degrees angular error. The comparison with ICP shows that the ellipsoid fitting based algorithm is faster but less accurate. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-08-31 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0065853 |
URI | http://hdl.handle.net/2429/28037 |
Degree |
Master of Applied Science - MASc |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2010-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2010_fall_prananta_thomas.pdf [ 15.51MB ]
- Metadata
- JSON: 24-1.0065853.json
- JSON-LD: 24-1.0065853-ld.json
- RDF/XML (Pretty): 24-1.0065853-rdf.xml
- RDF/JSON: 24-1.0065853-rdf.json
- Turtle: 24-1.0065853-turtle.txt
- N-Triples: 24-1.0065853-rdf-ntriples.txt
- Original Record: 24-1.0065853-source.json
- Full Text
- 24-1.0065853-fulltext.txt
- Citation
- 24-1.0065853.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}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0065853/manifest