UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 SKEL ETAL KINEMATIC MEASUREMENT USING MODELIMAGE REGISTRATION AND MECHANICAL CONSTRAINTS By SHANG MU A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORID A IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2010 1 PAGE 2 2010 Shang Mu 2 PAGE 3 ACKNOWLEDGMENTS I wish to express sincere appreciation to Dr. Banks, my advisor, who gave me enormous support throughout my quest. I thank J.D. Yamokoski for his preliminary work, discussion on the design of the soft ware framework, and general support on my work. I thank Nick Dunbar for building the bi plane calibration cube and for his help in the biplane experiment of natural joints and in various other situations. I thank Brody Huval for writing the biplane calibration program with Nicks help, Roshni Marballi, Michael DesRosiers, and Kyle Houston for performing measurement s, and Brian Maas for running statistics. Special thanks go to my father, who fostered my enthusiasm in science. 3 PAGE 4 TABLE OF CONTENTS page ACKNOWLEDG MENTS..................................................................................................3 LIST OF TABLES............................................................................................................6 LIST OF FI GURES..........................................................................................................7 LIST OF ABBR EVIATIONS.............................................................................................8 ABSTRACT .....................................................................................................................9 CHAPTER 1 INTRODUC TION....................................................................................................11 Backgroun d.............................................................................................................11 Research Ov erview.................................................................................................14 A Quantitative Evaluation of the Influences on the Accuracy from Various Effects...........................................................................................................14 Articular Constraints.........................................................................................15 Passive Elastic Joint Constrai nts......................................................................15 2 INFLUENCE OF BONE MODEL, NUMBER OF IMAGE PLANES, AND OBSERVERS ON THE ACCURACY OF JOINT KINEMATICS FROM RADIOGRAPHIC MODELIM AGE REGISTRATION..............................................18 Introducti on.............................................................................................................18 Methods ..................................................................................................................20 Result s....................................................................................................................23 Discussio n..............................................................................................................25 3 ARTICULAR CO NSTRAINT S.................................................................................34 Introducti on.............................................................................................................34 Methods ..................................................................................................................36 Penetration Comput ation Algor ithm..................................................................36 Image Comparis on Metr ic................................................................................38 Optimizers and Co st Func tions ........................................................................38 Experimental Setup..........................................................................................40 Result s....................................................................................................................42 Discussio n..............................................................................................................42 Acknowled gements.................................................................................................45 4 LIGAMENTLIKE JOIN T CONSTRAI NTS...............................................................49 4 PAGE 5 Introducti on .............................................................................................................49 Methods ..................................................................................................................50 Data se ts..........................................................................................................50 Registration method.........................................................................................51 Performance ev aluati on....................................................................................54 Result s....................................................................................................................55 Discussio n..............................................................................................................56 5 CONCLUS ION........................................................................................................62 APPENDIX: OPTIMIZER PARAMETER VALUES ........................................................64 LIST OF RE FERENCES...............................................................................................65 BIOGRAPHICAL SKETCH ............................................................................................69 5 PAGE 6 LIST OF TABLES Table page 21 Actual range of motion in the 10 image pairs per joint, measured using biplane bead registrations .....................................................................................30 23 RMS differences between bead a nd bone models for all observers and singleand biplane measuremen ts....................................................................30 24 Average difference of each observer fr om the mean of all observers, for all image data, including singleand bi plane measurements, and bead and bone models .......................................................................................................31 25 RMS difference between obs ervers for biplane registration of bead models.....31 26 RMS differences between biplane bead model measurements and singleplane bead model measurem ents......................................................................31 27 RMS differences between biplane bead model measurements and biplane bone model meas urement s................................................................................32 28 RMS differences between biplane bead model measurements and singleplane bone model measurem ents......................................................................32 31 Success rates and RMS registration errors for three methods with and without the contact function penalty or cons traint...............................................46 41 Success rates and the relative pose RMS deviations for two datasets with and without contact and/or ligament penalty. .....................................................59 A1 Particle swarm pattern s earch parameter values ................................................64 A2 RALG optimizer par ameter va lues. .....................................................................64 A3 LevenbergMarquardt nonlinear least squares optimizer parameter values.......64 6 PAGE 7 LIST OF FIGURES Figure page 11 Models of bones or implants can be registered with their radiographic projections to provide accurate m easures of dynamic 3D motion.......................17 12 Basic 3D2D registration process for m easuring 3D poses from 2D image(s)....17 21 Biplane registration of elbow bone models. ........................................................33 31 Knee implant models matched to the fluoroscopic image to measure their spatial pos es.......................................................................................................47 32 Illustration of inside/outs ide determination for articulating polygonal surfaces...47 33 Standard deviation of the errors in the measured relative joint position and joint angl es.........................................................................................................48 34 Mean errors in the measured relati ve joint position and joint angl es..................48 41 The symmetric view problem involv es nearly equivalent projections of symmetric obj ects...............................................................................................60 42 Coordinate system definition for t he femoral and tibial implant models..............60 43 Silhouette of illregister ed tibia component model s............................................61 44 An image from Datase t B...................................................................................61 7 PAGE 8 LIST OF ABBREVIATIONS CT Computed tomography. DOF Degree of freedom. LSQ Least squares (optimization). MR Magnetic resonance (imaging). Pswarm Particle swarm pattern search algorithm. RMS Root mean square. 8 PAGE 9 Abstract of Dissertation Pr esented to the Graduate School of the University of Fl orida in Partial Fulf illment of the Requirements for t he Degree of Doctor of Philosophy SKELETAL KINEMATIC MEASUREMENT USING MODELIMAGE REGISTRATION AND MECHANICAL CONSTRAINTS By Shang Mu December 2010 Chair: Scott A. Banks Major: Mechanical Engineering Skeletal kinematic measurement usi ng modelimage registration has been practiced for more than fifteen years and prov en useful elucidating the dynamic function of human joints with and without arthroplasty However, technical complexity and the time consuming nature of the analysis have restricted utilization to a few research groups, and improvements on the methods have been few. To facilitate the use of modelimage registration for skeletal ki nematics studies and to enhance research collaboration, we have been developing Jo intTrack, an opensource, expandable software platform for radiographic modelimage registration. Using this software, in the first part of this dissertation, we eval uate the influence of bone geometry, number of image planes, and the human operator on the accuracy of joint kinematics measurements, and identif y limitations and weaknesses in the methods. In the later parts of this dissertation, we investigate new techniques to improve measurement accuracy and reduce human labor. Compared to biplane (stereo) measurement, singleplane setup dominates the field due to physical constraint or equipment availability, but is intrinsically less accurate to a degree that limits its application, and delivers many difficulties to the registration process. Improvement on 9 PAGE 10 10 the accuracy in singleplane measurement was achieved by incorporating contact constraint into the registration process. Further incorporation of ligamentlike force models productively improved registration success rate, decreasing the need for human intervention in initial setup and/or latestage va lidation and correction. PAGE 11 CHA PTER 1 INTRODUCTION Background Register: 5.b. v.trans. To position with precision, in order to ensure an exact correspondenc e of parts ; to align. The Ox ford English Dictionary ( http://dictionary.oed.com/ accessed November 2010). Modelimage registration (Fig. 11) for the measurement of dynamic skeletal motion has been used for orthopaedic research since 19 91 (Banks et al., 1991). In this technique, a 3D model of a bone or implant is projected ont o a digitized radiographic image and its 3D pose is adj usted until the projection matches the image. The technique gained popularity due to its noninvasive nature, accuracy, and capability of obtaining measurements for dynamic in vivo (i.e. real life) ac tivities. Measurement precision on the order of 1mm translation or 1 rotation provides information that is useful to improve understanding of dynamic joint instability (e.g. Yamaguchi et al., 2009), improvement of joint replacement desi gns (e.g. Banks et al., 1997), and dynamic simulation of joint mechani cs (e.g. Zhao et al., 2007). In this technique, Xray images of a mo ving joint are first obtained from systems capable of continuously capt uring images, known generally as fluoroscopy systems. Given the 3D model of an object (a bone or an implant) and the radiographic projection parameters, a computer simu lated perspective projection of the 3D model is generated and compared with the actual xray images us ing an image similarity measure. Given a single frame radiographic image and a 3D model, a search pr ocedure, often a numerical optimization routine, is then used to find the best set of 3D model pose parameters that maximize the similarity measure (Fig. 12). This procedure is repeated for each bone or implant and for each frame to obtain the comp lete temporal kinematics of the objects. 11 PAGE 12 Numerous groups dev eloped measur ement tools based upon modelimage registration and applied these in studies of natural and prosthetic joints. In these reports, generation of perspective projection images mostly have been standard computer graphics routines using surface models of bone or implants (Banks and Hodge, 1996; Zuffi et al., 1999; Valstar et al., 2001; Mahfouz et al., 2003; Li et al., 2004; Yamazaki et al., 2004; Hirokawa et al., 2008; Prins et al., 2010), and recent reports have used digital reconstructed radiographs (DRR) from CT volumes (You et al., 2001; Tang et al., 2004). Image comparison metr ics have been based on correspondence of image edges (Banks and Hodge, 1996; Valsta r et al., 2001; Yamazaki et al., 2004; Hirokawa et al., 2008), wei ghted edge and grayscale correlation (Mahfouz et al., 2003), general normalized correlation (You et al., 2001), and gradient correlation (Tang et al., 2004). A different class of methods using 3D distance maps, based on efficient raytracing techniques, also have been reported ( Lavallee and Szeliski, 1995; Zuffi et al., 1999). A wide range of optimizat ion algorithms and numerical solvers have been utilized for modelimage registration, including Lev enbergMarquardt nonli near least squares algorithm (Zuffi et al., 1999), combined NelderMead downhill si mplex and simulated annealing method (Mahfouz et al ., 2003), feasible sequenti al quadratic programming (FSQP) (Valstar et al., 2001), DurandKern er method (Yamazaki et al., 2004), simple searches within lowdimens ion precomputed shapelibraries (Banks and Hodge, 1996; Hirokawa et al., 2008), and interactive manual matching by a human operator using solidmodeling software (Li et al., 2004). The ma jority of these reports provide careful evaluations of measurement accuracy, typically fo r a specific joint of interest or specific research application. However, none of these reports have quantified measurement 12 PAGE 13 accuracy with multiple factors and observers (human operators) us i ng consistent image sets on multiple natural joints. Kinematic measurement accu racy for modeimage registration techniques varies little with the details of the image processing and optimization algorithms, but is affected significantly by singleplane or biplane im aging techniques. Singleplane imaging is monocular vision, which is fundamentally insensitive for translations along the viewing axis (Zaxis in our standard coordinate syst em). Biplane imaging techniques capture images from two views simultaneously, provid ing stereo visualization and more uniform and precise measurements than singleplane imaging. However, biplane systems often have very limited space within which an activi ty can be performed, they double the xray exposure and are not generally available in hospitals. Sing leplane imaging systems are found in every hospital, and so are readily available for dynamic imaging. However, there have been few reports of methods that improve singleplane registration precision (Yamazaki et al., 2004). The robustness of kinematic measurement techniques rarely has been reported. Mahfouz et al. (2003) reported 50% registra tion success for knee replacement implants when the initial guess was 16 mm from the inplane soluti on, 128mm from the outofplane solution and 16 from the rotation so lution. Their registration success rate increased to over 90% when the initial guesses were within 4mm, 32mm and 4, respectively. The fact remains that all reported techniques for skeletal kinematic analysis rely on an operatorprovided guess to start the registrati on process for every image in the sequence. The laborintensive nature of these measurements is a significant factor limiting their use in re search and preventing their clinical use. 13 PAGE 14 Image quality and object shape (bone or implant) strongly affect modelimage registration success rates. For this reason, some bones are extremely difficult to measure. Patella kinematics, as an exampl e, are hard to measure due to the round shape, small size and the poor image contrast of the patella, and have only recently been measured using a biplane xray system (Bey et al., 2008). The bones of the foot/ankle, the talus and calcaneus in particu lar, also present measurement challenges due to their oblong rounded shapes and overlapping joints (poor contrast images), yet there is strong interest to quantify foot/ankle motion in pathologic feet during walking. In this dissertation, we evaluated the performance of a basic modelimage registration method that has been used by our group and identify limitations and weaknesses in the method. With these limitations clearly identified, we present two approaches to enhance accuracy and success rates for 3Dto2D modelimage registration. Research Overview A Quantitative Evaluation of the Influenc es on the Accurac y from Various Effects Biplane images were taken from four different natural joints in vitro Registration was performed by multiple human operators for four different conditions: biplane images with models of metallic beads implan ted in the bones, biplane images with 3D bone models, singleplane image s with bead models and singl eplane images with 3D bone models. Accuracy for each condition wa s quantified and analyz ed with statistical models. We specifically evaluated the infl uence of registration model (bones versus beads), number of image planes (biplane vers us singleplane), and the effect of human observers on the accuracy of joint ki nematics from radiographic modelimage registration, and the results are docum ented and discussed in Chapter 2. 14 PAGE 15 Articular Constraints Contact is the most important aspect of human joints, but is rarely used in modelimage registration methods to improve accura cy (Hirokawa et al., 2008, Prins et al., 2010). Hirokawa et al. (2008) in corporated contact into the registration algorithm by solving simultaneous equations for the contac t point(s). Because this method requires the articulating surfaces to be described in para metric equations, different joints need to be handled differently, and it is difficult to apply the method to complex contacts, especially natural joints wit hout prostheses. Prins et al (2010) reported a technique similar to the one proposed in this dissertati on. However, neither of these two methods can handle deformation, wear, or inaccuracy of the 3D surface models. The method proposed in Chapter 3 addresses these issues nat urally. In brief, t he method integrates contact constraints into the modelimage registration process to improve measurement in the outofplane direction in singlepla ne registration. In addition, the contact constraint dramatically shrinks the configur ation space of the optimization problem and therefore increases the rate of successful registrations in certain cases. Passive Elastic Joint Constraints The motions within human joints are c onstrained by the ligaments and softtissues surrounding each joint. These softtissues, even when damaged or diseased, define a passive envelope of possible motions. Kno wledge of this passive envelope might be useful to constrain the model image registration process to evaluate only physiologically feasible joint configurations. In Chapter 4, we incorporate in to the modelimage registration process passive joint constrai nts modeled as restoring ligament forces. These constraints reduced the need for accu rate initial guesses and significantly increased the rate of successful registrations. In particular, false extrema due to object 15 PAGE 16 symmetry were effectively avoided using suc h constraints. These ligamentlike constraints also have the potent ial reduce the requirements fo r careful human review of each registration result. We are not aware of any existing reports of similar methods. 16 PAGE 17 17 Figure 11. Models of bones or implants can be registered with their radiographic projections to provide accurate measures of dynamic 3D motion. Figure 12. Basic 3D2D registration process for measuring 3D pos es from 2D image(s). PAGE 18 CHA PTER 2 INFLUENCE OF BONE MODEL, NUMBER OF IMAGE PLANES, AND OBSERVERS ON THE ACCURACY OF J OINT KINE MATICS FROM RADIOGRAPHIC MODELIMAGE REGISTRATION Introduction Modelimage registration for in vivo skeletal kinematic measurement has become a very popular research tool since it was first reported i n 1991 (Banks and Hodge, 1992). For example, an ISI Web of Knowledge keyword search for the terms fluoro* and kinematics yields 106 relevant peerrevi ewed papers in archival journals from January 2008 to October 2010. There are a wide range of techniques implemented to perform these measurements, bu t they all share the commo n elements of (1) one or more radiographic projections with measured geometry, (2) threedimensional models of one or more moving bones or implants, and (3) a procedure to register the models with the radiographic projections by com puter synthesis of candidate projections, comparison with a metric and iterative num erical refinement (e.g. Banks and Hodge, 1996; You et al., 2001; Kaptein et al., 2003; Mahfouz et al., 2003). These techniques have been applied to quantify in vivo the mo vement of, for example, total knee replacements (e.g. Banks et al., 1997; Stiehl et al., 1999; Zuffi et al., 1999; Kaptein et al., 2003; Yamazaki et al., 2004), natural knees (e.g. Komistek et al., 2003; Li et al., 2007; Tashman et al., 2007; Morooka et al., 2008), natural and replaced shoulders (Kon et al., 2008; Nishinaka et al., 2008), hi p replacements (Komistek et al., 2002; Tanino et al., 2008), natural and replaced ankles (Yamaguchi et al., 2009), elbows (Matsuki et al., 2010) and spinal vertebrae (W ang et al., 2008; McDonald et al., 2010). Several studies have reported careful evaluat ions of measurement accuracy with knee implants (e.g. Banks and Hodge, 1996; Kaptei n et al., 2003; Mahfouz et al., 2003; 18 PAGE 19 Kaptein et al., 2007), yet there remain many issues associated wit h modelimage registration measurement techniques for jo ints without implants that have not been addressed. Important questions remain conc erning the relative accuracy of singleprojection and stereo measur ements, the influence of bone shape on measurement accuracy, and how multiple users affect measurement precision. Modelimage registration based measurements have t he potential to provide useful and clinically relevant information if sufficiently well developed. Technical challenges for a clinical implementat ion will include the need to automate the measurement, to characterize the measurement applied to different joints with and without implants, and to minimize t he influence of a human operator on the measurement procedure. Alt hough not yet a fully automatic and robust measurement capability, it remains worthwhile to quantif y present technical capabilities to inform research efforts and to guide priorities for futu re developments. In th is study, we sought to answer the following three questions: 1. How much does joint kinematics measur ement accuracy decrease using a singleplane imaging technique compared to a stereo imaging technique? 2. How is joint kinematics measurement a ccuracy affected by the shape model, i.e. models of bones versus models of discrete marker clusters? 3. How is joint kinematics measurement accuracy affected by the user? Answers to these questions put into pers pective reports currently appearing in the research literature and suggest how development efforts should be prioritized to reach a clinically useful in vivo skeleta l kinematic measurement capability. 19 PAGE 20 Methods A stereo flatpanel radiographic imaging syst em was used for this study (Allura Xper FD20/20, Philips Medical Systems B.V. ). All images were obtained at 15 frames per second with 3ms exposures. Beam energy and current were adjusted for each case prior to recording dynamic trials. The photogrammetric geometry of the two image systems was determined using a cubic threedimensional (3D) calibration object consisting of plastic plates with embedded me tallic spheres. The 3D sphere locations were determined by merging reconstructed bead locations from a pair of CT scans acquired with the calibration object in two or thogonal orientations relative to the scan direction. Repeated calibrations yielded imageplane bead reconstruction residuals of 0.3mm. One lower extremity and one forequarter cadaveric specimen were properly obtained and prepared for the study. Three to fi ve metallic beads (2mm Pb) were placed in the scapula, proximal and distal humerus, radius, ulna, femur, patella, proximal and distal tibia, talus and calcaneus. No ski n or softtissue was removed from the specimens. CT scans with 0.5mm slice thi ckness were acquired for the specimens (Sensation 16, Siemens) and the 3D surface geometries of the bones and the implanted beads were reconstructed using opensource software (Yushkevich et al., 2006). The 3D models were aligned with standard joint reference systems (Morooka et al., 2008; Nishinaka et al., 2008; Yamaguchi et al., 2009; Matsuki et al., 2010). The knee, lower shank and foot, and forequarter specimens were mounted to radiolucent (xray transparent) fixtures such that each join t could be taken through a range of motion during radiographic observation. 20 PAGE 21 Image detectors were aligned wit h each join t so that one plane was a lateral or anteroposterior projection as would commonly be used for a singleplane observation. The other detector was aligned to obtain a cl ear view of the moving joint with an angular separation of 5060 from the first detector. Image sequences were obtained for each joint as it was manually moved through a range of motion that included primary and secondary rotations of the joint. For each jo int a group of ten image pairs were chosen to cover the joint range of motion and thes e images were placed in a randomly ordered sequence. Modelimage registration was performed us ing custom opensour ce software with capability for both singleplane and stereo image analysis (JointTrack, www.sourceforge.net/projects/jointtrack ). The user adjusted grayscale parameters and edge enhancement (Canny, 1986) for individual images or the entire sequence. The user selected a model and manually pos itioned it in 3D space using the mouse and keyboard while viewing synthetic radiographic projections and a 3D viewing window showing the relative position of all the bones in the joint (Fig. 21). With this initial pose as a starting point, the modelimage regi stration was iteratively refined using a simulated annealing global optim ization routine (Kirkpatrick et al., 1983; Goffe et al., 1994) and a costfunction summing the crosscorrelations of the grayscale and edgedetected (Canny, 1986) versions of the actual and synt hesized radiographic projections. At least three, and in most cases f our, human users performed the modelimage registration measurements for each image sequence. Users ranged from orthopaedic surgeons who were very familiar with the anatomy and experienced using the software, 21 PAGE 22 to undergraduate engineering stu dent s with little anatomic familiarity and only basic instruction in how to use the measurement software. Users performed the modelimage measurement process independently four times for each image sequence with the following conditions: 1. Stereo images with 3D models of the bones with the metallic beads. 2. Stereo images with 3D bone models without beads. 3. Singleplane images with 3D model s of the bones with the metallic beads. 4. Singleplane images wit h 3D bone models without beads. The resulting 3D bone kinematics were used to compute joint kinematics in the anatomic reference frames for the patellofemo ral, tibiofemoral, tibiotalar, subtalar, glenohumeral, radiohumeral, and radioulnar joints. Mean values were compared using an analysis of variance model (The R Pr oject for Statistical Computing, www.rproject.org ) without interaction between factors: ),0(2 __ ....N yy y yyylframe kobserver boneorbeadj planesni ijkl (Eq. 21) where y is any measured pose parameter, the bar on top means taking the average value (e.g. ....y is the total mean of all the measured values), and ),0(2N is a normally distributed random error. Average ki nematics determined from biplane images registered with 3D bead models were taken as reference values, and relative measurement errors were determined for the other cases. Measurement precision (the reciprocal of the variance) within each co mbination of conditions (or factors) was determined from the computed differences from reference values, and factor effects were estimated using a generalized linear model with the gamma distribution: sijk 2 offset i n planes 2 j bead or bone 2 observer k 2 (Eq. 22) ),(~/12ksijk 22 PAGE 23 where is the sample variance within each combination and each is the variance due to a specific factor. Statistical signific ance was defined as a proba bility (p) of 0.05 or less. Rootmeansquare errors due to each fa ctor were constructe d from the mean and estimated variance values. Based on clinic ian feedback, measurem ent errors greater than 1mm or 1 are considered noteworthy. 2 ijks2Results The actual ranges of joint motion measur ed from the 10 image pai rs of each joint are presented in Table 21. All joints dem onstrated joint rotations greater than 25 about one or more axes.except the talocalcaneal joint. Joint kinematics measurement accuracy wa s significantly worse for all singleplane measurements compared to all biplane measurements (Table 22). Five of eight joints showed differences of greater than 1mm in outofplane (Zaxis) translations comparing singleplane measur ements to biplane measurements. Inplane translation measures were significantly different from th e biplane measures in all eight joints, and in 5of8 cases these joints had significant outofplane translation differences. Only the talocalcaneal joint showed significant large (>1 ) differences in the zaxis rotation (plantar/dorsiflexion). In general, the best singleplane results were obtained for the knee and elbow joints, with the shoulder and foo t/ankle joints having larger differences from the reference measures. Model type, 3D bones versus metallic beads, had a very strong influence on measurement accuracy (Table 23). Seven out of eight joints registered with 3D bone models showed differences greater than 1 (differences > 2 in six of eight joints) in the outofplane (Xand Yaxis) rotations compared to the reference data using implanted 23 PAGE 24 metallic be ads. Inplane (Zaxis) rotation differ ences, where statistically significant, were less than 1.2 for all joints. Significant differences in translations were less than 1mm for all joints, except for the outofplane (Zaxi s) translation of the talocalcaneal joint. For most joints there were not large differences (>1mm or >1 ) due to the observer performing the measurement (Table 24). The exception to this finding was the talocalcaneal joint, where all six degrees of freedom showed significant differences of greater than 1mm or 1 Three of eight joints showed in terobserver differences greater than 1 for longaxis (Ya xis) rotations. Direct comparisons of specific datasets pr ovides additional useful information. We see sometimes large (> 2mm or >2 ) RMS differences between observers performing the 3D pose measurement for the bipl ane bead dataset (Table 25). The largest average differences were observed for the ta localcaneal joint and for the Yaxis (longaxis) rotations. The tibiofemoral joint (knee) was most consistently measured. Singleplane measures using the bead models showed large (> 2mm or >2 ) RMS differences for the outofplane translati ons (Z) and rotations (X and Y), and were largest for the talocalcaneal joint (Table 26). The tibiofemoral joint showed the smallest differences contrasting the singl eand biplane bead measurements. Biplane and singleplane measurements using 3D bon e models also showed large (> 2mm or >2 ) RMS differences for the outofpl ane translations (Z) and rotations (X and Y), and were largest for the talocalcaneal joint (Tables 27 and 28). The RMS differences for the biplane bone measurem ents were generally smaller than those for the singleplane bone measures for the outofplane transla tions (Z) and rotations (X and Y), and similar for the inplane translations (X and Y) and rotations (Z). 24 PAGE 25 Discussion Radiograph derived joint kinematic m easurements have become a common and important research tool for joint biomechan ics studies. Use of these techniques is rapidly expanding with the introduction of free opensource and commercial software. As the users of these tools sh ifts from the original developer s to those who simply want to use the tools to perform measurements, it is increasin gly impor tant to establish expectations for measurement quality acro ss the broad range of applications. In this study we sought to quantify joint kinemati cs measurement accuracy across a range of joints, imaging setups, and users. We confir med that singleplane measures are not as accurate as biplane measures, especially for outofplane translations. We showed that outofplane rotations, in particular, are less accurate using bone models for registration compared to clusters of small spherical metallic markers. Long bones showed larger errors for measuring rotations about the long axis, except the knee, where the articular morphology has a rich shape compared to the proximal radius or humerus. Finally, we showed that interobserver differences can be significant, especia lly for outofplane rotations and for specific challenging joints. Beyond the specific findings, these data for different combinations of image planes, object models and observers provide a baseline assessment of registration performance wit h which future developments can be compared. The study has several limitations that ar e important to recognize. First, we used only one software tool and empl oyed a global numerical optimizer to perform the registrations. We did not refine the global op timizer results with a gradientbased routine afterward. We believe these results are representative of a very basic manualplus25 PAGE 26 numerical registration procedure, and serve to highlight areas where further refinement of registrati on techniques and numerical algo rithms might best be targeted for specific joints. The results do not suggest that ot her reported registration techniques are necessarily inferior or superior. Second, the study images were generated using a biplane imaging system that pulsed the two xray tubes in an evenly alternating fashion. Since the images were not simultaneously acqui red there will be errors introduced in the biplane kinematic measurement. We sought to minimize this source of error by slowly moving the test specimens during image captur e. Finally, the upper extremity specimen was not attached to the thorax, so the images lacked the ribs and chest anatomy present in vivo. Tang et al. (Tang et al., 2004) assessed t he accuracy of measuring patellofemoral joint kinematics using a cadaver setup with singleand biplane imaging, and performed registration with bone models and beadbased models. When comparing the singleplane measures to the biplane measures using beadbased models, they found RMS differences of 1.6 for rotations, 1.6mm for inplane translations, and 1.3mm for outofplane translations. We found RMS differences of 1.2 for rotations, 0. 8mm for inplane translations, and 5.9mm for outofplane transla tions (Table 26). When comparing the singleplane measures to the biplane measures using projected 3D bone models, they found RMS differences of 2.1 for rotations, 1.8 mm for inplane translations, and 2.0mm for outofplane translations. We f ound RMS differences averaging 2.4 for rotations, 0.9mm for inplane translations, and 5.1mm for outofplane translations (Table 28). Our results are mostly comparable to Tang et al., except our outofplane translation errors are much larger. Tang et al. used approximately double t he number of metallic 26 PAGE 27 beads in each bone for the beadmodel registrations, and used digitalreconstructedradiograph (DRR) projections for the bonemodel registrations. These technica l differences likely account for Tangs gr eater outofplane translation accuracy. Bey et al. (2008) compared biplane meas urements of patellofemoral kinematics and found 0.4mm RMS differences for all translations and 0.9 differences for all rotations comparing beadand bonemodel based measures. We found average differences of 0.8mm for all translations and 2.2 for all rotations. Again, use of a DRR projection method for bonem odel registration appears to provide superior biplane bonemodel registration results. Morooka et al. (Morooka et al., 2007) reported an a ssessment of modelimage registration using synthetic images and found bestcase r egistration errors for the tibiofemoral joint of 0.5 mm for inplane translations, 2 mm for outofplane translations, and 0.5 for all rotations. Tsai et al. (Tsai et al., 2010) recently repor ted the development of an imageonly registration method for tracking the knee bones using fluoroscopy. They report experimentally measured accuracy of 0.8mm for inplane translations, 3.1mm for outofplane translations and 1.1 for all rotations. Comparing our singleplane bone model results to the reference data (Table 28) we found RMS differences of 0.7mm for inplane translations, 2.8mm fo r outofplane translations, and an average of 1.3 for all rotations; strikingly similar to Tsais results. These results confirm that even simple methods with user supplied initial guesses can provide clinically useful kinematic measures for 5 of the 6 degrees of freedom for the tibiofemoral joint. Our results confirm that additional work is requi red to improve the outofplane translation measurement, and to further autom ate the registration process. 27 PAGE 28 Bey et al. (2006) validated a biplane m odelbased registration technique for the shoulder (glenohumeral joint) in com parison to a beadbased biplane technique. Comparing bonemodel and beadbased registrations, they found average RMS errors of 0.385mm and 0.25 for the scapula, and 0.374mm and 0.47 for the humerus. Assuming these errors are uncorrelated, RMS errors for measuring scapulohumeral kinematics would be approx imately 0.8mm and 0.75. Comparing the biplane bone model results with the biplane bead m odel results, we found RMS differences averaging 0.8mm and 2.0 for scapulohumeral kinematics. Bey et al. used a similar number of beads in each bone, but used a DRR projection method for the bone model registrations. This more realistic model pr ojection method likely accounts for the higher rotational accuracy. There appear to be no reports of accu racy assessment for singleplane modelimage registration with the shoulder, elbow or hindfoot. Our data will provide basic observations on measurement performance wit h these joints. Measurement results for the elbow showed translation errors of less than 2mm, but large longaxis (Yaxis) rotation errors were observed. This finding should not be surprising as both the radius and ulna present as long cylindric objects with few shape features to change the projection with longaxis rotation (Fig 21). Singleplane studies will benefit from visualization of the entire forearm to the wrist, wher e the radius and ulna are not cylindrical. Measurement results for hi ndfoot kinematics (tibiotalar and talocalcaneal joints) showed typical translation errors in exce ss of 1mm and typical outofplane rotation errors greater than 3 (Table 22). These measure are poor compared to the knee joint, 28 PAGE 29 and are large relative to the range of join t rotations (Table 21) Future efforts to measure hindfoot kinematics with singleplane registrati on techniques will need to incorporate surface interactions and, per haps, ligamentlike constraints to improve measurement accuracy. There were statistically significant differences in the kinematics measured by different observers for every joint and for most degrees of freed om. The talocalcaneal joint was notable for these interobserver differences. These findings underscore the importance of developing and using more capabl e optimization routines and registration procedures to eliminate the e ffect of the initial pose esti mate on the optimized pose. This study provides a comprehensive a ssessment of modelimage registration performance for eight joints considering t he number of image planes, the registration model and interobserver variability. Thes e observations highlight where simple registration methods currently provide meas ures useful for research or clinical application, and point out many cases wher e technical improvements are required. These data will provide a useful reference for assessing the utility of new registration techniques for quantifying 3D bone kinematics. 29 PAGE 30 Table 21. Actual range of motion in t he 10 image pairs per joint, measured using biplane bead registrations. Joint X translation ( mm ) Y translation ( mm ) Z translation ( mm ) X Rotation ( de g) Y Rotation ( de g) Z rotation ( de g) Patellofemoral 65 31 2 82 10 6 Tibiofemoral 16 3 4 124 2 22 Talocalcaneal 3 6 4 6 4 4 Tibiotalar 2 1 1 25 7 4 R a d i o u l n a r21461 11 3 Radiohumeral 1 4 2 12 10 139 Elbow (ulna/humerus)2 1 4 6 11 136 Scapulohumeral 3 10 7 32 8 33 6 Table 22. RMS differences between bip lane and singleplane measurements for all observers, and bead and bone models. Joint X translation ( mm ) Y translation ( mm ) Z translation ( mm ) X Rotation ( de g) Y Rotation ( de g) Z rotation ( de g) Patellofemoral (n=3) 0.9*# 0.3# 3.7*Tibiofemoral (n=3) 0.6* 2.3*#0.3*0.4*Talocalcaneal (n=4) 3.8*# 3.4*# 4.2*3.3*1.2*Tibiotalar (n=4) 0.3* 1.0*# 3.3*0.6*Radioulnar (n=4) 1.4* 0.4* 1.4*0.9*5.7*0.8*Radiohumeral (n=4) 0.2* 0.3* 0.5*Elbow (ulna/humerus) (n=4) 1.4* 0.4* 1.4*0.9*5.7*0.8*Scapulohumeral (n=3) 0.9* 1.4* 10.1*#2.3*1.1** indicates a statistically significant difference in mean values compared to the reference measures. # indicates a statistically significant difference in sa mple variance compared to the reference measures. indicates a statistically significant differenc e with reference measurements was not found. Shaded cells highlight differences greater than 1mm or 1. Table 23. RMS differences between bead and bone models for all observers and singleand biplane measurements. Joint X translation ( mm ) Y translation ( mm ) Z translation ( mm ) X Rotation ( de g) Y Rotation ( de g) Z rotation ( de g) P a t e l l o f e m o r a l ( n = 3 )Tibiofemoral (n=3) 0.4#0.9*0.6* 1.2*# 0.5*#Talocalcaneal (n=4) 0.5#1.2*2.8# 2.5*#Tibiotalar (n=4)0.4*0.9*3.1*# 4.7*# 1.2*#Radioulnar (n=4) 0.9* 0.5#3.3* 0.3*R a d i o h u m e r a l ( n = 4 )1.3* 2.6#Elbow (ulna/humerus) (n=4) 0.9* 0.5#3.3* 0.3*Scapulohumeral (n=3) 0.2*1.1* 2.0*# 0.4*# 30 PAGE 31 Table 24. Average difference of each observer from the m ean of all observers, for all image data, including singleand bi plane measurements, and bead and bone models. Joint X translation ( mm ) Y translation ( mm ) Z translation ( mm ) X Rotation ( de g) Y Rotation ( de g) Z rotation ( de g) Patellofemoral (n=3) 0.1*0.6*0.5*0.3*Tibiofemoral (n=3) 0.2*0.1*0.8*0.2*0.2*Talocalcaneal (n=4) 1.7*3.0*1.5*2.2*2.8*1.2*Tibiotalar (n=4) 0.1*0.3*0.8*0.6*0.1*Radioulnar (n=4) 0.4*1.0*0.5*1.1*Radiohumeral (n=4) 0.2*0.2*0.5*0.8*0.1*Elbow (ulna/humerus) (n=4)0.4*1.0*0.5*1.1*Scapulohumeral (n=3) 0.1*0.2*2.0*0.3*0.2* Table 25. RMS difference bet ween observers for biplane r egistration of bead models. Joint X translation ( mm ) Y translation ( mm ) Z translation ( mm ) X Rotation ( de g) Y Rotation ( de g) Z rotation ( de g) Patellofemoral (n=3) 0.3 0.3 0.6 1.1 1.2 0.8 Tibiofemoral (n=3) 0.4 0.4 0.4 0.3 0.3 0.4 Talocalcaneal (n=4) 3.6 6.8 3.8 4.8 6.5 5.5 Tibiotalar (n=4) 0.7 0.3 0.6 0.9 1.5 1.0 Radioulnar (n=4) 1.2 0.5 1.7 0.5 3.6 0.5 Radiohumeral (n=4) 0.5 0.4 1.2 0.5 2.4 0.8 Elbow (ulna/humerus) (n=4)1.2 0.5 1.7 0.5 3.6 0.5 Scapulohumeral (n=3) 0.4 0.5 0.7 0.8 1.6 0.4 Table 26. RMS differences between bip lane bead model measurements and singleplane bead model measurements. Joint X translation ( mm ) Y translation ( mm ) Z translation ( mm ) X Rotation ( de g) Y Rotation ( de g) Z rotation ( de g) Patellofemoral (n=3) 1.1 0.5 5.9 1.1 1.3 1.3 Tibiofemoral (n=3) 0.6 0.9 3.9 0.5 0.4 0.6 Talocalcaneal (n=4) 4.3 3.4 5.8 3.5 3.6 3.0 Tibiotalar (n=4) 0.8 1.5 4.1 3.6 1.9 1.1 Radioulnar (n=4) 1.8 1.2 3.5 3.1 6.4 1.1 Radiohumeral (n=4) 0.7 0.5 3.3 2.8 3.1 0.8 Elbow (ulna/humerus) (n=4)1.8 1.2 3.5 3.1 6.4 1.1 Scapulohumeral (n=3) 1.7 2.1 11.0 2.0 1.5 0.5 31 PAGE 32 Table 27. RMS differences between bip lane bead model measurements and biplane bone model measurements. Joint X translation ( mm ) Y translation ( mm ) Z translation ( mm ) X Rotation ( de g) Y Rotation ( de g) Z rotation ( de g) Patellofemoral (n=3) 0.4 0.4 1.4 2.2 3.2 1.2 Tibiofemoral (n=3) 0.7 0.5 1.0 0.7 1.5 0.9 Talocalcaneal (n=4) 3.7 6.8 3.5 8.5 4.5 5.6 Tibiotalar (n=4) 0.8 0.6 1.4 4.1 4.8 1.7 Radioulnar (n=4) 1.6 1.0 3.3 1.6 7.7 0.7 Radiohumeral (n=4) 0.7 0.7 4.5 1.4 5.8 0.8 Elbow (ulna/humerus) (n=4)1.6 1.0 3.3 1.6 7.7 0.7 Scapulohumeral (n=3) 0.5 0.6 1.2 1.5 3.8 0.6 Table 28. RMS differences between bip lane bead model measurements and singleplane bone model measurements. Joint X translation ( mm ) Y translation ( mm ) Z translation ( mm ) X Rotation ( de g) Y Rotation ( de g) Z rotation ( de g) Patellofemoral (n=3) 1.3 0.6 5.1 3.3 2.6 1.3 Tibiofemoral (n=3) 0.6 0.7 2.8 1.1 1.9 0.8 Talocalcaneal (n=4) 4.3 3.4 3.2 7.3 7.6 2.9 Tibiotalar (n=4) 0.8 0.7 1.6 6.6 6.1 2.0 Radioulnar (n=4) 1.5 1.3 3.3 2.0 6.2 1.0 Radiohumeral (n=4) 0.7 0.7 2.0 3.1 3.7 1.0 Elbow (ulna/humerus) (n=4)1.5 1.3 3.3 2.0 6.2 1.0 Scapulohumeral (n=3) 1.1 1.9 13.3 4.1 8.4 2.0 32 PAGE 33 33 Figure 21. Biplane registra tion of elbow bone models. PAGE 34 CHA PTER 3 ARTICULAR CONSTRAINTS Introduction 3Dto2D modelimage registration per formed with radiographic images is an important tool for measuring invivo joint kinematics. It has been proven useful for the study of a number of c linical and research topics in a va riety of joints, such as implant performance (Banks et al., 2005), pathological or natural joint mechanics (Yamaguchi et al., 2009; Matsuki et al., 2010), and contact model development (Fregly et al., 2005). Measurement uncertainty usually is on the order of 1mm or 1 for 5 of 6 degrees of freedom characterizing joint kinematics from a singlepl ane projection (e.g. Banks and Hodge, 1996; Zuffi et al., 1999). However, the same level of accuracy is not achieved in the outofplane direction. Biplane or stereoradiographic system s achieve better and more uniform accuracy, but they limit the vi ewing volume in which an activity can be performed, double subject xray exposures, increase cost, and are not available at many institutions. Thus, there remains st rong incentive to improve the accuracy of singleplane measures for research and, ultima tely, clinical joint kinematic assessments. Hirokawa et al. (2008) incorporated an articular surface model into their registration method to impr ove outofplane measurement accuracy for total knee replacements. They added a system of ten cons traint equations in fourteen variables to describe the points of contact between t he tibia and femur and found numerically a solution that best satisfied their imageregistration costfuncti on and articular constraint equations. Validation experi ments with a robot manipul ating knee replacement components reduced outofplane translation e rrors from approximately 5mm to less than 1mm. 34 PAGE 35 Prins et al. (Prins et al., 2010) impr oved outofplane measurement accuracy by detecting collision bet ween the femoral component and the tibial insert in total knee replacements. In their phantom study eval uation, accuracy was improved from 2.1mm to 0.1mm RMS translation in the outofplane direction, with no loss of accuracy in the other degrees of freedom when using accurate 3D implant models In their method, pose estimation is first performed with the tr aditional imagematchingonly method, and then a correction step with the inclusion of penet ration detection is carri ed out to correct the outofplane translation. In this chapter, we present a different approach to detect and compute a measure of articular contact or penetration. We in corporate an articular surface proximity measure with several optimization algor ithms and compare the success rate and accuracy compared to a reference dataset. Compared to the methods of Hirokawa et al. or Prins et al., this method prov ides several potential advantages: Surface penetration can be quantified as dept h, area, quasivolume, or can be represented as an elastic foundation (Win kler, 1867). Surface separations also can be calculated within the same framework. Thus, a wide range of penalty functions or optimization constraints c an be configured from the same basic method. Because explicit surface penetration/separation distances are computed, physically intuitive gaps or deadbands can be incorporated to account for model errors. Thus, model inaccuracy, surface def ormation and wear are simply handled. The pose estimation process is carried out in a single run. An initialization step using a traditional imageonly method is not needed. 35 PAGE 36 Methods In the traditional method of 3Dto2D model image registration for joint kinematics measurement, the pose of the each bone or im plant is estimated separately. Given the 3D model of an object (a bone or an im plant) and the radiographic projection parameters, a computer simu lated perspective projection of the 3D model can be generated and compared with the actual xray images using an image similarity measure. Given a single fram e radiographic image and a 3D model, a search procedure, often a numerical optimization routine, is then used to find the best set of 3D model pose parameters that maximize the similarity measure. This procedure is repeated for each bone or implant and for each frame to obtain the complete temporal kinematics of the bones (Figure 31). Similar to Prins et al. (2010), we combi ne a penetration computation routine with a numerical optimizatio n routine and simultaneously opt imize the image similarity measure and surface penetration. The 12 pose parameters of t he articulating bones or implants are optimiz ed together rather than independently. Penetration Computation Algorithm Triangle mesh surface models are used due to their si mplicity, universal availability and errorfree exchangeability between differ ent 3D modeling software packages. Natural bone models are typically obtained by CT or MRI scans and in the form of 3D images. The 3D images can be segmented and transformed into triangle meshes. Implant models are usually 3D CAD models and can be easily exported to triangle mesh models. Surface penetration/separation is com puted based on a pointwise closest distance search using a kdtree data st ructure (Bentley, 1975). For each point on one 36 PAGE 37 surface, a query is conducted to determine t he closest point on the other surface. The relative orientation of the two triangles and the outwardfacing normal determines penetration or separati on and the distance co ntribution to the final penetration score (Fig. 32). The kdtree data st ructure is a space subdivision algorithm specialized for spatial search and proximity query (Bentley, 1975). Given a swarm of spatial points, a single search for the point that is nearest/furthest to a given coordinate has a computational complexity of O(log(n)). For a pair of potentially contacting models in our application, the kdtree needs to be construc ted only once for only one of meshes. If one surface is designated the treemodel and the other the nontreemodel, the pseudocode of the algorithm follows: Preprocessing step: Tree_model_centroids = Centroids( triangles_in_the_tree_model ) KDTree = Construct_KDTree( Tree_model_centroids ) in the treemodels local coordinate system. Contact computation within each iteration of the optimization problem (when estimated poses of the two models are given): Transform the nontreemodel into th e treemodels coordinate system. FOR each triangle T in the nontree model FOR every vertex P of triangle T C_tree = KDTree_FindClosestPoint( P ) T_tree = Triangle_corresponding_to C_tree n = Outward_Normal( T_tree ) v = P C_tree (i.e. the vector connecting C_tree to P) IF Dot_product(n, v) < 0 Conclude P is inside the surface of the treemodel (Fig. 32, P1) ELSE Conclude P is outside the surf ace of the treemodel (Fig. 32, P2) ENDIF ENDFOR Calculate penetration score contribution from T based on the location of the vertices. ENDFOR 37 PAGE 38 Penetration depth, area, quasivolume or separation distance are readily computed using this approach. On a PC workstat ion with Intel Xeon 2.4GHz CPU and Windows XP SP3, the contact computation function takes roughly 0.2s to complete with a pair of 3D models of 5000 triangles each. Image Comparison Metric The image similarity measure is a measure of distances between the object contour in the actual fluoro scopic image and the cont our of the synthet ically projected object. First, edge pixels are identified us ing a Canny edge detec tor (Canny, 1986) in the actual fluoroscopic image and in the synthetic image. An image mask is created by dilating the edge image of the synthetic image by 25 pixels. The mask is then applied to the edge image of the actual xray image to keep only the nearby edge pixels. For each remaining edge point, the closest (in the sens e of image coordinates) contour pixel in the artificial projection image is identified and the distance recorded. A second closest point search is performed to find the distance of this point to its nearest edge pixel in the edge image of the actual xray image. This distance forms the basis of the contour similarity measure. The final image similarity measure is the sum of squares of all such distance values. Optimizers and Cost Functions The surface penetration measure can be incorporated into the registration procedure by adding the penetration score as a penalty to the image similarity measure or as an optimization constraint. Three representative combinations of optimization routines and registration metrics were eval uated. No preoptimization step was needed for any of the three methods. 38 PAGE 39 Particle sw arm pattern search (Pswarm) : Particle swarm pattern search (Vaz and Vicente, 2007) combines a stochastic particle swarm global optimization routine with the search/poll structure of local pa ttern search methods. The algorithm only supports boxbound constraints, so the articular surface penetration function only could be integrated as an augmented penalty function with the image similarity metric. The scalar penetration measure was squared and added to the image similarity measure to form the final cost function for optimization: i ixdxp xCost2 2)()()( (Eq. 31) where x is the 12 pose parameters of the two articulating objects, is the penetration measure, is the image comparison me tric as described earlier (sum of squares of edge pixel dis t ances), and the weighting factor )( xp2)( xdi was determined empirically (50 when p is penetration depth). The swarm si ze was set to 20, and default values for other parameters were used (cogniti al parameter = social parameter = 0.5, initial inertia = 0.9, final inertia = 0.4). RALG optimizer (OpenOpt, version 0.29, National Academy of Sciences of Ukraine, http://openopt.org/ralg ): RALG is a nonlinear/ nonsmooth problem solver based on Shors ralgorithm (Shor, 1985), and is a gradientbased deterministic local optimizer. A penalty approach exactly the same as used for the particle swarm optimization was tested. In a second approach, the penetrati on function was integrated into the optimization as a user supplied cons traint function. When supplying the penetration function as a constraint, it was found that scaling in the magnitude of the constraint function affected success rate, and the function value was scaled with a weight of 1e6. To provide the first derivativ es of the cost function and the constraint 39 PAGE 40 function to the optimizer, numerical forward differentiation was performed with a differentiation step size of 0.3mm or 0.3. Nonlinear least squares optimization (LSQ) : The LevenbergMarquardt nonlinear least squares algorithm takes as input a vector of errors. As a result, the penetration function could be added as one or more additional elements in the error vector. This is a penalty approach. The vector of image comparison errors was taken as the original pointwise cont our distances before taking sum of squares. The penetration penalty was appended to the error vector with a weight such that taking the sum of squares would result in the exact cost function used with the ot her two optimizers. Numerical forward differentiation was carried out with a differentiation step size of 0.3mm or 0.3 to generate the requi red Jacobian matrix. A full list of the optimizer parameter values used can be found in Appendix A. Experimental Setup Thirty biplane image pairs, reverse engi neered 3D knee implant models, and biplane kinematic measurements were obtained from P rins et al. (2010). The same dataset was used to evaluate our methods and their biplane measurements were taken as the reference measurements. Singleplane registration was performed using each of the three methods as described earlier, with and withou t the surface penetration function. A standard set of initial poses was generated by adding uniform random per turbations of 2mm for inplane translations, 3mm for outofplane translations, 3 to the reference measurements. This set of perturbed poses wa s used as the initial guess for all singleplane registrations. Modelimage registra tion was performed for each model in each 40 PAGE 41 frame. No correction or repeat registra tion was performed on any of the matching results. For the deterministic RALG and l eastsquares optimizers, one single run was performed. For the stochastic particle swarm pattern search, 10 run s were performed using different random seeds for the initial particles. Registration algorithm performance wa s quantified using success rate and accuracy. The success rate was defined as the percentage of successful registrations on a framebyframe basis. A frame qualif ied as an unsuccessful registration if, compared to the reference m easurements, the estimated abs olute pose of either the femur or tibia implant model had an error of larger than 1 mm of inplane translation, 2 in X and Z rotation, or 3 in Y rotation. (For the tibi a component, Y rotation corresponds to the rotation about the long ax is of the tibia, which is more difficult to measure accurately.) Accuracy was parameterized as mean, standard deviation and rootmeansquare (RMS) errors in the measured joint kinematics compared to the reference kinematics. When computing the error statis tics, only poses in successfully registered frames were included in the population. The RALG and particle swarm optimizers were allowed to run for 1000 cost function evaluations. The nonlinear least s quares optimizer typically converged within 500 cost function evaluations. In the imageonly registrations without contact functions, the femur and tibia components were regist ered separately. To allow for a fair comparison with the contactenabled methods, maximum iteration limits were set at 500 cost function evaluations, and in the case of the particle swarm optimizer, swarm size was also reduced by half to 10. 41 PAGE 42 All registrations were performed using custom opensource software (JointTrack, www.sourceforge.net/projects/jointtrack ). The penetration computation algorithm was implement ed as a JointTrack component in C++ using a third party kdtree implementation (CGAL, Computational Geometry Algorithms Library, http://www.cgal.org/). Regist ration routines and other al gorithms were implemented or connected to JointTrack in either C++ or Py thon, using available third party libraries when possible. Results Inclusion of surface penetration penalties or constraints improved the success rate and/or the outofplane translation accuracy for all three methods (Table 31, Figs. 33 and 34). Inclusion of surface interactions had no or little effect on the success rate for the nonlinear leastsquares method, but re sulted in 7% improvement for the RALG method. Significant improvements in outofplane translation precision were observed for each of the three methods when using a penetration function (Ftest, p < 1e10 for all three methods). Xaxis translations (inpl ane) were halved when contact constraints were included. Inclusion of surface penetrati on had little effect on the other four degrees of freedom, changing RMS errors by a maximum of 0.1mm or 0.2. Discussion Modelimage registration techniques prov ide important kinem atic measurements for joint biomechanics studies. These tec hniques often are applied f or singleplane radiographic image sequences, w here the measurement uncerta inties are much greater for outofplane translations. When joint surface geometries are known, it is possible to include additional information or constraints to the modelimage registration process to improve the outofplane m easurement accuracy. We showed that modelimage 42 PAGE 43 registration with quantified penetration resulted in significantly improved accuracy in the outofplane translation compared to imageonly registration. Inplane translation accuracy improved as a direct conseque nc e of improved outofplane accuracy. Accuracy in the rotational degrees of freedom was comparable between contactenhanced methods and the standard methods. Si milar accuracy was observed for the three numerical methods evaluated. Evaluation of the surf aceenhanced modelimage registration techniques was carried out using phantom images and revers eengineered implant CAD models. These conditions may not capture all of the important feat ures of in vivo images, and therefore may not give a true indication of the utility of these new methods in vivo. However, these new methods were evaluated using t he same image dataset and implant models as Prins et al. (2010), so a rigorous quantitative comparison of these related techniques is possible. Hirokawas method (Hirokawa et al., 2008) a ttempt to satisfy a contact constraint for either or both condyles, but does not provide directly a measure of surface separation that would be useful for uncertain or elastic surfaces. The collision detection algorithm used in Prins et al. (2010) is a type of boundi ng box tree algorithm. This Boolean collision detection algorithm quickly finds intersecting lines between two surfaces but requires additio nal special data structures and algorithms to obtain the area or volume of penetration. Our cont act algorithm, based on the kdtree data structure, permits fast co mputation of penetration area and quasivolume as well as surface separation distances, which c ould be used to implement ligamentlike separation constraints in the future. Hirokawa et al. report ed a reduction in outofplane 43 PAGE 44 translation errors from approxim ately 5mm to 1mm usi ng contac t constraints. Prins et al. performed a twostep registration procedure and reported a reduction in outofplane translation errors from 2.1mm RMS to 0.1m m RMS. Our singlestep registration results with and without the contact c onstraint gave outofplane regi stration errors averaging 2.5mm and 0.5mm RMS, respectively, comparable to the previous reports. The RALG method performed best in terms of success rate of the three methods incorporating penetration detection. The cons traintfunctionbased RALG optimization resulted in performance similar to penaltybas ed RALG optimization, but required trialanderror scaling of the constraintfuncti on to achieve good performance. Scaling sensitivity could have resulted from initial guesses being outside the feasible region defined by the constraints. Constrained optimization algorithms often require initial guesses to be within the feasible region to func tion correctly, but such a condition is not typically guaranteed in our application. The nonlinear least squares method wit h edgebased image metric provided the worst success rate but was insensitive to parameter changes, especially changes to contact penalty weighting. Appending 10,000 r eplicates of the co mputed penetration to the error vector did not produce a different success rate than appending 100 replicates. In addition, the least squares algorithm does not use or require a boxbound constraint on the search region. The other two algorithms, especially the particle swarm algorithm, demonstrated significantly lower success ra tes when the boxbound constraint was relaxed. The least squares method also converged quickly, typically in less than 300 function evaluations, while the other two met hods required at least twice the number of function evaluations to achieve stable success rate performance. 44 PAGE 45 The particle swarm pattern search al gorithm, a stochastic global optimizer, achieved a better success rate than the det erministic nonlinear least squares method. Howev er, the particle swarm algorithm is very sensitive to contact penalty weighting and to the bounds supplied for the search regi on. The reported success rate was the bestcase value, and the success rate often fell to lower than 50% when different sets of parameters were used. For bestcase success rate, the boxbound constraint was set to exactly the same bounds as were used to generate the perturbed st arting positions. It also should be noted that resu lts from global optimizers typi cally are refined with a local gradientbased optimizer upon completion. We did not include this refinement step, suggesting the results with particle swarm pa ttern search could be further enhanced. The present data set utilized a posteriors tabilized implant with a central femoral cam and tibial post, and the choice of cont act measure did not produce significantly different results. However, on a posteriorcruciateretaining implant design lacking a camandpost type constraint, re gistration using the quasivolu me measure significantly outperformed the depth measure in terms of success rate. This will be an important consideration for use of penetration penalties for registration of natural joints, where there is wide variation in articular congrui ty and lower accuracy in defining the joint surface models. It is reasonable to speculat e similar registration penalty functions could be useful with anatomic joints, but additiona l evaluations need to be conducted to prove this conjecture. Acknowledgements The authors thank Andre Prins and Bart Kaptein for generously providing the biplane dat aset and reference measurements. 45 PAGE 46 Table 31. Success rates and RMS registra tion errors for three methods with and without the contact function penalt y or constraint. Success tx ty tz rx ry rz Method Condition rate (mm)(mm)(mm)(deg)(deg) (deg) no contact 80% 0.7 0.2 2.1 0.4 1.4 0.7 LSQ contact 80% 0.3 0.2 0.6 0.5 1.6 0.7 no contact 88% 0.8 0.2 2.4 0.2 0.6 0.6 Pswarm contact 89% 0.4 0.3 0.5 0.2 0.5 0.8 no contact 90% 0.9 0.2 3.1 0.3 0.8 0.5 RALG penalty 97% 0.4 0.2 0.3 0.2 0.9 0.6 constraint 97% 0.4 0.2 0.6 0.3 0.8 0.7 tx and ty represent inplane translations, tz represents outofplane translations; rx, ry and rz represent rotations about the inplane and outofplane axes, respectively. Registration in a frame is unsuccessful if the estimated absolute pose of either the femur or tibia implant model had an error of larger than 1 mm of inplane translation, 2 in X and Z rotation, or 3 in Y rotation. 46 PAGE 47 Figure 31. Knee implant models matched to the fluoroscopic image to measure their spatial pos es. Figure 32. Illustration of inside/outside determination for ar ticulating polygonal surfaces. For clarity, the 3D problem is shown as a simplified 2D diagram, where line segments are used to represent triangl es and their midpoints are used to represent centroids of triangles. Nearest triangle Outward normal Nearest triangle Outward normal (kdtree of Treemodel triangle centroids) P1 P2 Nontreemodel 47 PAGE 48 48 Error in joint position and joint angles: standard deviation0.0 1.0 2.0 3.0 xyzx_rotationy_rotationz_rotation MeasureError SD (mm or degree) LSQ: traditional LSQ: contact Pswarm: traditional Pswarm: contact RALG: traditional RALG: contact penalty RALG: contact constraint Figure 33. Standard deviation of the errors in the measured relative joint position and joint angles Error in joint position and joint angles: mean1.0 0.0 1.0 2.0 3.0 xyzx_rotationy_rotationz_rotation MeasureMean Error (mm or degree) LSQ: traditional LSQ: contact Pswarm: traditional Pswarm: contact RALG: traditional RALG: contact penalty RALG: contact constraint Figure 34. Mean errors in the measured relative joint position and joint angles. PAGE 49 CHA PTER 4 LIGAMENTLIKE JOINT CONSTRAINTS Introduction 3Dto2D modelimage registration per formed with radiographic images is an important tool for measuring invivo joint kinematics. It has been proven useful for the study of a number of c linical and research topics in a variety of joints, including knee (Fregly et al., 2005), ankle (Yamaguchi et al., 2009), hip (Tanino et al., 2008), elbow (Matsuki et al., 2010), and shoulder joints (Nishinaka et al., 2008). In a typical setup, radiographic images of a patients joint are continuous ly taken during a motion trial. Poses of the bone or implants in each frame are then measured from the images using the method of 3Dto2D modelimage registration: Given the 3D model of an object (a bone or an implant) and the radiographic projection parameters, a computer simulated perspective projecti on of the 3D model can be generated and compared with the actual xray images usin g an image similarity measure. Given a single frame radiographic image and a 3D model, a search pr ocedure, often a numerical optimization routine, is then used to find the best set of 3D model pose parameters that minimize (or maximize) this similarity measure. The complete temporal kinematics of the bone is obtained by performing this proc edure for all image frames in the trial (Figure 41, Right). A phenomenon often present in the single pl ane 2D3D registration task is the symmetric view problem (Figure 41). In such situations, a choice must be made according to the objects relationship with other nonsymmetric objects whose orientations are much easier to determine. This represents an extreme case of the local minima problem where there are two equa lly good choices based on the registration 49 PAGE 50 criteria. The general registration problem is not simple, as typical costfunctions exhibit many local minima and/or long narrow valleys, which depend on the shape and pose of the object, the image similari ty measure, features in the image and image quality. Optimization procedures generally are not able to find the correct solu tion without a very good initial alignment provided by a human operator. W hen a symmetric view problem is present, even human operators sometimes fail to give correct initial guesses. Another problem intrinsic to singleplan e monocular measurements is the high uncertainty in outofplane translation. Hirokawa et al. (2008), Prins et al. (2010) and our group have introduced methods to incor porate articular surface proximity measurements to the traditional imageonly registration methods, and achieved improved outofplane measurement accuracy. These methods, however, are only effective when a relatively good initial guess is given or after a preliminary registration step. Human joints are not simply bones. Soft tissues, especially ligaments, play an important role in constraining the motion of a joint and provide a useful framework for additional registration constraints. In th is study, we exploit this knowledge and incorporate ligamentlike constraints in t he registration algorithm and characterize registration performance improvements compared to techniques without such constraints. Methods Data Sets Two total knee arthroplasty data sets were used to evaluate the performance of our methods. The first data set (D ataset A) consists of 47 singleplane fluoroscopic images of a patients knee in a stairup motion. The second data set (Dataset B) is 18 50 PAGE 51 images in a walking motion. The total knee implants in the two datasets were cemented posterior cruciate ligament re taining prostheses of differ ent designs from different manufacturers, and their corresponding 3D CAD models were obtained from the manufacturers. The joint kinematics in both datasets were previous ly measured using custom software by experienced human operators, using a traditi onal imagesimilarityonly registration method. Signific ant operator supervision was required for the measurement process, and measurement error in the outo fplane direction was visually corrected by the operator using side views of the 3D models. Results from dataset A were previously used for joint contact and wear modeling (Fr egly et al., 2005), and results from dataset B were crossvalidated with readings from force sensors instrumented in the tibial component of the knee implant (Zhao et al., 2007). These kinematic measurements were used as the reference kinemati cs for evaluation of new methods. Registration Method The new method minimized a combination of image similarity measure, surface penetration, and ligament force. The surfac e penetration measure and ligament force were incor porated into the registration proc edure as penalties to the image similarity measure. The twelve pose parameters of the femoral and tibial components were optimized together rather t han independently for each object. Surface penetration: Penetration volume (mm3) between the femur component and the tibial polyethylene insert was estimated using the method described in Chapter 3. For dataset A, the retriev ed tibial insert showed creep on its contact surface, so a dead band of 1mm depth was specified in co mputing the penetration. Penetration 51 PAGE 52 volume to the 0.7 power was used as the penalty score to reduce the cost functions sensitivity to deep penetration. Ligament force : Ligament forces limit the knee joints range of motion. We first defined bounds that represent the normal range of moti on of the knee joint. Ligament forces are modeled as linear or torsional springs with linear forcedisplacement relationships. The springs are active onl y when the estimated poses produce joint translations or joint rotations exceeding the predefined bounds. The bounds are defined as 30 mm femoral posterior translation to 10 mm anterior translation 5 to 5 mm in mediallateral translation, 20 to 20 in internalexternal rotation, 3 to 3 in varusvalgus angle, and 10 hyperextension to 150 flexion. For dataset A, no active bounds are defined in the relative superiorinferior position between femur and tibia. For dataset B measured forces showed both femoral condy les always to be in contact, so a superiorinferior virtual spring was constant ly loaded, pulling the femur and the tibia together. Image metric : The scalar image similarity measure was the sum of normalized crosscorrelations between grayscale and edge images (Canny, 1986), comparing the actual fluoroscopic image and the computer synthesized radiographic projection. The optimization routine we used is t he RALG optimizer (OpenOpt, version 0.29, National Academy of Sciences of Ukraine, http://openopt.org/ralg ). RALG is a nonlinear / nonsmooth problem solver based on Shors ralgorithm (Shor, 1985). To provide the first order partial derivatives of the cost function to the optimizer, numerical forward differentiation was performed with a different iation step size of 0.3mm for inplane translations, 0.9mm for outofplane translation, 0.6 for y rotation, and 0.3 for x and z 52 PAGE 53 rotations. The differences in these step si zes resulted from the scaling of pose parameters to adjust the cost functions sen sitivity in different directions. For Dataset A, the 12 degrees of fr eedom being optimized consisted of the absolute poses of the femur and tibia component s. A single image similarity score was computed in each cost function evaluatio n using a computed projection of both components. In the combined cost function the image metric has unit weighting, the ligament force penalties are summed across a ll directions with a penaltydisplacement relation of 0.002 mm1 or 0.002 deg1, and the penetration penalty has a weight of 0.0002. The image metric has a t heoretical range of 1 to +1, but typically takes a value between 0.3 and 0 when a fair initial guess is provided. Dataset B was more challenging due to the presence of the transmitting and receiving antennae being placed at the anterior tibia in direct apposition to the projected tibial implant. For Dataset B, the absolute pose of the fe moral component and the tibial components pose relative to the femoral component were optimized. A femoral image score and a tibial image score were combined to form a total image score with a relative weighting of 2:1. Compared to Dataset A, weights of the penalty functions were increased by a factor of 10, i.e. 0.02 mm1 or 0.02 deg1 for ligament penalties and 0.002 for the penetration penalty. In addition, the registration pr ocess was divided into two stages for Dataset B: the first stage did not include in the cost function the surface penetration measure or the loaded pseudoligament in the prox imaldistal direction; a second stage optimization started from the results of the fi rst stage, and used all the surface penetration and ligamentlike penalties. 53 PAGE 54 Performance Evaluation Singleplane registration was performed using three methods: imagemetriconly registration, registration with image me tric and surface penetration penalty, and registration with image metric surface penetration metric, and ligamentlike force metric. Three sets of initial poses were generated by adding unifo rm random perturbations to the reference measurements. These sets of perturbed poses were used as initial guesses for all registrations. Due to the shape of total knee implants, the femoral component pose is typically easier and more accurately measured than the tibial component pose; it is also easier for a human operator to align the femoral mo del with the actual image to give initial guesses to the optimization proc edure. Therefore, for the femoral component, the size of the perturbations was set at 3mm for inpl ane translations, 20mm for outofplane translations, and 5 for rotations. For the tibia com ponent, only inplane translations and inplane rotation are easy to determine, so the perturbation size was 3mm for inplane translations, 50mm for outofplane translation, 5 for inplane rotation, and 10 for outofplane rotations. Modelimage registration using each met hod was performed for each model in each frame using three sets of initial guesses. No correction or re peat registration was performed on any of the matching results. T he optimization routine was allowed to run for 70 iterations (in case of Dataset B, 40 and 30 iterations for the two stages respectively), which was roughly 1200 cost function evaluations including evaluations to compute the derivatives. In the image only registrations without contact or ligament 54 PAGE 55 functions, the femur and tibia implants we re registered separately, and maximum iteration lim its were set at 35 evaluations. Registration algorithm performance was quantified using success rate, which was defined as the percentage of successful regi strations on a framebyframe basis. A frame qualified as an unsuccessful registration if, comparing to the reference measurements, the estimated abs olute pose of either the femur or tibia implant model had a deviation of larger than 1.5mm of inplane translation, 3 in X and Z rotation, or 5 in Y rotation. For the tibia component, Y rotation corresponds to the rotation about the long axis of the tibia, which is typically very difficult to measure accurately (Figure 42). Relative joint kinematics were com puted from the measured femur and tibia poses, and their deviations from the va lues computed from the reference measurements were quantified using r ootmeansquare (RMS) values. Only measurements from successfully registered frames were included in the population. Results Inclusion of the ligamentlike force penal ty and articular surface penetration metrics resulted in 24% and 46% improvement on the success rate for Dataset A and Dataset B, respectively, compared to imageonly registration (Table 41). Inclusion of the ligamentlike force penalty and articular su rface penetration metrics resulted in 31% and 70% improvement on the success rate for Dataset A and Dataset B, respectively, compared to imagepluscontac t registration. Compared to the results of the new method, errors in the relative mediallater al position (z position) were an order of magnitude larger for the traditional imageonly registration. Large deviations in the x and y positions resulted from large errors in t he outofplane translation measurements, due 55 PAGE 56 to the fact that the knee join ts mediallateral axis was not perfectly aligned wit h the direction of the xray projection in either dataset. Discussion Modelimage registration techniques prov ide important kinem atic measurements for joint biomechanics studies. These tec hniques have the common issue of reduced success rate when the initial gues ses are not carefully given, or when the images are of bad quality lacking distinct features. We showed that modelimage registration augmented with quantified ligamentlike forc es and articular surface penetration resulted in significantly improved success rate compared to imageonly registration. We used two datasets with distinct registration challenges, and showed significant improvements in registration performance with inclusion of ligamentlike force penalties. Images in Dataset A have high contrast and the cement was not visible (Figure 41 and 43). However, there are certain properties of the tibia components that sometimes trap the optimization into a r egion far from the real soluti on (e.g. erroneous results in Figure 43). The new method was able to increas e the rate of successful registrations. Ligamentlike force penalties were not active when the estimated poses were within the natural limits of the knee join t, and hence would not affect the accuracy of successful measurements. Their primary role was to guide the solution search procedure into the region where other techniqu es, especially the surface contact augmented techniques, effectively function and produce accurate measurements. In Dataset B the image quality was poor, th e external antennae partially block the tibial component projection, the image is low re solution (which creates difficulties for the edge detection algorithm), and an edge line just below the top part of the true tibial component silhouette creates a local minimum in the image similarity metric very close 56 PAGE 57 to the true solution (Figure 44) The main objective in this case was to preve nt the tibia model from getting trapped in an y of the nearby strong local minima. In the first stage optimization without surfacepenetration penalties, the ligament like penalties helped to avoid local minima relatively far from the true solution. In the second stage, a virtual ligament force was added to pull the tibia m odel against the femu r implant and escape spurious local minima. To prevent this fo rce from erroneously c orrecting the femur pose rather than correcting the tibia pose, t he tibial components relative pose to the femur rather than its absolute pose was optimized. By doing so, a change in only the femurs pose parameters does not cause a change in the relative pose between the femur and tibia, and therefore does not in crease the surface penet ration or ligament force penalties. This transformation essentially shapes the cost function and its numerically estimated partial derivatives in favor of the direction search mechanism of the optimizer. However, a side effect is the pose parameters for the femur implant become harder for the optimizer to esti mate because any change would affect the perspective projection of both the femur and tibia models. The image matching score due to the femur model projection were ther efore given heavier we ight than the tibia projection to alleviate this issue. Compared to Dataset A, a much hi gher weight from the penalty functions was also used. By utilizing these strategies, we significantly reduced the chance of the tibial model get ting stuck in the local minimum. Conceptually, one may argue surface penetra tion alone could also correct these problems with local minima. However, we showed that the surface penetration measure actually degraded performanc e without a very good initial guess (Table 41, 57 PAGE 58 image+contact cases). Although joint limit constraints over lap wit h contact constraints to a certain degree, each has s eparate and complementary roles. This study is the first step to explore different forcetype penalties for modelimage registration. Many potential alternatives have not yet been explored. For example, we did not evaluate the ligament force constraint without the surface penalty measure. This will be an important alternative, because it is much easier to specify motion constraints with ligamentlike penalties than it is to measure or obtai n accurate surface geometry to enable surface penetration constraint s. We did not evaluate diffe rent stiffness values or varying stiffness values for each degree of freedom. A sensitivity study would provide guidance as to how best to use ligamentlike constraints for specific types of imaging studies (i.e. knees will likely be different from shoulders). We have shown in this study that a ligament like constraint model could be utilized to improve accuracy and the registration su ccess rate when crude initial guesses are given, and to significantly reduce erroneous local solutions when using a carefully chosen optimization strategy. The ultimate goal of all these efforts is to provide a measurement technique that requires minimal user intervention. Inclusion of ligamentlike registration penalties significantly r educes the need for accurate initial pose estimates, but does not eliminate the need for human assistance or supervision for robust modelimage registration for skeletal kinematic measurements. 58 PAGE 59 Table 41. Success rates and the relative pose RMS deviations for two datasets with and without contact and/or ligament penalty. Position (mm) Orientation ( ) Dataset Method Success rate x y z x y z image 55% 4.7 1.6 15.9 0.9 2.4 1.2 A image+contact 48% 4.8 1.6 16.3 0.8 2.3 1.2 image+contact+ligament 79% 0.8 0.4 1.7 0.7 1.3 1.0 image 28% 1.0 1.1 12.6 1.3 2.8 1.1 B image+contact 4% 2.4 1.8 29.5 2.1 3.9 0.7 image+contact+ligament 74% 1.1 0.4 3.8 0.4 1.9 1.7 59 PAGE 60 Figure 41. The symmetric view problem in v olves nearly equivalent projections of symmetric objects. Figure 42. Coordinate system definition fo r the femoral and tibial implant models. 60 PAGE 61 61 Figure 43. Silhouette of illregistered tibia component models. The right side error is a case of semisymmetrical view. Figure 44. An image from Dataset B. Left: Edges identified by a Canny edge detector were superimposed on top of the im age. Right: Silhouet te of the tibia component model whose top contour wa s erroneously registered with the wrong edge. PAGE 62 CHA PTER 5 CONCLUSION 3Dto2D modelimage registration is increa singly used by researchers worldwide for skeletal kinematic measurement. The met hod has proven utilit y in many situations where alternatives for quality kinematics are few. However, there remain a number of technical issues requiring improvement, and t he key obstacle to wider adoption remains the significant human intervention required to pr ovide initial guesses for registration. We investigated several methods to address these issues, and demonstrated better accuracy and success rates. Typical modelimage registration techniques estimate the pose of each object separately. We showed significant improvement by estimating all objects poses together using additional domainspecific knowledge. Increasing the dimensionality of the problem may seem counterintuitive, but in doing so we divided the problem into a coarse largeerror rejection problem and a fine precision measurement problem, which makes better use of available information. We showed that an image similarity measure alone is useful for all conditions but insu fficient to robustly solve any of them. Future work should address how to obtain initial pose estimates when no humanprovided initial guess is given. Human beings learn, identify and reproduce the rough position and orientation of objects easily and quickly. A machine learning approach might be appropriate for this purpose. In fact, a former lab member performed preliminary work with the SIFT (Scale In variant Feature Trans form) classifier and obtained some promising result s. Whether this approach pr ovides the required accuracy in all directions for our existing methods to function effectively remains unknown. 62 PAGE 63 63 We showed the success rate with any method is less than 100%. Further work is necessary to address this issue as operator verification and correction still requires a significant amount of human inte rvention. As shown in Chapt ers 2 and 4, image quality, image modality, and object geometric proper ties have a dramatic effect on measurement accuracy and success rate. Advanced image analysis techniques and utilization of artificial digitally re constructed radiographs (DRRs) are needed. The new methods we presented here prov ide improved accuracy and robustness, making 3Dto2D modelimage registration an increasingly superior technique compared to traditional 2D im age measurements or other invasive or skinmarkerbased techniques for skeletal kinematic measurem ents. These new techniques and tools are now readily available for use by re searchers in our Orthopedic Biomechanics Laboratory at the University of Florida, as well as researchers worldwide through the free distribution of the JointTrack program. PAGE 64 APPENDIX OPTIMIZER PARAMETER VALUES Table A1. Particle swarm patte rn search parameter values. Parameter Description Value s swarm size 20 mu cognitial parameter 0.5 nu social parameter 0.5 maxvfactor maximum velocity factor 0.5 Maxiter maximum number of iterations 50 Maxf maximum of function evaluations 1000 iweight initial inertial weight 0.9 fweight final inertial weight 0.4 tol x tolerance for stopping criteria 1.0E05 idelta pattern search mesh size expansion factor 2 ddelta pattern search mesh size reduction factor 0.5 Table A2. RALG optimizer parameter values. Parameter Value alpha 2.0 h0 1.0 nh 3 q1 0.9 q2 1.1 hmult 0.5 S 0 dilationType plain difference doBackwardSearch True xTolerance 0.01 For the meaning of the parameter s, please refer to Shor (1985) and Kappel and Kuntsevich (2000). Table A3. LevenbergMarquardt nonlinear least squares opt imizer parameter values. Parameter Description Value ftol relative error desired in the sum of squares 1.5e08 xtol relative error desired in the solution 1.5e08 gtol orthogonality desired between the function vector and the columns of the Jacobian 0.0 maxfev maximum of function evaluations 1000 64 PAGE 65 LIST OF REFERENCES Banks, S., Hodge, W., 1992. 3D knee pros thesis kinematics from single plane fluoroscopy. Orthopaedic Transactions 16(2), 530. Banks, S.A., Fregly, B.J., Boniforti, F., Reinschmidt, C., Romagnoli, S., 2005. Comparing in vivo kinemati cs of unicondylar and biuni condylar knee replacements. Knee Surgery, Sports Traumatology, Arthro scopy 13(7), 551556. Banks, S.A., Hodge, W.A., 1996. Accurate measurement of thr eedimensional knee replacement kinematics using singlepl ane fluoroscopy. IEEE Transactions on Biomedical Engineer ing 43(6), 638649. Banks, S.A., Markovich, G.D., Hodge, W.A. 1997. In vivo kinematics of cruciateretaining and substituting knee arthroplasties The Journal of Arthroplasty 12(3), 297304. Banks, S.A., Riley, P.O., Spector, C., Hodge W.A., 1991. In vivo bearing motion with meniscal bearing TKR. Ort hopaedic Transactions 15(2), 544. Bentley, J.L., 1975. Multid imensional binary search trees used for associative searching. Communications of the ACM 18(9), 509517. Bey, M.J., Kline, S.K., Tash man, S., Zauel, R., 2008. Accu racy of biplane xray imaging combined with modelbased tracking for measur ing invivo patellofemoral joint motion. Journal of Orthopaedic Sur gery and Research 3, 38. Canny, J., 1986. A Computational Approach to EdgeDetection. IEEE Transactions on Pattern Analysis and Machine Intelligence 8(6), 679698. Fregly, B.J., Sawyer, W.G., Harman, M.K., Banks, S.A ., 2005. Computational wear prediction of a total knee repl acement from in vivo kinema tics. Journal of Biomechanics 38(2), 305314. Goffe, W.L., Ferrier, G.D., Rogers, J., 1994. Gl obal Optimization of Statistical Functions with Simulated Annealing. Journal of Econometrics 60(12), 6599. Hirokawa, S., Abrar Hossain, M., Kihara, Y., Ariyoshi, S., 2008. A 3D kinematic estimation of knee prosthesis using Xray pr ojection images: clinical assessment of the improved algorithm for fluoroscopy images Medical & Biolog ical Engineering & Computing 46(12), 12531262. Kappel, F., Kuntsevich, A. V., 2000. An implementation of Shor's ralgorithm. Computational Optimization and Applications 15(2), 193205. 65 PAGE 66 Kaptein, B.L., Valstar, E.R ., Stoel, B.C., Reiber, H.C., N e lissen, R.G., 2007. Clinical validation of modelbased RSA for a total knee prosthesis. Clinical Orthopaedics and Related Research(464), 205209. Kaptein, B.L., Valstar, E.R ., Stoel, B.C., Rozing, P.M., Reiber, J.H.C., 2003. A new modelbased RSA method validated using CAD models and models from reversed engineering. Journal of Biomechanics 36(6), 873882. Kirkpatrick, S., Gelatt, C.D., Vecchi, M.P., 1983. Optimiza tion by Simulated Annealing. Science 220(4598), 671680. Komistek, R.D., Dennis, D.A., Mahfouz, M., 2003. In vivo fluoroscopic analysis of the normal human knee. Clinica l Orthopaedics and Related Research(410), 6981. Komistek, R.D., Dennis, D.A ., Ochoa, J.A., Haas, B.D., Hammill, C., 2002. In vivo comparison of hip separation after metalonmetal or metalonpolyethylene total hip arthroplasty. Journal of Bone and Joint Su rgeryAmerican Volume 84A(10), 18361841. Kon, Y., Nishinaka, N., Gamada, K., Tsutsu i, H., Banks, S.A., 2008. The influence of handheld weight on the scapulohumeral rhythm Journal of Shoulder and Elbow Surgery 17(6), 943946. Lavallee, S., Szeliski, R., 1995. Recoveri ng the position and ori entation of freeform objects from image contours using 3D dist ance maps. IEEE Transactions on Pattern Analysis and Machine Intelligence 17(4), 378390. Li, G., Papannagari, R., Nha, K.W., DeFrate, L.E., Gill, T.J., Ru bash, H.E., 2007. The coupled motion of the femur and patella duri ng in vivo weightbearing knee flexion. Journal of Biomechanical EngineeringTransactions of the Asme 129(6), 937943. Li, G., Wuerz, T.H., DeFrate, L.E., 2004. F easibility of using or thogonal fluoroscopic images to measure in vivo joint kinemati cs. Journal of Biomechanical Engineering 126(2), 314318. Mahfouz, M.R., Hoff, W.A., Komistek, R.D., Dennis, D. A., 2003. A robust method for registration of threedimens ional knee implant models to twodimensional fluoroscopy images. IEEE Transactions on M edical Imaging 22(12), 15611574. Matsuki, K.O., Matsuki, K., Mu, S., Sasho, T., Nakagawa, K., Ochiai, N., Takahashi, K., Banks, S.A., 2010. In vivo 3D kinematics of normal fore arms: Analysis of dynamic forearm rotation. Clinical Biomechani cs doi:10.1016/j.clinbi omech.2010.07.006. McDonald, C.P., Bachison, C.C., Chang, V ., Bartol, S.W., Bey, M.J., 2010. Threedimensional dynamic in vivo motion of the cervical spine: assessment of measurement accuracy and preliminary findings Spine Journal 10(6), 497504. 66 PAGE 67 Morooka, T., Hamai, S., Miura, H., Shimoto, T., Higak i, H ., Fregly, B.J., Iwamoto, Y., Banks, S.A., 2007. Can magnetic resonanc e imagingderived bone models be used for accurate motion measurement with singlepl ane threedimensional shape registration? Journal of Orthopaedic Research 25(7), 867872. Morooka, T., Hamai, S., Miura, H., Shimoto, T., Higaki, H ., Fregly, B.J., Iwamoto, Y., Banks, S.A., 2008. Dynamic activity dependenc e of in vivo normal knee kinematics. Journal of Orthopaedic Research 26(4), 428434. Nishinaka, N., Tsutsui, H., Mihara, K., Suzuki K., Makiuchi, D., Kon, Y., Wright, T.W., Moser, M.W., Gamada, K., Sugimoto, H., B anks, S.A., 2008. Determination of in vivo glenohumeral translation using fluoroscopy and shapematching techniques. Journal of Shoulder and Elbow Su rgery 17(2), 319322. Prins, A.H., Kaptein, B.L., Stoel, B.C., Reiber, J.H. Valstar, E.R., 2010. Detecting femurinsert collisions to improve precision of fluoroscopic knee arthroplasty analysis. Journal of Biomechanics 43(4), 694700. Shor, N.Z., 1985. Minimization methods for nondifferentiable functions, Springer Series in Computational Mathematics, vo l. 3. SpringerVerlag: Berlin. Stiehl, J.B., Dennis, D.A., Komistek, R.D., Crane, H.S., 1999. In vivo determination of condylar liftoff and screwhome in a mobilebearing total knee arthroplasty. Journal of Arthroplasty 14(3), 293299. Tang, T.S., MacIntyre, N.J., G ill, H.S., Fellows, R.A., Hill, N. A., Wilson, D.R., Ellis, R.E., 2004. Accurate assessment of patellar tra cking using fiducial and intensitybased fluoroscopic techniques. Med Image Anal 8(3), 343351. Tanino, H., Ito, H., Harman, M.K., Matsuno, T., Hodge, W.A. Banks, S.A., 2008. An in vivo model for Intraoperative assessment of impingement and dislocation in total hip arthroplasty. Journal of Arthroplasty 23(5), 714720. Tashman, S., Kolowich, P., Collon, D., Anderson, K., Anders t, W., 2007. Dynamic function of the ACLreconstructed k nee during running. Clinical Orthopaedics and Related Research(454), 6673. Tsai, T.Y., Lu, T.W., Chen, C. M., Kuo, M.Y., Hsu, H.C., 2010. A volumetric modelbased 2D to 3D registration method for measuri ng kinematics of natural knees with singleplane fluoroscopy. Medical Physics 37(3), 12731284. Valstar, E.R., de Jong, F.W., Vrooman, H.A., Rozing, P.M. Reiber, J.H., 2001. Modelbased Roentgen stereophotogrammetry of orthopaedic implants. Journal of Biomechanics 34(6), 715722. 67 PAGE 68 68 Vaz, A.I.F., Vicente, L.N., 2007. A particl e swarm pattern search method for bound constrained global optimization. Journal of Global Optimization 39(2), 197219. Wang, S., Passias, P., Li, G., Li, G., Wood, K., 2008. Measurement of vertebral kinematics using noninvasive image matching method Validat ion and application. Spine 33(11), E355E361. Yamaguchi, S., Gamada, K., Sasho, T., Kato H., Sonoda, M., Banks, S.A., 2009. In vivo kinematics of anterior cruciate ligam ent deficient knees during pivot and squat activities. Clinical Biomechanics (Bristol, Avon) 24(1), 7176. Yamaguchi, S., Sasho, T., Kato, H., Ku royanagi, Y., Banks, S.A., 2009. Ankle and Subtalar Kinematics During DorsiflexionPlantarflexion Activities. Foot & Ankle International 30(4), 361366. Yamazaki, T., Watanabe, T., Nakajima, Y., S ugamoto, K., Tomita, T., Yoshikawa, H., Tamura, S., 2004. Improvement of depth pos ition in 2D/3D registration of knee implants using singleplane fluoroscopy. IEEE Transactions on Medical Imaging 23(5), 602612. You, B.M., Siy, P., Anders t, W., Tashman, S., 2001. In vivo measurement of 3D skeletal kinematics from sequences of bi plane radiographs: application to knee kinematics. IEEE Transactions on Medical Imaging 20(6), 514525. Yushkevich, P.A., Piven, J., Hazlett, H.C., Smith, R.G., Ho, S., G ee, J.C., Gerig, G., 2006. Userguided 3D active contour s egmentation of anatom ical structures: Significantly improved efficiency and reliability. Neuroimage 31(3), 11161128. Zhao, D., Banks, S.A., D'Lima D.D., Colwell, C.W., Fregly, B.J., 2007. In vivo medial and lateral tibial loads during dynamic and high flexion activities. Journal of Orthopaedic Research 25(5), 593602. Zuffi, S., Leardini, A., Catani F., Fantozzi, S., Cappello A., 1999. A modelbased method for the reconstruction of total knee replacement kinematics. IEEE Transactions on Medical Imaging 18(10), 981991. PAGE 69 BIOGRAPHICAL SKETCH Shang Mu received his bachelors degree from Beihang University (Beijing University of Aeronautics and Astronautics), Be ijing, China, in 2005. He joined the University of Florida thereafter. Apart from his pursuit in engi neering, his biggest interest lies in physics. 69 