UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Robust bone detection in ultrasound using combined strain imaging and envelope signal power detection Hussain, Mohammad Arafat 2015

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2015_september_hussain_mohammadarafat.pdf [ 14.2MB ]
Metadata
JSON: 24-1.0166292.json
JSON-LD: 24-1.0166292-ld.json
RDF/XML (Pretty): 24-1.0166292-rdf.xml
RDF/JSON: 24-1.0166292-rdf.json
Turtle: 24-1.0166292-turtle.txt
N-Triples: 24-1.0166292-rdf-ntriples.txt
Original Record: 24-1.0166292-source.json
Full Text
24-1.0166292-fulltext.txt
Citation
24-1.0166292.ris

Full Text

Robust Bone Detection in UltrasoundUsing Combined Strain Imaging andEnvelope Signal Power DetectionbyMohammad Arafat HussainB.Sc., Bangladesh University of Engineering and Technology, 2011M.Sc., Bangladesh University of Engineering and Technology, 2013A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENT FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinTHE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES(Biomedical Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)June 2015c© Mohammad Arafat Hussain, 2015AbstractUsing ultrasound for tool navigation in an orthopaedic surgery requires bone lo-calization in ultrasound, a procedure which remains challenging despite encourag-ing advances in its practice. Current methods, e.g., the local image phase-basedfeature analysis, have shown promising results but continue to rely on delicateparameter selection processes and to be prone to error at confounding soft tissueinterfaces of similar appearance to bone interfaces. In addition, the 3-dimensionalphase-features-based automatic bone segmentation method is found to be timeinefficient (at ∼2 minutes).We have proposed a different approach to bone segmentation by combiningultrasound strain imaging and envelope power detection methods. After an es-timation of the strain and envelope power maps, we subsequently fused thesemaps into a single combined map that corresponds robustly to the actual boneboundaries. This study has achieved a marked reduction in false positive boneresponses at the soft tissue interfaces. We also incorporated the depth-dependentcumulative power of the envelope into the elastographic data as well as incor-porated an echo de-correlation measurement-based weight to fuse the strain andenvelope map. We also employed a data driven scheme to detect the presence ofany bone discontinuity in the scanned ultrasound image and introduced a multi-variate non-parametric Gaussian mixture regression to be used over the maximumintensity points of the fused map. Finally, we developed a simple yet effectivemeans to perform 3-dimensional bone surface extractions using a surface growingapproach that is seeded from the 2-dimensional bone contours.We employed mean absolute error calculations between the actual and es-timated bone boundaries to show the extent of the false positives created; ourmethods showed an average improvement in the mean absolute error of 20% onboth the 2- and 3-dimensional finite-element-models, and of 18% and 23%, re-spectively, on the 2- and 3-dimensional experimental phantom data, when com-pared with that of the local phase-based methods. Validation on the 2- and3-dimensional clinical in vivo data also demonstrates, respectively, an averageimprovement in the mean absolute fitting error of 55% and an 18 times improve-ment in the computation time.iiPrefaceContributions:This research was conducted under the supervision of my Principal-supervisorProf. Rafeef Abugharbieh, PhD and Co-supervisor Prof. Antony J. Hodgson,PhD in UBC’s Biomedical Signal and Image Computing Laboratory (Director:Prof. Rafeef Abugharbieh, PhD). My role included developing the research ques-tion, hypotheses and methodological design, completing the ethics application,conducting the literature review, preparing the experimental setups, organizingtechnology requirements, recruiting participants, collecting and analyzing thedata, and drafting this written document.Publications:A version of the chapter 2 has been published. Hussain, M.A., Hodgson, A.J., andAbugharbieh, R. “Robust Bone Detection in Ultrasound Using Combined StrainImaging and Envelope Signal Power Detection”, Medical Image Computing andComputer-Assisted Intervention (MICCAI), Cambridge-USA, September, Pages:356–363, 2014. I conducted all the testing and wrote the manuscript. My super-visors provided the valuable suggestions on the methodology, and reviewed andcorrected the manuscript.Some part of the chapter 3 has been accepted for publication. Hussain, M.A.,Pierre, G., Hodgson, A.J., and Abugharbieh, R. “Automatic Bone Segmentationin Ultrasound using Combined Ultrasound Strain Imaging and Envelope SignalPower”, The 2015 Computer Assisted Orthopedic Surgery (CAOS) Annual Meet-ing, Vancouver, Canada, June 17–20, 2015. I conducted all the testing andwrote the manuscript. My supervisors provided the valuable suggestions on themethodology, and reviewed and corrected the manuscript.Finally, a combination of the chapters 4 and 5 will be submitted as a full journalarticle which is now under preparation.iiiEthics:This research study met the criteria for a Minimal Risks Human Ethics applica-tion, and thus an expedited review by the UBC Clinical Research Ethics Board(CREB) was conducted. UBC CREB number is H15–00154.Conflict of Interest:The researchers and members of the thesis committee report no conflict of inter-est.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxii1 Introduction and Background . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation and Problem Statement . . . . . . . . . . . . . . . . . 11.2 Medical Imaging in Orthopedic Surgery . . . . . . . . . . . . . . . 31.2.1 Pre-operative Imaging . . . . . . . . . . . . . . . . . . . . 41.2.2 Intra-operative Imaging . . . . . . . . . . . . . . . . . . . 91.3 The Basic Physics of Ultrasound Imaging . . . . . . . . . . . . . . 111.3.1 Sound Properties . . . . . . . . . . . . . . . . . . . . . . . 121.3.2 The Generation and Reception of Ultrasound Pulses . . . . 131.3.3 Ultrasound-tissue Interaction . . . . . . . . . . . . . . . . 181.3.4 Ultrasound Beamforming . . . . . . . . . . . . . . . . . . . 201.3.5 Ultrasound B-mode Image Formation . . . . . . . . . . . . 211.4 Challenges in Ultrasound Image Interpretation . . . . . . . . . . . 23v1.4.1 Speckle Noise . . . . . . . . . . . . . . . . . . . . . . . . . 241.4.2 Reverberation . . . . . . . . . . . . . . . . . . . . . . . . . 241.4.3 Attenuation/Shadowing . . . . . . . . . . . . . . . . . . . 261.4.4 Mirror Images . . . . . . . . . . . . . . . . . . . . . . . . . 261.4.5 Anisotropy . . . . . . . . . . . . . . . . . . . . . . . . . . . 261.5 Ultrasound-based CAOS Systems . . . . . . . . . . . . . . . . . . 261.6 Current Ultrasound-based Bone Segmentation Methods . . . . . . 291.6.1 B-mode Imaging-based Methods . . . . . . . . . . . . . . . 291.6.2 Strain-imaging-based Methods . . . . . . . . . . . . . . . . 351.7 Thesis Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . 371.8 Thesis Organization . . . . . . . . . . . . . . . . . . . . . . . . . . 372 Bone Detection in Ultrasound Using Combined Strain-imagingand Envelope Power . . . . . . . . . . . . . . . . . . . . . . . . . . . 402.1 Ultrasound Signal Model . . . . . . . . . . . . . . . . . . . . . . . 402.2 Modified Strain Map (MSM) . . . . . . . . . . . . . . . . . . . . . 412.2.1 Displacement Estimation . . . . . . . . . . . . . . . . . . . 412.2.2 Strain Estimation and Processing . . . . . . . . . . . . . . 432.3 Modified Envelope Map (MEM) . . . . . . . . . . . . . . . . . . . 432.4 Fused Map Estimation . . . . . . . . . . . . . . . . . . . . . . . . 442.5 Forming Bone Contour Using Regression . . . . . . . . . . . . . . 452.6 Validation Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . 452.6.1 Simulated Phantom . . . . . . . . . . . . . . . . . . . . . . 452.6.2 Experimental Physical Phantom . . . . . . . . . . . . . . . 462.6.3 In Vivo Data . . . . . . . . . . . . . . . . . . . . . . . . . 472.7 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 472.7.1 FEM Results . . . . . . . . . . . . . . . . . . . . . . . . . 482.7.2 Experimental Phantom Results . . . . . . . . . . . . . . . 492.7.3 In Vivo Results . . . . . . . . . . . . . . . . . . . . . . . . 512.8 Advantages and Limitations . . . . . . . . . . . . . . . . . . . . . 563 Improved 2D Bone Contour Extraction Using Data-driven Pa-rameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59vi3.1 Improved Modified Envelope Map (MEM) and Modified StrainMap (MSM) Estimation . . . . . . . . . . . . . . . . . . . . . . . 603.2 Data-driven Weight Estimation . . . . . . . . . . . . . . . . . . . 613.3 Bone Discontinuity Detection . . . . . . . . . . . . . . . . . . . . 633.4 Gaussian Mixture Regression (GMR) . . . . . . . . . . . . . . . . 643.5 Validation Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . 663.5.1 Simulated Phantom . . . . . . . . . . . . . . . . . . . . . . 663.5.2 Experimental Phantom . . . . . . . . . . . . . . . . . . . . 673.5.3 In Vivo Data . . . . . . . . . . . . . . . . . . . . . . . . . 683.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 683.6.1 FEM Results . . . . . . . . . . . . . . . . . . . . . . . . . 693.6.2 Experimental Phantom Results . . . . . . . . . . . . . . . 713.6.3 In Vivo Results . . . . . . . . . . . . . . . . . . . . . . . . 723.6.4 Quantitative Comparison . . . . . . . . . . . . . . . . . . . 793.7 Advantages and Limitations . . . . . . . . . . . . . . . . . . . . . 834 Strain-initialized Bone Surface Detection in 3D Ultrasound . . 854.1 Ultrasound Scanning Protocol . . . . . . . . . . . . . . . . . . . . 854.2 3D Bone Surface Extraction Using Surface-growing . . . . . . . . 864.3 Validation Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . 874.3.1 Simulated Phantom . . . . . . . . . . . . . . . . . . . . . . 874.3.2 Experimental Phantom . . . . . . . . . . . . . . . . . . . . 884.3.3 In Vivo Data . . . . . . . . . . . . . . . . . . . . . . . . . 884.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 884.4.1 FEM Results . . . . . . . . . . . . . . . . . . . . . . . . . 904.4.2 Experimental Phantom Results . . . . . . . . . . . . . . . 904.4.3 In Vivo Results . . . . . . . . . . . . . . . . . . . . . . . . 904.4.4 Quantitative Validation . . . . . . . . . . . . . . . . . . . 944.5 Advantages and Limitations . . . . . . . . . . . . . . . . . . . . . 965 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 975.1 The Impact and Potential Clinical Application . . . . . . . . . . . 975.2 Thesis Contributions . . . . . . . . . . . . . . . . . . . . . . . . . 985.3 Remaining Challenges and Future Work . . . . . . . . . . . . . . 99viiBibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101Appendix A Clinical Evaluation Protocol Form . . . . . . . . . . . 112Appendix B Clinical Evaluation Patient Consent Form . . . . . . 120viiiList of Tables1.1 Acoustic impedances of different body tissues and organs [1]. . . . 183.1 Mean absolute error (MAE) Analysis of the individual effects of us-ing the depth-dependent cumulative power (at ρavg ≈ 0.90), data-driven weight (at ρavg ≈ 0.75 and 0.60), and GMR on the in vivodata. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 803.2 Mean absolute fitting error (mm) comparisons of the 2DSE and2DOPS methods using the in vivo data. . . . . . . . . . . . . . . 803.3 F-Test for the variances in Fig. 3.14. . . . . . . . . . . . . . . . . 814.1 Computation time in second for a volume size of 364×110×50 voxels. 95ixList of Figures1.1 A bird eye view of the ultrasound guided surgical tool navigationsystem in a CAOS procedure. Major system parts are shown assuper-blocks with numbered-dotted-boxes. . . . . . . . . . . . . . 31.2 X-rays are part of the electromagnetic spectrum, with wavelengthsshorter than visible light. Different applications use different partsof the X-ray spectrum [2]. . . . . . . . . . . . . . . . . . . . . . . 51.3 (a) An X-ray CT scanner [3], (b) an X-ray CT scanner with coverremoved to show internal components. Legend- T: X-ray tube, D:X-ray detectors, X: X-ray beam, and R: Gantry rotation [4]. (c) Aschematic image of a CT scanner with X-ray source and detectorsshown in three position [5]. . . . . . . . . . . . . . . . . . . . . . . 61.4 A schematic diagram of multi-slice helical CT acquisition process [6]. 71.5 (a) An illustration of a real MRI scanner [7], and (b) a cutawayschematic image of an MRI scanner [8]. . . . . . . . . . . . . . . . 81.6 (a) A typical fluoroscopic image of the human pelvic region [9],(b) fluoroscopic burn from long exposure [10], and (c) a diagramshowing the components of a fluoroscopic imaging chain [11]. . . . 91.7 (a) A modern medical ultrasound scanner [12], (b) a modern portableultrasound scanner [13], and (c) a schematic diagram of the partsof an ultrasound scanner. . . . . . . . . . . . . . . . . . . . . . . . 111.8 Sound frequency ranges for infrasound, audible (sound) and ul-trasound waves and the corresponding mammals that can hearthem [14]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121.9 A comparison of the resolution and penetration of different ultra-sound transducer frequencies [15]. . . . . . . . . . . . . . . . . . . 13x1.10 The (a) direct and (b) converse piezoelectric effect. Mechanicaldeformation and consequent oscillation caused by an electrical fieldapplied to certain material can produce a sound of high frequency[16]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141.11 (a) Different components of a convex type ultrasound probe, and(b) a graphical representation of the pulse transmission and recep-tion of a piezoelectric crystal [17]. . . . . . . . . . . . . . . . . . . 151.12 Illustration of the (a) function of the backing material, (b) functionof the acoustic matching layer, and (c) function of the acousticlens [17]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161.13 Electronic steering and focusing of ultrasound beam [18]. . . . . . 171.14 Schematic representation of ultrasound pulse generation [19]. . . . 171.15 Different types of ultrasound wave-tissue interactions [1]. . . . . . 191.16 Reflection and refraction of the ultrasound beams [20]. . . . . . . 191.17 Ultrasound image acquisition process [21]. . . . . . . . . . . . . . 211.18 Block diagram showing the B-mode image formation and pre-display processing. . . . . . . . . . . . . . . . . . . . . . . . . . . 211.19 Schematic image showing different tissue layers in a bony anatomy[22]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 241.20 (a) The reverberation artifact in a human ilium scan [23], (b) theshadowing artifact in a human shoulder scan [24], and (c) the Mir-ror image artefact. Persistent haematoma lesion anterior to thetibia also appears deep to the surface of the tibia [25]. . . . . . . . 25xi1.21 Illustration of the anisotropy artifacts in the ultrasound bone scan-ning. (a) A schematic diagram showing the causes of the anisotropyartifact. If the bone surface is parallel to the transducer face,then most of the echo power reflects back to the transducer (bluedotted-arrows), if the bone surface is partially (< 450) inclined tothe transducer face, then a partial echo power reflects back to thetransducer (green dotted-arrows), and if the bone surface is al-most perpendicular (> 450) to the transducer face, then almost noecho power reflects back to the transducer (red dotted-arrows). (b)and (c) show two clinical example images [26] showing anisotropyartifacts. Blue, green and red arrows correspond to regions that re-flects full, partial and almost no echo power from the bone surfaceto the transducer, respectively. . . . . . . . . . . . . . . . . . . . . 251.22 (a) Intra-operative 2D C-arm fluoroscopy, (b) 2D fluoroscopy scanshowing the fractured distal radius and inserted K-wire, (c) 2Dfluoroscopy scan showing the fractured distal radius and metal T-plate used during the surgery. The images are taken during adistal radius surgery at the Vancouver General Hospital (VGH),Vancouver, BC, Canada . . . . . . . . . . . . . . . . . . . . . . . 271.23 Typical isodose curves (in µGy/min) for a mobile C-arm fluo-roscopy system [27]. . . . . . . . . . . . . . . . . . . . . . . . . . . 281.24 Illustration of the 2D bone segmentation performance of the re-flected power imaging-based method developed by Xu et al. [28]using the bone phantom and in vivo data. (a)-(d) are the B-modeimages, and (e)-(h) are the reflected power images. . . . . . . . . 311.25 Illustration of the presence of false positives as well as the falsebone discontinuity in the estimated phase symmetry images usingthe in vivo data. (a)-(d) are the B-mode images with expert delin-eated bone boundaries overlaid. (e)-(h) are the phase symmetryimages correspond to images shown in (a)-(d), respectively. Thefalse positive intensities are shown with yellow arrows while falsebone discontinuities are shown with white arrows. . . . . . . . . . 34xii1.26 Illustration of the 2D bone segmentation performance of methodsdeveloped by Xu et al. [28] using the bone phantom and in vivodata. (a)-(d) are the B-mode images, and (e)-(h) are the strainimages. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 362.1 Illustration of the experimental phantom. (a) The phantom wascreated using a Sawbones radio-opaque hemi-pelvis, cut along redline. (b) The separated portion (shown with red arrow in (a)) ofthe pelvis is embedded in polyvinyl chloride (PVC) gel. . . . . . . 462.2 Illustration of bone boundary detection using the FEM phantom.(a) B-mode image, (b) ideal strain image, (c) strain image gener-ated by the AM method, (d) envelop power map, (e) fused map, (f)maximum intensity points on the fused map, (g) estimated boneboundary by the APS method, (h) estimated bone boundary bythe SEP method, (i) MAE analysis of the SEP and APS methodsat different SNRs, and (j) weight analysis at 10dB SNR (weight atminimum error is shown with red arrow). . . . . . . . . . . . . . . 482.3 Illustration of the bone boundary detection using the experimentalphantom. (a) B-mode image, (b) a CT projection along which USscanning is performed, (c) strain image, (d) maximum intensitypoints overlaid on the fused map, (e) estimated bone boundaryby the APS method, (f) estimated bone boundary by the SEPmethod, and (g) MAE analysis of the SEP and APS methods atdifferent SNRs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 502.4 Illustration of the bone boundary detection using the in vivo volunteer-I data set. (a) shows the scanned region on the volunteer (shownwith red rectangle), (b) is the B-mode image (bone surface is in-dicated with red arrows), (c) is the strain image, (d) shows themaximum intensity points overlaid on the fused map, (e) showsthe detected bone boundary by the APS method, and (f) showsthe detected bone boundary by the SEP method. . . . . . . . . . 52xiii2.5 Illustration of the bone boundary detection using the in vivo volunteer-II data set. (a) shows the scanned region on the volunteer (shownwith red rectangle), (b) is the B-mode image (bone surface is in-dicated with red arrows), (c) is the strain image, (d) shows themaximum intensity points overlaid on the fused map, (e) showsthe detected bone boundary by the APS method, and (f) showsthe detected bone boundary by the SEP method. . . . . . . . . . 532.6 Illustration of the bone boundary detection using the in vivo volunteer-III data set. (a) shows the scanned region on the volunteer (shownwith red rectangle), (b) is the B-mode image (bone surface is in-dicated with red arrows), (c) is the strain image, (d) shows themaximum intensity points overlaid on the fused map, (e) showsthe detected bone boundary by the APS method, and (f) showsthe detected bone boundary by the SEP method. . . . . . . . . . 542.7 Illustration of the bone boundary detection using the in vivo volunteer-IV data set. (a) shows the scanned region on the volunteer (shownwith red rectangle), (b) is the B-mode image (bone surface is in-dicated with red arrows), (c) is the strain image, (d) shows themaximum intensity points overlaid on the fused map, (e) showsthe detected bone boundary by the APS method, and (f) showsthe detected bone boundary by the SEP method. . . . . . . . . . 552.8 Illustration of the bone boundary detection using the in vivo volunteer-V data set. (a) shows the scanned region on the volunteer (shownwith red rectangle), (b) is the B-mode image (bone surface is in-dicated with red arrows), (c) is the strain image, (d) shows themaximum intensity points overlaid on the fused map, (e) showsthe detected bone boundary by the APS method, and (f) showsthe detected bone boundary by the SEP method. . . . . . . . . . 57xiv3.1 (a) Mean absolute error (MAE) analysis for different values of λ at10dB SNR (weight at minimum error is shown with red arrow), (b)Normalized cross-correlation coefficients with respect to appliedstrains for FEM-I phantom (see Section 3.5 for detail), and (c)Normalized cross-correlation coefficient dependent linear weight λ. 633.2 The alignment cost for varying cluster numbers. . . . . . . . . . . 653.3 FEM simulation phantoms. (a) FEM-I: ideal elastogram of a ho-mogeneous soft tissue of stiffness 10kPa. (b) FEM-II: normal distalradius bone and soft tissue region of stiffness 10GPa and 10kPa,respectively, and (c) corresponding ideal elastogram. (d) FEM-III: fractured distal radius bone and soft tissue region of stiffness10GPa and 10kPa, respectively, and (e) corresponding ideal elas-togram. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 673.4 Illustration of 2D bone boundary detection using the FEM-I phan-tom. (a) B-mode image, (b) strain image generated by the pro-posed approach, (c) MEM1, (d) Gaussian mixtures representationover fused map, and estimated bone boundary by the (e) 2DSEand (f) 2DOPS methods (arrows show leakage). Bone boundaryis shown after bottom up ray-casting for the 2DOPS method. . . . 693.5 Illustration of 2D bone boundary detection using the FEM-II phan-tom. (a) B-mode image, (b) strain image generated by the pro-posed approach, (c) MEM1, (d) Gaussian mixtures representationover fused map, and estimated bone boundary by the (e) 2DSEand (f) 2DOPS methods (arrow shows leakage). Bone boundaryis shown after bottom up ray-casting for the 2DOPS method. . . . 703.6 Illustration of 2D bone boundary detection using the experimentalphantom. (a) B-mode image, (b) strain image generated by theproposed approach, (c) MEM1, (d) Gaussian mixtures representa-tion over fused map, (e) estimated bone boundary by the 2DSEmethod, and (f) estimated bone boundary by the 2DOPS method(bone boundary is shown after bottom up ray-casting). . . . . . . 71xv3.7 Illustration of the proposed strain-imaging performance. (a)-(d)show the B-mode images where overlaid red-dotted boxes indicatethe bone locations. (e)-(h) show the strain images generated bythe traditional method adopted in [29] (red arrows indicate thetrue bone locations while yellow arrows indicate artifacts associ-ated with soft tissue layers), and (i)-(l) show the strain imagesgenerated by the proposed method (red arrows indicate the truebone locations). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 723.8 Illustration of 2D bone boundary detection using the in vivo volunteer-I data set. (a) shows bone region scanned on the volunteer (withred rectangle), (b) is the strain image generated by the proposedapproach, (c) is the MEM1, and (d) is the fused map with theGaussian mixtures overlaid. Finally, in (e), orthopaedic expertdelineated bone boundary, and detected bone boundary by the2DOPS and 2DSE methods are overlaid on the B-mode image(bone boundary produced by the 2DOPS method is shown afterbottom up ray-casting). . . . . . . . . . . . . . . . . . . . . . . . 733.9 Illustration of 2D bone boundary detection using the in vivo volunteer-II data set. (a) shows bone region scanned on the volunteer (withred rectangle), (b) is the strain image generated by the proposedapproach, (c) is the MEM1, and (d) is the fused map with theGaussian mixtures overlaid. Finally, in (e), orthopaedic expertdelineated bone boundary, and detected bone boundary by the2DOPS and 2DSE methods are overlaid on the B-mode image(bone boundary produced by the 2DOPS method are shown afterbottom up ray-casting). . . . . . . . . . . . . . . . . . . . . . . . 74xvi3.10 Illustration of 2D bone boundary detection using the in vivo volunteer-III data set. (a) shows bone region scanned on the volunteer (withred rectangle), (b) is the strain image generated by the proposedapproach, (c) is the MEM1, and (d) is the fused map with theGaussian mixtures overlaid. Finally, in (e), orthopaedic expertdelineated bone boundary, and detected bone boundary by the2DOPS and 2DSE methods are overlaid on the B-mode image(bone boundary produced by the 2DOPS method are shown afterbottom up ray-casting). . . . . . . . . . . . . . . . . . . . . . . . 763.11 Illustration of 2D bone boundary detection using the in vivo volunteer-IV data set. (a) shows bone region scanned on the volunteer (withred rectangle), (b) is the strain image generated by the proposedapproach, (c) is the MEM1, and (d) is the fused map with theGaussian mixtures overlaid. Finally, in (e), orthopaedic expertdelineated bone boundary, and detected bone boundary by the2DOPS and 2DSE methods are overlaid on the B-mode image(bone boundary produced by the 2DOPS method are shown afterbottom up ray-casting). . . . . . . . . . . . . . . . . . . . . . . . 773.12 Illustration of 2D bone boundary detection using the in vivo volunteer-V data set. (a) shows bone region scanned on the volunteer (withred rectangle), (b) is the strain image generated by the proposedapproach, (c) is the MEM1, and (d) is the fused map with theGaussian mixtures overlaid. Finally, in (e), orthopaedic expertdelineated bone boundary, and detected bone boundary by the2DOPS and 2DSE methods are overlaid on the B-mode image(bone boundary produced by the 2DOPS method are shown afterbottom up ray-casting). . . . . . . . . . . . . . . . . . . . . . . . 783.13 MAE (mm) comparisons of the 2DSE and 2DOPS methods atdifferent SNR values using (a) the FEM-II, (b) FEM-III, and (c)experimental phantoms. . . . . . . . . . . . . . . . . . . . . . . . 81xvii3.14 Distributions of the estimated bone boundary points by the 2DOPSand 2DSE methods are plotted with respect to the expert delin-eated bone contours for volunteers-I to -V. The ‘mean ± standard-deviation’ of each distribution is shown inside the correspondingplots. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 823.15 Mean absolute fitting error (mm) comparisons of the bone contoursproduced by using the Gaussian mixture regression over the ray-casted points for the 2DOPS method and the contours producedby the original 2DOPS and proposed methods. . . . . . . . . . . . 834.1 Illustration of the FEM results. (a) B-mode image, (b) estimated3D phase symmetry image by the 3DOPS method, (c) bone surfaceafter bottom up ray-casting in (b) (arrows show leakage), and (d)estimated bone surface by the proposed method. . . . . . . . . . . 894.2 Illustration of the experimental phantom results. (a) B-mode im-age, (b) estimated 3D phase symmetry image, (c) bone surfaceafter bottom up ray-casting in (b) (arrows show leakage), and (d)estimated bone surface by the proposed method. . . . . . . . . . . 904.3 Illustration of the 3D bone surface detection using the in vivovolunteer-I data set. (a) shows bone region scanned on the volun-teer (with red rectangle), (b) is the B-mode image, (c) shows theestimated 3D phase symmetry image by the 3DOPS method, (d)shows the bone surface after bottom up ray-casting in (c), and (e)shows the estimated bone surface by the 3DSE method. . . . . . . 914.4 Illustration of the 3D bone surface detection using the in vivovolunteer-I data set. (a) shows bone region scanned on the volun-teer (with red rectangle), (b) is the B-mode image, (c) shows theestimated 3D phase symmetry image by the 3DOPS method, (d)shows the bone surface after bottom up ray-casting in (c), and (e)shows the estimated bone surface by the 3DSE method. . . . . . . 93xviii4.5 Illustration of the 3D bone surface detection using the in vivovolunteer-I data set. (a) shows bone region scanned on the volun-teer (with red rectangle), (b) is the B-mode image, (c) shows theestimated 3D phase symmetry image by the 3DOPS method, (d)shows the bone surface after bottom up ray-casting in (c), and (e)shows the estimated bone surface by the 3DSE method. . . . . . . 944.6 MAE analysis of the 3DSE and 3DOPS methods at different SNRsusing (a) the FEM and (b) experimental phantom data. . . . . . . 954.7 Diagrams showing the computation time breakdown for differentcomponents of the (a) 3DOPS and (b) 3DSE methods. . . . . . . 96xixList of AbbreviationsAbbreviation DefinitionUS UltrasoundRF Radio-frequencyU.S. United StatesCAOS Computer Assisted Orthopedic SurgeryOR Operating Room1D One Dimensional2D Two Dimensional3D Three DimensionalFEM Finite Element ModelMAE Mean Absolute ErrorCT Computed TomographyMRI Magnetic Resonance ImagingMIS Minimally Invasive SurgeryVGH Vancouver General HospitalNSERC Natural Sciences and Engineering Research Council of CanadaRMS Root Mean SquareGMM Gaussian Mixture ModelsGMR Gaussian Mixture RegressionNCC Normalized Cross CorrelationROI Region of InterestSNR Signal-to-Noise RatioMSM Modified Strain MapMEM Modified Envelope MapFM Fused MapPVC Polyvinyl ChlorideCHHM Centre for Hip Health and MobilityGPU Graphics Processing UnitCREB Clinical Research Ethics BoardxxAcknowledgementsAfter the Almighty, it is my great pleasure to express gratitude to the people whomade this thesis possible. Foremost, I would like to express my deepest grati-tude to my supervisors, Prof. Rafeef Abugharbieh, PhD, and Prof. Antony J.Hodgson, PhD, whose expertise, understanding and patience, added significantlyto my graduate experience. This study could have never been done without mysupervisors’ motivation, guidance and inspiration. I would also like to thankProf. Pierre Guy, MD, MBA, who sacrificed his valuable time to demonstrateme the real trauma surgeries in the operating room at the Vancouver generalhospital (VGH), and also to delineate the bone surfaces on the B-mode imagesthat ultimately are considered as the ground truth in our clinical data-basedvalidation.I also like to remember my mentor Prof. Md. Kamrul Hasan of BangladeshUniversity of Engineering and Technology (BUET), who supervised my B.Sc. andM.Sc. theses. He taught me how to work with the medical ultrasound imagingsystem that ultimately helped me cent percent while I worked on this project.I would like to acknowledge the Natural Sciences and Engineering ResearchCouncil of Canada (NSERC) grants no.: 03–4685 and 03–4813 (UBC) that sup-ported this project. I would like to thank the Centre for Hip Health and Mobilityfor providing the lab facilities and the Institute for Computing, Information andCognitive Systems for program support.I would also like to thank my parents, my wife, and my two and half yearsold daughter for their support. Moving to a different continent to pursue anothergraduate degree was a huge step for me, and their encouragement and supportwere invaluable to me as I progressed through this phase of my life. Withoutthem, I would not be the person I am today. Specially, I feel greatly in debt tomy wife who has been raising our child all alone while allowing me to be seventhousand miles apart. This huge sacrifice of her was a great inspiration for me.Finally, I must say that working with brilliant researchers in the BiomedicalSignal and Image Computing Laboratory (BiSICL) at the University of BritishColumbia was my great honor. They provided a friendly environment for me tolearn and grow.xxiDedicationTo my familyxxiiChapter 1Introduction and Background1.1 Motivation and Problem StatementOrthopedic surgery focuses on the diagnosis, treatment, rehabilitation and pre-vention of diseases of the bones, joints, ligament, muscles, tendons and nerves.In these surgeries, orthopedic surgeons deal with various musculoskeletal diseasesand conditions, such as: (a) bone fractures and dislocations; (b) torn ligaments,sprains and strains; (c) tendon injuries, pulled muscles and bursitis; (d) ruptureddisks, sciatica, lower back pain and scoliosis; (e) abnormalities of the fingers andtoes and growth abnormalities; (f) surgical management of degenerative jointdisease; (g) knock knees, bow legs, bunions and hammer toes; (h) arthritis andosteoporosis; (i) bone tumors, muscular dystrophy and cerebral palsy; and, (j)club feet and unequal leg length. In 2005, at least one of these musculoskelatalconditions was reported by 107.67 million adults in the United States (the U.S.),representing nearly one in two persons aged 18 and older [30]. These clinicalconditions often require high risk surgery to improve the quality of the patients’lives after recovery. Moreover, the cost (direct and indirect) of these conditionsis approximately $849 billion annually in the U.S. alone [30].Typically, during an orthopedic surgery, surgical tools and mechanical im-plants are primarily guided or visualized by X-ray-based fluoroscopy (discussedin detail in Section 1.2.2.1). Some surgeries (e.g., bone grafting, knee arthro-plasty, pelvic fracture realignment) take, on average, over one and a half hours1during which the patient and the operating room (OR) staff are exposed to anapproximately 11.2 minutes of fluoroscopic radiation [31]. This exposure timedoes not include the time required to re-position the fluoroscopic C-arm; there-fore, operating the C-arm is one of the most time consuming steps of a fracturerealignment surgery. Minimizing the use of fluoroscopy in these procedures wouldnot only decrease the radiation exposure to the patients and staff but could alsofeasibly decrease the overall time required in an OR.In order to decrease the radiation exposure, surgery time and overall orthope-dic surgical cost, an emerging interest has recently been taking place towards thecomputer assisted orthopedic surgery (CAOS) procedure [32]. The term, CAOS,stands for the approaches that aim to improve the surgical field’s visibility and in-crease the accuracy of the surgical procedures by employing a navigation systemwhen carrying out surgical action. These goals are achieved by linking the bonyanatomy being operated on with a virtual representation, such as an image dataset [32]. Surgical navigation systems that use modern tracking technology are in-troduced and classified according to the chosen virtual representation of the sur-gical object, i.e, image-free and image-based (pre-operative and intra-operative)technology. Within the latter class in particular, computed tomography (CT)-and fluoroscopy-based systems have made their way successfully into the CAOSOR [32].However, since CT scans and fluoroscopy are significant sources of radiation,as we discussed briefly above, there has been an emerging interest in safer, non-ionizing real-time imaging alternatives such as the ultrasound (US) [33]. Gener-ally, an ultrasound guided surgical system has four major parts:1. Pre-operative imaging (US, CT, etc.) and segmentation2. Intra-operative ultrasound imaging and segmentation3. Pre- and intra-operative segmented image registration4. Tracking of the ultrasound probe and surgical tools for the surgical proce-dureThese major system parts are shown in Fig. 1.1 with numbered super blocks(dotted boxes). These numbers correspond to the aforementioned parts listed. In2this system, ‘bone segmentation in ultrasound images’ is one of the steps (shownin Fig. 1.1 with the smaller red block) which is the most crucial, challengingand error prone. Note that the ultrasound can be used as both the pre- andintra-operative imaging modality. If ultrasound is also used as the pre-operativeimaging system, then our thesis work will also include the segmentation block ofSuper block ‘1’ in Fig. 1.1.Figure 1.1: A bird eye view of the ultrasound guided surgical tool navigation system in a CAOSprocedure. Major system parts are shown as super-blocks with numbered-dotted-boxes.1.2 Medical Imaging in Orthopedic SurgeryIn this section, we describe the medical imaging modalities currently being em-ployed pre- and intra-operatively. Pre-operative imaging modalities are typicallyused to diagnose orthopedic clinical conditions and/or in planning surgery proce-dures, while intra-operative imaging modalities are used to assess fractured bonerealignment performance, and to visualize and/or navigate surgical implants,tools, and other associated materials. For the sake of discussion, in this section,we are considering X-ray, X-ray CT, and magnetic resonance imaging (MRI) asthe pre-operative imaging modalities, while fluoroscopy and medical ultrasoundwill be the intra-operative imaging modalities. However, these modalities areoften used interchangeably or solo depending on the requirements and surgical3plans. For example, we have already discussed, in Section 1.1, the fact that theultrasound imaging system can be used for both the pre- and intra-operativepurposes.1.2.1 Pre-operative Imaging1.2.1.1 The X-rayThe most common and widely used medical imaging modality that is typicallyused for orthopedic diagnostic purposes is the X-ray. When a patient is admit-ted to a hospital with a suspected clinical orthopedic condition, standard X-rayimages are first taken. X-radiation (composed of X-rays) is a type of electro-magnetic radiation. Most X-rays have a wavelength ranging from 0.01 to 10nm,corresponding to frequencies in the range of from 30PHz to 30EHz (3×1016Hz to3 × 1019Hz), and energies in the range 100eV to 100keV. X-ray wavelengths areshorter than those of ultraviolate rays and typically longer than those of gammarays. Different applications use different parts of the X-ray spectrum, as shownin Fig. 1.2. In many languages, X-radiation is referred to as Ro¨ntgen radiation,named after Wilhelm Ro¨ntgen, who won the Nobel Prize in Physics as its inventorin 1901 [34].Use of X-ray-based imaging for diagnostic purposes is typically known asradiography. A radiograph is an X-ray image obtained by placing a part of thepatient in front of an X-ray detector and then illuminating it with a short X-raypulse. Bones contain a large amount of calcium, which, due to its relatively highatomic number, efficiently absorbs X-rays. This reduces the number of X-raysreaching the detector in the shadow from the bones, making them clearly visibleon the radiograph. Radiographs are useful in the detection of pathology in theskeletal system as well as for detecting some disease processes in soft tissue. Onenotable example is the chest X-ray, which can be used to identify any orthopedicconditions.X-ray photons carry sufficient energy to ionize atoms and disrupt molecularbonds. This makes the process a type of ionizing radiation, and thus harmfulto living tissue. A very high radiation dose over a short period of time causesradiation sickness, while lower doses can provide an increased risk of radiation-4Figure 1.2: X-rays are part of the electromagnetic spectrum, with wavelengths shorter thanvisible light. Different applications use different parts of the X-ray spectrum [2].induced cancer. In medical imaging this increased cancer risk is generally greatlyoutweighed by the benefits of the examination. The ionizing capability of X-rays can be utilized in cancer treatment to kill malignant cells through radiationtherapy.1.2.1.2 Computed Tomography (CT) ScanX-ray computed tomography is a technology that uses computer-processed X-raysto produce tomographic images of specific areas of a scanned object, allowing theuser to noninvasively view the inside of the object. As the X-ray CT scan is themost common form of CT used in medicine and various other contexts, the term“computed tomography” alone is often used to refer to the X-ray CT scan. Digitalgeometry processing is used to generate a three-dimensional image of the inside ofan object from a large series of two-dimensional radiographic images taken arounda single axis of rotation [35]. Medical imaging is the most common application ofCT. Its cross-sectional images are used for diagnostic and therapeutic purposesin various medical disciplines, especially in the orthopedics.X-ray slice data is generated using an X-ray source that rotates around the5Figure 1.3: (a) An X-ray CT scanner [3], (b) an X-ray CT scanner with cover removed toshow internal components. Legend- T: X-ray tube, D: X-ray detectors, X: X-ray beam, and R:Gantry rotation [4]. (c) A schematic image of a CT scanner with X-ray source and detectorsshown in three position [5].object; X-ray sensors are positioned on the opposite side of the circle from the X-ray source (shown in Fig. 1.3(b)). Initial machines would rotate the X-ray sourceand detectors around a stationary object as shown in Fig. 1.3(c). Following acomplete rotation, the object would be moved along its axis, and the next rotationinitialized. Newer machines permit a continuous rotation with the object to beimaged by sliding it slowly and smoothly through the X-ray ring. These arecalled helical or spiral CT machines. A subsequent development of the helicalCT is the multi-slice (or multi-detector) CT; rather than employing a single rowof detectors, the device uses multiple rows of detectors, simultaneously capturingmultiple cross-sections. A schematic diagram of multi-slice helical CT acquisition6process is shown in Fig. 1.4.Figure 1.4: A schematic diagram of multi-slice helical CT acquisition process [6].The CT is often used to image complex fractures, especially ones surroundingthe joints, because of its ability to reconstruct the area of interest in multipleplanes. Fractures, ligamentous injuries and dislocations can easily be recognizedwith a 0.2mm resolution [36]. However, CT imaging is a major source of radiation.A typical plain film X-ray involves a radiation dose of 0.01 to 0.15mGy, whilea typical CT can involve doses of 10–20mGy for specific organs, and this canreach as high as 80mGy in certain specialized CT scans [37]. One gray (Gy) isequivalent to the absorption of one joule of radiation energy per one kilogram ofmatter. The radiation dose reported in the gray or mGy unit is proportional tothe amount of energy that the irradiated body part is expected to absorb, and thedegree of physical effect (such as the creation of DNA double strand breaks) onthe cells’ chemical bonds by X-ray radiation is proportional to that energy [37].1.2.1.3 Magnetic Resonance Imaging (MRI)Magnetic resonance imaging is a medical imaging technique used in radiologyto investigate the anatomy and physiology of the body during both health anddisease states. MRI scanners employ magnetic fields and radio waves to form im-7Figure 1.5: (a) An illustration of a real MRI scanner [7], and (b) a cutaway schematic imageof an MRI scanner [8].ages of the body. The technique is widely used in hospitals for medical diagnosisand, the staging of diseases and to perform follow-ups without causing exposureto ionizing radiation. In orthopedics, an MRI may be used to examine bones,joints, and such soft tissues as cartilage, muscles, and tendons for injuries, thepresence of structural abnormalities or certain other conditions, such as tumors,inflammatory diseases, congenital abnormalities, osteonecrosis, bone marrow dis-eases, and herniation or the degeneration of discs of the spinal cord. The MRImay be used to assess the results of corrective orthopedic procedures. Joint dete-rioration resulting from arthritis may be monitored by using magnetic resonanceimaging.To perform a study, the patient is positioned within an MRI scanner whichforms a strong magnetic field around the area to be imaged. In most medicalapplications, protons (hydrogen atoms) in tissues containing water molecules areused to create a signal that is processed to form an image of the body. First,energy from an oscillating magnetic field is temporarily applied to the patient atthe appropriate resonance frequency. The excited hydrogen atoms emit a radiofrequency signal which is measured by a receiving coil. The radio signal can bemade to encode position information by varying the main magnetic field usinggradient coils. As these coils are rapidly switched on and off, they create thecharacteristic repetitive noise of the MRI scan. The contrast between different8tissues is determined by the rate at which the excited atoms return to theirequilibrium state. A real life MRI scanner is shown in Fig. 1.5(a), while acutaway schematic image of an MRI scanner is shown with detailed legends inFig. 1.5(b).Although an MRI system is completely radiation free, due to the use in thedevice of a strong magnet, special precautions must be taken to perform MRIscans on patients processing certain implanted devices such as pacemakers orcochlear implants. The MRI technologist needs to obtain information regardingthe implanted device, such as its make and model number, to determine if it issafe to perform an MRI. Patients with internal metal objects, such as surgicalclips, plates, screws or wire mesh, might not be eligible for MRI examinations.1.2.2 Intra-operative Imaging1.2.2.1 FluoroscopyFigure 1.6: (a) A typical fluoroscopic image of the human pelvic region [9], (b) fluoroscopicburn from long exposure [10], and (c) a diagram showing the components of a fluoroscopicimaging chain [11].Fluoroscopy is a method providing real-time X-ray imaging that is especiallyuseful for guiding a variety of orthopedic diagnostic and interventional proce-dures. The ability of fluoroscopy to display motion is provided by a continuousseries of images produced at a maximum rate of 25–30 complete images per sec-9ond. This is similar to the means by which conventional television or video trans-mits images. The acquired images are transmitted to a digital monitor such thatthe body part and its motion can be seen in detail. Fluoroscopy, as an imagingtool, enables physicians to observe numerous body systems, including the skele-tal, digestive, urinary, respiratory, and reproductive systems. Fluoroscopy maybe performed to evaluate the immediate condition of bones, muscles, and joints.A typical fluoroscopic image of the human pelvic region is shown in Fig. 1.6(a).Fluoroscopy and radiography share some of the same imaging chain compo-nents, but difference between them exist. The primary difference is that theradiation exposure rate is much lower for fluoroscopy as compared with that ofradiography. A fluoroscopy of an average-sized adult abdomen typically requiresapproximately 45mGy/min [11]. For an abdominal radiograph, the entrance skinexposure to the patient is approximately 3mGy with an exposure time of 200msecfor an exposure rate of 900mGy/min, which is 20 times higher than the rate fora fluoroscopy. However, the total exposure for a radiograph is much lower thana typical fluoroscopic examination because the fluoroscopic exposure time is ex-tended [11]. For 10 minutes of abdominal fluoroscopy, the total patient exposureis approximately 450mGy, as compared with 3mGy for a radiograph. Long flu-oroscopic exposure causes skin and muscle burn; an example of this is shown inFig. 1.6(b).1.2.2.2 UltrasonographyUltrasound imaging also known as ultrasonography or simply sonography, useshigh-frequency sound waves to view the internal regions of the body. It is the mosteconomical and widely available medical imaging modality. Because ultrasoundimages are captured in real-time, they can also show the movements of the body’sinternal organs as well as the blood’s passage through the blood vessels. Unlike X-ray imaging, no ionizing radiation exposure is associated with ultrasound imaging.Two modern medical ultrasound scanners are shown in Figs. 1.7(a) and (b).Typical diagnostic sonographic scanners operate in the frequency range of 1 to18MHz, although frequencies of up to 50–100MHz have been used experimentallyin a technique known as biomicroscopy in special regions of the body, includingthe eye’s anterior chamber [38].10Figure 1.7: (a) A modern medical ultrasound scanner [12], (b) a modern portable ultrasoundscanner [13], and (c) a schematic diagram of the parts of an ultrasound scanner.A schematic diagram of the parts of an ultrasound scanner is shown in Fig.1.7(c). In an ultrasound examination, a transducer both sends sound wavesand receives their echoing waves. When the transducer is pressed against theskin, it directs small pulses of inaudible, high-frequency sound waves into thebody. As the sound waves bounce off internal organs, fluids and tissues, thesensitive microphone in the transducer records tiny changes in the sound’s pitchand direction. These signature waves are instantly measured and displayed by acomputer, which in turn creates a real-time picture of the information found onthe monitor. One or more frames of moving pictures are typically captured asstill images. Small loops of moving “real time” images may also be saved.Since this thesis work is based on the medical ultrasound, technical detailsof ultrasound data acquisition and its subsequent processing are described inSection 1.3.1.3 The Basic Physics of Ultrasound ImagingUltrasound has been used to visualize the living body’s international recesses forover half a century. An Austrian neurologist, Dr. Karl Theo Dussik was thefirst to apply ultrasound as a medical diagnostic tool to image the brain [39]. At11present, ultrasound is one of the most widely used imaging technologies found inmedicine. It is portable, free of radiation risk, and relatively inexpensive whencompared with other imaging modalities, such as the MRI and CT scans. Fur-thermore, ultrasound images are tomographic, i.e., they offer a “cross-sectional”view of anatomical structures. Here, we describe some of the fundamental prin-ciples and physics underlying ultrasound technology.1.3.1 Sound PropertiesFigure 1.8: Sound frequency ranges for infrasound, audible (sound) and ultrasound waves andthe corresponding mammals that can hear them [14].Sound is a vibration that is transmitted through a medium (e.g., air) that canbe detected by the human ear. Sound specifically refers to such information as thehuman ear is able to perceive, or “hear”. All vibrations, including sound, have afrequency. Frequency is a measure of the number of times something “vibrates”per second. The unit of frequency is Hertz (Hz), interpreted as “vibrations persecond”. The human ear can hear frequencies of between about 20 to 20,000Hz.Therefore, “sound” refers to vibrations within this frequency range. Any vibra-tion below 20Hz is called infrasound, while those beyond 20,000Hz are known asultrasound (US). Frequency ranges for infrasound, sound and ultrasound wavesare shown in Fig. 1.8. The US wavelengths and frequencies are inversely related,i.e., high frequency US has a short wavelength, and vice versa. Medical ultra-12sound devices use sound waves in the range of 1–20MHz. The proper selectionof transducer frequency is an important criterium for providing optimal imageresolution in diagnostic and procedural ultrasound. High-frequency ultrasoundwaves (of short wavelength) generate images of high axial resolution. Increas-ing the number of waves of compression and rarefaction for a given distance canhelp to more accurately discriminate between two separate structures along theaxial plane of wave propagation. However, at a given distance, high-frequencywaves are more attenuated than those of lower frequency; thus, they are mainlysuitable for imaging superficial structures [40]. Conversely, low-frequency waves(of long wavelength) offer images which are of lower resolution but can penetrateto deeper structures due to their lower degree of attenuation. A comparison ofthe resolution and penetration of different ultrasound transducer frequencies isshown in Fig. 1.9.Figure 1.9: A comparison of the resolution and penetration of different ultrasound transducerfrequencies [15].1.3.2 The Generation and Reception of Ultrasound PulsesUltrasound transducers (or probes) contain multiple numbers of a special typeof crystal known as piezoelectric crystals, which are interconnected electronicallyand vibrate in response to the potential applied across them. This phenomenon is13referred to as the piezoelectric effect. In 1880, the French physicists Pierre Curie,and his elder brother Jacques Curie, discovered the piezoelectric effect in certaincrystals [15]. They also demonstrated the reverse piezoelectric effect, i.e., theapplication of electrical current to the quartz resulting in quartz vibrations [41].Latter, Paul Langevin, a student of Pierre Curie, developed piezoelectric mate-rials, which can generate and receive mechanical vibrations with high frequency,and are therefore ultrasonic. These vibrating mechanical sound waves create al-ternating areas of compression and rarefaction when propagated through bodytissues.Figure 1.10: The (a) direct and (b) converse piezoelectric effect. Mechanical deformation andconsequent oscillation caused by an electrical field applied to certain material can produce asound of high frequency [16].The piezoelectric effect (or, direct piezoelectric effect) is a phenomenon ex-hibited by the generation of an electric charge in response to a mechanical force(squeeze or stretch) applied on certain materials (see Fig. 1.10(a)). Conversely,mechanical deformation can be produced when an electric field is applied tosuch material, also known as the converse piezoelectric effect (see Fig. 1.10(b)).Both natural and human-made materials, including quartz crystals and ceramicmaterials, can demonstrate piezoelectric properties. By stacking piezoelectricelements into layers in a transducer (see Fig. 1.11(a)), the latter can convertelectric energy into mechanical oscillations in pulse transmissions; these mechan-ical oscillations are then converted into electrical energy in echo reception (see14Figure 1.11: (a) Different components of a convex type ultrasound probe, and (b) a graphicalrepresentation of the pulse transmission and reception of a piezoelectric crystal [17].Fig. 1.11(b)). The stacking of the piezoelectric elements are called transducerarrays or elements (shown as blue array components in Fig. 1.11(a)).Some other essential components of the pulse generation system are the back-ing material, acoustic matching layer, and acoustic lens. The backing materialis located behind the piezoelectric element to prevent excessive vibration (seeFigs. 1.11(a) and 1.12(a)). It reduces the latter so that the element can generateultrasonic waves with shorter pulse lengths, improving axial resolution in images.On the other hand, ultrasonic waves transmitted from the piezoelectric elementare reflected off a target because there is a big difference in acoustic impedancebetween the piezoelectric element and the object. To avoid this phenomenon, an15Figure 1.12: Illustration of the (a) function of the backing material, (b) function of the acousticmatching layer, and (c) function of the acoustic lens [17].intermediate material called the acoustic matching layer is inserted between thetwo (see Figs. 1.11(a) and 1.12(b)) so that ultrasonic waves can efficiently enterthe object. The acoustic lens is a rubber like material attached to the tip of theprobe (shown in green in Fig. 1.11(a) and in grey in Fig. 1.12(c)). Ultrasonicwaves transmitted from the probe would spread and travel like light. The acous-tic lens prevents the ultrasonic waves from spreading and focuses them in theslice direction to improve the image resolution. Further focusing and steering ofultrasound beam can be done by using delays in the pulsing of elements (see Fig.1.13). This type of arrangement is called a phased or steered array.Ultrasound waves are generated in pulses (intermittent trains of pressure) thatcommonly consist of two or three sound cycles of the same frequency (see Fig.1.14). The pulse repetition frequency (PRF) is the number of pulses emitted by16Figure 1.13: Electronic steering and focusing of ultrasound beam [18].the transducer per unit of time. Ultrasound waves must be emitted in pulses withsufficient time between them to allow the signal to reach the target of interestand then be reflected back to the transducer as an echo before the next pulse isgenerated. The PRF for medical imaging devices ranges from 1 to 10kHz [19].Figure 1.14: Schematic representation of ultrasound pulse generation [19].17Body Tissue Acoustic Impedance (106Rayls)Air 0.0004Lung 0.18Fat 1.34Liver 1.65Blood 1.65Kidney 1.63Muscle 1.71Bone 7.8Table 1.1: Acoustic impedances of different body tissues and organs [1].1.3.3 Ultrasound-tissue InteractionAs ultrasound waves travel through tissues, they are partly transmitted to deeperstructures, partly reflected back to the transducer as echoes, partly scattered,and partly transformed into heat [40]. The echoes reflected back to the trans-ducer produce the image. The amount of echo returned after interacting witha tissue particle is determined by a tissue property called acoustic impedance.This is the intrinsic physical property of a medium defined as the density ofthe medium times the velocity of the ultrasound wave propagation through themedium. Ultrasound wave propagation speed is the speed at which sound cantravel through a medium and is typically 1540m/sec for soft tissue. The speed isdetermined solely by the medium’s characteristics, particularly those of densityand stiffness [18]. Organs containing air (e.g., the lungs) have the lowest acousticimpedance, while dense organs (e.g., bone) have very high-acoustic impedance(see Table 1.1). The intensity of a reflected echo is proportional to the differencein acoustic impedance between the two mediums [19]. If two tissues have iden-tical acoustic impedance, no echo is generated. Interfaces between soft tissuesof similar acoustic impedance usually generate low-intensity echoes. Conversely,interfaces between soft tissue and bone generate very strong echoes due to a largeacoustic impedance gradient [42].When an incident pulse encounters a large and smooth interface of two differ-ent types of tissues possessing different degrees of acoustic impedance, the type18Figure 1.15: Different types of ultrasound wave-tissue interactions [1].Figure 1.16: Reflection and refraction of the ultrasound beams [20].of reflection from this interface is called “specular reflection”, and the echo inten-sity generated is proportional to the difference between their acoustic impedancevalues (see Fig. 1.15(a)). A soft tissue-needle interface at the point where aneedle is inserted “in-plane” (where the needle and transducer faces are parallel)is a good example of specular reflection. If the incident ultrasound beam reachesthe linear interface at 900, almost all of the generated echo will travel back tothe transducer. However, if the angle of incidence with the specular boundary isless than 900, the echo will not return to the transducer, but rather be reflectedat an angle equal to the angle of incidence (i.e., 6 Θi = 6 Θr in Fig. 1.16). The19returning echo will potentially miss the transducer and not be detected [19].If the ultrasound pulse encounters reflectors whose dimensions are smallerthan the ultrasound wavelength, or when the pulse encounters a rough, irregulartissue interface, scattering occurs (see Figs. 1.15(b) and (c)). In this case, echoesreflected through a wide range of angles result in a reduction in echo intensity.However, some echoes return to the transducer regardless of the angle of theincident pulse. Most biologic tissues are filled with tiny scattering structures.The speckle signal that provides organs’ visible texture is a result of scatteredechoes produced within the volume of the incident ultrasound pulse [19].As ultrasound pulses travel through tissue along their depth, their intensity isreduced or attenuated. This attenuation is the result of reflection and scattering,and also of friction-like losses [40]. These losses result from the induced oscillatorytissue motion produced by the pulse, which causes a conversion of energy fromthe original mechanical form into heat. This energy loss to localized heatingis referred to as absorption and is the most important contributor to ultrasoundattenuation [19]. Longer path length and higher frequency waves result in greaterattenuation. Attenuation also varies among body tissues, with the highest degreefor any given frequency found in bone, a lesser degree in muscle and solid organs,and the lowest in blood.1.3.4 Ultrasound BeamformingUltrasound receiving beamforming or image formation is usually performed bysequential insonification of the medium using focused beams. Each focused beamallows the reconstruction of one image line (known as the scan-line or A-line). Atypical 2D image is made of a good number of A-lines (e.g., 64∼512). There area number of advanced beamforming techniques available, including synthetic re-ceive aperture, parallel beamforming, spatial compounding, adaptive beamform-ing technique [43]. However, a simple image acquisition process is illustrated inFig. 1.17. Typically, a subset of the total transducer elements transmits the pulseand receives the reflected echoes; then a delay-and-sum operation is performedover the received pulse sequences and forms the beamformed A-line. This opera-tion repeats by sliding the subset elements along the transducer arrays, with an20overlap from one end to other of the array.Figure 1.17: Ultrasound image acquisition process [21].1.3.5 Ultrasound B-mode Image FormationFigure 1.18: Block diagram showing the B-mode image formation and pre-display processing.21There are many different ways a ultrasound probe can make observations.These methods are called “modes”. Several modes ultrasound employs are uti-lized in medical imaging [44], including the A-, B-, C-, and M-modes, in additionto Doppler mode, pulse inversion mode, and harmonic mode. Since this thesiswork covers B-mode imaging, we discuss only this mode here. In the B-mode (or,brightness mode) ultrasound, a linear array of transducers simultaneously scansa plane through the body that can be viewed on a screen as a two-dimensionalimage. More commonly known as 2D mode now. We show a block diagramthat illustrates the basic B-mode image formation as well as pre-display imageprocessing procedures in Fig. 1.18.After beamforming the ultrasound radio-frequency (RF) signals are digitallyfiltered to remove the high-frequency noise resulting from the digital instruments.Then, the Hilbert transform is used [45] to detect the envelope from the RF signal.Since the envelope signal may contain very sporadic higher intensity points, it issubsequently log-compressed to decrease the dynamic range of the image. Thisin turn improves the visualization and quantification of the anatomical struc-tures in the image [29]. It is further processed by histogram equalization and/orspeckle reduction procedures. The B-mode ultrasound images usually have agrainy look, unique to ultrasound images. The grainy pattern arises from themultiple scattering of the waves originating in the scattereres in the tissue andtheir constructive-destructive interference. These grains are called “speckles” inthe jargon of the field. To enhance the quality of the B-mode images, image pro-cessing algorithms are used to reduce the number of speckles. The last and veryimportant stage of B-mode image formation is scan conversion. From a clinicalperspective, it is very important that the geometry of the B-mode image showthe exact geometry of the tissue being imaged. Since an image does not takeinto account the geometry of the transducer or the manner in which the RF datawas acquired, the generated image, although containing information about thestructure of the tissue in question, is not geometrically correct. Scan conversionis the process by which the data is mapped to the actual geometry of the tissue.221.4 Challenges in Ultrasound Image Interpreta-tionThe manual and automatic interpretation of an ultrasound image may be chal-lenged by the number of artifacts present in image. Artifacts are errors in theimages. They are normally caused by physical processes that affect the ultrasoundbeam and that in someway alter the basic assumptions the operator makes aboutthe beam. To understand artifacts, we need to consider the basic assumptionsthat are made in producing an ultrasound image:1. Sound waves travel in straight lines.2. Reflections are transmitted from structures along the central axis of thebeam.3. The intensity of a reflection corresponds to the reflector’s scattering strength.4. Sound travels at a constant speed of 1540m/sec.5. Sound travels directly to the reflector and back.These assumptions do not always hold true. There are numerous cases involvingexceptions to these assumptions. Although some of these artifacts may actuallyprovide useful information or allow for novel interpretations, the majority arepotential pitfalls that may result in confusion in the manual and automatic in-terpretation of the data. Since this thesis work focuses on bone imaging, we onlydiscuss some major artifacts that are usually encountered in association withwell-known speckle noise during ultrasound bone scanning. However, prior tothat discussion, we show a schematic image showing different tissue layers in abony anatomy in Fig. 1.19. Bone is the strongest reflector of US pulses amongvarious tissue types through which US pulses usually propagate. A typical bonesegmentation algorithm for US images targets to catch the US echo of higherintensity resulting from the periosteum-bone layer. Thus, a segmented bone sur-face actually represents the outer surface of the bone, right below the periosteumlayer.23Figure 1.19: Schematic image showing different tissue layers in a bony anatomy [22].1.4.1 Speckle NoiseSpeckle is a granular “noise” that inherently exists in and degrades the qualityof medical ultrasound images. The natural tissue surfaces are extremely roughon the wavelength scale. Images obtained from these surfaces by coherent imag-ing systems such as ultrasound suffer from a common phenomena called speckle.Speckle is primarily due to the interference of the returning wave at the transduceraperture. The origin of this noise is seen if we model our reflectivity function asan array of scatterers. Because of the finite resolution, at any given time, we arereceiving from a distribution of scatterers within the resolution cell. These scat-tered signals add coherently; that is, they add constructively and destructivelydepending on the relative phases of each scattered waveform. Speckle noise re-sults from these patterns of constructive and destructive interference, shown asbright and dark dots in the image. All the ultrasound images in Figs. 1.20(a)-(c) and 1.21(b), (c) contain common speckle noise with other specific types ofartifacts.1.4.2 ReverberationReverberation is an artifact that appears as multiple, equally spaced lines dis-tributed along a scan-line. Reverberation is caused by the sound bouncing backand forth between tissue boundaries and then returning to the receiver. A B-modeimage of a human ilium scan with associated reverberation artifacts is shown inFig 1.20(a).24Figure 1.20: (a) The reverberation artifact in a human ilium scan [23], (b) the shadowing artifactin a human shoulder scan [24], and (c) the Mirror image artefact. Persistent haematoma lesionanterior to the tibia also appears deep to the surface of the tibia [25].Figure 1.21: Illustration of the anisotropy artifacts in the ultrasound bone scanning. (a) Aschematic diagram showing the causes of the anisotropy artifact. If the bone surface is parallelto the transducer face, then most of the echo power reflects back to the transducer (blue dotted-arrows), if the bone surface is partially (< 450) inclined to the transducer face, then a partialecho power reflects back to the transducer (green dotted-arrows), and if the bone surface isalmost perpendicular (> 450) to the transducer face, then almost no echo power reflects backto the transducer (red dotted-arrows). (b) and (c) show two clinical example images [26]showing anisotropy artifacts. Blue, green and red arrows correspond to regions that reflectsfull, partial and almost no echo power from the bone surface to the transducer, respectively.251.4.3 Attenuation/ShadowingTissues deeper than strongly attenuating objects, such as bone or calcification,appear darker since most of the beam power reflects back from the hard surface,resulting in a lower intensity of the reflected beam beneath the surface of theobject. In the scan of the human shoulder in Fig. 1.20(b), the bottom side showsdecreased beam intensity because of attenuation in the bone surface.1.4.4 Mirror ImagesSound can bounce off a strong, smooth reflector such as a tendon. The surfaceacts as a mirror and reflects the pulse to another tissue interface. The ultrasoundsystem believes the second interface is beyond the first surface, and this is whereit appears on the scan. In Fig. 1.20(c), a B-mode image of a human tibia with apersistent haematoma lesion is shown with associated mirror artifacts.1.4.5 AnisotropyAnisotropy is a well-known property in medical ultrasound imaging describinga different resulting echogenicity of tissues, such as tendons, when the angle ofthe transducer is changed. Tendon fibers appear hyperechoic (bright) when thetransducer face is parallel with the tendon, but can appear hypoechoic (darker)when the transducer face is close to perpendicular. This can be a source ofinterpretation error for inexperienced practitioners. Two clinical sample images,showing anisotropy artifacts, are shown in Figs. 1.21(b) and (c). Blue, green andred arrows correspond to regions that reflect full, partial and almost nonexistentecho power from the bone surface to the transducer, respectively.1.5 Ultrasound-based CAOS SystemsNavigation systems in orthopaedic surgical procedures remain reliant on intra-operative fluoroscopy for the localizing of bone fragments and the guidance ofsurgical tools because of their relatively clear depiction of bone surfaces, frac-tures, implants and surgical tools [46]. However, with such 2D projection imag-26ing, surgeons typically encounter considerable difficulties in accurately localizingbone fragments in 3D space, and in assessing the adequacy and accuracy of re-duced fractures. Some sample images employing 2D fluoroscopy are shown inFig. 1.22. Surgical procedures in this case are thus highly dependent on the ex-Figure 1.22: (a) Intra-operative 2D C-arm fluoroscopy, (b) 2D fluoroscopy scan showing thefractured distal radius and inserted K-wire, (c) 2D fluoroscopy scan showing the fractured distalradius and metal T-plate used during the surgery. The images are taken during a distal radiussurgery at the Vancouver General Hospital (VGH), Vancouver, BC, Canadaperience of the surgeon and are prone to trial and error. Recently introduced 3Dfluoroscopy [47,48,49] provides 3D information about the anatomical area but isnearly twice as expensive as conventional 2D fluoroscopy. Also, in extremely obeseindividuals, image quality can be suboptimal for navigation purposes. Further-more the accuracy of 3D fluoroscopy depends on the rigid relationship betweenthe reference arc and the navigated anatomy. Current radiopaque retractors mustbe removed before image acquisition. Replacing the retractor could affect the ac-curacy of the navigation system [50]. Nonetheless, a recent study [51] reportsthe 92% accuracy of 6,617 screws insertions in 1,105 patients in a fluoroscopyguided surgery procedure. However, fluoroscopy involves significant radiation ex-posure which is potentially harmful to both patients and surgical teams, withthe latter enduring repeated exposure on a regular basis. We present an isodosecurve showing the radiation pattern for a mobile C-arm fluoroscopy system inFig. 1.23. Here we see that a significant amount of radiation exposure is presentat around a 100cm distance from a fluoroscopic X-ray tube, within the range ofwhich, surgical staff usually work.Recently, ultrasound (US) has been proposed as a safer, more promising al-ternative [52,53,54,55], and has been shown to be very effective in intra-operative27Figure 1.23: Typical isodose curves (in µGy/min) for a mobile C-arm fluoroscopy system [27].navigation. However, US can be used in two ways:1. Using US only intra-operatively for navigation procedures, or along withfluoroscopy, which reduces the use of the fluoroscopy, thus reducing theradiation exposure.2. Using US pre-operatively in diagnosis and surgical planning, and intra-operatively for navigation alone, such that there is no radiation exposure.To use US for surgical tool navigation procedures, intra-operative bone segmenta-tion in US images needs to be (a) automatic, (b) accurate, and (c) fast. However,as we have already discussed in the previous section, US is notorious for being avery noisy modality that is comparatively difficult to interpret, as well as beingmore suitable for the imaging of soft tissue interfaces. Hence, the method cannot28so reliably interpret data from such hard tissue localities as bone surfaces. Re-cent works on using phase-based features for US bone segmentations [56,57,58,59]have been most promising to date but, nonetheless, remain somewhat conduciveto false positives where soft tissue interfaces can generate phase-based responsesthat are extremely similar to those at bone surfaces and can thus be misclassi-fied [29].1.6 Current Ultrasound-based Bone Segmenta-tion Methods1.6.1 B-mode Imaging-based MethodsThroughout the last two decades, a number of manual, semi-automatic and au-tomatic methods for bone segmentation in US have been developed based onB-mode image information. These methods can be sub-categorized as image in-tensity and local gradient-based, a-priori -bone-appearance-based, active-shape-based, and image-phase-feature-based methods. In this section, we provide abrief overview on some of these methods, their achievement, and limitations.1.6.1.1 Image Intensity-based Manual MethodsIn the late 1990s and early 2000s, manual identification was the most commonimage intensity-based approach for identifying bone surface points from US data.Tonetti et al. [60] manually selected bone contour points from B-mode images ina percutaneous iliosacral-screwing surgery which achieved sufficient accuracy forthe placement of screws in a clinical validation procedure. The degree of overallerror resulting from the calibration, segmentation and registration was 2.57mm.Nevertheless, the manual segmentation was excessively time-consuming, resultingin an overall surgery time of ∼50min. In a human cadaver study, Barratt etal. [61], manually segmented pelvis and femur bone surfaces from B-mode imagesand registered them to corresponding CT data sets, where they reported anaverage root mean square (RMS) target registration error of 1.6mm. Here also,the manual segmentation process was excessively time-consuming, resulting in29an overall surgery time of ∼90min. Beek et al. [62, 63] also employed manualbone segmentation and, in one of their studies [63], used 5 manually selectedseed pixels from the B-mode image of a scaphoid phantom bone to obtain a pointcloud which was then used during the registration algorithm. The reported CTto US registration error was 0.54mm (max error: 0.68mm; standard deviation:0.07mm). However, the authors noted their segmentation non-real-time as being∼3min. Therefore, despite manual segmentation approaches obtaining a sub-millimetric accuracy, a major limitation of these methods was the time requiredfor manual segmentation.1.6.1.2 Image Intensity-based Automatic MethodsA number of groups attempted to automate the US image segmentation pro-cess for various applications [64]. Several traditional automatic methods usedultrasound image intensity and local gradient information for bone segmenta-tion [65, 66, 67]. Kowal et al. [66] used depth weighted thresholding, followed bymorphological operations and connected component labeling, in their bone seg-mentation framework. Validation on the porcine and bovine specimens revealeda mean accuracy of 0.42mm and a 0.8s average processing time. However, theauthors of this article used very well-chosen US image sets that had more promi-nent bone intensities than those of the soft tissue regions. Thus, a simple B-modeimage thresholding operation was able to extract the bone surface. However, therobustness of their method was not tested on challenging clinical image sets wherereverberation and anisotropy artifacts are prominent. Daanen et al. [65] devel-oped an automatic segmentation method for the delineation of bone in US images.Validation on clinical patient data shows a maximum mean error of 8.8 pixels anda minimum mean error of 4.545 pixels, with a pixel size of 0.112mm×0.109mm.The time required to delineate one image was reported to be less than 4s, and adataset of 69 images took them about 4 minutes to perform. However, it remainsdifficult to manage the sensitivity of intensity and gradient-based techniques toUS artifacts, machine settings and algorithm parameters [33].Few dynamic programming-based automatic bone segmentation methods arealso reported in the literature [68,69]. Still the bone detection accuracy is shownto be approximately 3mm with respect to the expert delineated bone surfaces. In30addition, these methods are also reported to be time inefficient, e.g., the methodof [69] which is implemented in C++ requires ∼5mins to extract 3D bone surfacesin a US volume of size 200×100×100 voxels.Figure 1.24: Illustration of the 2D bone segmentation performance of the reflected powerimaging-based method developed by Xu et al. [28] using the bone phantom and in vivo data.(a)-(d) are the B-mode images, and (e)-(h) are the reflected power images.Wen et al. [28] developed a reflected-power-imaging-based bone surface vi-sualization approach for 2D ultrasound images. Their qualitative results on aphysical bone phantom and human hand phalange are illustrated in Fig. 1.24.The reflected-power-imaging-based method depicts a band of region where thebone surface might be found but does not specify any contour for the actualbone (see Figs. 1.24(e)-(h)). Therefore, this method can be used to generate aprobabilistic bone region but is not suitable for specific bone contour localization.In addition, since reflected power imaging depends on the image intensity, thisapproach is prone to signal dropout in places where bone intensity is not muchmore prominent than that of the soft tissue region. This scenario occurs wherethe bone surface is not perpendicular to the ultrasound beam’s direction. A sim-ilar example is shown with yellow arrows in Figs. 1.24(d) and (h). Furthermore,Xu et al. [28] did not include any quantitative validation of their methods.311.6.1.3 A-priori -Bone-Appearance-based MethodsIn order to render the problems of image intensity-based method more tractable,some researchers have attempted to incorporate a priori bone appearance infor-mation into their frameworks [65,68,70]. Note that we are considering probabilis-tic bone appearance- [70] and active-shape-based [71] methods as a-priori -bone-appearance-based methods. These methods mathematically model the bone sur-face region. Combining such models with intensity and gradient information gavepromising results. Daneen et al. [65] reported a maximum mean bone localizationerror of 8.8 pixels (i.e., 0.8mm). However, fractured bone surfaces in orthopaedicsurgery applications, as well as reduced bones secured with internal fixation de-vices, do not have continuously smooth surfaces and often significantly violateprior assumptions regarding bone shape. Methods [71,72] using the active shapemodels also make a priori assumptions of the targeted bone structures, which alsofail for fractured bones. Building reasonable models of all possible fracture sce-narios into the system is not practical. Furthermore, methods using the evolvingcontour (i.e., active shape) encounter difficulties near the fractured regions wherethey begin to “leak” into the soft tissue regions. In addition, intensity-based aswell as active-shape-based bone segmentation methods were prone to a high levelof speckle noise, reverberation, anisotropy, and signal dropout [59,73].1.6.1.4 Multi-modal Registration-based MethodsSome research groups developed alternate methods by combining segmentationtechniques with the multi-modal registration of US to CT images [53, 54, 74, 75].Amin et al. [54] combined information from three perspectives: ultrasound boneintensity, shadowing beneath the bone surface, and a spatial prior obtained byprocessing a CT volume. Rather than explicitly segmenting the bone surfacefrom the US image, these three sources of information were used to obtain anidea of regions which are likely to contain bone surfaces. During the registrationprocess, the regions were refined to select data based on its consistency with the3D shape of the bone. The overall accuracy of the system is directly relatedto the accuracy of the initial registration estimate which also depends on theexperience of the surgeon and the quality of the image. However, they failed to32report any accuracy results, but, rather, reported a registration error of 1.94mmfor average translations and of 0.900 for average rotations. Brendel et al. [74]first segmented 3D bone surfaces from a CT image. Then, a part of the bonesurface, which should be visible in the CT image, was segmented from the USdata. The segmented bone surface in the CT image was subsequently registeredto the US image, which resulted in 0.5mm of surface fitting error. Ionescu etal. [75] employed a similar approach; the main difference was that they alsosegmented the US images using classical image segmentation techniques basedmainly on linear filtering or mathematical morphology. This segmentation wasthen updated using US-CT registration and a reported registration error of 2mmfor average translation and 20 for the average rotations.1.6.1.5 Image-phase-feature-based MethodsIn the last decade, Hacihaliloglu et al. [33,73,59,56,57,76,77] used a log-Gabor fil-tering based ridge detection approach to delineate bone in 2D and 3D US images,while reporting the achievement of sub-millimetric accuracy in bone localizations(with a mean absolute error <1mm) on their in vivo data. These methods exploitthe phase congruency features present in an ultrasound image. The underlyingprinciple is that phase is constant and congruent across all the scales at thelocation of which the human visual system would be perceived locally as a one-dimensional feature such as an edge. It has been demonstrated [76] that phaseinformation can provide more significant information within an image than am-plitude information; therefore, the human brain can perceive more enhanced fea-tures in phase space. Phase congruency [78] permits feature detections indepen-dent of the actual feature type. However, these image-local-phase-feature-basedbone segmentation methods use pre-determined (manually or data-driven) angu-lar orientations of bone in the corresponding US images along which log-Gaborfiltering is performed. Since soft tissue interfaces near the transducer face showssimilar intensity profiles as well as the fact that they lie along similar angularorientations to those of the bone contour, log-Gabor filtering eventually producesconfounding phase-symmetry profiles for both the bone and soft-tissue interfacesthat sometimes cause confusion when deciding on the actual bone contour. Theauthors further attempted to address this issue by adopting bottom-up ray cast-33ing and image cropping, which, however, may not have been a reliable approachfor use in a critical surgery environment. In addition, phase-based methods mayalso fail to extract bone boundaries in places where B-mode intensities are notprominent with respect to the background. These limitations are demonstratedin Fig. 1.25. We show the presence of such false positives (shown with yellowarrows) as well as the discontinuity in the bone boundaries (shown with whitearrows) in Figs. 1.25(e)-(h).(a) (b) (d)B-modePhase-symmetry (c)(e) (f ) (h)(g)In Vivo Data 1 In Vivo Data 2 In Vivo Data 3 In Vivo Data 4Figure 1.25: Illustration of the presence of false positives as well as the false bone discontinuityin the estimated phase symmetry images using the in vivo data. (a)-(d) are the B-modeimages with expert delineated bone boundaries overlaid. (e)-(h) are the phase symmetry imagescorrespond to images shown in (a)-(d), respectively. The false positive intensities are shownwith yellow arrows while false bone discontinuities are shown with white arrows.1.6.1.6 DiscussionFrom this brief discussion on the B-mode imaging-based bone segmentation meth-ods, we can conclude that, although the manual segmentation methods achievedbetter bone localization accuracy, they were too time-consuming to be widelyutilized in the OR. On the other hand, image-intensity-based automatic meth-ods are prone to displaying image artifacts (e.g., reverberation, anisotropy andsignal dropout/shadowing), machine settings and algorithm parameters. Someautomatic intensity-based methods, i.e., dynamic programming-based methods,34though, achieved moderate accuracy in bone localizations but at the cost of re-quiring very high computation time. A-priori -bone-appearance-based segmenta-tion methods showed promise in bone localizations, but failed to predict the shapeof the fractured bones. A few multi-modal registration-based methods for bonesegmentation are also reported, but nonetheless depend on ultrasound image in-tensity information, which is prone to displaying image artifacts similar to thoseof the image-intensity-based methods. Image-phase-feature-based bone segmen-tation methods addressed the limitations of the previous methods and achievedsub-millimetric accuracy in bone localization. These methods are also based onthe B-mode imaging. However, as we already have shown in Fig. 1.25, false pos-itive bone responses are produced by these methods at the soft tissue interfaceswhich sometimes result in confusion in the equipment’s determining the actualbone contour. In addition, the automatic 3D-phase-feature-based method [59]is found to be time inefficient (at ∼2mins), although a graphics-processing-unit(GPU)-based implementation of a version of this method, which requires man-ual selection of all the parameters, takes 0.8s for its execution [79]. Therefore,our thesis has focused on developing an effective method for performing bonesegmentations (i.e., an elastography/strain imaging-based method) that wouldexploit the inherent bone localizing information which is not usually revealedin the typical B-mode images and, thus, is not susceptible to the developmentof image artifacts. In addition, this thesis focuses on developing real-time ornear-real-time bone segmentation methods.1.6.2 Strain-imaging-based MethodsElastography (strain imaging) is an emerging medical diagnostic tool for captur-ing the mechanical properties (e.g. stiffness) of biological tissue [45,80,81,82,83].Elastography has shown promise in the detection of breast and prostate tu-mours [45], liver cirrhosis [80], vascular plaques [83], and ligament-bone inser-tions [84]. Because bones are much stiffer than the surrounding tissue, strainimaging is a good candidate for bone surface visualization. However, the onlymethod that investigated the potential of using elastography in 2D bone segmen-tation in US prior to us was that of Xu et al. [28], who investigated a strain-35Figure 1.26: Illustration of the 2D bone segmentation performance of methods developed byXu et al. [28] using the bone phantom and in vivo data. (a)-(d) are the B-mode images, and(e)-(h) are the strain images.imaging-based technique in [28] for bone surface visualization in 2D ultrasoundimages. Their qualitative results on a physical bone phantom and human handphalange are shown in Fig. 1.26. Their strain-imaging-based bone detection re-sults revealed a low degree of contrast in human data between the hard bone andsoft tissue regions (see Figs. 1.26(g) and (h)), although it exhibited some boneoutlines of bones in the phantom images (see Figs. 1.26(e) and (f)). Based onthe results of this method, a region can be predicted where the bone surface canbe found but where no contour for the actual bone can be specified. The exactbone extraction process from these fuzzy strain images is not further mentionedin the article. A quantitative validation of this method is also not reported.Moreover, problems arose from signal windowing as window-based elastographymethods (as used by Xu et al. [28]) have limitations associated with the sizeof the window segments. A significant amount of noise in the strain image caneasily be introduced with the choice of smaller window sizes and/or large overlapbetween successive windows [83]. In addition, in the window-based methods, theestimated strain for a particular window is considered to correspond to the centreof the window. This might result in a bone localization error of 12×window length(>1mm), which is larger than that of our target (sub-millimetric) accuracy.361.7 Thesis ObjectivesThe context of this thesis work was ‘US guided navigation in CAOS’, whichrequires robust, accurate and rapid bone segmentation in US images. Althoughdiverse, previously reported methods achieved satisfactory accuracy, they oftenfailed to vigorously extract bone surfaces due to the presence of various types ofimage-intensity-based artifacts. Therefore, the objective of this thesis work wasto develop more robust, accurate and automatic bone segmentation techniques forultrasound images than found in the state-of-the-art methods. This research wasintended to further advance the larger goal of developing a novel 3D ultrasound-based CAOS system. The goals of this thesis include the following:• Developing new image processing methods for automatically extractingbone surfaces from 2D and 3D ultrasound images that are substantiallymore robust to false positives than current phase symmetry methods.• Demonstrating comparable or better accuracy in bone localization on invivo data than found in the state-of-the-art methods.• Maintaining or reducing the computation time required for state-of-the-artmethods.We believe that the outcomes of this thesis work would (a) improve the per-formance of the ultrasound guided navigation procedure by providing better as-sessments and placements of the bone location, which might in turn decreasesurgery time; (b) decrease the amount of radiation exposure to patients and/orstaff during orthopedic surgeries; (c) promote minimally invasive surgery (MIS)by minimizing the extent of soft tissue incisions; and (d) decrease the cost andimprove the quality of health care by replacing fluoroscopy at key points in or-thopedic diagnoses and treatments.1.8 Thesis OrganizationThis thesis is separated into three parts. In Chapter 2, we present our preliminarywork [29], where we investigate the potential of combined US strain imaging and37envelope power signaling to delineate bone boundaries in US images. Our methoduses real-time strain imaging [82] based on an analytic minimization of regularizedcost functions to delineate bone from the tissue stiffness map and, consequently,we have achieved a marked reduction in false positive bone responses at the softtissue interfaces. However, our method relies on assumptions based on empiricalobservations.In Chapter 3, we present our subsequent study in which we mitigate the limi-tations of this preliminary work presented in Chapter 2 by introducing automaticparameter selection processes that replace all of the intuitive assumptions andthus, make the method more robust and self-tuning. We improve the methodin [29] by incorporating the depth-dependent cumulative power of the envelopeinto the pre- and post-compression data, which we demonstrate as reducing thesoft-tissue-related artifacts in the resulting strain images. We further improvethe method by incorporating a data driven weight into the technique, estimatedfrom the echo de-correlation measurements performed between the pre- and post-compression RF frames, and which controls the contribution of the strain imageinto the resulting fused map. We also use a local statistics-based bone discon-tinuity detection scheme. Finally, we introduce Gaussian mixture regression, amore flexible multivariate non-parametric regression model that better preservesthe curvature features in the bone boundary than the linear regression modelemployed in our preliminary work.Although our studies presented in Chapters 2 and 3 show great promise,they are limited only to 2D US images, while a near-real time 3D elastographytechnique using a 3D probe has been presented that can generate a 3D strainimage [81]; this technique has not been validated on clinical in vivo data, andwe have found it to be very sensitive to echo de-correlations while acquiring apost-compression US volume which limits its applicability. Therefore, in Chapter4, we present our simple yet effective technique to a 3D bone surface extraction,using a surface-growing approach that is seeded from 2D bone contours estimatedfrom combined ultrasound strain imaging and envelope power. The 2D seed bonecontour grows into the 3D surface by minimizing a combined intensity similarityand a voxel proximity-based cost function.The performance of these algorithms is evaluated using finite-element-model38(FEM) simulation phantoms and, experimental phantoms, as well as in vivodata, and compared with two state-of-the-art image-phase-features-based boneboundary estimation algorithms.39Chapter 2Bone Detection in UltrasoundUsing Combined Strain-imagingand Envelope PowerIn this chapter, we discuss our preliminary approach of bone segmentation bycombining US strain-imaging and envelope power detection at each RF sample.Our method uses real-time strain-imaging [82] based on an analytic minimizationof regularized cost functions to delineate bone from the tissue stiffness map andthus, we achieved a marked reduction of the false positive bone responses at thesoft tissue interfaces.2.1 Ultrasound Signal ModelThe one-dimensional (1D) time domain models of the back-scattered ultrasoundRF signals before and after compression can be written as [83]I1(t) = s(t) ∗ p(t) + ν1(t), (2.1)I2(t) = s(t− t0)∗ p(t) + ν2(t), (2.2)where I1(t) and I2(t) denote the pre- and post-compression RF echo signals,respectively, s(t) denotes the 1D US scattering function, p(t) denotes the pointspread function,  denotes the compression factor caused by mechanical deforming40pressure to the medium, ν1(t) and ν2(t) are the uncorrelated random noise profilesand * sign is used to denote the convolution operation. The post-compressionRF echo is assumed to be the scaled and shifted replica of the pre-compressionRF echo in the absence of noise. According to [45], the strain S is related to thecompression factor 1/ as S = 1−  where  ≤ 1 and S << 1.2.2 Modified Strain Map (MSM)The analytic minimization (AM)-based strain estimation method has been de-scribed in detail in [82]. In this section, we provide a brief overview of thismethod.2.2.1 Displacement Estimation2.2.1.1 Cost FormationTo estimate tissue displacement from the pre- and post-compression RF echoframes I1 and I2, respectively, of size m×n (considering in the discrete domain),the integer axial and lateral displacements of a seed scan-line are estimated us-ing dynamic programming [85]. Using these estimates as a priori information,displacements for all other scan-lines are estimated in the sub-sample level usingthe regularized cost function asCj(∆a1, ..,∆am,∆l1, ..,∆lm) =m∑i=1{[I1(i, j)− I2(i+ ai + ∆ai, j + li +∆li)]2+ α(ai + ∆ai − ai−1 +∆ai−1)2 + βa(li + ∆li − li−1 +∆li−1)2+ β′l(li +∆li − li,j−1)2}, (2.3)where α, βa and β′l are the tuneable parameters, ai and ∆ai are the integer andsub-sample axial displacements, respectively, and li and ∆li are the integer andsub-sample lateral displacements, respectively, of the ith sample in the jth scan-line. Now the post-compression data term in Eqn. (2.3) can be written in the2D Taylor expansion around (i+ ai, j + li) asI2(i+ ai + ∆ai, j + li + ∆li)≈ I2(i+ ai+, j + li) + ∆aiI ′2,a + ∆liI ′2,l, (2.4)41where I ′2,a and I ′2,l are the derivatives of the I2 in the axial and lateral directions,respectively.2.2.1.2 Cost OptimizationNow the optimal solution for (∆ai,∆li) can be estimated from Eqns. (2.3) and(2.4) by setting (∂Cj)/(∂∆ai) = 0 and (∂Cj)/(∂∆li) = 0 for i = 1...m. If4d = [∆a1∆l1∆a2∆l2...∆am∆lm]T and initial estimates d = [a1l1a2l2...amlm]T ,then after partial differentiating Eqn. (2.3) and equalize it to zero, it can bewritten as(I ′22 +D1 +D2)∆d = I ′2e−D1d, (2.5)D1 =α 0 −α 0 0 0 ... 00 βa 0 −βa 0 0 ... 0−α 0 2α 0 −α 0 ... 0: : : : : : : :0 0 0 ... 0 −βa 0 βawhere D2 = diag(0, β′l, 0, β′l, ..., 0, β′l) is a diagonal matrix of size 2m × 2m, e =[e1e1...emem]T , ei = I1(i, j) − I2(i + ai, j + li), I ′22 = diag(J ′2(1)...J ′2(m)) is asymmetric diagonal matrix of size 2m× 2m withJ ′2(i) =[I ′22,a I ′2,aI ′2,lI ′2,aI ′2,l I ′22,l]2.2.1.3 Practical Assumptions IncorporationTo cope with the practical scenarios, some assumptions are made in Eqn. (2.5):1. The first assumption incorporates depth-dependent US attenuation,2. The second assumption ensures very small lateral displacements (i.e. dif-ference between consecutive lateral displacements is up to 4 scan-lines),3. And the third assumption ensures zero displacement at the transducer face.42Moreover, for robust displacement estimates, iteratively reweighted least squares(IRLS) is used. Consequently, Eqn. (2.5) is modified along with facilitatinginverse gradient estimation as(WI ′21 + ZD1 + ZD2)∆d = WI ′1e− ZD1d+ s, (2.6)where W = diag(0, w(r1), w(r1)...w(rm), w(rm)), w(ri) = 11+(ri/T )2 is the IRLSCauchy weight, ri = I1(i, j) − [I2(i + di, j + ai) + ∆diI ′2,a + ∆aiI ′2,l], T is atunable threshold, s is the bias vector to ensure small lateral displacements,Z = diag(ζ1, ζ1, ...ζm, ζm), and ζi = [exp(1540×102atfolog(10)20fs×106 )]−i. Here, at denotesthe frequency dependent attenuation coefficient and is 0.6dB/cm/MHz for thesoft tissue [40], and fo and fs are the transducer operating and sampling frequen-cies (in MHz), respectively. Finally, Eqn. (2.6) is solved to calculate a regularizeddisplacement map ∆d.2.2.2 Strain Estimation and ProcessingA Kalman filter-based method is adopted to estimate the final strain S fromthe estimated displacement map [82]. We use this estimated strain map S as aprimary indicator to know the bone boundary location and expect to see lowerstrain values along the bone boundary. For robust bone boundary estimation,we combine the modified strain map with the envelope power map. Therefore,to produce a strain map with a dynamic range roughly equivalent to that of theenvelope power map, we estimate the MSM asMSM = Sˆmax[Sˆ], (2.7)where Sˆ =| (−S)+ | median(−S) || and median is the median operator.2.3 Modified Envelope Map (MEM)The performance of strain-imaging is prone to echo de-correlation [83]. In aclinical settings, a radiologist, inexperienced in elastography may find it hardto produce a good strain image to visualize bone. Therefore, we intend to fuse43bone appearance information from another reliable source, i.e., envelope of theRF data, into the strain image-based information. An envelope map can begenerated from a single RF frame and hence it is free of echo de-correlation-based artifacts. Moreover, an envelope map has better dynamic range in termsof its intensity contents than a log-compressed envelope (i.e., B-mode) map [68].This higher dynamic range helps better distinguishing between the intensitiescorrespond to the bone and soft tissue regions. We further increase the dynamicrange of a typical envelope map by estimating its power. However, as the strainis mapped with respect to the pre-compression RF frame, we use the same pre-compression RF frame to generate the modified envelope map (MEM) so thatit can be spatially matched to the MSM. The envelope of a RF scan-line atcolumn j is calculated as Ej(i) =| H[I1,j(i)] | [45], where H denotes the Hilberttransform. This transformation converts a periodic signal into its analytic formthat contains amplitude and phase information separately. In easy words, wecan get a demodulated signal (envelope) from its carrier (raw RF) after thistransformation. Finally, we estimate the MEM asMEM = E2max[E2]. (2.8)2.4 Fused Map EstimationOur next step involves fusing the MSM and MEM into a single fused map (FM).We estimate FM to complement bone information from two different sources andtypes of information. Strain image is based on the tissue stiffness informationwhile envelope contains the echo intensity-based information. As we already dis-cussed that the performance of the strain-imaging is prone to echo de-correlation,envelope power map complements with the higher echo intensity information re-sulting from bone surface. On the other hand, we already discussed in Chapter 2that intensity-based bone segmentation methods have limitations resulting fromhigh levels of speckle noise, reverberation and signal drop out, therefore, strainimage complements with the higher stiffness information resulting from the boneregion. The FM is estimated asFM = λ×MSM+ (1− λ)×MEM, (2.9)44where λ is the weight. To choose a suitable value of λ, we analyse the effect ofthis choice on the accuracy of the bone boundary localization in a FEM-basedsimulation (see Section 2.6.1) of a noisy environment (i.e., 10dB signal-to-noiseratio (SNR)). For this MAE analysis, we deliberately chose SNR to be lowerthan the typically encountered values (approximately 18∼22dB) in RF data afterbeamforming [86] to simulate a worse case scenario. We observe that the meanMAE is minimized for λ = 0.5 (see Fig. 2.2(j)). While this value for λ may notbe optimal under all possible imaging situations, the results seem to be relativelyinsensitive to this choice, so we opted to set λ = 0.5 for this study.2.5 Forming Bone Contour Using RegressionIn the fused map, the maximum intensity points are expected to correspondto the actual bone contour since the maximum intensity points in the MSMand MEM individually are likely to represent the bone surface. Therefore, afterestimating the FM, we detect the maximum intensity point along each scan-lineof the FM. Then we estimate the final bone boundary using a local linear-fitover the detected maximum intensity points. The linear fit window length ischosen to be 3mm. The choice of the fitting window length is important sincea small window will not adequately smooth the estimated surface while a largerwindow could fail to adequately represent real curvature in the bone contour. Wehave empirically selected a value of 3mm, as this seemed to work well with ourpilot data; however, in future work, it would be desirable to determine this valueusing a more principled approach, whether by incorporating a priori informationbased on knowledge of the scan target or derived directly from the image beingprocessed.2.6 Validation Setup2.6.1 Simulated PhantomWe built a 40mm×40mm finite-element-model (FEM) phantom using the AN-SYS analysis software (ANSYS, Inc., Canonsburg, PA) and then generated an45ultrasound simulation of the model using Field II [87]. Our phantom mimickedan US scan of the human distal radius bone with a total number of nodes of55,180. The stiffness of the homogeneous soft tissue and bone region were set to10kPa and 10GPa, respectively as previously reported in the literature [88]. Ourphantom was compressed from the top using a larger-width planar compressor infree-hand fashion. An ultrasonic transducer of center frequency, f0 = 5MHz andband-width = 50% was used to simulate the phantom scan from the top. Thetotal number of scan lines was set to 128. We set an applied pressure level thatcorresponds to 1% average strain. In addition, we did not model out-of-planemotion.(a) (b)Figure 2.1: Illustration of the experimental phantom. (a) The phantom was created using aSawbones radio-opaque hemi-pelvis, cut along red line. (b) The separated portion (shown withred arrow in (a)) of the pelvis is embedded in polyvinyl chloride (PVC) gel.2.6.2 Experimental Physical PhantomWe constructed an experimental phantom using a radio-opaque Sawbones hemi-pelvis (Sawbones, Pacific Research Laboratories, Inc., Vashon Island, WA), modelnumber 1297-22 (see Fig. 2.1). A portion of the pelvis was suspended in polyvinylchloride (PVC) gel and placed in an acrylic tube. A high resolution periph-eral quantitative computed tomography (CT) machine, model HR pQCT XtremeCT (Scanco USA, Inc., Wayne, PA) was used to acquire a single 482×482×402(lateral×axial×elevational) voxel volume. The resolution of the acquired volumewas 0.25mm×0.25mm×0.25mm. In addition, US was acquired using a SonixRP46(Ultrasonix Medical Corporation, Richmond, BC) scanner integrated with a L14-5W/60 probe operating at 10MHz in the Centre for Hip Health and Mobility,Vancouver, BC, Canada.2.6.3 In Vivo DataWe acquired five sets of data from five volunteers (volunteer-I: 25-year old male;volunteer-II: 26-year old male; volunteer-III: 24-year old male, volunteer-IV: 25-year old male; volunteer-V: 33-year old male) after obtaining prior consent. Alldata were acquired with free-hand compression. US was acquired using a SonixRP(Ultrasonix Medical Corporation, Richmond, BC) scanner integrated with a L14-5W/60 probe operating at 10MHz in the Centre for Hip Health and Mobility,Vancouver, BC, Canada. The study was approved by the UBC clinical researchethics board.2.7 ResultsWe provide comparative results of our proposed strain and envelope power basedmethod (SEP) with the automatic adaptive parameterization in local phase feature-based bone segmentation in US (APS) [57] method using the FEM phantom,experimental phantom and in vivo data. The APS method selects parametersautomatically for the Log-Gabor filters based on properties estimated directlyfrom the specific US image that is being analyzed. We calculate mean absoluteerror (MAE) which is defined asMAE = 1NN∑k=1| A(k)−G(k) |, (2.10)where N is the approximated number of columns on which the bone bound-ary spans, A is the matrix containing the actual bone boundary points (groundtruth), and G is the matrix containing the estimated bone boundary points. Asground truth, we use the estimated bone boundaries from the ideal elastogramand CT image for the FEM and experimental phantom, respectively. Note thatto estimate MAE for the SEP method, we use the bone boundary points beforeapplying the regression.47  010203040Height (mm)    0 10 20 30 40Width (mm)010203040Height (mm)0 10 20 30 40Width (mm)  (a) (b) (c) (d)B-mode Ideal Strain Strain by AM Method Envelop PowerAPS SEP 0.3 0.4 0.5 0.6 0.700.51.01.52.02.5Weight (λ)MAE (mm)Weight Analysis (h)00.51.01.52.0  10 15 20 25 30 35 40 SNR (dB)MAE (mm) SEPAPSMAE Analysis(g)3.0Fused Map Maximum Intensity Points0 10 20 30 40Width (mm) 0 10 20 30 40Width (mm)(e) (f )(i) (j)Figure 2.2: Illustration of bone boundary detection using the FEM phantom. (a) B-mode image,(b) ideal strain image, (c) strain image generated by the AM method, (d) envelop power map,(e) fused map, (f) maximum intensity points on the fused map, (g) estimated bone boundaryby the APS method, (h) estimated bone boundary by the SEP method, (i) MAE analysis ofthe SEP and APS methods at different SNRs, and (j) weight analysis at 10dB SNR (weight atminimum error is shown with red arrow).2.7.1 FEM ResultsWe provide comparative qualitative and quantitative results of the SEP and APSmethods using the FEM data. In Figs. 2.2(a), (b), (c), (d) and (e), we showthe B-mode, ideal strain, estimated strain S, envelop power, and fused maps,respectively. In addition, in Fig. 2.2(f), we show the maximum intensity points on48the fused map. Note that we are dealing with anatomy that includes bone, hence,the RF intensities fall sharply (compared to soft tissue) beneath the bone surface.In such a situation, the regularization terms in the cost equation become largerthan the signal part in the region beneath the bone surface and over-smoothingtends to occur [82]. In Figs. 2.2(g) and (h), we show the bone boundariesdetected by the APS and SEP methods, respectively. The figures demonstratethat the APS method produces several false positive bone responses (indicatedby white arrows) due to the presence of ridge-like features in the B-mode image.In contrast, the bone boundary produced by the SEP method (Fig. 2.2(h)) isapparently free of such artifacts.We also compare the quantitative performance of the SEP and APS methodsin terms of MAE in Fig. 2.2(i) with four different signal-to-noise ratio (SNR)simulations (40, 30, 20, and 10dB) with 100 realizations each. We employed meanabsolute error calculations between the actual and estimated bone boundaries toshow the extent of the false positives created by the comparing methods. As canbe seen in this figure, up to 20dB SNR, the APS method produces greater meanMAE than that of the SEP method, though the mean MAE of both the methodsare almost flat and standard deviations are close to zero. In addition, the meanMAE of the APS method at 10dB SNR is higher though the standard deviationof the SEP method in this case somewhat higher. Therefore, we can say that theSEP method shows better performance in terms of reduced false positive boneresponses than that of the APS method.2.7.2 Experimental Phantom ResultsFigure 2.3 demonstrates the qualitative and quantitative performance comparisonof the SEP and APS methods using the experimental phantom. In Fig. 2.3(a), weshow the B-mode image. Figure 2.3(b) shows the CT projection slice along whichthe US scanning is performed and fiducials in Fig. 2.3(a) corresponds to fiducialsin Fig. 2.3(b) (shown with yellow arrows) in the same order. The estimatedstrain image is shown in Fig. 2.3(c). We show the detected bone boundaries bythe APS and SEP methods in Figs. 2.3(e) and (f), respectively. We can see fromFig. 2.3(e) that the bone boundary image estimated by the APS method produces49   (a) (b)MAE AnalysisMAE (mm)B-modeSEPAPSCT Image  Strain by AM Method(c)(e) (f )10 15 20 25 30 35 40SNR (dB)00.51.51.02.0  SEPAPSMAE (mm)(g)(d)Maximum Intensity points on FMFigure 2.3: Illustration of the bone boundary detection using the experimental phantom. (a)B-mode image, (b) a CT projection along which US scanning is performed, (c) strain image,(d) maximum intensity points overlaid on the fused map, (e) estimated bone boundary by theAPS method, (f) estimated bone boundary by the SEP method, and (g) MAE analysis of theSEP and APS methods at different SNRs.false positive responses in soft tissue regions and fails to estimate bone boundarywhere US back-scatter response is very weak (shown with white circles in Figs.2.3(a) and (e)). In contrast, the bone boundary estimated by the SEP methodhas no discontinuity as well as produces no noticeable false positive response.We also show the quantitative performance of both the methods in Fig. 2.3(g).50Here too, we use four different SNRs (40, 30, 20, and 10dB) with 100 realizationsfor each of them to analyze the MAE. Since the quality of the ultrasound machine-acquired RF data is dependent on user-defined parameters, the signals may havedifferent signal to noise ratios in different situations. Therefore, we have useda range of SNRs to evaluate the robustness of our proposed SEP scheme. Wecan see from this figure that the mean MAE values of the SEP method are lowerthan that of the APS method at all SNRs though the standard deviations of theSEP method are slightly higher. Therefore, from here too, we can say that theSEP method shows better performance in terms of reduced false positive boneresponses than that of the APS method.2.7.3 In Vivo Results2.7.3.1 Case 1Figure 2.4 demonstrates the qualitative performance comparison of the SEP andAPS methods using the volunteer-I data set. The scanned region, B-mode im-age, and strain image are shown in Figs. 2.4(a), (b), and (c), respectively. Themaximum intensity points overlaid on the fused map are shown in Fig. 2.4(d).The estimated bone boundaries by the APS and SEP methods are shown inFigs. 2.4(e) and (f), respectively. We see from Fig. 2.4(e) that the APS methodproduced false positive response for bone in the soft tissue interfaces (shown withwhite arrows). In contrast, the bone boundary estimated by the SEP method isfree from any false positive related artifact and better match the shapes visiblein the corresponding B-mode image (see Figs. 2.4(b) and (f)).2.7.3.2 Case 2We show the qualitative performance comparison of the SEP and APS methodsusing the volunteer-II data set in Fig. 2.5. The scanned region, B-mode im-age, and strain image are shown in Figs. 2.5(a), (b), and (c), respectively. Themaximum intensity points overlaid on the fused map are shown in Fig. 2.5(d).The estimated bone boundaries by the APS and SEP methods are shown inFigs. 2.5(e) and (f), respectively. Note that for volunteer-II, only half of the totalscan-lines are considered for both methods since almost half of the transducer51    (a)Transducer Position B-mode Strain ImageAPS SEP(a) (b) (c)(d) (e)Maximum Intensity Points on FM(f )Figure 2.4: Illustration of the bone boundary detection using the in vivo volunteer-I dataset. (a) shows the scanned region on the volunteer (shown with red rectangle), (b) is the B-mode image (bone surface is indicated with red arrows), (c) is the strain image, (d) shows themaximum intensity points overlaid on the fused map, (e) shows the detected bone boundaryby the APS method, and (f) shows the detected bone boundary by the SEP method.face laterally was in the air. As can be seen in Fig. 2.5(c) that the estimatedstrain image does not indicate the bone surface clearly, however, the bone in-tensity is quite higher than that of the soft tissue region (see Fig. 2.5(b)). So,the envelope power map is expected to complement with the bone indicating in-formation. We see from Fig. 2.5(e) that the APS method failed to produce acontinuous bone boundary. In addition, it produced false positive response forbone in the soft tissue interfaces (shown with white arrow). In contrast, the boneboundary estimated by the SEP method is free from any false positive related52     (a) (a) (b) (c)(d) (e) (f )Transducer Position B-mode Strain ImageAPS SEPMaximum Intensity Points on FMFigure 2.5: Illustration of the bone boundary detection using the in vivo volunteer-II dataset. (a) shows the scanned region on the volunteer (shown with red rectangle), (b) is the B-mode image (bone surface is indicated with red arrows), (c) is the strain image, (d) shows themaximum intensity points overlaid on the fused map, (e) shows the detected bone boundaryby the APS method, and (f) shows the detected bone boundary by the SEP method.artifact and better match the shapes visible in the corresponding B-mode image(see Figs. 2.5(b) and (f)).2.7.3.3 Case 3Figure 2.6 demonstrates the qualitative performance comparison of the SEP andAPS methods using the volunteer-III data set. The scanned region, B-modeimage, and strain image are shown in Figs. 2.6(a), (b), and (c), respectively. Themaximum intensity points overlaid on the fused map are shown in Fig. 2.6(d).53  (a) (a) (b)(d) (e)(c)(f )Transducer Position B-mode Strain ImageAPS SEPMaximum Intensity Points on FMFigure 2.6: Illustration of the bone boundary detection using the in vivo volunteer-III dataset. (a) shows the scanned region on the volunteer (shown with red rectangle), (b) is the B-mode image (bone surface is indicated with red arrows), (c) is the strain image, (d) shows themaximum intensity points overlaid on the fused map, (e) shows the detected bone boundaryby the APS method, and (f) shows the detected bone boundary by the SEP method.The estimated bone boundaries by the APS and SEP methods are shown inFigs. 2.6(e) and (f), respectively. We see from Fig. 2.6(e) that the APS methodfailed to produce a continuous bone boundary and also produced false positiveresponses at the soft tissue interfaces (shown with white arrow). In contrast, thebone boundary estimated by the SEP method better match the shapes visible inthe corresponding B-mode image (see Figs. 2.6(b) and (f)). In addition, it doesnot produce any false positive artifact.542.7.3.4 Case 4(a) (b) (c)(d) (e) (f )Transducer Position B-mode Strain ImageAPS SEPMaximum Intensity Points on FMFigure 2.7: Illustration of the bone boundary detection using the in vivo volunteer-IV dataset. (a) shows the scanned region on the volunteer (shown with red rectangle), (b) is the B-mode image (bone surface is indicated with red arrows), (c) is the strain image, (d) shows themaximum intensity points overlaid on the fused map, (e) shows the detected bone boundaryby the APS method, and (f) shows the detected bone boundary by the SEP method.We show the qualitative performance comparison of the SEP and APS meth-ods using the volunteer-IV data set in Fig. 2.7. The scanned region, B-modeimage, and strain image are shown in Figs. 2.7(a), (b), and (c), respectively. Themaximum intensity points overlaid on the fused map are shown in Fig. 2.7(d).The estimated bone boundaries by the APS and SEP methods are shown inFigs. 2.7(e) and (f), respectively. Here, the estimated strain image satisfactorilyindicates the bone surface (see Fig. 2.7(c)). In addition, the B-mode intensityfor the bone surface is also quite prominent than that for the soft tissue region.55Therefore, the envelope power map is expected to complement with the bone in-dicating information. We see from Fig. 2.7(e) that the estimated bone boundaryby the APS method is not continuous. In addition, it produced false positiveresponse for bone in the soft tissue interfaces (shown with white arrows). Incontrast, the bone boundary estimated by the SEP method is free from any falsepositive related artifact and better match the shapes visible in the correspondingB-mode image (see Figs. 2.7(b) and (f)).2.7.3.5 Case 5Finally, in Fig. 2.8, we demonstrate the qualitative performance comparisonof the SEP and APS methods using the volunteer-V data set. The scannedregion, B-mode image, and strain image are shown in Figs. 2.8(a), (b), and (c),respectively. The maximum intensity points overlaid on the fused map are shownin Fig. 2.8(d). The estimated bone boundaries by the APS and SEP methods areshown in Figs. 2.8(e) and (f), respectively. Here also, we see that the estimatedbone boundary by the APS method is not continuous (see Fig. 2.8(e)). Inaddition, it produced false positive response for bone in the soft tissue interfaces(shown with white arrows). In contrast, similar to previous cases, the boneboundary estimated by the SEP method better match the shapes visible in thecorresponding B-mode image (see Figs. 2.8(b) and (f)).2.8 Advantages and LimitationsWe discussed a novel method for robust bone boundary localization based onthe fusion of strain-imaging and envelope signal power detection. We combinedreal-time strain-imaging based on analytic minimization of regularized cost func-tions with an envelope power map of the pre-compression RF frame. Our resultsdemonstrated reduced bone localization error through better elastogram reso-lution, adoption of envelope power map of higher dynamic range, addressingthe false positive bone response, and exploiting the smoothing feature of the AMmethod to better delineate the bone boundary in strain image. We demonstratedour improved performance on a wide range of validation data including a simu-lated FEM phantom, a physical experimental phantom, and in vivo data of the56(a) (b) (c)(d) (e) (f )Transducer Position B-mode Strain ImageAPS SEPMaximum Intensity Points on FMFigure 2.8: Illustration of the bone boundary detection using the in vivo volunteer-V dataset. (a) shows the scanned region on the volunteer (shown with red rectangle), (b) is the B-mode image (bone surface is indicated with red arrows), (c) is the strain image, (d) shows themaximum intensity points overlaid on the fused map, (e) shows the detected bone boundaryby the APS method, and (f) shows the detected bone boundary by the SEP method.human with reported improvements of approximately 13% and 15% in terms ofMAE in the FEM and experimental phantom tests, respectively, when comparedwith current state-of-the-art.Although this study demonstrates great promise in 2D bone segmentation,however, the empirical selection of the value for λ may result in great challengesin dealing with the data from patients having clinical conditions. In addition,proper choice of the linear fitting window length is also critical. Therefore, a data-driven analysis and choice of these parameter values are of great importance. This57method also lacks a feature to identify any discontinuity in the bone boundary.58Chapter 3Improved 2D Bone ContourExtraction Using Data-drivenParametersIn this chapter, we significantly improve our preliminary work described in theprevious chapter and also mitigate the limitations associated with it by intro-ducing automatic parameter selection processes that replace all of the intuitiveassumptions and thus, make the method more robust and self-tuning. We improveour previous 2D strain-imaging-based bone segmentation approach discussed inChapter 2 by incorporating a depth-dependent cumulative power of the enve-lope into the elastographic data as well as incorporating an echo de-correlationmeasure-based weight to fuse the strain and envelope map. We also use a datadriven scheme to detect the presence of any bone discontinuity in the scanned USimage and introduce a multivariate non-parametric Gaussian mixture regressionthat to be used over the maximum intensity points of the fused map.593.1 Improved Modified Envelope Map (MEM)and Modified Strain Map (MSM) Estima-tionWe already discussed in Section 2.2 that we use a real-time strain-imaging [82]based on analytic minimization (AM) of regularized cost functions. However, asdiscussed in [82], an inherent limitation of that method is that over-smoothing ofthe strain image occurs in places where the RF echo is weaker, which results inincrease in the share of the regularization term in the cost Cj defined asCj(∆a1, ..,∆am,∆l1, ..,∆lm) =m∑i=1{[I1(i, j)− I2(i+ ai + ∆ai, j + li +∆li)]2+ α(ai + ∆ai − ai−1 +∆ai−1)2 + βa(li + ∆li − li−1 +∆li−1)2+ β′l(li +∆li − li,j−1)2}, (3.1)where α, βa and β′l are the tuneable parameters, ai and ∆ai are the integerand sub-sample axial displacements, respectively, and li and ∆li are the integerand sub-sample lateral displacements, respectively, of the ith sample in the jthscan-line.In our subsequent work, we exploit this limitation to improve the bone de-lineation accuracy further. Since the RF echo and its envelope is already weakerbeneath the bone surface due to very high US beam reflection at the bone tissueinterface, we further modify the envelope map by reducing the intensity associ-ated with the soft tissue region (area between the bone and the transducer face)using the axially cumulative power and produce modified envelope maps, MEM1and MEM2 from I1 and I2, respectively, asMEMr(i, j) =Er(i, j)Pr(i,j)max[Er(i, j)Pr(i,j)], (3.2)where r corresponds to 1 and 2 denoting the pre- and post-compression images,respectively, Er is the envelope estimated using the Hilbert transform of Ir [45],and Pr(i, j) =∑ipo=1 E2r (po, j). Finally, we estimate strain (S) using MEM1 andMEM2 which results an over-smoothing in both the top and bottom regions of60the bone boundary. For robust bone boundary estimation, we combine S withMEM1. Therefore, in order to get an equivalent dynamic range of S with that ofthe MEM1, we estimate a modified strain map MSM = Sˆ/max[Sˆ], whereSˆ(i, j) ={−S(i, j); if S(i, j) < 0S(i, j); otherwise(3.3)3.2 Data-driven Weight EstimationHere, we recap again from Section 2.4 on the estimation of the fused map usingthe MSM and MEM1 asFM = λ×MSM + (1− λ)×MEM1, (3.4)where, λ controls the contribution of MSM into the fused map. In this section,we target to automate the process of choosing the value of λ.To observe the effect of different values of λ on the bone boundary estimationaccuracy, we use FEM-II simulation phantom (see Section 3.5 for detail). We usethe actual and estimated bone boundaries of FEM-II to estimate mean absoluteerror (MAE) at different values of λ which is shown in Fig. 3.1(a). For this MAEanalysis, we use pre- and post-compression RF frames having 10dB signal-to-noiseratio (SNR), where the SNR is deliberately chosen to be lower than the typicallyencountered values (approximately 18∼22dB) in RF data after beamforming [86]to simulate a worse case scenario. We also use 1% strain data to ensure that theecho de-correlation is low enough. We see from Fig. 3.1(a) that the lowest MAEis at λ = 0.5. However, by the increase of echo de-correlation, estimated strainimage becomes more noisy. Free hand elastography is prone to echo de-correlationand by the increase of applied pressure, echo de-correlation increases significantlythat leads to noisy strain map [83]. To demonstrate this effect, we use FEM-I simulation phantom (see Section 3.5 for detail) and estimate the normalizedcross-correlation (NCC) between pre- and post-compression RF frame of FEM-Iat different applied strains (i.e., 1%-8% applied strains). The mean NCC peaksat different applied strains are plotted in Fig. 3.1(b). We can see from this figurethat the peaks of the NCC coefficients deteriorate and thus, the quality of strainimages deteriorate with the increase of applied strain.61We incorporate the echo de-correlation measure between the pre- and post-compression RF frames to automatically choose a suitable value for λ. To esti-mate the degree of echo de-correlation, we consider a smaller (i.e., 10mm×10mm)region-of-interest (ROI) closer to the transducer face in the pre- and post-compressionecho frames. We choose axial length of the ROI to be 10mm so that three 1DRF signal windows of length L = 5mm each with 50% axial overlap can be usedto estimate cross-correlation coefficients. We also choose L to be 5mm since thecross-correlation between the pre- and post-compression RF windows of length3∼5mm (with transducer center frequency 5∼10MHz) produces satisfactory lagto estimate tissue displacement [80]. We use lateral length of the ROI to be10mm, however can be chosen up to the width of the RF frame. After definingthe ROI, we stretch the post-compression signal windows in the ROI by a factorα, i.e., Iα(i) = I2(αi). By neglecting the noise term (i.e., να(i) = 0 in Eqn. (2.2)),we define the NCC ρα(k) (≤ 1) as [83]ρα(k)=∑Li=1 I1(i) · Iα(i+ k)√∑Li=1{I1(i)}2∑Li=1{Iα(i)}2. (3.5)Equation (3.5) becomes maximum for α =  with the assumption of p(αt) ∼= p(t).After estimating ρα(k) for all the scan-lines inside the ROI, the mean NCC peakis estimated asρavg =1MM∑w=1max{ρα(k)}w, (3.6)where M is the total number of 1D signal windows inside the ROI. Note that tofind individual peak among NCC coefficients of each window in the sub-samplelevel, we use cosine interpolation [83].Now, based on the value of ρavg, we consider three cases: (i) if ρavg is equal orgreater than 0.9, we can consider echo de-correlation between the pre- and post-compression echo frames is significantly low [89]. In that case, λ = 0.5 works thebest as evident from Fig. 3.1(a). (ii) If ρavg is less than 0.5, it is obvious thatthe degree of echo de-correlation is greater than the correlation. Then we chooseλ = 0 (i.e., no contribution from the noisy strain image). And (iii) if the valueof ρavg varies between 0.5 and 0.9, we assume a linear relation of λ with ρavg and621 2 3 4 5 6 7 800.10.20.30.40.50.60.70.80.91Strain (%)NCC PeaksNCC Peak Analysis(a) (b)Weight (λ)0.3 0.4 0.5 0.6 0.700.51.01.52.02.5Weight (λ)MAE (mm)Weight Analysis 3.0Average NCC Coe!cient (        )NCC Dependent Weightρavg(c)0.5 0.6 0.7 0.8 0.90.050.100.150.200.300.350.400.450.500.00.25Figure 3.1: (a) Mean absolute error (MAE) analysis for different values of λ at 10dB SNR(weight at minimum error is shown with red arrow), (b) Normalized cross-correlation coeffi-cients with respect to applied strains for FEM-I phantom (see Section 3.5 for detail), and (c)Normalized cross-correlation coefficient dependent linear weight λ.thus choose the value for λ from a linear function of ρavg (see Fig. 3.1(c)) definedasλ = 54(ρavg − 0.5). (3.7)It might seem that the data-driven weight has a number of fix parameters, e.g.,the upper and lower bound of the weight and ρavg. However, these parametersare chosen based on the state-of-the-art practice (i.e., the method in [89] runsin the commercial US scanners manufactured by the Ultrasonix medical corpo-ration, Richmond, BC) and vigorous observations made on the bone localizationaccuracy using different types of data. In addition, our sole target is to controlthe contribution of a strain image into the fused map based on the fact thathow good that strain image is. Therefore, we believe that our data-driven weightserves this purpose irrespective of data types, and its efficacy is quantitativelydemonstrated in the results section (see Section 3.6.3).However, after estimating the fused map, the location of maximum intensitypoint along each scan-line of the map is then used as the initial bone boundaryY(= [y1, y2, ..., yn]), where yj is the axial sample number at jth scan-line in theFM. We also use Grubbs test [90] to discard any outlier in Y.3.3 Bone Discontinuity DetectionTo detect possible discontinuity in the bone surface, we estimate the lease-square-error-based gradient ∇Y = dY˜/dX [91], where X(= [x1, x2, ..., xn]), Y˜ is the63bone boundary after regression, and xj represents the location of yj in the lateraldirection. Although the selection of proper window length for regression is crucial,our primary target here is to find a comparably larger jump in the values of ∇Yin the location of bone discontinuity. Therefore, we use larger regression windowof length 5mm so that the noisy estimates can be better minimized as well as anoticeable jump in the location of bone discontinuity is preserved in the valuesof ∇Y. Then we use the local statistics of ∇Y such that any point on ∇Yexceeds µ + 3σ (where µ and σ are the mean and standard deviation of ∇Y) isflagged as a bone discontinuity. We choose µ+ 3σ to cover a confidence intervalof 0.997 assuming that ∇Y is normally distributed. Each continuous portion ofY is considered as a segment and after finding all the segments, we use Gaussianmixture regression over each of the segments separately.3.4 Gaussian Mixture Regression (GMR)In this work, we use GMR that overcomes the following drawbacks of linearregression which we used in our previous work [29]:1. Choosing a suitable length for regression window is ambiguous since smallerwindows cannot produce smoother bone boundary if Y is noisy. On theother hand, comparatively larger regression window causes loss of localcurvature information in Y.2. Parametric linear regression models are often too rigid to model generalnonlinear patterns in data space and thus, a more flexible model of non-parametric regression is required [92].Now let fX,Y (x, y) =∑Kk=1 pikφ(x, y;µk,Σk) be the joint density function ofthe maximum intensity points identified in the FM. K is the total number ofclusters present in Y as estimated from a spectral clustering analysis [93]. Thespectral clustering approach is based on an eigenvalue analysis of the affinitymatrix of the interrogated point clouds. Eigenvalues 1 or close to 1 typically in-dicate a natural estimate of the number of clusters to be used to model the data.However, if the groups are not clearly separated, once noise is introduced, the641 2 3 4 5 6 7 8 90.90.910.920.930.940.950.960.970.980.99Cluster Number (K  )Cost ( J )iFigure 3.2: The alignment cost for varying cluster numbers.values start to deviate from 1, thus the criterion of choice becomes tricky. Thus,authors in [93] adopted an approach which relies on the structure of the eigenvec-tors. This idea is based on the fact that each eigenvector column corresponding toan eigenvalue 1 will have only a single non-zero entry. In addition, the eigenvectorcolumns corresponding to eigenvalues 1 form a canonical coordinate system in anideal case. Thus, for a noisy data set, the algorithm in [93] forms a cost func-tion with respect to different possible cluster numbers Ki, where for each Ki, theeigenvector columns are repeatedly rearranged such that the rearranged first Kicolumns form a canonical coordinate system and contain only a single non-zeroentry in each of first Ki columns. Let Z be the matrix of size n ×Ki obtainedafter rotating (rearranging) the eigenvector matrix, and Mp = maxq Zpq. Thenthe cost function is defined as [93]J =n∑p=1Ki∑q=1Z2pqM2p. (3.8)The cluster number Ki for which the cost J (Eqn. (3.8)) becomes minimum istaken as the final cluster number K. For example, we plot the cost values fordifferent cluster number Ki using the in vivo volunteer-I data set (see Section3.6.3) in Fig. 3.2. Here we see that the cost is minimum for Ki = 3.However, pik, µk = [µkX ;µkY ] and Σk = [ΣkX ΣkXY ; ΣkY X ΣkY Y ] in the jointdensity function are the probability, mean and covariance matrix, respectively, of65the data points in the kth cluster and are estimated from an expectation maxi-mization algorithm [94], and φ(x, y;µk,Σk) is the probability density function ofpoints inside kth cluster defined asφ(x, y;µk,Σk)= 12pi√|Σk|exp[−12(D − µk)TΣ−1k (D − µk)], (3.9)where D = [X;Y]. Then, from fX,Y (x, y), our GMR equation can be derivedwhich takes the form as [92]b(x) = E[Y|X = x] =K∑k=1wk(x)dk(x), (3.10)where dk(x) = µkY + ΣkY XΣ−1kX(x − µkX), and µkX and ΣkX are the mean andcovariance matrix of the marginal density function (fX) of X, in the kth datacluster. In addition, wj(x) is the mixing weight defined aswk(x) =pikφ(x;µkX ,ΣkX)∑Kk=1 pikφ(x;µkX ,ΣkX). (3.11)GMR covers a spectrum of regression models of varying flexibility, ranging fromthe multivariate nonparametric kernel regression (K = N ;N>1) to the classicallinear regression model (K = 1) [92].3.5 Validation Setup3.5.1 Simulated PhantomWe built three 40mm×40mm (axial×lateral) FEM phantoms using the analysissoftware ANSYS (ANSYS, Inc., Canonsburg, PA) and the ultrasound simula-tion software Field II [87]. The first phantom was generated to show the echode-correlation effect with the applied mechanical stress where total number ofnodes were 54,405. Our second and third phantoms mimicked 2D US scans ofa normal human distal radius bone (FEM-II) and a fractured human distal ra-dius bone (FEM-III) with a total number of 55,180 nodes. The stiffness of the66  10kPa10GPa010203040 0 10 20 30 40 0 10 20 30 400 10 20 30 40Height (mm)Width (mm) Width (mm)Width (mm)(b) (c) (d)10kPa10GPa0 10 20 30 40Width (mm) (e)00.0020.0040.0060.0080.010.0120.0140.0160.0180.020 10 20 30 40Width (mm)(a)FEM - I FEM - II FEM - IIIFigure 3.3: FEM simulation phantoms. (a) FEM-I: ideal elastogram of a homogeneous softtissue of stiffness 10kPa. (b) FEM-II: normal distal radius bone and soft tissue region ofstiffness 10GPa and 10kPa, respectively, and (c) corresponding ideal elastogram. (d) FEM-III:fractured distal radius bone and soft tissue region of stiffness 10GPa and 10kPa, respectively,and (e) corresponding ideal elastogram.homogeneous soft tissue and bone regions were set to 10kPa and 10GPa, respec-tively, (see Fig. 3.3) as previously reported [88]. Our phantom was compressedin the lateral direction from the top using a planar compressor that was widerthan the phantom. We set an applied pressure level that corresponds to 1%average strain. An ultrasonic transducer of center frequency, f0 = 5MHz andband-width = 50% was used to simulate the phantom scan from the top. Thetotal number of scan lines was set to 128. The resulting pixel size of the 2D datawas 0.015mm×0.3mm. Figures 3.3(a), (c) and (e) show the ideal elastograms ofFEM-I, FEM-II and FEM-III, respectively, at 1% applied strain.3.5.2 Experimental PhantomWe constructed an experimental phantom using a radio-opaque Sawbones hemi-pelvis (Sawbones, Pacific Research Laboratories, Inc., Vashon Island, WA), modelnumber 1297-22. A portion of the pelvis was suspended in polyvinyl chloride(PVC) gel (see Fig. 2.1(a)) and placed in an acrylic tube (see Fig. 2.1(b)). Ahigh resolution peripheral quantitative CT machine, model HR pQCT XtremeCT (Scanco USA, Inc., Wayne, PA) was used to acquire a single 482×482×402(lateral×axial×elevation) volume with a resolution of 0.25mm×0.25mm×0.25mm.The US images were acquired using a SonixRP (Ultrasonix Medical Corporation,Richmond, BC) scanner in the Center for Hip Health and Mobility, Vancouver67Coastal Health Authority, Vancouver, BC, Canada. We used a L14-5W/60 lineararray probe operating at 10MHz for collecting data for the 2D implementations.3.5.3 In Vivo DataWe acquired five sets of data from five volunteers (volunteer-I: 25-year old male;volunteer-II: 33-year old male; volunteer-III: 26-year old male; volunteer-IV: 24-year old male; volunteer-V: 27-year old male) after obtaining prior consent. Alldata were acquired with free-hand compression. The US images were acquiredusing a SonixRP (Ultrasonix Medical Corporation, Richmond, BC) scanner inthe Center for Hip Health and Mobility, Vancouver Coastal Health Authority,Vancouver, BC, Canada. Here also, we used a L14-5W/60 linear array probeoperating at 10MHz. The study was approved by UBC clinical research ethicsboard.3.6 ResultsWe provide comparative results of our proposed 2D strain and envelope power-based bone detection method (2DSE) with the previously reported adaptivelyparameterized local phase feature-based (2DOPS) [57] using the FEM phantoms,experimental phantom and in vivo data. The 2DOPS method selects parametersautomatically for the Log-Gabor filters based on properties estimated directlyfrom the specific US image that is being analyzed. We calculate mean absoluteerror (MAE) which is defined asMAE = 1PP∑p=1| A(p)−G(p) |, (3.12)where P represent the number of data points spanned by the bone boundary inthe lateral directions, A is a matrix containing the ‘ground truth’ bone boundarypoints, and G is a matrix containing the estimated bone boundary points. Asground truth, we used the bone geometry defined for Field II simulation for theFEM, the estimated bone boundary from the computed tomography (CT) imagefor the experimental phantom, and expert delineation of bone boundaries on68B-mode images for the in vivo data. We use fiducial-based CT and US bonealignment before estimating the MAE for the experimental phantom. Note thatwe use bottom-up ray casting in the phase-symmetry images for the 2DOPSmethod to get G as suggested in [57]. Also note that we use the initial boneboundary estimates Y to estimate the MAE for the 2DSE method for a farecomparison. Also note that since we use expert delineated bone contours asground truth for the in vivo data, therefore, the estimated error using Eqn. 3.12for the in vivo data is called mean absolute fitting error.        0.10.20.30.40.50.60.70.80.910.10.20.30.40.50.60.70.80.9100102030B-Mode Strain Image MEM GM Over Fused Map2DOPSHeight (mm)Height (mm)Width (mm) Width (mm)(a) (b) (c) (d)(e) (f)0.520 3039.5 10 Width (mm)3 37 Width (mm) Width (mm) Width (mm)20 30103 37 20 30103 37 20 30103 3720 30103 37 20 30103 371020300.539.512DSEFigure 3.4: Illustration of 2D bone boundary detection using the FEM-I phantom. (a) B-modeimage, (b) strain image generated by the proposed approach, (c) MEM1, (d) Gaussian mixturesrepresentation over fused map, and estimated bone boundary by the (e) 2DSE and (f) 2DOPSmethods (arrows show leakage). Bone boundary is shown after bottom up ray-casting for the2DOPS method.3.6.1 FEM ResultsWe provide comparative qualitative results of the 2DSE and 2DOPS methodsusing the FEM data. In Figs. 3.4(a), (b), (c) and (d), we show the B-modeimage, estimated strain image, MEM1, and fused map, respectively, for FEM-69        0.10.20.30.40.50.60.70.80.910.10.20.30.40.50.60.70.80.9100B-Mode Strain Image2DOPS2DSEGM Over Fused Map102030Height (mm)0.539.5102030Height (mm)0.539.5Width (mm)20 30103 37 Width (mm)20 30103 37 Width (mm)20 30103 37 Width (mm)20 30103 37Width (mm)20 30103 37 Width (mm)20 30103 37MEM1(a) (b) (c) (d)(e) (f)Figure 3.5: Illustration of 2D bone boundary detection using the FEM-II phantom. (a) B-modeimage, (b) strain image generated by the proposed approach, (c) MEM1, (d) Gaussian mixturesrepresentation over fused map, and estimated bone boundary by the (e) 2DSE and (f) 2DOPSmethods (arrow shows leakage). Bone boundary is shown after bottom up ray-casting for the2DOPS method.II. To show the efficacy of our proposed method in case of discontinuous boneboundary localization, we show results for FEM-III in Fig. 3.5. In Figs. 3.5(a),(b), (c) and (d), we show the B-mode, estimated strain S, MEM1, and fusedmaps, respectively. In Figs. 3.4(e) and 3.5(e), we show the bone boundariesdetected by the 2DSE method for the FEM-II and FEM-III, respectively. Inaddition, we show the bone boundaries detected by the 2DOPS method for theFEM-II and FEM-III in Figs. 3.4(f) and 3.5(f), respectively. Since the 2DOPSmethod depends on ray-casting to extract final bone boundary from the phasesymmetry images, these methods become prone to “leaking” through weak areasin the phase symmetry profiles. As evident we see in Figs. 3.4(f) and 3.5(f)that the ray-casting leaks through weak spots in the phase symmetry profiles ofthe bone and catches intensities associated with the soft tissue layer (indicatedby yellow arrows in Figs. 3.4(f) and 3.5(f)). In contrast, the bone boundaries70  Strain Image 2DOPS00.10.20.30.40.50.60.70.80.91B-mode00.10.20.30.40.50.60.70.80.911.510203038.5Height (mm)Width (mm) Width (mm)1.510203038.52DSEGM Over Fused Map10 20 30 40 50 57.52.5 Width (mm)Height (mm)10 20 30 40 50 57.52.5 10 20 30 40 50 57.52.510 20 30 40 50 57.52.5 Width (mm) 10 20 30 40 50 57.52.5 Width (mm) 10 20 30 40 50 57.52.5 Width (mm)MEM1(a) (b) (c)(d) (e) (f)Figure 3.6: Illustration of 2D bone boundary detection using the experimental phantom. (a)B-mode image, (b) strain image generated by the proposed approach, (c) MEM1, (d) Gaussianmixtures representation over fused map, (e) estimated bone boundary by the 2DSE method,and (f) estimated bone boundary by the 2DOPS method (bone boundary is shown after bottomup ray-casting).produced by the 2DSE method (see Figs. 3.4(e) and 3.5(e)) appear to be free ofsuch artefacts.3.6.2 Experimental Phantom ResultsFigure 3.6 demonstrates the qualitative performance comparison of the 2DSEand 2DOPS methods using the experimental phantom. In Figs. 3.6(a), (b),(c) and (d), we show the B-mode image, strain image, MEM1, and fused map,respectively. We show the detected bone boundaries by the 2DSE and 2DOPSmethods in Figs. 3.6(e) and (f), respectively. We can see from Fig. 3.6(f) that thebone boundary estimated by the 2DOPS method is discontinuous. In contrast,the bone boundary estimated by the 2DSE method is continuous and smooth(see Fig. 3.6(f)).7100.10.20.30.40.50.60.70.80.91Volunteer - I Volunteer - II Volunteer - III Volunteer - IVB-mode00.10.20.30.40.50.60.70.80.9100.10.20.30.40.50.60.70.80.91       Strain (traditional)     Strain (proposed)(a) (b) (c) (d)(e) (f) (g) (h)(i) (j) (k) (l)Figure 3.7: Illustration of the proposed strain-imaging performance. (a)-(d) show the B-modeimages where overlaid red-dotted boxes indicate the bone locations. (e)-(h) show the strainimages generated by the traditional method adopted in [29] (red arrows indicate the true bonelocations while yellow arrows indicate artifacts associated with soft tissue layers), and (i)-(l)show the strain images generated by the proposed method (red arrows indicate the true bonelocations).3.6.3 In Vivo ResultsIn Fig. 3.7, we show the strain estimation performance of the proposed methodand that used in [29]. We see from Figs. 3.7(e)-(h) that the method used in [29]produces confounding artifacts in the soft tissue layers (indicated with yellowarrows). In contrast, the strain images produced by our proposed method isfree of such artifacts as evident in Figs. 3.7(i)-(l). Next we compare the boneboundary detection performance of the 2DSE and 2DOPS methods for the five invivo volunteer data cases. An Orthopedic expert delineated the bone boundarieson the B-mode images which we consider as the ground truths for a comparativeanalysis of the 2DSE and 2DOPS methods.723.6.3.1 Case 1Scanned Bone Area Strain ImageGM over Fused Map Expert, 2DSE, 2DOPSExpert Delineation2DSE2DOPSMEM1(a) (b) (c)(d) (e)Figure 3.8: Illustration of 2D bone boundary detection using the in vivo volunteer-I data set.(a) shows bone region scanned on the volunteer (with red rectangle), (b) is the strain imagegenerated by the proposed approach, (c) is the MEM1, and (d) is the fused map with theGaussian mixtures overlaid. Finally, in (e), orthopaedic expert delineated bone boundary, anddetected bone boundary by the 2DOPS and 2DSE methods are overlaid on the B-mode image(bone boundary produced by the 2DOPS method is shown after bottom up ray-casting).We show the qualitative performance comparison of the 2DSE and 2DOPSmethods using the volunteer-I data set in Fig. 3.8. The scanned region, strainimage, envelope power map, and the fused map are shown in Figs. 3.8(a), (b),(c), and (d), respectively. The expert delineated bone contour, and the estimatedbone boundaries by the 2DSE and 2DOPS methods are overlaid on the B-modeimage as shown in Fig. 3.8(e). We see from Fig. 3.8(e) that the bone boundary73produced by the 2DOPS method satisfactorily close to the expert delineatedcontour though there are some sporadic points (shown with red arrows) resultedfrom the false positive responses. However, the bone boundary estimated by the2DSE method better matches the shape as marked by the expert in the B-modeimages (see Figs. 3.8(e)). In addition, the 2DSE method does not produce falsepositives at soft tissue interfaces.3.6.3.2 Case 2Scanned Bone Area Strain ImageGM over Fused Map Expert, 2DSE, 2DOPSMEM1(a) (b) (c)(d) (e)Expert Delineation2DSE2DOPSFigure 3.9: Illustration of 2D bone boundary detection using the in vivo volunteer-II data set.(a) shows bone region scanned on the volunteer (with red rectangle), (b) is the strain imagegenerated by the proposed approach, (c) is the MEM1, and (d) is the fused map with theGaussian mixtures overlaid. Finally, in (e), orthopaedic expert delineated bone boundary, anddetected bone boundary by the 2DOPS and 2DSE methods are overlaid on the B-mode image(bone boundary produced by the 2DOPS method are shown after bottom up ray-casting).74Figure 3.9 demonstrates the qualitative bone localization performance com-parison of the 2DSE and 2DOPS methods using the volunteer-II data set. Thescanned region, strain image, envelope power map, and the fused map are shownin Figs. 3.9(a), (b), (c), and (d), respectively. The expert delineated bone con-tour, and the estimated bone boundaries by the 2DSE and 2DOPS methods areoverlaid on the B-mode image as shown in Fig. 3.9(e). We see from Fig. 3.9(e)that the bone boundary produced by the 2DOPS method varies noticeably insome places. This variation results from the false positive bone responses. Incontrast, the bone boundary estimated by the 2DSE method does not produceany false positive-based artifact and better matches the shape as marked by theexpert in the B-mode images (see Fig. 3.9(e)).3.6.3.3 Case 3We show the qualitative performance comparison of the 2DSE and 2DOPS meth-ods using the volunteer-III data set in Fig. 3.10. The scanned region, strainimage, envelope power map, and the fused map are shown in Figs. 3.10(a), (b),(c), and (d), respectively. The expert delineated bone contour, and the estimatedbone boundaries by the 2DSE and 2DOPS methods are overlaid on the B-modeimage as shown in Fig. 3.10(e). We see from Fig. 3.10(e) that the bone boundaryproduced by the 2DOPS method varies sporadically in some places due to thefalse positive responses. In this case, our 2DSE method also performs poorly asit could not align the produced bone contour with the expert delineated contouron the right side of the image. However, in other places, the 2DSE method per-forms better and the produced bone contour matches the shape as marked by theexpert in the B-mode images (see Figs. 3.10(e)).3.6.3.4 Case 4Figure 3.11 demonstrates the qualitative bone localization performance compar-ison of the 2DSE and 2DOPS methods using the volunteer-IV data set. Thescanned region, strain image, envelope power map, and the fused map are shownin Figs. 3.11(a), (b), (c), and (d), respectively. The expert delineated bone con-tour, and the estimated bone boundaries by the 2DSE and 2DOPS methods are75Scanned Bone Area Strain ImageGM over Fused Map Expert, 2DSE, 2DOPSMEM1(a) (b) (c)(d) (e)Expert Delineation2DSE2DOPSFigure 3.10: Illustration of 2D bone boundary detection using the in vivo volunteer-III dataset. (a) shows bone region scanned on the volunteer (with red rectangle), (b) is the strainimage generated by the proposed approach, (c) is the MEM1, and (d) is the fused map with theGaussian mixtures overlaid. Finally, in (e), orthopaedic expert delineated bone boundary, anddetected bone boundary by the 2DOPS and 2DSE methods are overlaid on the B-mode image(bone boundary produced by the 2DOPS method are shown after bottom up ray-casting).overlaid on the B-mode image as shown in Fig. 3.11(e). We see from Fig. 3.11(e)that the bone boundary produced by the 2DOPS method varies with respect tothe expert delineated contour. In contrast, the bone boundary estimated by the2DSE method better matches the expert delineated contour (see Figs. 3.11(e)).In addition, the 2DSE method does not produce false positives at soft tissueinterfaces.76Scanned Bone Area Strain ImageGM over Fused Map Expert, 2DSE, 2DOPSMEM1(a) (b) (c)(d) (e)Expert Delineation2DSE2DOPSFigure 3.11: Illustration of 2D bone boundary detection using the in vivo volunteer-IV dataset. (a) shows bone region scanned on the volunteer (with red rectangle), (b) is the strainimage generated by the proposed approach, (c) is the MEM1, and (d) is the fused map with theGaussian mixtures overlaid. Finally, in (e), orthopaedic expert delineated bone boundary, anddetected bone boundary by the 2DOPS and 2DSE methods are overlaid on the B-mode image(bone boundary produced by the 2DOPS method are shown after bottom up ray-casting).3.6.3.5 Case 5Finally, we show the qualitative performance comparison of the 2DSE and 2DOPSmethods using the volunteer-V data set in Fig. 3.12. The scanned region, strainimage, envelope power map, and the fused map are shown in Figs. 3.12(a), (b),(c), and (d), respectively. The expert delineated bone contour, and the estimatedbone boundaries by the 2DSE and 2DOPS methods are overlaid on the B-modeimage as shown in Fig. 3.12(e). We see from Fig. 3.12(e) that the bone boundary77Scanned Bone Area Strain ImageGM over Fused Map Expert, 2DSE, 2DOPSMEM1(a) (b) (c)(d) (e)Expert Delineation2DSE2DOPSFigure 3.12: Illustration of 2D bone boundary detection using the in vivo volunteer-V dataset. (a) shows bone region scanned on the volunteer (with red rectangle), (b) is the strainimage generated by the proposed approach, (c) is the MEM1, and (d) is the fused map with theGaussian mixtures overlaid. Finally, in (e), orthopaedic expert delineated bone boundary, anddetected bone boundary by the 2DOPS and 2DSE methods are overlaid on the B-mode image(bone boundary produced by the 2DOPS method are shown after bottom up ray-casting).produced by the 2DOPS method varies sporadically in some places due to thefalse positive responses. In contrary, the 2DSE method performs better and theproduced bone contour matches very closely with the shape as marked by theexpert in the B-mode images (see Figs. 3.12(e)).783.6.4 Quantitative ComparisonAt first we show the individual effect of using the depth-dependent cumulativepower, data-driven weight and Gaussian mixture regression in terms of MAEusing the in vivo data in Table 3.1. To show the effect of using the depth-dependent cumulative power, we choose data sets having satisfactory correlation(i.e., ρavg ≈ 0.90), and used fixed weight (i.e., λ = 0.5) and simple linear regressionfor both the 2DSE and SEP methods (SEP method is described in Chapter 2).While we show the effect of using the data-driven weight, we do not use thedepth-dependent cumulative power and use simple linear regression for both the2DSE and SEP methods. In addition, we use fixed weight (i.e., λ = 0.5) for theSEP method. Finally, to show the effect of using the GMR, we use GMR andlinear regression separately for the 2DSE method. From Table 3.1, we see thatthe use of the depth-dependent cumulative power improves the MAE performanceby approximately 13% on an average, the use of the data-driven weight improvesthe MAE performance by approximately 20% on an average, and the use of theGMR improves the MAE performance by approximately 16% on an average.We also compare the quantitative performance of the 2DSE and 2DOPS meth-ods using the FEM-II, III, and experimental phantoms in terms of MAE in Fig.3.13 with four different SNR simulations (40, 30, 20, and 10dB) with 100 real-izations each. In real clinical settings, SNR of the RF data after beamformingis reported to be approximately 18∼22dB [86]. Note that we manually excludethe leakage points in the bone boundaries produced by the 2DOPS method forestimating MAE. As can be seen in Figs. 3.13(a) and (b) for the FEM-II and III,respectively, that the 2DOPS method produces greater mean MAE than that ofthe 2DSE method although both the methods show a trend of increasing meanMAE with the decrease of SNR. We also see from Fig. 3.13(c) for the experi-mental phantom that the mean MAE values for the 2DSE method are lower thanthat of the 2DOPS method at all SNR values.We also compare the quantitative performance of the 2DSE and 2DOPS meth-ods using the in vivo data in terms of mean absolute fitting error in Table 3.2.Here also, we see that our 2DSE method improves the mean absolute fitting er-ror performance by approximately 55% on an average than that of the 2DOPSmethod.79Table 3.1: Mean absolute error (MAE) Analysis of the individual effects of using the depth-dependent cumulative power (at ρavg ≈ 0.90), data-driven weight (at ρavg ≈ 0.75 and 0.60),and GMR on the in vivo data.Volunteer MAE at MAE at MAE at MAE forρavg ≈ 0.90 ρavg ≈ 0.75 ρavg ≈ 0.60 Regression(mm) (mm) (mm) (mm)2DSE SEP 2DSE SEP 2DSE SEP GMR LinearI 0.36 0.43 0.39 0.48 0.44 0.59 0.32 0.36II 0.43 0.50 0.48 0.59 0.53 0.72 0.35 0.43III 0.58 0.66 0.64 0.71 0.65 0.82 0.50 0.58IV 0.39 0.45 0.45 0.53 0.50 0.65 0.29 0.39V 0.36 0.40 0.40 0.48 0.43 0.58 0.32 0.36Table 3.2: Mean absolute fitting error (mm) comparisons of the 2DSE and 2DOPS methodsusing the in vivo data.Volunteer 2DSE 2DOPS(mm) (mm)I 0.32 0.65II 0.35 0.92III 0.50 0.85IV 0.29 0.80V 0.32 0.718000.51.01.52.0  10 15 20 25 30 35 40 SNR (dB)MAE (mm)(a)  2DSE2DOPS00.51.01.52.010 15 20 25 30 35 40 SNR (dB)2DSE2DOPS(b)FEM-II FEM-III10 15 20 25 30 35 40SNR (dB)00.51.51.02.0  2DSE2DOPS(c)Experimental PhantomFigure 3.13: MAE (mm) comparisons of the 2DSE and 2DOPS methods at different SNR valuesusing (a) the FEM-II, (b) FEM-III, and (c) experimental phantoms.Table 3.3: F-Test for the variances in Fig. 3.14.2DOPS 2DSEMean 5.118 0.896Variance 8.36967 0.72568Observations 5 5Degrees of Freedom 4 4F 11.5336P (F ≤ f) one-tail 0.01808F Critical one-tail 6.38823In addition to the mean absolute fitting error, we also show in Fig 3.14 thenormalized distributions of the estimated bone boundary points by the 2DOPSand 2DSE methods with respect to the expert delineated bone contours usingthe in vivo volunteer data sets. Although the mean absolute fitting errors forboth the methods are below 1mm (Table 3.2), the spread of the signed distancesproduced by the 2DSE method appears to be lower than those for the 2DOPSmethod as evident from the standard deviation values of the distributions. Tofurther support this claim, we perform an F-test on the standard deviations ofthe distributions produced by the 2DOPS and 2DSE methods as shown in Table3.3. We see in this table that F > F Critical one-tail, which rejects the nullhypothesis: variances for the 2DOPS and 2DSE methods are equal. Thus, thevariance for the 2DOPS method is greater than that of the 2DSE method.8100.20.40.60.81−2 0 2 4 6 8 10 12 1400.20.40.60.81−4 −3 −2 −1 0 1 2 3 4−2 0 2 4 6 8 10 12 14 −4 −3 −2 −1 0 1 2 3 4 −6 −4 −2 0 2 4 6 8 10−6 −4 −2 0 2 4 6 8 10−2 0 2 4 6 8 1000.20.40.60.81−2 0 2 4 6 8 1000.20.40.60.81−5 −4 −3 −2 −1 0 1 2 3 4−5 −4 −3 −2 −1 0 1 2 3 41.21.21.21.2Distance (mm)Distance (mm)Distance (mm)Distance (mm) Distance (mm)Volunteer - I Volunteer - II Volunteer - IIIVolunteer - IV Volunteer - V2DOPS2DOPS2DSE2DSE(a)(b)(c)(d)(e)(g)(h)(f)(i)(j)Point DistributionPoint DistributionPoint DistributionPoint Distribution  0.17+__ 1.70 +__+__+__+__+_+__0.46 1.30 0.50 3.000.60 2.50 0.60 2.400.300.30+__+__+__0.01 0.30 1.00 1.370.55 0.86 0.40 1.30Figure 3.14: Distributions of the estimated bone boundary points by the 2DOPS and 2DSEmethods are plotted with respect to the expert delineated bone contours for volunteers-I to -V.The ‘mean ± standard-deviation’ of each distribution is shown inside the corresponding plots.Note that we use Gaussian mixture regression over the maximum intensitypoints in the fused map, while the bone boundaries produced by the 2DOPSmethod are found by bottom-up ray casting. To evaluate whether we have intro-duced a bias by comparing the result following regression with the 2DSE method82with the ray-cast points found in 2DOPS, we also calculated the mean absoluteerrors for the bone contours produced by applying Gaussian mixture regressionto the ray-cast points from the 2DOPS method. Fig. 3.15 shows that the MAEproduced by 2DSE is still approximately 28% lower than that produced by GMRapplied to the ray-cast points found in the 2DOPS method.Figure 3.15: Mean absolute fitting error (mm) comparisons of the bone contours produced byusing the Gaussian mixture regression over the ray-casted points for the 2DOPS method andthe contours produced by the original 2DOPS and proposed methods.3.7 Advantages and LimitationsIn this chapter, we improved our previous 2D bone segmentation approach (de-scribed in Chapter 2) by incorporating (i) a depth-dependent cumulative powerof the envelope into the elastographic data, (ii) an automatic weight estimatedfrom the echo de-correlation measure between the pre- and post-compression RFframes to fuse the strain and envelope power map, (iii) a local statistics-basedbone discontinuity detection scheme, and (iv) a non-parametric Gaussian mix-ture regression to produce smoother yet curvature preserved 2D bone contour.We demonstrated our improved performance on a wide range of validation dataincluding three sets of FEM data, an experimental phantom, and in vivo human83data. We achieved an improvements of approximately 19%-20% in the FEM, and18% in the experimental phantom in terms of MAE when compared with currentphase-feature-based state-of-the-art [57]. In addition, we achieved an improve-ments of approximately 55% in terms of mean absolute fitting error using the invivo data when compared with current phase-feature-based state-of-the-art [57].Although this study automates the bone segmentation approach that we de-veloped based on combined ultrasound strain-imaging and envelope signal power,however, this approach is limited to 2D images only. 2D methods do not takeadvantage of correlations between adjacent images (i.e., along the axis perpendic-ular to the scan plane direction) and are therefore subject to spatial compoundingerrors as well as errors due to beam thickness effects [56]. Therefore, our nextfocus is on developing a 3D bone segmentation approach by extending our 2Dmethod. A near-real time 3D elastography technique using 3D probe has beenpresented by Treece et al. [81] that can generate a 3D strain image, however, thistechnique has not been validated on clinical in vivo data and we found it to bevery sensitive to echo de-correlation while acquiring post-compression US volumewhich limits its applicability.84Chapter 4Strain-initialized Bone SurfaceDetection in 3D UltrasoundIn this chapter, we discuss our 3D bone surface extraction method using a surface-growing approach that is seeded from the 2D bone contours estimated from com-bined ultrasound strain-imaging and envelope power which we already describedin Chapter 3.4.1 Ultrasound Scanning ProtocolOur scanning protocol has two steps: (1) acquiring 2D pre- and post- compres-sion RF frames, I1 and I2, respectively, (of size m × n) where the US beamdirection is perpendicular to the bone surface, and (2) acquiring 3D US RF data(consist z number of frames) as the transducer sweeps over the target. To im-plement this protocol, we use a 4D linear array transducer that initially collectsa number of 2D US RF frames (axial×lateral) with free-hand compression whilethe transducer array is positioned in the middle of its sweep (i.e., in the probe’scentral plane). Once we are satisfied with the quality of the elastography data,we start acquiring the 3D US data by switching the transducer to the sweepingmode without removing it from the skin surface. The transducer sweeps in theelevation direction. In this step, the user can decrease the field-of-view anglealong the elevation direction if scanning a bone with a narrow cross-section (e.g.,85humerus, radius, ulna, etc.).4.2 3D Bone Surface Extraction Using Surface-growingWe modify the 3D B-mode (B) image prior to surface-growing in order to high-light the bone intensities more compared to that of soft tissue region. We es-timate the axially cumulative sum of intensities for each voxel as P (i, j, k) =∑ipo=1 B(po, j, k)2, where k denote the voxel index in the elevation direction. Us-ing P (i, j, k), we estimate the modified B-mode volume asBˆ(i, j, k) = B(i, j, k)P (i,j,k)max[B(i, j, k)P (i,j,k)] . (4.1)We assume that the 2D seed points bo(io, jo, ko) corresponds to voxels at (io, jo, z/2),where ko ≈ z/2. Then the bone surface expands frame-wise along the elevationdirection by minimizing the cost function defined asC(i, jo, ko ± 1) =jo+Ly∑jv=jo−Ly[Ei(i, jo, ko ± 1) + κEs(i, jo, ko ± 1)]× e−|jo−jv |, (4.2)where Ei and Es are the energy functions associated with the new candidatevoxel, κ is the mixing weight, and Ly is the lateral voxel number upon which theexponential weight spans to regularize C in the y direction. Energy Ei ensuresthat the new surface voxel at (jo, ko±1) carries the highest intensity value amongthe voxels in the interrogated scan-line, and is defined asEi(i, jo, ko ± 1)= 1− [Bˆ(i, jo, ko ± 1)− µt]2 − [Bˆ(i, jo, ko ± 1)− µb]2max[1− [Bˆ(i, jo, ko ± 1)− µt]2 − [Bˆ(i, jo, ko ± 1)− µb]2], (4.3)whereµt =1it − 1it−1∑a=1Bˆ(a, jo, ko ± 1), (4.4)86µb =1m− ibm∑a=ib+1Bˆ(a, jo, ko ± 1), (4.5)[it, ib] is the axial search range for a new surface candidate centered at io. Inaddition, the energy function Es ensures a smooth transition between the seedand new candidate voxel, and is defined asEs(i, jo, ko ± 1) =(|io − i|/|io − it|)2max[(|io − i|/|io − it|)2]. (4.6)Finally, we estimate the location of a new bone surface voxel from Eq. (4.2)for the frame at ko ± 1 asio |ko±1= argmini {C(i, jo, ko ± 1)}. (4.7)We update io |ko±f as the new seed voxel, proceed to the next frame at ko±(f+1),where f = 1, 2, 3.., and continue growing until min[C] rises to a value that makesit unlikely that the next voxel still represents bone. In our analyses, we set thisthreshold empirically at a value of 0.8; we plan to revisit this choice in future.4.3 Validation Setup4.3.1 Simulated PhantomWe built a 40mm×40mm×20mm (axial×lateral×elevation) FEM phantom usingthe analysis software ANSYS (ANSYS, Inc., Canonsburg, PA) and the ultrasoundsimulation software Field II [87]. The phantom mimicked a 3D US scan of thehuman distal radius bone with a total number of 2,207,200 nodes. The stiffnessof the homogeneous soft tissue and bone regions were set to 10kPa and 10GPa,respectively, as previously reported [88]. Our phantom was compressed in thelateral direction from the top using a planar compressor that was wider thanthe phantom. We set an applied pressure level that corresponds to 1% averagestrain. At first we acquired 2D elastographic data along lateral direction. Thenwe acquired 3D ultrasound data using a phased array transducer that swept boththe sides of the 2D frame along the elevation direction. An ultrasonic transducerof center frequency, f0 = 5MHz and band-width = 50% was used to simulate the87phantom scan from the top. The total number of scan lines was set to 128 and64 in the lateral and elevation directions, respectively. The resulting voxel size ofthe 3D data was 0.015mm×0.3mm×0.3mm. The B-mode image of the phantomis shown in Fig. 4.1(a).4.3.2 Experimental PhantomWe constructed an experimental phantom using a radio-opaque Sawbones hemi-pelvis (Sawbones, Pacific Research Laboratories, Inc., Vashon Island, WA), modelnumber 1297-22. A portion of the pelvis was suspended in polyvinyl chloride(PVC) gel (see Fig. 2.1(a)) and placed in an acrylic tube (see Fig. 2.1(b)). Ahigh resolution peripheral quantitative CT machine, model HR pQCT XtremeCT (Scanco USA, Inc., Wayne, PA) was used to acquire a single 482×482×402(lateral×axial×elevation) volume with a resolution of 0.25mm×0.25mm×0.25mm.The US images were acquired using a SonixRP (Ultrasonix Medical Corporation,Richmond, BC) scanner in the Center for Hip Health and Mobility, VancouverCoastal Health Authority, Vancouver, BC, Canada. We used a 4DL14-5/38 linear4D array probe operating at 5MHz for collecting data for the 3D implementations.4.3.3 In Vivo DataWe acquired three sets of in vivo data from three volunteers (volunteer-I: 27-yearold male; volunteer-II: 26-year old male; volunteer-III: 29-year old male) afterobtaining prior consent. All data were acquired with free-hand compression.The US images were acquired using a SonixRP (Ultrasonix Medical Corporation,Richmond, BC) scanner in the Center for Hip Health and Mobility, VancouverCoastal Health Authority, Vancouver, BC, Canada. We used a 4DL14-5/38 linear4D array probe operating at 5MHz for collecting data for the 3D implementations,respectively. The study was approved by the UBC clinical research ethics board.4.4 ResultsWe provide comparative results of our proposed 3D bone detection (3DSE) meth-ods with the previously reported volume-specific parameter optimized 3D local88(d)B-mode(a)3DOPS Ray-casting in 3DOPS 3DSE(b) (c)Figure 4.1: Illustration of the FEM results. (a) B-mode image, (b) estimated 3D phase symme-try image by the 3DOPS method, (c) bone surface after bottom up ray-casting in (b) (arrowsshow leakage), and (d) estimated bone surface by the proposed method.phase feature-based bone segmentation (3DOPS) [59] method using the FEMphantoms, experimental phantom and in vivo data. The 3DOPS method selectsparameters automatically for the Log-Gabor filters based on properties estimateddirectly from the specific US image that is being analyzed. We calculate meanabsolute error (MAE) which is defined asMAE = 1P ×QP∑p=1Q∑q=1| A(p, q)−G(p, q) |, (4.8)where P and Q represent the number of data points spanned by the bone bound-ary in the lateral and elevation directions, respectively, A is a matrix containingthe ‘ground truth’ bone boundary points, and G is a matrix containing the es-timated bone boundary points. As ground truth, we used the bone geometrydefined for Field II simulation for the FEM, the estimated bone boundary fromthe computed tomography (CT) image for the experimental phantom, and expertdelineation of bone boundaries on B-mode images for the in vivo data. We usefiducial-based CT and US bone alignment before estimating the MAE for theexperimental phantom. Note that we use bottom-up ray casting in the phase-symmetry images for the 2DOPS and 3DOPS methods to get G as suggestedin [59]. In our surface-growing technique, we set parameters as io − it ≈ 5mm,Ly = 3 and κ = 0.5.894.4.1 FEM ResultsWe provide comparative qualitative results of our proposed 3DSE and the 3DOPSmethods using the FEM data in Fig. 4.1. Note that the 3DOPS method usesbottom up ray-casting in the 3D phase symmetry image to produce final bonesurface; this can produce “leakage” through weak areas of the phase symmetryvolumes. As an evident, we see that the ray-casting leaks through weak spots inthe phase symmetry image (shown with white arrows in Fig. 4.1(c)). In contrast,the bone boundary produced by the proposed method (Fig. 4.1(d)) appears tobe free of such artifacts.(d)B-mode(a)3DOPS Ray-casting in 3DOPS 3DSE(b) (c)Figure 4.2: Illustration of the experimental phantom results. (a) B-mode image, (b) estimated3D phase symmetry image, (c) bone surface after bottom up ray-casting in (b) (arrows showleakage), and (d) estimated bone surface by the proposed method.4.4.2 Experimental Phantom ResultsFigure 4.2 demonstrates the performance comparison of the proposed 3DSE and3DOPS methods using the experimental phantom. Here also we see that theray-casting leaks through weak spots (indicated with white arrows in Fig. 4.1(c))in the 3D phase symmetry image (see Fig. 4.1(b)). In contrast, the bone surfaceestimated by the proposed method has no discontinuity (see Fig. 4.1(d)).4.4.3 In Vivo ResultsWe show the bone surface detection performance of the proposed 3DSE and3DOPS methods on three cases of volunteer in vivo data.904.4.3.1 Case 1B-mode 3DOPSRay-casting in 3DOPS 3DSEScanned Bone Area(a) (b) (c)(d) (e)Figure 4.3: Illustration of the 3D bone surface detection using the in vivo volunteer-I data set.(a) shows bone region scanned on the volunteer (with red rectangle), (b) is the B-mode image,(c) shows the estimated 3D phase symmetry image by the 3DOPS method, (d) shows the bonesurface after bottom up ray-casting in (c), and (e) shows the estimated bone surface by the3DSE method.We show the qualitative performance comparison of the 3DSE and 3DOPSmethods using the volunteer-I data set in Fig. 4.3. The scanned region, B-modeimage, 3D phase symmetry image, ray-casted version of the phase symmetryimage, and bone surface produced by the 3DSE method are shown in Figs. 4.3(a),(b), (c), (d), and (e), respectively. We see in Fig. 4.3(c) that the 3DOPS methodproduces false positive responses in the soft tissue interfaces (shown with white91arrows). We also see in Figs. 4.3(d) that the ray-casting results in leaks throughweak responses on the bone surface and instead detects false positive soft tissueinterfaces (shown with corresponding white arrows in Figs. 4.3(c) and (d)). Incontrast, the bone surfaces estimated by the 3DSE method better match theshapes visible in the corresponding B-mode images (see Figs. 4.3(b) and (e)).In addition, our method successfully extracts the bone surfaces that are notperpendicular to the US beam direction (shown with corresponding green arrowsin Figs. 4.3(b) and (e)). In contrast, the 3DOPS method fails to extract thesecurved bone regions.4.4.3.2 Case 2Figure 4.4 demonstrates the qualitative performance comparison of the 3DSE and3DOPS methods using the volunteer-II data set. The scanned region, B-modeimage, 3D phase symmetry image, ray-casted version of the phase symmetryimage, and bone surface produced by the 3DSE method are shown in Figs. 4.4(a),(b), (c), (d), and (e), respectively. We see in Fig. 4.4(c) that the 3DOPS methodproduces false positive responses in the soft tissue interfaces (shown with whitearrows). We also see in Figs. 4.4(d) that the ray-casting results in leaks throughweak responses on the bone surface and instead detects false positive soft tissueinterfaces (shown with corresponding white arrows in Figs. 4.4(c) and (d)). Incontrast, the bone surfaces estimated by the 3DSE method better match theshapes visible in the corresponding B-mode images (see Figs. 4.4(b) and (e)).4.4.3.3 Case 3Finally, in Fig. 4.5, we demonstrate the qualitative performance comparisonof the 3DSE and 3DOPS methods using the volunteer-II data set. The scannedregion, B-mode image, 3D phase symmetry image, ray-casted version of the phasesymmetry image, and bone surface produced by the 3DSE method are shown inFigs. 4.5(a), (b), (c), (d), and (e), respectively. We see in Fig. 4.5(c) that the3DOPS method produces false positive responses in the soft tissue interfaces(shown with white arrows). We also see in Figs. 4.5(d) that the ray-castingresults in leaks through weak responses on the bone surface and instead detects92B-mode 3DOPSRay-casting in 3DOPS 3DSEScanned Bone Area(a) (b) (c)(d) (e)Figure 4.4: Illustration of the 3D bone surface detection using the in vivo volunteer-I data set.(a) shows bone region scanned on the volunteer (with red rectangle), (b) is the B-mode image,(c) shows the estimated 3D phase symmetry image by the 3DOPS method, (d) shows the bonesurface after bottom up ray-casting in (c), and (e) shows the estimated bone surface by the3DSE method.false positive soft tissue interfaces (shown with corresponding white arrows inFigs. 4.5(c) and (d)). In contrast, the bone surfaces estimated by the 3DSEmethod better match the shapes visible in the corresponding B-mode images(see Figs. 4.5(b) and (e)). In addition, like case 1, our method successfullyextracts the bone surfaces that are not perpendicular to the US beam direction(shown with corresponding green arrows in Figs. 4.5(b) and (e)). In contrast,the 3DOPS method fails to extract these curved bone regions. However, very fewfalse positives are produced by the 3DSE method in this case as indicated with93B-mode 3DOPSRay-casting in 3DOPS 3DSEScanned Bone Area(a) (b) (c)(d) (e)Figure 4.5: Illustration of the 3D bone surface detection using the in vivo volunteer-I data set.(a) shows bone region scanned on the volunteer (with red rectangle), (b) is the B-mode image,(c) shows the estimated 3D phase symmetry image by the 3DOPS method, (d) shows the bonesurface after bottom up ray-casting in (c), and (e) shows the estimated bone surface by the3DSE method.the white arrow in Fig. 4.5(e).4.4.4 Quantitative ValidationWe estimate the MAE for the FEM and experimental phantom with four differentsignal-to-noise ratio (SNR) simulations (40, 30, 20, and 10dB) with 100 realiza-tions each. In real clinical settings, SNR of the RF data after beamforming isreported to be approximately 18∼22dB [86]. We also estimate the Matlab(R) exe-cution time for the in vivo data running in Intel(R) Xeon(R) CPU E3 @ 3.20GHz,9410 15 20 25 30 35 40SNR (dB)Experimental Phantom00.51.51.02.0  3DSE3DOPS(b)10 15 20 25 30 35 40MAE (mm)SNR (dB)FEM Phantom  3DSE3DOPS(a)00.51.51.02.0Figure 4.6: MAE analysis of the 3DSE and 3DOPS methods at different SNRs using (a) theFEM and (b) experimental phantom data.Table 4.1: Computation time in second for a volume size of 364×110×50 voxels.Volunteer 3DSE 3DOPSI 5.5 103.0II 5.9 105.6III 6.0 108.7Memory: 8GB in Table 4.1. Note that we use the bone surfaces produced by the3DOPS method after ray-casting in the MAE analysis after manually excludingleakage regions. Nonetheless, we see in Figs. 4.6(a) and (b) that for all SNRvalues, the mean MAE for the our method is lower than that of the 3DOPSmethod. In addition, the Table 4.1 shows that the execution times for the 3DSEmethod are approximately 18 times smaller than those of the 3DOPS method forall three in vivo data sets.We also show computation time diagrams in Fig. 4.7 that illustrate the runtimes of different components of the 3DOPS and 3DSE methods. Althougha GPU-based implementation of a version of the 3DOPS method is reportedwhich takes 0.8sec for its execution [79], this time refers solely to the core log-Gabor convolution and presumes that all parameters are already pre-selected.In the full 3DOPS method [59], all the parameters are automatically estimatedusing a process in which the core log-Gabor convolution, implemented in Mat-95lab(R), is performed iteratively over 50 times to select the optimum parameters.These repeated iterations are performed in the (a) filter scale selection, (b) fil-ter orientation refinement, and (c) filter angular band-width selection blocks inFig. 4.7(a). If the GPU-based implementation were incorporated in the current3DOPS method, the overall run-time, including automatic optimization, wouldbe approximately 1min, compared with ∼5s for the 3DSE method. If both com-ponents of the proposed 3DSE method were implemented in GPU, we anticipatethat we could reduce the run time to near real-time.Figure 4.7: Diagrams showing the computation time breakdown for different components of the(a) 3DOPS and (b) 3DSE methods.4.5 Advantages and LimitationsIn this chapter, we extend our 2D bone contour into a 3D bone surface by usinga surface-growing-based approach. We demonstrated our improved performanceon a wide range of validation data including a FEM, an experimental phantom,and in vivo human data with reported improvements of approximately 20% and23% in terms of MAE and 18 times in terms of computation time, respectively,when compared with current state-of-the-art.However, we used some empirical parameters in our 3D implementation, i.e.,io − it, Ly and κ. We plan to investigate this promising study further to makethese parameters data driven so that this method becomes fully automatic.96Chapter 5Conclusions5.1 The Impact and Potential Clinical Applica-tionWe have proposed novel methods for robust bone boundary localizations based onthe fusion of strain imaging and envelope signal power detection for 2D US imagesas well as surface-growing-based bone segmentation in 3D US volumes. We havecombined real-time strain imaging based on an analytic minimization of regular-ized cost functions with an envelope power map of the pre-compression RF frame.Our results have demonstrated reduced bone localization error through a betterelastogram resolution, adoption of an envelope power map of higher dynamicrange, addressing the false positive bone response, and exploiting the smoothingfeature of the AM method to better delineate the bone boundary in strain images.We have also developed an effective way to perform 3D bone surface extractionusing a surface growing approach that is seeded from 2D bone contours estimatedfrom combined ultrasound strain imaging and envelope power. We have demon-strated our improved performance on a wide range of validation data, includingthree sets of FEM data, an experimental phantom, and in vivo human data. Wehave achieved improvements of approximately 20% on both the 2D and 3D FEM,18% on the 2D experimental phantom, and 23% on the 3D experimental phantomin terms of MAE when compared with current phase-features-based method [57].In addition, we have achieved improvements of approximately 55% on the 2D in97vivo data in terms of mean absolute fitting error when compared with currentphase-features-based state-of-the-art technique [57]. We have also performed 3Dbone segmentation in in vivo data approximately 18 times faster than that of thecurrent state-of-the-art method [59].Our laboratory believes that 3D US has great potential to assist orthopaediccare, especially during surgeries, if the anatomical structures of interest can belocalized and visualized with sufficient accuracy and clarity and in a highly auto-mated and rapid manner [95]. Our laboratory has validated its 3D-local-phase-feature-based US segmentation technique on scans obtained from 29 trauma pa-tients with distal radius and pelvic ring fractures, and achieved an average surfacefitting error of 0.62±0.42mm for pelvic patients and 0.21±0.14mm for distal ra-dius patients. As we have demonstrated in this thesis, our method has achievedbetter accuracy in bone localizations than the state-of-the-art image-phase-basedmethod, and we hope that our method will likewise perform at a higher qualitylevel than the former in more challenging situations. Since our method is basedon the tissue stiffness estimation, it is expected to be robust in bone localizationin a challenging situation when ligament/tendon appears in between the trans-ducer face and bone surface. However, our method might produce false positivesin a rare situation when ligament is ossified and low strain profiles are producedalong the ligament-soft tissue layers similar to bone-soft tissue layers.Accordingly, we believe that our proposed method may have numerous clinicalapplications. The method can be used to assess bone fractures, diagnose bonedeformities, and guide computer-assisted orthopedic surgery. However, becauseit requires applying some pressure to the probe to generate the strain image, itwould be difficult to apply in situations in which the patient is awake and unableto tolerate pressure to an injury site.5.2 Thesis ContributionsThe goal of our work was to contribute some effective methods for robust bonesegmentation in ultrasound images that would ultimately address some of thechallenging limitations of the previously reported bone segmentation methods.This research was intended to further advance the larger goal of developing a98novel 3D ultrasound-based CAOS system. To this end, we have contributed thefollowing to this project:1. We have shown the potential of combined US strain-imaging and envelopepower signaling to delineate bone boundaries in US images. Our methoduses real-time strain imaging [82] based on an analytic minimization ofregularized cost functions to delineate bone from tissue stiffness maps and,thus, we have achieved a marked reduction in the number of of false positivebone responses at the soft tissue interface.2. We have introduced automatic parameter selection processes to make our2D bone segmentation method more robust and self-tuning. We have in-troduced the depth-dependent cumulative power of the envelope into thepre- and post-compression data, which we reveal to reduce the number ofsoft tissue related artifacts in the resulting strain images.3. We have incorporated a data-driven weight into our method, estimated fromthe echo de-correlation measurements as being between the pre- and post-compression RF frames, which controls the contribution of strain images inthe resulting fused map. We have also employed local statistics-based bonediscontinuity detection scheme.4. We have introduced Gaussian mixture regression, a more flexible multivari-ate non-parametric regression model that better preserves the curvaturefeatures at the bone boundary than the parametric linear regression model.5. We have developed a method for 3D bone surface extractions using a surfacegrowing approach that is seeded from the 2D bone contours estimated fromcombined ultrasound strain imaging and envelope power. The 2D seedbone contour grows into the 3D surface by minimizing a combined intensitysimilarity and a voxel proximity-based cost function.5.3 Remaining Challenges and Future WorkThe bone segmentation algorithm presented in this thesis is one of the compo-nents of an ultrasound guided CAOS system to assist surgeons during orthopedic99surgery. We have demonstrated in this thesis that the proposed segmentationmethods are accurate, rapid and robust. However, the slowest portion of themethod is, currently, the 3D segmentation runtime. Therefore, we plan to im-plement the 3D bone segmentation algorithm in C++ and run it on a graphicsprocessing unit (GPU) to accelerate the runtime. The current runtime for 3Dbone segmentation in MATLAB(R) is ∼5 seconds and we hope to achieve animprovement factor of at least 500x after the code is optimized. While we havedemonstrated excellent robustness and accuracy in our initial phantom and shortclinical study, there remains room for improvement to cope with the challengingclinical data sets.Although all of the parameters in our 2D bone segmentation method aredata-driven, our 3D bone segmentation method still has three free parameters:the search region io − it, the lateral regularization distance, Ly, and the weightof the fusing energy function, κ. In our future work, we plan to investigate thepossible approaches to having these parameters become automatically selected.Finally, we plan to perform two advanced validation tests of our proposedsegmentation methods before they are considered to be fully capable of beingincorporated into surgery environments. Thus,1. We plan to perform a laboratory-based simple fracture realignment pro-cedure, using a bovine long bone (e.g., femur) specimen, and assess theaccuracy of the bone localization as well as its fracture realignment perfor-mance.2. Latter, we plan to perform a surgical tool navigation performance procedureutilising a human/animal cadaver, and to assess the accuracy of the boneand surgical tool tip localization.In the long run, we hope to see this work employed in the ultrasound guidedCAOS system, with the goal of minimizing radiation time and the degree ofpossible surgical error during orthopedic surgeries.100Bibliography[1] P. Marhofer and V. W. Chan, “Ultrasound-guided regional anesthesia: cur-rent concepts and future trends,” Anesthesia & Analgesia, vol. 104, no. 5,pp. 1265–1269, 2007.[2] “X-ray (wikipedia),” http://en.wikipedia.org/wiki/X-ray, accessed: 2015-05-15.[3] “Nutek medical centre,” http://www.nutekmedicalcentremumbai.com/, ac-cessed: 2015-05-15.[4] “X-ray computed tomography (wikipedia),” http://en.wikipedia.org/wiki/X-ray computed tomography, accessed: 2015-05-15.[5] “Computerised axial tomography (CAT or CT),” http://www.schoolphysics.co.uk/age16-19/Atomic%20physics/X%20rays/text/CAT scanning/index.html, accessed: 2015-05-15.[6] “New paradigm in GPU-based monte carlo simulations for X-ray CT imagingdoses,” http://www.rpi.edu/dept/radsafe/public html/GPU project/indexGPU.html, accessed: 2015-05-15.[7] “Philips,” http://www.healthcare.philips.com/main/products/mri/systems/achieva3t/, accessed: 2015-05-15.[8] “Online image arcade,” http://imgarcade.com/1/magnetic-resonance-imaging-diagram/, accessed: 2015-05-15.[9] “Fluoroscopy & digital photospot,” http://www.upstate.edu/radiology/education/rsna/fluoro/fluoro.php, accessed: 2015-05-15.101[10] “Fluoroscopy (wikipedia),” http://en.wikipedia.org/wiki/Fluoroscopy, ac-cessed: 2015-05-15.[11] B. A. Schueler, “The AAPM/RSNA physics tutorial for residents generaloverview of fluoroscopic imaging,” Radiographics, vol. 20, no. 4, pp. 1115–1126, 2000.[12] “Color doppler ultrasound scanner,” http://www.lhdz.org/vet/ultrasound-op480.htm, accessed: 2015-05-15.[13] “MAB001 digital portable cheapest ultrasound scanner,” http://www.mayamed.cn/product/1298666144-213160187/MAB001 Digitalportable cheapest ultrasound scanner.html, accessed: 2015-05-15.[14] “How ultrasound imaging works explained simply,” http://www.howequipmentworks.com/ultrasound basics/, accessed: 2015-05-15.[15] C. Otto, “Principles of echocardiographic image acquisition and Doppleranalysis,” Otto CM: Textbook of clinical echocardiography. WB Saunders,Philadelphia, pp. 1–3, 2000.[16] “Overview of the smart materials technology,” http://resources.edb.gov.hk/physics/articleIE/smartmaterials/SmartMaterials e.htm, accessed: 2015-05-15.[17] “Basic principle of medical ultrasonic probes (transducer),” http://www.ndk.com/en/sensor/ultrasonic/basic02.html, accessed: 2015-05-15.[18] J. E. Aldrich, “Basic physics of ultrasound imaging,” Critical care medicine,vol. 35, no. 5, pp. S131–S137, 2007.[19] S. N. Narouze, Atlas of ultrasound-guided procedures in interventional painmanagement. Springer Science & Business Media, 2010.[20] “Ultrasound physics,” http://www.nysora.com/mobile/regional-anesthesia/foundations-of-us-guided-nerve-blocks-techniques/3084-ultrasound-physics.html, accessed: 2015-05-15.102[21] J. Bercoff, Ultrafast ultrasound imaging. INTECH Open Access Publisher,2011.[22] “Studyblue,” https://www.studyblue.com/notes/note/n/hubs/deck/7055337, accessed: 2015-05-29.[23] “What does this image of the hip display?” http://www.aium.org/soundwaves/article.aspx?aId=666&iId=20130912, accessed: 2015-05-15.[24] “Clinical implications of the echo enhancement artifact in volume sonogra-phy of the uterus,” http://www.jultrasoundmed.org/content/25/11/1431/F5.expansion, accessed: 2015-05-15.[25] “Artefacts on musculoskeletal ultrasound,” https://theultrasoundsite.co.uk/artefacts-musculoskeletal-ultrasound/, accessed: 2015-05-15.[26] “Pathology of the shoulder,” http://shoulderultrasound.com/973/1036.html, accessed: 2015-05-15.[27] “Basic physics of digital radiography/the image receptor,”http://en.wikibooks.org/wiki/Basic Physics of Digital Radiography/The Image Receptor, accessed: 2015-05-15.[28] W. Xu and S. T. Salcudean, “Enhancement of bone surface visualizationusing ultrasound radio-frequency signals,” in Ultrasonics Symposium, 2007.IEEE. IEEE, 2007, pp. 2535–2538.[29] M. A. Hussain, A. Hodgson, and R. Abugharbieh, “Robust bone detectionin ultrasound using combined strain imaging and envelope signal power de-tection,” in Medical Image Computing and Computer-Assisted Intervention–MICCAI 2014. Springer, 2014, pp. 356–363.[30] “How great is the burden of musculoskeletal disease?” http://www.aaos.org/news/aaosnow/jan08/youraaos15.asp, accessed: 2015-05-15.[31] K. Lefaivre, A. Starr, B. Barker, S. Overturf, and C. Reinert, “Early experi-ence with reduction of displaced disruption of the pelvic ring using a pelvic103reduction frame,” Journal of Bone & Joint Surgery, British Volume, vol. 91,no. 9, pp. 1201–1207, 2009.[32] L. P. Nolte and T. Beutler, “Basic principles of CAOS,” Injury, vol. 35,no. 1, pp. 6–16, 2004.[33] I. Hacihaliloglu, A. Brounstein, P. Guy, A. Hodgson, and R. Abugharbieh,“3D ultrasound-CT registration in orthopaedic trauma using GMM regis-tration with optimized particle simulation-based data reduction,” in Med-ical Image Computing and Computer-Assisted Intervention–MICCAI 2012.Springer, 2012, pp. 82–89.[34] R. A. Novelline and L. F. Squire, Squire’s fundamentals of radiology. LaEditorial, UPR, 2004.[35] G. T. Herman, Fundamentals of computerized tomography: image recon-struction from projections. Springer Science & Business Media, 2009.[36] K. A. Buckwalter, J. Rydberg, K. K. Kopecky, K. Crow, and E. L.Yang, “Musculoskeletal imaging with multislice CT,” American Journal ofRoentgenology, vol. 176, no. 4, pp. 979–986, 2001.[37] S. E. Polo and S. P. Jackson, “Dynamics of DNA damage response proteinsat DNA breaks: a focus on protein modifications,” Genes & development,vol. 25, no. 5, pp. 409–433, 2011.[38] J. M. Liebmann, “Ultrasound biomicroscopy of the eye,” Journal of Glau-coma, vol. 4, no. 4, p. 299, 1995.[39] I. Edler and K. Lindstro¨m, “The history of echocardiography,” Ultrasoundin medicine & biology, vol. 30, no. 12, pp. 1565–1644, 2004.[40] M. K. Hasan, M. A. Hussain, S. R. Ara, S. Y. Lee, and S. K. Alam, “Usingnearest neighbors for accurate estimation of ultrasonic attenuation in thespectral domain,” Ultrasonics, Ferroelectrics, and Frequency Control, IEEETransactions on, vol. 60, no. 6, pp. 1098–1114, 2013.104[41] A. Weyman, “Physical principles of ultrasound,” Principles and practice ofechocardiography, pp. 3–25, 1994.[42] G. Kossoff, “Basic physics and imaging characteristics of ultrasound,” Worldjournal of surgery, vol. 24, no. 2, pp. 134–142, 2000.[43] P. Kruizinga, F. Mastik, N. de Jong, A. F. van der Steen, and G. van Soest,“Plane-wave ultrasound beamforming using a nonuniform fast Fourier trans-form,” Ultrasonics, Ferroelectrics, and Frequency Control, IEEE Transac-tions on, vol. 59, no. 12, 2012.[44] R. S. Cobbold, Foundations of biomedical ultrasound. Oxford UniversityPress on Demand, 2007.[45] M. A. Hussain, S. K. Alam, S. Y. Lee, and M. K. Hasan, “Robust strain-estimation algorithm using combined radiofrequency and envelope cross-correlation with diffusion filtering,” Ultrasonic Imaging, vol. 34, no. 2, pp.93–109, 2012.[46] M. D. Kraus, G. Krischak, P. Keppler, F. T. Gebhard, and U. H. Schuetz,“Can computer-assisted surgery reduce the effective dose for spinal fu-sion and sacroiliac screw insertion?” Clinical Orthopaedics and RelatedResearch R©, vol. 468, no. 9, pp. 2419–2429, 2010.[47] U. Sto¨ckle, B. Ko¨nig, A. Scha¨ffler, T. Zschernack, and N. Haas, “Clinicalexperience with the siremobil Iso-C (3D) imaging system in pelvic surgery,”Der Unfallchirurg, vol. 109, no. 1, pp. 30–40, 2006.[48] U. Linsenmaier, C. Rock, E. Euler, S. Wirth, R. Brandl, D. Kotsianos,W. Mutschler, and K. J. Pfeifer, “Three-dimensional CT with a modified C-arm image intensifier: feasibility 1,” Radiology, vol. 224, no. 1, pp. 286–292,2002.[49] L. Nolte, H. Walti, H. Fujita, G. Zheng, W. Seissler, and J. Hey, “A novelmobile C-arm system for use in image-guided surgery: A feasibility study,”in Proc of the International Conference on Computer Assisted OrthopaedicSurgery, 2002.105[50] J. S. Hott, V. R. Deshmukh, J. D. Klopfenstein, V. K. Sonntag, C. A.Dickman, R. F. Spetzler, and S. M. Papadopoulos, “Intraoperative Iso-CC-arm navigation in craniospinal surgery: the first 60 cases,” Neurosurgery,vol. 54, no. 5, pp. 1131–1137, 2004.[51] I. D. Gelalis, N. K. Paschos, E. E. Pakos, A. N. Politis, C. M. Arnaoutoglou,A. C. Karageorgos, A. Ploumis, and T. A. Xenakis, “Accuracy of pediclescrew placement: a systematic review of prospective in vivo studies compar-ing free hand, fluoroscopy guidance and navigation techniques,” EuropeanSpine Journal, vol. 21, no. 2, pp. 247–255, 2012.[52] L. Carrat, J. Tonetti, P. Merloz, and J. Troccaza, “Percutaneous computerassisted iliosacral screwing: clinical validation,” in Medical Image Computingand Computer-Assisted Intervention–MICCAI 2000. Springer, 2000, pp.1229–1237.[53] B. Brendel, S. Winter, A. Rick, M. Stockheim, and H. Ermert, “Registra-tion of 3D CT and ultrasound datasets of the spine using bone structures,”Computer Aided Surgery, vol. 7, no. 3, pp. 146–155, 2002.[54] D. V. Amin, T. Kanade, A. M. Digioia, and B. Jaramaz, “Ultrasound regis-tration of the bone surface for surgical navigation,” Computer Aided Surgery,vol. 8, no. 1, pp. 1–16, 2003.[55] T. K. Chen, P. Abolmaesumi, D. R. Pichora, and R. E. Ellis, “A system forultrasound-guided computer-assisted orthopaedic surgery,” Computer AidedSurgery, vol. 10, no. 5-6, pp. 281–292, 2005.[56] I. Hacihaliloglu, R. Abugharbieh, A. Hodgson, and R. Rohling, “Bone seg-mentation and fracture detection in ultrasound using 3D local phase fea-tures,” in Medical Image Computing and Computer-Assisted Intervention–MICCAI 2008. Springer, 2008, pp. 287–295.[57] I. Hacihaliloglu, R. Abugharbieh, A. J. Hodgson, and R. N. Rohling, “Auto-matic adaptive parameterization in local phase feature-based bone segmen-tation in ultrasound,” Ultrasound in medicine & biology, vol. 37, no. 10, pp.1689–1703, 2011.106[58] A. Brounstein, I. Hacihaliloglu, P. Guy, A. Hodgson, and R. Abugharbieh,“Towards real-time 3D US to CT bone image registration using phase andcurvature feature based GMM matching,” in Medical Image Computing andComputer-Assisted Intervention–MICCAI 2011. Springer, 2011, pp. 235–242.[59] I. Hacihaliloglu, P. Guy, A. J. Hodgson, and R. Abugharbieh, “Volume-specific parameter optimization of 3D local phase features for improved ex-traction of bone surfaces in ultrasound,” The International Journal of Med-ical Robotics and Computer Assisted Surgery, vol. 10, no. 4, pp. 461–473,2014.[60] J. Tonetti, L. Carrat, S. Blendea, P. Merloz, J. Troccaz, S. Lavalle´e, andJ.-P. Chirossel, “Clinical results of percutaneous pelvic surgery: Computerassisted surgery using ultrasound compared to standard fluoroscopy,” Com-puter Aided Surgery, vol. 6, no. 4, pp. 204–211, 2001.[61] D. C. Barratt, G. P. Penney, C. S. Chan, M. Slomczykowski, T. J. Carter,P. J. Edwards, and D. J. Hawkes, “Self-calibrating 3D-ultrasound-based boneregistration for minimally invasive orthopedic surgery,” Medical Imaging,IEEE Transactions on, vol. 25, no. 3, pp. 312–323, 2006.[62] M. Beek, P. Abolmaesumi, S. Luenam, R. W. Sellens, and D. R. Pichora,“Ultrasound-guided percutaneous scaphoid pinning: Operator variabilityand comparison with traditional fluoroscopic procedure,” in Medical ImageComputing and Computer-Assisted Intervention–MICCAI 2006. Springer,2006, pp. 536–543.[63] M. Beek, P. Abolmaesumi, S. Luenam, R. E. Ellis, R. W. Sellens, and D. R.Pichora, “Validation of a new surgical procedure for percutaneous scaphoidfixation using intra-operative ultrasound,” Medical image analysis, vol. 12,no. 2, pp. 152–162, 2008.[64] J. A. Noble and D. Boukerroui, “Ultrasound image segmentation: a survey,”Medical Imaging, IEEE Transactions on, vol. 25, no. 8, pp. 987–1010, 2006.107[65] V. Daanen, J. Tonetti, and J. Troccaz, “A fully automated method for thedelineation of osseous interface in ultrasound images,” in Medical ImageComputing and Computer-Assisted Intervention–MICCAI 2004. Springer,2004, pp. 549–557.[66] J. Kowal, C. Amstutz, F. Langlotz, H. Talib, and M. G. Ballester, “Auto-mated bone contour detection in ultrasound B-mode images for minimally in-vasive registration in computer-assisted surgery An in vitro evaluation,” TheInternational Journal of Medical Robotics and Computer Assisted Surgery,vol. 3, no. 4, pp. 341–348, 2007.[67] A. Kryvanos, “Computer assisted surgery for fracture reduction and defor-mity correction of the pelvis and long bones,” 2003.[68] P. Foroughi, E. Boctor, M. J. Swartz, R. H. Taylor, and G. Fichtinger, “Ul-trasound bone segmentation using dynamic programming,” in UltrasonicsSymposium, 2007. IEEE. IEEE, 2007, pp. 2523–2526.[69] K. A. Patwardhan, K. Cao, D. Mills, and R. Thiele, “Automated bone andjoint-region segmentation in volumetric ultrasound,” in Biomedical Imaging(ISBI), 2012 9th IEEE International Symposium on. IEEE, 2012, pp. 1327–1330.[70] A. K. Jain and R. H. Taylor, “Understanding bone responses in B-modeultrasound images and automatic bone surface extraction using a bayesianprobabilistic framework,” in Medical Imaging 2004. International Societyfor Optics and Photonics, 2004, pp. 131–142.[71] P. He and J. Zheng, “Segmentation of tibia bone in ultrasound images usingactive shape models,” in Engineering in Medicine and Biology Society, 2001.Proceedings of the 23rd Annual International Conference of the IEEE, vol. 3.IEEE, 2001, pp. 2712–2715.[72] A. Alfiansyah, R. Streichenberger, P. Kilian, M. Bellemare, and O. Coulon,“Automatic segmentation of hip bone surface in ultrasound images using anactive contour,” International Journal of Computer Assisted Radiology andSurgery, vol. 1, p. 496, 2006.108[73] I. Hacihaliloglu, R. Abugharbieh, A. J. Hodgson, and P. Guy, “Automaticbone bontour detection in 3D B-mode ultrasound images using optimizedphase symmetry features - A clinical evaluation evaluation for pelvic frac-ture,” Journal of Bone & Joint Surgery, British Volume, vol. 94, no. SUPPXLIV, pp. 64–64, 2012.[74] B. Brendel, S. Winter, A. Rick, M. Stockheim, and H. Ermert, “Bone reg-istration with 3D CT and ultrasound data sets,” in International CongressSeries, vol. 1256. Elsevier, 2003, pp. 426–432.[75] G. Ionescu, S. Lavalle´e, and J. Demongeot, “Automated registration of ultra-sound with CT images: Application to computer assisted prostate radiother-apy and orthopedics,” in Medical Image Computing and Computer-AssistedIntervention–MICCAI99. Springer, 1999, pp. 768–777.[76] I. Hacihaliloglu, R. Abugharbieh, A. J. Hodgson, and R. N. Rohling, “En-hancement of bone surface visualization from 3D ultrasound based on localphase information,” in Ultrasonics Symposium, 2006. IEEE. IEEE, 2006,pp. 21–24.[77] I. Hacihaliloglu, R. Abugharbieh, A. Hodgson, and R. Rohling, “Automaticdata-driven parameterization for phase-based bone localization in US usinglog-gabor filters,” in Advances in Visual Computing. Springer, 2009, pp.944–954.[78] P. Kovesi, “Image features from phase congruency,” Videre: Journal of com-puter vision research, vol. 1, no. 3, pp. 1–26, 1999.[79] A. Amir-Khalili, R. Abugharbieh, and A. J. Hodgson, “Using graphicsprocessing units to achieve real-time bone surface extraction from volumetricmedical ultrasound image data using local phase features,” in ComputerAssisted Orthopaedic Surgery (CAOS), 2013, conference. [Online]. Available:http://bisicl.ece.ubc.ca/papers/CAOS2013.pdf[80] T. Varghese, J. Ophir, E. Konofagou, F. Kallel, and R. Righetti, “Tradeoffsin elastographic imaging,” Ultrasonic Imaging, vol. 23, no. 4, pp. 216–248,2001.109[81] G. M. Treece, J. E. Lindop, A. H. Gee, and R. W. Prager, “Freehand ul-trasound elastography with a 3D probe,” Ultrasound in medicine & biology,vol. 34, no. 3, pp. 463–474, 2008.[82] H. Rivaz, E. M. Boctor, M. A. Choti, and G. D. Hager, “Real-time regu-larized ultrasound elastography,” Medical Imaging, IEEE Transactions on,vol. 30, no. 4, pp. 928–945, 2011.[83] M. A. Hussain, E. M. Abu Anas, S. K. Alam, S. Y. Lee, and M. K. Hasan,“Direct and gradient-based average strain estimation by using weighted near-est neighbor cross-correlation peaks,” Ultrasonics, Ferroelectrics, and Fre-quency Control, IEEE Transactions on, vol. 59, no. 8, pp. 1713–1728, 2012.[84] J. P. Spalazzi, J. Gallina, S. D. Fung-Kee-Fung, E. E. Konofagou, and H. H.Lu, “Elastographic imaging of strain distribution in the anterior cruciate lig-ament and at the ligament-bone insertions,” Journal of orthopaedic research,vol. 24, no. 10, pp. 2001–2010, 2006.[85] H. Rivaz, E. Boctor, P. Foroughi, R. Zellars, G. Fichtinger, and G. Hager,“Ultrasound elastography: a dynamic programming approach,” MedicalImaging, IEEE Transactions on, vol. 27, no. 10, pp. 1373–1377, 2008.[86] W. Boonleelakul, U. Techavipoo, D. Worasawate, R. Keinprasit, P. Pinun-sottikul, N. Sugino, and P. Thajchayapong, “Compression of ultrasound RFdata using quantization and decimation,” in Biomedical Engineering Inter-national Conference (BMEiCON), 2013 6th. IEEE, 2013, pp. 1–4.[87] J. A. Jensen, “Field: A program for simulating ultrasound systems,” in 10thNordicbaltic conference on Biomedical Imaging, PART 1, vol. 4. Citeseer,1996, pp. 351–353.[88] W. Pistoia, B. Van Rietbergen, E.-M. Lochmu¨ller, C. Lill, F. Eckstein, andP. Ru¨egsegger, “Estimation of distal radius failure load with micro-finiteelement analysis models based on three-dimensional peripheral quantitativecomputed tomography images,” Bone, vol. 30, no. 6, pp. 842–848, 2002.110[89] R. Zahiri-Azar and S. E. Salcudean, “Motion estimation in ultrasound im-ages using time domain cross correlation with prior estimates,” BiomedicalEngineering, IEEE Transactions on, vol. 53, no. 10, pp. 1990–2000, 2006.[90] R. Nass, L. S. Farhy, J. Liu, S. S. Pezzoli, M. L. Johnson, B. D. Gaylinn,and M. O. Thorner, “Age-dependent decline in acyl-ghrelin concentrationsand reduced association of acyl-ghrelin and growth hormone in healthy olderadults,” The Journal of Clinical Endocrinology & Metabolism, vol. 99, no. 2,pp. 602–608, 2013.[91] F. Kallel and J. Ophir, “A least-squares strain estimator for elastography,”Ultrasonic imaging, vol. 19, no. 3, pp. 195–208, 1997.[92] H. G. Sung, “Gaussian mixture regression and classification,” Ph.D. disser-tation, Rice University, 2004.[93] L. Zelnik-Manor and P. Perona, “Self-tuning spectral clustering,” in Ad-vances in neural information processing systems, 2004, pp. 1601–1608.[94] D. A. Cohn, Z. Ghahramani, and M. I. Jordan, “Active learning with sta-tistical models,” Journal of artificial intelligence research, 1996.[95] I. Hacihaliloglu, P. Guy, A. J. Hodgson, and R. Abugharbieh, “Automaticextraction of bone surfaces from 3D ultrasound images in orthopaedic traumacases,” nternational Journal of Computer Assisted Radiology and Surgery(IJCARS), 2015.111Appendix AClinical Evaluation Protocol Form112113114115116117118119Appendix BClinical Evaluation Patient ConsentForm120121122123124125126127

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0166292/manifest

Comment

Related Items