Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Registration of preoperative CT to intraoperative ultrasound via a statistical wrist model for scaphoid… Anas, Emran Mohammad Abu 2016

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

Item Metadata

Download

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

Full Text

Registration of preoperative CT tointraoperative ultrasound via astatistical wrist model for scaphoidfracture fixationbyEmran Mohammad Abu AnasB.Sc., Bangladesh University of Engineering and Technology, 2009M.Sc., Bangladesh University of Engineering and Technology, 2011A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Electrical and Computer Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)December 2016c© Emran Mohammad Abu Anas 2016AbstractScaphoid fracture is the most probable outcome of wrist injury and it oftenoccurs due to sudden fall on an outstretched arm. To fix an acute non-displaced fracture, a volar percutaneous surgical procedure is highly recom-mended as it provides faster healing and better biomechanical outcome tothe recovered wrist. Conventionally, this surgical procedure is performedunder X-ray based fluoroscopic guidance, where surgeons need to mentallydetermine a trajectory of the drilling path based on a series of 2D projectionimages. In addition to challenges associated with mapping 2D informationto a 3D space, the process involves exposure to ionizing radiation. Ultra-sound (US) has been suggested as an alternate; US has many advantagesincluding its non-ionizing nature and real-time 3D acquisition capability.US images are, however, difficult to interpret as they are often corruptedby significant amounts of noise or artifact, in addition, the appearance ofthe bone surfaces in an US image contains only a limited view of the truesurfaces.In this thesis, I propose techniques to enable ultrasound guidance inscaphoid fracture fixation by augmenting intraoperative US images withpreoperative computed tomography (CT) images via a statistical anatomicalmodel of the wrist. One of the major contributions is the development ofa multi-object statistical wrist shape+scale+pose model from a group ofsubjects at wide range of wrist positions. The developed model is thenused to register with the preoperative CT to obtain the shapes and sizes ofthe wrist bones. The intraoperative procedure starts with a novel US boneenhancement technique that takes advantage of an adaptive wavelet filterbank to accurately highlight the bone responses in US. The improved boneenhancement in turn enables a registration of the statistical pose model tointraoperative US to estimate the optimal scaphoid screw axis for guidingthe surgical procedure. In addition to this sequential registration technique,I propose a joint registration technique that allows a simultaneous fusion ofthe US and CT data for an improved registration output. We conduct acadaver experiment to determine the accuracy of the registration process,and compare the results with the ground truth.iiPrefaceThis thesis primarily consists of four published manuscripts, resulting fromthe collaboration between multiple researchers. All these publications havebeen modified accordingly to present a consistent thesis. Ethical approvalfor conducting this research has been provided by the Clinical ResearchEthics Board, certificate numbers: H13-01891 and 6007480.A version of Chapters 1 and 5 has been accepted for publication in:• Emran Mohammad Abu Anas, Alexander Seitel, Abtin Rasoulian,Paul St John, Tamas Ungi, Andras Lasso, Kathryn Darras, DavidWilson, Victoria A. Lessoway, Gabor Fichtinger, Michelle Zec, DavidPichora, Parvin Mousavi, Robert N. Rohling and Purang Abolmae-sumi. Registration of a statistical model to intraoperative ultrasoundfor scaphoid screw fixation, International Journal of Computer As-sisted Radiology and Surgery, 1-9, 2016.The author contributed by developing, implementing and evaluating themethod. Dr. Seitel and Dr. Rasoulian helped in developing and imple-menting the method. Dr. Seitel, P. St John, Dr. Ungi, Dr. Lasso, Prof.Fichtinger, Dr. Zec, Dr. Pichora and Prof. Mousavi were involved in con-ducting the cadaver experiment used in this study. Prof. Wilson and K.Darras helped in acquiring and analyzing the CT images. V. A. Lessowayhelped in developing the ultrasound scanning protocol. Profs. Mousavi,Rohling and Abolmaesumi provided important feedback and suggestions inimproving the methodology.All co-authors contributed to the editing of the manuscript.A version of Chapters 2 and 3 has been accepted for publication in:• Emran Mohammad Abu Anas, Abtin Rasoulian, Alexander Seitel,Kathryn Darras, David Wilson, Paul St John, David Pichora, ParvinMousavi, Robert N. Rohling and Purang Abolmaesumi. AutomaticSegmentation of Wrist Bones in CT Using a Statistical Wrist Shape+PoseiiiPrefaceModel, IEEE Transactions on Medical Imaging, 13(11): 2025-2034,2016.The author contributed by developing, implementing and evaluating themethod. Dr. Seitel and Dr. Rasoulian helped in developing and imple-menting the method. P. St John, Dr. Pichora, Prof. Wilson and K. Darrashelped in acquiring and analyzing the CT images. Profs. Mousavi, Rohlingand Abolmaesumi provided important feedback and suggestions in improv-ing the methodology.All co-authors contributed to the editing of the manuscript.A version of Chapters 2 and 3 has been accepted for publication in:• Emran Mohammad Abu Anas, Abtin Rasoulian, Paul St John, DavidPichora, Robert N. Rohling and Purang Abolmaesumi. A Statisti-cal Shape+Pose Model for Segmentation of Wrist CT Images, SPIEMedical Imaging, 90340T, 2014.The author contributed by developing, implementing and evaluating themethod. Dr. Rasoulian helped in developing and implementing the method.P. St John and Dr. Pichora helped in acquiring and analyzing the CT im-ages. Profs. Rohling and Abolmaesumi provided important feedback andsuggestions in improving the methodology.All co-authors contributed to the editing of the manuscript.A version of Chapter 4 has been accepted for publication in:• Emran Mohammad Abu Anas, Alexander Seitel, Abtin Rasoulian,Paul St John, David Pichora, Kathryn Darras, David Wilson, VictoriaA. Lessoway, Ilker Hacihaliloglu, Parvin Mousavi, Robert N. Rohlingand Purang Abolmaesumi. Bone enhancement in ultrasound usinglocal spectrum variations for guiding percutaneous scaphoid fracturefixation procedures, International Journal of Computer Assisted Ra-diology and Surgery, 10(6):959-969, 2015.The author contributed by developing, implementing and evaluating themethod. Prof. Hacihaliloglu, Dr. Seitel and Dr. Rasoulian helped indeveloping and implementing the method. P. St John and Dr. Pichorawere involved in phantom development used in this study. Prof. Wilsonand K. Darras helped by providing idea in scaphoid fracture fixation systemdevelopment. V. A. Lessoway helped in developing the ultrasound scanningprotocol. Profs. Mousavi, Rohling and Abolmaesumi provided importantivPrefacefeedback and suggestions in improving the methodology.All co-authors contributed to the editing of the manuscript.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Clinical background . . . . . . . . . . . . . . . . . . . . . . . 11.1.1 Diagnosis . . . . . . . . . . . . . . . . . . . . . . . . . 11.1.2 Types . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.1.3 Treatment . . . . . . . . . . . . . . . . . . . . . . . . 21.2 US guidance in computer assisted interventions . . . . . . . . 51.2.1 CT+US fusion . . . . . . . . . . . . . . . . . . . . . . 51.2.2 SM+US fusion . . . . . . . . . . . . . . . . . . . . . . 71.3 Objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71.4 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . 101.5 Thesis outline . . . . . . . . . . . . . . . . . . . . . . . . . . 122 Novel statistical shape+scale+pose model for wrist . . . . 142.1 Introduction and background . . . . . . . . . . . . . . . . . . 142.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.2.1 Digital wrist database . . . . . . . . . . . . . . . . . . 162.2.2 Development of statistical shape model . . . . . . . . 172.2.3 Statistical scale model . . . . . . . . . . . . . . . . . . 192.2.4 Development of statistical pose model . . . . . . . . . 212.2.5 Model representation . . . . . . . . . . . . . . . . . . 23viTable of Contents2.2.6 Model registration . . . . . . . . . . . . . . . . . . . . 242.3 Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 242.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 252.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 262.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 293 Automatic segmentation of wrist bones using statistical modelregistration to CT . . . . . . . . . . . . . . . . . . . . . . . . . 313.1 Introduction and background . . . . . . . . . . . . . . . . . . 313.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 333.2.1 Initial alignment . . . . . . . . . . . . . . . . . . . . . 333.2.2 Model registration . . . . . . . . . . . . . . . . . . . . 353.2.3 Non-rigid registration . . . . . . . . . . . . . . . . . . 363.3 Experiments and evaluation . . . . . . . . . . . . . . . . . . 363.3.1 Evaluating parameters . . . . . . . . . . . . . . . . . 373.3.2 Parameter selection . . . . . . . . . . . . . . . . . . . 373.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 383.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 403.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 454 Bone enhancement in US using local spectrum variation . 464.1 Introduction and background . . . . . . . . . . . . . . . . . . 464.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 484.2.1 Spectrum variation of US bone responses . . . . . . . 484.2.2 Phase symmetry estimation . . . . . . . . . . . . . . 494.2.3 Bone enhancement . . . . . . . . . . . . . . . . . . . 524.3 Experiments and evaluation . . . . . . . . . . . . . . . . . . 574.3.1 Phantom . . . . . . . . . . . . . . . . . . . . . . . . . 574.3.2 Ex vivo . . . . . . . . . . . . . . . . . . . . . . . . . . 584.3.3 In vivo . . . . . . . . . . . . . . . . . . . . . . . . . . 584.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 594.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 634.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 635 Scaphoid screw axis estimation based on US guidance . . 655.1 Introduction and background . . . . . . . . . . . . . . . . . . 655.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 675.2.1 Initial alignment . . . . . . . . . . . . . . . . . . . . . 675.2.2 Pose optimization . . . . . . . . . . . . . . . . . . . . 675.2.3 Scaphoid screw axis estimation . . . . . . . . . . . . . 68viiTable of Contents5.3 Experiments and evaluation . . . . . . . . . . . . . . . . . . 695.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 725.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 735.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 766 Joint atlas-based registration of multiple bones in CT andUS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 776.1 Introduction and background . . . . . . . . . . . . . . . . . . 776.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 786.3 Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 806.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 806.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 826.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 837 Conclusion and future work . . . . . . . . . . . . . . . . . . . 847.1 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . 847.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . 86Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88viiiList of Tables1.1 Summary of CT+US-based methods. The top and bottomparts of the table present the intensity- and volume-basedapproaches, respectively. . . . . . . . . . . . . . . . . . . . . . 81.2 Summary of SM+US-based methods. The top and bottomparts of the table present the intensity- and volume-basedapproaches, respectively. . . . . . . . . . . . . . . . . . . . . . 92.1 Average and standard deviation of the mSDEs achieved byregistering the model to surfaces of the training data in theLOSO experiment. At each wrist position, the mSDE is com-puted for each training sample, and the average mSDE iscomputed across all training samples at that particular wristposition. Note that the sample size is different at differentwrist positions. Flexion (F), Extension (E), ulnar deviation(UD) and radial deviation (RD). . . . . . . . . . . . . . . . . 282.2 Average and standard deviation of the mSDEs achieved byregistering the model to surfaces of the training data in theLOPO experiment. The sample sizes for this table are sameas those in Table 2.1 . . . . . . . . . . . . . . . . . . . . . . . 293.1 Statistical measures of the mSDEs and JIs of nine wrist bonescomputed from 66 CT images (Radius= Rad, Scaphoid= Sca,Capitate= Cap, Hamate= Ham, Lunate= Lun, Pisiform=Pis, Trapezium= Trzm, Trapezoid= Trzd, Triquetrum= Triq.) 425.1 Results of the proposed approach with respect to some definedcriterion. ‘MD’, ‘VA’, ‘Y’ and ‘N’ indicate minimum distance,volar accessibility, yes and no answers, respectively. . . . . . 736.1 Results of the proposed approach with respect to some definedcriterion. ‘MD’, ‘VA’, ‘Y’ and ‘N’ indicate minimum distance,volar accessibility, yes and no answers, respectively. . . . . . 81ixList of Figures1.1 A schematic diagram showing the wrist bones of an examplewrist. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2 A plain X-ray of a fractured wrist (from Wikimedia Commons). 31.3 Different types of scaphoid fractures. . . . . . . . . . . . . . . 31.4 An example demonstration of the scaphoid bone response inUS. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61.5 The proposed workflow for scaphoid screw axis estimation.Pre-operatively, the statistical wrist model is registered to thepreoperative CT at a neutral wrist position for estimation ofthe size and shape of the wrist bones. Intra-operatively, a3D US volume is acquired at extension wrist position andthe wrist model is registered to the bone surfaces extractedfrom the US volume. In this step, only pose optimizationis performed and shape and scale parameters are applied ascalculated from the CT volume. Finally, the scaphoid screwposition is obtained by maximizing the screw length within asafe zone inside the registered scaphoid bone. . . . . . . . . . 101.6 The proposed workflow for joint registration approach. In-stead of sequential registration, the statistical wrist model isjointly registered to US and CT. Note that the shape andscale of the wrist bones are determined from the joint opti-mization. Poses of the bones are estimated indepedently tocapture the pose differences across the US and CT acquisitions. 11xList of Figures2.1 A schematic diagram showing the example wrist alignment,mean shape generation and similarity transformation calcu-lation. Correspondences among the wrist examples are es-tablished using a group-wise registration-based framework.Subsequently, the wrist examples are aligned in a commoncoordinate and mean shape is generated across all exampleshapes for each wrist bone. Similarity transformation Tn,l iscomputed for l-th bone of n-th example from the mean shapeto the l-th bone of the n-th example. . . . . . . . . . . . . . . 192.2 A schematic diagram showing the process of developing thetwo inter- and intra-subject pose models. I concatenate thepose features from all wrist bones at different wrist positionsto form a joint pose feature. Initially, PCA extracts the inter-subject pose variations. After that an instance created fromthe inter-subject pose model is used to develop the intra-subject pose model that represents the pose variations acrossdifferent wrist positions. . . . . . . . . . . . . . . . . . . . . . 202.3 Compactness of the developed (a) shape, (b) scale, (c) inter-subject pose and (d) intra-subject pose models. . . . . . . . . 262.4 Variations of the wrist model by changing the first mode ofthe (a) shape, (b) scale, (c) inter-subject pose and (d) intra-subject pose models. In each figure, the white and greensurfaces represent two model instances at extreme values ofthe corresponding mode. . . . . . . . . . . . . . . . . . . . . . 273.1 Flow chart of whole registration procedure. The CT volumeis thresholded to create a rough wrist structure. It is usedto initially align the wrist model. An iterative approach isfollowed to refine the wrist structure and estimate the modelcoefficients. In the end, non-rigid registration approach isused to correct the residual errors. . . . . . . . . . . . . . . . 34xiList of Figures3.2 Influence of the number of repetitions (JR) and the numberof points (N) in the target point cloud on the registration ac-curacy and computation time. Note that when I vary JR tocompute the affected indices, I keep N constant at 20000 andvice versa I keep JR constant at 6. The mSDE is computedacross all vertices of the registered model. For each JR andN , mSDEs and JIs between the resulting wrist models andthe reference segmentations are averaged over all 15 CT vol-umes. In contrast, the computation time is computed for anarbitrary CT volume for each JR and N . (a) Average mSDEand JI vs. JR. (b) Computation time vs. JR. (c) AveragemSDE and JI vs. N . (d) Computation time vs. N . . . . . . . 393.3 Dependence of the mSDE of the average voxel size in theCT volumes. The mean SDE is computed across all verticesof the registered model to a CT, and the average voxel sizeis calculated as the mean of the voxel dimensions along x-, y- and z-axes in that CT volume. Results from 66 CTvolumes are used in this figure and each point represents thecoordinate of (average voxel size, mSDE) for a particular CTvolume. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 403.4 Segmented wrist bone boundaries superimposed on exampleCT slices at different planes. Each color represents a partic-ular bone and their correspondences are given in the figure.(a-c) Example segmented axial, sagittal and coronal CT slicesof a healthy individual at neutral wrist position. (d-e) Twoexample axial CT slices showing segmentations of fracturedtrapezium (d) and lunate (e) bones. (f) Example segmentedsagittal CT slice of a cadaver at extension position. . . . . . . 413.5 Effect of the non-rigid registration step on the registrationerror. I compute mSDE across all vertices of each registeredsurface. Two types of segmentations are performed here. Thefirst one is my proposed approach and in the second one Iavoid the non-rigid registration step from my proposed seg-mentation technique. The median value of the mSDE across66 CT samples is marked by the red color. The solid bluerepresents the range between 25 and 75 percentiles. The redcrosses indicate outliers. . . . . . . . . . . . . . . . . . . . . . 43xiiList of Figures3.6 Error distributions across all vertices of the registered model.Segmented bone surfaces of all wrist bones of an arbitrary CTvolume at palmar (a) and dorsal (b) views. The color rep-resents the error in mm of the segmentation with respect togold standard. (c-e) Axial, sagittal and coronal slices demon-strating the errors in the corresponding CT planes marked in(a). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 444.1 Simulated US bone responses and their spectrums. (a-c)Three simulated bone responses, (d-f) corresponding 2D Fouriertransformations. . . . . . . . . . . . . . . . . . . . . . . . . . 504.2 Orientations in 2D and 3D frequency spectrums. (a) An ex-ample orientation in 2D frequency spectrum. (b) A 3D fre-quency spectrum is divided into different cones, and each seg-mented cone is further partitioned into different orientations. 514.3 The flowchart of the proposed 2D bone enhancement ap-proach. (1) The corresponding Fourier transform P (ω, φ) inthe polar domain is computed from the US image. (2) Esti-mation of the orientational filter parameters and orientationalfilter design. (3) Log-Gabor filter design. (4) The local PSis estimated from the designed filter and the US image. (5)After that, I obtain BR image from PS using the approachproposed by Foroughi et al. [25]. . . . . . . . . . . . . . . . . 524.4 Utilization of the spectrum variation in the filter design. (a)The variation of spectrum strength over the polar angle. Themaxima and neighboring two minima are detected from thespectrum strength. (b) The variation of spectrum strengthover the angular frequency. The spectrum is divided intodifferent parts and the lower and upper cut-off frequenciesare determined from each part. . . . . . . . . . . . . . . . . . 544.5 Experimental setup of the proposed approach. (a) A cadaverand (b) a volunteer wrist are kept fixed in our custom-builtwrist holder. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 584.6 Comparative bone enhancement results in the phantom ex-periment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 604.7 Results of the proposed, 2DLPS and 3DLPS methods for theex vivo experiment. Example axial and sagittal US framesare shown in (a), (f). The corresponding bone enhancementsare demonstrated below. The differences in the enhancementare prominent in the surfaces marked by arrows. . . . . . . . 61xiiiList of Figures4.8 Results of the proposed, 2DLPS and 3DLPS methods for thein vivo experiment. Example axial and sagittal US framesare shown in (a), (f). The corresponding bone enhancementsare demonstrated below. The differences in the enhancementare prominent in the surfaces marked by arrows. . . . . . . . 625.1 An example of manual landmark selection across the wristmodel instance (a) and the extracted bone surface in US (b). 685.2 Evaluation of the proposed approach. (a) Effect of the UStransducer’s pressure on the fiducial-based registration of man-ually segmented CT to US. The registered segmented CT(marked by white) is overlaid on an US image, and the regis-tered surface does not align with the bone responses in US. Tominimize the effect of the US transducer’s pressure, I manu-ally correct the registration and the manually registered seg-mented CT (marked by yellow) corresponds better to the USbone responses. (b) Several evaluation criteria are measuredfrom the estimated screw position and the reference wristbone surfaces. . . . . . . . . . . . . . . . . . . . . . . . . . . . 705.3 Example results of the proposed approach. (a) Example reg-istration result of the statistical model to US. (b) Screw tra-jectory estimated by my method. The trajectory is not ob-structed by the surrounding trapezium bone. (c) An exam-ple of a registered model overlaid on an US frame. I noticea misalignment between the registered model and the boneresponse marked by a dashed arrow. The source of this er-ror mainly includes the inaccurate segmentation in CT im-ages. (d) Segmentation result of corresponding preoperativeCT. The mentioned inaccurate segmentation is marked by adashed arrow. . . . . . . . . . . . . . . . . . . . . . . . . . . . 745.4 Possible sources of errors. (a) A limited field of view for thescaphoid bone. (b) Minimal bone response is observed for thetrapezium bone. The minimal bone responses mainly occurdue to non-perpendicular scanning of the trapezium bone. . . 756.1 A work-flow of the proposed joint registration approach. . . . 78xivList of Figures6.2 Results of the proposed joint registration approach and com-parison with the sequential registration technique for the scaphoidbone. (a-b) Example US registration results for the proposedjoint and sequential registration methods, respectively. (c)Corresponding CT slice to show the source of error of thesequential registration method. . . . . . . . . . . . . . . . . . 81xvAcknowledgementsI would like to express my sincere gratitude to my supervisor Prof. Pu-rang Abolmaesumi for his continuous support during my PhD study andrelated research. I always enjoy working under his guidance in all the timeof research and writing of this thesis. I also feel fortunate and proud to getopportunity to work with my supervisor. In addition to an excellent mentorand advisor, my supervisor is an excellent human being who always asks meabout the difficulties in the abroad life. I also would like to express my kindregards to my supervisor for his recommendations that helped me to winthe GSI awards for four times in my PhD career.My sincere thanks also goes to Prof. Robert Rohling not only for his consis-tent support in my PhD life, but also for steering me in the right direction.I also thank Robert for assistance with the methodology and for commentsthat greatly improved the manuscript.I would also like to thank Prof. Parvin Mousavi for her extensive supportin method development and manuscript preparation.Special thanks to Alexander Seitel for his insightful feedbacks and sharinghis knowledge in image analysis, and to Abtin Rasoulian for helping withprogrammming.I am very grateful to Paul St John, Tamas Ungi, Andras Lasso, Alireza,Gabor Fichtinger, Michelle Zec, David Pichora and Prof. Parvin Mousaviat Queen’s University and Kingston General Hospital for conducting a totalof 13 cadaver experiments. Without their precious support it would not bepossible to conduct this research.I am also thankful to Vickie Lessoway, Kathryn Darras, David Wilson andIlker Hacihaliloglu who were far more than collaboration partners.I am also thankful to all the former and current members of the RCL labat UBC: Mehran, Delaram, Manyou, Shekoofeh, Mehdi, Siavash, Saman,Amin, Mohammad Najafi, Phillipe and Hussam. I also wish to thank allother friends who shared their joy with me.I thankfully acknowledge the fellowships and grants that supported my re-search during the course of my studies: the University of British Columbia,the National Science and Engineering Council of Canada (NSERC), and thexviAcknowledgementsCanadian Institutes of Health Research (CIHR).I offer my gratitude to my beloved parents for providing me unconditionalsupport during my PhD career. Last but not least I would like to thankmy wife, Sabrina, who has been with me in every step with her support,encouragement, quiet patience and unwavering love.xviiChapter 1Introduction1.1 Clinical backgroundThe wrist injury is one of the most common joint injuries, often occursdue to sudden fall on an outstretched hand. A statistical analysis on theincidents of wrist injuries indicates more than 250,000 reported cases eachyear only in North America [62, 70]. Wrist injury can occur to any bone,however, the scaphoid bone of roughly 25×10×10 mm in an adult, is themost frequently (about 90%) fractured bone due to its special anatomicalposition in the wrist [64].The wrist anatomy consists (Fig. 1.1) of eight small intricately shapedcarpal bones, two forearm bones (radius and ulna) and five metacarpalbones. The carpal bones are arranged in two rows, the lower row is calledas proximal row and the upper one is known as distal one. The four carpalbones (scaphoid, lunate, pisiform, triquetrum) in the proximal row connectto the forearm, and the other four carpal bones (trapezium, trapezoid, cap-itate, hamate) in the distal row link to the metacarpal bones. Among all ofthe carpal bones, only scaphoid bone links the proximal and distal rows to-gether that leads this bone as the most commonly fractured carpal bone [1].To fix a scaphoid fracture, the required clinical courses consist of fracturediagnosis (Sec. 1.1.1), classification based on the diagnosis (Sec. 1.1.2) andcorresponding medical treatment (Sec. 1.1.3).1.1.1 DiagnosisScaphoid fractures are usually diagnosed by taking X-rays of the wrist(Fig. 1.2). However, the fractures may not be visible in initial (taken infirst week) X-rays. Considering this effect, if a patient has significant ten-derness directly over the scaphoid bone, he/she will be suspected of hav-ing a scaphoid fracture. In such questionable cases, computed tomography(CT)/Magnetic resonance imaging (MRI) scans are also used to help thediagnosis process.11.1. Clinical backgroundScaphoidRadiusTrapeziumTrapezoidUlnaPisiformTriquetrumCapitateHamate      1stmetacarpal      2ndmetacarpal      3rdmetacarpal      4thmetacarpal      5thmetacarpalLunateFigure 1.1: A schematic diagram showing the wrist bones of an examplewrist.1.1.2 TypesScaphoid fractures can be broadly categorized into two groups: displacedand non-displaced fractures. If any fracture develops a displacement be-tween the fractured parts of the scaphoid bone greater than 1 mm, then thefracture is defined as displaced, otherwise as non-displaced [68]. In addition,each type of fracture can be grouped into three sub-categories based on thelocation of the fracture in the scaphoid bone: distal pole fracture (the partclose to the upper carpal row), waist fracture (middle part) and proximalpole fracture (the part close to the radius bone). Examples of these threetypes of fractures are presented in Fig. 1.3.1.1.3 TreatmentEarly treatment (whether it is displaced or non-displaced) of this fractureincludes immobilization in a plaster cast, which prevents a further degra-dation of the fracture. The typical healing time is 10-12 weeks, however, it21.1. Clinical backgroundFigure 1.2: A plain X-ray of a fractured wrist (from Wikimedia Commons).(a) distal (b) waist (c) proximalFigure 1.3: Different types of scaphoid fractures.can be longer especially for a fracture located at the proximal pole of thescaphoid bone [64]. The reason behind such complication is the retrogradenature of the blood supply in the scaphoid bone, i.e., the artery enters onlyvia the distal pole of scaphoid bone and then travels to the proximal pole.31.1. Clinical backgroundTherefore, a fracture at the proximal pole does not allow the proximal frag-ment to receive nutrients for the recovery, as a result, the fractured partmight even die. For a better outcome and a faster recovery of such compli-cated cases, open or percutaneous surgical procedure is usually suggested,where a surgical screw is inserted along the longest axis of the fracturedscaphoid bone.Open surgical approachThe open surgical approach is usually recommended for displaced fracturecases as it allows a direct view of the targeted anatomy and provides aconfirmation of the anatomic reduction (i.e., close the fracture) before thescrew placement. However, as an invasive procedure, it is subject to damagesto blood supply vessels, surrounding soft tissue, ligaments and joints [27].Percutaneous surgical approachThe minimally invasive percutaneous surgical approach is usually suggestedfor non-displaced scaphoid fracture cases at the proximal pole. A tiny in-cision (1-2 mm long) is made on the palmar/dorsal side of the wrist, andthe screw is guided along its desired trajectory [101] usually under 2D fluo-roscopic guidance. The main limitation of the fluoroscopy-based navigationincludes the 2D projection view of the fluoroscopic images that make itchallenging to visualize the complex 3D anatomy of the carpal bones duringsurgery. Though the amount of dose is negligible in this kind of peripheralprocedures, the staff members frequently working in the operating roommay be affected by the X-ray radiation as a long term effect. Further-more, the close proximity of other carpal bones in the wrist joint interfereswith the acquisition of clear fluoroscopic images of the scaphoid. To over-come these difficulties, while reducing the number of fluoroscopic images,several computer-assisted fluoroscopic navigation systems have been devel-oped [57, 95]. These systems reduced the number of drilling attempts andX-ray exposure in the procedure, but did not improve the accuracy of screwplacement, likely because the complex 3D geometry of the scaphoid wasdifficult to interpret in 2D radiographic projections [87].To address the drawbacks of fluoroscopy, other guidance, e.g., cone beamcomputed tomography (CBCT) [87] or ultrasound (US) [7] has been sug-gested in the literature. A 3D guidance solution based on CBCT was pre-sented in [87], where volume rendering was used for the purpose of planningand navigation of the screw placement. However, due to the poor CBCT41.2. US guidance in computer assisted interventionscontrast in the narrow inter-bone spaces of the wrist, the visualization ofscaphoid surface in the rendered volume may be difficult. Beek et al. [7] pro-posed an US-based guidance approach for percutaneous scaphoid fracturefixation that uses preoperative CT data of the patient to derive a patientspecific model. Initially, a manual segmentation of the scaphoid surface inCT was performed, which was later registered to the semi-automatically seg-mented US bone surfaces. The main advantages of intraoperative US overfluoroscopy include the real-time 3D imaging, easy accessibility, low costand radiation free imaging. However, due to the smooth, featureless surfaceof the scaphoid bone, the registration of the segmented scaphoid surface tothe partial visible bone surfaces in US was challenging.1.2 US guidance in computer assistedinterventionsApart from the scaphoid fracture fixation, a significant number of previ-ous studies investigated the use of intraoperative US in different computerassisted applications (e.g., femur, pelvis, lumbar spine). Though there arecertain advantages in using intraoperative US, this data are often corruptedby significant amount of speckle noise and artifacts. In addition, US pro-vides only a partial view of the bony anatomy due to mismatches of accousticimpedances between the soft-tissue and bone surfaces (an example is shownin Fig. 1.4). Therefore, it has been suggested to fuse US data with anatom-ical information from CT (Sec. 1.2.1) or statistical model (SM) (Sec. 1.2.2)to obtain 3D contextual information of the bony anatomies during surgicalintervention.1.2.1 CT+US fusionTable 1.1 presents a summary of the CT+US-based methods suggested fordifferent anatomical applications [83]. These works can be broadly cate-gorized into two groups: surface-based (top part of Table 1.1) and volume-based (bottom part of Table 1.1) techniques. Among surface-based methods,the most common fusion approach includes a preoperative manual/semi-automatic CT segmentation, followed by the intraoperative registration ofthe segmented CT surface to the segmented bone surface in US. The bonesurfaces in US could be segmented either manually (e.g., Tonetti et al. [90])or automatically (e.g., Ionescuet al. [45]) or using the prior CT segmentationas a guidance (Amin et al. [2]). For the volume-based CT+US fusion, the51.2. US guidance in computer assisted interventionsTrue bone surfaceUS bone responseFigure 1.4: An example demonstration of the scaphoid bone response in US.key advantage is to avoid the CT or US segmentation (Brendel et al. [13])or both (Penney et al. [72]). Apart from these possibilities, there were alsodifferent registration examples reported in the literature, e.g., by simulatingUS images from CT and subsequently matching it to the US data (Wein etal. [97]). However, the volume-based techniques have more computationalcomplexities compared to the surface-based approaches.For both surface- and volume-based approaches, registration between theCT and US was achieved by optimizing the pose (rigid-body) parameters.For registration of multiple bones, either the ensemble of bone structures wasconsidered as a single rigid object (Brendel et al. [13]), or each bone wasindependently registered (Yan et al. [102]), or a biomechanical model wasused to regularize a multi-object registration (Gill et al. [30]). As the posesof the wrist bones change significantly while moving a wrist from one toanother position, expectedly a single rigid object registration does not pro-vide a good registration accuracy. In contrast, an independent registrationallows each of the wrist bones to move differently during the registration.However, this kind of registration does not ensure non-overlap between theneighboring bones, especially for the wrist anatomy with narrow inter-bonespaces between the bones. A bio-mechanical model instead addresses theaforementioned problems by constraining the transformation parameters ofthe bony anatomies during the registration. The main challenge of develop-ing a biomechanical model is the availability of all of the connecting complex61.3. Objectiveligament information.Although CT images can provide improved anatomical representation toUS, a significant challenge with the CT+US fusion is that the pose registra-tion of the joint bones does not always guarantee acceptable pose estimation,e.g., the registered bones could overlap with each other. In addition, thereis a need for extensive segmentation of the CT images in the surface-basedtechniques.1.2.2 SM+US fusionThe key reason behind introduction of the SM in US-CT fusion is to avoidthe extensive manual CT segmentation; in addition, a statistically accept-able shape and pose optimization is expected from the registration. Table 1.2presents a summary of the SM+US-based methods (top part: surface-basedand bottom part: volume-based). The SMs are known as powerful toolsin shape and pose analysis. The models are developed in order to capturethe main modes of features variations across population, and they can beused to register to unseen observation. Most of the developed models in theliterature include only the shape statistics; a few are based on multi-objectstatistical shape+pose analysis (e.g., Rasoulian et al. [81]). The multi-objectshape+pose model [12] was introduced to include the pose statistics in theintraoperative registration, as a results, an improved pose optimization wasachieved [81]. Similar to the CT-US fusions, the US images could be seg-mented manually or automatically. Either the SM could be registered di-rectly to US (Barratt et al. [5]) or simulated US images generated from thestatistical atlas could be registered to US (Khallaghi et al. [50]). A short-coming of the SM+US-based approach is the lack of quantitative measuresthat describe how faithful the statistical anatomical model is to the trueanatomy of the patient.1.3 ObjectiveIn an attempt to address the individual limitations of the CT+US andSM+US fusions, I suggest a registration framework, where preoperative CTimages are registered to intraoperative US via a statistical anatomical modelwithout the need for segmentation, to provide guidance of the screw inser-tion procedure for non-displaced percutaenous scaphoid fracture fixation.Fig. 1.5 shows an overview of my proposed SM+CT+US approach. In thissequential registration framework, a statistical anatomical model is initially71.3. ObjectiveTable 1.1: Summary of CT+US-based methods. The top and bottom partsof the table present the intensity- and volume-based approaches, respec-tively.Reference Year Anatomy ExperimentsIonescu et al. [45] 1999 Pedicle vertebra, 1×synthetic bone,ilio-sacral joint 1×cadaveric specimenTonetti et al. [90] 2001 Ilio-sacral joint 4×patientsAmin et al. [2] 2003 Pelvis 1×synthetic bone,1×patientBarratt et al. [5] 2008 Femur, pelvis 3×cadaveric specimenBeek et al. [7] 2008 Scaphoid 57×synthetic bonesMoore et al. [66] 2009 Lumbar spine 1×synthetic bone,1×cadaveric specimenBrounstein 2011 Pelvis 1×synthetic bone,et al. [14] 1× patientRasoulian 2012 Lumbar spine 5×synthetic bones,et al. [77] 1×cadaveric specimenGoncalves 2015 Femur 1×synthetic boneet al. [34]Wein et al. [98] 2015 Femur, tibia 1×cadaveric specimenand fibulaBrendel et al. [13] 2002 Lumbar spine 1×cadaveric specimenPenney et al. [72] 2006 Femur, pelvis 3×cadaveric specimenWein et al. [99] 2007 Abdomen 25×patientsWein et al. [97] 2008 Abdomen 10×patientsWinter et al. [100] 2008 Lumbar spine 5×patients,12× vertebraeGill et al. [30, 31] 2009, Lumbar spine 1×,6×synthetic bones,2012 1×cadaveric specimenYan et al. [102] 2011 Lumbar spine 1×phantom,1×cadaveric specimenHacihaliloglu 2013 Pelvis 1×synthetic bones,et al. [42] 10×patientsregistered to preoperative CT to obtain the shape and scale, and subse-quently to intraoperative US to derive the pose of multiple bones. Since CTis usually acquired in the diagnosis step to make decision whether or not to81.3. ObjectiveTable 1.2: Summary of SM+US-based methods. The top and bottom partsof the table present the intensity- and volume-based approaches, respec-tively.Reference Year Anatomy ExperimentsStindel et al. [88] 2002 Distal femur 1× dry bone,1× cadaveric specimenChan et al. [16] 2004 Femur, pelvis 3× cadaveric specimenKilian et al. [51] 2008 Distal femur 1× dry bone,1× cadaveric specimenBarratt et al. [6] 2006 Femur, pelvis 3× cadaveric specimenForoughi et al. [26] 2008 Pelvis 1× dry bone,2× cadaveric specimenSchumann et al. [84] 2012 Pelvis 1× synthetic bone,1× cadaveric specimenSchumann et al. [85] 2012 Pelvis 2× synthetic bones,2× dry bonesRasoulian et al. [81] 2015 Lumbar spine 13× patientsKhallaghi et al. [50] 2010 Lumbar spine 3× synthetic bonesGhanavati et al. [29] 2010 Pelvis 2× cadaveric specimenGhanavati et al. [28] 2011 Pelvis 5× synthetic bones,1× dry bonego for operative treatment for the scaphoid fractures [82], no additional CTacquistion is required for my approach. For the US registration, I developan automatic bone enhancement method to enable the SM registration toUS for determining the optimal scaphoid screw axis.In addition to the aforementioned two-step sequential registration ap-proach, I also propose a joint model registration (Fig. 1.6) to CT and USas an alternate way to derive the patient-specific 3D models of the wristbones. Instead of estimation of the shape and scale of the wrist bones onlyfrom preoperative CT, within the proposed joint registration framework,these estimations are based on both CT and US images to acheive a betterregistration accuracy.Finally, it is important to note that the proposed approaches (both se-quential and joint) do not interfere with established clinical work-flow andsetup, while they propose significant reduction in ionizing radiation andimproved patient safety.91.4. ContributionsCT registration(shape + scale+ pose)CT volumeStatistical modelUS volumeBone enhancementBone enhanced volumeUS registration(pose)Preoperative registered modelIntraoperative registered modelScaphoid screw axis estimationIntraoperativeWrist at neutralScrew overlapped on scaphoidCT acquisitionUS acquisitionWrist at extensionPreoperativeFigure 1.5: The proposed workflow for scaphoid screw axis estimation. Pre-operatively, the statistical wrist model is registered to the preoperative CTat a neutral wrist position for estimation of the size and shape of the wristbones. Intra-operatively, a 3D US volume is acquired at extension wristposition and the wrist model is registered to the bone surfaces extractedfrom the US volume. In this step, only pose optimization is performed andshape and scale parameters are applied as calculated from the CT volume.Finally, the scaphoid screw position is obtained by maximizing the screwlength within a safe zone inside the registered scaphoid bone.1.4 ContributionsThe thesis is an attempt to enable US guidance in a percutaneous non-displaced scaphoid fracture fixation for estimation of the optimal scaphoidscrew axis. It is interesting to note that the proposed work-flow could also beused in other US-based orthopaedic applications corresponding to differentanatomies of human body (e.g. lumbar spine). In the course of achievingthis objective, the following contributions were made:• Development of a multi-object statistical shape+scale+pose model101.4. ContributionsCT volumeStatistical modelUS volumeBone enhancementBone enhanced volumeJoint CT & US registrationRegistered model to USScaphoid screw axis estimationIntraoperativeWrist at neutralScrew overlapped on scaphoidCT acquisitionUS acquisitionWrist at extensionPreoperativeRegistered model to CTFigure 1.6: The proposed workflow for joint registration approach. Insteadof sequential registration, the statistical wrist model is jointly registered toUS and CT. Note that the shape and scale of the wrist bones are determinedfrom the joint optimization. Poses of the bones are estimated indepedentlyto capture the pose differences across the US and CT acquisitions.that captures the main modes of shape, scale and pose variations acrossa training population at different wrist positions.• Registration of the developed statistical model to preoperative CTfor automatic segmentation of the wrist bones in CT images. Theregistration allows us to estimate the shape and size of the wrist bonesthat will be used in the intraoperative US registration.• Development of an automatic bone enhancement technique in US thatenhances the blurry, disconnected, noisy bone responses in US com-pared to those responses from other anatomical interfaces. The en-hancement allows a better interpretation of the bone responses in theintraoperative US that in turns enables registration of the model toUS.• Registration of the statistical model to intraoperative US to obtain111.5. Thesis outlinethe poses of the wrist bones during surgical navigation for estimationof the scaphoid screw axis.• Joint registration of the statistical model to preoperative CT and in-traoperative US, which allows a simultaneous fusion of the model toCT and US to derive the patient-specific 3D models of the wrist bones,even where the bone poses substantially vary between the CT and US.1.5 Thesis outlineThe rest of this thesis is subdivided into seven chapters as outlined below:CHAPTER 2: NOVEL STATISTICAL SHAPE+SCALE+POSE MODELFOR WRISTA statistical model captures different feature variations across population sothat the model is capable to generate an unseen observation based on thestatistics. In this chapter, I introduce a statistical wrist shape+scale+posemodel that covers the main modes of shape, scale and pose variations acrossa training population at different wrist positions. To establish the correspon-dences across the training shapes, the wrist bone surfaces are jointly alignedusing a group-wise registration framework based on a Gaussian MixtureModel. Next, the development of shape model starts with the generation ofthe mean shape, and subsequent calculation of the shape residual at tangentspace. Principal component analysis is then used to determine the majormodes of shape variations. The scale and pose models development startswith the computation of the similarity transformation between each wristbone of the mean shape to the corresponding bone of each training sam-ple, followed by mapping them at tangent space. For the pose statistics,the variations in poses not only across the population but also across differ-ent wrist positions are incorporated in two pose models. An intra-subjectpose model is developed by utilizing the similarity transforms at all wristpositions across the population. Further, an inter-subject pose model is con-structed to model the pose variations across different wrist positions.CHAPTER 3: AUTOMATIC SEGMENTATION OF WRIST BONES US-ING STATISTICAL MODEL REGISTRATION TO CTSegmentation of the wrist bones in CT images has been frequently used indifferent clinical applications including arthritis evaluation, bone age assess-ment and image-guided interventions. The major challenges for development121.5. Thesis outlineof an automatic segmentation technique include non-uniformity and spongytextures of the wrist bone tissue as well as narrow inter-bone spaces. In thischapter, for automatic segmentation of the wrist bones in CT images, thedeveloped model (in the previous chapter) is registered to the edge pointcloud extracted from the CT volume through an expectation maximizationbased probabilistic approach. Residual registration errors are corrected atthe end by application of a non-rigid registration technique.CHAPTER 4: BONE ENHANCEMENT IN US USING LOCAL SPEC-TRUM VARIATIONThe main challenge of the US based computer assisted interventions is thedifficulty in interpreting the US images due to the low-contrast and noisy na-ture of the data. In this chapter, I propose bone enhancement methods thatexploit local spectrum features of the ultrasound images. These featuresare utilized to design a set of quadrature band-pass filters and subsequentlyestimate the local phase symmetry, where high symmetry is expected at thebone locations. I also incorporate the shadow information below the bonesurfaces to further enhance the bone responses.CHAPTER 5: SCAPHOID SCREW AXIS ESTIMATION BASED ON USGUIDANCEIn this chapter, the extracted bone surfaces in US are used to register thedeveloped model to US volumes, allowing to derive a patient specific 3Dmodel for guiding scaphoid fracture fixation. The registered model is thenused to determine clinically important intervention parameters, includingthe screw length and the trajectory of screw insertion in the scaphoid bone.CHAPTER 6: JOINT ATLAS-BASED REGISTRATION OF MULTIPLEBONES IN CT AND USIn this chapter, an alternate way is presented to derive the patient-specific3D models of the wrist bones by jointly registering the model to CT and US.The registration is facilitated by the joint optimization of the atlas shapeand scale coefficients using information from both imaging modalities.CHAPTER 7: CONCLUSIONA short summary about the thesis is presented in this chapter, followed bya discussion of the proposed approach and possible future extensions.13Chapter 2Novel statisticalshape+scale+pose model forwrist2.1 Introduction and backgroundStatistical models (SMs) are widely recognized as powerful tools in severalclinical applications such as image segmentation [21], kinematic analysis [18]or prosthesis development. The purpose of a SM is to capture differentfeature (e.g., shape, pose) variations across a set of training examples. Asone of the objectives in this thesis is to develop a SM of a wrist joint,therefore, it is desirable to capture different feature variations across a setof wrist training examples at wide range of wrist positions.A study on the literature of SM development indicates that majorityof works were based on capturing the shape differences of a single-objectacross population. The first attempt [21] to develop a single-object SM re-quired manual landmark identification for establishing the point correspon-dences across training shapes. However, manual landmark identification is atime consuming and subjective task. To address this problem, a number ofworks [43, 48] had been published based on surface matching algorithms likeIterative Closest Point (ICP) [9] or Softassign Procrustes [76]. They chosean arbitrary shape from a set of training shapes as a reference and foundThis chapter is adapted from the following three papers: 1) Anas E. M. A., RasoulianA., John P. St., Pichora D., Rohling R. N., and Abolmaesumi P., A Statistical Shape+PoseModel for Segmentation of Wrist CT Images, In SPIE Medical Imaging, 90340T, 2014,and 2) Anas E. M. A., Rasoulian A., Seitel A., Darras K., Wilson D., John P. St., PichoraD., Mousavi P., Rohling R. N., and Abolmaesumi P., Automatic Segmentation of WristBones in CT Using a Statistical Wrist Shape+Pose Model, IEEE Transactions on MedicalImaging, 13(11): 2025-2034, 2016, and 3) Anas, E. M. A., Seitel A., Rasoulian A., John P.St., Ungi T., Lasso A., Darras K., Wilson D., Lessoway V., Fichtinger G., Zec M., PichoraD., Mousavi P., Rohling R. N., and Abolmaesumi P., Registration of a statistical model tointraoperative ultrasound for scaphoid screw fixation, International Journal of ComputerAssisted Radiology and Surgery, 1-9, 2016.142.1. Introduction and backgroundthe optimal similarity transformation from one shape to the other. However,bias was introduced in the generation of SM due to the choice of referenceshape. To solve this problem, group-wise registration techniques were pro-posed [4, 10, 78] that assumed both the mean shape and its transformationto the training shapes are unknown.Though the above mentioned works demonstrated quite good perfor-mance in different applications, they could not be directly applicable to mywrist fracture fixation problem as only shape variations are not sufficientto capture the feature variations of the wrist joint across population. Forsuch a complex anatomical structure, multi-object joint shape+pose mod-els [11, 12, 73] are usually used. The purpose of a multi-object shape+posemodel is to capture the shape differences of the joint bones, as well as tomodel the relative position variations of the wrist bones across population.Bossa et al. [11, 12] presented techniques for multi-object shape+pose mod-els by concatenating the shape and pose features in a long vector, and thenapplying standard multivariate statistical tools, such as Principal Compo-nent Analysis (PCA) or Principal Geodesic Analysis (PGA), to identifycommon modes of variation. However, such statistical tools could only beused to determine the variations at a specific anatomical position in a pop-ulation, and are not directly applicable for determining the statistics over alarge range of wrist motions.To model the complex variations of the wrist joint for relatively largemotions, a few techniques have been proposed in the literature. Recently,Van de Giessen et al. [93] estimated the coordinate differences of pairs ofpoints on adjacent bone surfaces of a healthy wrist, and developed a LocalStatistical Motion Model (LSMM) based on the distribution of these differ-ences as a function of the wrist position. However, this method only modeledthe inter-bone distances and it did not include the shape or pose variations.Chen et al. [18, 19] proposed a joint statistical model to include the shapeand pose variations obtained from a number of 3D CT data sets of a groupof subjects at various wrist poses, where the required correspondences weredetermined using a pair-wise registration approach. However, their methodincluded pose variation only due to the large wrist motion, i.e., it discardedthe relative pose differences due to the population variations.In this chapter, I present a statistical shape+scale+pose model of thewrist, capturing the shape, scale and pose variations of the joint across apopulation and over a wide range of wrist motions. I use a wrist database [65]that includes segmented 3D surface models of carpal bones from 30 subjects(both arms), including 15 females (ages 21-28) and 15 males (ages 22-34)at different targeted wrist positions. Initially, I use a group-wise Gaus-152.2. Methodssian Mixture Model (GMM)-based registration technique to establish thecorrespondences among the shape points of the wrist bones across differ-ent example sets. Next step is to generate the mean shape of each carpalbone, followed by the development of the shape model [80]. After that Idetermine the similarity transformations for each wrist bone at all wrist po-sitions with respect to the mean shape of that carpal bone. These similaritytransformations can be broken down into two transformations: A transfor-mation containing only the scaling factor and a transformation related tothe rigid-body movement. One of the contributions in my developed modelis to separate the scale statistics from the pose statistics unlike the previoustechniques [11, 12, 80], which allows us to use the shape and size of the wristbones estimated from the CT data in the US registration procedure. Thescaling factors are utilized to develop a multi-object scale model, and thesimilarity transformations related to the rigid-body movement are used todevelop inter-subject and intra-subject pose model, another contribution inmy model. I refer to the variation of the wrist poses (for all wrist bones) in asubject population at a particular wrist posture as the “inter-subject pose”variation and the variation of the wrist poses within an individual subjectas a result of changing the wrist position as “intra-subject pose” variation.2.2 MethodsThe statistical model in this section is constructed such that it covers themain modes of shape, scale and pose variations across a training populationat different wrist positions. The publicly available wrist database developedby Moore et al. [65] is used as a training data set. A brief description of thewrist database is presented in Sec. 2.2.1. The information in the databasehas been used during the development of different parts of the statisticalmodel (shape: Sec. 2.2.2, scale: Sec. 2.2.3, and pose: Sec. 2.2.4). Sec. 2.2.5presents the steps of generating a wrist instance from a given set of shape,scale and pose modes. Finally, I describe (Sec. 2.2.6) how I register thedeveloped model to a target 3D point cloud.2.2.1 Digital wrist databaseTo develop the statistical wrist model I have used the wrist database pub-lished by Moore et al. [65] that includes segmented 3D surface models ofwrist bones at neutral wrist position from 30 healthy persons (no history ofwrist trauma or signs of wrist pathology or obvious degenerative changes)with both arms, including 15 females (ages 21-28) and 15 males (ages 22-34),162.2. Methodswhere it was found that the bones of the female volunteers were smaller onaverage than those of the males (by an average of 33.3%). The surface mod-els were generated from CT segmentation, where a GE helical CT scanner(GE Medical, Milwaukee, WI) was used to scan both wrists of each vol-unteer simultaneously. 60-80 slices were acquired from each individual with1.0 mm slice thickness. The typical field-of-view includes 1-2 cm of the distalradius and ulna, the entire carpus, and 1-2 cm of the proximal portion of themetacarpals. The CT of each wrist at neutral position was segmented usinga combination of commercially available and custom-written software. Allthe left wrists were converted to right wrists by inverting the X componentof the bone surface models.The data set contains CT volumes at different wrist positions that areacquired using i) the incremental orthogonal protocol, and ii) the combinedmotion protocol. The incremental orthogonal protocol (5 males, 5 females)included scans at eight targeted positions: neutral position, 30◦ and 60◦ offlexion, 30◦ and 60◦ of extension, 20◦ and 40◦ of ulnar deviation, and 20◦of radial deviation. The combined motion protocol (10 males, 10 females)included scans at nine targeted wrist positions: neutral position, 40◦ offlexion, 40◦ of extension, 10◦ of radial deviation, 30◦ of ulnar deviation,and combined motions of 40◦ flexion and 30◦ ulnar deviation, 40◦ extensionand 30◦ ulnar deviation, 40◦ extension and 10◦ radial deviation, and 40◦flexion and 10◦ radial deviation. The positions of the wrist bones at eachnon-neutral wrist position are calculated with respect to those at the neutralwrist position using rigid body kinematics provided in the wrist database. Inaddition to these shape and pose information, the database provided originand three principal axes of radius-based coordinate system (RCS) for eachwrist example set.2.2.2 Development of statistical shape modelAssume my training set consists of Ns example shapes of an ensemble ofL wrist bones, and the shape of each bone is represented by a set of Msurface points. In this work, I have Ns = 60 (30 wrists, both arms), L =9 (eight carpal bones+radius) and M = 1502. The radius bone is notcompletely visible in the CT volumes and its covered portion also variesacross the population. I account for this variation during model developmentby determining the portion of radius that can be consistently extracted fromthe training samples. To accomplish this, I calculate the ratio between thevolume of the radius segmentation and the average volume of all other carpalbones for each case. The radius segmentations are then cut along the plane172.2. Methodsperpendicular to the radius long axis of the RCS coordinate system suchthat the resulting new ratio for each training sample is close to the smallestratio found in the data set. The resulting segmentation is then used inthe consequent model development. I was not able to include the ulna andmetacarpal bones because a uniquely defined coordinate system was notavailable in the wrist database.Fig. 2.1 shows a schematic diagram of different processing steps of myshape model development approach. Note that all example shapes are ob-tained at neutral wrist position and it is expected that shapes will not changedue to changes in the wrist position. Therefore, generating the shape statis-tics only at the neutral wrist position is sufficient to describe the shapestatistics at any wrist position. To generate the shape model, I have used ourpreviously published multi-object shape model development technique [79].In short, a group-wise Gaussian mixture model (GMM)-based registrationtechnique [78] is used to establish the correspondences across the trainingshapes of each wrist bone. After establishment of the correspondences, allthe wrist instances are aligned using the RCS coordinates provided in thedatabase, where the origins of the RCS coordinates are expected to be con-sistently located at same anatomical location across the training examples.After aligning all the example shapes of each wrist bone, the mean shapefor l-th wrist bone is obtained across the population and it is representedas Πsl (M × 3 matrix). Here, the right superscript ‘s’ stands for shapestatistics. Then, I use Principal Geodesic Analysis (PGA) [89] to exponen-tially map 3D coordinates of the correspondences to the tangent space ofthe mean shape and perform Principal Component Analysis (PCA) in thetangent space. PGA is an extension of PCA [89] which is typically usedin non-linear analysis for pose statistics. In this work, I chose to use PGAover PCA to maintain consistency with the computation of pose statistics.There is no disadvantage in using PGA over PCA for shape analysis in thiscontext. Next, the joint shape features are formed and the result of PGA areCs principal components in the tangent space, where the c-th mode princi-pal eigenvector is vsc (= (vsc,1T,vsc,2T, · · · ,vsc,LT)T) and c-th mode eigenvalueis λsc, and 1 ≤ c ≤ Cs. Shape Principal Geodesics (PGs) are obtained forthe l-th wrist bone by exponential mapping, expΠsl(Vsc,l), where Vsc,l is aM×3 matrix generated by rearranging the 3M×1 vector vsc,l. A new shapeinstance Φsl (M × 3 matrix) of the l-th wrist bone can then be described asΦsl (Θs) = expΠsl(Cˆs∑c=1θsc Vsc,l), (2.1)182.2. Methodswhere θsc represents the weight of the c-th shape mode, Cˆs (≤ Cs) is thetotal number of shape coefficients used for instance generation, and Θs =(θs1, θs2, · · · , θsCˆs)T are the shape coefficients.1…1…1…..1…1…1…..Correspondence establishmentAligning1…2Mean shapegeneration12..Wrist examplesAt neutralWrist examplesAt neutralWrist examples aligned to one common coordinateMean shapesSimilarity transformationneutralFigure 2.1: A schematic diagram showing the example wrist alignment,mean shape generation and similarity transformation calculation. Corre-spondences among the wrist examples are established using a group-wiseregistration-based framework. Subsequently, the wrist examples are alignedin a common coordinate and mean shape is generated across all exampleshapes for each wrist bone. Similarity transformation Tn,l is computed forl-th bone of n-th example from the mean shape to the l-th bone of the n-thexample.2.2.3 Statistical scale modelThe scale model development starts with the computation of the similaritytransform Tn,l (4×4 matrix) between the l-th wrist bone of the mean shape(Πsl ) to the corresponding bone of the n-th training sample at a neutralwrist position using the generalized procrustes analysis. The 4 × 4 matrixTn,l can be represented as,Tn,l =[kn,lRn,l tn,l0 1]. (2.2)192.2. Methods(a) Inter-pose model developmentx1y1z1rxryrzh(7)(7L)(7LK)(7LK)(7L)................PCA across nSimilaritytransformationMean across kMapping attangent spaceFeaturegenerationMean acrossnand k Mapping attangent spaceFeaturegenerationx1y1z1rxryrzh(7)(7L)(7L)........PCA across kSimilaritytransformation1. Inter-pose feature generation2. Inter-pose statistics1. Intra-pose feature generation2. Intra-pose statistics(b) Intra-pose model development: similarity transformation of l-th carpal bone of n-th sample at k-th wrist position : mean transformation of l-th carpal bone: residual transformation: pose feature: c-th joint inter-pose eigenvector of the l-th carpal bone at the k-th wrist position: c-th mode inter-pose coefficient: a pose instance of the wrist model of the l-th carpal bone at k-th wrist position: total number of inter-pose modes: mean transformation of l-th carpal bone: residual transformation: pose feature: c-th joint intra-pose eigenvector of the l-th carpal bone: c-th mode intra-pose coefficient: a pose instance of the wrist model of the l-th carpal bone: total number of intra-pose modesFigure 2.2: A schematic diagram showing the process of developing thetwo inter- and intra-subject pose models. I concatenate the pose featuresfrom all wrist bones at different wrist positions to form a joint pose fea-ture. Initially, PCA extracts the inter-subject pose variations. After thatan instance created from the inter-subject pose model is used to develop theintra-subject pose model that represents the pose variations across differentwrist positions.Here, kn,l, Rn,l and tn,l represent the scaling factor (a scalar value), the rigid-body rotation matrix (a 3×3 matrix) and the rigid-body translational vector(a 3 × 1 vector), respectively, for the l-th wrist bone of the n-th trainingsample. 0 indicates a zero-valued 1×3 row vector. Tn,l can be broken downinto two transformations: the scaling and the rigid-body transformations,i.e., Tn,l = T2n,l ×T1n,l, whereT1n,l =kn,l 0 0 00 kn,l 0 00 0 kn,l 00 0 0 1 , (2.3)T2n,l =[Rn,l tn,l0 1]. (2.4)202.2. MethodsIn the scale model development, I model the variations of kn,l across n. Themean pil of the scaling factors for the l-th wrist bone is determined, andsubsequently the residual scale factor at the tangent space for the l-th wristbone of the n-th training sample is computed as hn,l = log(kn,lpil).PCA is then applied to develop the scale model, consisting of a set ofeigenvectors vkc = (vkc,1, vkc,2, · · · , vkc,L)T and a set of eigenvalues λkc , wherethe total number of eigenvectors is Ck. The right superscript ‘k’ indicatesthe parameters related to the scale model. A scale instance for the l-th wristbone can be generated asφkl (Θk) = pilCˆk∏c=1exp(θkc vkc,l), (2.5)where θkc is the c-th scale mode coefficient, Cˆk (≤ Ck) is the total number ofscale coefficients used for instance generation, and Θk = (θk1 , θk2 , · · · , θkCˆk)Tare the scale coefficients.2.2.4 Development of statistical pose modelIn the following, I present the steps of developing the inter-subject pose andintra-subject pose models. An illustration of generating the pose statisticsis shown in Fig. 2.2.Similarity transformations and their meanThe pose model development starts with the calculated similarity transfor-mation T2n,l (in previous section) between the l-th wrist bone of the meanshape (Πsl ) to the corresponding bone of the n−th training sample at neu-tral wrist position. To include similarity transformations from other wristpositions, I have used the neutral to non-neutral transformations providedin the wrist database, i.e., Tn,l,k = (k1Tn,l)(T2n,l), where 2 ≤ k ≤ K, and Kis the total number of wrist positions. Here, k1Tn,l indicates the rigid-bodytransformation matrix that maps surface points of the l-th wrist bone ofthe n-th subject from neutral to the k-th non-neutral position. The meantransformation Πpl is computed from Tn,l,k across the population n and thewrist position k for each wrist bone. The right hand superscript ‘p’ denotesthe pose related variables to distinguish from Πsl .212.2. MethodsPose feature generationUsing the mean transformation, the residuals in the tangent space, Rpn,l,k =log(Πpl−1Tn,l,k) are calculated. The 4× 4 matrix Rpn,l,k can be representedas,Rpn,l,k =0 −rz ry xrz 0 −rx y−ry rx 0 z0 0 0 0 (2.6)For each Rpn,l,k, a pose feature vector upn,l,k is defined [11] and this vector isrepresented by six free parameters in the tangent space, i.e.,(x, y, z, rx, ry, rz)T . Here, x, y, z are related to the translational parameters;and rx, ry, rz represent the Rodrigues parameters [18]. Therefore, one set ofjoint pose feature can be expressed as upn,k = (upn,1,kT,upn,2,kT, · · · ,upn,L,kT)T, where 1 ≤ n ≤ Np and 1 ≤ k ≤ K. Np and K represent the total numberof pose samples and total wrist positions, respectively.Inter-subject pose statisticsFor pose model development, I have used data from the combined motionprotocol that includes nine (K = 9) targeted wrist positions of 40 (Np = 40)samples. I concatenate the pose features from all different wrist positions(Fig. 2.2(b)), i.e., upn = (upn,1T,upn,2T, · · · ,upn,KT)T. The newly generatedlong concatenated vector includes pose information of all wrist bones at allnine wrist positions of a particular training sample. Subsequently, PCAis applied to upn across n to extract Cp principal components representingthe inter-subject pose variations, vpc = (vpc,1T,vpc,2T, · · · ,vpc,KT)T, where,vpc,k = (vpc,1,kT,vpc,2,kT, · · · ,vpc,L,kT)T. The corresponding eigenvalues fromthe PCA analysis are λpc . Inter-subject pose PGs are generated by theexponential mapping, Πpl exp (Vpc,l,k). Given a 6 × 1 vector vpc,l,k, a 4 × 4matrix Vpc,l,k is generated using eq. (2.6). The similarity transform Φpl,krepresenting a new pose instance for the l-th bone at the k-th wrist positioncan be expressed asΦpl,k(Θp) = ΠplCˆp∏c=1exp(θpc Vpc,l,k), (2.7)where θpc denotes the weight of the c-th inter-subject pose mode, Cˆp isthe total number of inter-subject pose modes used for instance generation,222.2. Methodsand the inter-subject pose coefficient vector Θp = (θp1, θp2, · · · , θpCˆp)T. Bychanging the value of k, I can generate the pose of any wrist bone at any ofthe K wrist positions.Intra-subject pose statisticsTo model the positions of wrist bones at different wrist positions, I furtheranalyze the obtained pose instances. Analyses are carried out on Φpl,k, whereI compute the residual transformation and the subsequent pose features asdetailed in the previous sections. In short, Πp′l is the mean computed fromΦpl,k across k, and the residuals are expressed as Rp′l,k = log(Πp′l−1Φpl,k) inthe tangent space. I use the right hand subscript ‘p′’ to indicate the intra-subject pose coefficient related variables only. Subsequently, I compute thepose feature up′l,k corresponding to Rp′l,k, and generate the joint pose feature asup′k = (up′1,kT,up′2,kT, · · · ,up′L,kT)T. PCA is applied to up′k across k to computeintra-subject pose PGs across different wrist positions. The eigenvector isvp′c = (vp′c,1T,vp′c,2T, · · · ,vp′c,LT)T, and the corresponding eigenvalues are λp′c ,where 1 ≤ c ≤ Cp′ and Cp′ is the total number of intra-subject pose PGs.A pose instance for the l-th bone can be calculated asΦp′l (Θp,Θp′) = Πp′lCˆp′∏c=1exp(θp′c Vp′c,l). (2.8)Here, Vp′c,l is a 4×4 matrix generated from a 7×1 vector vp′c,l. θp′c represents c-th intra-subject pose coefficient and Cˆp′is the total number of intra-subjectpose modes used for instance generation. The intra-subject pose coefficientvector is Θp′= (θp′1 , θp′2 , · · · , θp′Cˆp′)T. Note that Φp′l is function of Θp′ as wellas Θp.2.2.5 Model representationGiven the shape (Θs), scale (Θk) and two sets of pose (Θp,Θp′) coefficients,a wrist instance for the l-th wrist bone can be generated asYl(Θs,Θk,Θp,Θp′) = (2.9)φkl (Θk)Φsl (Θs)× {Rp′l (Θp,Θp′)}T + IM × {tp′l (Θp,Θp′)}T ,where Rp′l and tp′l are the rotation matrix (3×3 matrix) and the translationalvector (3 × 1 vector), respectively, which are determined from the 4 × 4232.3. Evaluationsimilarity transformation matrix Φp′l . IM is a vector of ones of size of M×1.A complete wrist instance is then determined from L wrist bones usingl = 1, · · · , L.2.2.6 Model registrationThe model registration to a target point cloud is performed using a GMM-based technique [80] in which the target point cloud (number of target pointsN) is assumed to be an observation generated by a GMM, and the Expec-tation Maximization (EM) framework is applied to optimize the model pa-rameters. In the expectation step of the EM framework, the probability ofhow likely a point in the model generates another point from the target iscomputed.P (tml |zn) =exp(−12 ||zn−yml (Θs,Θk,Θp,Θp′)σ ||2){∑Ll=1∑Mm=1 exp(−12 ||zn−yml (Θs,Θk,Θp,Θp′ )σ ||2)+(2piσ2)1.5(w1− w)(MLN)},(2.10)where zn represents the n-th point of the target point cloud,yml (Θs,Θk,Θp,Θp′) is the m-th point of the developed model belonging tothe l-th wrist bone, L is the total number of wrist bones used, M is the totalnumber of points per bone in the model, N is the total number of pointsin zn and w is a factor to account for outliers. Θs, Θk, Θp and Θp′repre-sent the shape, scale, inter-subject pose and intra-subject pose coefficients,respectively. σ is the standard deviation of the Gaussian components whichis chosen to be large initially and its value is decreased in each iteration. Inthe maximization step, the following objective function is minimized:Q =L∑l=1M,N∑m,n=1P (tlm|zn)||zn − yml (Θs,Θp,Θp′)||2, (2.11)The cost function is minimized alternately with respect to the model co-efficients using the Quasi-Newton method that transform the model into asurface registration closest to the actual target point cloud.2.3 EvaluationFor evaluation of the developed statistical model, a compactness analysiswas performed for each (e.g. shape or scale) model to measure the variations242.4. Resultscaptured across the number of principal modes. Compactness is the abilityof the model to use a minimal set of parameters, and it is measured as afunction of mode c. The compactness at mode c is the ratio of sum of thevariances of first c modes to the sum of the variances of all modes.In addition, the capability of the model to register to any wrist positionof an unseen example was determined in two leave-one-out experiments:leave-one-subject-out (LOSO) and leave-one-pose-out (LOPO) experiments.In LOSO test, the statistical wrist model was developed using all but onesubject of the training set, such that the training data set consisted of a totalof 58 cases (29 left and 29 right wrists) and the test case contained both leftand right wrists. It was important to consider both sides of the wrist in thisanalysis because the corresponding left and right bone surfaces are rathersimilar and may bias the generalization result of the model. In contrast, thestatistical wrist model in LOPO experiment was developed using all but onewrist pose (including both arms) of the training set.In both experiments, the model was registered to each excluded subject(both arms) at each wrist position, and the mean surface distance error(mSDE) between the registered and the target point cloud was calculated. Icalculated the surface distance error (SDE) for each point of the registeredsurface by computing its Euclidean distance to the corresponding point ofthe reference surface.2.4 ResultsThe compactness analysis of the shape, scale and pose models is presented inFigs. 2.3(a-d). It indicates the differences among the number of the modesof the models to capture a specific percentage of variation. For example, a95% of the variations can be covered by 50 shape modes, 2 scale modes, 8inter- and 2 intra-subject pose modes.Figs. 2.4(a-d) demonstrate the variations of the wrist model due to thechanges of first principal mode of the shape, scale, inter-subject and intra-subject pose models, respectively. For each principal mode, the model in-stance is generated at two different values of the corresponding coefficient,i.e., −1.5√λc and 1.5√λc, where λc indicates c-th principal mode eigenvalueof the corresponding model.Tables 2.1 and 2.2 show the ability of the shape+scale+pose model toregister to different wrist positions using LOSO and LOPO experiments,respectively. For LOSO experiment, we obtain mSDEs in the range between0.29 mm (0.28 mm for LOPO) for the neutral wrist position to 0.39 mm252.5. Discussion10 20 30 40 50102030405060708090100Number of shape modesShape compactness1 2 3 4 5 6 7 8 9949596979899100Number of scale modesScale compactness5 10 15 20 25 30 35 40 45 50556065707580859095100Number of inter−subject pose modesInter−subject pose compactness1 2 3 4 5 6 7 865707580859095100Number of intra−subject pose modesIntra−subject pose compactness(a) (b)(c) (d)Figure 2.3: Compactness of the developed (a) shape, (b) scale, (c) inter-subject pose and (d) intra-subject pose models.(0.41 mm for LOPO) for the wrist positioned in 40◦ ulnar deviation withthe numbers of modes as chosen earlier.2.5 DiscussionThe compactness analysis in Fig. 2.4 indicates the scale model as the mostcompact one among the four models. It is expected as only one scale featureis used in the scale model development, where more number of shape andpose features have been used in the shape and pose model development,respectively.The intra-subject pose model represents the pose variations of the wristsbones due to the changes of the wrist positions, e.g., the first principal mode(Fig. 2.4(d)) in the intra-subject pose model represents pose variation alongextension-flexion directions. I have also investigated the second mode and262.5. Discussion(a) Shape (b) Scale (c) Inter-subject pose (d) Intra-subject poseFigure 2.4: Variations of the wrist model by changing the first mode of the(a) shape, (b) scale, (c) inter-subject pose and (d) intra-subject pose models.In each figure, the white and green surfaces represent two model instancesat extreme values of the corresponding mode.it represents the variations along the radial-ulnar deviations.A comparison between the inter-subject and intra-subject pose modelsindicates more movement of the wrist bones in the intra-subject pose modelcompared to that in the inter-subject pose model (Fig. 2.4(c) vs. Fig. 2.4(d)).The ability of the model to be registered to target point clouds at differ-ent wrist positions of different unseen examples have been demonstrated inLOSO and LOPO experiments on the model training data where a mSDElower than 0.41 mm could be achieved. Note that mSDEs are higher forsome of the cases in the experiments. This can be mainly attributed tothe initial misalignment as my current initialization is based on the a wristinstance at neutral position. However, if a priori information on the targetwrist position is provided, the initialization can be performed with a model272.5. DiscussionTable 2.1: Average and standard deviation of the mSDEs achieved by reg-istering the model to surfaces of the training data in the LOSO experiment.At each wrist position, the mSDE is computed for each training sample, andthe average mSDE is computed across all training samples at that particularwrist position. Note that the sample size is different at different wrist posi-tions. Flexion (F), Extension (E), ulnar deviation (UD) and radial deviation(RD).Wrist positions average mSDE (mm) sample sizeNeutral 0.29±0.10 6030◦ F 0.29±0.09 2040◦ F 0.30±0.11 4060◦ F 0.32±0.12 2030◦ E 0.31±0.12 2040◦ E 0.30±0.11 4060◦ E 0.35±0.14 2020◦ UD 0.28±0.13 2030◦ UD 0.35±0.13 4040◦ UD 0.39±0.14 2010◦ RD 0.33±0.13 4020◦ RD 0.33±0.12 2040◦ F and 30◦ UD 0.37±0.14 4040◦ E and 30◦ UD 0.33±0.14 4040◦ E and 10◦ RD 0.31±0.10 4040◦ F and 10◦ RD 0.34±0.12 40instance at the approximate target pose for better initial alignment.282.6. SummaryTable 2.2: Average and standard deviation of the mSDEs achieved by regis-tering the model to surfaces of the training data in the LOPO experiment.The sample sizes for this table are same as those in Table 2.1Wrist positions average mSDE (mm)Neutral 0.28±0.0930◦ F 0.30±0.1040◦ F 0.32±0.1160◦ F 0.32±0.1330◦ E 0.33±0.1140◦ E 0.30±0.1060◦ E 0.34±0.1320◦ UD 0.29±0.1330◦ UD 0.37±0.1440◦ UD 0.41±0.1510◦ RD 0.35±0.1320◦ RD 0.34±0.1240◦ F and 30◦ UD 0.40±0.1440◦ E and 30◦ UD 0.32±0.1440◦ E and 10◦ RD 0.33±0.1040◦ F and 10◦ RD 0.32±0.122.6 SummaryIn this chapter, I have presented a novel statistical wrist shape+scale+posemodel that captures shape, scale and pose variations of 30 subjects (botharm) at 9 different wrist positions. A group-wise registration frameworkhas been used for a better correspondence establishment during the modeldevelopment. The pose model captures pose variations not only due topopulation variation but also across different wrist positions of a particularindividual. The model has been registered to target wrists of different wrist292.6. Summarypositions using LOSO and LOPO tests, and the results are promising enoughto use the model in different applications.30Chapter 3Automatic segmentation ofwrist bones using statisticalmodel registration to CT3.1 Introduction and backgroundIn the last decade, an increasing interest has been noticed about wrist bonesegmentation in Computed Tomography (CT) volumes due to its potentialapplications to various diagnostic procedures [19, 44, 52, 103]. For evalua-tion of arthritis, the size of the bone is an important indicator as the wristbones are eroded gradually during the course of the disease [52]. A seg-mentation of the wrist bones can help radiologists estimate the bone sizeaccurately. Moreover, wrist bone segmentation has been used in bone ageassessment [103] that is important to evaluate the stage of skeletal matu-ration of a child. The assessment is often done by determining the bonegrowth of the wrist joint [103], and a segmentation of the wrist joint helpsto assess the bone age. In addition, wrist bone segmentation has been usedto simulate the complex wrist motion, which helps in understanding the nor-mal wrist function and allows for detecting kinematic abnormalities such asthe scaphoid-lunate dissociation [19, 93]. Furthermore, automatic segmenta-tion of the wrist bones can be used in developing image-guided interventiontechniques for wrist fracture fixation such as guiding a percutaneous non-displaced scaphoid fracture fixation, where the segmentation of the wristbones in CT could be used to augment a 3D anatomical model of the wristin 3D intraoperative ultrasound [7].The anatomy of the wrist bones poses some challenges to their segmen-tation in CT volumes. The inter-bone distances are narrow relative to theThis chapter is adapted from Anas E. M. A., Rasoulian A., Seitel A., Darras K.,Wilson D., John P. St., Pichora D., Mousavi P., Rohling R. N., and Abolmaesumi P.,Automatic Segmentation of Wrist Bones in CT Using a Statistical Wrist Shape+PoseModel, IEEE Transactions on Medical Imaging, 13(11): 2025-2034, 2016.313.1. Introduction and backgroundtypical slice thickness (1.0-2.0 mm) of conventional CT imaging [86]. Thismakes it difficult to correctly delineate the bone boundaries. In addition,the outer layer (cortical bone) of the bone tissue is found denser than theenclosed spongy bone and as a result, the cortical bone appears brighterthan the spongy bone in CT images. As a result, strong gradients are ap-peared on both sides of the cortical surface. Moreover, due to the closespacing of wrist bones and inherent blurring in CT imaging, the inter-bonespaces often appear brighter than the image background which leads to apoor boundary contrast. All these challenges observed in CT imaging makethe segmentation of wrist bones in CT images a challenging task [86]. Ac-curate segmentation of the wrist to date involves a lot of manual work andis a rather time consuming procedure. A robust and accurate automaticsegmentation approach would thus be of great benefit for the involved radi-ologists.There is a significant number of works published in the literature fa-cilitating segmentation of wrist bones in CT images, though most of themare semi-automatic. Sebastian et al. [86] provided a comparative analysis ofvarious fundamental (e.g., intensity based, edged based, region growing anddeformable) methods for segmentation of wrist bones in CT images. Theauthors showed that these basic techniques are insufficient for this particularsegmentation task, and they proposed an improved segmentation techniquenamed ‘skeletally coupled deformable model’ to segment the wrist bones inCT images. They used a curve evolution approach from manually initial-ized seeds, where growth was controlled by skeletally-mediated competitionbetween neighboring regions. Duryea et al. [23] reported a semi-automatedsegmentation tool using a combination of edge-tracking and active contourmethods to segment the wrist bones for a study of rheumatoid arthritisprogression. The reported technique at first segmented the center slice ofa particular wrist bone by using an edge tracking method. The segmen-tation at the center slice was then considered as an initial contour to sub-sequent slices and an active contour method was utilized to segment thebones. In addition, they used a gradient threshold refinement algorithmfor further improving the segmentations. Another interactive segmentationframework [17] combined the grow-cut [94] segmentation method with rigidimage registration to segment and align the wrist bones simultaneously fromdifferent subjects or at different poses from the same subject. However, theabove mentioned semi-automatic methods required expert users to achieveacceptable segmentation results.For development of automatic segmentation techniques, the use of statis-tical models has been suggested in the literature [19, 21, 54, 71, 79]. A sta-323.2. Methodstistical model allows incorporation of a priori information about the shapesand poses of objects that can then be used during the segmentation process.However, few works have been published for segmentation of complex hu-man joints, consisting of multiple structures such as the wrist [19] and thespine [79].Most recently, Chen et al. [19] automatically segmented the radius, ulnaand eight carpal bones from a set of wrist CT volumes using a successivelydeveloped statistical pose model, where the required correspondences aredetermined using a pair-wise registration approach. The statistical wristpose model captured the major two (corresponding to radial-ulnar deviationand flexion-extension) modes of pose variations across a training set at 5different wrist positions. A set of 125 CT volumes from 25 subjects at 5different wrist positions were used for development of the pose model andevaluation of the segmentation through a leave-one-out-subject experiment(with all corresponding poses). The authors reported approximately 20% ofunsuccessful segmentations mainly due to the misalignment of symmetricalulna and trapezoid bones.In this chapter, I register my developed statistical shape+scale+posemodel to target CT volumes to enable fully automatic segmentation ofwrist bones in the target CTs. The results of the segmentation will besubsequently used in the intraoperative registration step (chapter 5). Theregistration of the multi-object statistical model to a target CT volume isperformed using a three-step registration technique that tackles challengespresent in CT wrist segmentation such as non-uniformity of the bone tis-sue and narrow low contrast inter-bone spaces. The proposed integratedsegmentation framework is evaluated on a total of 66 CT volumes.3.2 MethodsSegmentation of the wrist bones in a CT volume can be seen as a registrationproblem of the wrist model to CT, and it can be broken down into an initialalignment (Sec. 3.2.1), a multi-step wrist model registration step (Sec. 3.2.2)and a non-rigid registration (Sec. 3.2.3) at the end. An overview of mymethod registering the wrist model to an unseen CT volume is shown inFig. 3.1.3.2.1 Initial alignmentAn initial rigid registration of the wrist model and the target CT volume isnecessary to prevent the following model-based registration from converging333.2. MethodsRepetitive model registrationThresholdingWrist modelTarget wrist point cloud1) Initial rigid alignmentRegistered Model (initial)2.i) Model optimizationRegistered model (iteration j)2.ii) Target wrist refinementRefined  target wrist point cloudj > total number of repetitions?3) Non-rigid registrationRegistered model (final) overlaid on CTnoyesCT volumeFigure 3.1: Flow chart of whole registration procedure. The CT volume isthresholded to create a rough wrist structure. It is used to initially alignthe wrist model. An iterative approach is followed to refine the wrist struc-ture and estimate the model coefficients. In the end, non-rigid registrationapproach is used to correct the residual errors.into a local minimum. For this purpose, I roughly segment the boundary ofthe wrist in the CT volume by applying Otsu’s thresholding technique [69].Fig. 3.1 shows a wrist surface resulting from the thresholding operation, andit is evident that the shapes of the wrist bones are quite visible. The ini-tial registration is then performed by aligning the corresponding principalaxes of the mean model and the segmented target CT surface. As coordi-nate axes order or direction might not be determined consistently across themean model and the target point cloud, I choose the best initial alignmentfrom the set of all 23 = 8 (2 possible directions at 3 possible orders) possiblealignments. For each possible axes configuration, I first transfer the meanmodel to the target space by aligning the coordinate axes using GeneralizedProcrustes analysis, and then perform a probabilistic rigid-body (includingscale) registration [67] between the transformed mean model and the targetpoint cloud for fine registration. The alignment resulting in the smallestσ (please see Eqn. 2.10 for definition of σ) value of the probabilistic rigidregistration which represents how well the registered model corresponds to343.2. Methodsthe target point cloud is chosen as the initial alignment.Note that in the absence of a priori information on the wrist positionin the target CT, I use the mean model for initializing the registration. Ifapproximate pose information is known, the registration can be initiatedwith a model instance at that pose for better initial alignment.3.2.2 Model registrationFor model registration two main processing steps are repeated multipletimes. First, the model coefficients are optimized with respect to the targetpoint cloud as extracted from CT (i). The number of shape (Cˆsj ) and pose(Cˆpj and Cˆp′j ) coefficients of the model are evenly increased throughout theregistration loop to avoid the optimization converging in local minima dueto the initially coarse segmentation of the bone surface in the CT volumeor a poor initial registration. In other words, the gradual increase of thenumber of coefficients during registration can be seen as a gradual coarseto fine registration. At first, I emphasis on coarse registration, i.e., quicklyconverging to a close solution, that can be further refined using higher ordersof eigenmodes. This approach also helps with improving computation timeas it reduces the search space for registration parameters. The optimizedmodel is then used in the second processing step to refine the target pointcloud (ii). This process is repeated until the user defined maximal numberof repetitions (JR) is reached.i) Model optimization Model optimization is performed using a GMM-based registration technique [80] in which the target point cloud is assumedto be an observation generated by a GMM, and the Expectation Maximiza-tion (EM) framework is applied to optimize the model parameters. Thedetails of the registration was presented in the previous chapter.ii) Refinement of target point cloud The optimized model is usedto then refine the target point cloud after each iteration especially in theinter-bone spaces where the initial CT segmentation could not provide re-liable target points. A gradient-based approach is used for the refinementprocedure. The gradients of the target CT volume are computed along theoutward normal directions for each point of the partially registered modelobtained after each step. Since the intensity of the CT image decreases fromhigh to low values along the outward normal direction of the desired bone353.3. Experiments and evaluationboundaries, the points corresponding to the strongest negative gradient val-ues across the outward normal direction of the partially registered modelconstitute the refined target point cloud.3.2.3 Non-rigid registrationIn the final step I use a probabilistic non-rigid registration technique [67]to correct the residual registration errors that cannot be captured by theproposed model. In short, the non-rigid registration considers one pointset as GMM centroids and the other as the data points. The registrationis achieved by moving the GMM centroids coherently to the data point bymaximizing the likelihood function. To avoid overlapping segmentations ofneighboring bones, I introduce an iterative refinement of the target pointcloud into the non-rigid registration approach. After each iteration, I checkfor occurring intersections between the registered bone surfaces. Then, Iupdate the target point cloud by choosing the refined target points at thenext strongest negative gradient in the volume along the inverted surfacenormal at the corresponding point.Finally, the registered model to CT for the l-th wrist bone can be ex-pressed asYl,CT (ΘsCT ,ΘkCT ,Sl,CT ,ΘpCT ,Θp′CT ) = (3.1)φkl (ΘkCT )Φsl (ΘsCT ) + Sl,CT } × {Rp′l (ΘpCT ,Θp′CT )}T+IM × {tp′l (ΘpCT ,Θp′CT )}T ,Here, ΘsCT , ΘkCT , ΘpCT and Θp′CT represent the shape, scale, inter-subjectand intra-subject model coefficients, respectively, estimated from the sta-tistical model registration to CT. Sl,CT represents the residual registrationerror estimated using the non-rigid registration technique.3.3 Experiments and evaluationFor evaluation and selection of tuning parameters of the presented model-based segmentation approach, I collected a total of 81 wrist CT volumesfrom two different hospitals: Vancouver and Kingston General Hospitals,using different CT machines. Note that these 81 CT volumes were com-pletely independent from the data [65] (mesh models of the wrist bones) weused for the model development in the previous chapter. To obtain a refer-ence for evaluation, all 81 CT volumes were semi-automatically segmented363.3. Experiments and evaluationusing the interactive segmentation framework of the Medical Imaging In-teraction Toolkit (MITK) [59]. Three biomedical engineering students wereappointed to perform the semi-automatic segmentations, and these studentswere trained by an expert postdoctoral fellow with many years of experiencein image guided interventions. In addition, to minimize inter-observer vari-ability the STAPLE algorithm [96] was used to combine the independentsegmentation data. Moreover, to select the trade-off variables of the seg-mentation technique, the 81 CT volumes were arbitrarily categorized intotwo groups. The cross-validation set was chosen to consist of 15 CT volumesand was used to select the trade-off variables of the algorithm. The testingdata set contains the remaining 66 CT volumes including the fractured casesand scans at extension position, and was used to measure the accuracy ofthe CT segmentation technique. A description about the evaluating param-eters and optimization of the tuning parameters are presented in the nexttwo sections.3.3.1 Evaluating parametersTo evaluate the segmentation quality quantitatively, I calculate the meansurface distance error (mSDE) and the Jaccard index (JI) [46] from eachsegmentation and its reference. SDE is defined as the distance betweeneach point of the registered surface and its closest neighboring point in thereference CT segmentation, and mSDE is the average of SDEs across allvertices. The JI, very similar to Dice coefficeint, measures the similaritybetween two sample sets, and is defined as the size of their intersectiondivided by the size of their union. The coefficient is defined between [0,1],where 1 indicates a perfect overlap of the compared segmentations.In addition, to evaluate the computational complexity of the proposedsegmentation method, I measure the computation time of complete segmen-tation algorithm as well as its different processing steps. The computationtime is recorded from unoptimized MATLAB code running on an Intel Corei3-2310M CPU at 2.17 GHz.3.3.2 Parameter selectionThe proposed algorithm for segmenting wrist bones in unseen CT volumesdepends on several parameters that need to be defined to achieve sufficientsegmentation quality. The main parameter dependencies of the approachare the number of points in the target point cloud (N), and the number ofrepetitions in the model registration cycle (JR). To analyze the influence of373.4. Resultsthese parameters on the segmentation outcome, I use the cross-validationdata set consisting of 15 CT volumes and perform the segmentation algo-rithm for each of these volumes while varying N and JR. Fig. 3.2 shows thedependency of both segmentation accuracy and computation time on theseparameters indicating an increased accuracy and an increased computationtime with both increasing N and JR. Note that either N or JR are keptconstant at high values (N=20000 and JR=6) while varying the respectiveother parameter. Given that the increase in accuracy is neglectable abovecertain values of JR and N I chose JR = 3 and N = 12000 for my evaluation.The total number of model coefficients are chosen such that they cover95% of the shape and pose variations in the model. From the compactnessanalysis (in the previous chapter) I get that I need a total of 50 shape, 2 scale,8 interpose and 2 intra-subject pose coefficients to cover 95% variations. Inthe model registration process I evenly increase the number of coefficientsthroughout the JR = 3 repetitions such that I finally get 17 shape, 2 scale,3 inter-subject pose and 2 intra-subject pose coefficients for repetition 1; 34shape, 2 scale, 6 inter-subject pose and 2 intra-subject pose coefficients forrepetition 2; and 50 shape, 2 scale, 8 inter-subject pose and 2 intra-subjectpose coefficients for repetition 3. Note that I will not increase the number ofscale and intra-subject pose coefficients in each step, as I only have 2 scaleand intra-subject pose coefficients capturing 95% of the variation.3.4 ResultsQuantitative results of the model-based wrist segmentation on the 66 CTvolumes are summarized in Table 3.1. An average mSDE of 0.33±0.11 mmbetween the resulting wrist model and the reference segmentation averagedover all 66 CT volumes and all 9 wrist bones could be achieved. Median,Haussdorf and RMS error were recorded at 0.32 mm, 1.19 mm and 0.36 mm,respectively. The average similarity between the calculated segmentationand the reference segmentation is recorded with a JI of 0.87. Most of thetarget CT volumes could be segmented with an average mSDE smaller thanthe average voxel size (Fig. 3.3). Fig. 3.4 visualizes the segmentation resultsof some example slices.Applying the non-rigid registration step leads to an improvement of31.87% of the average mSDE (Fig. 3.5). Figs. 3.6(a-b) demonstrate thequality of my segmentations with respect to the ground truth segmentationat palmar and dorsal views for all wrist bones. The colors on the bone383.4. Resultsaverage mSDE (mm)average JIcomputation time (min)average JIcomputation time (min)(a) (c)average mSDE (mm)(b) (d)1 2 3 4 50.320.340.360.380.40.420.650.70.750.80.850.92000 4000 6000 8000 10000 12000 14000 160000.340.360.380.40.420.70.750.80.851 2 3 4 51012141618202000 4000 6000 8000 10000 12000 14000 160001213141516171819Figure 3.2: Influence of the number of repetitions (JR) and the numberof points (N) in the target point cloud on the registration accuracy andcomputation time. Note that when I vary JR to compute the affected indices,I keep N constant at 20000 and vice versa I keep JR constant at 6. ThemSDE is computed across all vertices of the registered model. For eachJR and N , mSDEs and JIs between the resulting wrist models and thereference segmentations are averaged over all 15 CT volumes. In contrast,the computation time is computed for an arbitrary CT volume for each JRand N . (a) Average mSDE and JI vs. JR. (b) Computation time vs. JR.(c) Average mSDE and JI vs. N . (d) Computation time vs. N .surfaces represent the SDEs with respect to the ground truth segmentation.Higher errors are mostly found on the capitate and hamate bone surfaces.The average overall computation time was recorded as 17:58 min where0:25 min were used for the initialization step, 14:46 min for the model reg-istration step and 3:12 min for the non-rigid registration step.393.5. DiscussionAverage voxel size (mm)mSDE across all bones (mm)voxel size = mSDE 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.60.150.20.250.30.350.40.450.50.55Figure 3.3: Dependence of the mSDE of the average voxel size in the CTvolumes. The mean SDE is computed across all vertices of the registeredmodel to a CT, and the average voxel size is calculated as the mean of thevoxel dimensions along x-, y- and z-axes in that CT volume. Results from 66CT volumes are used in this figure and each point represents the coordinateof (average voxel size, mSDE) for a particular CT volume.3.5 DiscussionThese results are comparable with results achieved by the state-of-the-artapproach proposed by Chen et al. [19] that reported a mean segmentationaccuracy of 0.34 mm achieved on a subset of the data labeled by the al-gorithm as ‘successfully registered’. The authors pointed out limitations oftheir method especially for symmetrical wrist bones such as the trapezoid.Introduction of the inter-subject pose variations in my proposed statisti-cal model of the wrist allows us to avoid misalignment of those symmetricalcarpal bones. In addition, my developed model is based on a group-wise reg-istration framework compared to the pairwise registration framework usedin Chen et al. [19]. A group-wise registration approach establishes bettercorrespondences across the population compared to a pair-wise one as thelatter is biased by the template selection [5]. Improved correspondencesin turn allow my statistical model to capture anatomical variations more403.5. Discussion(a) (b) (c)ScaphoidTriquetrumLunateCapitateTrapeziumHamateTrapezoidRadius(d) (e) (f)axialsagittalcoronalaxialaxialsagittalFigure 3.4: Segmented wrist bone boundaries superimposed on example CTslices at different planes. Each color represents a particular bone and theircorrespondences are given in the figure. (a-c) Example segmented axial,sagittal and coronal CT slices of a healthy individual at neutral wrist posi-tion. (d-e) Two example axial CT slices showing segmentations of fracturedtrapezium (d) and lunate (e) bones. (f) Example segmented sagittal CTslice of a cadaver at extension position.accurately.The statistical model used for the segmentation has been developed usinga training set of healthy individuals. Since most isolated wrist fractures arenot associated with significant carpal disruption and pose perturbation, Ianticipate that the statistics captured in the model is also applicable for caseswith such fractures. To further investigate the potential of my approach forsegmenting fractured CT cases, I applied the approach on the two fracturedcases in our CT dataset. Figs. 3.4(d-e) show the promising performance ofthe approach on those cases.The proposed segmentation technique could also be applied to segmentthe wrist bones at non-neutral positions. Example segmentation results isdemonstrated in Fig. 3.4(f). Overall, I could achieve a mSDE of 0.31 ± 0.11mm and a JI of 0.88 ± 0.05, calculated from four non-neutral test cases.The proposed approach was able to segment the individual wrist bonesand could successfully delineate the bone surfaces despite the narrow inter-bone distances for most of the data sets (Fig. 3.4). However, it had dif-413.5. DiscussionTable 3.1: Statistical measures of the mSDEs and JIs of nine wrist bonescomputed from 66 CT images (Radius= Rad, Scaphoid= Sca, Capitate=Cap, Hamate= Ham, Lunate= Lun, Pisiform= Pis, Trapezium= Trzm,Trapezoid= Trzd, Triquetrum= Triq.)mSDE (mm) Rad Sca Cap Ham LunMedian 0.34 0.33 0.30 0.33 0.40Haussdorf (maximum) 0.74 0.50 0.50 0.54 1.19Average 0.34 0.33 0.30 0.33 0.40Standard deviation 0.13 0.12 0.09 0.10 0.18RMS 0.37 0.33 0.33 0.38 0.44mSDE (mm) Pis Trzm Trzd Triq AllMedian 0.31 0.28 0.42 0.32 0.32Haussdorf (maximum) 0.52 0.52 0.63 0.79 1.19Average 0.30 0.30 0.37 0.33 0.33Standard deviation 0.10 0.12 0.10 0.11 0.11RMS 0.33 0.33 0.36 0.38 0.36JI Rad Sca Cap Ham LunMedian 0.87 0.89 0.88 0.86 0.84Average 0.88 0.89 0.89 0.85 0.84Standard deviation 0.06 0.05 0.05 0.04 0.08RMS 0.88 0.89 0.88 0.84 0.82JI Pis Trzm Trzd Triq AllMedian 0.84 0.87 0.85 0.87 0.87Average 0.85 0.88 0.83 0.82 0.87Standard deviation 0.05 0.04 0.04 0.05 0.05RMS 0.86 0.87 0.85 0.87 0.86423.5. Discussionwith non-rigid without non-rigid mSDE across all bones (mm)0.20.40.60.811.21.41.61.82Figure 3.5: Effect of the non-rigid registration step on the registration error.I compute mSDE across all vertices of each registered surface. Two types ofsegmentations are performed here. The first one is my proposed approachand in the second one I avoid the non-rigid registration step from my pro-posed segmentation technique. The median value of the mSDE across 66CT samples is marked by the red color. The solid blue represents the rangebetween 25 and 75 percentiles. The red crosses indicate outliers.ficulties in finding the correct bone outline in 11 CT volumes, where wegot mSDE greater than average+standard deviation (0.44 mm). The mostcommon reason for this segmentation inaccuracy is the higher CT slice thick-ness which results in weaker gradients along the bone boundaries due to thepartial volume effect. Though the iterative target point cloud refinement ap-proach in my method is robust to overcome this effect, however, especially inthe areas where neighboring bones are very close to each other, bone edgesget merged due to the partial volume effect which leads to locally wrongassignments of points to the target point cloud. I could show that with in-creasing image resolution this issue becomes less prominent and registrationaccuracy increases (Fig. 3.3)). For applications with high accuracy require-ments, acquisition of high resolution CT data should be advised. Moreover,some of my test cases showed high curvatures in the bone boundaries that433.5. Discussion00.050.10.150.20.250.30.350.40.450.5(b)(a)(c) (d) (e)axial sagittal coronalaxialcoronalsagittalFigure 3.6: Error distributions across all vertices of the registered model.Segmented bone surfaces of all wrist bones of an arbitrary CT volume atpalmar (a) and dorsal (b) views. The color represents the error in mm of thesegmentation with respect to gold standard. (c-e) Axial, sagittal and coronalslices demonstrating the errors in the corresponding CT planes marked in(a).were not covered by the model, and could thus not be segmented correctly(Figs. 3.6(c-e)). Including more of those distinct cases in the model trainingprocess is likely to improve segmentation performance for similar such cases.For example, we can include these 81 CTs along with the wrist database [65]to develop a statistical model that possibly includes more variations in thedeveloped model. Note that these 81 CTs were acquired only at one wristposition. Therefore, only the shape and size models could possibly obtainbenefits from this extension of dataset.The performance of my model-based wrist segmentation approach de-pends on various factors that will be discussed in the following. The currentinitialization of the model is dependent on the Otsu threshold segmentationof the CT volume. Poor segmentations due to noise or image artifacts pro-duced e.g. by implants or screws might result in initial misalignment since443.6. Summarythe computed principal axes do not correspond to the anatomical geometryof the wrist anymore. Defining a scanning protocol that keeps the wrist ina defined fixed orientation for each scan would reduce these issues.The target point cloud is currently subject to a repetitive refinementwhich results in an increased computation time. Inclusion of the boundaryappearance in e.g. an appearance model could be interesting for furtherimprovement of segmentation performance.Achieving accurate segmentations currently relies on the number of repe-titions in the model registration step (Sec. 3.2.2.ii) and the number of pointsin the target point cloud to correctly align the individual bones of the wristmodel to their respective bone surfaces in the CT volume. A careful selec-tion of the number of repetitions and the number of points in the targetpoint cloud has been performed to achieve a balance between segmentationaccuracy and computational performance (Fig. 3.2).The gold-standard segmentation of the 81 CT volumes was achievedthrough semi-automatic segmentations by three biomedical engineering stu-dents under the supervision of an expert postdoctoral fellow. Furthermore,we used the STAPLE algorithm [96] to minimize the inter-observer variabil-ity. It should be noted that a randomly selected subset of the segmentationswere supplied to an expert radiologist for evaluation of the reference seg-mentations, who confirmed the quality of the gold-standard segmentations.3.6 SummaryIn this chapter, I have presented an application of the statistical wrist modelto automatic wrist bone segmentation in CT images. Segmentation of wristbones in an unseen target CT image is performed by applying this modelin a multi-step registration procedure. The main registration step is basedon a probabilistic approach robust to noise, outliers and missing points thatallows the model to repetitively register to a target point cloud representingthe bone surfaces in the CT image. Residual registration errors have beencorrected by including a final non-rigid registration step. I applied my seg-mentation approach to 66 unseen CT volumes and could achieve sub-voxelaccuracy (average mSDE of 0.33 mm) and an average JI of 0.87. In conclu-sion, the promising segmentation result underlines its potential for differentclinical applications.45Chapter 4Bone enhancement in USusing local spectrumvariation4.1 Introduction and backgroundA comprehensive investigation has been performed by a significant number ofrecent past studies regarding potential of US to guide orthopaedic surgeries,and for percutanous pain management [2, 20, 53, 72, 91, 102]. However, themain challenge of the US-based approaches is the detection/enhancement ofthe noisy, disconnected, low-contrast and blurry bone surface in US image.The previous works ondetection/enhancement of the US bone responses can be broadly catego-rized into two groups: intensity-based [2, 25, 72, 98, 102] and phase-basedapproaches [38–41]. Intensity-based approaches use the intensities of thebone responses as well as the shadow below the bone surfaces to enhancethe bone features. In contrast, phase-based approaches mainly rely on theimplicit features of the bone responses, e.g., local symmetry of the boneresponses.Penney et al. [72] proposed to enhance the pixel response of the bonyparts in the US image by defining a measure based on the pixel intensityand the ratio of the number of zero pixels between the pixel of interest andthe last pixel of the corresponding column to the total number of pixelsbetween them. Both measurements were expected to have high values forthe bony part because of the large change in acoustic impedances betweensoft-tissue and bone. Amin et al. [2] used a directional edge detector todetect the transition between the bone shadow region and the bone surfaceThis chapter is adapted from Anas E. M. A., Seitel A., Rasoulian A., John P. St.,Pichora D., Darras K., Wilson D., Lessoway V. A., Hacihaliloglu I., Mousavi P., Rohling R.N., and Abolmaesumi P., Bone enhancement in ultrasound using local spectrum variationsfor guiding percutaneous scaphoid fracture fixation procedures. International Journal ofComputer Assisted Radiology and Surgery, 10(6):959-969, 2015.464.1. Introduction and backgroundreflection. They found the pixels greater than an empirically set bone surfacethreshold in each column along the upward direction. The output of thedirectional edge was convolved with a Gaussian kernel to ensure that theedge detector region captured the entire bone surface reflection. Yan etal. [102] used a backward tracing method that was similar to the directionaledge detector, except the idea was based on a fan-shaped sector of the phasedarray probe. Kowal et al. [53] presented an automatic US bone contourdetection approach that was subdivided into two main steps. Firstly theyidentified a sub-image range that most likely contained the bony part inthe US image. Then, they searched for the bone contour within the sub-image region. Foroughi et al. [25] modelled the shadow below a pixel thatbelongs to the bony part of the US image as a Gaussian weighting function,and subsequently calculated a bone probability map in an US image. Theythen proposed a cost function that includes continuity and smoothness ofthe bone surface, and dynamic programming was used for finding the globalminima in order to detect the bone contour. Most recently, Wein et al. [98]presented an intensity-based automatic bone contour detection approachbased on a bone-specific feature descriptor. Initially, non-vertical ridge-like responses in US were enhanced using an orientation-dependent filterapplied to the input US image. Next, an ultrasound confidence map [49]was used for enhancement of the expected bone locations while suppressingthe responses from other anatomical interfaces. However, the main challengeof all of the described intensity-based bone contour detection approaches isthe sensitivity of their methods to US artifacts and machine settings [39, 41].For example, the assumption of high intensity values along the bone regionis not always true, particularly when the US transducer is not positionedperpendicular to the bone surface. This issue is often noticed for curved wristbone surfaces, which in turn leads to low contrast, blurred and disconnectedbone features [47].To address the above mentioned challenges, local phase-based US imageenhancement approaches were recently suggested. They were reported forbone enhancement [38–41] as well as localizing soft tissue interfaces [35, 74,75]. In the local phase-based approach, the phase symmetry (PS) informa-tion was utilized to detect the ridge-like features in US images, using a setof quadrature filters at different scales and orientations. Hacihaliloglu etal. [38] used Log-Gabor filters to obtain the PS information from a 2D USimage that showed high values at the expected bone locations. The sameauthors extended their 2D approach to 3D US by designing 3D Log-Gaborfilters [40]. Most recently, Hacihaliloglu et al. [41] proposed a Local PhaseTensor metric for enhancing bone features in US images. The tensor was474.2. Methodscomputed using an isotropic band-pass (Log-Gabor) filtered image, its gra-dient and Hessian. The filter parameters (scale, bandwidth, orientation andangular bandwidth of the filter) were estimated from the image content byincorporating the principal curvature computed from the Hessian matrixand directional filter banks in a phase scale-space framework [40]. The filterdesign considered an isotropic scale and bandwidth of the filters at all orien-tations, i.e., the frequency responses of the filters are the same throughoutall orientations. However, the frequency response of a curved surface in USshows variations over different orientations. As a result, this approach maynot detect/enhance this kind of surface without accounting for the spectrumvariations over different orientations.In this chapter, I propose improved techniques for enhancing bone sur-faces in both 2D and 3D US images based on an empirical wavelet approachthat exploits the local spectrum variation of the ultrasound image. Em-pirical wavelet filters have been proposed in the literature to extract thedifferent modes of a signal according to its contained information [32, 33].I propose to account for the local spectrum variations to be able to designa filter bank for subsequent local PS estimation. Bone shadow informationis also utilized to better differentiate the bone surface from the soft tissue,nerves and other anatomies in US images. For qualitative evaluation, phan-tom, ex vivo and in vivo data are acquired for the feasibility experiments,and the results are compared to two previously published methods [38, 40].4.2 MethodsThis section starts with how the bone responses in US give rise to anisotropicvariation in the frequency spectrum (Sec. 4.2.1), followed by the introduc-tion of quadrature Log-Gabor filters for 2D and 3D local phase symmetryestimation in Sec. 4.2.2. Finally I will present bone enhancement techniquesfor 2D and 3D US images (Sec. 4.2.3). Though my thesis is based on us-ing 3D US as guidance, I also propose a 2D bone enhancement techniqueparticularly for different 2D US applications.4.2.1 Spectrum variation of US bone responsesA bone surface in an US image produces highest intensity values when theUS beam falls perpendicularly to that bone surface, otherwise, the responsefrom that bone surface will be weaker. Since an US signal cannot prop-agate through the bony structure due to the large mismatches in acousticimpedances between the bone and neighboring tissues, shadows are observed484.2. Methodsbelow bone surfaces. The bone response above the shadow has a certainwidth, reaching up to 4 mm [47]. In addition, blurred and disconnectedbone responses are also often observed.To understand the variation of frequency responses over the spectrum,I present some examples of 2D Fourier spectrums of idealized bone-like re-sponses. Figs. 4.1(a-c) show three images representing a horizontal line, aninclined line and a curved structure. The corresponding 2D Fourier magni-tude spectrums are demonstrated in Figs. 4.1(d-f) on a logarithmic scale. Inotice that a horizontal structure in the image domain produces a responsealong the vertical axis in the Fourier domain. In addition, the width of theresponse is different in the horizontal and vertical directions, and the widthalong the vertical direction is significantly greater than that along the hor-izontal direction. An inclined line such as in Fig. 4.1(b), also produces aninclined line-like structure in the Fourier domain. It is worth noting that theFourier spectrum of Fig. 4.1(b) is rotated with respect to that of Fig. 4.1(c)by the same angle difference between the inclined and the horizontal lines.Here, most of the energy is concentrated along the inclined major axis. Fi-nally, a curved structure in the spatial domain also creates a curved-likestructure in the Fourier domain, and the zero crossings of the lobes are notlocated along the trajectory of a circle. Moreover, the energy is more spreadout in the Fourier domain compared to the line-like structure.The above examples reveal that the energy of the Fourier spectrum dis-tributes non-uniformly all over the region. Furthermore, the width of thelobe is different in different directions. Though the presented examples arebased on 2D US images, the non-uniform distribution of the spectrum vari-ations also holds true for a 3D US volume.4.2.2 Phase symmetry estimationThe local phase symmetry (PS) estimation starts with dividing a 2D/3Dfrequency spectrum into different orientations (example orientations inFig. 4.2(a-b), and subsequent application of the Log-Gabor filters at eachorientation. The estimation of PS using a set of Log-Gabor filters are pre-sented in the next two sections for 2D and 3D US images.494.2. Methods−1 −0.5 0 0.5 1−1−0.500.51−1 −0.5 0 0.5 1−1−0.500.51−1 −0.5 0 0.5 1−1−0.500.51/pi /pi/pi/pi/pi/pi(a) (b) (c)(d) (e) (f)Figure 4.1: Simulated US bone responses and their spectrums. (a-c) Threesimulated bone responses, (d-f) corresponding 2D Fourier transformations.2DTo divide the spectrum into different orientations, a set of orientationalfilters are used, where the frequency response of each filter is defined as:O(φ) = exp(−(φ− φ0)22σφ), (4.1)where, φ represents the azimuthal angle in polar coordinate, φ0 indicatesthe center of the orientation, and σφ represent the span/bandwidth of theorientation (Fig. 4.2(a)).After applying the orientational filter, a set of Log-gabor filters are usedat each orientation. Mathematically, the frequency response of a 2D Log-Gabor filter is defined as below:R(ω) = exp(−(ln(ωω0))22(ln(κ))2), (4.2)where ω (0 ≤ ω ≤ √2pi) represents the angular frequency, ω0 representsthe peak tuning frequency, and 0 < κ < 1 is related to the octave band-width. Finally, the frequency response of a band-pass Log-Gabor filter at aparticular orientation can be expressed as: F (ω, θ) = R(ω)O(φ).504.2. Methods(b)(a)θ0ϕ0coneorientationxyzxy orientationϕ0Figure 4.2: Orientations in 2D and 3D frequency spectrums. (a) An exampleorientation in 2D frequency spectrum. (b) A 3D frequency spectrum isdivided into different cones, and each segmented cone is further partitionedinto different orientations.3DThe frequency response of an orientational filter is defined as a multiplicationof azimuthal (Φ(φ)) and polar (Θ(θ)) orientational filters:O(φ, θ) = Φ(φ)×Θ(θ) = exp(−(φ− φ0)22σφ)× exp(−(θ − θ0)22σθ), (4.3)where the azimuthal angle φ (0 ≤ φ ≤ 2pi) measures the angle in the xy-plane from the positive x-axis in counter-clockwise direction, and the po-lar angle θ (0 ≤ θ ≤ pi) indicates the angle from the positive z-axis. φ0and θ0 indicate the center of the orientation, and σφ and σθ represent thespan/bandwidth of the orientation (Fig. 4.2(a)). The purpose of the polarorientational filter is to divide the 3D spectrum into different cones, andthe azimuthal orientational filter further divides each cone into differentsub-spectrums/orientations (Fig. 4.2(a)).After selection of a particular orientation, the frequency response of aband-pass Log-Gabor filter at a particular orientation can be expressed as:F (ω, φ, θ) = R(ω)O(φ, θ), where R(ω) is defined before in Eqn. 4.7.514.2. Methods4.2.3 Bone enhancementMy proposed approaches are primarily based on accounting the local spec-trum variation by using filters with appropriate parameters at different ori-entations.2D bone enhancementIn the following, I will present the different steps of my proposed boneenhancement method which are illustrated in the flowchart shown in Fig. 4.3.2D US frame( )(1) Polar FT (2) Orientational filter(3) Log-Gabor filter(4) PS using eqn. (4.8)(5) Bone enhancement( )Figure 4.3: The flowchart of the proposed 2D bone enhancement approach.(1) The corresponding Fourier transform P (ω, φ) in the polar domain iscomputed from the US image. (2) Estimation of the orientational filterparameters and orientational filter design. (3) Log-Gabor filter design. (4)The local PS is estimated from the designed filter and the US image. (5)After that, I obtain BR image from PS using the approach proposed byForoughi et al. [25].524.2. Methods(1) Calculation of the polar Fourier transform For designing thefilters, I use the polar Fourier transform instead of the conventional Fouriertransform computed over a rectangular grid, because the parameters of thequadrature filters are directly related to the coordinates of the polar domain.Let U(x, y) be the input US image and P (ω, φ) be the corresponding Fouriertransform in polar domain. Nw and Nφ represent the number of radius andazimuthal bins in the polar Fourier transform.(2) Orientational filter For estimation of the parameters of the orienta-tional filters O(φ), I average the logarithmic magnitude of P (ω, φ) over ω,p(φ) =1Nω∑ωlog(|P (ω, φ)|), (4.4)where a higher value in p(φ) indicates the presence of high energy at thatparticular orientation (Fig. 4.1(i)). To segment p(φ), I use the idea proposedin [33]. The segmented portions of p(φ) correspond to different orientationsof the filter, i.e., the total number of orientations is equal to the numberof segmented portions. Each separate portion is centered around a specificangular value φo (a maximum), and two neighboring supports/boundaries(minima or sufficient low magnitude), represented as φ+o and φ−o (examplesare shown in Fig. 4.4(a)). These values are subsequently used for orienta-tional filter design as:Oo(φ) = exp(− (φ− φo)22(φ+o − φ−o )). (4.5)(3) Log-Gabor filter In contrast to the previously published approaches[38, 39, 41], the parameters of the Log-Gabor filters for each different ori-entation are computed separately. For each orientation o, I compute theaverage of the logarithmic magnitude of P (ω, φ) over φ within that orienta-tion range (φ−o to φ+o ),qo(ω) =1Noφ∑φlog(|P (ω, φ)|), (4.6)where Noφ is the number of angle bins within φ−o and φ+o . Fig. 4.4(b) showsan example qo(ω). To estimate the parameters of the Log-Gabor filter, Isegment qo(ω). Let the neighboring radial supports (cut-off frequencies) fora particular orientation o be ω−s,o and ω+s,o (Fig. 4.4(b)), and the total number534.2. Methods0.5 1 1.5 2 2.56.86.977.17.27.37.47.57.610−210−1100−70−60−50−40−30−20−100(a) (b)Figure 4.4: Utilization of the spectrum variation in the filter design. (a)The variation of spectrum strength over the polar angle. The maxima andneighboring two minima are detected from the spectrum strength. (b) Thevariation of spectrum strength over the angular frequency. The spectrum isdivided into different parts and the lower and upper cut-off frequencies aredetermined from each part.of segmented portions of qo(ω) be equal to the total number of scales in thatparticular orientation.The estimated cut-off frequencies for each scale of the Log-Gabor filters(ω−s,o and ω+s,o) are used to compute the Log-Gabor filter parameters: ωm,ns,0 =√ω−s,o ω+s,o and κm,ns = exp(−0.25 log2(ω+s,oω−s,o) √2 ln 2). The Log-Gabor filterresponse at orientation o and scale s can be represented as:Rs,o(ω) = exp− (ln( ωωm,ns,0 ))22(ln(κm,ns ))2 . (4.7)(4) Phase symmetry (PS) estimation After designing of Rs,o(ω) andOo(φ), I can compute the local PS according to the following equation:PS(x, y) =∑s∑o(|P(s,o)(x, y)| − |Q(s,o)(x, y)| − To)∑s∑o√P 2(s,o)(x, y) +Q2(s,o)(x, y) + ε, (4.8)where P(s,o)(x, y) and Q(s,o)(x, y) are the quadrature outputs of the designedfilter at scale s and orientation o, i.e., P(s,o)(x, y) = real(F−1(F(s,o)(ω, φ)) ∗U(x, y)) and Q(s,o)(x, y) = imag(F−1(F(s,o)(ω, φ)) ∗ U(x, y)). Here, F−1(.)544.2. Methodsindicates the inverse Fourier transform, * represents the convolution op-eration and F(s,o)(ω, φ) = Rs,o(ω)Oo(φ). ε (=0.0001) is a small constantintroduced to avoid division by zero error. To is an orientation-dependentnoise compensation term and it represents the maximum response that couldbe generated from the noise alone in the signal.(5) Bone response (BR) estimation Since PS is expected to have highmagnitude at symmetry locations, it also demonstrates high values at otherfeatures (e.g., soft tissue-muscle interface) in the US images [38]. Since myinterest is mainly to enhance the bony anatomy in the image, I use theshadow information of each bone in the US image. To estimate the shadowresponses, I have followed the idea proposed by Foroughi et al. [25]. In short,the shadow below an image pixel is estimated by weighted summation of theintensity values of all pixels beneath as follows:SM(x, y) =∑Hj=xG(j, y)U(j, y)∑Hj=xG(j, y), (4.9)where, SM(x, y) is the estimated shadow map (SM) for pixel (x, y) of the USimage, x and y indicate the row and column indices of the US image, andH isthe total number of rows of U(x, y). G(.) represents the Gaussian weightingfunction that models the transition of high intensity pixels near to the bonesurface to the dark pixels deeper under the bone. The product of SM withPS is defined as the bone response (BR) (BR(x, y) = PS(x, y)SM(x, y))value in my work. Since, both PS(x, y) and SM(x, y) vary from 0 to 1,BR(x, y) also has a range from 0 to 1.3D bone enhancementFor 3D bone enhancement, I initially aim to estimate the parameters for theorientational filter (φ0, θ0, σφ, σθ) for each orientation (1). The estimationof the Log-Gabor filter parameters for each orientation is presented in (2).The phase symmetry estimation and subsequent bone surface extraction aredescribed afterward (3).(1) Parameters for orientational filter The first step is to computespherical Fourier transform (FT) P (ω, φ, θ) from a given US volume. To doso, the 3D conventional FT in rectangular coordinates is calculated, followedby transforming them into spherical coordinates. For segmentation of thespectrum into different cones, I compute the strength of the spectrum along554.2. Methodsthe polar angle coordinate as: u(θ) =∑√3piω=0∑2piφ=0 log(|P (ω, φ, θ)|). Thelocations θm of the maxima of u(θ) are detected, where m = 1, 2, ...,M ,and M is the total number of detected maxima. For each θm, the detectedleft and right minima are represented as −θm and +θm, and the differencebetween these two minima positions is estimated as: σθ,m =+ θm −− θm.Note that each detected maxima corresponds to a cone in the 3D frequencyspectrum, i.e., the total number of cones is M , the center and the bandwidthof the m-th cone are θm and σθ,m, respectively.Subsequently, each segmented cone is further divided into different orient-ations/sub-spectrums. To do so, the strength of the spectrum is calculatedalong the azimuthal angle within a particular cone (say, m-th cone). Themaxima and corresponding two neighboring minima are estimated from thespectrum strength as before. Then, the center φmn and the bandwidth σmφ,nof the n-th sub-spectrum within m-th cone are calculated.(2) Parameters for Log-Gabor filters For estimation of the Log-Gaborfilter parameters at each orientation, the spectrum strength is calculatedwithin that orientation as: wm,n(ω) =∑θ∑φ20 log(|P (ω, φ, θ)|) dB. Asegmentation of wm,n(ω) is performed to estimate the parameters of theLog-Gabor filters at different scales. The lower ωm,ns,l and upper ωm,ns,u cut-offfrequencies for a scale s are determined from wm,n(ω), where, 1 ≤ s ≤ Sm,n,Sm,n is the total number of scales at n-th orientation within m-th cone. Theright subscripts ‘l’ and ‘u’ indicate the lower and upper cut-off frequencies.The parameters of the Log-Gabor filters (ω0 and κ) can be directly calculatedfrom the lower and upper cut-off frequencies as: ωm,ns,0 =√ωm,ns,l ωm,ns,u andκm,ns = exp(−0.25 log2(ωm,nu,lωm,ns,l) √2 ln 2).(3) Bone surface extraction The above estimated filter parameters arethen utilized to compute the frequency responses of the orientational andLog-Gabor filters (same as 2D bone enhancement). These filters are subse-quently used in 3D phase symmetry estimation [40]. A shadow map is thenestimated for each voxel by weighted summation of the intensity values ofall voxels beneath [25]. The product of the shadow map with the phasesymmetry is defined as the bone response (BR) in this work, which has arange from 0 to 1. To construct a target bone surface, I use a simple thresh-olding with a threshold of Tbone on the BR volume to detect the bones inthe US volume. An optimized selection of the threshold is not possible due564.3. Experiments and evaluationto a smaller sample size (13) in this thesis, therefore, an empirical thresholdvalue of 0.15 is chosen.4.3 Experiments and evaluationFor evaluation of the proposed bone enhancement approaches, I performedphantom (Sec. 4.3.1), ex vivo (Sec. 4.3.2) and in vivo (Sec. 4.3.3) experi-ments. A motorized transducer (Ultrasonix 4D L14-5/38, Ultrasonix, Rich-mond, BC, Canada) with a field-of-view of 30◦, an angular difference of0.24◦, and a center frequency of 10 MHz was used for 3D US acquisition forall experiments. A qualitative comparison of my techniques were performedin this study with two state-of-the-art approaches: a 2D local phase sym-metry (2DLPS) [37] and a 3D local phase symmetry (3DLPS) [40] methods.In addition, I recorded the run-times of the comparing bone enhancementtechniques from unoptimized MATLABTM (Mathworks, Natick, MA, USA)code on an Intel Core i7-2600M CPU at 3.40 GHz for an US volume of sizeof 57.3×36.45×32.7 mm3 with a pixel spacing of 0.4 mm in all dimensions.4.3.1 PhantomFor construction of the wrist phantom, we printed (3D) a wrist bone surfacemodel generated from a cadaver wrist CT segmentation. In addition, webrushed on the printed surface with some anti-skid powder (Rust-OleumEpoxy Shield) to make the acoustic impedance of the printed bony surfacesequal to that of the true human bone surfaces. To embed the printed bonein soft-tissue, we followed the soft-tissue preparation procedures in [58]. Thephantom construction started with addition of 50.4 g of dry gelatin and 64.5g of glass beads to 460 cc distilled water in a beaker. In addition, 16.33 gof dry agar were added to 800 cc distilled water in another beaker. The twomixtures were heated in double boilers; after the two materials had clarified(usually clarification had occurred when the temperature had risen to the90◦C range), 460 cc of the molten gelatin and 720 cc of the molten agar weremixed together in another beaker. The volume per cents are approximately40% and 60%, respectively, for all agar/gelatin materials reported here. Bypartially immersing the beaker in cold water and stirring, the mixture wascooled to 50◦C, and 18 g of Germall-plus were dissolved into it. The mixturewas then cooled to 36◦C and 3.0 cc of formalin solution was mixed in (notethat formalin solution is 37% formaldehyde). At that point, the moltenagar/gelatin was ready to be used as a soft-tissue model. Next, a portion ofthe molten agar/gelatin was poured into a mould, followed by the placing574.3. Experiments and evaluationthe printed wrist bone on it. Finally, the rest molten agar/gelatin waspoured on the printed bone for cooling and congealing overnight to roomtemperature. Note that two phantoms were constructed for my experimentusing two different printed wrist models.4.3.2 Ex vivoFig. 5.2(a) shows the general setup for US acquisition of 13 cadaver wrists.After thawing, each forearm specimen was mounted onto a custom-builtwrist holder as shown in Fig. 5.2(a). The US probe was positioned to acquirebone responses mainly from the scaphoid, part of the radius, trapezium,lunate and capitate bones. During the US scan, the probe was held free-hand on the wrist in a stationary position for at least two complete 3Dimaging sweeps.There may be some microscopic changes in the soft-tissue due to freez-ing and defrosting of the cadaver wrists, that lead to changes in the tissuedensities and corresponding acoustic impedances. However, we expect aminimal change in the reflection coefficient (in turn US image intensity) ofa particular anatomical location, as it is equal to the normalized differencebetween the acoustic impedances of neighboring interfaces.(a) (b)Figure 4.5: Experimental setup of the proposed approach. (a) A cadaverand (b) a volunteer wrist are kept fixed in our custom-built wrist holder.4.3.3 In vivoA total of 12 young volunteer wrists were scanned in the in vivo experiment.Consent was obtained from each volunteer participated in the experiment.The same wrist holder in the previous section was used to keep the wrist584.4. Resultsfixed at extension position. Fig. 5.2(b) demonstrates the general setup forUS acquisition of a volunteer wrist.It is interesting to note that the bone responses in the phantom dataappear better compared to those in the ex vivo and in vivo data. Becausethe soft-tissue model in the phantom contains almost uniform substances, incontrast, the soft-tissue in the ex vivo and in vivo wrists consist of muscle,fat, ligaments, nerves and other anatomical elements. Therefore, the USimages for the ex vivo and in vivo wrists demonstrate higher intensities inthe soft-tisue region (along with for the bony anatomies) due to mismatchesof acoustic impedances between the soft-tissue constituents (muscle or fat).In addition to demonstrating our bone enhancement results for the phan-tom, ex vivo and in vivo data, I also provide a comparison with the 2DLPS [37]and 3DLPS [40] methods. Note that only qualitative comparison is providedin this chapter as an extensive effort is needed for a quantitative compari-son to collect the gold-standard manual segmentation of the bone surfacesin a series of US images by a number of expert sonographers. However,for the ex vivo data I provide an indirect quantitative comparison (insteadof direct comparison using manual segmentation) in the next chapter byevaluating our statistical model registration to the extracted bone surfacewith gold-standard CT. As the registration accuracy is directly dependenton the extracted bone surface by a bone enhancement method, a compar-ison of the registration accuracies in turn indicates a comparison amongthe bone enhancement methods. Unfortunately, for the phantom and invivo data we cannot compute the registration accuracy due to unavailabilityof corresponding gold-standard CTs. The details about the calculation ofregistration accuracies will be discussed in the next chapter.4.4 ResultsFigs. 4.6-4.8 demonstrate the bone enhancement of my 2D and 3D techniquesfor the phantom, ex vivo and in vivo data, respectively. Comparative resultsby the 2DLPS [37] and 3DLPS [40] methods are also shown in these figures.As mentioned earlier, the bone appearance in the phantom is bettercompared to those in the ex vivo and in vivo data, as a result, a better boneenhancement is expected for the phantom data for all of the comparingmethods. A volume rendering of the bone enhancement, therefore, providesa clear representation of the wrist bone surfaces (Fig. 4.6). The comparativeresults for the phantom do not indicate significant noticeable differencesbetween the competing methods except at the pisiform bone locations. In594.4. Results(a) Our 2D (b) Our 3D(c) 2DLPS (d) 3DLPSPisiformFigure 4.6: Comparative bone enhancement results in the phantom experi-ment.addition, the 2D-based methods produce more amount of outliers comparedto the 3D-based methods.In contrast, for the ex vivo and in vivo data, my methods demon-strate better performance compared to the corresponding state-of-the-artapproaches (Figs. 4.7-4.8). Note that the 2D methods are applied across theaxial slices to obtain the bone enhancement of the given US volume, there-fore, a better bone enhancement is expected across the axial slices thanacross the other directions.604.4. Results2DLPS3DLPSOur 2DUS frameOur 3D(a)(d)(c)(b)(e)(f)(i)(h)(g)(j)Axial SagittalFigure 4.7: Results of the proposed, 2DLPS and 3DLPS methods for theex vivo experiment. Example axial and sagittal US frames are shown in(a), (f). The corresponding bone enhancements are demonstrated below.The differences in the enhancement are prominent in the surfaces markedby arrows.The recorded run-times are 4, 11, 3 and 10 seconds for the proposed 2D,proposed 3D, 2DLPS and 3DLPS methods, respectively.614.4. Results2DLPS3DLPSOur 2DUS frameOur 3DAxial Sagittal(a)(d)(c)(b)(e)(f)(i)(h)(g)(j)Figure 4.8: Results of the proposed, 2DLPS and 3DLPS methods for thein vivo experiment. Example axial and sagittal US frames are shown in(a), (f). The corresponding bone enhancements are demonstrated below.The differences in the enhancement are prominent in the surfaces markedby arrows.624.5. Discussion4.5 DiscussionA comparison of four bone enhancement approaches indicates the advantagesof the proposed methods (especially the 3D one) in detection of the curvedbone edges and low contrast bone surfaces compared to two state-of-the-artapproaches. I also noticed that both 2D-based methods have demonstratedworse performance along the sagittal directions.In phantom construction, the soft-tissue has been modelled using an uni-form chemical substance, therefore, there were minimal changes of acousticimpedances in the soft-tissue region. As a result, I have not noticed anyhigher intensities in the soft-tissue region. However, in the in vivo and exvivo data, I obtained higher intensities in the soft-tissue regions as the acous-tic impedances are not uniform within the soft-tissue regions, that mainlylead to less accuracies in the in vivo and ex vivo experiments.The volume rendering of the bone enhancement in Fig. 4.6 demonstratedthe partial wrist bone surfaces available in the US data, i.e., only a portion(around 30%) of true wrist bone surfaces is available in the US data. Suchsmall amount of available information in US data also indicates the impor-tance of a better bone enhancement technique as a less accurate methodfurther limits the information.One of the drawbacks of my method is the segmentation of the localfrequency responses (e.g., p(φ), qo(ω)). It may detect false minima in theresponse or if two maxima are too close, then it may detect them as one singlemaximum. As a result, the local spectrum will not be accurately segmentedand subsequent local PS estimation will be less accurate. Another limitationis enhancement of the symmetrical noise. This type of noise mainly appearsas scattered objects (marked by circles in Figs. 4.7-4.8 in the bone enhancedvolumes. In addition, the shadow map used in this work was not able toreduce the non-bony responses substantially.4.6 SummaryI have presented bone enhancement methods in US based on the utilizationof local spectrum variation. The introduction of the spectrum variationin the filter design allows us to estimate the local phase symmetry moreaccurately, subsequently better enhancing the expected bone locations. Ihave applied my technique to 2 wrist phantoms, 13 cadaver wrists, and 14healthy volunteer wrists, the comparative results demonstrate the improvedperformance of my detection compared to two state-of-the-art approaches634.6. Summaryespecially for the ex vivo and in vivo data.64Chapter 5Scaphoid screw axisestimation based on USguidance5.1 Introduction and backgroundAcute non-displaced scaphoid fracture is conventionally fixed using a percu-taneous surgical procedure by placing a surgical screw along the longitudinalaxis of the scaphoid bone. The main objective of the screw insertion is toallow the fractured parts to move as a rigid-body. To achieve improvedbiomechanical and clinical results of the recovered wrist, it is advocated toinsert the screw along the central axis of the scaphoid bone, use a longerscrew length [22] and place the screw in the central one-third of the proxi-mal pole [92]. In addition, the optimal screw axis has to be volar accessible,i.e., the elongated screw should not be obstructed by the trapezium bone.In fact, the anatomical position of the trapezium limits the potential entrypoints of the screw to a narrow window [55, 61]. To overcome this problem,alternatives were suggested by either keeping the wrist at ulnar deviationand extension position during the surgery [3] or partial resection of thetrapezium [56] or placement of the screw through the trapezium [55, 60].In order to estimate the optimal screw placement a number works havebeen published in the literature. Menapace et al. [63] suggested to place thescrew in a continuous column along the longitudinal axis of the scaphoidbone that keeps intact the vascular supply and the articular surfaces ofthe scaphoid bone. They defined a safe zone for the entry point of theguide wire based on the bone vascularity and morphology of the scaphoidbone, however, optimal screw length and placement of the screw within theThis chapter is adapted from Anas, E. M. A., Seitel A., Rasoulian A., John P. St.,Ungi T., Lasso A., Darras K., Wilson D., Lessoway V., Fichtinger G., Zec M., Pichora D.,Mousavi P., Rohling R. N., and Abolmaesumi P., Registration of a statistical model tointraoperative ultrasound for scaphoid screw fixation, International Journal of ComputerAssisted Radiology and Surgery, 1-9, 2016.655.1. Introduction and backgroundsafe zone were not investigated [55]. Levitz and Ring [56] estimated theparameters of the position of a 3.0-mm screw in the scaphoid bone usingCT images, however, they restricted themselves to the coronal, sagittal andaxial planes of the scanner, leads to a limited number of solutions [55].Recently, Leventhal et al. [55] investigated the optimal screw placement inthe scaphoid bone using the 3D digital model of the scaphoid bone. Theydefined a safe zone based on the thickness of the inner cortical surface ofthe scaphoid bone and the radius of the cannulated screw. In addition,they proposed two methods for estimation of the optimal screw axis: (1)maximum screw length (MSL) within the safe zone and (2) a cylinder best-fit (CYL) to the safe zone. Finally, the obtained screw axes were evaluatedin regard to the volar accessibility, screw length, and location of the screwaxis.The above mentioned works are based on the assumption of availabilityof gold-standard patient specific 3D model of the bony anatomies duringthe surgical navigation, which is practically almost impossible. Current sur-gical procedures use fluoroscopy as the intraoperative imaging modality toderive an estimate of the patient specific 3D model, however, fluoroscopyprovides only projection view of the 3D anatomy. US has been suggested asan alternative modality for a real-time 3D guidance, unfortunately, it pro-vides only partial view of the bony anatomies. To improve interpretationof US images, it has been proposed to fuse this data with anatomical infor-mation from computed tomography (CT) [7]. Beek et al. [7] explored thefeasibility of accurate screw insertions by registering manually segmentedsurfaces of the scaphoid bone from preoperative wrist CT to intraoperativeUS images. The main drawbacks of their method include the challengingregistration given the smooth featureless surface of the scaphoid bone andthe mandatory, labour intensive manual bone segmentation.In this chapter, I propose to fuse a multi-object statisticalshape+scale+pose model with the preoperative CT and the intraoperativeUS to more accurately estimate patient specific 3D model of the wrist bones.The fusion consists of two steps, the first step of model registration to CTwas described in Chapter 3. The second step of registration of the model toUS is presented in this chapter. As only the wrist pose during CT and USacquisitions are different, the shape+scale information estimated from CTsegmentation (Chapter 3) are used in this chapter. Therefore, the statisticalmodel registration to US is formulated by optimizing only the pose modes.The target point cloud corresponding to US is generated from the boneenhancement (described in Chapter 4) of the US volume. The registeredwrist bones are finally used to estimate the optimal trajectory for scaphoid665.2. Methodsscrew insertion. I determine the accuracy of the scaphoid screw axes in 13cadaver experiments.5.2 MethodsEstimation of the scaphoid screw axis is accomplished by registering thewrist model to an intra-operatively acquired US volume. Initially the modelis aligned to the bone surfaces extracted from the bone enhanced volume(Sec. 5.2.1), next, the pose coefficients of the model are optimized in Sec. 5.2.2.The registered model is then used to estimate the screw axis (Sec. 5.2.3).5.2.1 Initial alignmentGiven an US volume, the bone enhancement of the US volume is per-formed (the 3D bone enhancement presented in the previous chapter). Thebone enhanced volume is then thresholded using the empirically determinedthreshold value of 15% of the maximum intensity value of the enhanced vol-ume to extract the wrist bone surfaces. Note that a small variation in thementioned threshold value can be compensated by the probabilistic-basedregistration technique (Section 5.2.2) that is robust to outliers and miss-ing points. For initialization of the model in the US volume, I have useda landmark-based approach where at least five anatomical correspondinglandmarks (e.g., Fig. 5.1) are chosen manually across the model and theextracted bone surface.5.2.2 Pose optimizationAfter initial alignment, the model is registered to the extracted bone sur-faces using the probabilistic-based registration presented in chapter 2. Notethat in this registration step I only optimize the pose coefficients; the shapeand scale coefficients are re-used from the previous CT registration. Theregistered model to US for the l-th wrist bone can be expressed asYl,US(ΘsCT ,ΘkCT ,Sl,CT ,ΘpUS ,Θp′US) = (5.1){φkl (ΘkCT )Φsl (ΘsCT ) + Sl,CT } × {Rp′l (ΘpUS ,Θp′US)}T+IM × {tp′l (ΘpUS ,Θp′US)}T ,Here, ΘpUS and Θp′US represent the inter-subject and intra-subject modelcoefficients, respectively.675.2. MethodsFigure 5.1: An example of manual landmark selection across the wrist modelinstance (a) and the extracted bone surface in US (b).5.2.3 Scaphoid screw axis estimationTo estimate the optimal screw axis of the scaphoid bone for the fracturefixation, I have used the maximum screw length (MSL) technique developedby Leventhal et al. [55]. In the MSL method, a safe zone is defined toensure a safe screw placement without penetrating the outer cortical bonesurface. Assume the thickness of the cortical bone surface is xc mm, thescrew radius is xr mm and an additional safety offset is xs mm. Then,the safe zone is estimated by computing the surface located (xc + xr +xs) mm inside the original scaphoid bone surface. In this work, I considerxc = 0.35, xs = 0.25 [55] and xr = 1.50. xr is taken as the maximum oftypical screw radius values (1.20-1.50 mm). The safety offset xs is introducedhere to account for the registration inaccuracy. A registration techniquewith greater registration error should uses a greater offset to avoid anyintersection of the screw with the outer cortical bone surface. However, itleads to a safe zone with smaller area, in turns a reduced screw length; wherea greater screw length is clinically suggested for a better recovered wrist.The screw axis for scaphoid fracture fixation is then determined as thelongest possible axis within the safe zone considering the combination ofeach pair point set in the inner safe zone surface. Note that the screw axiscomputation is completely independent from the optimization/registration685.3. Experiments and evaluationof the model to the US or CT. In addition to the screw axis computation, Ielongate the screw to the distal part of the scaphoid bone to check the volaraccessibility, i.e., whether the elongated screw intersects the trapezium boneor not. If the screw is not volar-accessible, then the estimated screw canbe still drilled through the scaphoid bone either by partial resections of thetrapezium or by placing the screw through the trapezium [55, 60].5.3 Experiments and evaluationFor evaluation of the proposed approach, I performed a cadaver study of13 cadaver wrists. Note that these cadavers were also used in the previouschapter for bone enhancement experiment. In addition to acquire the USimages from each cadaver wrist, the US transducer was tracked with anoptical tracking system (Polaris Spectra, Northern Digital Inc (NDI), Wa-terloo, ON, Canada). For ground-truth generation from the wrist at theextension position and for preoperative wrist CT at the neutral position,CT volumes were acquired for each cadaver wrist using the GE Lightspeedequipment in Kingston General Hospital (KGH) Imaging Services depart-ment. Six fiducial markers were attached to the wrist holder in order toobtain a transformation between the ground-truth CT and US coordinateframes.The reference CT volume was manually segmented using the interactivesegmentation framework of the Medical Imaging Interaction Toolkit [59] toobtain a ground-truth surface representation of the wrist bones as well asthe fiducial positions on the wrist holder. Using the segmented fiducialpositions in CT and the recorded fiducial positions with respect to the ref-erence frame, the segmented bone surfaces in CT were transformed to theUS coordinate. Due to the US transducer’s pressure on the wrist, the wristbones were mostly moved down during the US scanning as there was nosupport below the wrist joint. Therefore, the transformed segmented CTbones superimposed on the US volume appear above compared to the ac-tual position. An example US slice along with the transformed segmentedCT bones is demonstrated in Fig. 5.2(a). This error, which was mostly alongthe direction of the US scanning axis, is minimized by manually registeringthe manually segmented CT bone surfaces to the US bone responses (shownin Fig. 5.2(a)). The manual correction is performed using 3DSlicer [24] bytranslating the volume rendered CT such that it visually fits the bone sur-faces in the US volume. I start the correction with a translation along the USscanning direction, followed by fine translations along the other directions.695.3. Experiments and evaluation(a)Before correctionAfter correctionScaphoidTrapeziumScrew lengthScaphoid lengthVolar accessibilityCentral one-third of proximal poleMinimum distance between scaphoid and screw surfaces(b)Figure 5.2: Evaluation of the proposed approach. (a) Effect of the US trans-ducer’s pressure on the fiducial-based registration of manually segmented CTto US. The registered segmented CT (marked by white) is overlaid on an USimage, and the registered surface does not align with the bone responses inUS. To minimize the effect of the US transducer’s pressure, I manually cor-rect the registration and the manually registered segmented CT (marked byyellow) corresponds better to the US bone responses. (b) Several evaluationcriteria are measured from the estimated screw position and the referencewrist bone surfaces.705.3. Experiments and evaluationUsing the generated ground truth, I evaluated my proposed method ac-cording to the following parameters. Fig. 5.2(b) shows a schematic diagramfor calculation of the parameters from the estimated screw position overlaidon the generated ground truth.i) Mean surface distance error: The surface distance error (SDE) isdefined for each surface point (of the registered scaphoid bone) by comput-ing its Euclidean distance to the closest neighboring point in the referencesurface. The mean surface distance error (mSDE) between the registeredand the reference scaphoid surfaces is the average of SDEs across all surfacepoints. As mentioned in the previous chapter, we also report the mSDEsusing other state-of-the-art bone enhancement techniques (2DLPS [37] and3DLPS [40] methods). To use other methods, the whole proposed work-flowremains same except the bone enhancement is performed by each of thecomparing methods.ii) Volar accessibility: Volar accessibility is achieved if the calculatedscrew axis does not intersect the neighboring trapezium bone in neither themodel nor the ground truth CT space (Fig. 5.2(b)).iii) Screw length: The screw length (Lscrew) estimation is evaluated bycalculating its ratio (Lratio) to the maximally allowed screw length definedin the CT reference as the scaphoid length (Lscap) minus 4 mm [22]. Lscrewis defined as the length of the screw in mm along its main axis and Lscap ismeasured as the length of the longest axis in mm within the scaphoid boneof the reference surface (Fig. 5.2(b)).iv) Cortical surface intersection: I calculate the minimum distance(Fig. 5.2(c)) between the scaphoid and cylindrical screw surfaces, and com-pare this distance with the measured standard thickness of the cortical sur-face (0.35 mm). If any intersection occurs, then I estimate the maximumdistance between the intersection areas and the number is reported as anegative number.v) Proximal position: I check the relative location of the screw withrespect to the central one-third of the proximal pole of the scaphoid bone(Fig. 5.2(c)). According to Trumble’s criteria [92], part of the screw shouldbe located in the central one-third of the proximal pole for improved ratesof healing.715.4. Resultsvi) Computation time: To evaluate the computational complexity ofthe proposed segmentation method, I measure the overall computation timeof the algorithm as well as its different processing steps. The computationtime is recorded from unoptimized MATLABTM (Mathworks, Natick, MA,USA) code running on an Intel Core i7-2600M CPU at 3.40 GHz.5.4 ResultsThe results obtained from 13 cadaver wrists based on the above mentionedcriterion are presented below:i) Mean surface distance error: Fig. 5.3 demonstrates an example reg-istered atlas superimposed on the US volume for our method. Table 5.1presents mSDEs for our method in mm for all 13 cadaver scaphoid bones.The mSDE averages to 0.8 ± 0.2 mm. The corresponding obtained mS-DEs for the 2DLPS and 3DLPS methods are 1.1 ± 0.3 mm and 0.9 ± 0.2mm, respectively. The reported Hausdorff distance errors are 1.9 ± 0.3 mm,2.7 ± 0.4 mm and 2.4 ± 0.4 mm for our, 2DLPS and 3DLPS methods,respectively.ii) Volar accessibility: The results of the volar accessibility for both ofthe registered model and the reference surface are demonstrated in Table 5.1.92% successful volar accessibility was obtained in the cadaver wrists usingmy model-estimated drill axes. An example case is shown in Fig. 5.3(b).iii) Screw length: The screw lengths estimated using my method for 13cadaver wrists are summarized in Table 5.1. The screw length was estimatedas 23.8 ± 2.6 mm, corresponding scaphoid length: 29.5 ± 2.8 mm meaningthat the estimated lengths are in between 90% and 98% of the maximumallowable screw lengths (average 94%).iv) Cortical surface intersection: The minimum distance between thescaphoid and cylindrical screw surfaces for each cadaver averaged to 0.5 mm(Table 5.1). Out of 13 cases, one case has value less than 0.35 mm (standardcortical surface thickness) and another has a negative value indicating anintersection between the screw and the scaphoid surfaces.725.5. DiscussionTable 5.1: Results of the proposed approach with respect to some definedcriterion. ‘MD’, ‘VA’, ‘Y’ and ‘N’ indicate minimum distance, volar acces-sibility, yes and no answers, respectively.Case mSDE VA by Lscrew Lratio MD Prox.(mm) my ref (mm) (mm) pos1 0.6 Y Y 24.0 93.2 0.6 Y2 0.8 Y Y 24.9 95.1 0.4 Y3 0.6 N N 24.2 93.5 1.2 Y4 0.7 Y N 20.1 92.1 0.6 Y5 1.1 Y Y 21.7 93.9 0.2 Y6 0.8 Y Y 27.4 92.9 0.5 Y7 0.8 Y Y 27.4 93.0 0.6 Y8 0.6 Y Y 20.3 98.2 0.7 Y9 1.2 N N 19.8 90.9 -0.9 Y10 0.6 Y Y 25.2 94.7 0.6 Y11 0.5 Y Y 24.4 93.0 1.1 Y12 0.7 Y Y 24.8 92.8 0.7 Y13 0.8 Y Y 25.2 93.5 0.7 YAvg. 0.8 92.3% 23.8 93.6 0.5 100%v) Proximal position: Table 5.1 shows the results of the relative lo-cations of the screw positions with respect to the central one-third of theproximal pole and my estimated screws have 100% accuracy in locating thescaphoid screw in the central one-third of the proximal pole.vi) Computation time: The average computation time is recorded as1:34 min for the intra-operative 3D US registration step.5.5 DiscussionWe have applied our technique to 13 cadaver wrists, and obtained an averagemSDE of 0.8 mm and an average mxSDE of 1.9 mm between the registeredand reference scaphoid bone surfaces. Though our mxSDE improvementof 0.5 mm compared to the state-of-the-art approach (3DLPS) is small inabsolute magnitude, the achieved improvement is significant at about 25%of the clinical surgical accuracy (2 mm). The rank-sum test comparing themxSDEs of our and 3DLPS methods rejected the null hypothesis of equal735.5. Discussion(a)ScaphoidTrapezium(c) US scanning    direction(b)(d)Figure 5.3: Example results of the proposed approach. (a) Example reg-istration result of the statistical model to US. (b) Screw trajectory esti-mated by my method. The trajectory is not obstructed by the surroundingtrapezium bone. (c) An example of a registered model overlaid on an USframe. I notice a misalignment between the registered model and the boneresponse marked by a dashed arrow. The source of this error mainly in-cludes the inaccurate segmentation in CT images. (d) Segmentation resultof corresponding preoperative CT. The mentioned inaccurate segmentationis marked by a dashed arrow.medians with p=1.7×10−5.Although the results in this feasibility study are promising, in three outof 13 cases the proposed approach fails to meet the clinical guidelines. In twofailed cases (case no. 5 and 9), my estimated screw surfaces intersected thecortical scaphoid surfaces. In another case (case no. 4), my method couldnot correctly determine the volar accessibility due to incorrect registration745.5. Discussionto the trapezium bone. The above mentioned violations mainly occur dueto the registration inaccuracies and the major source of error is the absenceof surrounding bone responses in US. For example, in case no. 5 and 9,the field of views of the US scans did not contain the radius and scaphoidbones, respectively, and the estimated screw in those cases intersected thecortical surface of the scaphoid bone. An example of improper field-of-view is demonstrated in Fig. 5.4(a) (case no. 9). For the case no. 4, aminimal bone response was observed for the trapezium bone (Fig. 5.4(b)).A less registration accuracy has been obtained for this bone that leads tothe failure of the volar accessibility test. Therefore, one possible way toimprove the accuracy is to ensure inclusion of all of the four surroundingwrist bones (part of radius, trapezium, lunate and part of capitate) in the3D US acquisition. Another way to improve the accuracy is to consider agreater safety offset xs. However, a greater offset reduces the area enclosedby the safe zone that leads to a reduced screw length.(a) (b)part of scaphoid     is missingminimal bone responses        for trapeziumFigure 5.4: Possible sources of errors. (a) A limited field of view for thescaphoid bone. (b) Minimal bone response is observed for the trapeziumbone. The minimal bone responses mainly occur due to non-perpendicularscanning of the trapezium bone.Although the registration errors in CT registration step are quite small(mean error of around 0.3 mm), those registration inaccuracies propagate tothe US registration result as the shape and scale coefficients are not updatedin the US registration step. Such a case is demonstrated in Fig. 5.3(c),where the registered model does not align well with the corresponding USbone responses. This misalignment mainly occurs due to the inaccuratesegmentation in CT images (Fig. 5.3(d)). A manual correction (at least forthe scaphoid and trapezium bones) after the pre-operative CT registrationstep may improve the overall accuracy.755.6. SummaryThe time for intra-operative registration is currently not yet suited fora real-time scenario. However, in a related research project for US spineregistration [15], I have demonstrated an efficient C++ implementation us-ing multi-core processing that has the potential to improve run time to anorder of seconds. Especially, an efficient implementation of the calculationof the time consuming matrix exponential, used several times during modeloptimization, is expected to significantly decrease the current run time.Future work should include further evaluation of this technology withinthe clinical workflow in a prospective patient study. The major challengesin the clinical application may include the tracking errors of drilling instru-ments and unwanted movements of the wrist bones during drilling. Fur-thermore, the post-operative confirmation of the screw position inside thescaphoid bone is also a challenge in my ultrasound-based guidance. A post-operative fluoroscopy may be used for the confirmation. Future work alsoincludes introducing an improved screw axis estimation by jointly optimiz-ing most important clinical guidelines. In addition, to minimize the bonemovements due to the US transducer’s pressure during US acquisition, abetter US wrist imaging setup can be used. Moreover, an automatic initial-ization of the wrist model to the US volume is planned to be integrated intothe model registration workflow.5.6 SummaryI have proposed a model-based approach to estimate the screw insertion axisfor a volar percutaneous scaphoid fracture fixation using 3D US volumes.A statistical shape+scale+pose model was used that allows incorporationof pre-operative CT data in the shape and scale estimation process, andoptimizes the model pose with respect to the information available from theintra-operative 3D US scan. I have applied the proposed approach to 13cadaver wrists and achieved a mSDE 0.8 ± 0.2 mm for the scaphoid bone.An average 94% of maximum allowable screw length is obtained based onmeasurements from gold standard CT. Also, I obtained an average 92%successful volar accessibility. In 10 out of 13 cases, the estimated screwtrajectory meets all important clinical criteria. In conclusion, the developedapproach has demonstrated promising results in scaphoid screw positionestimation using 3D US-based guidance.76Chapter 6Joint atlas-based registrationof multiple bones in CT andUS6.1 Introduction and backgroundUltrasound has become popular as an intraoperative imaging modality incomputer-assisted orthopaedic surgical procedures due to its certain advan-tages including non-ionizing nature, real-time 3D acquisition capability, easyaccessibility and low cost over the fluoroscopic navigation. US images are,however, difficult to interpret as they are often corrupted by significantamount of noise or artifact, and bone responses can appear blurry and dis-connected. To improve interpretation of US images, it has been proposed tofuse this data with anatomical information from CT [7, 102] or a statisticalmodel [81].Prior works on fusion with US mostly focus on registration to either CTor the statistical model. In the previous chapter, a sequential registrationof preoperative CT to intraoperative US via a statistical model has beenpresented to address the individual limitations of CT-US and statisticalmodel-US fusions. In short, a statistical anatomical model was initially reg-istered to preoperative CT to obtain the shape and size, and subsequentlyto intraoperative US to derive the pose of multiple bones. As the shapeof the bones were optimized only in the preoperative registration using CTinformation, that approach did not take advantage of additional shape in-formation in readily available US images; further, any error in the shape ofthe registered bones was propagated to the final registration results.To avoid the aforementioned problems in the sequential registration, ajoint registration of a statistical shape+pose spine model to preoperative CTand intraoperative US images was presented recently [8] to simultaneouslyalign the information from the two datasets (CT and US). That method776.2. Methodsshowed significant promise in improving the overall alignment of the modelto both imaging modalities; however, an outstanding challenge was that posestatistics of the bones also captured the modes of variation that account fortheir scale. As a result, optimization of the pose also altered the model’sscale, which was not desirable.In this chapter, instead of using a statistical shape+pose model, a sta-tistical shape+scale+pose model is used to enable a joint fusion to CT andUS data of multiple bone structures. The joint optimization of the shapeand scale coefficients allows for flexibility to modify each parameter inde-pendently. The patient position differences in the acquired CT and USvolumes are handled by independent pose optimizations to CT and US. Ievaluate the proposed registration framework on a wrist dataset consistingof 13 pairs of CT and US images of wrist bones. In addition, the registra-tion accuracy of this technique is compared with those obtained using thesequential registration approach.6.2 MethodsDetails of my proposed joint registration work-flow is shown in Fig. 6.1. Ioptimize the shape+scale coefficients using both CT and US, and indepen-dently optimize the pose coefficients to register to CT and US. The steps ofthe proposed approach are described below:CT aquisitionUS aquisitionTarget generation & initializationBone enhancement & initializationSoft correspondence assignmentShape+scale joint optimizationCT pose optimizationUS pose optimizationCT volumeUS volumeUS target point cloudCT target point cloudModel instance registered to CTModel instance registered to USBoth converged?Statistical shape+scale+pose modelModel instantiation & transformationModel instantiation & transformation(1b)(2)(4a)(4b)(3)(5a)(5b)NoIterave framework(1a)Figure 6.1: A work-flow of the proposed joint registration approach.786.2. Methods1. Target generation and initialization: The joint registration startswith the mean model, i.e. values of zero for modes of variations ofshape, scale and pose. The number of assumed variation modes forshape, scale and pose are selected such that they capture 95% of thetraining samples.(a) A point cloud is obtained from the CT volume by applying a thresh-old. The initialization of the model instance to CT is performed au-tomatically by aligning the principal axes between the target and themean model instance.(b) The bone responses in a 3D US volume are enhanced using theUS bone enhancement technique presented in Chapter 4. For initial-ization, a manual landmark-based technique is followed to align themean model instance to US.2. Soft correspondence assignment: An Expectation-Maximization(EM) framework is followed in the model registration. The key ideawas described in Chapter 2. The expectation step consists of compu-tation of the soft correspondences PI(ylm|zn,I) between the m-th modelinstance point of l-th bone (ylm) and n-th target point of I-th imagingmodality (zn,I), where I ∈ {CT,US}.3. Shape+scale joint optimization: The joint objective function opti-mized in the maximization step can be written as: Q =∑I∈{CT,US}QI,where the cost function for I-th imaging modality is defined as:QI =L∑l=1M,NI∑m,n=1PI(ylm|zn,I)||zn,I − yml (Θs,Θk,ΘpI ,Θp′I )||2, (6.1)where Θs, Θk, Θp and Θp′represent the shape, scale, inter-subjectpose and intra-subject pose coefficients, respectively. M , NI and Lrepresent the total number of points in the model per bone, totalnumber of points in the I-th target point cloud and the total numberof bones in the atlas, respectively. A Quasi-Newton method is followedto optimize the joint cost function (Q) with respect to the shape, scaleand inter-subject pose coefficients.4. Pose optimization: The modality-specific pose coefficients are opti-mized using the individual cost function (QI).5. Model instantiation: After estimation of the shape, scale and posecoefficients, I update the model instance, and repeat steps (2-5) untilthe iterative EM framework reaches convergence.796.3. Evaluation6.3 EvaluationThe cadaver experiment in the previous chapter is used here to evaluatethe proposed joint registration technique. In short, a study consisting of 13cadavers was performed to acquire 3D US at an extension position focusingmainly on the scaphoid and its surrounding bones. In addition, two CTs atextension and neutral wrist positions for each cadaver wrist were acquired.The first CT provided the ground truth at the extension wrist position,while the latter was considered as preoperative data. An optical tracker andsix fiducial markers were used for tracking the US transducer.To generate the ground truth for US registration, the wrist CT at ex-tension position was segmented manually using Medical Imaging InteractionToolkit (MITK). Subsequent fiducial-based registration along with manualcorrection were used to transform the segmented CT to US coordinates. Thestudy was approved by an institutional research ethics board.To evaluate the proposed approach, we followed the same evaluatingcriteria defined in the previous chapter. In addition, I compare my jointregistration approach with the sequential registration technique describedin the previous chapter.6.4 ResultsThe results obtained from 13 cadaver wrists based on the previously de-scribed criterion are presented below:i) Mean surface distance error: Table 6.1 presents mSDEs in mm forall 13 cadaver scaphoid bones. Using my proposed joint registration ap-proach, the mSDE of 0.7 ± 0.2 mm between the registered and ground-truthscaphoid bone surfaces was obtained, while the sequential registration-basedmethod leads to an mSDE of 0.8 ± 0.2 mm. The corresponding mxSDEs are1.7 ± 0.3 mm and 1.9 ± 0.3 mm, respectively. Hence, a good improvementof mxSDE has been achieved.Fig. 6.2(a-b) demonstrates an example of improved results using the jointregistration vs. the sequential registration methods. The comparative re-sults indicate that the sequential registration cannot correct the errors dueto the CT image artifacts. For example, the shape error in Fig. 6.2(b) resultsfrom the partial volume artifact in CT registration (Fig. 6.2(c)). As my jointregistration technique is based on using complementary information in CTand US for the shape+scale estimation, it can overcome the effect of the CT806.4. ResultsJoint registration Sequential registration(b)(a) (c)Figure 6.2: Results of the proposed joint registration approach and com-parison with the sequential registration technique for the scaphoid bone.(a-b) Example US registration results for the proposed joint and sequentialregistration methods, respectively. (c) Corresponding CT slice to show thesource of error of the sequential registration method.image artifacts by utilizing the corresponding bone responses in US.Table 6.1: Results of the proposed approach with respect to some definedcriterion. ‘MD’, ‘VA’, ‘Y’ and ‘N’ indicate minimum distance, volar acces-sibility, yes and no answers, respectively.Case mSDE VA by Lscrew Lratio MD Prox.(mm) my ref (mm) (mm) pos1 0.6 Y Y 24.2 93.8 0.6 Y2 0.7 Y Y 24.7 94.3 0.5 Y3 0.6 N N 24.5 94.7 1.1 Y4 0.7 Y N 20.2 92.6 0.7 Y5 1.1 Y Y 22.0 95.2 0.2 Y6 0.7 Y Y 27.5 93.2 0.6 Y7 0.7 Y Y 27.2 92.3 0.6 Y8 0.6 Y Y 20.1 97.2 0.6 Y9 1.1 N N 19.9 91.4 -0.8 Y10 0.6 Y Y 25.4 95.5 0.7 Y11 0.5 Y Y 24.8 94.5 1.3 Y12 0.5 Y Y 24.7 92.4 0.7 Y13 0.7 Y Y 25.3 93.9 0.6 YAvg. 0.7 92.3% 23.9 93.9 0.6 100%816.5. Discussionii) Volar accessibility: 92% successful volar accessibility was obtained inthe cadaver wrists using the joint registration technique, which is comparableto the sequential registration result.iii) Screw length: The screw lengths estimated using my method for 13cadaver wrists are summarized in Table 6.1. The screw length was estimatedas 23.9 ± 2.6 mm, and the estimated length ratio to the maximum possiblelength has an average of 94%. These results are also comparable to thesequential registration result.iv) Cortical surface intersection: The minimum distance between thescaphoid and cylindrical screw surfaces for each cadaver averaged to 0.6 mm(Table 6.1). Out of 13, two cases have less than the average cortical surfacethickness (0.35 mm).v) Proximal position: Table 6.1 shows 100% accuracy in locating thescaphoid screw in the central one-third of the proximal pole.vi) Computation time: The computation time is recorded as 9:16 min.This computation is higher than the seuqnetial registration.6.5 DiscussionWe observe an improvement in the registration accuracies (mSDE and mxSDE)for the joint registration technique compared to the sequential registrationapproach. An example graphical comparison has been demonstrated inFig. 6.2, where a preoperative shape registration error due to partial volumeeffect (Fig. 6.2(b)) propagates to the US registration result. In contrast, ourjoint registration technique overcomes the effect of weaker gradients in theCT images by utilizing the corresponding bone responses in US (Fig. 6.2(a)).For clinical evaluations, unfortunately any improvement has not beenobserved for the joint registration approach. Though these clinical evalua-tions do not indicate improvement of the joint registration technique overthe sequential registration approach within this limited number of test ex-amples, we believe that in a large number of test examples we may noticethe clinical improvement for the joint registration method.The main limitation of the proposed joint-registration method is its highcomputational complexity. Though tolerable for preoperative CT registra-tion, the computational complexity of the proposed method is still too high826.6. Summaryto allow real-time implementation for intraoperative procedures. One wayto improve the computational complexity is to perform registration only toCT preoperatively, followed by starting the intraoperative joint registrationapproach with the preoperatively optimized model coefficients, instead ofthe mean shape, for faster convergence.6.6 SummaryI presented a joint registration framework that takes advantage of a sta-tistical model to align CT and US data of multi-bone structures. My keycontribution involved the concurrent optimization of the shape and scalecoefficients, which lead to accurate fusion of the two modalities despite thesubstantial pose differences between them. The cadaver experiment hasdemonstrated improved performance of the joint registration technique overthe sequential registration method.83Chapter 7Conclusion and future workUltrasound has been emerged as a potential intraoperative imaging modalityin computer-assisted orthopaedic surgical applications. However, the limitedview of the bony anatomies and corruption of noise/artifacts in the ultra-sound images make interpretation of the ultrasound data difficult. An at-tempt to address this problem in percutaenous non-displaced scaphoid frac-ture fixation is presented in this thesis by augmenting the ultrasound datawith anatomical information from CT and statistical anatomical model. Thestatistical model is generated by capturing the shape, scale and pose varia-tions of the wrist bones across population at wide range of wrist positions.An improved registration of the model to ultrasound is made possible byenhancing the bone responses in ultrasound. The registered model is finallyused for estimation of the optimal scaphoid screw axis to guide the surgicalprocedure. This study demonstrates the potential of the proposed techniqueto be included in a clinical ultrasound-based percutaenous scaphoid fracturefixation procedure.In this thesis, novel techniques were presented for fusion of a statistical modelwith preoperative CT and intraoperative ultrasound images. In Chapter 2a multi-object statistical model is developed to capture the main modes ofshape, scale and pose variations across population at wide range of wristpositions. A registration of this model to CT images is presented in Chap-ter 3 for an automatic segmentation of the wrist bones in CT images. Toimprove the interpretation of the ultrasound, bone enhancement techniquesare proposed in Chapter 4 by utilizing the local spectrum variation of theultrsound images. The optimal scaphoid screw axis is estimated based onthe model registration to intraoperative ultrasound in Chapter 5. Finally, asimultaneous optimization of the model shape and scale coefficients is pre-sented in Chapter 6 based on the CT and US data as an alternate way toderive the patient-specific 3D models of the wrist bones.7.1 ContributionsThe contributions of this thesis are summarized as follows:847.1. Contributions• A novel statistical shape+scale+pose model is developed to capturethe main modes of shape, scale and pose variations across populationat wide range of wrist positions. A group-wise registration frameworkis used for an improved correspondence establishment across the train-ing shapes. The pose model consists of two models, inter-subject andintra-subject pose models, representing the wrist pose variations acrosspopulation and wrist positions, respectively. The presented model alsoseparates the scale statistics from the pose one, allowing an improvedfusion of the model to CT and US. I demonstrate that the multi-objectstatistical wrist model is capable to register to target wrists at differentwrist positions.• A novel technique is presented for an automatic segmentation of thewrist bones in CT images. Segmentation of the wrist bones is chal-lenging due to blurry narrow inter-bone spaces of the wrist bones. Totackle such challenges, the proposed segmentation technique allowsthe model to register to the target CT by optimizing the model co-efficients. In addition, a non-rigid registration is used at the end tocorrect the shape residuals that cannot be captured by the proposedmodel. The presented technique is demonstrated to segment the wristbones in healthy and unhealthy subjects.• A novel bone enhancement technique is presented based on local spec-trum variation of the US images to design a set of Log-Gabor filters forimproved bone enhancement in US images. The technique to utilizethe local spectrum variation is presented for a 2D US image, followedby extension of the technique to 3D US volume. Shadow informationis also used for enhancing the expected bone locations rather than thesoft-tissue interfaces. The presented techniques are capable to enhancethe low-contrast and curvy bone surfaces in the US images.• A novel technique is developed for estimation of the scaphoid screwaxis based on intraoperative US to guide a scaphoid fracture fixation.The estimated screw trajectory ensures placement of the screw alongthe longest axis of the scaphoid bone as well as within the central onethird of the proximal pole of scaphoid bone.• A framework is developed for joint registration of a statistical wristmodel to CT and US with the aim to improve the interpretation ofthe wrist bone anatomy. The joint registration allows a more accurateoverlay of the anatomical information of the wrist anatomy in both CT857.2. Future workand US compared to the approach where either US or CT informationwas used for registration.7.2 Future workThe thesis is an attempt for construction of a multi-object statistical wristmodel, and its registration to preoperative CT and intraoperative US forscaphoid screw axis estimation to guide a non-displaced scaphoid fracturefixation. A number of interesting areas of research can be suggested asfollows:• The presented shape model is less compact compared to the scale andpose models. The main reason is the higher number of shape featuresused here for shape representation. A more effective shape represen-tation, e.g. spherical harmonic (SPHARM) shape analysis could bebe used instead. In addition, the concept of non-linear PCA could beinvestigated for a more compact model development. Addressing thisproblem can possibly further improve accuracy and robustness of theregistration of the model to CT/US.• The proposed model is developed from a limited number of trainingsamples. Including more number of examples in the training set possi-bly introduces more shape, scale and pose variations in the developedmodel.• The main computational complexity of the registration procedure in-cludes the expectation and maximization steps of the model optimiza-tion step. Computation time can be improved by applying an effi-cient implementation of the registration algorithm using the fast Gausstransform [36] and porting the rather non-optimized Matlab code toan efficient C/C++ implementation. Efficient parallel implementationcan be also applied to the bone enhancement for further reduction ofthe runtime.• The proposed segmentation technique has been evaluated using 66 CTimages; most of them are normal and a few of them have fractures. Oneof the aims is to validate our method using more test CTs includingmore different kinds of abnormalities (e.g. implants or screws).• The current registration technique is based on a probabilistic ideathat assumes a constant amount of outlier in the target point cloud.867.2. Future workHowever, the extracted target point cloud from the US data does notgenerate same amount of outliers for all of the visible wrist bones in theultrasound scanning. An adaptive estimation of the outlier for eachwrist bone may possibly further improve the registration accuracy.• Future work includes development of a post-filtering approach (e.g.,machine learning concept could be used) on the bone enhanced volumeto remove the scattered outliers. Moreover, an automatic initializationof the wrist model to the US volume is planned to be integrated intothe model registration workflow.• We would also like to make the US images of 13 cadavers availableto the community, and want to collect the manual segmentation forthe quantitative measures. This would also allow us to apply recentlydeveloped deep neural networks to our data.• The proposed statistical model could be applied to segment the wristbones in MRI images. Though CT is usually preferred for complicateddiagnosis, MRI is taken for some cases to diagnose the fracture.• Future work should include further evaluation of this technology withinthe clinical workflow in a prospective patient study. In addition, Iwould also like to include feedback from surgeons, e.g., what kind ofinformation from the registration and screw axis computation will beuseful for the surgical guidance. The major challenges in the clinicalapplication may include the tracking errors of drilling instruments andunwanted movements of the wrist bones during drilling. In addition,it is not always possible to maintain a controlled experimental condi-tion (as like in our cadaver study) in a real surgical procedure, e.g.,movement of wrist bones during surgical procedure. To account forthis effect, I plan to include fluoroscopy in the procedure to maintainquality assurance during surgical experiments. However, I expect thatmuch smaller number of fluoroscopy shots will be required than theconventional surgical procedure.• The post-operative confirmation of the screw position inside the scaphoidbone is also a challenge in my ultrasound-based guidance. A post-operative fluoroscopy may be used for the confirmation. Future workalso includes introducing an improved screw axis estimation by jointlyoptimizing most important clinical guidelines. In addition, to mini-mize the bone movements due to the US transducer’s pressure duringUS acquisition, a better US wrist imaging setup can be used.87Bibliography[1] Scaphoid fractures. American Society for Surgery of the Hand,13(2):217–225, 1995.[2] Devin V Amin, Takeo Kanade, Anthony M Digioia, and BranislavJaramaz. Ultrasound registration of the bone surface for surgical nav-igation. Computer Aided Surgery, 8(1):1–16, 2003.[3] R Arora, M Gschwentner, D Krappinger, M Lutz, M Blauth, andM Gabl. Fixation of nondisplaced scaphoid fractures: making treat-ment cost effective. Archives of Orthopaedic and Trauma Surgery,127(1):39–46, 2007.[4] Brian Avants and James C Gee. Geodesic estimation for large defor-mation anatomical shape averaging and interpolation. Neuroimage,23:S139–S150, 2004.[5] Dean C Barratt, Carolyn SK Chan, Philip J Edwards, Graeme PPenney, Mike Slomczykowski, Timothy J Carter, and David J Hawkes.Instantiation and registration of statistical shape models of the femurand pelvis using 3D ultrasound imaging. Medical Image Analysis,12(3):358–374, 2008.[6] Dean C Barratt, Graeme P Penney, Carolyn SK Chan, Mike Slom-czykowski, Timothy J Carter, Philip J Edwards, and David J Hawkes.Self-calibrating 3D-ultrasound-based bone registration for minimallyinvasive orthopedic surgery. IEEE Transactions on Medical Imaging,25(3):312–323, 2006.[7] Maarten Beek, Purang Abolmaesumi, Suriya Luenam, Randy E Ellis,Richard W Sellens, and David R Pichora. Validation of a new surgi-cal procedure for percutaneous scaphoid fixation using intra-operativeultrasound. Medical Image Analysis, 12(2):152–162, 2008.[8] Delaram Behnami, Alexander Seitel, Abtin Rasoulian, Emran Mo-hammad Abu Anas, Victoria Lessoway, Jill Osborn, Robert Rohling,88Bibliographyand Purang Abolmaesumi. Joint registration of ultrasound, CT anda shape+ pose statistical model of the lumbar spine for guiding anes-thesia. International Journal of Computer Assisted Radiology andSurgery, pages 1–9, 2016.[9] Paul J Besl and Neil D McKay. Method for registration of 3D shapes.In Robotics-DL Tentative, pages 586–606. International Society forOptics and Photonics, 1992.[10] Kanwal K Bhatia, Joseph V Hajnal, Basant K Puri, A David Edwards,and Daniel Rueckert. Consistent groupwise non-rigid registration foratlas construction. In IEEE International Symposium on BiomedicalImaging: Nano to Macro, pages 908–911, 2004.[11] Matias Bossa, Salvador Olmos, et al. Statistical linear models inprocrustes shape space. In 1st MICCAI Workshop on MathematicalFoundations of Computational Anatomy: Geometrical, Statistical andRegistration Methods for Modeling Biological Shape Variability, pages102–111, 2006.[12] MN Bossa and Salvador Olmos. Multi-object statistical pose+shapemodels. In 4th IEEE International Symposium on Biomedical Imaging:From Nano to Macro, pages 1204–1207. IEEE, 2007.[13] Bernhard Brendel, Susanne Winter, Andreas Rick, Martin Stockheim,and Helmut Ermert. Registration of 3D CT and ultrasound datasetsof the spine using bone structures. Computer Aided Surgery, 7(3):146–155, 2002.[14] Anna Brounstein, Ilker Hacihaliloglu, Pierre Guy, Antony Hodgson,and Rafeef Abugharbieh. Towards real-time 3D US to CT bone imageregistration using phase and curvature feature based gmm matching.In Medical Image Computing and Computer-Assisted Intervention–MICCAI 2011, pages 235–242. Springer, 2011.[15] Mikael Brudfors, Alexander Seitel, Abtin Rasoulian, Andras Lasso,Victoria A Lessoway, Jill Osborn, Atsuto Maki, Robert N Rohling, andPurang Abolmaesumi. Towards real-time, tracker-less 3D ultrasoundguidance for spine anaesthesia. International Journal of ComputerAssisted Radiology and Surgery, pages 1–11, 2015.[16] Carolyn SK Chan, Dean C Barratt, Philip J Edwards, Graeme P Pen-ney, Mike Slomczykowski, Timothy J Carter, and David J Hawkes.89BibliographyCadaver validation of the use of ultrasound for 3D model instantiationof bony anatomy in image guided orthopaedic surgery. In Medical Im-age Computing and Computer-Assisted Intervention–MICCAI 2004,pages 397–404. Springer, 2004.[17] Xin Chen, Jim Graham, and Charles Hutchinson. Integrated frame-workfor simultaneous segmentation and registration of carpal bones.In 18th IEEE International Conference on Image Processing (ICIP),pages 433–436, 2011.[18] Xin Chen, Jim Graham, Charles Hutchinson, and Lindsay Muir. Auto-matic inference and measurement of 3D carpal bone kinematics fromsingle view fluoroscopic sequences. IEEE Transactions on MedicalImaging, 32(2):317–328, 2013.[19] Xin Chen, Jim Graham, Charles Hutchinson, and Lindsay Muir. Au-tomatic generation of statistical pose and shape models for articulatedjoints. IEEE Transactions on Medical Imaging, 33(2):372–383, 2014.[20] Collin Clarke, John Moore, Christopher Wedlake, Donald Lee,Su Ganapathy, Maher Salbalbal, Timothy Wilson, Terry Peters, andDaniel Bainbridge. Virtual reality imaging with real-time ultrasoundguidance for facet joint injection: a proof of concept. Anesthesia &Analgesia, 110(5):1461–1463, 2010.[21] T.F. Cootes, C.J. Taylor, D.H. Cooper, and J. Graham. Active shapemodels-their training and application. Computer Vision and ImageUnderstanding, 61(1):38 – 59, 1995.[22] Seth D Dodds, Manohar M Panjabi, and Joseph F Slade. Screw fixa-tion of scaphoid fractures: a biomechanical assessment of screw lengthand screw augmentation. The Journal of Hand Surgery, 31(3):405–413,2006.[23] J Duryea, M Magalnick, S Alli, L Yao, M Wilson, and R Goldbach-Mansky. Semiautomated three-dimensional segmentation software toquantify carpal bone volume changes on wrist CT scans for arthritisassessment. Medical Physics, 35:2321, 2008.[24] Andriy Fedorov, Reinhard Beichel, Jayashree Kalpathy-Cramer,Julien Finet, Jean-Christophe Fillion-Robin, Sonia Pujol, Christian90BibliographyBauer, Dominique Jennings, Fiona Fennessy, Milan Sonka, John Bu-atti, Stephen Aylward, James V. Miller, Steve Pieper, and Ron Kiki-nis. 3D slicer as an image computing platform for the quantita-tive imaging network. Magnetic Resonance Imaging, 30(9):1323–1341,2012.[25] P. Foroughi, E. Boctor, and M. Swartz. 2D ultrasound bone segmen-tation using dynamic programming. In IEEE Ultrasonics Symposium,pages 2523–2526, 2007.[26] Pezhman Foroughi, Danny Song, Gouthami Chintalapani, Russell HTaylor, and Gabor Fichtinger. Localization of pelvic anatomical co-ordinate system using US/atlas registration for total hip replacement.In Medical Image Computing and Computer-Assisted Intervention–MICCAI 2008, pages 871–879. Springer, 2008.[27] Ryan M Garcia and David S Ruch. Management of scaphoid fracturesin the athlete: open and percutaneous fixation. Sports Medicine andArthroscopy Review, 22(1):22–28, 2014.[28] Sahar Ghanavati, Parvin Mousavi, Gabor Fichtinger, and PurangAbolmaesumi. Phantom validation for ultrasound to statistical shapemodel registration of human pelvis. In SPIE Medical Imaging, pages79642U–79642U. International Society for Optics and Photonics, 2011.[29] Sahar Ghanavati, Parvin Mousavi, Gabor Fichtinger, Pezhman For-oughi, and Purang Abolmaesumi. Multi-slice to volume registration ofultrasound data to a statistical atlas of human pelvis. In SPIE MedicalImaging, pages 76250O–76250O. International Society for Optics andPhotonics, 2010.[30] Sean Gill, Purang Abolmaesumi, Gabor Fichtinger, JonathanBoisvert, David Pichora, Dan Borshneck, and Parvin Mousavi. Biome-chanically constrained groupwise ultrasound to CT registration of thelumbar spine. Medical Image Analysis, 16(3):662–674, 2012.[31] Sean Gill, Parvin Mousavi, Gabor Fichtinger, David Pichora, and Pu-rang Abolmaesumi. Group-wise registration of ultrasound to CT im-ages of human vertebrae. In SPIE Medical Imaging, pages 72611O–72611O. International Society for Optics and Photonics, 2009.[32] Jerome Gilles. Empirical wavelet transform. IEEE Transactions onSignal Processing, 61(16):3999–4010, 2013.91Bibliography[33] Je´roˆme Gilles, Giang Tran, and Stanley Osher. 2D empirical trans-forms. wavelets, ridgelets, and curvelets revisited. SIAM Journal onImaging Sciences, 7(1):157–186, 2014.[34] Paulo Jorge Sequeira Gonc¸alves, Pedro MB Torres, F Santos,R Anto´nio, N Catarino, and JMM Martins. A vision system forrobotic ultrasound guided orthopaedic surgery. Journal of Intelligent& Robotic Systems, 77(2):327–339, 2015.[35] Vicente Grau and J Alison Noble. Adaptive multiscale ultrasoundcompounding using phase information. In Medical Image Computingand Computer-Assisted Intervention–MICCAI 2005, pages 589–596.Springer, 2005.[36] Leslie Greengard and John Strain. The fast gauss transform. SIAMJournal on Scientific and Statistical Computing, 12(1):79–94, 1991.[37] I Hacihaliloglu, R Abugharbieh, AJ Hodgson, and P Gug. Automatedbone contour detection in 3D B-mode ultrasound images using opti-mised phase symmetry features-a clinical evaluation for pelvic frac-tures. Journal of Bone & Joint Surgery, British Volume, 94(SUPPXLIV):64–64, 2012.[38] Ilker Hacihaliloglu, Rafeef Abugharbieh, Antony J Hodgson, andRobert N Rohling. Bone surface localization in ultrasound using imagephase-based features. Ultrasound in Medicine & Biology, 35(9):1475–1487, 2009.[39] Ilker Hacihaliloglu, Rafeef Abugharbieh, Antony J Hodgson, andRobert N Rohling. Automatic adaptive parameterization in localphase feature-based bone segmentation in ultrasound. Ultrasound inMedicine & Biology, 37(10):1689–1703, 2011.[40] Ilker Hacihaliloglu, Rafeef Abugharbieh, Antony J Hodgson, Robert NRohling, and Pierre Guy. Automatic bone localization and fracturedetection from volumetric ultrasound images using 3D local phase fea-tures. Ultrasound in Medicine & Biology, 38(1):128–144, 2012.[41] Ilker Hacihaliloglu, Abtin Rasoulian, Robert N Rohling, and PurangAbolmaesumi. Local phase tensor features for 3D ultrasound to sta-tistical shape+pose spine model registration. IEEE Transactions onMedical Imaging, 33(11):2167–2179, 2014.92Bibliography[42] Ilker Hacihaliloglu, David R Wilson, Michael Gilbart, Michael A Hunt,and Purang Abolmaesumi. Non-iterative partial view 3D ultrasoundto CT registration in ultrasound-guided computer-assisted orthopedicsurgery. International Journal of Computer Assisted Radiology andSurgery, 8(2):157–168, 2013.[43] Peter Heinze, Dietmar Meister, Rudolf Kober, Jo¨rg Raczkowsky, andHeinz Worn. Atlas-based segmentation of pathological knee joints.Studies in Health Technology and Informatics, pages 198–203, 2002.[44] Yan Chai Hum. Segmentation of hand bone for bone age assessment.Springer, 2013.[45] Gelu Ionescu, Ste´phane Lavalle´e, and Jacques Demongeot. Automatedregistration of ultrasound with ct images: application to computer as-sisted prostate radiotherapy and orthopedics. In Medical Image Com-puting and Computer-Assisted Intervention–MICCAI99, pages 768–777. Springer, 1999.[46] Paul Jaccard. The distribution of the flora in the alpine zone. Newphytologist, 11(2):37–50, 1912.[47] Ameet K Jain and Russell H Taylor. Understanding bone responsesin B-mode ultrasound images and automatic bone surface extractionusing a bayesian probabilistic framework. In Medical Imaging, pages131–142. International Society for Optics and Photonics, 2004.[48] Klas Josephson, Anders Ericsson, and Johan Karlsson. Segmentationof medical images using three-dimensional active shape models. In14th Scandinavian Conference on Image Analysis, SCIA, pages 719–728. Springer, 2005.[49] Athanasios Karamalis, Wolfgang Wein, Tassilo Klein, and NassirNavab. Ultrasound confidence maps using random walks. MedicalImage Analysis, 16(6):1101–1112, 2012.[50] Siavash Khallaghi, Parvin Mousavi, Ren Hui Gong, Sean Gill,Jonathan Boisvert, Gabor Fichtinger, David Pichora, Dan Borsch-neck, and Purang Abolmaesumi. Registration of a statistical shapemodel of the lumbar spine to 3D ultrasound images. In Medical ImageComputing and Computer-Assisted Intervention–MICCAI 2010, pages68–75. Springer, 2010.93Bibliography[51] P Kilian, C Plaskos, S Parratte, J-NA Argenson, E Stindel, J Tonetti,and S Lavallee. New visualization tools: computer vision and ul-trasound for mis navigation. The International Journal of MedicalRobotics and Computer Assisted Surgery, 4(1):23–31, 2008.[52] Martin Koch, Alexander G Schwing, Dorin Comaniciu, and MarcPollefeys. Fully automatic segmentation of wrist bones for arthri-tis patients. In 2011 IEEE International Symposium on BiomedicalImaging: From Nano to Macro, pages 636–640. IEEE, 2011.[53] Jens Kowal, Christoph Amstutz, Frank Langlotz, Haydar Talib,and Miguel Gonzalez Ballester. Automated bone contour detectionin ultrasound b-mode images for minimally invasive registration incomputer-assisted surgery: an in vitro evaluation. The InternationalJournal of Medical Robotics and Computer Assisted Surgery, 3(4):341–348, 2007.[54] Hans Lamecker, Martin Seebass, Hans-Christian Hege, and Peter Deu-flhard. A 3D statistical shape model of the pelvic bone for segmenta-tion. In Medical Imaging 2004, pages 1341–1351. International Societyfor Optics and Photonics, 2004.[55] Evan L Leventhal, Scott W Wolfe, Eric F Walsh, and Joseph J Crisco.A computational approach to the “optimal” screw axis location andorientation in the scaphoid bone. The Journal of Hand Surgery,34(4):677–684, 2009.[56] Seth Levitz and David Ring. Retrograde (volar) scaphoid screwinsertion–a quantitative computed tomographic analysis. The Journalof Hand Surgery, 30(3):543–548, 2005.[57] Philippe A Liverneaux, Anes Gherissi, and Matthieu Beustes Ste-fanelli. Kirschner wire placement in scaphoid bones using fluoroscopicnavigation: a cadaver study comparing conventional techniques withnavigation. The International Journal of Medical Robotics and Com-puter Assisted Surgery, 4(2):165–173, 2008.[58] Ernest L Madsen, Maritza A Hobson, Hairong Shi, Tomy Varghese,and Gary R Frank. Tissue-mimicking agar/gelatin materials for usein heterogeneous elastography phantoms. Physics in Medicine andBiology, 50(23):5597, 2005.94Bibliography[59] Daniel Maleike, Marco Nolden, H-P Meinzer, and Ivo Wolf. Interactivesegmentation framework of the medical imaging interaction toolkit.Computer Methods and Programs in Biomedicine, 96(1):72–83, 2009.[60] G Meermans and Frederik Verstreken. Percutaneous transtrapezialfixation of acute scaphoid fractures. Journal of Hand Surgery (Euro-pean Volume), 33(6):791–796, 2008.[61] Geert Meermans, Francis Van Glabbeek, Marc J Braem, Roger P vanRiet, Guy Hubens, and Frederik Verstreken. Comparison of two percu-taneous volar approaches for screw fixation of scaphoid waist fractures.The Journal of Bone & Joint Surgery, 96(16):1369–1376, 2014.[62] Stephen W Meldon and Stephen W Hargarten. Ligamentous injuries ofthe wrist. The Journal of Emergency Medicine, 13(2):217–225, 1995.[63] Kurt A Menapace, Lawrence Larabee, Steven P Arnoczky, Vivek SNeginhal, A George Dass, and Lawrence M Ross. Anatomic placementof the herbert-whipple screw in scaphoid fractures: a cadaver study.The Journal of Hand Surgery, 26(5):883–892, 2001.[64] AB Mink van der Molen, JW Groothoff, GJP Visser, PH Robinson,and WH Eisma. Time off work due to scaphoid fractures and othercarpal injuries in the netherlands in the period 1990 to 1993. TheJournal of Hand Surgery: British & European Volume, 24(2):193–198,1999.[65] Douglas C Moore, Joseph J Crisco, Theodore G Trafton, and Evan LLeventhal. A digital database of wrist bone anatomy and carpal kine-matics. Journal of Biomechanics, 40(11):2537–2542, 2007.[66] John Moore, Colin Clarke, Daniel Bainbridge, Chris Wedlake, AndrewWiles, Danielle Pace, and Terry Peters. Image guidance for spinalfacet injections using tracked ultrasound. In Medical Image Computingand Computer-Assisted Intervention–MICCAI 2009, pages 516–523.Springer, 2009.[67] Andriy Myronenko and Xubo Song. Point set registration: Coher-ent point drift. IEEE Transactions on Pattern Analysis and MachineIntelligence, 32(12):2262–2275, 2010.[68] Ryan Nishihara. The dilemmas of a scaphoid fracture: a difficultdiagnosis for primary care physicians. Hospital Physician, 36(3):24–42, 2000.95Bibliography[69] Nobuyuki Otsu. A threshold selection method from gray-level his-tograms. Automatica, 11(285-296):23–27, 1975.[70] Victoria S Pao and James Chang. Scaphoid nonunion: diagnosisand treatment. Plastic and Reconstructive Surgery, 112(6):1666–1677,2003.[71] Vladimir Pekar, Michael R Kaus, Cristian Lorenz, Steven Lobregt,Roel Truyen, and Juergen Weese. Shape-model-based adaptation of3D deformable meshes for segmentation of medical images. In MedicalImaging 2001, pages 281–289. International Society for Optics andPhotonics, 2001.[72] Graeme P Penney, Dean C Barratt, Carolyn SK Chan, Mike Slom-czykowski, Timothy J Carter, Philip J Edwards, and David J Hawkes.Cadaver validation of intensity-based ultrasound to CT registration.Medical Image Analysis, 10(3):385–395, 2006.[73] Stephen M Pizer, Ja-Yeon Jeong, Conglin Lu, Keith Muller, andSarang Joshi. Estimating the statistics of multi-object anatomic geom-etry using inter-object relationships. In Deep Structure, Singularities,and Computer Vision, pages 60–71. Springer, 2005.[74] Bahbibi Rahmatullah, Aris T Papageorghiou, and J Alison Noble.Integration of local and global features for anatomical object detectionin ultrasound. In Medical Image Computing and Computer-AssistedIntervention–MICCAI 2012, pages 402–409. Springer, 2012.[75] Kashif Rajpoot, Vicente Grau, and J Alison Noble. Local-phase based3D boundary detection using monogenic signal and its application toreal-time 3D echocardiography images. In IEEE International Sympo-sium on Biomedical Imaging: From Nano to Macro, 2009. ISBI’09.,pages 783–786. IEEE, 2009.[76] Anand Rangarajan, Haili Chui, and Fred L Bookstein. The softassignprocrustes matching algorithm. In Information Processing in MedicalImaging, pages 29–42. Springer, 1997.[77] Abtin Rasoulian, Purang Abolmaesumi, and Parvin Mousavi. Feature-based multibody rigid registration of ct and ultrasound images of lum-bar spine. Medical Physics, 39(6):3154–3166, 2012.96Bibliography[78] Abtin Rasoulian, Robert Rohling, and Purang Abolmaesumi. Group-wise registration of point sets for statistical shape models. IEEETransactions on Medical Imaging, 31(11):2025–2034, 2012.[79] Abtin Rasoulian, Robert Rohling, and Purang Abolmaesumi. Lum-bar spine segmentation using a statistical multi-vertebrae anatom-ical shape+pose model. IEEE Transactions on Medical Imaging,32(10):1890–1900, 2013.[80] Abtin Rasoulian, Robert N Rohling, and Purang Abolmaesumi. Astatistical multi-vertebrae shape+ pose model for segmentation of CTimages. In SPIE Medical Imaging, pages 86710P–86710P, 2013.[81] Abtin Rasoulian, Alexander Seitel, Jill Osborn, Samira Sojoudi,Saman Nouranian, Victoria A Lessoway, Robert N Rohling, and Pu-rang Abolmaesumi. Ultrasound-guided spinal injections: a feasibilitystudy of a guidance system. International Journal of Computer As-sisted Radiology and Surgery, 10(9):1417–1425, 2015.[82] Steven J Rhemrev, Daan Ootes, Frank JP Beeres, Sven AG Meylaerts,and Inger B Schipper. Current methods of diagnosis and treatment ofscaphoid fractures. Int J Emerg Med, 4(4), 2011.[83] Steffen Schumann. State of the art of ultrasound-based registration incomputer assisted orthopedic interventions. In Computational Radi-ology for Orthopaedic Interventions, pages 271–297. Springer, 2016.[84] Steffen Schumann, Lutz-P Nolte, and Guoyan Zheng. Compensationof sound speed deviations in 3D B-mode ultrasound for intraoperativedetermination of the anterior pelvic plane. IEEE Transactions onInformation Technology in Biomedicine, 16(1):88–97, 2012.[85] Steffen Schumann, Lutz-P Nolte, and Guoyan Zheng. Determinationof pelvic orientation from sparse ultrasound data for tha operated inthe lateral position. The International Journal of Medical Roboticsand Computer Assisted Surgery, 8(1):107–113, 2012.[86] Thomas B Sebastian, Hu¨seyin Tek, Joseph J Crisco, and Benjamin BKimia. Segmentation of carpal bones from CT images using skeletallycoupled deformable models. Medical Image Analysis, 7(1):21–45, 2003.[87] Erin J Smith, Hisham A Al-Sanawi, Braden Gammon, Paul J St John,David R Pichora, and Randy E Ellis. Volume slicing of cone-beam97Bibliographycomputed tomography images for navigation of percutaneous scaphoidfixation. International Journal of Computer Assisted Radiology andSurgery, 7(3):433–444, 2012.[88] E Stindel, JL Briard, Ph Merloz, S Plaweski, F Dubrana, C Lefevre,and J Troccaz. Bone morphing: 3D morphological data for total kneearthroplasty. Computer Aided Surgery, 7(3):156–168, 2002.[89] Martin Styner, Kevin Gorczowski, Tom Fletcher, Ja Yeon Jeong,Stephen M Pizer, and Guido Gerig. Statistics of pose and shape inmulti-object complexes using principal geodesic analysis. In MedicalImaging and Augmented Reality, pages 1–8. Springer, 2006.[90] Jerome Tonetti, Lionel Carrat, Sorin Blendea, Philippe Merloz, Jo-celyne Troccaz, Stephane Lavalle´e, and Jean-Paul Chirossel. Clinicalresults of percutaneous pelvic surgery. computer assisted surgery us-ing ultrasound compared to standard fluoroscopy. Computer AidedSurgery, 6(4):204–211, 2001.[91] Denis Tran and Robert N Rohling. Automatic detection of lumbaranatomy in ultrasound images of human subjects. IEEE Transactionson Biomedical Engineering, 57(9):2248–2256, 2010.[92] Thomas E Trumble, Todd Clarke, and Hans J Kreder. Non-union ofthe scaphoid. treatment with cannulated screws compared with treat-ment with herbert screws. The Journal of Bone & Joint Surgery,78(12):1829–37, 1996.[93] Martijn van de Giessen, Mahyar Foumani, Frans M Vos, Simon DStrackee, Mario Maas, Lucas J Van Vliet, Cornelis A Grimbergen,and Geert J Streekstra. A 4D statistical model of wrist bone mo-tion patterns. IEEE Transactions on Medical Imaging, 31(3):613–625,2012.[94] Vladimir Vezhnevets and Vadim Konouchine. Growcut: Interactivemulti-label n-d image segmentation by cellular automata. In Proc. ofGraphicon, pages 150–156, 2005.[95] Eric Walsh, Joseph J Crisco, and Scott W Wolfe. Computer-assistednavigation of volar percutaneous scaphoid placement. The Journal ofHand Surgery, 34(9):1722–1728, 2009.98Bibliography[96] Simon K Warfield, Kelly H Zou, and William M Wells. Simultaneoustruth and performance level estimation (staple): an algorithm for thevalidation of image segmentation. Medical Imaging, IEEE Transac-tions on, 23(7):903–921, 2004.[97] Wolfgang Wein, Shelby Brunke, Ali Khamene, Matthew R Callstrom,and Nassir Navab. Automatic CT-ultrasound registration for diagnos-tic imaging and image-guided intervention. Medical Image Analysis,12(5):577–585, 2008.[98] Wolfgang Wein, Athanasios Karamalis, Adrian Baumgartner, andNassir Navab. Automatic bone detection and soft tissue awareultrasound-CT registration for computer-aided orthopedic surgery.International Journal of Computer Assisted Radiology and Surgery,10(6):971–979, 2015.[99] Wolfgang Wein, Ali Khamene, Dirk-Andre´ Clevert, Oliver Kutter, andNassir Navab. Simulation and fully automatic multimodal registrationof medical ultrasound. In Medical Image Computing and Computer-Assisted Intervention–MICCAI 2007, pages 136–143. Springer, 2007.[100] Susanne Winter, Bernhard Brendel, Ioannis Pechlivanis, KirstenSchmieder, and Christian Igel. Registration of CT and intraoperative3D ultrasound images of the spine using evolutionary and gradient-based methods. IEEE Transactions on Evolutionary Computation,12(3):284–296, 2008.[101] Scott W. Wolfe. Scaphoid fractures and non-union: Wrist fracturesand treatment. 2009.[102] Charles XB Yan, Benoˆıt Goulet, Julie Pelletier, Sean Jy-Shyang Chen,Donatella Tampieri, and D Louis Collins. Towards accurate, robustand practical ultrasound-CT registration of vertebrae for image-guidedspine surgery. International Journal of Computer Assisted Radiologyand Surgery, 6(4):523–537, 2011.[103] Aifeng Zhang, Arkadiusz Gertych, Sylwia Kurkowska-Pospiech,Brent J Liu, and HK Huang. Carpal bone analysis in bone age assess-ment. In SPIE Medical Imaging, pages 61450Z–61450Z, 2006.99

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-0340302/manifest

Comment

Related Items