Towards Real-Time Registration of Ultrasound and CT in Computer Aided Orthopaedic Surgery Applications by Anna B Marie Brounstein B.A.Sc., The University of British Columbia, 2012 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF Master of Applied Science in The Faculty of Graduate Studies (Biomedical Engineering) The University Of British Columbia February 2012 © Anna B Marie Brounstein 2012 Abstract Pelvic fractures are serious injuries that are most commonly caused by mo- tor vehicle accidents and affect people of all ages. Surgeries to realign the pelvis and fix the bone fragments with screws have inherent risks and rely on cumbersome intra-operative radioscopic imaging methods. Ultrasound (US) is emerging as a desirable imaging modality to replace fluoroscopy as an intra-operative tool for pelvic fracture surgery because it is safe, portable and inexpensive. Despite the many advantages of US, it suffers from speckle noise, a limited field of view and a low signal-to-noise ratio. Therefore, we must find a way to efficiently process and utilize ultrasound data so that it can be used to effectively visualize bone. In the past decade, there has been much research focused on fusing US with pre-operative Computed To- mography (CT) to be used in an intra-operative guidance system; however, current methods are either too slow or not robust enough to use in a clinical setting. We propose a method to automatically extract bone features in US and CT volumes and register them using a fast point-based method. We use lo- cal phase features to estimate the bone surfaces from B-mode US volumes. We simplify the bone surface using particle simulation, which we optimize using the hierarchical Barnes-Hut algorithm. To ensure the point cloud best represents the bone surface, we reinforce them with high curvature features. We then represent the point clouds using Gaussian Mixture Models (GMMs) and find the correspondence between them by minimizing a measure of dis- tance between the GMMs. We have validated our proposed algorithm on a phantom pelvis and clinical data acquired from pelvic fracture patients. We demonstrate a registration runtime of 1.4 seconds and registration error of 0.769 mm. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Acronyms . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii 1 Introduction and Background . . . . . . . . . . . . . . . . . . 1 1.1 Motivation and Problem Statement . . . . . . . . . . . . . . 1 1.2 Thesis Objectives . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Medical Imaging in Pelvic Fracture Surgery . . . . . . . . . . 2 1.3.1 Pre-operative Imaging . . . . . . . . . . . . . . . . . . 3 1.3.2 Intra-operative Imaging . . . . . . . . . . . . . . . . . 4 1.4 Current Pelvic Fracture Surgery . . . . . . . . . . . . . . . . 9 1.5 State-of-the-Art US-CT Registration in Orthopaedic Surgery 16 1.5.1 Intensity-Based Registration Methods . . . . . . . . . 16 1.5.2 Point-Based Registration Methods . . . . . . . . . . . 19 1.6 Thesis Organization . . . . . . . . . . . . . . . . . . . . . . . 21 2 Bone Localization from US and CT . . . . . . . . . . . . . . 24 2.1 Feature Extraction and Point Cloud Creation . . . . . . . . . 25 2.1.1 Local Image Phase Features . . . . . . . . . . . . . . 25 iii Table of Contents 2.1.2 Automatically Optimized Parameters for Log-Gabor Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.1.3 Bone Surface Extraction from CT Volumes . . . . . . 35 2.2 Point Cloud Simplification . . . . . . . . . . . . . . . . . . . 35 2.2.1 Point Cloud Simplification using Particle Simulation . 36 2.2.2 Optimization of Force Calculations . . . . . . . . . . 38 2.3 Point Cloud Reinforcement Using High Curvature Features . 41 3 Fast 3D Correspondence of US and CT Volumes and Its Application in Pelvis Fracture Surgery . . . . . . . . . . . . 51 3.1 Registration of US to CT Volumes using Gaussian Mixture Model Matching . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.1.1 The Application of Gaussian Mixture Models to Point Clouds . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.1.2 Minimizing the Difference between Two Guassian Mix- ture Models . . . . . . . . . . . . . . . . . . . . . . . 55 3.1.3 The effect of Point Cloud Size on Registration Run- time. . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.2 Extension of Methods to US-US Stitching . . . . . . . . . . . 56 3.3 US to CT Registration in a Pelvic Fracture Surgery Guidance System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.3.1 Introducing ultrasound (US) into Pelvic Fracture Surg- eries . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4 Validation of Proposed Registration Algorithm . . . . . . . 66 4.1 Registration Requirements for Proposed Pelvic Fracture Sur- gical System . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.2 Experimental Setup and Results for Phantom Study . . . . . 67 4.2.1 Phantom Creation and Equipment . . . . . . . . . . . 67 4.2.2 Point Cloud Size and the Effect on Registration Run- time and Accuracy . . . . . . . . . . . . . . . . . . . 70 4.2.3 Particle Simulation’s Effect on Registration Robust- ness . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 iv Table of Contents 4.2.4 Registration Results and the Contribution of Curva- ture . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 4.2.5 A Comparison of GMM Registration and Other Point- Based Methods . . . . . . . . . . . . . . . . . . . . . 80 4.2.6 Validation of the US Stitching Adaptation of Regis- tration Algorithm . . . . . . . . . . . . . . . . . . . . 85 4.3 Clinical Evaluation of GMM Registration . . . . . . . . . . . 86 5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 5.1 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 5.2 Thesis Contributions . . . . . . . . . . . . . . . . . . . . . . 92 5.3 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 v List of Tables 1.1 Summary of state-of-the-art orthopaedic registration methods. 22 4.1 Quantitative results demonstrating the effect of point cloud size; specifically 300 and 500 points per dataset. The test was run on 49 US-CT phantom dataset pairs and point clouds were simplified using particle simulation. . . . . . . . . . . . . 72 4.2 Quantitative results of test comparing registration results for tests using particle simulation (PSim) and näıve down sam- pling. The test was run on 49 US-CT phantom dataset pairs and correspond to the boxplots in Figure 4.5. . . . . . . . . . 74 4.3 Quantitative results of test comparing registration results for tests using PSim and down sampling. The test was run on 49 US-CT phantom dataset pairs and correspond to the boxplots in Figure 4.7. . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.4 Comparison of GMM registration, iterative closest point (ICP) and coherent point drift (CPD). For the purposes of this ta- ble, the fail rate was defined as 2 mm. The test was repeated for each of the 49 US datasets described in Section 4.2 and was repeated 10 times for each range tested. . . . . . . . . . . 84 4.5 Comparison of GMM registration, ICP and CPD. Registra- tion algorithms were tests on 49 US-CT data pairs (unlike Table 4.4). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.6 Quantitative results of test comparing registration results for US stitching. Tests were run on nine pairs of US data; the qualitative results can be seen in Figure 4.11. . . . . . . . . 86 vi List of Tables 4.7 Quantitative results for clinical tests comparing automati- cally and empirically optimized parameters phase symme- try (PSym). One non-optmized three dimensional (3D) PSym test failed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 vii List of Figures 1.1 Inlet and outlet fluoroscopy images are taken of the patient prior to surgery. Image courtesy of Wiss et al. [1] . . . . . . . 3 1.2 Example inlet and outlet views of the pelvis. The inlet view (left) shows the pelvic ring while the outlet view (right) pro- vides a view of the sacrum. Image courtesy of the UBC Di- vision of Orthopedic Trauma. . . . . . . . . . . . . . . . . . . 4 1.3 CT volumes are taken pre-operatively to help diagnose frac- tures and plan surgery. The red arrows in the (a) 3D volume and the (b) single two dimensional (2D) slice show the loca- tion of the sacral fracture. . . . . . . . . . . . . . . . . . . . . 5 1.4 Pelvic fracture surgery set up. An inlet view (a) and an outlet view (b) are taken prior to starting surgery to detect any movement in the pelvic fragments. Images courtesy of the UBC Division of Orthopedic Trauma. . . . . . . . . . . . . . 7 1.5 Fluoroscopic image of a pelvis. The black arrows indicate the iliocortical density and the red screw shows the intended screw placement. Image courtesy of Wiss et al. [1] . . . . . . 8 1.6 O-arms can be used to acquire both 2D fluoroscopic images and 3D CT images. (a) shows fluoroscopic shots of a pa- tient’s illium and (b) shows a 2D slice of a CT volume of a patient’s spine and hip. Picture courtesy of Medtronic (www.medtronic.com). . . . . . . . . . . . . . . . . . . . . . . 10 viii List of Figures 1.7 Frontal slice of a sacrum fracture showing K-wire placements. The yellow structures are sacral nerves and the red and blue structures are blood vessels. Both the nerves and blood ves- sels must be avoided during K-wire and cannulated screw insertion. Image courtesy of Wiss [1]. . . . . . . . . . . . . . . 12 1.8 Frontal (a) and sagittal (b) view of the nerves in the sacrum. Image courtesy of Wiss [1]. . . . . . . . . . . . . . . . . . . . 13 1.9 Drilling the K-wire requires several fluoroscopic images. (a) is the anterior-posterior (AP) view of the K-wire being drilled. When the screw reaches the sacral foramen, (b), the surgeon must take a lateral shot (c) to ensure the screw does not come out of the sacrum anteriorly. The wire must therefore be behind the iliocortical density, seen marked as a white line in (d). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.10 Tools used in pelvic fracture surgery: (a) Low profile recon- struction plates are used bent to shape the bone surface and are attached with screws. (b) Cannulated screw on a K-wire. (c) Direct measuring device used to measure insertion depth of K-wire. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.1 Overview of the proposed algorithm’s pipeline. The US Bone Segmentation block is described in Section 2.1.1. The CT Bone Segmentation block is described in Section 2.1.3. The Point Cloud Simplification blocks are described in Section 2.2 and the Reinforcement blocks in Section 2.3. The Registra- tion block is described in Chapter 3. . . . . . . . . . . . . . . 26 2.2 The CT volume is cropped to the region of interest (ROI) and thresholded. The bone surface is extracted using ray-tracing, see Section 2.1.3. The bone in the US volume is localized using PSym and the final bone surface is determined using ray-tracing, Section 2.1.1 . . . . . . . . . . . . . . . . . . . . . 27 ix List of Figures 2.3 Once the bone surface is extracted from both the US and CT data, the points on the surface are converted to a list and simplified using PSim, see Section 2.2. The resulting point cloud is reinforced with high curvature features, Section 2.3. . 28 2.4 Example of a 2D Log-Gabor radial and angular filter compo- nents. (a) and (b) show the radial component at two different scales, as described by Equation 2.3. The angular components of the Log-Gabor filter are seen in (c) - (h) at six different orientations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.5 Three filter angles are demonstrated below to capture three different bone surfaces in the US volume. The first row is a slice from the B-Mode US volume. The second row is a slice from three different angular components of the filter. The bottom row is the angular filter response for each corre- sponding filter. . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.6 An example of the US surface extraction steps used in our proposed algorithm. The top row is a single 2D slice taken from the volume in the bottom row. The B-Mode US can be seen in (a) and (d). The PSym of the B-Mode volume is seen in (b) and (e). Finally, the outline can be seen in (c) and (f). 32 2.7 The filter’s orientation must be steered to cross different seg- ments of the bone surface to capture the entire bone surface. The red arrows in each B-Mode image show the surface being captured with the adjacent angular filter. . . . . . . . . . . . 34 2.8 PSim is used to decrease the size of PS without losing salient information. (a) shows the original dense point cloud, PS . (b) shows the randomly sampled point cloud. FIn is shown in (c) and a close up is shown in (d). (e) shows the relaxed points projected back on the surface described by PS . The final point cloud, Ps, is shown in (f) after three iterations. . . 39 2.9 An example of the algorithm used to build a 2D Barnes-Hut tree. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 2.10 Continuation of example in Figure 2.9. . . . . . . . . . . . . . 43 x List of Figures 2.11 2D example calculating Fs,C for θ = 1.5. . . . . . . . . . . . . 44 2.12 2D example Calculating Fs,B for θ = 1.5. . . . . . . . . . . . 45 2.13 Pelvic bone anatomy and the effect of high curvature rein- forcement - We are primarily interested in imaging the iliac spine and the iliac crest (a). While the simplified point cloud can register with respect to the pelvic table (indicated by the outlined arrow in (b)), high curvature features are needed to register correctly along to the iliac crest (indicated by the solid arrow in (b)). . . . . . . . . . . . . . . . . . . . . . . . . 48 2.14 Gaussian curvature reinforcement: (a) shows a 3D rendering of an US volume acquired from the iliac spine and (b) corre- sponding parametrized surface. The dark background is not part of the surface and is therefore not included in the K calculation. (c) shows the high curvature features, K, on the parametrized surface. (d) is the final simplified point cloud (blue) with reinforced curvature features (red). . . . . . . . . 50 3.1 The proposed registration pipeline using GMM. The inputs of the registration algorithm is the simplified point clouds de- scribing the US and CT bone surfaces as described in Chapter 2. The output is the transformed US volume. . . . . . . . . . 53 3.2 Example Illustrating GMMs - An example of a randomly gen- erated 2D GMM comprised of three components with mean vectors, µi, at [−3,−5], [1, 2] and [5,−4]. . . . . . . . . . . . 54 3.3 Flowchart describing the image processing pipeline for US stitching using our proposed algorithm. . . . . . . . . . . . . 58 3.4 Flowchart of patient’s pre-operative care. . . . . . . . . . . . 61 3.5 Flowchart of proposed operating room (OR) preparation. . . 62 3.6 Flowchart of pelvic reduction procedure in proposed system. . 63 3.7 Flowchart of K-wire insertion procedure in proposed system. 64 3.8 Flowchart of screw placement procedure. . . . . . . . . . . . . 65 xi List of Figures 4.1 The phantom was created with a Sawbone radio-opaque hemi pelvis, cut along red line (a). The final phantom embedded in PVC gel is seen in (b). Sawbone graphic courtesy of Pacific Research Laboratories Inc. . . . . . . . . . . . . . . . . . . . . 68 4.2 A single slice of a US volume (a) and corresponding CT slice (b) of the phantom. The red arrow (a) shows the reverbera- tion caused by the fiducial. The blue arrow (a) shows an air bubble. Notice that the air bubble does not have reverbera- tions and does not appear in the the CT image (b). . . . . . . 71 4.3 To justify the use of particle simulation, we must compare the registration error and overall runtime. As the optimal number of points per point cloud is not known for registration attempts that utilise particle simulation, the test was re-run for different point cloud sizes. PSym was omitted from the runtime calculation, as it is not dependent on the number of points in a point cloud. . . . . . . . . . . . . . . . . . . . . . 73 4.4 Comparison of point cloud simplification methods - The full point cloud extracted from an US bone surface is seen in (a). (b) shows the point cloud simplified with K-means clustering, (c) hierarchical clustering and (d) particle simulation. . . . . 75 4.5 Quantitative results of registration validation, comparing the target registration error (TRE) and surface registration er- ror (SRE) of registration performed on point clouds simpli- fied with particle simulation and down sampling. Note the im- proved registration results for datasets simplified using PSim. The test was run on 49 US-CT phantom dataset pairs and correspond to the results seen in Table 4.2.3 . . . . . . . . . . 76 4.6 The runtimes of force calculations comparing the Barnes-Hut algorithm and exhaustive force calculations. Data was run on a single point cloud for various desired output point cloud sizes. The test was repeated 10 times; the values above are the mean runtime. . . . . . . . . . . . . . . . . . . . . . . . . 78 xii List of Figures 4.7 Quantitative results evaluating the effect of reinforcing a point cloud with high curvature features, K. Note the improvement of the TRE in cases that include K features. The test was run on 49 US-CT phantom dataset pairs and correspond to the results seen in Table 4.2.4. . . . . . . . . . . . . . . . . . 79 4.8 Qualitative results of GMM Registration with PSim and Cur- vature Features. Three tests are shown, with both 3D and a single corresponding 2D slice. (a) is the CT data, (b) is the B-Mode US, prior to registration, (c) is the fused result of the GMM registration and (d) is fused rendering of the extracted bone surfaces. . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 4.9 Comparison of GMM, ICP and CPD registration. The typ- ically acceptable error for a number of orthopaedic surgery applications is between 2 mm (blue line) and 4 mm (red line) for the entire system [2]. This test was run on the 49 phan- tom US point sets previously extracted and correspond to the data seen in Figure 4.10 and Table 4.4. For each dataset, a translation and rotation was applied in a random direction to create the model dataset. . . . . . . . . . . . . . . . . . . . . 82 4.10 Comparison of the initial vs. final TRE values for GMM, ICP and CPD. The error limits in orthopaedic surgery are between 2 mm (blue line) and 4 mm (red line) for the entire system [2]. These results correspond to Figure 4.9 and Table 4.4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.11 Qualitative results showing two of the nine US tests used to validate the proposed algorithm using US Stitching. Quanti- tative results can be seen in Figure 4.12 and Table 4.6. . . . . 87 4.12 Quantitative results comparing the runtime and error of US stitching for points clouds of size 300 and 500. Nine pairs of US volumes were acquired from the phantom for this test and correspond to the qualitative results seen in Figure 4.11 and quantitative results seen in Table 4.6. . . . . . . . . . . . . . 88 xiii List of Figures 4.13 Results of clinical US to CT registration. The top row is the registration results using automatically optimized phase symmetry parameters to localize the US bone surface. The bottom row uses empirical parameters to calculate the phase symmetry. Note that the bottom row is misregistered. Quan- titative results can be seen in Table 4.3. . . . . . . . . . . . . 89 xiv List of Acronyms Abbreviation Definition 2D Two Dimensional 3D Three Dimensional AP Anterior-Posterior CAOS Computer Assisted Orthopaedic Surgery CPD Coherent Point Drift CPU Central Processing Unit CR Correlation Ratio CT Computed Tomography CUDA Compute Unified Device Architecture EMG Electromyogram GMM Gaussian Mixture Models GPU Graphical Processing Unit HU Hounsfield Unit ICP Iterative Closest Point ITK Insight Toolkit NCC Normalized Cross Correlation OR Operating Room PDF Probability Density Function PSim Particle Simulation Psym Phase Symmetry RMS Root-Mean-Square ROI Region of Interest SNR Signal-to-Noise Ratio SRE Surface Registration Error xv List of Acronyms SSM Statistical Shape Model TRE Target Registration Error US Ultrasound VGH Vancouver General Hospital xvi Acknowledgements I would like to thank my supervisors Rafeef Abugharbieh and Antony Hodg- son for all of their time and dedication. Thank you to Ilker Hacihaliloglu for the help he has given me throughout my time here. I would also like to thank my parents for their unwavering support and James for always believing in me. xvii Chapter 1 Introduction and Background 1.1 Motivation and Problem Statement Pelvic fractures are serious injuries that often require high risk surgery to re- align and fix bone fragments. These fractures are most commonly caused by car accidents and sports injuries; however, a simple fall can cause a fracture in the pelvis of an elderly patient. While many elderly patients are treated non-operatively, majority of patients are operated on to improve the quality of life after recovery by realigning and fixing the bone fragments. There are considerable risks, however, involved in such a surgery. The purpose of our work is to minimize the risk to the patient due to error, reduce the exposure to radiation and potentially improve the realignment of the bone fragments. The treatment almost invariably requires reduction surgery to re-align the bone fragments and stabilize the pelvis. Pelvic fractures only account for 3% of all fractures; however, they are usually severe and are an indi- cator of a major trauma with a wide range of complications that must be diagnosed and treated early. Over 40% of patients who require pelvic sta- bilization result in long-term complications, usually neurological, urological and non-specific pain [3]. Computer Assisted Orthopaedic Surgery (CAOS) systems implemented to assist surgeons navigate complex local anatomy during pelvic surgeries could help minimize the complications by providing guidance and insight into the pelvic region during surgery [3]. Pelvic frac- ture surgeries take over one and a half hours, on average, during which the patient and the OR staff are exposed to approximately 11.2 minutes of fluo- 1 Chapter 1. Introduction and Background roscopic radiation [4]. This exposure time does not include the time required to reposition the fluoroscopic C-arm; therefore operating the C-arm is one of the most time consuming steps of a pelvic fracture surgery. Minimizing the fluoroscopy in these procedures would not only decrease the radiation exposure to the patient and staff but could also feasibly decrease the overall time required in the OR. 1.2 Thesis Objectives The aim of this thesis is to provide an algorithm to enable us to integrate pre-operative CT volumes of pelvic fracture patients in a guidance and nav- igation system to assist surgeons in realigning bone fragments and inserting fixation implants. We propose a robust and accurate registration algorithm, which will be a key component in this framework for orthopaedic surgery. This registration algorithm must be sufficiently fast, i.e. not impede the surgeon during surgery. In this thesis, we propose and demonstrate a fully automatic registration algorithm for fusing US and CT volumes. We vali- date this proposed algorithm on phantom and clinical data, using CT as the gold-standard. This work is ultimately intended to contribute to an implementation of the CAOS system described in this thesis. We propose to introduce 3D US into pelvic fracture procedures to minimize the radiation exposure for the OR staff and patients, and register the US scans to pre-operative CT volumes to provide a high quality visualization and guidance system for pelvic fracture surgery. It should be noted that we intend to track bone fragments, rather than the gap between the fragments, known as fractures. 1.3 Medical Imaging in Pelvic Fracture Surgery In this section we describe the medical imaging modalities currently used to diagnose pelvic fractures and view the pelvis during surgery. 2 Chapter 1. Introduction and Background 1.3.1 Pre-Operative Imaging When a patient is admitted to a hospital with a suspected pelvic fracture, standard AP view X-rays are taken. If a fracture is observed in these images, inlet and outlet view X-rays are taken to further delineate the injury of the patient. An inlet image is the superior view of the pelvis and is acquired by positioning the X-ray source at 45◦ from the body towards the caudal (inferior) direction. The inlet image is taken to visualize the pelvic ring. The outlet view is taken at 45◦ in the cephalad (superior) direction [5] and shows the sacrum, pubic rami and obturator foramen [3]. An example of the inlet and outlet image directions of a patient in supine position are seen in Figure 1.1 and inlet/outlet fluoroscopic images of the pelvic can be seen in Figure 1.2. Life-threatening injuries can be detected from the AP views, and additional details are gained from the inlet and outlet views. It is also advantageous to use the inlet and outlet views prior to surgery, as they are the projections viewed intra-operatively. If the patient is stable and does not require immediate treatment, they will then undergo further diagnostic imaging. Figure 1.1: Inlet and outlet fluoroscopy images are taken of the patient prior to surgery. Image courtesy of Wiss et al. [1] 3 Chapter 1. Introduction and Background Figure 1.2: Example inlet and outlet views of the pelvis. The inlet view (left) shows the pelvic ring while the outlet view (right) provides a view of the sacrum. Image courtesy of the UBC Division of Orthopedic Trauma. If a fracture is observed in the 2D X-ray, the patient will undergo a CT scan. The surgeon uses this CT information to decide on surgical or non-operative care. If the surgical route is chosen, this CT will be used for pre-operative planing. The CT volume is often rendered and displayed in axial cuts with increased use of sagittal, coronal views as well as a 3D reconstructions during the surgery. Figure 1.3 shows an example CT volume of a pelvis with a sacral fracture, indicated by a red arrow. A pelvic fracture will often be accompanied by other injuries. Radio- graphic imaging methods are used by surgeons to diagnose fractures and other injuries, as many other organs are present in the pelvis. Modalities such as CT, US and electromyogram (EMG) may be used to assist the di- agnosis of soft-tissue trauma [3]; however US is not currently being used to specifically assist fracture care. We do not use this pre-operative US in our proposed method. 1.3.2 Intra-Operative Imaging The pre-operative CT is displayed for the surgeons during surgery, although the standard procedure is to only use this data to indicate the location of the 4 Chapter 1. Introduction and Background Figure 1.3: CT volumes are taken pre-operatively to help diagnose fractures and plan surgery. The red arrows in the (a) 3D volume and the (b) single 2D slice show the location of the sacral fracture. 5 Chapter 1. Introduction and Background bone fragments. The wait time between patient intake and surgery may be anywhere between one day to three weeks after the CT, unless immediate attention is needed. During this time, the pelvic fragments are liable to move and the CT will no longer reflect the current position of the bone fragments. Surgeons use fluoroscopic C-arms to assess alignment during surgery to track how well fragments are reduced and to guide fixation implant place- ment. In most ORs, a radiology technologist operates the C-arm during surgery to display the results of the realignment and screw placement. Sev- eral fluoroscopic images are taken from multiple angles to visualize the loca- tion of the screws; an experienced surgeon needs as many as 20 fluoroscopy images to visualize a single screw. Surgeons have identified anatomical land- marks on radiographic images to assist the safe placement of implants. For example, Figure 1.5 shows the iliocortical density (indicated by the black arrows), an anatomical feature that demonstrates the target depth of the iliosacral screw placement. Subtle features like these are used by surgeons to place the K-wires and pins and often require several images from different angles to locate [1]. If a patients receives over 2 Gys, they are at risk of skin injuries due to hyperemia of capillaries, known as Erythema. For absorbed dosages over 10 Gys, the patient is recommended to follow up with a doctor to ensure Erythema does not occur. These injuries are not immediate and can take anywhere between two to ten weeks to appear [6]. While the risk of cancer and other long-term illness is not high for the patient, the surgeons are ex- posed to significantly more ionizing radiation over the course of their career. The surgeons and staff are more likely to get cataracts, cancer and other diseases as a results of their work conditions. In particular, surgeons and staff are more susceptible to brain cancers such as glioma because they do not wear head shields [7]. Some ORs are equipped with O-arms — portable CT machines dis- tributed and sold by Medtronics that are placed around the patient during surgery. These machines are cumbersome and expose the patient to similar doses of radiation are produced by first generation CTs. However, they pro- 6 Chapter 1. Introduction and Background Figure 1.4: Pelvic fracture surgery set up. An inlet view (a) and an outlet view (b) are taken prior to starting surgery to detect any movement in the pelvic fragments. Images courtesy of the UBC Division of Orthopedic Trauma. 7 Chapter 1. Introduction and Background Figure 1.5: Fluoroscopic image of a pelvis. The black arrows indicate the iliocortical density and the red screw shows the intended screw placement. Image courtesy of Wiss et al. [1] 8 Chapter 1. Introduction and Background duce a relatively high quality volume of the pelvis which allows the surgeon to better visualize the bone alignment and screw placement. This machine is not available in most ORs due to the high capital purchase cost of the machine. The O-arm is also capable of creating 2D fluoroscopic image, an example of intra-operative CTs and a 2D fluoroscopy image obtained with an O-arm can be seen in Figure 1.6. Both C-arm and O-arm images have been linked to surgical navigation systems. These systems have shown limited uptake in trauma related to their cost but are important as they show the phenomenon where fracture fragments move between the pre-operative CT and surgery. 1.4 Current Pelvic Fracture Surgery When a patient arrives at a hospital with a suspected pelvic fracture, the patient first undergoes an X-ray to determine if (1) the patient has a frac- tured bone and (2) if a CT should be taken to further diagnose the patient. If the patient is in critical condition and is not hemodynamically stable, he/she will be operated on immediately, bypassing the CT. Otherwise, an orthopaedic surgeon will use the CT to find the fracture location and deter- mine if surgery is necessary. If surgery is necessary, the surgeon will then plan the surgery; i.e. the surgery to realign the bone fragments and possibly insert screws and plates. If the patient is in critical condition, they may require immediate action; however, most patients will wait at least one day. The patient’s wait time depends on the judgement of the surgeon, among other factors some surgeons will operate within a few days to minimize the patient’s recovery time while other surgeons will wait up to three weeks to allow the patient to heal after their trauma to minimize complications. During this time, the pelvic fragments will move with respect to both the patient and other bone fragments, therefore the CT will be inaccurate when the surgery begins, and cannot be used for guidance without being updated. Once the patient enters the OR, he/she are placed under general anaes- thetic by an anaesthesiologist while the OR is being prepared for surgery. This preparation involves setting up tools, equipment and the sterile work 9 Chapter 1. Introduction and Background Figure 1.6: O-arms can be used to acquire both 2D fluoroscopic images and 3D CT images. (a) shows fluoroscopic shots of a patient’s illium and (b) shows a 2D slice of a CT volume of a patient’s spine and hip. Picture courtesy of Medtronic (www.medtronic.com). 10 Chapter 1. Introduction and Background environment. The CT volumes are loaded on displays in the OR and the surgeons will often study these images before the surgery to do any last minute planning. This preparation step can take up to one hour between the time the patient enters the OR and the time the surgeon begins the surgery. In general, a patient is placed in supine position for closed surgery and prone position for open surgery, although this is not necessarily true for all patients. The first step of the surgery is to reduce the pelvis, i.e. to realign the pelvic fragments. Traction is commonly used to align the pelvis, either being applied by an assistant or a weight. At this point, the surgeon will send a request for a fluoroscopic C-arm which must be operated by a radiology technologist. Every person in the operating room must wear a lead vest and thyroid shield and will continue to wear it until the end of the surgery. These lead vests are very heavy and are not only an inconvenience, but can potentially cause back or muscular injuries with long term use over the career of the surgeon and staff [7]. The technologist will take fluoroscopic shots of the adjusted pelvis from several different angles and the surgeon will judge if the pelvis is adequately aligned. If the pelvis needs further alignment, they will repeat the process of reducing and imaging the pelvis until it is aligned. This is an iterative process and will likely happen several times. The next step is to insert K-wires. These are long wires that are drilled into the pelvis to temporarily secure the pelvis and act as guides for can- nulated screws. The surgeon holds the drill and K-wire at the intended insertion point and asks the radiology technologist to take several fluoro- scopic shots around the tip of the tool to ensure the K-wire’s trajectory and tip location will fix the bone fragments together. Figure 1.7 shows an illus- tration of the K-wire placement with respect to major nerves and vascular structures. This takes several radiographic shots from different angles to de- termine the trajectory of the wire. The surgeon then drills the K-wire into place and checks the placement using another set of fluoroscopic images. The iliocortical density is an anatomical feature observed in fluoroscopic images that surgeons use to ensure the screw does exit the pelvis in the anterior 11 Chapter 1. Introduction and Background Figure 1.7: Frontal slice of a sacrum fracture showing K-wire placements. The yellow structures are sacral nerves and the red and blue structures are blood vessels. Both the nerves and blood vessels must be avoided during K-wire and cannulated screw insertion. Image courtesy of Wiss [1]. direction, Figure 1.9. These steps are seen in Figure 3.3.1. The surgeon will repeat this until he or she is satisfied that enough wires are placed to fix the bone fragments. The surgeon will use a direct measuring device to measure the K-wire’s depth in the bone, as seen in Figure 1.10. Cannulated screws (Figure 1.10) are over-drilled on the K-wires to prop- erly fix the bone and the K-wires are removed. If necessary, the surgeon will also install reconstruction plates on the surface of the pelvis to provide fur- ther fixation for the patient. The surgeon evaluates the work with a last set of fluoroscopic images and either finishes the surgery or repeats the K-wire or screw insertion steps. The surgery described above requires on average 12 minutes of X-ray exposure time. Our proposed system will not completely eliminate the need of this radiographic modality, but it will minimize its necessity. 12 Chapter 1. Introduction and Background Figure 1.8: Frontal (a) and sagittal (b) view of the nerves in the sacrum. Image courtesy of Wiss [1]. 13 Chapter 1. Introduction and Background Figure 1.9: Drilling the K-wire requires several fluoroscopic images. (a) is the AP view of the K-wire being drilled. When the screw reaches the sacral foramen, (b), the surgeon must take a lateral shot (c) to ensure the screw does not come out of the sacrum anteriorly. The wire must therefore be behind the iliocortical density, seen marked as a white line in (d). 14 Chapter 1. Introduction and Background Figure 1.10: Tools used in pelvic fracture surgery: (a) Low profile recon- struction plates are used bent to shape the bone surface and are attached with screws. (b) Cannulated screw on a K-wire. (c) Direct measuring device used to measure insertion depth of K-wire. 15 Chapter 1. Introduction and Background 1.5 State-of-the-Art US-CT Registration in Orthopaedic Surgery This section describes research efforts in US-to-CT registration and other relevant research in image registration. The research in this area generally falls into one of two categories: intensity-based registration or feature-based point-to-point registration. Both methods have merits, and we will outline the current research of both methods. Registration is almost exclusively rigid because the shape and size of the bone is not expected to change. For certain cases, (e.g. the spine, bone fragments), the US and CT volumes are registered using articulated methods, as the bone segments move with respect to each other. Most 2D-to-3D registration techniques are used to register fluoroscopy to CT or MRI images. Unfortunately, many of these papers focus on methods using DRR (digitally reconstructed radiographs) which is not applicable to 2D US or other non-projected images [8–10]. 1.5.1 Intensity-Based Registration Methods One leading group in US-to-CT registration for orthopaedic surgery is at Ruhr-Universität Bochum, lead by B. Brendel and S. Winter; they focus on intra-operative US imaging for spinal surgeries [11–14]. Brendel and Winter propose to register pre-operative CT volumes of the spine to 3D US volumes of vertebra. Their research focuses on intensity-based registration algorithms, with minimal pre-preprocessing of US images. Winter et al. do not segment the US bone surfaces to avoid lengthy image processing steps and issues of automaticity and robustness. Instead they use adaptive depth gain compensation on 3D US images of vertebrae to emphasize the bone surfaces. The CT volumes are pre-operatively extracted by thresholding the CT image and then identifying the surface visible from the skin, i.e. the bone surfaces that provide a response in B-Mode US [13]. Brendel et al. register the US images using a set of points defined by the extracted CT bone surface and optimize the translation ⇀ s and the rotation R using the 16 Chapter 1. Introduction and Background optimization criterion, for the number of points N , and the gray value u( ⇀ p ) for each 3D coordinate of point ⇀ pn: o = 1 N N∑ n=1 u(R · ⇀pn + ⇀s ) (1.1) This method is used for 3D US to 3D CT registration in [13] and further described in [11, 14]. However it is susceptible to errors due to the compar- ison of a well defined bone surface in CT volumes due to a blurred partial surface in US volumes. Brendel et al. validate this method in [12] and Win- ter et al. further improves this method in [15] by proposing four optimization methods: a gradient-based method and three separate evolutionary-based methods with an average registration runtime of five seconds. Winter pro- poses an US acquisition protocol to improve accuracy of this registration method [16]. Their method shows very good results, with the best method successfully converging 100% but the runtime is not reported. Winter et al. imply, however, their method with the best accuracy is not suitable for real-time applications. Penney et al. introduced a new registration method, that creates two intermediate volumes from the US and CT datasets and registers them us- ing Normalized Cross Correlation (NCC). The intermediate volumes are described by probability density functions (PDF)s of the bone surface from each modality. This method showed promising results, but it requires many test volumes to create the intermediate volumes and therefore it is not fea- sible in a clinical setting [17]. US simulation from CT datasets was proposed by Wien et al. for the purpose of registration in [19] and [20]. Their main contribution is simulat- ing the physical principles of US based on tissues detected in pre-operative CT volumes and optimizing the runtime as previous attempts at US sim- ulation required several weeks to run [21, 22] down to seconds [23]. The tissue density can be estimated from the Hounsfield Unit (HU) in the CT volumes and since the tissue density is directly proportional to the acoustic impedance, the transmitive and reflective properties of the tissue can be 17 Chapter 1. Introduction and Background calculated. The similarity metric used for the registration is the correlation ratio (CR), which is similar to NCC without the assumption of a linear map- ping between the intensities of the two volumes to account for inaccuracies in the US simulation due to unknown large-scale reflections and homogeneity signal. The direction and type of probe is factored in to find the simulated US intensity from these transmitivity and reflectivity properties. The av- erage fiducial Registration Error (FRE) and TRE are 10.4 mm and 9 mm for rigid registration with a registration runtime of approximately twenty seconds (note the FRE and TRE are due to the algorithm being tested on soft tissue). In 2008, Shams et al. [24] propose a more complex simulation model that incorporates Lambertian scattering into the simulated CT vol- ume and requires 2 days to simulate a B-Mode ultrasound image with 128 RF scans. Khallaghi et al. optimized Wein’s US simulation and validated their implementation on the lumbar spine [25] by creating an US simula- tion from a statistical shape model (SSM) and registering it to a 3D US volume. They report an accuracy of 3.25mm and a simulation runtime of 8 ms. Kutter et al. [23] optimize Wein’s and Sham’s simulation method on the Graphical Processing Unit (GPU) and validated it for image fusion and simulator training, respectively. They show up to a 300x speedup factor using Wein’s simulation method and reduce the run time of Sham’s method from two days to 6 ms. Note this speed up is only for the US simulation pre-processing, not the registration algorithm. The registration runtime is still over 20 minutes, as they use an intensity based method. In the past five years, over twenty-five papers proposing registration op- timization on the GPU have been published; however, few of these focus on rigid registration methods or US applications [26]. Reichl et al. [27] incorporate speckle noise, acoustic absorption and blurring into the simula- tion framework and implement the algorithm using Compute Unified Device Architecture (CUDA), requiring 8.7 ms to process a 512x512 CT slice. It should be noted that these papers apply the proposed method to soft-tissue and use deformable registration methods, and therefore these results cannot be directly compared to our work. We mention these papers; however, be- cause recent work has applied this US simulation technique to orthopaedic 18 Chapter 1. Introduction and Background registration problems. Gill et al. [28] validated Wein’s simulation model on the lumbar spine of a sheep cadaver by registering the US volumes to the simulated US volumes using the CR similarity metric used by Wein et al. [19, 20]. They imposed biomechanical constraints to prevent registra- tion results that are anatomically impossible and report a TRE of 1.44 mm and 1.25 mm for a phantom and cadaver, respectively. Their success rate is high for the phantom study at 98.5%, but the success rate is lower at 87% for the cadaver study, likely due to the introduction of soft tissue. While their results are promising, their registration runtime average is 43 minutes when implemented with the Insight Toolkit (ITK) and run on the central processing unit (CPU) and is still 4 minutes long when implemented on a GPU. Ghanavati et al. created a statistical atlas from 110 male pelvis CTs generated using principal component analysis and simulated US from this model. They used both rigid and deformable Linear Correlation of Linear Combinations to register the simulated US with several 2D US slices. While the results were quick to obtain, the TRE was reported to be 9.15 mm and 7.21 mm in their cadaver study [18]. These intensity-based registration methods give accurate results, how- ever their computational runtime makes them impractical to be used during surgery. Despite the advances made in the US simulation pre-processing steps by Gill et al., the actual registration of orthopaedic images takes up to twenty minutes [28]. 1.5.2 Point-Based Registration Methods A relatively successful point-based registration method called ICP has been used in various orthopaedic applications [29, 30]. While the results of ICP are fast and often accurate, the robustness of the algorithm suffer due to local minima. To overcome this robustness issue, Penney et al. [30] introduce a perturbance into one of the point sets to move the point clouds out of a small local minima. They show a TRE of 1.17 mm on a femur phantom without any failed trials; however, the test was only done on phantoms 19 Chapter 1. Introduction and Background without soft tissue and required manual segmentation of the bone surfaces to ensure adequate point-to-point correspondence. Carrat et al. [31] and Tonetti et al. [32] used ICP to register intra-operative US to CT to guide the placement of iliosacral screws during pelvic surgery and decrease the X-ray exposure time from 1.07 minutes to 0.35 minutes with no errors in screw placement. However their method requires manual US segmentation and increases the operating time by over 25 minutes [31]. Amin et al. improved ICP by extracting likely bone surfaces from the US volume and weighting the points by the boundary of the bone shadow, the bone surface reflection and the initial pixel proximity of the US and CT [33], however there method was limited to 2D. Barratt et al. take a different approach to their registration algorithm and eliminate the need for pre-operative CT volumes. They create a sta- tistical model of a female pelvis using the CT data from 10 patients and a model of a non-gender specific pelvis using 16 patients, both male and female. They adopt a statistical deformation model as described by Rueck- ert et al. [34] to construct a pelvic model by deforming training sets to a template. Rueckert et al. use non-rigid registration to create the statistical model, yet they rigidly register intra-operative free-hand 3D US by mini- mizing the root-mean-square (RMS) distance between the pelvis model and a manually extracted point cloud from the US bone surface, a method com- monly known as ICP [29]. Barratt et al. report an RMS surface error of 3.5 mm, which is too great for many orthopaedic applications [2]. In an older paper, Barrat et al. also propose a self-calibrating system for registering pre-operative CT to US using intra-operative infrared light emitting diode fiducials and register these to contrast-filled markers attached to the bone during the CT. This method was tested on three cadaver pelves and six cadaver femurs and required 90 minutes to calibrate the cadaver prior to surgery. The registration worked well with femur volumes, resulting in a mean TRE of 1.3 mm and no failed registrations (failure defined as TRE ≥ 5.0 mm). The pelves, however, had an average TRE of 2.2 mm and failed in sixty-three of 735 tests [35]. Similar to Ghavanti et al. above, Foroughi et al. created a statistical 20 Chapter 1. Introduction and Background shape model from the pelves of healthy patients and registered the model to US volumes using ICP to find the model’s real-world coordinates. They reported results of approximately 1.6 mm for dry bone experiments [36]. The Unscented Kalman filter was introduced by Moghari et al. [37, 38] to register US volumes of the scaphoid to CT data. This point-based method shows a significant improvement over ICP in terms of accuracy and robustness. Moghari et al. report 100% convergence when 20% of data points are outliers, compared to 13% convergence when the same test is run using ICP. This method, however, requires manual point cloud selection from the bone surfaces and can potentially have a runtime of over 20 seconds in Matlab, regardless of the number of points, due to the random nature of the algorithm. These point-based registration methods, with the exception of the Un- scented Kalman Filter method, have a very short runtime that is reasonable for intra-operative use in CAOS systems. These same methods; however, suffer from poor accuracy and failure to converge. The Unscented Kalman Filter registration method, on the other hand, provides reliably accurate re- sults but has a runtime which makes it unusable for our proposed method. We propose a registration algorithm which will address these perfor- mance gaps. A summary of the papers described in this section can be seen in Table 1.5.2. Intensity-based methods are, in general, very accurate and robust; however, their runtime makes them impractical for clinical use. Point-based registration methods are generally fast, but are susceptible to falling into local minima. In the next chapter, we propose an algorithm which increases the accuracy and robustness of a fast, point-based method. 1.6 Thesis Organization This thesis is separated into three parts. The first part, Chapter 2, describes the image pre-processing required in our propose algorithm, i.e., the bone lo- calization and point cloud creation and simplification. Chapter 3 describes the implementation of the registration algorithm and how it will be inte- grated into pelvic fracture surgeries. We outline the system requirements 21 Chapter 1. Introduction and Background Table 1.1: Summary of state-of-the-art orthopaedic registration methods. Author Ref. Summary Category Limitations Anatomy Brendel and Winter et al. [11– 15] Maximizing the binary bone surface from the CT image voxel-wise multiplication of the overlaid US volume. Intesity Susceptable to errors due to blurred US sur- face. Vertebra Barratt et al. [29] Create a statistical model of the pelvis from CT scans . The registration is accomplished with ICP. Point Requires manual seg- mentation of the US volume Pelvis Foroughi et al. [36] Statistical shape model registered using ICP Point Error was reported as 1.6 mm for dry bone experiments Pelvis Barratt et al. [35] Register using fiducial implanted in patient during surgery setup. Point Calibration of fiducials requires approximately 90 minutes to setup. Femur Penney et al. [30] ICP registration with introduced peturbance to move point clouds out of local minima Point While this is an im- provement over ICP, it still not sufficiently ro- bust. Femur Carrat et al. [31] ICP Point US bone surface is ex- tracted using manual segmentation Pelvis Tonetti et al. [32] ICP Point US bone surface is ex- tracted using manual segmentation Pelvis Amin et al. [33] Register point clouds weighted by the US intensity value using ICP. Point This method was lim- ited to 2D US and re- quired 30 US images to complete registration, taking approximately 5 minutes to obtain Pelvis Moghari et al. add in more [37] Unscented Kalman Fil- ter registration Point Results are accurate and robust, however the run time is variable due to the random na- ture of the algorithm. Pelvis and Scaphoid Penney et al. [17] Intermediate probabil- ity volumes are created from both the US and CT data sets and reg- istered using NCC. Intensity This method requires many tests data sets to create the intermediate volumes to create prob- ability volumes. Femur and pelvis Wien et al., Sham et al., Kutter et al. and Gill et al. [19, 20, 23– 25, 27, 28] US volumes are simu- lated from the CT data sets and registered us- ing an intensity based method. Intesity While the US sim- ulation has been optimized to run on the GPU, the reg- istration step still approximately 20 min- utes, as it is intensity based. Registration took 90 seconds. Liver and ver- tebra * Ghanavati et al. [18] Statistical atlas of male pelvis from CT volumes and created a simulated US vol- ume. Registered using Linear Combination of Linear Correlations Intensity This method used 110 CT volumes to create the statistical atlas and the reported error was over 7 mm Pelvis 22 Chapter 1. Introduction and Background and validate our registration algorithm and test an US stitching adapta- tion of our propose method in Chapter 4. Our validation is done on both phantom and clinical datasets. 23 Chapter 2 Bone Localization from US and CT In this thesis, we propose to fuse an US volume with a CT volume using a point based registration method. We will show this algorithm to be fast, accurate and robust; however, for this method to work, we must first rep- resent both the US and the CT volumes as a set of points. Our proposed algorithm does not require point correspondence, i.e. the points need not be in the same anatomical location. As our method is targeted for orthopaedic surgery, we can take advantage of the distinct bone surfaces in both US and CT volumes and register the volumes rigidly. US deficiencies make this a very challenging task. In US volumes, only the hard tissue surface nearest the probe is visible. This surface is blurred and noisy and appears as a wide ridge, unlike CT, where the bone surface is generally clean and is often represented as a step. Therefore, we must find a way to to represent both volumes. In this section, we propose a method to create two point sets from the US and CT volumes by localizing the bone surfaces and representing the data with a small point cloud. Figure 2.1 is a pipeline describing the relationship between the surface extraction, point cloud simplification and registration. Figure 2.2 shows the image pre-processing pipeline and Figure 2.3 shows the point cloud simplification pipeline. Simplification is necessary to reduce the registration runtime, as this property is primarily dependent on the size of the point clouds. Simplifica- tion provides a trade-off between accuracy and speed, and will be further discussed in Chapter 4. The runtime of the registration algorithm, which 24 Chapter 2. Bone Localization from US and CT is described in the next chapter, decreases as the number of points in each point cloud decreases, as indicated by the complexity of the registration algorithm O(nm), where n and m are the number of points in each point cloud. If we were to assign a point for each voxel in the bone surfaces that have been extracted from both the US and CT volumes, the resulting point cloud would contain in the order of twenty thousand points. Using the en- tire point cloud would result in the registration lasting minutes instead of seconds. Therefore, we must decrease the number of points in each point cloud while keeping the salient information of the surface. 2.1 Feature Extraction and Point Cloud Creation 2.1.1 Local Image Phase Features Hard tissues appear as highly reflective surfaces in US and therefore reliably appear in US volumes as a bright,ridge-like structure. Therefore, we propose to use the bone surfaces in both the US and CT to register the volumes. The challenge remains to extract the bone surface from the soft tissue features in US. We propose to use a method developed by Hacihaliloglu et al. to automatically extract the bone surface using local phase features [39]. This section describes PSym, which is a method of detecting symmetry about an axis using quadrature filters, and how it is applied to US volumes. Features can be extracted from quadrature filters and be used to find the local phase features. The Log-Gabor filter was chosen as it provides high frequency specificity due to its arbitrary band-width. While other methods of extracting local phase information have been suggested, such as the monogenic signal [40], the Log-Gabor was proven to be best suited for US volumes [41] because of its steerable orientation. Rajpoot et al. report a runtime improvement of 7x [40], however the results are of poorer quality and are insufficient to localize the bone surfaces for registration. The Log- Gabor filter for a 2D signal can be defined as follows: LG(ω) = e − [ (log(ω/ωo) 2 2(log(κ/ωo)) + (φ−φo)2 2σφ ] (2.1) 25 Chapter 2. Bone Localization from US and CT Figure 2.1: Overview of the proposed algorithm’s pipeline. The US Bone Segmentation block is described in Section 2.1.1. The CT Bone Segmen- tation block is described in Section 2.1.3. The Point Cloud Simplification blocks are described in Section 2.2 and the Reinforcement blocks in Section 2.3. The Registration block is described in Chapter 3. 26 Chapter 2. Bone Localization from US and CT Figure 2.2: The CT volume is cropped to the ROI and thresholded. The bone surface is extracted using ray-tracing, see Section 2.1.3. The bone in the US volume is localized using PSym and the final bone surface is determined using ray-tracing, Section 2.1.1 . 27 Chapter 2. Bone Localization from US and CT Figure 2.3: Once the bone surface is extracted from both the US and CT data, the points on the surface are converted to a list and simplified using PSim, see Section 2.2. The resulting point cloud is reinforced with high curvature features, Section 2.3. 28 Chapter 2. Bone Localization from US and CT In Equation 2.1, above, ωo is the center frequency of the Log-Gabor filter and κ defines the bandwidth. σφ can be used to define the angu- lar bandwidth, ∆Ω = 2σφ √ 2log2. However, as shown by Hacihaliloglu et al. in [39, 42], 2D methods ignore correlations between adjacent slices and are therefore subject to large spatial compounding errors as well as errors associated with the US beam thickness effects. LG(ω) = e − [ (log(ω/ωo) 2 2(log(κ/ωo)) + α(φi,θi) 2 2σ2α ] (2.2) The 3D extension, above, includes the term α(φi, θ) 2 which is the filter’s angle between the azimuth (φ) and elevation (θ). For both the 2D and 3D cases, the Log-Gabor filters are built into banks with individual components separated by a scaling factor of δ, based on a minimum wavelength λmin and multiple m. Therefore, ωo can be determined for each filter component as: ωo = 2δm−1 λmin (2.3) These properties of the Log-Gabor filter are important, as they allow us to selectively steer the filters to obtain directional information about the bone surfaces. These properties can be seen in Figure 2.4. The effect of the filter orientation can be see in Figure 2.5. In regions of high symmetry, i.e. around ridge-like objects, the even com- ponents of the local phase features greatly outweigh the odd components. This can be used to find bone surfaces in US images, as these surfaces ap- pear as ridges across the surface. The local phase features extracted with the Log-Gabor filters, above, can be used to calculate the even and odd components of the image and therefore calculate the PSym seen in Figure 2.4 [43, 44] and similarily in 3D, [39]: PS(x,y, z) = ∑ r ∑ mb[|erm(x,y, z)| − |orm(x,y, z)|]− T c∑ r ∑ m √ erm(x,y, z)2 + orm(x,y, z)2 + ε (2.4) erm(x, y, z) and orm(x, y, z) are the responses to the even and odd filters 29 Chapter 2. Bone Localization from US and CT Figure 2.4: Example of a 2D Log-Gabor radial and angular filter compo- nents. (a) and (b) show the radial component at two different scales, as described by Equation 2.3. The angular components of the Log-Gabor filter are seen in (c) - (h) at six different orientations. 30 Chapter 2. Bone Localization from US and CT Figure 2.5: Three filter angles are demonstrated below to capture three different bone surfaces in the US volume. The first row is a slice from the B-Mode US volume. The second row is a slice from three different angular components of the filter. The bottom row is the angular filter response for each corresponding filter. 31 Chapter 2. Bone Localization from US and CT Figure 2.6: An example of the US surface extraction steps used in our proposed algorithm. The top row is a single 2D slice taken from the volume in the bottom row. The B-Mode US can be seen in (a) and (d). The PSym of the B-Mode volume is seen in (b) and (e). Finally, the outline can be seen in (c) and (f). over the number of scales m and orientations r. The threshold component, Tr, is added to PSym to account for noise in the US image. Hacihaliloglu et al. [44] determined that Tr can remain constant for all US images without degrading the results in both phantom and in-vivo images. To avoid divi- sion by zero, a very small constant, ε (typically 10−4), is included in the denominator of 2.4. The outline of the bone is calculated as the maximum PSym along the direction of the probe; an example can be seen in Figure 2.6. 2.1.2 Automatically Optimized Parameters for Log-Gabor Filter The Log-Gabor filter has been shown to be very effective for localizing bone surfaces in US images, however one of the challenges remains to find au- 32 Chapter 2. Bone Localization from US and CT tomatically optimized parameters for all cases. We will show in Chapter 4 that empirically set parameters are sufficient for phantom data that lack soft- tissue interfaces; however, these parameters need to be optimized for each image that contains soft-tissue. Hacihaliloglu et al. developed a method to automatically obtain the parameters for the Log-Gabor filter and we outline their methods below. A more detailed explanation is found in [45]. The Log-Gabor’s radial bandwidth, determined by κ/ωo, is related to the speckle size in the US volume. Using autocorrelation as a measure of speckle size, we can calculate this ratio as seen in Equation 2.5, where F is the full- width at half-maximum of the autocorrelation and r is the pixel resolution [45]: κ ωo = e−.25 √ 2ln(2)F×r (2.5) The Log-Gabor filters are very sensitive to speckle noise for small scales, λ. The minimum filter scale, λmin is obtained by determining the “rid- geness” of the B-mode image by applying a line-enhancing filter based on eigenvalue analysis of the Hessian matrix, because the bone surfaces gener- ally appear as ridge-like structures in B-Mode US images. The orientation of the filters must be chosen to best capture the bone surface. Figure 2.7 shows an example of a filter bank that captures all regions of the bone surface by steering the angular filter. The radon transform is used to determine the filter angles, where the peaks correspond to the mean value of the projection angles. These initial angles are used to determine the filter scale and are re-calculated for the ridge strength image obtained above. The filter’s angular bandwidth, σφ, corresponds to the standard devi- ation of the angular spreading function. This parameter is important, as it determines the angular filter’s selectivity. σφ is determined to be the peak kurtosis value of the radon transform of the ridge strength image over multiple angular bandwidth values because a greater kurtosis indicates the variance is due to infrequent high-intensity voxel values. This indicates that only the strongest response, i.e., bone surface, is found. The kurtosis is 33 Chapter 2. Bone Localization from US and CT Figure 2.7: The filter’s orientation must be steered to cross different seg- ments of the bone surface to capture the entire bone surface. The red arrows in each B-Mode image show the surface being captured with the adjacent angular filter. 34 Chapter 2. Bone Localization from US and CT calculated from the radon transform of multiple ridge strength images; the angular bandwidth is chosen to correspond to the largest kurtosis. 2.1.3 Bone Surface Extraction from CT Volumes In contrast to US, extracting the bone surface from the CT volumes is comparatively easy. CT volumes of bone tissue can be reliably thresholded, in particular the pelvis can be consistently chosen to be between 110-320 HU, (200 HU is commonly accepted as a realistic value of segmenting the pelvic bone), according to Sugano et al [46]. Our system assumes the US probe is tracked and the pelvis is anatomically digitized prior to the start of the surgery, therefore, the position of the US probe is known with respect to the CT dataset. In our application, the location of the probe is known with respect to the CT volume. Therefore we know the location of the US probe with respect the the pelvis and can estimate the ROI to crop the CT and the direction of the probe to extract the visible bone surface in the thresholded CT. Given this information, we use ray-tracing to extract the bone surface from the thresholded CT volume visible in the US volume. This results in a binary volume with a one-voxel thick outline, similar to the US volume extracted from the PSym volume. Both of these bone surfaces are comprised of approximately twenty thou- sand voxels. The next section details how we extract a small point cloud from this surface to be used in the registration algorithm proposed in Chap- ter 3. For the remainder of this section, we will refer only to a point cloud P = {pi ∈ R} for n points, p, and not distinguish between US and CT data, as the steps taken to simplify each bone surface is the same. 2.2 Point Cloud Simplification The outlines of the bone surfaces, extracted in the previous section, must be simplified before being registered. This section describes PSim and how we have optimized the algorithm using a hierarchical tree to estimate the force calculations. The user can input the number of points in the algorithm and 35 Chapter 2. Bone Localization from US and CT we explore the effect of point cloud size on the runtime of the registration algorithm in Section 4.2.2. PSim is a method of simplifying the point cloud that can be analogized as finding the equilibrium of point forces that are restricted to a 3D surface. 2.2.1 Point Cloud Simplification Using Particle Simulation One of the simplest and most time-effective methods of reducing the number of points in a large point cloud is simple down-sampling. While the results were shown to be robust with ten out of ten tests converging, [47], there is a possibility that the registration will fail if a large surface area lies in one plane, which could be omitted in this näıve down-sampling method. Therefore, another method must be found to simplify P . We propose a method, called PSim, to decrease the number of points in P , while retaining the salient features needed to register the US and CT point clouds [48]. This method redistributes a random selection of point to better represent a surface, analogous to the diffusion of particles in a volume. This method selects a random sample of n points from the original point cloud, PS , and iteratively relaxes the new subsampled point cloud, Ps, by calculating the forces exerted on pi by every other point in Ps and moving each point along force Fi. Finally, pi must be projected back onto the bone surface. Ps is initialized by randomly selecting n points from PS , allowing the user to specify the size of the simplified point cloud instead of the size being dependent on characteristics of PS and the US volume. PS and the initialized Ps can be seen in figures 2.8 (a) and (b). Once Ps is initialized, the forces must be calculated for each point in the point cloud. We have implemented the force, Fs, to have an internal and an external force component, FIn and FEx respectively, Equation 2.6. If two points are forced to overlap, i.e. share the same coordinates, the PSim cannot recover and separate the points; therefore a small random jitter force, 36 Chapter 2. Bone Localization from US and CT Fj was added to Fs. Fs = FIn + FEx + Fj (2.6) FIn is a linear force calculation based on the distance between two points, r and a force constant, k. Pauly et al. suggest a constant value of k to be acceptable for most cases [48]. To simplify the calculations of FIn, only points within r voxels contribute to FIn. Equation 2.7 shows FIn. Figures 2.8 (c) and (d) show FIn superimposed on Ps. FIn = n∑ m=1 k(r − ‖p− pi‖) · ( p− pi ‖p− pi‖ ) (2.7) PSim has been demonstrated on point clouds which encompass a con- tinuous surface [48], such as a sphere or torus, however it is not ideal for a surface with discrete edges, as is the case with the bone surfaces extracted in 2.1.1. This is caused by FIn forcing points towards low density areas, i.e. the edge of the surface. To prevent this edge crowding effect, we have added FEx, which acts to push points away from the edge. As Fj is only necessary to prevent overlapping points when FIn and FEx are non-zero, Fs is zero when FIn is zero. FEx was implemented as a trough function; a force was imposed to push the points away from the edge of the US surface. An exhaustive calculation to find FIn is used by Pauly et al. [48] and has an undesirable complexity of O(n2). Therefore, we propose to optimize the force calculation using a hierarchical octree, initially proposed by Barnes et al. [49], that is typically used for N-Body simulations in astro-physics and modelling computer networks [50], which will be be described in Section 2.2.2. After the force has been calculated, each point is moved according to Fs. This step is called the relaxation of the point cloud. Finally, the relaxed points must be projected back onto the surface, Figure 2.8 (e). As suggested by Pauly et al., the point is projected onto the nearest point in PS . This technique is sufficient because of the relatively high particle density of PS . These steps are iteratively repeated until an equilibrium is reached, Figure 37 Chapter 2. Bone Localization from US and CT 2.8 (f). This algorithm is outlined in Algorithm 1. Algorithm 1 Particle Simulation Algorithm Input: Desired size of final point cloud, n Point cloud, PS , size > n Output: Simplified point cloud, Ps, of size n 1: procedure PSim(PS , n) 2: Randomly select n points from PS 3: repeat 4: Calculate force, Fs, for each point in Ps 5: Relax each point along Fs 6: Project points onto surface of PS 7: until equilibrium reached 8: end procedure 2.2.2 Optimization of Force Calculations Pauly et al. use an exhaustive search and only include significant forces, determined by r, which in their implementation of PSim has complexity O(n2). As PSim requires the recalculation of the force calculation for each iteration, we propose to optimize this step. Barnes et al. proposed an algorithm in 1986 (now called the Barnes-Hut algorithm) which reduces the force calculations to O(nlogn) [49]. This algorithm is widely used in astro- physics for the gravitational N -Body simulation problem [50–53]; however, we have adapted it to calculate repulsive forces for points in Ps. The Barnes-Hut algorithm creates a hierarchical tree to represent the point cloud volume using octnodes, i.e. a tree where each node has eight children, each representing equal volumes. The root node represents the entire volume containing Ps, each child node in the first layer represents one eighth of the volume, hence the term octnode. The tree is started by inserting the first point into the root node. One-by-one, other points from the point cloud Ps are added to the tree. The position of a point in the tree is determined by its spacial coordinate in the volume and any prior insertions or points that are within a neighbourhood of the current point. If the node is empty, a point, pi from Ps, is assigned to that node. Whenever the point 38 Chapter 2. Bone Localization from US and CT Figure 2.8: PSim is used to decrease the size of PS without losing salient information. (a) shows the original dense point cloud, PS . (b) shows the randomly sampled point cloud. FIn is shown in (c) and a close up is shown in (d). (e) shows the relaxed points projected back on the surface described by PS . The final point cloud, Ps, is shown in (f) after three iterations. 39 Chapter 2. Bone Localization from US and CT encounters a node with another point already assigned, the node is branched into eight nodes (i.e. the octant is subdivided into smaller octants), and both points are moved down to the appropriate children nodes. This branching continues until the two points are assigned to distinct nodes, at which point pi has been placed in the tree. The tree is complete once every point in Ps has been inserted into the tree with at most one point assigned to each node. Once the tree is complete, Fs,i is calculated for each point. The tree is traversed, starting at the tree node. If the visited node, nj , is an external (leaf) node, the force of pj acting on pi is calculated and added to FIn. If nj is internal, the size-distance ratio q = sd is calculated, where s is the width of nj and d is the distance between pi and the center-of-mass of nj . If q is less than some constant, θ, we can assume that the particles in this branch are sufficiently far away from the point of interest and can be modeled as a single point. The force is calculated using the center-of-mass and total mass of the node. This avoids individual force calculations for points in children nodes of nj . If q > θ, however, the algorithm proceeds to nj ’s children nodes and recursively repeats this test, as the points are too close to the point of interest to be modelled as a single point. This algorithm is described in figures 2.9 to 2.12, using a quadtree as an example. The 2D Barnes-Hut tree is built by adding points one-by-one into a quadtree. As demonstrated in Figure 2.9 (a), the first point, A, is placed into the root of the tree and the center-of-mass and total mass for that node is calculated. In Figure 2.9 (b), Point B is inserted into the root node, but because there are now two points in the root, both points are inserted into the appropriate node, one level down as in Figure 2.9 (c). Points A and B are still in the same quadrant, therefore we repeat the previous step, but this time, A is inserted into the south-east child and B is inserted into the north-east child, Figure 2.10 (a). Now A and B are in different nodes, C is inserted into the appropriate child of the root, Figure 2.10 (b); the south east quadrant. C is the only point in this quadrant, so no further action must be made. For each step, the total mass and the center-of-mass are updated for each node visited. Note that this example is completed using a quadtree to 40 Chapter 2. Bone Localization from US and CT demonstrate the 2D case. Our implementation of the Barnes-Hut algorithm is in 3D, using octrees where each node has eight children instead of four. The building of the Barnes-Hut tree is described in Algorithm 2. Once the Barnes-Hut tree has been built, the point forces can be esti- mated by assuming the effect of several points as one, providing they are suf- ficiently far away. Continuing the example above, Figure 2.11 demonstrates how the Barnes-Hut algorithm estimates Fs,C . First, we start traversing the Barnes-Hut tree at the root, by determining if q = sd is less than θ. In this example, we use a θ equal to 1.5. In Figure 2.11 (a) q = 100/40 and is greater than θ, therefore we drop to the next level, Figure 2.11 (b). Repeating our previous step, we see that q < θ, therefore we can treat this node as an external node and the masses of A and B are treated as one. The force of this node acting on C is calculated, and added to Fs,C . There are no more points, therefore Fs,C is complete. The example is continued in Figure 2.12, which demonstrating the cal- culation of Fs,B. This force can be calculated starting with the root node, see Figure 2.12 (a), resulting in q > θ. Therefore, the root is treated as an internal node and we inspect its children. C is an external node, therefore we calculate the force of C acting on B and add it to Fs,B Figure 2.12 (b). This is also the case for A, therefore we add the force to Fs,B as well, see Figure 2.12 (c). Fs,A is calculated similarly to Fs,B. The force calculation can be seen in Algorithm 3. 2.3 Point Cloud Reinforcement Using High Curvature Features Our particular anatomical region of interest is the human pelvis, but US can only image the iliac spine and the pubis (Figure 2.13 (a)) because other regions of the pelvis are occluded by many layers of soft tissue. Simplify- ing the point cloud significantly reduces the runtime, but in doing so, we could potentially lose key salient points which contribute to an accurate registration result. For example, if we consider our particular region of in- 41 Chapter 2. Bone Localization from US and CT Figure 2.9: An example of the algorithm used to build a 2D Barnes-Hut tree. 42 Chapter 2. Bone Localization from US and CT Figure 2.10: Continuation of example in Figure 2.9. 43 Chapter 2. Bone Localization from US and CT Figure 2.11: 2D example calculating Fs,C for θ = 1.5. 44 Chapter 2. Bone Localization from US and CT Figure 2.12: 2D example Calculating Fs,B for θ = 1.5. 45 Chapter 2. Bone Localization from US and CT Algorithm 2 Barnes-Hut Tree Building Algorithm Input: Point cloud, Ps Output: Hierarchical Barnes-Hut tree, BH procedure BHtree(PS , FIn) 2: for n points do Insert point, pi, into root node repeat 4: if Node does not contain a point then Leave pi in node and update centre-of-mass and total mass 6: else if Node is internal then Update centre-of-mass and total mass of node. Recursively insert point pi into child octants 8: else if Node is external and contains point pk then while necessary do 10: Create children octnodes and insert pi and pk in ap- propriate children nodes end while 12: end if until the point is in an empty external node 14: end for end procedure 46 Chapter 2. Bone Localization from US and CT Algorithm 3 Barnes-Hut Tree Force Algorithm Input: Barnes-Hut tree, θ, Output: FIn procedure BH(BH) Start traversing tree at root node 3: for each point, pi do if external node then Calculate the force of node exerted on point pi 6: else if internal node then Calculate the distance between pi and the center-of-mass of the node, d Determine the size of the octant, s 9: Calculate q = s/d if q < θ then Treat internal node as single point force on pi 12: else Traverse node’s children end if 15: end if end for end procedure 47 Chapter 2. Bone Localization from US and CT Figure 2.13: Pelvic bone anatomy and the effect of high curvature reinforce- ment - We are primarily interested in imaging the iliac spine and the iliac crest (a). While the simplified point cloud can register with respect to the pelvic table (indicated by the outlined arrow in (b)), high curvature features are needed to register correctly along to the iliac crest (indicated by the solid arrow in (b)). terest, registering the volumes using the simplified point clouds calculated in Section 2.2 results is adequate registration with respect to the pelvic table (outlined arrow in Figure 2.13 (b)). The 3D US volumes tend to slide down the iliac crest, resulting in a failed registration. Therefore, we propose to take advantage of the high curvature regions around the iliac crest to correct this misregistration. We propose to measure the high curvature regions of the bone by using an intrinsic measurement of curvature, Gaussian curvature (K). Gaussian curvature, at any point on a parametrized surface, is the product of the principle curvatures κ1 and κ2 [54], where e through G are defined in the following paragraph: K = κ1κ2 = eg − f2 EG− F 2 (2.8) Our method calculates K from the coefficients of the first fundamental form for the surface s(u, v) [54]: ds2 = Edu2 + 2Fdudv +Gdv2 E = |∂x ∂u |2, F = ∂x ∂u ∂x ∂v ,G = |∂x ∂v |2 (2.9) 48 Chapter 2. Bone Localization from US and CT e, f and g are the coefficients of the second fundamental form: e = det(xuuxuxv)√ EG− F 2 , f = det(xuvxuxv)√ EG− F 2 , g = det(xvvxuxv)√ EG− F 2 (2.10) Note that special care must be taken around the edges of the parametrized surface of the US volume, because our choice of curvature features would detect the edge surfaces as high curvature points, when in fact, they do not represent an anatomical feature. Therefore, we do not use any points on the edge of the surface to reinforce the final point cloud. Even if the edge happens to lie on a true feature, it is not possible to distinguish from the US volume alone. Figure 2.14 demonstrates the point cloud reinforcement using high curvature features. The points with high curvature features determined in this section are combined with the output of the PSim. In the next chapter, we will use this combined point cloud to register the US and CT volumes using the point sets extracted in this chapter. 49 Chapter 2. Bone Localization from US and CT Figure 2.14: Gaussian curvature reinforcement: (a) shows a 3D render- ing of an US volume acquired from the iliac spine and (b) corresponding parametrized surface. The dark background is not part of the surface and is therefore not included in the K calculation. (c) shows the high curvature features, K, on the parametrized surface. (d) is the final simplified point cloud (blue) with reinforced curvature features (red). 50 Chapter 3 Fast 3D Correspondence of US and CT Volumes and Its Application in Pelvis Fracture Surgery As discussed in Chapter 1, fusing US with another modality is very difficult, even in orthopaedics, where bones appear as bright, rigid features. This is due to the poor signal-to-noise ratio (SNR), limited field of view and occlusion of deep tissues. Many methods have been proposed to fuse US and CT and they generally fall under two categories: point- and intensity- based registration. Point-based methods tend to be fast; however it can be challenging to represent the US and CT volumes with an appropriate set of points and standard methods tend to fall into local minima, making these methods unreliable. Intensity based methods have been shown to be accurate and robust; however, they are time intensive. For example, mutual information is a common intensity-based similarity metric used for multi-modal registration; however even optimized to be run on a GPU, registration take eleven seconds [28]. In Chapter 2 we concentrated our efforts on representing the bone surfaces and a simplified bone surface a simplified point cloud which was reinforced with high curvature features. In this Chapter, we will propose a fast point-based registration method that is robust to noise and variations in the point clouds representing our US and CT volumes. This method does not require a point-to-point correspondence; instead it matches a statistical model of each point cloud by minimizing the 51 Chapter 3. Fast 3D Correspondence of US and CT Volumes and Its Application in Pelvis Fracture Surgery distance between them. 3.1 Registration of US to CT Volumes Using Gaussian Mixture Model Matching Most CAOS registration research in the past five years has focused on intensity-based methods, as previous research had shown point-based regis- tration methods to be less robust. Any point-based method that requires an identification of corresponding points is likely to fail in US as finding a pair of features on the iliac spine in both the US and CT is very challenging due to noise and the inherent different natures of the modalities. We could potentially gain an advantage if we model the point cloud as a whole in- stead of requiring point-to-point correspondence. Therefore, we propose to represent the bone surface as a GMM and minimize the difference between the model from the two modalities. A method of registration using GMMs was proposed by Jian et al. in 2005 [55] and further described by Jian et al. in 2011 [56], but (to our knowledge) has not previously been applied to medical images. We propose to use GMM registration, as it faster than intensity-based methods and significantly more robust than other point- based methods. While GMM registration can be adapted for deformable registration, rigid registration is best suited to register the CT and US bone surfaces. 3.1.1 The Application of Gaussian Mixture Models to Point Clouds A GMM is a set of continuous density functions that statistically models data as a group of multi-dimensional Gaussian distributions. Each Gaussian distribution is called a component and is defined for the point set x ∈ Rd: φ(x|µi,Σi) = exp[−12(x− µi)TΣ−1i ](x− µi)√ (2pid)|det(|Σi)| (3.1) In Equation 3.1 above, d is the number of of dimensions, µ is the mean 52 Chapter 3. Fast 3D Correspondence of US and CT Volumes and Its Application in Pelvis Fracture Surgery Figure 3.1: The proposed registration pipeline using GMM. The inputs of the registration algorithm is the simplified point clouds describing the US and CT bone surfaces as described in Chapter 2. The output is the transformed US volume. 53 Chapter 3. Fast 3D Correspondence of US and CT Volumes and Its Application in Pelvis Fracture Surgery Figure 3.2: Example Illustrating GMMs - An example of a randomly gen- erated 2D GMM comprised of three components with mean vectors, µi, at [−3,−5], [1, 2] and [5,−4]. of the Gaussian and Σ is the covariance. Out data is 3D, therefore our Gaussian distributions are defined using d = 3. Each GMM, p(x), has k components which are weighted by w. Figure 3.2 is 2D example of a GMM with three components. p(x) = k∑ i=1 wiφ(x|µi,Σi) (3.2) To simplify the GMM construction, Jian et al. define the number of components to be the number of points in the point cloud, k = n and the mean vector, µ, to be the spatial location at each point. This allows the GMMs to be constructed using a Gaussian kernel and reduces the number of parameters needed to implement the registration. More details on the implementation of this algorithm can be found in [56]. 54 Chapter 3. Fast 3D Correspondence of US and CT Volumes and Its Application in Pelvis Fracture Surgery 3.1.2 Minimizing the Difference between Two Guassian Mixture Models The two point sets are registered by iteratively minimizing the distance between them. Jian et al. chose the L2 distance as the similarity metric as there exists a closed-form expression between the GMMs that is related to the estimator L2E, allowing for a registration algorithm that is inherently robust due to the statistical modelling of the point clouds [56]. The optimization algorithm minimizes the dL2 cost function for point sets modelled by g(x) and f(x): dL2(g, f,R, t) = ∫ (g − fR,t)2dx = ∫ (g2 − 2gfR,t + f2R,t)dx (3.3) The transformed model above is written as: fR,t(x) = m∑ i=1 wiφ(x|Rµi + t,R(Σi)RT ) (3.4) Taking into account Equation 3.4, Equation 3.3 can be reduced in the rigid case, because ∫ p(T (M, θ))2 = ∫ p(M), for the transformed model M rotated by θ, and ∫ g2 does not depend on any transformation. Therefore, only the cross-term fR,tg need be considered during optimization. The cross-term becomes the summation of n Gaussian functions at m points, or a discrete Gauss transform: gfR,t = m∑ i=1 n∑ j=1 w(f,i)w(g,j)φ(x|Rµ(f,i) +t−µ(g,j),RΣ(f,i)RT +Σ(g,j)) (3.5) Equation 3.5 is iteratively minimized until convergence or a certain num- ber of iterations have passed, specified by the user. Jian et al. set the maxi- mum number of iterations to eighty [56]. In our experiments, the algorithm is iterated approximately four times and never failed to converge. To our knowledge, this method has not been applied to 3D data apart 55 Chapter 3. Fast 3D Correspondence of US and CT Volumes and Its Application in Pelvis Fracture Surgery from initial validation completed on standard data sets by Jian et al. [56] or to any medical applications. In Chapter 4, we will demonstrate and validate this registration algo- rithm on both phantom and clinical data. 3.1.3 The Effect of Point Cloud Size on Registration Runtime. In Chapter 2 we discussed the localization of the bone surfaces and the point cloud creation including the simplification process used to select a subset of points; however, here we will discuss the size of the point cloud and how it affects the registration. Using a näıve approach to simplifying the point clouds, i.e. simple down- sampling, at least 800-1000 points are needed to represent the bone surfaces [47], as the surfaces are not well distributed and the registration tends to fall into local minima; even this can be insufficient if large surface areas are omitted. Introducing PSim significantly decreases the number of points needed to register the point clouds, which we will demonstrate in Chapter 4 using a phantom pelvis. The runtime grows in O(nm) when directly computing the Gauss transform. While the fast Gauss transform can be used to decrease the runtime to O(n+m) complexity, it is not effective for fewer than 1000 points due to computation overhead (Jian [56]). For the case of US stitching (see Section 3.2 below), even fewer points are required to accurately register the US volumes because of the inherent similarity of the point clouds, resolution and quality of bone surface extrac- tion. 3.2 Extension of Methods to US-US Stitching We have chosen to apply this registration algorithm to fuse US to CT vol- umes of the iliac crest; however it can easily be applied to other applications. Stitching is a specific type of registration which fuses datasets with differ- ent spatial information with varying overlap, instead of fusing two data sets 56 Chapter 3. Fast 3D Correspondence of US and CT Volumes and Its Application in Pelvis Fracture Surgery with different temporal or modal information (as previously discussed in this thesis). We propose to extend our method to stitch two US volumes; we need only change the image pre-processing bone localization step outlined in Chapter 2 to repeat the PSym calculation on the second US volume, seen in Figure 3.3. The main difference observed is the pre-processing runtime due to the necessity of calculating the PSym for each US volume as well as the initial alignment of the two datasets. 3.3 US to CT Registration in a Pelvic Fracture Surgery Guidance System The purpose of this thesis is to contribute to a larger orthopaedic guid- ance system, and part our work is the description of the system and its requirements. In Chapter 1, we described a standard pelvic fracture sur- gical procedure. In this section we will describe the proposed system and how US is to be incorporated into the surgical procedure. This thesis is fo- cused on the US component of the the system described in this section. The US volumes acquired by the surgeon are registered with the pre-operative CT and are used to update the displayed pelvis model using the algorithm proposed in Chapters 2 and 3. 3.3.1 Introducing US into Pelvic Fracture Surgeries The surgery described in Section 1.4 requires a large number of time-consuming fluoroscopic images obtained with a C-arm operated by a technician and the entire staff must wear lead protective lead shields, which are inconvenient and can lead to back pain. Readjusting the C-arm multiple times to image a single region of interest can take several minutes. We propose to develop a system which will use US imaging to not only minimize exposure to ra- diation, but also minimize the time needed to reduce the pelvis and place implants. We will discuss the potential use of US in pelvis fracture surgery in this section. Our system takes advantage of the pre-operative CT data 57 Chapter 3. Fast 3D Correspondence of US and CT Volumes and Its Application in Pelvis Fracture Surgery Figure 3.3: Flowchart describing the image processing pipeline for US stitch- ing using our proposed algorithm. 58 Chapter 3. Fast 3D Correspondence of US and CT Volumes and Its Application in Pelvis Fracture Surgery taken to diagnose the patient and plan the surgery, therefore we do not re- quire the patient to undergo another harmful and expensive procedure to acquire high SNR images to fuse with US datasets. The main objection of the proposed system is to estimate the bone frag- ment movement prior to and during surgery. We plan to accomplish this by tracking the near side of the pelvis with US datasets and overlaying it on a pre-operative CT volume. It might be necessary to stitch several US volumes together to create a larger ROI to fuse with the CT volume. The proposed pre-operative care is depicted in Figure 3.3.1. This is necessary to find the location of the patient’s pelvis with respect to the CT model. The system will ask the user to use the US to find the pelvis coordinates using anatomic digitization. These steps are described in Figure 3.3.1 and must be completed before the surgery is started. The three regions of the pelvis that are best captured using an US probe are the left and right iliac crests and the pubic symphysis and therefore can be used as anatomical landmarks to digitize the pelvic coordinates as proposed by Spencer et al. [57]. After the digitization, the surgeon will start the first incision and begin the surgery. The surgeon will then reduce the pelvic fragments, updating the system with the US probe after each attempt, replacing the previously used fluoroscopy, Figure 3.3.1. The system will use the US volumes to track the near side of a fragment and then overlay the CT to predict alignment at the fixation site. The system will update the model of the patient’s pelvis on surgeon’s display. Once the surgeon is satisfied with the bone alignment, he/she will position the first K-wire. The K-wire is inserted with a drill, which may be tracked to show the tool’s trajectory. It should be noted; however, fluoroscopy will likely still be needed to see the tool tip which will be occluded by bone or many layers of soft tissue. Once satisfied with the position, the surgeon will drill the K-wire and check the placement again, with fluoroscopy or US, if the conditions are suitable, Figure 3.3.1. When enough K-wires are inserted, the surgeon will inform the system that they are ready to insert a cannulated screw. The surgeon must input the type of screw (length, diameter, manufacturer, threading length, etc.). The 59 Chapter 3. Fast 3D Correspondence of US and CT Volumes and Its Application in Pelvis Fracture Surgery system will then update the display to show the tracked drill’s trajectory and the depth of the screw as it is being inserted. It is important for the surgeon to know the depth of the screw and the location of the tip, to avoid the L5, S1 and S2 nerves and the iliac vessels in the sacral region, Figure 1.8. The gastrointestinal tract, urinary track and the uterus (in female patients) are also located close the bony pelvis and are at risk during pelvis reduction and fixation. Misalignment of the screws can cause nerve or vascular damage and will often result in severe sequelae [3]. After the screws are over-drilled, the K-wires are removed and the work will be confirmed with both fluoroscopy and an US update. The surgeon then decides if they need to go back and further reduce or fix the patient, or if the surgery can be completed, Figure 3.3.1. We do not propose to eliminate fluoroscopy entirely, because of the limi- tations of US. It is important to incorporate such a system into the surgeon’s work flow without introducing unnecessary steps and equipment, if integra- tion is to be successful. We estimate less than twenty fluoroscopic views will be needed for the surgery. Three initial views will be needed to detect the initial bone fragment location and movement. Three views will be taken for each k-wire inserted and three more for each screw. Lastly, the three final views will be needed to asses the completed reduction and fixation before the surgery can be completed. Each fluoroscopic view exposes the patient to approximately 1-3 seconds of radiation. Therefore, we expect our system to reduce the exposure time from 12 minutes to approximately 1 minute [4]. 60 Chapter 3. Fast 3D Correspondence of US and CT Volumes and Its Application in Pelvis Fracture Surgery Assessment CT Move Patient to O.R. Plan Treatment Patient moved CT Pre- processing Patient moved Pre-Operation Is patient a good candidate for CAOS? DICOM CT volume No Patient Admitted End Process Diagnostic X- ray Patient moved Sterilize US Equiptment Is the patient a good candidate for surgery? Yes Use Alternative Method No Load CT Data to CAOS DICOM CT volume Yes Waits Patient moved Legend Patient Surgeon Staff Automatic Software Function Interactive Software Function To OR Preparation Figure 3.4: Flowchart of patient’s pre-operative care. 61 Chapter 3. Fast 3D Correspondence of US and CT Volumes and Its Application in Pelvis Fracture Surgery Legend Patient Surgeon Staff Automatic Software Function Interactive Software Function Acquire US Volumes Landmark Digitization Patient and Equipment is Prepared Preparation (in O.R.) Display Pelvis Model CAOS Software Started Patient is ready Ready for data input Prompt user Attach Optical Trackers to Bone Surface Optically tracked pointer Surgery Starts Anatomic Digitization Software Visualization From Pre-op. To Pelvic Reduction Procedure Figure 3.5: Flowchart of proposed OR preparation. 62 Chapter 3. Fast 3D Correspondence of US and CT Volumes and Its Application in Pelvis Fracture Surgery Pelvic Reduction Fracture Reduction Sufficient? No Repeat as necessary Update surface model Rescan Read Display Software Visualization To K-Wire Insertion Legend Patient Surgeon Staff Automatic Software Function Interactive Software Function From Digitization and Treatment Plan Yes Figure 3.6: Flowchart of pelvic reduction procedure in proposed system. 63 Chapter 3. Fast 3D Correspondence of US and CT Volumes and Its Application in Pelvis Fracture Surgery K-Wire Insertion To Screw Drilling Legend Patient Surgeon Staff Automatic Software Function Interactive Software Function From Pelvic Reduction Position K- wire Ready to drill Well Positioned? No Yes Insert K-wire Trajectory Correct? Images of k-wire positions No Remove K- wire Add more K-wires? Yes Yes Display K- wires using 2D Fluoroscopy Display the Final K-wire Position with 2D Fluoroscopy Display Tracked Drill Display showing K-wire placement No Software Visualization Figure 3.7: Flowchart of K-wire insertion procedure in proposed system. 64 Chapter 3. Fast 3D Correspondence of US and CT Volumes and Its Application in Pelvis Fracture Surgery Screw Drilling Legend Patient Surgeon Staff Automatic Software Function Interactive Software Function Confirm with 2D Fluoro Is Reduction and K-wire placements correct? Yes Finish Surgery Determine if you should go back to reduction or K-wire insertion No Choose Cannulated Screw Input Screw Choice Screw length and type Is the screw self drilling? Drive ScrewYes Over Drill No Display screw placement Confirm Reduction with Display Images of reduction and implant placement Software Visualization From K-Wire Insertion Figure 3.8: Flowchart of screw placement procedure. 65 Chapter 4 Validation of Proposed Registration Algorithm In the previous chapters, we described an US to CT registration algorithm and proposed a surgical system which incorporates US into pelvic fracture surgery. In this chapter we will list the requirements of a registration algo- rithm for the proposed system, validate our registration algorithm according to these requirements and demonstrate the validity of our algorithm choices by quantitatively testing the runtime, accuracy and robustness of our sys- tem. This validation will include a phantom study and a preliminary clinical study. 4.1 Registration Requirements for Proposed Pelvic Fracture Surgical System To ensure our registration algorithm can be introduced into the surgical work flow without impeding the surgeon we impose the following constraints on the US to CT registration algorithm to comply with the system’s require- ments: 1. A CAOS system should have an error no greater than 2 - 4 mm for fracture reduction applications [2]. 2. The surgeon can require up to 20 fluoroscopy images to place one K- wire, and this process can take up to 5 minutes. Therefore our system should not take longer than 5 minutes to perform the same task. We propose to use a scanned 3D linear array which requires approximately 66 Chapter 4. Validation of Proposed Registration Algorithm 5 seconds to acquire a single US volume (this is dependent on the depth). 3. The algorithm must also be robust, as we cannot rely on the surgeon monitoring the system. The surgeon will, however, acquire fluoro- scopic images to verify the alignment of the pelvic fragments. 4. The registration must be fully automatic and should not require the surgeon to segment or manually align either the US or CT and should not require any parameters to be set for the registration algorithm. The purpose of this chapter is to test the suitability of our proposed algorithm according to the requirements stated above. We created a pelvis phantom and test the runtime, accuracy and robustness as well as the effect of introducing PSim and high curvature features. Clinical datasets were also tested to ensure the algorithm meets these requirements for patient data and show the necessity of automatically optimized parameters for bone localization using PSym. 4.2 Experimental Setup and Results for Phantom Study 4.2.1 Phantom Creation and Equipment We created a phantom to test our proposed registration algorithm using a high-resolution radio-opaque Sawbone hemi-pelvis (Sawbone, Pacific Re- search Laboratories Inc., Vashon, Washington), model number 1297-22. We were only interested in validating our method on the iliac crest because this is the most accessible part of the pelvis in a patient [57]. We cut this portion of the pelvis and placed it in an acrylic tube to place in a high-resolution CT. This model does not take specimens longer than 12 cm or wider than 10 cm. The phantom was suspended in a PVC gel with sixty-one 1 mm metal fiducials suspended approximately 1 cm from the surface of the gel to ensure any US volume will contain between four and six fiducials. The pelvic segment used for our validation can be seen in Figure 4.1 (a) and 67 Chapter 4. Validation of Proposed Registration Algorithm Figure 4.1: The phantom was created with a Sawbone radio-opaque hemi pelvis, cut along red line (a). The final phantom embedded in PVC gel is seen in (b). Sawbone graphic courtesy of Pacific Research Laboratories Inc. the finished phantom submerged in PVC can be seen in (b). The finished phantom is 10 cm long and has a 10 cm diameter. The US volumes were taken with a General Electric Voluson 730 Expert Ultrasound Machine with a 3D RSP5-12 probe. The RSP5 probes are linear arrays that are mechanically swept in a 20◦ arc, alleviating the need for free- hand US, which is predominant in literature. The center image frequency is 7.5MHz. The US volumes are 152x198x148 voxels and have a resolution of .25mm3. The sweep time is dependent on the depth of the volume (in our case, 49.5 mm); but each scan takes approximately 5 seconds. 49 US scans were acquired along the iliac spine of the phantom. All phantom tests in this section (with the exception of Section 4.2.6 - Validation of the US Stitching Adaptation of Registration Algorithm) were run on all 49 datasets. Figure 4.1 is an example of the phantoms, the red line in (a) is where the phantom 68 Chapter 4. Validation of Proposed Registration Algorithm was cut to create the embedded phantom in (b). Note that bubbles formed in the phantom when the sawbone was submerged in the PVC gel; however they do not affect the PSym calculation and the fiducials are distinguishable from the air pockets because they reverberate in the US volume, see the red arrows Figure 4.2 (a). All the fiducials were sufficiently distanct from the bone surface to not affect the registration. A high resolution peripheral quantitative CT machine, model HRpQCT Xtreme CT (Scanco Medical, Wayne, Pennsylvania), was used to acquire a single 482x482x402 voxel volume with a resolution of .25mm2. Note that the surface of the phantom is radio-opaque, but the interior is radio-lucent, Figure 4.2 (b). The fiducials were omitted from registration and were only used for validation purposes. A full pelvis radiolucent phantom was used for the US stitching test, Sawbone model number 1301 (Pacific Research Laboratories, Inc., Vashon, Washington) with thirty-eight 1 mm metal fiducials fixed to the phantom surface. The phantom was submerged in a warm water bath and the US volumes were acquired with the same US machine used in the registration phantom validation, described above. Nine sets of two volumes were fused using the stitching adaptation of our proposed registration algorithm. To quantitatively validate the US stitching results, we calculated the TRE and SRE, similar to our CT-to-US validation results. The SRE is only calculated for the overlapping surface. The registration algorithm was tested with variable point cloud sizes. The point clouds were tested with and without curvature reinforcements and näıve sampling was compared to particle simulation simplification. We compare the runtime of the Barnes-Hut algorithm and the runtime of an exhaustive force calculation to justify the introduction of the Barnes-Hut algorithm into the PSim. We compare our registration method to two other point-based registration methods; ICP and CPD. ICP is the most commonly used point-based registration method; however we show it is not robust for medical images. CPD, while previously unused for CAOS systems, is a point-based method very similar to our proposed GMM registration method; however, the similarity metric used to compare the GMMs is different. We 69 Chapter 4. Validation of Proposed Registration Algorithm chose to compare our method to CPD to show our method is more robust for our application, albeit slightly slower. While the Unscented Kalman Filter registration algorithm was developed for orthpaedic US to CT registration, we do not compare our method to the Unscented Kalman filter, because the runtime of the algorithm was cited to be over 30 seconds [37] compared to our estimated 1-2 seconds runtime for GMM registration. The GMM algorithm (as well as the ICP and CPD algorithms used for comparison) were written in C++ and called using MATLAB. The PSim algorithm was implemented in MATLAB with the exception of the Barnes- Hut optimization and comparison of the exhaustive force calculation. All experiments were run on a 2.67GHz Intel(R) Core(TM) i7 CPU using 64-bit Windows 7 and 4 GB of RAM. Quantitative validation was done by calculating the SRE and TRE. The SRE measures the distance between the CT and registered US surfaces, as the name suggests, and was calculated as the root mean square of the euclidean distance of the model from the scene outline (the CT bone sur- face) was used to measure the SRE between the data sets. To prevent an artificially high SRE, only overlapping values were used to calculate the SRE. This overlap was determined to be any portion of the surface that lay directly in the direction of the probe as the second surface. The fiducial registration error for each data set was calculated as the average distance between fiducials in the CT volume and registered US volume. The fiducials were not used during registration. 4.2.2 Point Cloud Size and the Effect on Registration Runtime and Accuracy The size of the point cloud is the main consideration when minimizing the registration runtime for GMM registration. The complexity of the algo- rithm is O(nm) where n and m are the numbers of points in each point cloud. Jian et al. stated that point clouds smaller than 1000 points [56] result in a reasonable runtime, but they did not specify what they consider reasonable. In [47] we determined that 800 - 1000 points provided the best 70 Chapter 4. Validation of Proposed Registration Algorithm Figure 4.2: A single slice of a US volume (a) and corresponding CT slice (b) of the phantom. The red arrow (a) shows the reverberation caused by the fiducial. The blue arrow (a) shows an air bubble. Notice that the air bubble does not have reverberations and does not appear in the the CT image (b). 71 Chapter 4. Validation of Proposed Registration Algorithm Table 4.1: Quantitative results demonstrating the effect of point cloud size; specifically 300 and 500 points per dataset. The test was run on 49 US- CT phantom dataset pairs and point clouds were simplified using particle simulation. 500 Points 300 Points SRE Time SRE Time (mm) (s) (mm) (s) Mean 0.315 1.406 0.546 0.978 Std. Dev. 0.092 0.306 0.236 0.551 Median 0.289 1.337 0.519 0.887 Min. 0.151 1.100 0.057 0.254 Max. 0.654 2.861 1.281 2.428 trade-off between accuracy and runtime for our proposed application. In this paper, however, näıve down-sampling was used to minimize the point cloud. Including PSim in the algorithm reduces the size of the point cloud needed to accurately register the US and CT volumes. Figure 4.3 shows the registration and PSim runtime and compares it to the SRE. The SRE registration error settles to a reasonable error at approximately 300 points and the runtime increases drastically at 800 points, due to the point cloud projection step in the particle simulation algorithm. We tested the registra- tion of 300 and 500 points, the results can be seen in Table 4.2.2. While the mean accuracy did not fail in cases with 300 points, the standard deviation was significantly larger, indicating that using 500 points is more robust. 4.2.3 Particle Simulation’s Effect on Registration Robustness The method of point cloud simplification affects both the accuracy and ro- bustness of the registration algorithm. In [47], we show that it is possible to use a näıve down-sampling method; however the US volumes were of high quality and the validation was limited to 10 datasets. We tested this näıve down-sampling method and discovered that it is not suitable for all phantom volumes. This is due to large bone surfaces being omitted under certain con- 72 Chapter 4. Validation of Proposed Registration Algorithm Figure 4.3: To justify the use of particle simulation, we must compare the registration error and overall runtime. As the optimal number of points per point cloud is not known for registration attempts that utilise particle simulation, the test was re-run for different point cloud sizes. PSym was omitted from the runtime calculation, as it is not dependent on the number of points in a point cloud. 73 Chapter 4. Validation of Proposed Registration Algorithm Table 4.2: Quantitative results of test comparing registration results for tests using PSim and näıve down sampling. The test was run on 49 US-CT phantom dataset pairs and correspond to the boxplots in Figure 4.5. PSim Näıve SRE TRE Time SRE TRE Time (mm) (mm) (s) (mm) (mm) (s) Mean 0.315 0.769 1.406 0.277 6.837 1.010 Std. Dev. 0.092 0.289 0.306 0.069 5.351 0.476 Median 0.289 0.791 1.337 0.265 5.761 0.915 Min 0.151 0.250 1.100 0.181 0.803 0.252 Max 0.654 1.609 2.861 0.477 25.050 2.644 ditions, e.g. a bone surface lies in a plane which is omitted from the point cloud. PSim, however, ensures the bone surface is well represented using a smaller point cloud in comparison to down-sampling. When developing our proposed algorithm, we used qualitative inspection and comparison of runtimes to omit other methods of point cloud simplification Figure 4.4, such as K-means clustering (b) and hierarchical cluster (c). While K-means clustering produces a very well represented simplified surface, the run-time was very long (approximately 2 minutes in our implementation). Hierar- chical clustering is a simplification method that models the surface as a hierarchical tree and uses a top-down method of grouping points based on surface variation. Hierarchical clustering is very fast (> 1s); however, it does not represent the entire bone surface well, because this method emphasized points along the iliac crest, where there is greater variation, while group- ing the pelvic table into one cluster. The user cannot specify the number of points in hierarchical clustering, which made the resulting point cloud unreliable for registration. We demonstrate the registration robustness im- provement introduced by PSim by comparing the registration results of the 49 pairs of point clouds simplified with PSim and näıve down sampling. The SRE and the TRE can be seen in Table 4.2.3 and Figure 4.5. PSim keeps the registration SRE and TRE submillimetric at low point cloud sizes and provides a good trade-off in runtime. 74 Chapter 4. Validation of Proposed Registration Algorithm Figure 4.4: Comparison of point cloud simplification methods - The full point cloud extracted from an US bone surface is seen in (a). (b) shows the point cloud simplified with K-means clustering, (c) hierarchical clustering and (d) particle simulation. 75 Chapter 4. Validation of Proposed Registration Algorithm Figure 4.5: Quantitative results of registration validation, comparing the TRE and SRE of registration performed on point clouds simplified with particle simulation and down sampling. Note the improved registration re- sults for datasets simplified using PSim. The test was run on 49 US-CT phantom dataset pairs and correspond to the results seen in Table 4.2.3 76 Chapter 4. Validation of Proposed Registration Algorithm The Barnes-Hut algorithm was incorporated in the PSim algorithm to optimize the force calculation and minimize the runtime. We compare the runtime of the force calculation for different point cloud sizes to show the Barnes-Hut algorithm is more effective that an exhaustive calculation (see Figure 4.6). For 500 points, the runtime is approximately twice as fast, or .25 seconds faster. PSim calculates the force every iteration, so we see a runtime reduction of approximately 1 second, including the the time it takes to construct the Barnes-Hut tree. 4.2.4 Registration Results and the Contribution of Curvature We ran the phantom experiments with and without high curvature feature reinforcement, K, to determine their effect on the registration accuracy. The SRE error is similar between tests reinforced with K and cases without this K. Table 4.2.4 and Figure 4.7 shows the quantitative results of the regis- tration, comparing the SRE and TRE of the two tests. As expected, the fiducial registration error is larger for tests that did not include K reinforce- ments. Interestingly, the SRE is similar for both tests, which indicates that the registration algorithm matches the bone surfaces well but the surfaces can be relatively displaced along the iliac spine, as expected. The red arrows show the US fiducial and the blue arrows point to the fiducial in the CT volume; or where the fiducial should be in the case of the tests run without K. Notice that although the bone surfaces are well aligned, the fiducials are in the correct position only in the ’K’ case. These results support our proposal of introducing K reinforcement to improve the registration of the bone surface with respect to anatomical features as the results. From the results in Figure 4.7 and Table 4.2.4, we have shown that K-reinforced point clouds provide the most robust and accurate results. We present qualitative results of this algorithm in Figure 4.8. These re- sults show that including high curvature features in the set of points being matched with the GMM improves the accuracy of the final results. An in- dependent two-sample t-test was conduction to compare the impact of the 77 Chapter 4. Validation of Proposed Registration Algorithm Figure 4.6: The runtimes of force calculations comparing the Barnes-Hut algorithm and exhaustive force calculations. Data was run on a single point cloud for various desired output point cloud sizes. The test was repeated 10 times; the values above are the mean runtime. 78 Chapter 4. Validation of Proposed Registration Algorithm Figure 4.7: Quantitative results evaluating the effect of reinforcing a point cloud with high curvature features, K. Note the improvement of the TRE in cases that include K features. The test was run on 49 US-CT phantom dataset pairs and correspond to the results seen in Table 4.2.4. Table 4.3: Quantitative results of test comparing registration results for tests using PSim and down sampling. The test was run on 49 US-CT phantom dataset pairs and correspond to the boxplots in Figure 4.7. K Reinforcement noK Reinforcement SRE TRE Time SRE TRE Time (mm) (mm) (s) (mm) (mm) (s) Mean 0.315 0.769 1.406 0.340 1.690 1.132 Std. Dev. 0.092 0.289 0.306 0.127 1.653 0.245 Median 0.289 0.791 1.337 0.298 1.298 1.113 Min 0.151 0.250 1.100 0.167 0.354 0.806 Max 0.654 1.609 2.861 0.889 10.539 1.779 79 Chapter 4. Validation of Proposed Registration Algorithm reinforcement of high-curavture features on the accuracy of the results. The p-value was 0.0002 for the t-test conducted on the TRE results and 0.27 for the SRE results. There was no significant difference in the SRE between the tests that included K and those that did not; however, we can conclude that the TREs are different. This supports our claim that GMM registration matches the surfaces of the bones very well and that K points are necessary to accurately represent the bone surfaces. 4.2.5 A Comparison of GMM Registration and Other Point-Based Methods In this section we compare the registration results of two other point-matching registration methods, ICP and CPD. The same point clouds used for GMM registration in the above tests were used to test the ICP and CPD algorithms to eliminate the effect of the random point selection in the PSim step. To compare these methods, we tested the convergences over a range of rigid transformations up to 100 mm in a random direction and a 100◦ rotation around a random axis. We chose this range to ensure we tested the worst-case initial alignment caused by a poor anatomical digitization. Each of the 49 US datasets used in the previous sections was used to create the registration scenes in this comparison test. The registration model, i.e. the moving dataset in the registration, was made by imposing a specific rotation around a random axis and translation in a random direct on each scene. As the model obtained by transforming the scene, we would expect a perfect registration result to have a TRE equal to zero. The results are seen in Table 4.4. The TRE was calculated as the RMS of the distance between each pair of points in the registered model and scene data set. We observe the TRE to be 0 mm for GMM registration in Figure 4.9 for rotation and translation errors up to and including 70 mm and 70◦, with a slight increase in error up to 100 mm and 100◦. CPD; however, increases in error after a transformation of only 20 mm and 20◦, and ICP shows significant errors with a mere a 10 mm and 10◦ transformation. Figure 4.10 shows three scatter plots comparing the registration results of these 80 Chapter 4. Validation of Proposed Registration Algorithm Figure 4.8: Qualitative results of GMM Registration with PSim and Curva- ture Features. Three tests are shown, with both 3D and a single correspond- ing 2D slice. (a) is the CT data, (b) is the B-Mode US, prior to registration, (c) is the fused result of the GMM registration and (d) is fused rendering of the extracted bone surfaces. 81 Chapter 4. Validation of Proposed Registration Algorithm Figure 4.9: Comparison of GMM, ICP and CPD registration. The typically acceptable error for a number of orthopaedic surgery applications is between 2 mm (blue line) and 4 mm (red line) for the entire system [2]. This test was run on the 49 phantom US point sets previously extracted and correspond to the data seen in Figure 4.10 and Table 4.4. For each dataset, a translation and rotation was applied in a random direction to create the model dataset. tests. Note that GMM registration converges perfectly for initial errors of 3 mm and below, and few tests fail in the higher range. In contrast, ICP and CPD both fail increasingly more as the initial TRE grows. This is reflected in the results shown in Table 4.4. We tested these algorithms on the 49 US and CT volumes obtained pre- viously; the results are seen in Table 4.2.5. All methods were implemented in C++ and run on the computer described in this section. Both ICP and CPD have faster runtimes; however, neither are appropriate for registering these data sets, exemplified by the high SREs. All 49 tests fail for both ICP and CPD. It is known that ICP is susceptible to local minima, and therefore the results are likely due to initializations outside the region of at- 82 Chapter 4. Validation of Proposed Registration Algorithm Figure 4.10: Comparison of the initial vs. final TRE values for GMM, ICP and CPD. The error limits in orthopaedic surgery are between 2 mm (blue line) and 4 mm (red line) for the entire system [2]. These results correspond to Figure 4.9 and Table 4.4. 83 Chapter 4. Validation of Proposed Registration Algorithm Table 4.4: Comparison of GMM registration, ICP and CPD. For the pur- poses of this table, the fail rate was defined as 2 mm. The test was repeated for each of the 49 US datasets described in Section 4.2 and was repeated 10 times for each range tested. GMM ICP CPD TRE (mm) Fail Rate (%) TRE (mm) Fail Rate (%) TRE (mm) Fail Rate (%) ±1mm,±1◦ 0.00 0 0.17 0 0.00 0 ±5mm,±5◦ 0.00 0 0.89 0 0.06 0 ±10mm,±10◦ 0.00 0 1.77 44 0.26 0 ±20mm,±20◦ 0.00 0 4.98 92 0.95 0 ±30mm,±30◦ 0.00 0 9.30 99 2.07 63 ±40mm,±40◦ 0.00 0 12.1 100 3.65 89 ±50mm,±50◦ 0.00 0 14.0 100 5.72 97 ±60mm,±60◦ 0.00 0 14.2 86 6.91 85 ±70mm,±70◦ 0.09 2 16.4 90 10.0 90 ±80mm,±80◦ 0.45 12 20.0 100 13.6 99 ±90mm,±90◦ 0.94 25 21.8 100 16.6 100 ±100mm,±100◦ 2.00 52 22.8 99 20.1 100 84 Chapter 4. Validation of Proposed Registration Algorithm traction of the true global minimum. CPD has not previously been applied to orthopaedic US volumes. This is likely because it was originally devel- oped for analysing temporal-related registered images, such as those found in elastography [58]. The TRE was not calculated for this test, as the exact point correspondence was not known. Table 4.5: Comparison of GMM registration, ICP and CPD. Registration algorithms were tests on 49 US-CT data pairs (unlike Table 4.4). GMM Registration ICP CPD SRE Time SRE Time SRE Time (mm) (s) (mm) (s) (mm) (s) Mean 0.315 1.406 40.97 0.115 37.13 1.125 Std. Dev. 0.092 0.306 13.14 0.044 13.51 0.329 Median 0.289 1.337 41.65 0.099 41.02 1.263 Min 0.151 1.100 15.39 0.065 4.57 0.298 Max 0.654 2.861 79.10 0.246 56.86 1.579 4.2.6 Validation of the US Stitching Adaptation of Registration Algorithm Nine pairs of US volumes were stitched using our proposed registration algo- rithm. Due to the similarity in the US (i.e., same pre-processing technique and resolution), fewer points are needed to successfully register the volumes resulting in a very low registration runtime. This does not make up for the extra time needed to compute the PSym on both volumes, however. We chose to compare the registration error and runtime between 300 and 500 points in the point cloud because the plot in Figure 4.3 indicates the this range provides a good trade-off between accuracy and runtime. Figure 4.11 demonstrates the qualitative results of the US stitching tests. In US stitching, both the scene and model are of the same modality, and therefore are inherently similar in nature, resulting in the low SRE and TRE observed in Table 4.6 and Figure 4.12. It can be observed that registration attempts run using 300 points results in a TRE that is slightly higher, but are still 85 Chapter 4. Validation of Proposed Registration Algorithm Table 4.6: Quantitative results of test comparing registration results for US stitching. Tests were run on nine pairs of US data; the qualitative results can be seen in Figure 4.11. 300 Points 500 Points TRE SRE time TRE SRE time (mm) (mm) (s) (mm) (mm) (s) Average 1.20 0.45 0.89 1.06 0.44 1.04 Std. Dev 0.70 0.29 0.48 0.74 0.21 0.27 Median 1.02 0.32 0.80 0.77 0.33 0.96 Min 0.46 0.21 0.40 0.28 0.18 0.69 Max 3.02 1.08 2.09 3.14 0.87 1.44 well within an acceptable error. Therefore, we can decrease the runtime on average by 0.15 seconds by reducing the size of the point cloud by 200 points. 4.3 Clinical Evaluation of GMM Registration In this section, we will describe the clinical validation of our proposed US to CT registration. The goal of this validation was to determine if our algorithm met the requirements outlined in the introduction of this chapter for datasets that contained soft tissue. Soft tissue can obscure the bone surface and make it very difficult to localize the bone. A clinical study was done on twelve pelvic fracture patients at Vancouver General Hospital (VGH). Each of these twelve patients required a pre- operative CT scan to diagnose the fracture and plan the surgery. US volumes were taken in one of two places: in the OR prior to surgery when the patient was under general anaesthetic or post-surgery in the ward while the patient was recovering. For patients scanned post-surgery, only the non-fractured side was scanned to minimize patient’s discomfort. The same G.E. Voluson 730 Expert Ultrasound Machine and 3D RSP5- 12 probe used in the phantom study was used for this clinical study. 86 Chapter 4. Validation of Proposed Registration Algorithm Figure 4.11: Qualitative results showing two of the nine US tests used to validate the proposed algorithm using US Stitching. Quantitative results can be seen in Figure 4.12 and Table 4.6. 87 Chapter 4. Validation of Proposed Registration Algorithm Figure 4.12: Quantitative results comparing the runtime and error of US stitching for points clouds of size 300 and 500. Nine pairs of US volumes were acquired from the phantom for this test and correspond to the qualitative results seen in Figure 4.11 and quantitative results seen in Table 4.6. The quantitative results of the clinical tests can be seen in Table 4.3. Twenty-four tests were run, twelve for automatically optimized parameters and twelve for empirical parameters. Only one test failed. This failed test was run with empirical parameters. These results imply that not only that the algorithm is robust, but the accuracy of the registration results improve when the parameters are automatically optimized for clinical cases. While empirical parameters were sufficient for the phantom study in Section 4.2, the soft-tissue response in the clinical volumes degrades the PSym results and, therefore, the registration results. It should be noted that the runtime required to automatically optimize the parameters for PSym increased the PSym runtime by approximately ten fold. This will be addressed in Chapter 5, Section 5.3: Future Work. 88 Chapter 4. Validation of Proposed Registration Algorithm Figure 4.13: Results of clinical US to CT registration. The top row is the registration results using automatically optimized phase symmetry pa- rameters to localize the US bone surface. The bottom row uses empirical parameters to calculate the phase symmetry. Note that the bottom row is misregistered. Quantitative results can be seen in Table 4.3. Table 4.7: Quantitative results for clinical tests comparing automatically and empirically optimized parameters PSym. One non-optmized 3D PSym test failed. Optimized PSym non-Optimized PSym SRE SRE (mm) (mm) Average 0.58 0.89 Std. Dev. 0.39 1.02 Median 0.68 0.54 Min 0.02 0.03 Max 1.35 Fail 89 Chapter 5 Conclusions In this chapter we discuss the significance of the work and the conclusions of this thesis. We also discuss the future work required to continue moving towards integrating our proposed registration algorithm in the CAOS system discussed in this thesis. 5.1 Discussion In this thesis, we have proposed a CAOS system to use US during pelvic fracture surgery, in which US is fused with pre-operative CT volumes. The US volumes would be used to update the pelvic model during surgery, ne- cessitating a near-realtime registration algorithm. Doing so would minimize exposure to ionizing radiation for the surgeon, staff and patient due to intra- operative fluoroscopy. We have presented a novel approach to registering US and CT to con- tribute to this proposed system. We have applied a method of automatically localizing US bone surfaces. To ensure a fast runtime, we have simplified the extracted bone surface to a small point cloud containing 500 points by using particle simulation in this novel application. To further reduce the regis- tration runtime of this algorithm, we have optimized the force calculations in the PSim algorithm by integrating the Barnes-Hut algorithm. We have incorporated high curvature features, extracted using Gaussian Curvature, to reinforce salient anatomical features and thus prevent misregistration of the bone surfaces. The extracted point clouds were registered using GMM registration, a fast and robust point-based method, which to our knowledge has not previously been used for medical images. We have validated our proposed algorithm on 49 phantom data sets 90 Chapter 5. Conclusions and twelve clinical data sets taken from patients with pelvic fractures. We compared the accuracy and runtime for point sets created with and without particle simulation and high curvature feature reinforcements, and also com- pared 2D and 3D PSym. While the dimensionality of the PSym had little effect, we can conclude that high curvature feature reinforcement increases the accuracy of the results, and particle simulation increases the robust- ness of the algorithm and reduces the runtime of the overall algorithm by reducing the number of points required. We have demonstrated a mean reg- istration runtime of 1.4 seconds and a mean SRE of 0.769 mm. We observed a significantly larger average TRE (approximately 3x the mean SRE). This is because the fiducials used to measure the TRE were approximately 1-3 cm from the surface of the bone and therefore were more greatly affected by small rotational misalignments. While our runtime is not yet under 5 seconds for the entire algorithm, we have proposed a robust, accurate and fully automatic registration method that addresses three of the four require- ments stated in the beginning of Chapter 4. The fourth requirement will be discussed in Section5.3: Future Work. We demonstrated that PSim represents the bone surface better than näıve down-sampling because näıve down-sampling runs the risk of ommit- ting entire bone surfaces. It should be noted that the SRE is lower for nïıve down-sampling than registration completed with PSim. This is due to the random nature of the PSim and the low error and standard deviation of the SRE results. We have shown considerable improvement to the state-of-the-art US-to- CT for CAOS applications. Our accuracy and robustness is comparable to Moghari’s implementations of the Unscented Kalman Filter registration algorithm, however we have shown our runtime to be approximately 15x faster [37]. While Gill et al. have shown great results with implementing US on the GPU [28], their registration runtime is still approximately 20 minutes. We have managed to show a considerable improvement in runtime, even including the runtime of our PSym implementation. We demonstrated out method improved the accuracy robustnessof both CPD and ICP in Chapter 4, a significant improvement over the ICP registration results reported by 91 Chapter 5. Conclusions Penney et al., Carrat et al. and Tonettie et al. [30–32]. 5.2 Thesis Contributions The goal of our work was to contribute a method of multi-modal image reg- istration to a guidance and navigation system for pelvic ring injury surgery. This system will use US and CT imaging to assist a surgeon during surgery. To this end, we contributed the following to this project: 1. We outlined and described a proposed orthopaedic guidance system for pelvic fracture surgeries and how it is to be used. We also described the components necessary for the final system and its requirements. 2. To meet one of the requirements in the aforementioned system, we proposed the novel use of a 3D registration algorithm not previously implemented in medical imaging. This registration algorithm uses GMM to model two sets of points extracted from the US and CT vol- umes and computes the rigid transformation by iteratively minimizing a measure of distance between them. 3. We applied previous research in US localization to extract the bone surface from US volumes to be used in the above registration algo- rithm. We showed that this localization method can be effectively used to segment the bone surface in 3D US data sets. We also demon- strate that automatically optimized parameters significantly improve the results of the registration algorithm. 4. To use the bone surfaces extracted from the US and CT volumes, we proposed to use PSim to simplify the surface into a small point cloud. To our knowledge, PSim has never been applied to pre-processing for image registration or to medical data. 5. We reduced the PSim runtime by optimizing the force calculations in the above point cloud simplification method by implementing a Barnes-Hut hierarchical tree, previously applied to N -Body simula- tions in astrophysics. 92 Chapter 5. Conclusions 6. We extended our multi-modal registration method to demonstrate US- to-US stitching, which contributed to another requirement of the pro- posed surgical system. 5.3 Future Work The registration algorithm presented in this paper is one component of the proposed CAOS system to assist surgeons during pelvic fracture surgery. In this thesis we have shown the registration to be accurate, fast and robust. However the slowest step is currently the PSym runtime. PSym is an im- portant step in our algorithm and cannot be omitted, and therefore we are currently implementing the 3D PSym algorithm to be optimized to run on a GPU. This implementation with use CUDA to exploit the parallelizability of the algorithm. The current runtime for the 3D PSym in MATLAB is over twenty-five seconds, and initial results lead us to believe that we will see an improvement factor of at least 500x after the code is optimized. The proposed GMM registration algorithm takes less than two seconds to run, and future work will include adding intensity and weight information to the dataset to increase robustness to initial alignment, which will increase the runtime and may require the code to be optimized as well. This im- provement will likely also decrease the overall error in our method. While we have demonstrated excellent robustness and accuracy in our phantom study, there is room for improvement in clinical datasets. Our method has five free parameters: the force constant, k; in the Par- ticle simulation, θ in the Barnes-Huts describing the distance of a cluster of points compared to the size of the volume division; the number of points in the point cloud; the number of curvature reinforcement points; and the number of threshold value for the CT segmentation. Of these parameters, only the final size of the point cloud used in the GMM registration algorithm greatly affects the outcome. In the future, we plan to estimate the optimal number of points from the size of the US volume and the image resolution. In Chapter 4, Section 4.2.2, we showed 500 points to be the optimal trade-off between runtime and accuracy. The greatest contribution to the runtime of 93 Chapter 5. Conclusions the PSim algorithm is currently the projection step completed after every iteration. In the future, we propose to optimize this step by incorporating k-dimensional trees [59]. While we have demonstrated that GMM registration is robust, accurate and fast, incorporating intensity and structure information may improve the overall results. We may also benefit from investigating the method of segmenting the CT volume, to ensure it is robust for other anatomical regions of interest. It should be noted that any pre-processing done on the CTcan be done pre-operatively and therefore can have a longer runtime that the US pre-processing. Further work should be done to extract the bone surface from the PSym volume and perhaps incorporate regularization to ensure the bone surface is continuous. We have demonstrated that this work can be applied to US stitching; however, we have not conducted a clinical study and therefore this requires further investigation to validate our method for this purpose. In the valida- tion completed in this work, we only measured the the SRE of the overlapped area, which likely is artificially low. We also expect, providing the runtime is not too adversely affected, to incorporate the US stitching into the regis- tration algorithm to further improve the robustness of the system in terms of initial alignment and local minima due to the US probes’ small ROI and the similarity of anatomical features along the iliac crest. This registration work has been tested on clinical data, however, a full cadaver study must be performed before it can be fully incorporated into the CAOS system. This cadaver study should involve tracking the probe and possibly the surgical instrumentation for the K-wire drilling and cannulated screw over-drilling. Ultimately, we hope to see this work employed in the proposed CAOS system, with the goal of minimizing the radiation time and possible surgical errors during pelvic surgeries. Our proposed algorithm needs to be adapted into an articulated registration method to account for pelvic fragment move- ment. This information will need to be updated and reflected in the pelvic model displayed to the surgeon during surgery. 94 Bibliography [1] A. D. Wiss, Fractures. Master Techniques in Orthopaedic Surgery, Van Nuys, California: Lippincott-Raven, 1998. [2] R. Phillips, “The accuracy of surgical navigation for orthopaedic surgery,” in Current Orthopedics, vol. 21, pp. 180–192, 2007. [3] T. Rüedi, R. Buckley, and C. Moran, AO Principles of Fracture Man- agement, vol. 2-Specific Fractures. Boca Raton, FL: Thieme, 2 ed., 2007. [4] K. A. Lefaivre, A. J. Starr, B. P. Barker, and C. M. Reinert, “Early experience with reduction of displaced disruption of the pelvic ring using a pelvic reduction frame,” Journal of Bone and Joint Surgery, vol. 91-B, pp. 1201–1207, 2009. [5] N. Joseph, “Online radiography continuing education for radiologic x- rey technologist,” 2007. [6] H. J. M. D. W. L. Z. M. Balter, S., “Fluoroscopically guided interven- tional procedures: A review of radiation effects on patients’ skin and hair,” Radiology, vol. 254, pp. 326–341, Feb. 2010. [7] M. D. B. S. L. W. D. H. N. A. M. M. G. J. Klein, W.L., “Occupa- tional health hazards in the interventional laboratory: Time for a safer environment,” Radiology, vol. 250, pp. 538–544, Feb. 2009. [8] G. Penney, J. Weese, P. Batchelor, D. Hill, and D. Hawkes, “Valida- tion of a two- to three-dimensional registration algorithm for aligning preoperative ct images and intraoperative fluoroscopy images,” Medical Physics, vol. 28, no. 6, pp. 1024–1032, 2001. 95 Bibliography [9] D. Tomazevic, B. Likar, T. Slivnik, and F. Pernus, “3-d/2-d registration of ct and mr to x-ray images,” Medical Imaging, IEEE Transactions on, vol. 22, pp. 1407 –1416, Nov. 2003. [10] D. Russakoff, T. Rohlfing, K. Mori, D. Rueckert, A. Ho, J. Adler, J.R., and J. Maurer, C.R., “Fast generation of digitally reconstructed radio- graphs using attenuation fields with application to 2d-3d image reg- istration,” Medical Imaging, IEEE Transactions on, vol. 24, pp. 1441 –1454, Nov. 2005. [11] B. Brendel, S. Winter, a. Rick, M. Stockheim, and H. Ermert, “Reg- istration of 3D CT and ultrasound datasets of the spine using bone structures.,” Computer aided surgery : official journal of the Interna- tional Society for Computer Aided Surgery, vol. 7, pp. 146–55, Jan. 2002. [12] B. Brendel, J. Siepermann, S. Winter, and H. Ermert, “In vivo evalua- tion and in vitro accuracy measurements for an ultrasound-CT registra- tion algorithm,” International Congress Series, vol. 1281, pp. 583–588, May 2005. [13] S. Winter, B. Brendel, A. Rick, M. Stockheim, K. Schmieder, and H. Ermert, “Registration of bone surfaces, extracted from CT-datasets, with 3D ultrasound.,” Biomedizinische Technik. Biomedical engineer- ing, vol. 47, pp. 57–60, Jan. 2002. [14] C. Dekomien, M. Mildenstein, K. Hensel, S. Hold, and S. Winter, “Reg- istration of bone structures in 3D ultrasound and CT data: Comparison of different optimization strategies,” Advances in Medical Engineering, vol. 114, pp. 57–60, May 2007. [15] S. Winter, B. Brendel, I. Pechlivanis, K. Schmieder, and C. Igel, “Regis- tration of ct and intraoperative 3-d ultrasound images of the spine using evolutionary and gradient-based methods,” Evolutionary Computation, IEEE Transactions on, vol. 12, pp. 284 –296, june 2008. 96 Bibliography [16] S. Winter, I. Pechlivanis, C. Dekomien, C. Igel, and K. Schmieder, “Toward registration of 3d ultrasound and ct images of the spine in clinical praxis: Design and evaluation of a data acquisition protocol,” Ultrasound in Medicine and Biology, vol. 35, no. 11, pp. 1773–1782, 2009. [17] G. Penney, D. Barratt, C. Chan, M. Slomczykowski, T. Carter, P. Ed- wards, and D. Hawkes, “Cadaver validation of intensity-based ultra- sound to ct registration,” in Medical Image Computing and Computer- Assisted Intervention MICCAI 2005 (J. Duncan and G. Gerig, eds.), vol. 3750 of Lecture Notes in Computer Science, pp. 1000–1007, Springer Berlin / Heidelberg, 2005. [18] S. Ghanavati, P. Mousavi, G. Fichtinger, P. Foroughi, and P. Abol- maesumi, “Multi-slice to volume registration of ultrasound data to a statistical atlas of human pelvis,” Medical Imaging 2010: Visualiza- tion, Image-Guided Procedures, and Modeling, vol. 7625, no. 1, 2010. [19] W. Wein, A. Khamene, D. Clevert, O. Kutter, and N. Navab, “Sim- ulation and fully automatic multimodal registration of medical ultra- sound,” in Medical Image Computing and Computer-Assisted Interven- tion - MICCAI 2007 (N. Ayache, S. Ourselin, and A. Maeder, eds.), vol. 4791 of Lecture Notes in Computer Science, pp. 136–143, Springer Berlin / Heidelberg, 2007. [20] W. Wein, S. Brunke, A. Khamene, M. R. Callstrom, and N. Navab, “Automatic CT-ultrasound registration for diagnostic imaging and image-guided intervention,” Medical Image Analysis, vol. 12, pp. 577– 585, 2008. [21] J. Jensen, “A model for the propagation and scattering of ultrasound in tissue,” The Journal of the Acoustical Society of America, vol. 89, no. 1, pp. 182–190, 1991. [22] J. Jensen and I. Nikolov, “Fast simulation of ultrasound images,” in Ultrasonics Symposium, 2000 IEEE, vol. 2, pp. 1721 –1724, Oct. 2000. 97 Bibliography [23] O. Kutter, A. Karamalis, W. Wein, and N. Navab, “A gpu-based frame- work for simulation of medical ultrasound,” in SPIE Medical Imaging, (Orlando, Florida, USA), February 2009. [24] R. Shams, R. Hartley, and N. Navab, “Real-time simulation of med- ical ultrasound from ct images,” in Medical Image Computing and Computer-Assisted Intervention 2008, vol. 5242 of Lecture Notes in Computer Science, pp. 734–741, Springer Berlin / Heidelberg, 2008. [25] S. Khallaghi, P. Abolmaesumi, R. H. Gong, E. C. S. Chen, S. Gill, J. Boisvert, D. Pichora, D. Borschneck, G. Fichtinger, and P. Mousavi, “Gpu accelerated registration of a statistical shape model of the lumbar spine to 3d ultrasound images,” in SPIE Medical Imaging (K. H. Wong and D. R. H. III, eds.), vol. 7964, (Lake Buena Vista, Florida, USA), p. 79642W, SPIE, SPIE, 03/2011 2011. [26] O. Fluck, C. Vetter, W. Wein, A. Kamen, B. Preim, and R. Wester- mann, “A survey of medical image registration on graphics hardware,” Computer Methods and Programs in Biomedicine, Nov 2010. [27] T. Reichl, J. Passenger, O. Acosta, and O. Salvado, “Ultrasound goes gpu: real-time simulation using cuda,” in SPIE Medical Imaging, vol. 7261, February 2009. [28] S. Gill, P. Abolmaesumi, G. Fichtinger, J. Boisvert, D. Pichora, D. Bor- shneck, and P. Mousavi, “Biomechanically constrained groupwise ultra- sound to ct registration of the lumbar spine,” Medical Image Analysis, 2010. [29] D. Barratt, C. Chan, P. Edwards, G. Penney, M. Slomczykowski, T. Carter, and D. Hawkes, “Instantiation and registration of statistical shape models of the femur and pelvis using 3D ultrasound imaging.,” Medical image analysis, vol. 12, pp. 358–74, June 2008. [30] G. Penney, P. Edwards, A. King, J. Blackall, P. Batchelor, and D. Hawkes, “A stochastic iterative closest point algorithm (stochas- 98 Bibliography ticp),” in Medical Image Computing and Computer-Assisted Interven- tion MICCAI 2001 (W. Niessen and M. Viergever, eds.), vol. 2208 of Lecture Notes in Computer Science, pp. 762–769, Springer Berlin / Heidelberg, 2001. [31] L. Carrat, J. Tonetti, P. Merloz, and J. Troccaza, “Percutaneous com- puter assisted iliosacral screwing: Clinical validation,” in Medical Im- age Computing and Computer-Assisted Intervention MICCAI 2000 (S. Delp, A. DiGoia, and B. Jaramaz, eds.), vol. 1935 of Lecture Notes in Computer Science, pp. 97–140, 2000. [32] J. Tonetti, L. Carrat, S. Blendea, P. Merloz, J. Troccaz, S. Lavallée, and J.-P. Chirossel, “Clinical results of percutaneous pelvic surgery. computer assisted surgery using ultrasound compared to standard flu- oroscopy,” Computer Aided Surgery, vol. 6, no. 4, pp. 204–211, 2001. [33] D. V. Amin, T. Kanade, A. M. Digioia, and B. Jaramaz, “Ultrasound registration of the bone surface for surgical navigation,” Computer Aided Surgery, vol. 8, no. 1, pp. 1–16, 2003. [34] D. Rueckert, A. Frangi, and J. Schnabel, “Automatic construction of 3d statistical deformation models of the brain using nonrigid registration,” Medical Imaging, IEEE Transactions on, vol. 22, pp. 1014 –1025, aug. 2003. [35] D. Barratt, G. Penney, C. Chan, M. Slomczykowski, T. Carter, P. Ed- wards, and D. Hawkes, “Self-calibrating 3d-ultrasound-based bone reg- istration for minimally invasive orthopedic surgery,” Medical Imaging, IEEE Transactions on, vol. 25, pp. 312–323, march 2006. [36] P. Foroughi, D. Song, G. Chintalapani, R. Taylor, and G. Fichtinger, “Localization of pelvic anatomical coordinate system using us/atlas registration for total hip replacement,” vol. 5242, pp. 871–879. [37] M. Moghari and P. Abolmaesumi, “A novel incremental technique for ultrasound to ct bone surface registration using unscented kalman fil- 99 Bibliography tering,” in Medical Image Computing and Computer-Assisted Interven- tion MICCAI 2005 (J. Duncan and G. Gerig, eds.), vol. 3750 of Lecture Notes in Computer Science, pp. 197–204, Springer Berlin / Heidelberg, 2005. [38] M. Moghari and P. Abolmaesumi, “Point-based rigid-body registration using an unscented kalman filter,” in IEEE Transactions on Medical Imaging, vol. 26, pp. 1708–1728, Springer Berlin / Heidelberg, 2007. [39] I. Hacihaliloglu, R. Abugharbieh, A. Hodgson, and R. Rohling, “Bone segmentation and fracture detection in ultrasound using 3d local phase features,” in Medical Image Computing and Computer-Assisted Inter- vention MICCAI 2008 (D. Metaxas, L. Axel, G. Fichtinger, and G. Szkely, eds.), vol. 5241 of Lecture Notes in Computer Science, pp. 287–295, Springer Berlin / Heidelberg, 2008. [40] K. Rajpoot, V. Grau, and J. Noble, “Local-phase based 3d boundary detection using monogenic signal and its application to real-time 3- d echocardiography images,” in Biomedical Imaging: From Nano to Macro, 2009. ISBI ’09. IEEE International Symposium on, pp. 783– 786, 28 2009-july 1 2009. [41] I. Hacihaliloglu, Towards a Novel minimally Invasive 3D Ultrasound Imaging Based Computer Assisted Orthopaedic Surgery System for Bone Fracture Reduction. PhD thesis, The University of British Columbia, 2010. [42] I. Hacihaliloglu, R. Abugharbieh, A. Hodgson, and R. Rohling, “Au- tomatic adaptive parameterization in local phase feature-based bone segmentation in ultrasound.” 2011. [43] P. Kovesi, “Symmetry and asymmetry from local phase,” in Tenth Aus- tralian Joint Converence on Artificial Intelligence, pp. 2–4, 1997. [44] I. Hacihaliloglu, R. Abugharbieh, A. Hodgson, and R. Rohling, “Bone surface localization in ultrasound using image phase-based features,” 100 Bibliography Ultrasound in Medicine & Biology, vol. 35, no. 9, pp. 1475–1487, 2009. In Press. [45] I. Hacihaliloglu, R. Abugharbieh, A. Hodgson, and R. Rohling, “Au- tomatic data-driven parameterization for phase-based bone localiza- tion in us using log-gabor filters,” in Advances in Visual Computing (G. Bebis, R. Boyle, B. Parvin, D. Koracin, Y. Kuno, J. Wang, J.-X. Wang, J. Wang, R. Pajarola, P. Lindstrom, A. Hinkenjann, M. En- carnao, C. Silva, and D. Coming, eds.), vol. 5875 of Lecture Notes in Computer Science, pp. 944–954, Springer Berlin / Heidelberg, 2009. [46] N. Sugano, T. Sasama, Y. Nakajima, Y. Sato, T. Nishii, T. Iida, K. Nak- agawa, K. Ono, S. Nishihara, S. Tamura, K. Yonenobu, and T. Ochi, “Effects of CT threshold value to make a surface bone model on ac- curacy of shape-based registration in a CT-based navigation system for hip surgery,” International Congress Series, vol. 1230, pp. 319–324, June 2001. [47] A. Brounstein, I. Hacihaliloglu, P. Guy, A. Hodgson, and R. Abughar- bieh, “Towards real-time 3d us to ct bone image registration using phase and curvature feature based gmm matching,” in Medical Image Computing and Computer-Assisted Intervention 2011 (G. Fichtinger, A. Martel, and T. Peters, eds.), vol. 6891 of Lecture Notes in Computer Science, pp. 235–242, Springer Berlin / Heidelberg, 2011. [48] M. Pauly, M. Gross, and L. Kobbelt, “Efficient simplification of point- sampled surfaces,” in Visualization, 2002. VIS 2002. IEEE, pp. 163– 170, nov. 2002. [49] J. Barnes and P. Hut, “A hierarchical O(N log N) force-calculation algorithm,” Nature, vol. 324, pp. 446–449, Dec. 1986. [50] J. P. Singh, C. Holt, T. Totsuka, A. Gupta, and J. L. Hennessy, “Load balancing and data locality in adaptive hierarchical n-body methods: Barnes-hut, fast multipole, and radiosity,” Journal Of Parallel and Dis- tributed Computing, vol. 27, pp. 118–141, 1995. 101 Bibliography [51] J. K. Salmon, Parallel hierarchical N-body methods. PhD thesis, Pasadena, CA, USA, 1991. UMI Order No. GAX91-37285. [52] L. Nyland and J. Prins, “Fast n-body simulation with cuda,” Simula- tion, vol. 3, pp. 677–696, 2007. [53] J. Singh, J. Hennessy, and A. Gupta, “Scaling parallel programs for multiprocessors: methodology and examples,” Computer, vol. 26, pp. 42 –50, jul 1993. [54] A. Gray, E. Abbena, and S. Salamon, Modern Differential Geometry of Curves and Surfaces with Mathematica. Boca Raton, FL: CRC Press, 2 ed., 1997. [55] B. Jian and B. Vemuri, “A robust algorithm for point set registration using mixture of gaussians,” in Computer Vision, 2005. ICCV 2005. Tenth IEEE International Conference on, vol. 2, pp. 1246 – 1251, oct. 2005. [56] B. Jian and B. Vemuri, “Robust point set registration using gaussian mixture models,” Pattern Analysis and Machine Intelligence, IEEE Transactions on, vol. 33, pp. 1633–1645, aug. 2011. [57] J. Spencer, R. Day, K. Sloan, and R. Beaver, “Computer navigation of the acetabular component,” Journal of Bone and Joint Surgery, vol. 88- B, pp. 972–975, 2006. [58] B. Chuang, A. Myronenko, R. English, and A. Noble, “Interpreting ultrasound elastography: Image registration of breast cancer ultrasound elastography to histopathology images,” in Biomedical Imaging: From Nano to Macro, 2010 IEEE International Symposium on, pp. 181 –184, april 2010. [59] J. L. Bentley, “Multidimensional binary search trees used for associative searching,” Commun. ACM, vol. 18, pp. 509–517, September 1975. 102