LIBRARY Mlchlgan State University PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. DATE DUE DATE DUE DATE DUE MSU Is An Affirmative Action/Equal Opportunity Institution $deth TE! APPLICKTION 0F HELICAL AXIS THEORY TO THE STUDY OF HMIAN KINENAIICS: A KINENAIIGALEY DEFINED CENTER FOR BIOHECHANIGAL JOINTS BY Raymond Robert Brodeur A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Metallurgy, Mechanics and Material Science 1991 2% §‘: ‘0 \ Ne 1?) Q ABSTRACT THE APPLICATION OF NELICAL AXIS THEORY TO THE STUDY OF HDIAN KINEMATICS: A KINEMATICALEY DEFINED CENTER FOR BIONECNANICAL JOINTS By Raymond Robert Brodeur The forces acting within a biomechanical joint can be succinctly described by a joint reaction force and a muscle moment. The least invasive means of calculating these quantities is to use inverse dynamics. The muscle moment is used in clinical evaluations, therefore it is imperative that it be determined as accurately as possible; however, the muscle moment is dependent on the point used for the summation of the moments. Anatomically based definitions for a joint center have been suggested, but these are-not applicable for subjects with anatomical abnormalities. An alternative is to define the joint center kinematically. This can be done for two dimensional motion using the instantaneous center of rotation. The purpose of this dissertation is to determine a kinematically defined center for biomechanical joints for three-dimensional motion. A solution can be found using helical axis theory. The instantaneous helical axis (IHA) translates and rotates about a second order helical axis. The first and second order helical axes intersect at the center point of the IHA, the point of minimum acceleration on the IRA. The path of center points is the path of minimum velocity for the entire motion. An algorithm was developed for calculating the IRA center point for continuous three-dimensional motion, including solutions for two singularities. The equation for the velocity of a point relative to the IRA can be written as an overdetermined equation; this provides a means of locating the point of minimum velocity for the entire motion. This equation can be expanded to determine the least squares best-fit polynomial estimate for the center point path. Three methods were investigated for defining the kinematic joint center: the IRA center point algorithm; the point of minimum velocity; and a linear best-fit path of minimum velocity. Kinematic and kinetic data were collected for the right ankle joint of two male runners. An anatomically defined ankle joint center was compared to the three kinematic joint center definitions. The kinematic Center of motion was most often medial to the ankle anatomical center, near the medial maleoli. The linear best-fit path of minimum velocity provides a better means of calculating the center point path since it serves to smooth the IRA center points. The moments about the ankle were determined for each of the four joint center definitions. It was found that moments about the three kinematic joint centers were very consistent, and the inversion/eversion moments were significantly larger than the inversion/eversion moments about the anatomical center. ACKNOWLEDGMENTS The author wishes to express his appreciation to the following people and institutions for making the completion of his Doctor of Philosophy degree possible: To my major professor, Dr. Robert W. Soutas-Little, for his help, patience, and guidance, without which this project would not have been possible. To Palmer College of Chiropractic, for their financial support of my doctoral studies. To the Foundation for Chiropractic Education and Research, for their financial support of my doctoral studies. To Karthy Nair-Brodeur, my wife and closest friend, for her support and encouragement when I needed it most, and for her help in preparing this manuscript. TABLE OF CONTENTS LIST OF TABLES ................................................. LIST OF FIGURES ................................................ Chapter I. II. III. IV. INTRODUCTION ............................................ LITERATURE REVIEW ....................................... 2.1 2.2 2.3 2.4 Historical Background .............................. Modern Developments ................................ Biomechanical Applications of Helical Axis Theory .. A Biomechanical Joint Center ....................... MATHEMATICAL TOOLS ...................................... wwwww \IO‘U‘NH W on Definitions ........................................ Properties of Two Lines in Space ................... Determine the Angular Velocity of a Rigid Body ..... Locating the Instantaneous Helical Axis ............ Instantaneous Helical Axis Between Two Rigid Bodies and the Three Axes Theorem ......................... The Ruled Surfaces Generated by the Fixed and Moving Axodes ...................................... CENTER POINT OF THE INSTANTANEOUS HELICAL AXIS AND A KINEMATICALLY DEFINED COORDINATE SYSTEM ................. 4.1 4.2 4.3 Determining the Center Point of an Instantaneous Helical Axis ....................................... The Rate of Translation and Rotation of the Second Order Helical Axis ......; .......................... A Kinematically Defined Coordinate System .......... A KINEMATICALLY DEFINED JOINT CENTER .................... 5.1 5.2 A Discrete Definition For The Center Point of An Instantaneous Helical Axis ......................... 5.1a The Parallel Helical Axes Singularity ........ 5.1b The Zero Angular Velocity Singularity ........ 5.1c An Algorithm For Determining the Path of the Helical Axis Center Point Given a Series of Discrete Instantaneous Helical Axes .......... A Point of Minimum Velocity and Path of Minimum Velocity ........................................... 3O 30 36 37 40 40 43 45 48 VI. LICATION TO THE ANKLE JOINT .......................... Data Collection .................................... Angular Velocity Results ........................... Helical Axis Location .............................. S Helical Axes Center Point Location ................. 6 Ankle Joint Muscle Moments and the Work Done at the Ankle Joint ........................................ APP 6.1 6.2 Ankle Joint Biomechanics ........................... 6.3 6.4 6. 6. VII. CONCLUSIONS ............................................. APPENDIX I ERROR ANALYSIS ..................................... APPENDIX II NUMERICAL SIMULATION .............................. BIBLIOGRAPHY vi 91 97 104 113 TABLE 6.1 TABLE 6.2 TABLE 6.3 TABLE A2.1 LIST OF TABLES Length and standard deviation of vectors between targets of the foot (in cm.) in the shank coordinate system, during stance phase ........................ Helical axes centroid (HAC) for each subject and trial (the point in space that is closest to all the helical axes) ........................... ....... Minimum velocity center (MVC) for both subjects .. Root mean square error between the calculated center points and the true center point ............ 60 76 76 111 Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure U! U! 0‘ 0‘ 0‘0\ 0‘ 0‘ LIST OF FIGURES Two non-parallel and non-intersecting lines in space ...................... . .................. Foot, Shank, and Ankle IHA ....................... The ruled surfaces generated by the fixed and moving axodes .................................... Two consecutive instantaneous helical axes ....... Vectors relating two consecutive IHA's to points on a rigid body .................................. Two consecutive instantaneous helical axes separated by a finite amount of time ............. Two consecutive parallel instantaneous helical axes followed by a third, non-parallel instantaneous helical axis ....................... Camera set up, calibration space and force plate location ......................................... Target placement on the shank and foot for each subject ..................................... Bones and joints that define the ankle ........... Angular velocity of the ankle for subject A Figure 6.4a X axis (inversion(+)/eversion(-)) Figure 6.4b Y axis (plantar(+)/dorsi(-)) Figure 6.4c Z axis (medial(+)/lateral(-) rot.) Angular velocity of the ankle for subject 8 Figure 6.5a X axis (inversion(+)/eversion(-)) Figure 6.5b Y axis (p1antar(+)/dorsi(-)) ....... : Figure 6.5c Z axis (medial(+)/lateral(-) rot.) The shortest vector from the MAC to each helical axis, comparing the IHA to the IHA/IAA for subject A, trial 3 Figure 6.6a X direction .......................... Figure 6.6b Y direction .......................... Figure 6.6c Z direction .......................... X-Z (top) and Y-Z (bottom) projection of the IHA's (or IAA where appropriate) for subject A trial 1 .......................................... viii 20 25 28 31 33 41 42 54 55 58 6O 61 61 62 63 65 66 69 Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure Figure X-Z (top) IHA's (or trial 2 X-Z (top) IHA's (or trial 3 X-Z (top) IHA‘s (or trial 1 X-Z (top) IHA's (or trial 2 X-Z (top) IHA's (or trial 3 Path of Path of Path of Path of Path of Path of Helical Helical Helical Helical Helical Helical Ankle all 4 Ankle all 4 Total force Vectors for defining the center point position ... Model for generating the simulated data .......... minimum minimum minimum minimum axis axis axis axis axis Y-Z (bottom) projection of the where appropriate) for subject Y-Z (bottom) projection of the where appropriate) for subject Y-Z (bottom) projection of the where appropriate) for subject Y-Z (bottom) projection of the where appropriate) for subject Y-Z (bottom) projection of the where appropriate) for subject axis center center center center center center muscle moment ankle center definitions velocity for minimum velocity for velocity for velocity for velocity for point for point for point for point for point for point for for subject A, trial 1 using ix subject A trial subject A trial subject A trial subject B trial subject B trial minimum velocity for subject B trial subject A trial subject A trial subject A trial subject B trial subject B trial subject B trial work, sagittal plane work and colinear work for the ankle (Subject A, trial 1) .......................................... .......................................... .......................................... .......................................... .......................................... H 2 3 ................... muscle moment for subject B, trial 1 using ankle center definitions 71 72 73 74 77 77 78 78 79 79 80 80 81 81 82 82 86 87 88 103 105 Figure A2.2 Figure A2.3 Figure A2.4 Ruled surfaces generated by the fixed and moving axodes ........................................... 105 Angle between "true and uhoisy ................... 109 Angle between ”true and ”noisy’ substituting the IAA for the IHA for w < 2.0 rad/sec .............. 109 In the study of biokinematics and biokinetics, the primary objectives are to describe the motions of each body segment and to determine the forces involved in generating these motions. These objectives can be achieved by modeling the human body as a series of rigid segments, so that segments such as the foot, shank and thigh are each considered a separate rigid body. The kinematics and kinetics of each body segment are then determined separately. In order to describe the kinematics of one body segment relative to another, or relative to a fixed reference frame, it is necessary that the positions of at least three non-colinear points are known on all of the body segments being studied. The velocity and acceleration of each point can be determined by numerical differentiation. The angular velocity and angular acceleration can then be calculated for each body segment (Verstraete and Soutas-Little, 1990). The forces between two adjacent body segments are transmitted via the articular joint surfaces, and the muscles and ligaments connecting the two segments. These forces can be succinctly described by a joint reaction force and a joint muscle moment. The joint reaction force is the sum of the forces acting at a joint, including external forces, inertial forces and the muscle contraction forces between the two segments. The effective muscle moment acting at a joint is the net moment generated by the muscles and ligaments connecting the two body segments. The joint muscle moment is an extremely important factor in biomechanical analysis since it is often used as a comparison between "normal" and "abnormal" biomechanics and in pre- and post-operative comparisons (Prodromos, et a1, 1985, Wang, et a1, 1990). The least invasive means of calculating the joint reaction force and the joint muscle moment is to use inverse dynamics. External forces and moments acting on each body segment can be measured via a force platform, load cell, or pressure plate. The inertial properties of the body segments can be measured using methods such as those outlined by Clauser, et al (1969); or the body segments can be modeled mathematically and the inertial properties can be calculated directly from the idealized geometry of a mathematical model (Verstraete, 1988). The joint reaction force is found by summing the external forces and inertial forces acting on each body segment. The net muscle moment acting at a joint can then be calculated by summing the moments about a point that defines the joint center. Unfortunately, it is not clear as to what defines a joint center. In three dimensional kinetic analysis, the point in space used for the summation of the moments has a direct effect on the magnitude and direction of the calculated moments. This is due to the r X F term in the moment equation. Thus the joint muscle moment is highly dependent on the definition of the joint center. Many researchers have used anatomical definitions for the joint center (Procter and Paul, 1982, Prodromos, et a1, 1986, Verstraete, 1988, Wang, et al, 1990). The center of a joint is defined relative to anatomical landmarks and is based on assumptions of the motions which take place at the joint. One major disadvantage with an anatomically defined joint center is the question of how to modify the definition when dealing with subjects that have anatomical abnormalities. If the motion under consideration is planar, then the instantaneous center of rotation (ICR) is the natural definition for the joint center between two body segments. Unfortunately most biological joints cannot be modeled as having planar motion due to significant nonplanar axes of rotation. The purpose of this research is to define a biomechanical joint center for three dimensional motion. By considering the acceleration of a rigid body it will be shown that a unique point exists on the instantaneous helical axis, and that this point is the center of motion for a rigid body moving in three dimensions. The motion of a rigid body can be uniquely defined by the instantaneous helical axis (IHA). In general, the helical axis itself translates and rotates about a second order helical axis. The second order helical axis intersects the first order helical axis at a unique point, called the center point of the helical axis. It will be shown that the center point of an instantaneous helical axis is the point on the IHA where the acceleration is minimal. It will also be shown that the path of center points is the path of minimum velocity for the entire motion. Thus the center point path is the path of points where motion is minimal; hence the center point of each IHA is a naturally defined center of motion for three dimensional kinematics. In order to apply the helical axis center point to define the joint center for biomechanical joints, some adaptation is required for two singularities where the center point cannot be defined. The first singularity occurs when the angular acceleration vector is parallel to the angular velocity vector: the center point is not unique for this condition. The second singularity occurs when the angular velocity vector is zero: if the angular velocity vector is zero, then the instantaneous helical axis may not exist. The above two problems are addressed in Chapter 5. Helical axis theory is essential to this research, therefore, in Chapter 2, a portion of the literature review is dedicated to a brief history of helical axis theory, including recent developments. The biomechanical application of helical axis theory is also reviewed. Finally, the literature regarding the definition of the center of a biomechanical joint is reviewed. In the third chapter, the fundamental definitions of helical axis theory are presented. The equations of a line, the extension of those equations to a helical axis, and basic mathematical definitions are made. In the last part of Chapter 3, there is a discussion on the relative motion of two rigid bodies in space. The helical axis of one body moving relative to another body is defined and the three axis theorem is reviewed. I In Chapter-4, the center point of the instantaneous helical axis is defined. It is shown that the center point of an instantaneous helical axis is at the point of minimal acceleration on the helical axis. The center point is the center of motion for three dimensional kinematics and therefore is the obvious choice for defining the center of a three- dimensional joint. A natural result of the definition of the center point is a kinematically defined coordinate system. Some discussion is given to the physical interpretation of this naturally defined coordinate system. In Chapter 5, two singular cases are investigated wherein the center point is not defined. The solutions to these problems are found in higher order derivatives of the helical axis. The application of these theoretical solutions to develop a practical algorithm for defining the center point is then given. A least squares approach to the joint center problem is also given in this chapter. The least squares approach defines a point of minimum velocity for the entire motion. It is shown that the minimal velocity point is a weighted average of the helical axis center points. This point of minimum velocity can be used to define the joint center, if the motion at a joint is known to be about a fixed point such as the motion at the hip joint. It is also shown that the path of center points is the path of minimum velocity for the entire motion. The center point path can be modeled using a least squares, best-fit path of minimal velocity. The best-fit path of minimum velocity acts as a filter of the center points. Thus, in this chapter, three center point definitions were developed: the algorithm for calculating the center point for each IHA; the least squares point of minimum velocity; and the least squares best-fit linear path of minimum velocity. In Chapter 6, the above three center point definitions are applied to the ankle joint. The angular velocity vector results and the positions of the instantaneous helical axis are given. The center point of each helical axis is found using the algorithm developed in Chapter 5. The helical axis center points are compared to the point of miminum velocity and to the least squares best-fit linear path of minimum velocity. The moments about an anatomically defined ankle joint center are compared to the moments about the three kinematically defined center points. The total work done at the ankle is compared to the work done by the sagittal plane moment (the dorsi- plantar-flexion moment). The work done by the colinear force was also determined and compared to the total work done at the ankle. In Chapter 7, the theoretical and experimental results are summarized and the conclusions of this dissertation are given. In Appendix I an error analysis is done for the angular velocity vector equation; the location of the instantaneous helical axis equation; and for the equation used to determine the location of the center point of each helical axis. Appendix II shows the application of the theory developed in this dissertation to numerically simulated data. Random noise is introduced to simulated target points. The velocity, angular velocity, angular acceleration and helical axis positions are determined given only the noise induced data. The results from the noise induced simulated target data are then compared to the known results. CHAPTER II LITERATURE REVIEH 2.1 W Chalses (1830) is generally regarded as being the first to show that a finite or infinitesimal motion of a rigid body can be described by a translation and a rotation about a line in space. However, it has been reported by Rittershaus (1878) and Schoenflies (1893) that the screw axis was first described by Giulio Mozzi in 1763 or 1765 (see Hunt, 1967). Poinsot (1806) was the first to show that the forces acting on a rigid body can be reduced to a single force and a couple in a plane perpendicular to the force. The mid and late nineteenth century saw a rapid growth in the area of rigid body kinematics with the works of Poinsot (1851), Hamilton (1830,1845,1848), Mobius (1837,1838), and Plucker (1865,1866). In order to apply the helical axis theory to the study of mechanisms several means of expressing screw coordinates were developed. In 1865, Plucker defined his six coordinates of a line. In this publication (titled "On a New Geometry of Space"), he defined a line based geometry and applied it to optics, kinematics and kinetics. Hamilton developed quaternions, a vector and scalar based system that used traditional algebraic multiplication and complex number definitions to define what we today call the vector dot product and cross product. Hamilton's quaternions looked much like our modern vectors; the multiplication of two quaternions results in the equivalent of a dot product and a cross product in terms of modern vector terminology. In 1873, Clifford developed a mathematics based on a new type of complex number, where the multiplier w was used, and defined by w-w - 0. He termed this complex number system dual numbers, consisting of a real number and the dual quantity, w. He combined dual numbers with Hamilton's quaternion notation to develop biquaternions. In 1900, Ball published his monumental work on the theory of screws. Ball used a screw coordinate system based on six parameters, five of which must remain independant. Ball studied the infinitesimal displacements of a rigid body having between zero and six degrees of freedom. His major emphasis was on the restraining forces necessary to maintain the given degree of freedom and the motion that would be imparted to a body given the restraints, the initial position of the body and the direction of an impulse acting on the body. He is known for his geometrical intuition in his approach to kinematics. Study (1903) used Clifford's dual numbers to represent a screw as a dual vector. A dual vector was defined as two vectors; the first being a real vector that defined the direction of a line in space, the other being the moment of the line about the origin multiplied by the dual multiplier, w. Dual vectors are similar to Plucker coordinates, in that they both use the directional cosines of a line and the moment of a line about the origin as a means of describing a line in space. 2.2 Wilma Little work was done in the area of screw theory between the early twentieth century and the post World War Two era. Dimentberg (1948) adopted Study's notation and used it in the analysis of kinematic chains and in the finite displacement of general three dimensional mechanisms. In 1965, Dimentberg published "The Screw Calculus and Its Application in Mechanics", a work that summarized his substantial contribution to mechanism design. The publication also helped raise Study’s notation from obscurity. A good portion of the helical axis theory literature over the past thirty years has been the rediscovery and further advancement of work that was done in the late ninteenth century, but re-written in terms of modern vector and matrix notation. A. T. Yang (1963) used dual vectors and re-introduced the concept of dual quaternions to the analysis of spatial mechanisms. He combined the work of Clifford, Hamilton and Study to define vectors and line geometery in terms of dual quaternions. The result of this is a very concise and elegant notation for expressing the motions of a mechanism. Yang (1971) expanded dual notation to express the dynamics of a rigid body. He expressed the acceleration of a rigid body and the rate of change of momentum of the body as one dual vector equation. This lead Beggs (see the discussion of Yang, 1971) to comment that Yang's dual dynamic equation "... compares with a DNA molecule in information density!". Woo and Freudenstein (1970) and Yuan and Freudenstein (1971) precede their work by extensively reviewing the notation and findings of Plucker. Clifford, Study and Ball. The above two papers extended the application of screw coordinates to the kinematics and statics of mechanisms. The application of screw coordinates to rigid body motion and mechanisms provided the development of more efficient methods for studying the displacement, velocity, acceleration, and forces acting on spatial mechanisms. Veldkamp (1963, 1967a, 1967b) defined and described the concepts of instantaneous invariants and a canonical system of reference for rigid body motion. The instantaneous invariants are scalar parameters that 10 describe the motion of a rigid body and are independent of the coordinate system used to define the motion of the body. For first order motion (velocity), there are two instantaneous invariants, the first being the magnitude of the angular velocity, the second being the translational velocity parallel to the angular velocity vector. The canonical system of reference is a kinematically defined coordinate system that has the instantaneous helical axis as one line of the . coordinate system and the second order helical axis as a second line in the coordinate system. The third coordinate axis is the cross product of the first two. The canonical reference system can be an extremely useful tool in the analysis of spatial mechanisms. There have been a number of papers concerned with the finite helical axis (FHA) and the finite center of rotation for biomechanical joints (Hicks, 1953, Close, 1956, Van Langelaan, 1983, Lundberg, et a1, 1989, Blankevoort, et a1, 1990). Most of the finite helical axis research has been collected using stereoradiogrammetry and human cadavers and/or excised joints. The biomechanical researcher faces a number of dilemmas when investigating joint motion. Extremely precise data can be collected on cadavers (or excised joints), using radiopaque targets implanted into the bones that define the joint under investigation. However, by gaining a higher degree of accuracy the researcher sacrifices functional normality. 0n the other hand, when using living subjects and skin mounted targets, the function of the joint may be "normal" but there is usually a loss in accuracy due to tissue motion. v"?! ( .l‘a-‘ “‘v“. ‘Jfl‘r‘ .1. .ws—y «8' 11 The application of helical axis theory for describing biomechanical motion date as far back as Weber and Weber (1836) who made the first recorded attempt to measure the center of rotation of the knee joint (using a planar model for describing knee motion). In this century several researchers have applied the finite helical axis or its planar equivalent, the finite center of rotation, to biological joints (Manter, 1941; Hicks, 1953; Isman and Inman, 1969; Walker, et a1, 1972; Smidt, 1973; Blacharski, et a1, 1975; Procter and Paul, 1982; Lundberg, et al, 1989; and Blankevoort, et a1, 1990). The methods for determining the FHA vary widely between the papers listed above. Hicks (1953) and Isman and Inman (1969) used simple visual techniques for defining the axis of rotation. They defined the axis of rotation to be through points of minimal motion. For example, for talo- tibial motion, a point of minimum movement was visually determined on the medial and lateral sides; then the axis of rotation was simply the line between the points of minimal motion. More recently, stereoradiographic and biplanar radiographic methods have been developed for finite helical axis measurements (Selvik, 1974). Van Langelaan (1983) did an extensive radiographic study of the tarsal joints using ten shank-foot preparations. He described the FHA for the talotibial, talocalcaneal, cuboid calcaneal, naviculocalcaneal, and the talonavicular joints. Lundberg, et a1 (1989) studied the talotibial joint in vivo using eight volunteers. Tantalum balls (0.8mm diameter) were implanted into the tibia and talus of each subject. Weight bearing flexion-extension; pronation-supination; and medial-lateral rotation were performed at 10° increments. They found that all of the axes tended to cross near a point in the center of the talus. Although many studies have been done on the finite helical axes, few of the studies agree on the location and variation of the axes 12 'locations. Blankevoort, et al (1990) attribute this to several factors: 1) individual subject variation; ii) measurement errors; iii) variation due to different loading and variations in the kinematic paths forced on the body segments by the researchers. Blankevoort, et a1 expressed the opinion that the third factor is probably the factor of greatest importance. They, therefore, described in great detail the loading and kinematic paths of the femur and tibia in their work. However, the problem in general with biomechanical research is control over the repeatability of the experiment. Very little work has been done in the area of instantaneous helical axes for biomechanical joints. Fioretti, et a1 (1990) used simulated data to determine the accuracy of calculating the IHA. They simulated a cylinder rolling on a plane, introducing noise of 0.5mm and 4.0mm standard deviation to the known target locations. For the 0.5mm error, the location of the IHA could be determined to within 1.27mm using a readily available spline smoothing technique (Woltring, 1986). For the 4.0mm standard deviation noise, the IHA could be determined to within 5.09mm using the same smoothing technique. Fioretti et al also collected in 1129 data on the metacarpophalangeal joint and described the IHA path as it intersected the sagital plane of the metacarpal. Verstraete and Soutas-Little (1990) described a method for finding the angular velocity and angular acceleration from discrete position data. They used numerical differentiation to determine the velocity and acceleration of the position data, then used a least squares method to find the angular velocity and angular acceleration. Sommer and Buczek (1990a, 1990b) also described a least squares method for determining the angular velocity and angular acceleration vector from discrete target position data. They determined the accuracy of the angular velocity and acceleration on a fixed axis mechanism. 13 They found the average measured angular velocity and angular acceleration to be very close to the true angular velocity and acceleration. Karlsson et a1 (1990) compared bone mounted targets to skin mounted targets on the location of tibial-femoral instantaneous helical axodes. As would be expected, they found that the IHA's found from bone mounted target data had less variation than skin mounted targets. 2.4 Biome a1 Jo Cent In the study of biomechanics, the inverse dynamics method is used to determine the forces and moments acting at a joint. Body segments such as the foot, shank and thigh are modeled as rigid bodies. External forces acting on these bodies can be measured via a force platform, load cell, or pressure plate. The inertial properties of the body segments can be measured (Clauser, et al, 1969) or the body segments can be modeled mathematically and the inertial properties can be calculated directly from the idealized geometry of a mathematical model (Verstraete, 1988). Since the external forces and moments are known, the forces and moments acting at a joint can be calculated using the inverse dynamics method. The moment acting at a biomechanical joint is called the joint muscle moment. It is the effective moment acting at the joint due to the forces generated by the muscles and ligaments connecting the two body segments. The magnitude of the joint muscle moment is an extremely important factor in the analysis of a biomechanical joint. The muscle moment is often used as a comparison between "normal" and "abnormal" biomechanics and in pre- and post-operative comparisons (Prodromos, et a1, 1985, Wang, et a1, 1990). 14 Unfortunately, it is not clear as to what defines a joint center. In three dimensional kinetic analysis, the point in space used for the summation of the moments acting on a rigid body has a direct effect on the magnitude of the calculated moments. This is due to the r X F term in the moment equation. Thus the joint muscle moment is highly dependent on the definition of the joint center. If the motion under consideration is planar, then the instantaneous center of rotation (ICR) is the natural definition for the joint center between two body segments. Yamaguchi and Zajac (1989) used the instantaneous center of rotation to more accurately determine the effective moment arm of the quadriceps. Unfortunately most biological joints, including the knee joint (Karlsson, et a1, 1990, and Blankevoort, et a1, 1990), cannot be modeled as planar due to significant nonplanar axes of rotation. Many researchers have used anatomical definitions for the joint center (Procter and Paul, 1982, Verstraete, 1988, Prodromos, et alh 1986, Wang, et a1, 1990). The center of a joint is defined relative to anatomical landmarks and is based on assumptions of the motions which take place at the joint. One major disadvantage with an anatomically defined joint center is the question of how to modify the definition when dealing with subjects with anatomical abnormalities. The need for a three dimensional biomechanical joint center has been recognized by several researchers. Chao and An (1982) proposed that a joint center be defined by the midpoint of the common perpendicular between two consecutive helical axes. They proposed that the path of this midpoint could be used as a description of the motion taking place at a joint. Woltring (1990a) has proposed that the intersection of the first and second order helical axis be used to 15 define a joint center. It will be shown in this thesis that, in the limit, the above two definitions converge. Much of the debate regarding the definition of a joint center has taken place on Biomch-L, a biomechanics electronic bulletin board. Woltring (1990b) suggested a point at which the absolute velocity and acceleration are minimal. For any given instant in time, this point occurs at the center point of the IHA. Chao (1990) stated that the three dimensional joint center as defined by Woltring (1990b) should be utilized for the calculation of joint muscle moments. Chao indicated that the joint articulating surface contact point does not necessarily coincide with the center of motion of the joint. Joint muscle models require the moment arm between the line of action of a muscle and the center point of the joint in order to accurately model the joint kinematics and kinetics. Spoor, et a1 (1990) suggested that the moment arm be determined using the relationship between the work done by a muscle and the work done by the moment of the muscle. For planar motion, the work done by a muscle generating a force F is simply F ds where ds is the infinitesimal change in the tendon length. The work done by the muscle moment on the joint is M do where M is the moment generated by the muscle and d0 is the change in angle of the joint. Then: F ds - M d0 the effective moment arm "a" is: M ds a-_-_ F d9 This relationship allows the researcher to calculate the moment arm by relating the change in tendon position (ds) to the change in the joint angle (d0); the actual moment arm does not need to be measured. This method assumes that the motion is planar; an assumption that is not true 16 for most biological joints. The other obvious disadvantage to this method is that tendon length changes are difficult, if not impossible to measure in vivo. The purpose of this research is to define a biomechanical joint center based on the kinematics of the body segments that define a joint. Helical axis theory is essential to this research, therefore, a review is provided for basic definitions and mathematical tools that are used in this thesis. CHAPTER III MATHEMATICAL TOOLS 3-1 Wm Six independent parameters are needed to describe the motion of a rigid body in space. A line in space can be defined by four independent parameters. Equation 3.1 defines the projection of a line onto the x-z and y-z planes, respectively: 3.1a x - az + b 3.1b y - C: + d The parameters a,b,c,d are four independent parameters that can be used to define a line in three dimensional space. The rotation and translation along a line account for two more degrees of freedom; therefore the motion of a rigid body can always be expressed as a translation and a rotation about a line. The line about which a rigid body is translating and rotating is call the helical axis. The ratio of the translation to the rotation is called the piggh of the helical axis. Thus the motion of a rigid body can be described by the four parameters that define a line in space and two other parameters, the pitch and the rotation. Alternatively, the rotation and translation can be used instead of the rotation and the pitch, due to the relationship between the rotation, translation and pitch. The motion of a body can be described in terms of finite relative positional changes or instantaneous motion. The finite hellcal axis (FHA) describes the line in space about which a body translates and rotates in order to move from one position to some other position in space. However, it must be emphasised that the FHA describes the most 17 18 Snocinct means by which a body may move from one position to another position; it does not necessarily describe the aaag; means by which a body has moved. The motion along a helical axis is described in terms of the £315; and the pitch, where the twist is the magnitude of the rotation of the body. As the motion of the body becomes infinitesimally small, the helical axis is referred to as the laa;an§aaagaa_hallgal_a§la (IHA). The IHA is not only the most succinct means of describing the motion of the body, it is also the exact means by which the body moves from one position to another new position. Since the infinitesimal translations and rotations occur over an infinitesimal time, the ratio of the differential quantities to the infinitesimal time do not go to zero, but, in the limit, become the translational velocity and the angular velocity, respectively. The pitch of the instantaneous helical axis is the ratio of the translational velocity to the angular velocity. The 'gglaglng_mg§lgn of the body is the magnitude of the angular velocity. In a similar fashion, the forces acting on a body can be reduced to a force and a colinear moment. That is, any number of forces and moments acting on a rigid body can be reduced to a single force and moment acting along a line in space. This combination of a force and a colinear moment is called a grangh. The magnitude of the wrench is equal to the magnitude of the resultant force. The pitch of the wrench is equal to the ratio of the colinear moment to the force. The free translation of a body or a pure moment acting on a body cannot be reduced to acting on a single line in space. However, in terms of helical axis theory, a pure translation is an infinitesimally small rotation about a line at infinity; and a pure moment is an infintesimally small force acting on the body on a line at infinity. Thus. a helical axis is the sum of two rotations. One is the rotation 19 0f the body about the line of the helical axis; the other is a rotation about a line at infinity. Similarly for the wrench: it is the sum of two forces acting on the body. One is the force acting along the line of the wrench; the other is a force acting on the body at a line at infinity. Helical axis theory provides a succinct and concise means of describing the kinematics and kinetics of a rigid body and for describing the relative motions across an anatomical joint. In addition, the axis of rotation is essential for any type of muscle modeling (Spoor, et a1, 1990), since the moment created by the muscle force is dependent on the moment arm from the axis of rotation to the muscle tendon. For the inverse dynamics problem, the axis of rotation defines a line, but a unique point on that line is needed for the summation of moments, in order to accurately ascertain the net muscle moment acting across a joint. 3.2 Iggnartles of 129 Liaas lg §2ace This section reviews the basic geometry of two lines in space; defines a simple means by which it can be determined if the lines intersect, and defines a means for finding the point on one line that is closest to the second line. Figure 3.1 shows two non-parallel, non~intersecting lines in space; lines "a" and "b". A unit vector is known for each line, defining the direction of each of the lines (ea, 2b). A point is known on each line, where R3 and Rh are vectors from the origin to the known point on lines a and b, respectively. In order to define the point of intersection or the point on one line closest to the other line, it is convenient to use the two lines to 9| ”I Figure 3.1 Two non-parallel and non-intersecting lines in space. define a coordinate system as follows: 3.2 ed- egxfl; leaXebl 3.3 cc - ed X ca 3.4 ab/a-ab-na 21 The unit vectors (Qa,;c,;d) define an orthogonal coordinate system. Place the coordinate system at point q, the point on line "a" closest to line "b". The vector q is the vector from the origin to point q, as shown in Figure 3.1. Lines "a" and "b" intersect if: A 3.5 nb/a - ed - 0 If the two lines do not intersect then: 3.6 Rb/a - ed - d where d is the shortest distance between the two lines. The vector q can be found from the simple geometry defined by the two lines and the coordinate system defined above and shown in Figure 3.2: A 5:.ea 3'7 q - 3a + [Rb/a. ea - tan(0) lea If the two lines are parallel, then q cannot be defined. The dot product of the two lines can be defined using the above coordinate system. Let the wrench acting on a body lie on line "b"; then the force and moment acting on the body are: 3.8 F - F ;b M - M ;b Let the velocity of the body be expressed by an instantaneous helical axis on line "a". The veloctiy of the body is then: 3.9 V - V ea a - w ea The distance the body travels along the helical axis line "a" during the application of the force and moment is: 3.10 ds - V dt dd - u dt The work done by the wrench on the body is the dot product of the force and moment with the distance the body travels. In order to calculate the dot product, the force and moment must be transfered to a point on the line "a". The equivalent force and moment acting on point q is: 3.11a Fa - F cos(0) Ma - M cos(0) - F d sin(0) 22 3.11b FC - F sin(0) Mc - M sin(0) + F d cos(9) where Fa and Ma are the force and moment acting in the eadirection and Fc and Mc are the force and moment acting in the ;c direction. Note that there is no force or moment acting along the ed direction at point q. The work on the body is then: 3.12 {F V cos(0) - F d w sin(0) + M w cos(0)}dt If the work done on the body is zero and F,M,V,w are non-zero, then the wrench is said to be reciprocal to the twist. This occurs if the two lines of action are orthogonal and intersecting, but can also occur for an infinite number of combinations of F,M,V,w, and 0. The interesting note about this result is that screws do not necessarily have to be perpendicular to be orthogonal. The implications of this finding have been applied to the study of constraint problems in mechanism design and robotics (Ball, 1900, Ohwovoriole, 1981; Lipkin and Duffy, 1982). 3.5 Beta ' e vs id The purpose of this section is to review one method of calculating the angular velocity vector from discrete position data. Given a solid body in general three-dimensional motion and at least three non-colinear points rigidly attached to the body, the problem is to find the angular velocity for any instant in time. The discrete position data can be numerically differentiated to determine the velocity of each point. At time t, the position and velocity of the three points B,C,D are: 3.12a rB, rc, rD 3.12b vb, vc, VD 23 The position and velocity of any two of the points can be expressed relative to any other point. Thus, relative to point D, the position and velocity of points B and C are: 3.13a r B/D ' ‘B '_rD 3.13b rC/D - rc - rD 3.13c vD/D - VB - vb 3.13d vC/D - vb - vb The angular velocity of the body can be found for time t: v X v 3.14 w - B/D C/D v'B/D ' rC/D Unfortunately equation 3.14 does not provide a reliable means of determining the angular velocity, especially if the vector rC/Dis perpendicular with the velocity vector vB/D vB/D - rC/D This occurs if the angular velocity vector is parallel to the plane (making - 0). formed by the points B,C,D. Verstraete and Soutas-Little (1990) developed a means for determining the angular velocity using a least-squares approach. Sommer and Buczek (1990a, 1990b) independently developed the sameethod for defining the angular velocity and acceleration. The velocity of one point relative to another must always be perpendicular to the vector between the two points (for a rigid body). Thus: /D ' 'b/D The cross product can be written in terms of a skew symmetric matrix: 3.15 a X to 0 -r r w z x 3.16 uXr- - r 0 -t:x w --[£]u -r r 0 24 Where the "~" below the r indicates that [5] is a skew symmetric matrix. Thus equation 3.15 can be written as: 3.17 ' [EC/D] " ' 'b/D Equation 3.17 can be written in an overdetermined fashion by repeating the same process for all of the points on the body: V 53/0 B/D 3.18 - EC/D a - vC/D 50/3 vC/B The "~" indicates that each of the vectors r in equation 3.18 are written as (3x3) skew symmetric matrices as in equation 3.16. By solving 3.18 using the least square method, the angular velocity u can be obtained. The accuracy of the angular velocity increases with the number of targets attached to the rigid body. Once a is known, the IHA can be located relative to any of the known points. To find the vector from a point D to a point P on the IHA (where UP is the shortest perpendicular vector from D to the IHA): 3.19 B? - ” x 'b 2 |"l Thus the instantaneous screw axis has angular velocity u, translational velocity VD - u/|u|, and is located at: 3.20 6? - 66 + B? where "0" is the origin of the coordinate system, and P is a point on the IHA as defined by equation 3.19. Sections 3.5 and 3.6 were concerned with calculating the instantaneous helical axis for a rigid body relative to a fixed coordinate system. This section is concerned with calculating the IHA between two bodies. Figure 3.2 shows the foot with a foot coordinate Figure 3.2 Foot, shank, and ankle IHA. system attached at point B and the shank with a shank coordinate system attached at point A. The position and velocity of point A on the shank and point B on the foot are known; and the angular velocity of both the foot and shank (of and us), relative to the global XYZ coordinate system, are also known. 26 The position of point B can be defined as: 3.21 RB - RA + nB/A The velocity of B is: 3.22 VB_VA+vB/s+"san/A where Vh/s is the velocity of point B relative to the shank coordinate system. Therefore: 3.23 Vh/S - VB - (VA + as X nB/A) The angular velocity of the ankle joint, 63, is simply the difference between the foot and shank angular velocities: 3.24 U3 - uf - as The instantaneous helical axis of the foot relative to the shank (IHAf/S) can be found relative to point B on the foot using equation 3.19: 3.25 DB - “ x V55, 0 . w J .1 where D is the shortest vector from point B to the IHA of the foot B relative to the shank. The translational velocity of the foot along the IHAf/s is: 3.26 A A vparallel - (VB/A ' e'j)ej The triple axis theorem (Phillips and Hunt, 1964, Ball, 1900) states that given three bodies, U,V,W moving in space with relative instantaneous helical axes IHAU/v, IHAV/W' IHAw/U, the three relative IHA's share a common perpendicular line. The theorem can be easily proven given the above equations. Let the three bodies be the shank, the foot, and the fixed coordinate frame. IHAs is the IHA of the shank relative to the fixed coordinate frame. IHA is the IHA of the foot f relative to the fixed coordinate frame. IHAf/S relative to the shank. If IHAf and IHAs are known, then find the is the IHA of the foot location of IHAf/s. 27 In equation 3.23, the points A and B are any points on the shank and foot respectively. Therefore, without any loss in generality, the point A can be located on IHAS at the point on IHAs that is closest to IHAf. Similarly point B can be located on IHAf at the point on IHAf that is closest to IHAS. Then AB is the common perpendicular between IHAS and IHA Substituting equation 3.23 into equation 3.25, the f' shortest vector from point B to the IHAf/s can be found (DB): 3.27 D _ (uf- as) X [VB - (YA,+ 05X AB) (”f- as) . (”f- Us) _ -0fX VA - qu (05X AB) - uSX VB + 08X (05X AB) uf- "f - Zuf- as + 05- ”5 Point A lies on IHAs therefore VA is parallel to as, and point 8 lies on IHAf therefore VB is parallel to of. An examination of equation 3.27 reveals that the numerator is a vector parallel with vector AB. Thus vector DB lies on the line that passes through line segment AB. Therefore IHAf/s shares the same common perpendicular line as the common perpendicular between IHAf and IHAS. An alternative means of calculating the instantaneous helical axis of the foot relative to the shank is to define the foot targets relative to the shank coordinate system and calculate the velocity of the foot targets relative to the shank coordinate system using standard numerical differentiation techniques. The angular velocity of the foot can be determined from equation 3.18 and the location of the IHAf/s can be found from equation 3.19. Figure 3.3 shows a rigid body moving relative to a fixed coordinate system. At any given instant in time the rigid body is translating and rotating about the instantaneous helical axis. Relative to the fixed frame of reference the IHA is called the fixed axode (Skreiner, 1966). Relative to the moving coordinate system the IHA is called the moving axode. The IHA lies on the line in space that is shared by both the fixed and moving coordinate systems. The motion of a line in space generates a ruled surface. Thus as the motion of the rigid body generates a new IHA position, two ruled surfaces are formed by the changing IHA position. The motion of the fixed axode generates one Fixed Axode T Ruled Surface Figure 3.3 The ruled surfaces generated by the fixed and moving axodes. 29 ruled surface; the motion of the moving axode generates a second ruled surface. During the motion of the rigid body the two ruled surfaces roll relative to each other, always maintaining one line in contact; that line is coincident with the IHA. The shape of the two ruled surfaces are dependent on the motion of the body that generated the IHA. Kirson and Yang (1978) and Yang, et a1 (1975) have proposed that instantaneous invariants for describing ruled surface geometry be used to describe the motions of spatial mechanisms and to aid in the design and synthesis of mechanisms. CHAPTER IV CENTER POINT OF THE INSTANTANEOUS HELICAL AXIS AND A KINEMATTCALLY DEFINED COORDINATE SYSTEM Figure 4.1 shows two instantaneous helical axes, one at time t (IHA) and one at time t + dt (IHA'). In general these two lines are non-intersecting; however, they share a common perpendicular, which is also the shortest distance between IHA and IHA'. Call the unit vector along the common perpendicular 2k. The common perpendicular intersects the IHA at point C. The body is translating and rotating about the IHA, but the IHA is translating a distance ds and rotating by an angle d0 about an axis parallel to 2k located at the point C. This line is a helical axis for the instantaneous helical axis; it is called the second order helical axis. If there is no rotation involved in the change in the IHA (for example, in planar motion all IHA’s are parallel to one another) then there is no unique second order helical axis; it compares to a pure translation of a rigid body, but in this case it is the pure translation of the IHA. More will be discussed on the second order helical axis in section 4.2. Equation 4.1 gives the angular velocity of the body at time t + dt: 4.1 w'-u+cdt The sin and cos of angle d0 can be expressed using a, a, and dt: 30 31 IHA' ds Figure 4.1 Two consecutive instantaneous helical axes. 4.2a sin do - '“’ x ”'l - "” x “d” '0' l"'| MI I“ + adtl 4.2b cos d0 - u . ”I w + u-adt IWI IU'I IUI I" + Gdtl A coordinate system can be defined, as shown in Figure 4.1, by defining the following: A 4.3 e - u w ___ . Ivl 4.4 ek _ a X 0 In X u'I 4.5 eT - ek X ew The coordinate system is located at point C, which is the point of intersection of IHA and the common perpendicular between IHA and IHA'. The point C is called the center point of the IHA. It can be shown that 0 lies in the ow - eT plane. Examine: 32 A “36 aoeK-s Thus: 4.7 °';K"" UXB' _ ¢l°[UX(u+odt)] |o X 0" ID X u'| But a X a is zero, and a X adt is perpendicular to a, so then: 4.8 a - (a X adt) - 0 - 6 Thus 6 - O and the vector u lies entirely in the ew - eT plane. Define: 4.9 ab - (a - ew)ew 4.10 “'1' - (a . eT)e.r Let point D be any point on the body and let the shortest vector from D to the IHA be vector DP. Let D’P' be the shortest (perpendicular) vector from D’ to IHA'. Then: _ u’Xv, (u+adt)X(v+ dt) 4.11 D'P' - D - D 8” |u'|2 |u + adtl2 The vector 55' is the change in the position of D between time t and time t+dt, thus: 4.12 35' - dec From Figure 4.3, the vector PP' is: 4.13 fi'- fiI+fi'-‘D‘P —I -l —I —I A A A The vectors PP“, PPT' PPK are the components of PP in the ew, eT, eK directions, respectively. From Figure 4.2, the magnitude of CE is: _ IPP'TI 4.14 |CE| tan d0 Using equations 4.2, it can be shown that: 1 “2 + w awdt 4.15 - tan d9 w a dt T Figure 4.2 Vectors relating two consecutive IHA's to points on a rigid body. So: __ ”2 + w awdt A 4.16 |CE| - |PP'TI cw w aTdt A It can be shown that the acceleration of the point C in the eK direction is zero. The acceleration of the point P is: 4.17 8P - aD + a X DP + 9 X (a X DP) expressing 4.17 in component form gives: a V T DT 4.18 an - an“ + w a v w DT 4.19 aPT - aDT - w + w vDK a v w DK 4.20 aPK - aDK - - w vDT (a) 34 The acceleration of any other point, F, on the IHA is: 4.21 8F - aP + a X PF + 0 X (a X PF) but FF is parallel to both a and am, thus 11F and aP differ only in the acceleration in the eK direction. Thus: 4.22 aF - a? + «T x F? So: a v T DT 4'23 Fm - Pw - aDu) + w cw a v w DT 4.24 an- 8PT- EDT- to eT+vaKeT 4.25 aFK - a.PK + “T X PF Thus all points on the IHA have the same acceleration in the em and eT directions; they differ only in their relative eK accelerations. From Figure 4.3 the acceleration of the point E can be found: 4.26 aE-aP+aTXfi But PE - PPL, thus the acceleration of the point C is: 4.27 a.c - a.P + “T X (PP; - CE) or: |P_P-i-| _ 4.28 aC - aP + [ - |PPw|] GT eK tan d0 From equations 4.11 - 4.16, the following relationships are found: o X v + a X dt + a X v dt a X v 4.29 PP' - vbdt + D ‘D D - D |u + adtl2 |w|2 The quantities PP; and PPi and PPk can be found: a v dt 4.30 PP' - v dt + T DR w Dw ------§ |u + adtl w v + w a dt + a v dt v 4.31 PPi - vDTdt + ”K ”R “ DK + DK In + adtl2 |u| 35 w v + w a dt + a v dt - a v dt w v 4.32 PP' - v dt+ ”T ”T w DT TDw _ DT K DK 2 2 I0 X adtl '0' From equation 4.28, the acceleration of the point C in the eK direction is: aw VDX w + awdt 4.33 aAK - aDK - ( + w VDT) - aTIPPw] + aT[PPT] w a T Thus, using equation 4.30, 4.31 and using the common denominator: |o| |u + adt|2 - w3 + Zuzawdt (droping higher order differentials), the acceleration of C in the eK direction is: 3 2 w aDK + 2w awaDKdt 4.34 a - GK 2 [0| I0 + adt| 2 4 3 w avaK + 2w avaKdt + w vDT - 2w avaTdt |0| lu + adtl2 3 2 w aTvadt + w aTvDKdt |u| |u + «1:12 4 3 3 2 + w vDT - w aDK + w avaK + w avaTdt - w aDKdt + w avaKdt |u| |u + adt|2 This reduces to: 4.35 aCK The acceleration of the point C in the eK direction is zero. Hence to - order(dt) z 0 find the center point C simply requires finding the point on the IHA at which the acceleration in the 2k direction is zero. It should be noted that if there is no component of a perpendicular to m then the center point is not unique for the IHA. This singular condition will be discussed in section 5.2. The infinitesimal translation and rotation of the IHA are ds and d0. Figure 4.2 shows that ds - |CC'I - IPPkI. From equation 4.32: uvDT + waDTdt + avaTdt - aTvadt w vDT ds - vDKdt + , - In + adt|2 |u|2 The point D is any point on the rigid body. For convenience sake, let the point D be the center point of the ISA, then point D - point C, then, recalling that VCK - vCT - 0, the above equation reduces to: w a dt - a V dt ds _ CT T Cw J/ In + adtl2 The term lu + adtl - “2 + 2wawdt + aidt2 + aidt2 = ”2 + Zwawdt, and: w a - a V a a v 4.36 d: _ CT T Cw a fig: _ T Cw dt w2 + 2waudt w w2 From equation 4.15: w aTdt tand9 - d0 - 7—, u + w 0 dt (d Dividing through by dc gives: w a a 4.37 f:— T .. T dt “2 + w a dc w w Thus the IHA is translating and rotating about the second order helical axis with the translational velocity and rotational velocity of the IHA about the second order helical axis given by equations 4.36 and 4.37 respectively. These values are independent of the coordinate 37 System used to define the helical axes; thus, the quantities given by equations 4.36 and 4.37 are invariant and can be used to describe the motion of the rigid body. Recall in the introduction, it was stated that the angular velocity magnitude and the translational velocity along the IHA are invariant to the coordinate system. Thus, the first order invariants of motion are u and vs. The second order invariants of motion are do/dt and ds/dt. The motion of a line in space describes a ruled surface. For any given instant, the ratio of the translation of the line to the rotation of the line is called the distribution parameter. Thus, the distribution parameter of the IHA is equal to the pitch of the second order helical axis. The distribution parameter is: The distribution parameter for the ruled surface generated by the IHA is the pitch of the second order helical axis. 4.3 LWWM In section 4.1, a coordinate system, based on two consecutive instantaneous helical axes, was used in order to find the center point of a helical axis. That coordinate system is completely defined by the kinematics of a rigid body. Using the angular velocity and angular acceleration of a body, and knowing the velocity and acceleration of at least one point on the rigid body, a coordinate system can be defined: 38 4.39s cw - w/Iul x 0 X a 4.391) 8k - m 4.39c eT - ek X em The coordinate system is located at the center point of the IHA, a point where the magnitude of the acceleration of points on the IHA is at a minimum (the acceleration at the center point is zero in the QR direction). This coordinate system is called the IHA trihedron by Kirson (1975). The significance of the IHA trihedron is that if E - (2”,;T,;k) defines the coordinates system in the fixed frame of reference and E' - (35,;i,;£) defines the coordinate system relative to the moving coordinate system, it can be shown (Veldkamp, 1967) that E - E'. That is, the lines defining the coordinate system are coincident for both the fixed and moving body. Thus, for each instant in time, there exists a coordinate system whereby points in both bodies can be defined. This naturally defined coordinate system is called the canonical coordinate system (Veldkamp, 1967). This coordinate system can be a useful tool in the analysis of spatial mechanisms. There is a physical significance to each of the axes of the IHA trihedron. The em axis lies on the IHA; the rigid body translates and rotates about this line. The 2k axis lies on the second order helical axis. The first order helical axis (the IHA) is translating and rotating about this axis. The ;T axis is pointing in the direction toward which the angular velocity vector is changing. In terms of the linear velocity, em is the line all points in the body are translating parallel to. 39 If the motion of the body is planar it is still possible to define the IHA trihedron, but the origin of the trihedron has no unique point on the IHA; any point on the IHA can be used to define the origin of the trihedron. The directions of the trihedron are defined using the acceleration and the angular velocity (Kirson, 1975): 4.40a cw - UVw mob e, - fl / lu X a| 4.40c eT - ek X 2“ The vectors em and eT define a plane, with ek being the vector normal to that plane. The direction of ek can be used for defining the ew- eT plane. If the motion is planar, then ek will lie in a plane parallel to the motion of the body. For general three dimensional motion, the derivative of ek will indicate the direction the ew- eT plane is turning toward. CHAPTER V A EINEMATICALLY DEFINED JOINT CENTER 5-1 WWW MAW—MAM In Chapter 4, a center point for the instantaneous helical axis was defined. The IHA center point can be defined as long as the angular acceleration has a component perpendicular to the angular velocity. 'The purpose of this chapter is to modify the instantaneous definition for a discrete series of helical axes. The center point for each helical axis can be defined in an instantaneous sense, if the angular velocity and angular acceleration of a rigid body are known, and if the velocity and acceleration of at least one point is known. However, two singularities exist. The first singularity occurs when the angular acceleration vector is parallel with the angular velocity vector: there is no unique center point on the IHA. .The second singularity occurs if the angular velocity vector is zero. Both of these singularities have theoretical solutions that can be defined in the instantaneous sense, but due to the need for higher order derivatives, and the errors associated with numerical differentiation of discrete data, the solutions to these singularities will be expressed in a discrete sense. The center point of an instantaneous helical axis can be approximated using two consecutive instantaneous helical axes. Figure 5.1 shows two consecutive helical axes, one at time t-i, the other at time t-i+l, where i and i+l are a separated by a finite amount of time. The center point of the helical axes at time t-i is simply the point on h the it helical axis that is closest to the helical axis i+1. As long as the two helical axes are not parallel, a unique point exists on the ith axis that is closest to the i+l axis. 40 Figure 5.1 Two consecutive instantaneous helical axes separated by a finite amount of time. If the position of each helical axis is known, the center point of helical axis 1 can be found. The next two sections consider the two singularities for which the center point cannot be defined. 513W If the angular acceleration of a rigid body is parallel to the angular velocity, or if the angular acceleration is zero, then no unique center point can be defined for the instantaneous helical axis. Consider the infinitesimally separated helical axes shown in Figure 5.2. The first axis is at time t-i, the last axis is at t-i+2. The time between each of the successive IHA's is infinitesimally small. The 42 “’ 1+2 is ”1‘1 ., ci+1 C t: ”1 Figure 5.2 Two consecutive parallel instantaneous helical axes followed by a third, non-parallel instantaneous helical axis. first two axes are parallel, the third axis is not parallel to the other two. The angular acceleration is parallel to the angular velocity for helical axis 1 (or the angular acceleration may be zero). The second derivative of the angular velocity must have a component that is perpendicular to the angular velocity, otherwise the third IHA would be parallel to the first two IHA's. Although a center point as defined in chapter four cannot be defined here, a higher order center point can be defined using the second derivative of u since it has a component perpendicular to 0. Since a and u are parallel, the velocity and the acceleration of all points on the ith IHA are the same, but the derivative of the acceleration (sometimes called the super acceleration, .43 or jerk) is not the same for all points on the IHA; thus a unique point can be found where the super acceleration is a minimum. This point can be used to define a higher order center point on the IHA. This argument can be extended to any number of derivatives of the angular velocity vector, until one is found that has a component perpendicular to o. This definition allows a center point to be defined for any series of parallel helical axes, followed by at least one helical axis that is not parallel to the others. This definition is also "forward looking", that is, only subsequent events affect the current definition of the center point; past events have no effect. A discrete version of this argument can be defined given any number of discretely defined helical axes. If the helical axes in Figure 5.2 .were discretely separated rather than infinitesimally separated then the center point of the ith axis can be defined as the point on axis i that is closest to axis 1+2. This is the discrete version of using higher order derivatives to define the center point. anW When the angular velocity is zero the IHA cannot be defined. However, if the angular acceleration is non-zero when the angular velocity is zero, L’Hopital's rule can be used to define the position of the IHA. Consider the equation for defining the shortest vector from a point to the helical axis: a X v 5.1 DP - —P 0'” where v? is the velocity of point P and DP is the shortest (perpendicular) vector from point P to the IHA. Let the angular velocity be zero at the time t-t then in the limit, equation 5.1 is 1; zero over zero. Applying L’Hopital's rule to the numerator and denominator: 44 a X v: + a X :2 5-2 lim 1) - lim 2m” tmti t-tti If v? is zero, or if a and v? are parallel, the above equation is still zero over zero, thus L'Hopital's rule can be applied again: & X v + a X + a X + a X . 5.3 lim 1) - lim 1 P 2;“. 3’2“." a? 8”] t-ot:1 t:-~t:1 |_ If v? is zero, or if m and VP are parallel, then equation 5.3 reduces _/ to: a X 8? 5.4 lim D - tvti Thus, in the limit, if the angular velocity is zero only instantaneously, the IHA can still 62 defined in terms of the acceleration and angular acceleration. The direction of the IHA axis is simply the direction of the angular acceleration vector. The above definition can be called the instantaneous axis of acceleration (IAA), that is, it is the axis of minimum acceleration. Such an axis can only be defined if the angular velocity is zero, and therefore this axis can only exist for the brief instant when a is zero. By defining the helical axis in terms of the angular acceleration, the definition of the center point must be shifted up one derivative higher, as discussed in the previous section. That is, a point of minimal super acceleration must be found on the IAA/ As with the center point, this point can only be defined if the super angular acceleration has a component perpendicular to the angular acceleration. However, by analyzing the center point of each IHA using the discrete definition given in section 5.1, it is only necessary to define the line of the IHA; the center point of that line can then be found using the discrete definitions given in section 5.1 and 5.1a. A second approach to this problem is to simply interpolate the center point across a region where u is small. That is, given that the 45 angular velocity is small (approximately equal to the error in locating the angular velocity vector), and therefore the error in determining the helical axis is large, a more straightforward solution would be to simply assume that the center point of motion is a smooth and continuous function. Thus, given the path of the center point before and after an interval of small angular velocity, the path of the center point can be estimated using standard interpolation techniques. t ven a es 0 D e He Given a series of discrete instantaneous helical axes that represent the continuous motion of a rigid body, this section will define an algorithm for defining the center point for each IHA. In the preceding sections, suggestions were made for defining the helical axis center point given discrete helical axes. Any algorithm for defining the center point must consider the special case of lines that are nearly parallel, and the case when the angular velocity vector is approximately equal to the error in determining the angular velocity. Define the following: “min is the smallest allowable angular velocityJ 0min is the smallest allowable angle between adjacent IHA’s. Angular velocities below wmin are considered to be highly influenced by the error in determining the angular velocity; thus the center point for these IHA's will not be determined directly, but will be determined by interpolation. The angle 0min assumes that two lines that have an angle between them that is less than or equal to 0min can be considered to be parallel lines. The following algorithm provides the simplest means of defining the center point of each instantaneous helical axis: Define a subroutine CP(e.,R ,e ,R.,R. ) where: 1 i j j 1c e is 1 R1 is e is J R3 is Ric 46 a unit vector on line i a point on line i a unit vector on line j a point on line j Define: is the point on line i that is closest to line j if - the last time frame in the data set. The input to the subroutine CP is (e1,Ri,ej,Rj); the output is the vector Ric’ the point on line i that is closest to line j. The algorithm is as follows: 1) 2) 3) 4) 5) i-O i-i+1 j-i J-j+1 If (mi < wmin) Ric - undefined, goto 2 J If ((01 < omin) goto 3 Steps 1 through 4 If eio eJ < cos(0min) goto 3 find a center point for axis 1, if axis 1 call CP(ui,Ri,u3,R3,Ri;) ‘ can be accurately defined. If (i< if) goto 2 If (R1c - undefined) Then find the first time frame, k, If (i - l) --> where ch is defined. Set all prior values of Ric - ch If (1 > 1 AND R1c is defined) Then Find the time frame before and after frame i where the center points are defined, that is, find 47 Ric and ch (2 < i < k) Then linearly interpolate between the two known center points to define the center points for the time frames 2 < i < k If (i > 1 AND Rfc is not defined) Then If the center points could not be defined for lines i through f then define the center point so that R Ric where Ric is the last defined center point. ic - End of Algorithm The above algorithm allows a center point to be defined for each instant in time as long as at least one line in the data set is not parallel to the other lines. The advantage of this algorithm is that it smooths the center point for lines that are approximately parallel and for lines where the angular velocity is small. This algorithm assumes that the center point moves in a linear fashion. If lines are approximately parallel, then the center point can be defined by the first line with a sufficiently large difference in its orientation and an angular velocity larger than a specified minimum. This algorithm will be tested on numerically simulated data and on human ankle joint data. The above algorithm is "forward looking", that is, except when forced to interpolate, and except when the lines are parallel at the end of a data set, the past history of the motion of a body is not used to influence the calculation of the center point position. In the next section a "best fit" solution for the center point will be described. 48 It has been shown that the center point of the instantaneous screw axis is the point of minimum acceleration of all the points on the line of the screw axis. In order to determine the center point, the acceleration and angular acceleration of the body must be known. However, since the data collected is discrete, it is possible to use two consecutive screw axes and find the point on the screw axis at time 1 that is closest to the screw axis at time 2. The velocity equations can be written in an overdetermined manner, providing a new approach to solving for the center point between two helical axes, or for an entire motion. The velocity of a point in space can be expressed relative to a point on the screw axis, as shown below: 5'5 "‘k"'k"ik"sk where V1 is the velocity of the point i, vs is the velocity of the point i that is parallel to the screw axis, R is the vector from the point i to any point on the screw axis (hence the negative sign at the begining of equation 5.5), and k is the time t-k. Let the point i be a point that is instantaneously at the origin of the coordinate system, then R is a vector from the origin of the coordinate system to the screw axis. Then, at time t-k+l, equation 22 becomes: 5 '6 '”k+1x I"k+1 ' vol<+l ' vsk+l Equations 5.5 and 5.6 can be used to find the point in space that has the minimal velocity for both time k and time k+1 by writing the equations in an overdetermined sense and solving for R as shown below. 49 5 . - .. - 7 9k R vok vsk Ek+1 vok+1 ' vsk+1 where vb is the velocity of a point at the origin for time k and time k+1. The solution to this overdetermined set of equations yields R that has the minimum velocity in a least squares sense for the two time frames. This point will lie on the shortest perpendicular line segment between the two screw axes. If the screw axes intersect at a point, then R will be the point of intersection. If the two lines are parallel, then R will be the shortest vector from the origin that satisfies the above equation. The significance of the above equation is that it can be used to find the point of minimum velocity for the entire event of interest, as shown below. 5.8 -[9] R - [Vo - VS] ‘31 where [ 9 ] - 92 u -n and vo1 - vsl - v [V _ V ] _ 02 52 o s v - v on sn Equation 5.8 finds the point of minimum velocity (in the least squares sense), for the entire event. In the case of a body fixed at 50 °r“! point, this method would provide a good numerical approach for firuiing the point the body is rotating about. The hip is often modeled as being fixed at one point, and (Lundberg et al, 1989) suggested that the ankle joint may also be accurately modeled as rotating about a fixed point. A variation of equation 5.8 is to weight the equation with the magnitude of the angular velocity. This would minimize the effect of "outlying" screw axes due to a small angular velocity magnifying the error. Equation 5.8 can be modified further to determine the path of minimum velocity for the event. Let the center point be modeled as a linear function of time: n 5.9 C(t) - k§1 Ak Fk(t) where A.k is a (3x1) vector and Fk(t) is any function of time, t. Thus equation 5.9 becomes: 5.10 - 91F1(0) 91F2(0) --- 91Fm(0) A1 92F1(h) 92F2(h) "' 92Fm(h) A2 A3 93F1(2h) 93F (2h) --- Q3Fm(2h) QnFl(nh) ghF2(nh) ... gth(nh) A'm where h is the time increment of data collection (assumed to start at time t-O). Equation 5.10 reduces to equation 5.8 if F1(t) - l and all other Fk(t)-0, in which case A1 - R. The above equation acts as a filter for the path of the center point. This equation can be modified by weighting both sides with the magnitude of the angular velocity, thus reducing the effect of "outlying" screw axes caused by error magnification for small angular velocities. If a polynomial is used for Fk(t) then: 51 91 0 0 ... 91Fm(0) A1 ‘32 92“ 92h2 ' ' ' 92%“) A2 5.11 - 93 9321: «33(211)2 93(2h)“'1 A3 - [vo- vs] 2 m-l _ _n gnnh 9n(nh) 9n(nh) ‘ _Am _ where the polynomial is of degree m. Equation 5.11 fits a polynomial to the path where the velocity of the points on the path is minimal for the entire event. This path is equivalent to a least squares fit path of the center points for each IHA. Thus equations 5.10 and 5.11 act as filters to the center point path. It must be noted that the solutions to the above two equations are "best fit" for the entire event, that is, past time frames affect the solution for future time frames and vice- versa. The advantage of this result is that measurement error can be smoothed. The disadvantage is that the results may be over-smoothed. CHAPTER VI APPLICATION TO THE ANKLE JOINT 6-1 W The purpose of this section is to apply the theory developed in this thesis to kinematically define the ankle joint center. Ankle joint kinematic data was collected at the Biomechanics Evaluation Laboratory (BEL) at St. Lawrence Hospital, Lansing, MI. Using state of the art video technology, it is possible to track the position of a target in three dimensional space. Thus, with the use of a force platform, both kinetic and kinematic data can be recorded for the lower limb,/ Four NEC cameras were used to track the position of targets placed on the foot and shank of two male subjects. The four cameras operate at 60 Hz with a shutter speed of 1/1000 seconds. The camera shutters are synchronized so that an image of the object is simultaneously recorded by each camera pixel screen, assuring that there is no time lapse between the images on each camera. Given the target position from two camera pixel screens, the three dimensional position of the target can be reconstructed if the camera positions are known in space. The camera positions are determined using a calibration structure targeted so that the target positions on the structure are known to a high degree of accuracy. Using a direct linear transformation (DLT) the camera positions can be defined relative to the known target positions of the calibration structure. The procedure to calibrate a space is as follows: the cameras are positioned according to the protocol of the study; once the cameras are positioned, the calibration structure is placed in the center of the space of interest. The video image of each camera must "see" each target on the structure. Thus, by knowing the position of the targets 52 53 °r| the calibration structure relative to a coordinate system on the calibration structure, it is possible to reconstruct the three dimensional position of each camera relative to the calibration coordinate system. A minimum of six calibration targets are needed in order to reconstruct the position of each camera. The calibration structure used for this experiment had sixteen targets, allowing the equations to be written in an overdetermined sense, thus increasing the accuracy of the calibration. The video system is able to "see" a target, as long as the target has a retroreflective surface,’ A light source is placed near the camera lens so that the retroreflective surface of targets in front of the lens reflect a bright image onto the pixel screen within each video camera. Since the target image is very bright, the threshold of the video system can be adjusted so that background images are removed and only the target images remain. The video image of on-off pixels can then be imported to a computer (SUN 4/26OC work station). The center of each target is approximated by finding the centroid of the pixel image. Once the centroid of the pixel image is found for each target, the three dimensional position of every target can be determined using a direct linear transformation (as long as each target has been viewed by at least two video cameras). The kinetic data for this experiment was measured by an AMTI force platform. The voltage output of the force platform was imported to the SUN computer via an analog to digital device. The force plate input and the video input were monitored so that forces measured by the force plate, that were above a given threshold, triggered an event marker on the video files, thus allowing the video and force platform data to be synchronized. 54 camera 1 camera 2 z 75cm amera 3 camera 4 C Tap View side View Figure 6.1 Camera set up, calibration space and force plate location The camera set-up, the calibration space volume and the force plate location for this experiment are shown in Figure 6.1. The calibration structure was positioned so that the force platform origin was at the center of the calibration space. The calibration coordinate system origin was aligned with the force platform coordinate system origin so that the position of the force plate relative to the calibration coordinate system was known. The cameras were set-up to allow tracking of targets placed on the foot and the shank. The targets were attached to the right lower limb of the subjects using hypoallergenic tape. The foot targets were placed on the posterior inferior aspect of the shoe over the calcaneous; the posterior medial aspect of the shoe below the medial condyle; and the posterior latBl‘tal aspect of the shoe below the lateral condyle. The shank targets were placed on the proximal shank, below the quadriceps tendon insertion; the distal shank; and the posterior shank on the tendon of the gastrochnemius, below the belly of the muscle. The target attachments are shown in Figure 6.2. From the three targets on the foot the position of the anatomical center of the ankle joint can be estimated. The anatomical center of the ankle was estimated to lie anterior and superior to the midpoint of N 01 ?\ 0 “ a. O -natom cal center of the ankle join Figure 6.2 Target placement on the shank and foot for each subject 56 the vector 5 (from targets E and F of the foot, as shown in Figure 6.2). A foot coordinate system (xf,yf,zf) was defined using the direction of vector EF as the yf direction; the cross product of A DE and DF as the zf direction and the cross of yf and 2f for the xf direction. The ankle anatomical center was then defined as: 6.1 Ac — fi+§F/2+hxf+ vzf Where h and v are the horizontal and vertical distances from the center point of vector EF to the anatomical center of the ankle joint (see Figure 6.2). This definition is similar to that used by Verstraete (1988).J The location of this anatomical joint center was compared to the locations of the kinematically defined ankle joint centers. The results from three trials from two male runners were investigated. The results in the following section are defined relative to the right shank coordinate system, which is defined as shown in Figure 6.2. The 2 axis is parallel to the long axis of the shank, the y axis points medially and the x axis points in the posterior to anterior direction. The shank coordinate system is located on the distal shank target, as shown in Figure 6.2. The next section provides a brief review of ankle joint biomechanics. The following sections investigate the angular velocity, the direction of the angular velocity vector, the location of the helical axis, the location of the center point of the helical axis, the location of the point of minimum velocity and the path of minimum velocity for the ankle joint of each subject. 57 6.2 W In much of the biomechanical literature, the ankle joint is modeled as a simple hinge joint between the shank and foot. Unfortunately, the ankle is far more complex than this simple model. The ankle joint consists of five bones: the tibia, fibula, talus, calcaneous, and the navicular. There are two major joints: the talocrural joint, also called the upper ankle joint (the joint between the talus and the tibia and fibula) and the talocalcaneonavicular joint, also called the lower ankle joint (the joint between the talus, the calcaneous and the navicular bones). Figure 6.3 shows the bones and joints that define the ankle joint. There has been extensive research into the biomechanics of the ankle joint. Barnett and Napier (1952) and Hicks (1953) found that the talo- crural joint has a different axis of rotation for dorsiflexion and plantarflexion. They found that for dorsiflexion the axis lies on a line that begins on the medial side of the ankle and proximal to the ankle down to a more distal point on the lateral side of the joint. On plantarflexion the axis was on a line that begins on the medial side but distal to the malleoli, and moves toward a point on the lateral side that is more proximal. Other researchers have investigated the axes of rotation of the ankle joint for motions other than dorsi-plantar flexion. Close (1956) described rotation about the vertical axis during walking. McCullough and Burge (1980) and Van Langelaan (1983) described the axes of rotation of the ankle joint for several different ankle motions. In general, the talc-crural joint allows rotation in all directions (Van Langelaan, 1983, Lundberg, et a1 , 1989) and behaves more like a ball and socket joint than the hinge joint many try to model it after. Right Foot Right Foot Medial View D Lateral View Joints Makin u the Ankle Joint Bones of the Ankle Joint A Talocrural (talus, tibia, fibula) 1. Tibia B. Subtalar (talus and calcaneus) 2. Fibula C Talocalcaneonavicular (talus. 3. Talus calcaneous and navicular) 4. Calcaneus D Calcaneocuboid (calcaneus and 5. Navicular cuboid) 6. Cuboid Figure 6.3 Bones and joints that define the ankle. 6.3 W All data presented in this section and the following sections are presented relative to the shank coordinate system. Table 6.1 shows the lengths and standard deviations of the lengths of the vectors between each of the foot targets (vectors DE, BF, and EF, using the target names shown in Figure 6.2) for both subjects and each trial. The data listed in the table are after smoothing and after being transformed into the shank coordinate system. This table provides an approximate estimate 59 for the variation in the length between relative targets due to skin motion and errors in the measurement of the target positions. The data listed in Table 6.1 are for stance phase only. The largest standard deviation for the foot data is 0.249 cm; the smallest standard deviation is 0.046 cm. Most of the vector lengths vary between 0.1 cm and 0.2 cm. Nine of the twenty-seven vector lengths have standard deviations less. than 0.110 cm. Figures 6.4 and 6.5 show the angular velocity vector of the foot relative to the shank, for both subjects. The angular velocity vectors show consistent inter- and intra-subject patterns. Rotation about the shank coordinate system y axis is the dorsi/plantar flexion axis, as shown in Figure 6.2, with plantar flexion being positive (using the right hand rule). Rotation about the x axis is inversion-eversion of the foot with inversion being positive; and rotation about the z axis is medial and lateral rotation of the foot with medial rotation being positive. Soutas-Little, et al (1987) applied an Euler angle approach to describe ankle joint position, angular velocity and acceleration. Although the Euler angle axes are not orthogonal, it is still possible to compare the results shown in Figures 6.4 and 6.5 to the angular velocity results found using Euler angles. The angular velocities reported here agree with the findings of Soutas-Little and co-workers. At heel strike the foot is dorsiflexing, everting and laterally rotating. By midstance the angular velocity vector is near zero; between midstance and toe off the foot is plantar flexing, inverting and medially rotating. RADIANSPERSEOOND 60 Table 6.1 Length and standard deviation of vectors between targets of the foot in the shank coordinate system, during agance phase. Subject & 5E (cm) (cm) (cm) (cm) t dev th td ev 8.068 0.162 8.394 0.058 8.132 0.178 8.672 0.092 8.196 0.193 8.788 0.190 8.235 0.177 8.726 0.138 8.251 0.122 8.746 0.119 8.229 0.149 8.677 0.140 Foot: 'l'arvetg TIME IN SECONDS Figure 6.4a (cm) en b 9.961 9.957 10.027 10.288 10.131 10.190 X axis (inversion(+)/eversion(-)). (cm) std dev 0.102 0.105 0.109 0.249 0.206 0.068 - TRIALI - TRIAL! - TRIALS 0 rhuaww :0 T000" 5 61 .. o - TRIALI ‘ TRIAL2 - TRIAL3 0 HOOISWO RADIANSPERSEOOND 'bbhbbLGhAo—nueuouoo .. 0 TIME IN SECONDS Figure 6.4b Y axis (plantar(+)/dorsi(-)). - TRIALI - TRIAL: - TRMLII I I flJ<> HoeIStrIIo D finch HADIANSPERSECDND TIME IN SECONDS Figure 6.4c Z axis (medial(+)/lateral(-) rot ). Angular velocity of the ankle joint for Subject A RADIANSPERSEOOMJ RADIANSPERSEWM) ébbhbbbbbéo—nueoou O b b L TIME N SECONDS Figure 6.5a X axis (inversion(+)/eversion(-)). I... f" 0 TIME IN SECONDS Figure 6. 5b Y axis (plantar(+)/dorsi(-)). / mum TRIALZ TRIAL: I0 HoolSInko ‘0 TooOII " TRIAL! - TRIALZ - TRIALS 0 Ikdfiflw U T000" ,1 63 l — a _. g - - TRMLI 1 __ -TRW.2 - TRMLG I E o y I I I I0 Ibdaflu g .. __ V? a mo" 4 .. .3 __ '4 —— J -— nuemsecouos Figure 6.5c Z axis (medial(+)/lateral(-) rot.). Angular velocity of the ankle joint for subject B 6.4 e a s at on The helical axis location can be expressed in a number of ways. One means is to show the points of intersection of the axes with a given plane. However, the axes of the ankle joint vary to such a degree that there is difficulty defining a meaningful anatomical plane such that none of the axes lie parallel to the plane. If any of the axes of interest lie parallel or nearly parallel to the plane, the points of intersection with that plane can be largely dispersed. An alternative means of expressing the position of the helical axes closest to all the helical axes is to find the point in space that is 64 and plot the shortest (perpendicular) vector from this point to each helical axis. The point in space closest to all the helical axes was called the helical axes centroid (HAC). If all the axes cross near a point, the (shortest) vector from this point (the HAC) to each axis will provide a concise description of the dispersion of the helical axes. In this investigation the helical axes centroid (HAC) was found for the stance phase for each trial of both subjects. The shortest vector from the HAC to each IHA was then determined. These vectors provide a means by which the variation in the helical axis positions can be measured. If the vector from the HAC to each helical axes has a very small magnitude, then the motion in question is most likely motion about a fixed point. Figure 6.6 shows the shortest vector from the HAC to each IHA for subject A, trial 3. In Chapter 5 it was found that the IHA can be determined using the acceleration of a point and the angular acceleration of the body, if the angular velocity vector has zero magnitude. When this occurs, the axis is called the instantaneous axis of acceleration (IAA). The IAA can only be defined if w is zero; however, in this chapter, that definition is modified and is used to locate the IHA when w is small. Since the IAA definition assumes that the angular acceleration is large when the angular velocity is small, a definition of "small" angular velocity is necessary. From equation 5.4, if luj is zero, then the shortest vector from a point D to point P on the IAA is: 6.2 HIS-2:0 If the angular velocity is non-zero, then the acceleration of point D with respect to point P is: 6.3 aD-aP-aXD—P- uX(uXDP) 65 where a? must be parallel to a by the definition of equation 6.1. Thus the error in UP due to a non-zero angular velocity, is: '2 2 __ IUI _ Ifll 6.5 |DP| - |DP| + m |DP| - |DP|(1 + ) Therefore, if wz/a is small, equation 6.2 can be used to define the IHA. Figure 6.6 shows the shortest vector from the HAC to each helical axis for subject A, trial 3, comparing the results obtained from the standard definition for calculating the IHA position (equation 3.19) to the results obtained if the IAA (equation 6.2) is substituted for the IHA for wz/ a s 0.06. A clear improvement is obtained in the continuity of the position of the helical axes if the IAA is used to determine the IHA for small angular velocities. 8 6 -IIIAI 4 -Im: 0M“ 01000! 2 ° . . . I . . -fi' . . I. L."'.. \. -2 META”! Figure 6.6a The shortest vector from the HAC to each helical axis (X direction) comparing the IHA to the IHA/IAA for subect A, trial 3. 1.0 ..m, -IINIAA1 0mm ten 0 A I I I I I J!“ a '. I. I. I. I. I. I. .. -1.0 Peacemrmcca Figure 6.6b The shortest vector from the HAC to each helical axis (Y direction) comparing the IHA to the IHA/IAA for subject A, trial 3. Figure 6.6c The shortest vector from the HAC to each helical axis (Z direction) comparing the IHA to the IHA/IAA for subject A, trial 3. 67 Figures 6.7 through 6.12 show the helical axes of the foot relative to the shank coordinate system origin. The IAA is substituted for the IHA for small wz/ a. The helical axes are projected onto the X-Z (sagittal) plane and the Y-Z (frontal) plane. The projection of the helical axes begins at the point on each IHA that is closest to the helical axes centroid (HAC). The direction of the angular velocity vector is indicated by the arrow on each projection. The IHA labeled "1" in each of the Figures 6.7 - 6.12 describes the motion of the foot relative to the shank at heel contact. Midstance occurs at IHA 6 for subject A and at IHA 7 for subject B. Toe off occurs at IHA 12 for subject A and IHA 13 for subject B. A foot and shank have been drawn in to aid in visualizing the positions of the helical axes. The black rectangle in each of Figures 6.7 - 6.12 shows the location of the anatomically defined center of the ankle (from equation 6.1). In the X-Z (sagittal) plane, the axes are all nearly parallel to one another. The slope of the axes indicates that the rate of medial/lateral rotation and the rate of inversion/aversion are approximately equal for most of the stance phase of running. In the Y-Z (frontal) plane, the axes describe almost a complete circle. The pattern of the helical axes appears to have three stages, corresponding to heel contact, midstance, and toe off. At heel contact the dominant axial direction is that of lateral rotation and eversion. From heel contact to midstance (IHA 1-6 or 1-7 for subject A and IHA 1-7 or 1-8 for subject B) the direction of the axes change so that the dominant axis of rotation is increasingly about the dorsiflexion axis; with the foot dorsiflexing until midstance. At midstance there is a sudden change in the direction of the axes as the foot motion switches from dorsiflexion to plantar flexion. At midstance the helical axes are 68 suuflx that medial rotation and inversion are the dominant rotations, since the dorsi- plantar-flexion angular velocity is zero, or near zero, at midstance. From midstance to toe off the direction of the axes changes until they are nearly parallel with the plantarflexion (+Y) axis. It is obvious from Figures 6.7 through 6.12 that a planar model for ankle motion is inadequate as a description of the kinematics of the ankle during running. The axodes for subject A all cross near a single point (the HAC) in the ankle joint. The HAC lies posterior to the Z axis of the shank and slightly medial to the midline of the tibia. Table 6.2 lists the location of the helical axes centroid (HAG) for each trial of both subjects (relative to the origin of the shank coordinate system). The HAC's listed in Table 6.3 all substitute the IAA for the IHA for small wZ/a (typically < 0.10). There is relatively high intra-subject consistency for the location of the HAC. Subject B has a greater dispersion of the helical axes compared to subject A. In addition, the location of the HAC differs slightly for subject B; it is closer to the midline of the tibia, except for trial 3, where it is slightly medial to the tibial midline. For most of the trials, the helical axes lie within a sphere with a 1.5 cm radius from the HAC. Many of the helical axes lie within a sphere of 1 cm radius. Several researchers have investigated the finite helical axes of the ankle joint (Hicks. 1953; Close, 1956; Lundberg, et a1, 1989; Van Langelaan, 1983). However, it is not possible to use this literature to compare the findings of this thesis. Besides the differences between the finite helical axis and the instantatneous helical axis, the joint loading and kinematics differ substantially between the above studies and this study. 6 95.10 a no 7 n #7 I 12 -aa -10 z A 1 no" 2 34 1 q‘.\\ ll ‘ 5 C) ' : =: 3 -aa -10 IO a 9 no 11 12 IO - 5 ===—~ ve # Heel Contact 1 4 Midstance 6 Toe off 12 3 2 Figure 6.7 l X-Z (top) and Y-Z (bottom) projection of the IHA's (or IAA where appropriate) for subject A trial 1. " .—_==12 4 # Heel Contact l Midstance 6 3 Toe off 12 Figure 6.8 X-Z (top) and Y-Z (bottom) projection of the IHA's (or IAA where appropriate) for subject A trial 2. Ion _ I , . f ”I", V t # 557:?" V"; Heel Contact 1 s != Midstance 6 4 Toe off 12 3 2 ‘ Figure 6.9 X-Z (top) and Y-Z (bottom) projection of the IHA's (or IAA where appropriate) for subject A trial 3. 72 /“;o.13 , jg, \ ' J ~23 -I 1 2 . 3 4 s -28 10 n w-- 2 13 ml. ’ HM - ‘ ua“":«v ‘ . ”(fii' ' Event IHA # ( Heel Contact l Midstance 7 5 Toe off 13 4 3 2 Figure 6.10 X-Z (top) and Y-Z (bottom) projection of the IHA's (or IAA where appropriate) for subject B trial 1. » $34” “"“ 1" ve Heel Contact Midstance Toe off I 3 2 Figure 6.11 X-Z (top) and Y-Z (bottom) projection of the IHA's (or IAA where appropriate) for subject B trial 2. l3 / $12?” 5.1-r29 >1 the term (l/f)2 can be considered small compared to the (AR)f terms, thus the magnitude of the velocity error is: i+l i-l A1.4 |Av| - [Jinn - AR )|f s 12‘(|AR| + [AR|)f - ARf where AR is the maximum measurement error for the collected data. A1.2 W; The angular velocity of a rigid body can be calculated using the method of Verstraete and Soutas-Little (1990), as discussed in section 3.5. With this method, the least square angular velocity can be found if at least three non-colinear targets are attached to'a rigid body and the positions and velocities of these targets are known at any given time. The vector between any two targets on a rigid body and the relative velocity between any two targets are: A1.5 Dij - Ri - Rj (i f j, i < j) and Al 6 vij - v1 - v3 The velocity vector vij must be perpendicular to the vector Dij due to the rigid body assumption. The relative velocity, vij can be expressed in terms of the angular velocity vector and the relative position vector: ij - u x Dij 99 01': A1.8 [Pljlu-vlj I where [913] is a skew symmetric matrix. Equation A1.8 can be expanded so that all the vectors Dij and vij are included. This is an overdetermined set of equations, that can be solved in the least squares sense: Al.9 [p] u - [v] T T T T T T where [D] - [912... 91“ 923... PZn - Pn-1,n] and [v] _ [v1 v; VI v1 VI 1T 12"' ln 23"' Zn "' n-1.n where n is the number of targets on the body. The solution to the above set of equations is: A1.lO w - IIQITIQII‘IIPI Iv] The error in a due to errors in [P] and [v] can be found with the aid of the relationship (Stoer and Bulirsch, 1983): A1.ll (I + If1 a II - P) where I is the identity matrix and F is "small". Thus if D - D + AD where D is the true vector between each target and AD is the error in vector between each target, then: A1.12 {[9 + A9]T[D + AQII'l - (~TD + _TA9 + A~TP + A_TA9I’1 - I_Tg + _TA9 + AyTpI '1 . l_T9[I + (-T9)'1<~TA9 + A9T9>II'1 100 = II - <9T9>'1<9TA2 + A2T2>I<9T9r1 In equation Al.12 it is assumed that APTAQ a 0. Using the approximation from Al.12 the error in the angular velocity vector can be estimated: Al.13 u + A» - [I-(9T9)'1(EFA9 + Apr)](pr)‘1 [p + A9]T[v + Av] so A0 is approximately equal to: A114 Aw - <9T951<9Tm+<1f1331 - <9T9)‘1I9TA9+A2T21<2T9>'1<9Tv) The magnitude of the error (using the euclidean norm) is: T T T T T T Al.15 |Au| _<_ “21“”+ IAPVI . 12 “2| |va + IA!) 2| In vI IPTPI IPTPI IPTPI the term Av can be estimated by using equation (Al.4) (Av = AR f = AD f) and the term AD can be estimated by finding the maximum error in D by the standard deviation in the length of each Dij for the entire data collection. In this manner, an estimate can be obtained for the magnitude of the error in the angular velocity for each time frame. The right hand side of equation Al 15 can be approximated as follows: Al.16 DH f AD + ADM vM + 2ADvM _ AD[ _£_ + 3 v” I 2 2 ”M ”H OH where the subscript M denotes the maximum value of the vectors D and v, AD is the maximum error in the matrix [9]. Equation Al.16 is meant only 101 to provide an approximation of the maximum error in the angular velocity. Equations Al.15 and Al.16 reveal that the error in the angular velocity is directly dependent on the accuracy of the target position measurements and that the error is amplified by the frequency of data collection and the velocity of the rigid body. At low velocities the error in the angular velocity is controlled by the measurement error, for higher velocities the effect of the velocity term on the error is only significant if DM is small. By maximizing the distance between targets, the error in the angular velocity vector can be minimized. A1.3 Brr2r_EatLsats_ia_Qalsslatins_£he_§£rez_sai§ The instantaneous screw axis can be calulated from a point in space if the velocity at that point is known: Al.17 r - ” x ‘1 0'9 where P is the vector from point i to the IHA and v is the velocity of 1 point i. The error is then easily calculated: (u + An) x (vi + Avi) (u + Au) - (u + Au) Al.18 P + AP - Thus the error in determining the position of the ISA (dropping higher order terms) is: Al.19 AP - w x Av: + A” x Vi a a X AvI + Au x vi a - u + 2 w-Au u - u where it is assumed that |Aw| << lwl. Thus: lwIIAvil + IAwllvil 2 A1.20 [API 5 102 The error can be minimized if the target that is being used to find the ISA is close to the screw axis and therefore has a minimal velocity. Equation A1.20 also indicates that the position of the screw axis can be more accurately measured for higher angular velocities. A1.3 Brr2r_in_the_§enter_221nt_rgsitisn The position of the center point can be found given the direction of each IHA and a point on each IHA. Let 21 and 32 be the direction vectors for two consecutive IHA's. Let 31 and R2 be vectors from the origin to each of the IHA's, as shown in Figure Al.l. The angle between the two lines is found from the dot product of $1 and £2: A1.22 eloe2 - coso Define: A1.23 P - 22 - 21 A A A1.24 e3 - e1 X e2 |e1 X e2] A1.25 ea - e3 X e1 The point on line 1 that is closest to line 2 can be found as follows: A1.26 q - - P . e4 e1 + Poe e tan(0) 1 A l where q is a vector from the origin to the point on line 1 closest to line 2. The error in the angle 0 can be determined by differentiating equation A1.22: A1.27 -sin0 d0 - deloe2 + e1~de2 or in terms of A9 the error is: A1.28 A0 - '[A;1';g + ;10A;2] sinfi 103 The error in the position of q due only to errors in the location of 21 and R2 is determined by differentiating equation A1.26: . .e . .e A A A A A A1.29 dq - de - 552 4 1 4 e1 + dR2-e1 e1 - de-e1 e1 tano or: A A .e A A A .; A Al.30 Aq - A21 - AR1oe1e1 + cinch e1 + Akzcele1 - tinofi e1 The magnitude of the above error is: Al.31 |Aq| s ARI + ARI-e1 + ARzoel + Afi'ea + 22W. tano s 21ml .+ AR2 + AR] + AR; tan0 equation Al.3l indicates that the error in q increases if the angle between line one and two is small. From equation A1.28 it can be seen that if two lines are approximately parallel, the error in the orientation of each line is greatly magnified, making the error in the angle between the two lines much larger. 0 Figure A1.1 Vectors for defining the center point position. APPENDIX II NUMERICAL SIMULATION AZ-l W The main purpose of this section is to investigate the accuracy of finding the instantaneous helical axes and the center point of the helical axes given noisy measurement data. The model used in this simulation was designed so that the orientation of the instantaneous helical axes changed through time. The accuracy of assessing the helical axis position and orientation must also involve variations in the magnitude of the angular velocity vector, including the case of zero angular velocity. The data for this simulation was generated from the kinematic equations that describe the model shown in Figure A2.1 The model consists of a wheel on an axle, with the axle fixed to a pivot. The axle rotates about the pivot with an angular velocity 0(t). A coordinate system is defined so that the z axis is along the axis of the pivot, as shown in Figure A2.1; thus 0(t) lies along the z axis. The wheel rolls without slip on the x-y plane. The angular velocity of the wheel is a function of 0(t), the length of the axle, RA' and the radius of the wheel, R. The IHA’s generated by the motion of the wheel describe a conically shaped surface. Figure A2.2 shows the shape of the ruled surfaces generated by the fixed and moving axodes of the wheel IHA's. The motion of this model is constrained to motion about a fixed point; thus all IHA’s must pass through the fixed point (0,0,10) as shown in Figure A2.2. The motion of this model can be defined mathematically, and therefore the positions of points on the wheel can be used to generate data for simulating rigid body motion. This data can then be used to 104 105 0(t) Tz'Ta RA Y > T1 Figure A2.1 ml x v Model for generating the simulated data. Moving axode ruled (0,0,10) Figure A2.2 Ruled surfaces generated by the fixed and moving axodes. 106 determine the accuracy of the theory presented in this thesis. The exact positions of points located on the wheel were defined as functions of 0(t), RA, and R. The velocity, angular velocity, acceleration, and angular acceleration can be determined exactly, for a given 0(t), R , R. Random errors were introduced to the known point positions to simulate measurement errors and soft tissue motion. The following criteria were used to define the motion of this model: A2.1a 0(t) - 4.0 sin(2wt) (0.0 S t 5 0.98333 seconds) A2.1b RA - 20.0 cm A2.1c R - 10.0 cm A2.1d T1(0) - (20.0, 0.0, 2.0) A2.1e T2(0) - (20.0, -6.928, 14.0) A2.1f 13(0) - (20.0, 6.928, 14.0) Where T1(0) are the initial positions of each target at time t-O. Three targets are the minimum necessary for defining three dimensional rigid body motion; thus this simulation used three points distributed equally about the surface of the wheel as shown in Figure A2.1. These points were used to simulate targets. The position of each of the three simulated targets can be determined mathematically; thus the exact position for each simulated target can be generated for a given 0(t). The target positions were determined for the values of 0(t), RA’ R and the initial simulated target positions given in equation A2.1. A random error bounded within i2.0mm was introduced to the known position of each point. The following sections investigate the accuracy of the angular velocity, the angular acceleration, the position of the instantaneous helical axes, and the position of the center point of the helical axes, comparing the data with the 12.0mm random noise to the exact known solutions. 107 The simulated data consists of 60 "frames"; however, due to the ' effects of smoothing and differentiating numerical data, the beginning and end points of the data set are unreliable. Therefore, the first five frames and the last five frames of the simulated data will not be used in the comparison to the true data. In this section, the true angular velocity, angular acceleration, the location and position of the instantaneous helical axes, and the location of the center point of the helical axes will be compared to the solutions obtained from the noise induced data. The noise induced data was smoothed and differentiated using a quintic spline smoothing routine (Woltring, 1986). The velocity and acceleration were calculated for each of the three targets from the derivatives of the quintic spline. The angular velocity was determined using the method of Verstraete and Soutas-Little (1990). The true angular velocity was compared to the noisy data angular velocity. The differences in the measured angular velocity and the true angular velocity are small. One means of quantifying the difference between two vectors is to find the root mean square error (RMSE) between the vectors. The RMSE between the true angular velocity and the angular velocity calculated from the noise induced data is 0.218 rad/sec; the maximum magnitude of the angular velocity was 9.0 rad/sec. The RMSE between the true angular acceleration and the angular acceleration calculated from the noise induced data was 3.151 rad/secz; the magnitude of the angular acceleration was as high as 50 rad/secz. The two primary concerns for the helical axes are the position and ggigggagign of the lines of the axes. First, the orientation of the angular velocity vector calculated from the noise induced data was compared to the true angular velocity vector. The error in the 108 orientation between the true, and noise induced, data angular velocities were plotted against the magnitude of the true angular velocity, as ' shown in Figure A2.3 (the direction of the true angular velocity for |"truel - 0 was defined by the direction of the true angular acceleration vector). For an angular velocity above 4 radians per second, the error in the orientation between the true and error induced angular velocity vectors is less than 3 degrees. For angular velocities above 3 radians per second, the orientation error is less than 6 degrees. For 2 radians per second, the orientation error is less than 8 degrees; and for 1 radian per second, the error is less than 12 degrees. Thus, for large angular velocities ( > 4 rad/sec), lines can be considered parallel if they are within 3 degrees of each other. For small angular velocities ( < 2 rad/sec), it is not possible to determine the angle between lines any better than 12 degrees. The positions of the noise induced IHA's were compared to the true IHA position. The shortest vector from the origin to each IHA (for each instant in time (every 1/60 second)) was calculated, for both the true and the noise induced data. The differences in these two vectors were then used as a measure of the error in the IHA position. The RMSE in the position of these two vectors, and therefore in the IHA position, was found to be 0.4241 cm. In Chapter 5, a means was developed for calculating the IHA for an angular velocity of zero; this method was expanded in Chapter 6, equation 6.3, to determine the IHA for small angular velocities. The instantaneous axis of acceleration (IAA) was substituted for the IHA for the noise induced data (for w < 2.0 rad/sec), in order to determine the position and orientation of the IHA more accurately. When this was done, the RMSE between the true and noise induced IHA positions dropped DEREES DEGREES O‘HHAOOQOI 0"“...VIC. 109 ANGJLAR VEOCITY ”NUDE Figure A2 . 3 Angle between ”true and "noisy I I I I I I I l l I I l I 1 a a s 7 s s ANGULAR VELOCITY MAMlTUDE (RAD/SEC) Figure A2 . 4 Angle between "true and ”noisy' substituting the IAA for the IHA for w < 2.0 rad/sec. 110 to (3.2925 cm. Thus, a significant improvement in the IHA position was achieved by simply substituting the IAA for the IHA for small angular velocities. In addition to improving the position of the IHA, an improvement in the orientation of the IHA was also obtained by substituting the IAA for small angular velocities. Figure A2.4 shows the improvement in the orientation of the noise induced angular velocity direction relative to the true angular velocity direction. The error in the orientation improved dramatically, decreasing the error by as much as 10 degrees. The center point of the motion was calculated using three different methods. The first method used the algorithm in Chapter 5 to determine the center point of each IHA. For the algorithm, 0min and mm are in required. The minimum angular velocity magnitude was ”min - 0, due to the substitution of the IAA for the IHA for small angular velocity, no minimum angular velocity vector was needed. For the minimum angle between consecutive helical axes, several different minimum angles were investigated. It was found that the larger the angle for 9m the in greater the accuracey of the results. In the table below, the results are given for 0min - 15' and 30°. The center point was also calculated by determining the minimum velocity center using equation 5.8. The third method used to find the kinematic center of motion was a linear (1° polynomial) model to determine the path of minimal velocity (equation 5.11). The true center of motion for the simulated data is located at the point (0,0,10) in the coordinate system shown in Figure A2.1. The RMSE between this true center point and the center point(s) calculated from the above three methods are shown in Table A2.1. 111 Table A2.1 Root mean square error between the calculated center points and the true center point. Method Used to Determine the Center Point RMSE Algorithm 1.205 cm 0 - 15° min Algorithm 0.622 cm 0 - 30° min Min Vel Point 0.160 cm Min Vel Path 0.126 cm The error in the center point when using the linear path of minimum velocity or the point of minimum velocity are extremely small relative to the induced error in the original data. The error in the center point calculated using the algorithm of Chapter 5 was larger. The reason for this increase in error is the l/tan0 term in the equation for calculating the position of the center point on each IHA (see equation 3.7). For small angles, the error is magnified by the l/tanfl term. The error in locating the position of the IHA was previously shown to be 0.2925 cm. The algorithm in Chapter 5 requires that a minimum angle and a minimum angular velocity be entered by the user in order to determine: i) if two lines are "parallel" and ii) if the axis under consideration is stable enough to be used to define a center point. The angle used for this algorithm was 30°, a significantly high angle, but still having a 1/tan9 of 1.73. 112 Therefore the error in the center point will be almost double the error in the position of the IHA. The path of minimum velocity and the point of minimum velocity are more accurate than the results obtained from the algorithm. These two methods are therefore preferable for calculating the kinematic center of a biological joint. BIBLIOGRAPHY Ball, R. (1900). A Treatise on the Theory of Screws, Cambridge University Press. Blacharski, P. A., Somerset, J. H., Murray, D. G. (1975). A Three- Dimensional Study of the Kinematics of the Human Knee. Journal of Biomechanics, Vol. 8, pp 375-384. Blankevoort, L., Huiskes, R., Lange, A. D. (1990). Helical Axes of Passive Knee Joint Motions. Journal of Biomechanics, vol 23, pp 1219-1229. Chalses, M. (1830). Note sur les Proprietes Generales du Systeme de deux corps semblables entr’eux at places d’une maniere quelconque dans l’espace; et sur le deplacement fini ou infiniment petit d'un corps solide libre. Ferussac, Bulletin des Sciences Mathematiques, vol. 14, pp 321-326. Chao, E. Y. S. (1990). Issues Related to Joint Moment Calculation. A discussion on Biomch-L, a biomechanics electronic bulletin board. The discussion is retrievable by sending the command: send Biomch-L log 9011 (on one line) to LISTSERV@HEARN.BITNET, a BITNET electronic mail service. Chao, E. Y. S., An, K. N. (1982). Perspectives in Measurements and Modeling of Musculoskeletal Joint Dynamics. In Biomechanics: Principles and Applications, Huiskes, R., Van Campen D., and De Wijn, J., eds. pp 1-19. Clauser, C. E., McConville, J. T., Young, J. W. (1969). Weight, Volume and Center of Mass of Segments of the Human Body. AMRL Technical Report, Wright-Patterson Air Force Base, Ohio. Clifford, W.K. (1873), Preliminary Sketch of Bi-quaternions, Proceedings of the London Mathematical Society, vol 4, p 361. Close, J. R. (1956). Some Applications of the Functional Anatomy of the Ankle Joint. Journal of Bone and Joint Surgery (American), Vol. 38A, pp 761-781. Dimentberg, F.M. (1948). A General Method of Investigation of Finite Displacements of Three-Dimensional Mechanisms, and Certain Cases of Passive Couplings. Transactions of a Seminar on the Theory of Machines 113 114 amui Mechanisms, USSR Academy of Sciences Institute of Mechanical Engineering, No. 17 Dimentberg, F.M. (1965), The Screw Calculus and its Applications in Mechanincs, Izdatel'stvo, "Nauka", Moscow, USSR, (English translation: AD6890993, Clearinghouse for Federal and Scientific Technical Information). Fiorettei, S., Jetto, L., Leo, T. (1990). Reliable In Vivo Estimation of the Instantaneous Helical Axis in Human Segmental Movements. IEEE Transactions on Biomedical Engineering, vol. 37, pp 398-409. Hamilton, W. R. (1830). First Supplement to an Essay on the Theory of Systems of Rays. Trans. of the Royal Irish Academy, vol. 16, pp 4-62. Hamilton, W. R. (1845). On Some Additional Applications of the Theory of Algebraic Quaternions. Royal Irish Academy Proceedings, vol. 3. Hamilton, W. R. (1848). Some Applications of Quaternions to Questions Connected with the Rotation of a Solid Body. Royal Irish Academy Proc. vol. 4, pp 38-56. Hicks, J. H. (1953). The Mechanisms of the Foot.I-the joints. Journal of Anatomy, Vol. 87, pp 345-357. Hunt, K. H. (1967). Screw Axes and Mobility in Spatial Mechanisms via the Linear Complex. Journal of Mechanisms, vol. 3, pp 307-327. Isman, R. E., Inman, V. T. (1969). Anthropometric Studies of the Human Foot and Ankle. Bulletin of Prosthetic Research, Vol. 10/11, pp 97-129. Karlsson, J.0.M., Murphy, M.C., Mann, R. W. (1990). Using Axodes to Compare In Vivo Knee Kinematics Measured with Bone- and Skin-Mounted Markers. ASME Winter Annual Meeting, Symposium on Dynamics and Control of Biomechanical Systems: Musculoskeletal Modeling and Control, November 25-30, 1990. Kirson, Y., Yang, A. T. (1978). Instantaneous Invariants in Three Dimensional Kinematics. Trans. of the ASME. Journal of Applied Mechanics, Vol. 45, pp 409-414. Lundberg, A., Svensson, 0. K., Nemeth, G., Selvik, G. (1989). The Axis of Rotation of the Ankle Joint. Journal of Bone and Joint Surgery vol 7l-B, pp 94-99. Manter, J. T. (1941). Movements of the Subtalar and Transverse Tarsal Joints. Anat. Rec. Vol. 80, pp 397-410. Mobius, A. F. (1837). Lehrbuch der Statik, Leipzig [English abstract available in Ball, 1900]. Mobius, A. F. (1838). Ueber die Zusammensetzung unendlich kleiner Drehungen. Crelle's Journal, vol. 18, pp 189-212 [English abstract I ”fa—m L 115 available in Ball, 1900]. Plucker, J. (1865). On a new Geometry of Space. Phil. Transactions, vol. 155. PP 725-791. Plucker, J. (1866). Fundamental views Regarding Mechanics. Phil. Transactions, vol. 156, pp 361-380. Poinsot, L. (1806). Sur la Composition Des Moments et la Composition des Aires. Journal de l'Ecole Polytechnique, vol 6, pp 182-205 [English abstract available in Ball, 1900]. Poinsot, L. (1851). Theorie Nouvelle de la rotation des corps. Liouvelle's Journal Math, Vol 16, pp 9-129,289-336 [English abstract available in Ball, 1900]. Procter, P., Paul, J. P. (1982). Ankle Joint Biomechanics. Journal of Biomechanics, vol. 15, pp 627-634. Prodromos, C. C., Andriacchi, T. P., Galante, J. 0. (1985). A Relationship Between Gait and Clinical Changes Following High Tibial Osteotomy. Journal of Bone and Joint Surgery (American), Vol. 67A, pp 1188-1194. Selvik, G., Alberius, P., Aronson, A. S. (1982). A Roentgen Stereophotogrammetric System: Construction, Calibration, and Technical Accuracy. Acta Radiological Diagnosis, Vol. 24, pp 343-352. Smidt, G. L. (1973). Biomechanical Analysis of Knee Flexion and Extension. Journal of Biomechanics, Vol. 6, pp 79-92 Sommer, H. J., Buczek, F. L. (1990a). Experimental Determination of the Instant Screw Axis and Angular Acceleration Axis. IEEE 16th Annual N.E. Bioengineering Conf., pp 141-142. Sommer, H. J., Buczek, F. L. (1990b). Least Squares Estimation of the Instant Screw Axis and Angular Acceleration Axis. ASME Winter Annual Meeting Advances in Bioengineering, Dallas, Texas, November 25-30, pp 339-342. Spoor, C. W., Van Leeuwen, J. L., Meskers, G. M., Titulaer, A. F., Huson, A. (1990). Estimation of Instantaneous Moment Arms of Lower Leg Muscles. Journal of Biomechanics, vol. 23, pp 1247-1259. Study, E. (1903),Geometrie der Dynamen [Geometry of the Dynamics], Verlag Teubner, Leipzig. Van Langelaan, E. J. (1983). A Kinematical Analysis of the Tarsal Joints. Supplementum 204, Vol 54, ISBN:87-l6-06251-5, ISSN: 0300-8827, Acta Orthopaedica Scandinavia, Munksgaard, Copenhagen Veldkamp, G. R. (1963). Curvature Theory in Plane Kinematics. J. B. Wolters, Groningen. 116 ‘Veldkamp, G. R. (1967a). An Approach to Spherical Kinematics Using Tools Suggested by Plane Kinematics. Journal of Mechanisms, Vol. 2 pp 437-450. Veldkamp, G. R. (1967b). Canonical Systems and Instantaneous Invariants in Spatial Kinematics. Trans. of the ASME, Journal of Engineering for Industry, Vol. 89, pp 84-86. Verstraete, M. C. (1988). A Method for Computing the Three-Dimensional Forces and Moments in the Lower Limb During Locomotion. Ph. D Dissertation, Michigan State University, East Lansing, MI. Verstraete, M. C., Soutas-Little, R. W. (1990). A Method for Computing the Three-Dimensional Angular Velocity and Acceleration of a Body Segment from Three-Dimensional Position Data. Journal of Biomechanical Engineering, Vol. 112, pp 114-118. Walker, P. S., Shoji, H., Erkman, M. J. (1972). The Rotational Axis of the Knee and its Significance to Prosthesis Design. Clinical Orthopeadics and Related Research, Vol. 89, pp 160-167. Wang, J. W., Kuo, K. N., Andriacchi, T. P., Galante, J. O. (1990). The Influence of Walking Mechanics and Time on the Results of Proximal Tibial Osteotomy. Journal of Bone and Joint Surgery (American), Vol 72A, pp 905-909. , Weber, W., Weber, E. (1836). Mechanik der Menshlichen Gehwerkreuge, ~tie1 II, pp 161-202. Uber das Kniegelenk, Gottingen Winter, D. A. (1983). Moments of Force and Mechanical Power in Jogging. Journal of Biomechanics, Vol. 16, pp 91-97. Woltring, H. J. (1986). A Fortran Package for Generalized, Cross- Validatory Spline Smoothing and Differentiation. Adv. in Engineering Software, Vol. 8, pp 104-113. Woltring, H. J. (1990a). Definition and Calculus of Joint Angles, Axes, and Centers of Rotation from Noisy Position and Attitude Data. AIM Euroforum Meeting, Sevilla, Spain, December 13-15, p 1. Woltring, H. J. (1990b). RE: kinematics & kinetics. A discussion on Biomch-L, a biomechanics electronic bulletin board. The discussion is retrievable by sending the command: send Biomch-L log 9011 (on one line) to LISTSERV@HEARN.BITNET, a BITNET-electronic mail service. Woo, L., Freudenstein, F. (1970). Application of Line Geometry to Theoretical Kinematics and the Kinematic Analysis of Mechanical Systems. Journal of Mechanisms, Vol. 5, pp 417-461 Yamaguchi, G. T., Zajac, F. E. (1989). A Planar Model of the Knee Joint to Characterize the Knee Extensor Mechanism. Journal of Biomechanics, Vol. 22, pp 1-10. Yang, A.T. (1963), Application of Quaternion Algebra and Dual Numbers 117 to the Analysis of Spatial Mechanisms, Doctoral Dissertation, Columbia University, New York City, No. 64-2803 (University microfilm, Ann Arbor, MI). Yang, A. T. (1971). Inertial Force Analysis of Spatial Mechanisms. Journal of Engineering for Industry, Vol. 93, pp 27-33. Yang, A. T., Kirson, Y., Roth, B. (1975). On A Kinematic Curvature Theory For Ruled Surfaces. Proc. of the Fourth World Congress on the Theory of Machines and Mechanisms, Pub. by the ASME, pp 737-742. Yuan, M. S. C., Fruedenstein, F. (1971). Kinematic Analysis of Spatial Mechanisms by Means of Screw Coordinates. Part I - Screw Coordinates Journal of Engineering for Industry, vol. 93, pp 61-66.