- Email: [email protected]
A technique for obtaining spatial kinematic parameters of segments of biomechanical systems from cinematographic data
J Bwmechamcr Vol 13. pp. 535-547. c Per&mm Press Ltd. 1980. Prrnted m Great Bntam.
0021-9290/80m01-0535
$02.00/0
A TECHNIQUE FOR OBTAINING SPATIAL KINEMATIC PARAMETERS OF SEGMENTS OF BIOMECHANICAL SYSTEMS FROM CINEMATOGRAPHIC DATA* NORMANR. MILLER Assistant Professor of Mechanical Engineering and of Bioengineering, Department of Mechanical and Industrial Engineering, University of Illinois at Urbana-Champaign, Urbana, IL 61801, U.S.A.,
ROBERTSHAPIRO Assistant Professor of Physical Education, Department of Physical Education, Northern Illinois University, DeKalb, IL 60115, U.S.A.
and THOMASM. MCLAUGHLIN Assistant Professor of Physical Mucation, Department of Health, Physical Education and Recreation, Auburn University, Auburn AL 36830, U.S.A.
Abatract - Ibis paper describes a technique for establishing the three-dimensional kinematic parameters of a body using only the basic cinematographic
equipment available in most biomechanics
laboratories.
determined using only one camera, but accuracy is improved by utilizing data from two or more aimubaneous, syncbronixed film strips. In general, at least 3 points of known location must be available on each body. More points can be utitizcd to obtain a better solution in a Icaat squares aenae. Fewer than 3 poiatgaaeachbodyuencededwhantbebodiesarecaupladbyjomts~witbkno~cburctaiPtigrsina ma&me or approximately in biomecbanicnl systems. Tbe pammemrr wedad for biomcclrpaial anaIyais via matrix mctbods are obtainable dimctIy. The accuracy of the method is explored by means of a contr&d experiment u&g a tbree-dimcnaional mechanical linkage. In addition, the data was utilized to invcatigate the accuracy of several smoothing and differentiation techniques. Parameters can be
INTRODUCTION
et al. (1978) also used diode triades to define body
In recent years there has been increasing interest among workers in biomechanics in methods of obtaining the three-dimensional kinematic parameters of the various segments of the musculo-skeletal systems of humans or animals. A number of papers have described three-dimensional radiographic techniques (Chao and Morrey, 1978; Suh, 1974; Brown et al., 1976). These techniques are, of course, primarily suitable for determining static geometrical configurations. A number of authors have discussed various three-dimensional cinematographic technique-s (Bullock and Harley, 1972; Miller and Petake, 1973; Ayoub ef al., 1970; Erdman et al., 1975). Most of these authors have limited their discussion to techniques for finding the spatial coordinates of points from two photographs. Erdman et al. (1975) uses triades of light emitting diodes attached to body segments to determine the spatial parameters of the body via the method of the intersection of 3 spheres. Crowninshield
* Received 14 November 1978; in revisedjorm
7 May 1979.
positions. These techniques do not take into account the inevitable experimental error in photographically obtained point coordinates. This paper describes a rather general technique for obtaining the spatial location of segments of biomechanical systems. The emphasis is on simplicity of equipment. The technique also provides an estimate of the errors in positional parameters. THE RELATIONSHIPBETWEENCAMERASPACE AND OBJECTSPACE Every point in the portion of space within a camera’s view may be thought to have a unique mapping into the two-dimensional space represented by the photographic image. Unfortunately, the converse is not true, that is, a point on the image does not represent a unique point in space but rather a line in space. If multiple cameras are used to film the same scene, each camera will have its own unique mapping. If X, Y and 2 are the three-dimensional coordinates of a point in space and ui and vi are the corresponding photographic coordinates of the image of the point in the
NORMAN R. MILLER, ROEJEWSHAPIROand THOMASM. MCLAUGHLIN
536
photograph taken by the ith camera, then for the ith camera, there exists the following relationship:
L -cJ
Lz-zoJ
with C, = C/L,;C, = C/n,,; and the m,, are elements of the rotation matrix M. The systematic error terms, Aq and Ar, can be considered to be due to symmetrical and asymmetrical lens distortion and expressed as follows : Aq = q’(K,R’
-I- K,R4
+ KaR6
+ . ti.)
@a) + P,(R’ + 2q”) + 2P,q’r’ where C denotes camera principal distance; S scale factor; up and vi,, photo coordinates of the principle Ar = r’(K,R’ + K,R4 + KaR6 + . ..) point; [M], 3 x 3 rotation matrix ; and X0. Y,,,and (4b) ZO,three-dimensional coordinates of the camera per+ P,(R’ + 2r”) + 2P,q’r’, spective center. whereq’ = q - qO;r’ = r - to; R2 = q12 + r’l; K,,K2 This relationship is based on the basic theoretical and K s are coefficients of symmetrical lens distortion ; concept of photogrammetry; that is, that the photoand P, and P2 are coefficients of asymmetrical lens graph, being a perfect plane, is a central projection of distortion. the object space. In practice, the film coordinates are Equations (3) and (4) express the relationship not measured directly but instead are determined by between three-dimensional coordinates in space and the use of an optical comparator. If we assume that the the corresponding two-dimensional coordinates as coordinate systems of the photograph and the commeasured by an optical comparator. It is normal parator are parallel, we can make a transformation photogrammetric practice to determine the constants from the comparator coordinates to the photo coorin these equations by careful calibration of the cameras dinates and correct for linear film deformation, linear and comparators. In addition, the photogrammetric lens distortion and comparator errors as follows: camera must be constructed to image fiducial marks on the film to establish the image coordinate system. nt - u, - k(q + & - 40) (2) The specialized equipment required is extremely vi - vp = I& + Ar - ro), expensive. where J,, and I, are scale factors for each ax is ; q and r Abdel-Aziz and Karara (1971) regrouped the paraare the comparator coordinates corresponding to u meters of equations (3) to form the following and v; q. and r. are coordinates of the principal point relationships : in the comparator coordinate system ; and Aq and Ar L,X+L,Y+L,Z+L, _ are systematic errors in the comparator coordinates. q + Lpx + LlOY + LllZ + 1 Figure 1 illustrates the fundamental relationships of equations (1) and (2). Upon eliminating 1 from &X+L,Y+L,Z+Ls (5) equation (1) and substituting equations (2), we obtain r + Ar = b,x + L,oY + L,rZ + 1
4=
4+&-go
They solve for the 11 unknown parameters (L,, i = l-11) in quations (5) and as many of the parameters in equations (4) as desired by utilizing a two-pass m&Y Yo) m&Z-Zo) linear regression technique called the Direct Linear \ (3) Transformation (DLT) technique. The 11 parameters are referred to as the DLT parameters. This technique only provides estimates for the m21W-Xo)+ m22(Y-Yo) + m,,(Z-Z,) parameters of equations (3) and (4). but in practice it ‘mdX-X,) + MY-Y,) + m&Z-Z,)
mllW-X0) + m12(Y-Yo) + m,,(Z-Zo)
=-
cqmslW-Xo) + - +
r+Ar-r.
= _
c
Fig. 1. Relationship between spatial, photo and comparator coordinate systems.
Segments of biomcchanical systems
seems to perform in a very satisfactory fashion. The original application of the technique was for use with still cameras. Shapiro (1978) reports the use of the method with high speed cinematographic cameras to measure the acceleration of a ball in a free fall. Once the parameters in equations (3) or (5) and possibly equations (4) have been determined for two or more cameras, the equations can be used to determine the three-dimensional coordinates of any point which appears in at least two photographs. Each photograph provides data for one set of equations (5) in the three unknown coordinates, X, Y and 2. Thus, two photographs provide four equations in three unknowns which by simple algebraic rearrangement can be solved by least squares for the unknown coordinates.
537
in 4 x 4 (see Denavit and Hartenberg, notation as
1955) matrix
(6)
where 1
0
0
0
Tp,=
LOCATION OF A BODY IN SPACE WHICH CARRIES THREE POINTS WITH KNOWN COORDINATES
We have described a reasonably simple technique for determining the spatial coordinates of points falling in the field of view of two or more cameras. Let us assume that three points are attached to a body. Further, we assume that the position of the points are measured relative to a local coordinate system attached to the body. For simplicity, we assume that point A with local coordinates 0,O and 0 global coordinates (as measured photographically) X, Y, and 2, Jies at the origin of the local system. We assume point B (local coordinates xr,, 0 and 0, global coordinates X, Y, and 2,) lies on the local positive X axis and finally that point C(x, y, and 0 ; X, Y, and Z,) lies in the part of the local x, y plane where both x and y are greater than or equal to zero. Now if:5 G, are the three unit vectors in the global system andi’,j’ and G’are the unit vectors in the local system, i’s
ii =
G’ =i’ T,
I’
ZJG + (Y*- r,y + (z*-z,)y
(X, - x,); + (Y, - Y.h' + (Z, -
[(X*-X,,’
(X,-x,)i+ x
(Y,-Y,)j+
(Z,-ZJG
$1;’ x jq;
G’ XT’.
The position of the body can then be expressed by the origin of the local system (X,, Y., Z,) and the rotation matrix
J&I =
i.1 r.i’ 1 I 1
-_ ;.;’
G .T
cos (X, y)
cos (X, 2)
cos(Y,y)
cos(Y,z)
cos (Z, y)
cos (Z, z)
i.i’
T-G’ i.G’
G-i’
G.G’
.
The complete relationship between coordinates in the ground and local coordinate systems can be expressed
X0, Ye, Z. are the coordinates of the origin of the local coordinate system as measured in the global coordinate system. M,, can he expressed in terms of three Eulerian angles (see, for example, Greenwood, 1965). The set of angular rotations a about the Z axis, followed by fi about the new X axis and finahy y about the new Z axis will be used in this paper. The transformation matrix in terms of &se an&s is shown in the Appendix. Once the trrnrtanartion matrix has been calculated from the photographic data, one can extract the three Eukrian angles from the expanded version of the matrix.
A LEAST SQUARES TECHNIQUE FOR OBTAINING THE SIX SPATIAL PARAMETERS OF A RIGID BODY
The technique described above for obtaining the spatial coordinates of a body from photographic data suffers from two deficiencies. First, it is implicitly assumed that there is no error in the photographically obtained three-dimensional coordinates. Second, if more than three points are attached to the body, there is no way to take advantage of the additional information to improve the estimate of the spatial parameters. (One can, of course, calculate the-distance. between points on the body and in this manner eliminate points which are completely in error.) To overcome these problems, Jet us observe first of all that we are interested in the six spatial parameters of the body in question. Thus, let us combine equations (5) and (6), eliminating the global coordinates X, Y and Z. (We could just as well combine equations 3 and 6.) We write q+Aq=tj=NgfDg
and r+Ar=i=NrIDr
(7)
538
NORMANR. MILLER ROBERT SHAPIRO and THOMASM. MCLAUGHLIN
where the superscript 0 indicates evaluation using the current estimate of the six spatial parameters and AX0 = X0 - X$; AYo = Y. - Yg, etc. We can then write for the entire system of equations
where 1
Ng = L, + [0 Li Ls LJT,
’
of:,,
-
Y
ax0
z
hl% I
ar,
ay
af:,, a7
3ff12 _... ax0
1 Dg = 1 + (0 La Lro
af;,, af:,, -...-
’
kc
Y
AYo
Z
MO
and Nr and Dr are defined in a similar manner. We see that 4 and i are functions of the six spatial parameters, or
ami
6,
-j&f,,
Yo,
20,
a, A
YX
Aa
AB
63)
L
aft, _... ax0
aft,
L
AY
ay
where the subscript i indicates the local point on the
body and the subscript i indicates the camera The power ofthe technique now becomes apparent. ,Let us assume that we have only one camera (i = 1) and three points on the body (all visible to the camera). Then we can write quation (8) three times (i = 1,2 and 3), giving a total of six non-linear equations in the six unknown parameters, X0, Yo, Z,,, a, /I and y. These equations are independent and, given suitable initial estimates for the parameters can be solved iteratively by Newton’s method. Let us now examine the general case where we have more than three points on the body and more than one camera. We have one set of equations for each point visible in each camera’s image. Thus if all the points on the body are visible in all cameras, we will have 2(i x i) equations. (Notice that points need not be visible in all cameras to be used in the solution.) It is now possible to solve for a set of spatial parameters which minimizes the differences between the left- and right-hand sides of these equations in a least squares sense. Because of the highly non-linear nature of these equations, normal linear least squares does not apply. Many non-linear least squares algorithms exist, but it was decided to try the linearization technique described in Chapter 10 of Draper and Smith (1966) because of its close similarity to the iterative solution technique used in mechanism theory (see Uicker et al., 1964). Using a truncated Taylor series expansion, we write
+ aft* &#
+
da
art1 -A/?+-Ay ag
We now apply standard linear least squares techniques to this set of linearized equations; that is & z (gZ,)-iZ;F,,.
AN ESTIMATION OF THE ERROR IN THE SPATIAL PARAMEIEP ESTIMATE An
estimate of the variance-covariance parameters can be obtained from
$+Xo+$AY,,+...+dyAy, 0
0
(10)
We calculate X0 = AX: + X0, Y. = Y8 + A Yo, etc. and repeat the process using these new estimates of the spatial parameters in equation (9). The procedure continues until all elements of Ist are less than some arbitrary convergence criteria. Normally, the initial set of parameters used in quation (9) for the first frame in a set of fimstrips is obtained from the three coordinate technique described above. If this is not possible (say, only one camera was used), one can generally guess the parameters to an accuracy sufficiently close for convergence of the algorithm After the first set of parameters has been calculated, the starting values for the next frame are simply the final parameters from the last calculation. In practice, the linearization technique described above has been quite successful. Some problems with convergence might be experienced if the frame rate were slow relative to the velocity of the object. In such a situation, the initial values of the parameters provided by the preceeding calculation might not be sufficiently close to ensure convergence. Under these conditions, a more powerful algorithm such as described by Brown and Dennis (1972) would probably ensure convergence of quation (8).
aft1 av
and similarly, i&-f;,+
(9)
matrix of the
aft2 V(S) = (Z;Zo)-‘S2,
(11)
Segments of biomcchanical systems where (ZgTZO)-l is the same matrix used to evaluate equation (10) during the last iteration. S2 is an estimate of the variance of the residuals. It is calculated by summing the squares of the elements of vector Fb in equation (9) and dividing by the number of elements in F0 (that is, n) less six (the number of parameters estimated or p). For true linear regression, the joint confidence region for the six parameters would he an ellipsoid in six-dimensional space. This would be approximately true of this non-linear case (see pp. 64 and 273 of Draper and Smith, 1966). A rough estimate of the confidence band associated with each parameter can be calculated as
s.e. (&)J,
where Ci is the estimated value of one of the parameters ; t(u, 1 - (1/2)x,) is the “student-t” statistic with c==n - p degrees of freedom and at a risk level of a, Estimated s.e. ({J is the square root of the ith diagonal term of matrix (ZiZ,)- ’ in equation (11). This error estimate has a very small computational cost since the inverted matrix in equation (11) already exists. It can be used to assess the accuracy of the estimated parameters based solely on the data taken during an experiment Thus, it is most useful. In calculating this error estimate, the tacit assump tion is made that the DLT parameters are error free. Since, of course, this is not the case, these error estimates are changed in an unknown and complex manner. To obtain an idea of theefkt of the estimated DLT parameters, one can imagine 8 non-linear least squares approximation involving simultaneously the DLT parameters and the spatial parameters. The matrix Z, in equation (9) now includes partial derivatives with respect to each DLT parameter. In addition, vector F, will contain the residuals from fitting the DLT parameters for each camera. An equation analogous to equation (11) can be used to obtain an augmented variance-covariance matrix. Naturally, this is computationally expensive. Fortunately, as will be shown, the resulting error estimate is not significantly different from that obtained using equation (11). THE REDUCIION OF THE NUMBER OF PARAMETERS TO BE ESTIMATED BY THE USE OF GEOMETBK CONSTBAlNTS
Frequently in biomechanics, one is concerned with the measurement of the spatial parameters of a system of coupled bodies. It also frequently happens that one knows quite accurately the position of one body and wishes to find the position of a body attached to it by a joint. An example of this type occured in the study of tennis carried out by McLaughlin and Miller (1980). In that study, the position of a body consisting of the tennis racket and subject’s hand were relatively well known (having been measured by the methods de-
539
scribed above). It was desired to find the position of the forearm. Now one may consider the wrist to be approximately a Hooke or “Universal” joint having two degrees of freedom. Thus, to find the location of the forearm, one needs only estimate the two Hooke joint parameters. This can be done by a simple extension of the techniques described above. In this case, T*i = T,,@,(OT,,,
(12)
where T,,i is the 4 x 4 transformation matrix of the type described in equation (6) and used in equation (7). T,j is the known transformation from the fixed coordinate system to the joint in question. It may or may not have been found by the techniques described above. @At) is the 4 x 4 transformation for the joint type in question. For a Hooke joint, t represents two parameters consisting of the angles of the joint. T,! is a fixed known transformation relating the local coordinate system at the joint to the system in which the local coordinates on the body are known. Equation (12) can be substituted into equations (7) and a set of equations (8) can be written where the parameters in A,, andJ;1, sre only those of the joint. The Taylor series expansion will only involve the parameters for the joint. A{ in equation (9) wiii also only have as many elements as joint parameters. An iterative solution using equation (10) is used as above. The variance-covariance matrix of the estimpted joint vsri8bles can be found from an equation of the same form as equation (11). Finally, the confidence interval for the joint parameters can be calculated exactly as described above for the spetiai parameters. If one plans to use analysis techniques based upon 4 x 4 methods, such as the MAN program (see Williams, 1977), the estimates of the joint parameters will be sufficient. If, however, one wants the complete set of spatial parameters for the second body, one may use the estimated joint parameters in equation (12) and extract the six spatial parameters from the definition of the transformation matrix. One may estimate the variance-covariance matrix for the six spatial parameters by first constructing equation (9) using F, from the last stage of the evaluation of the joint parameters and using the six estimates of the spatial parameters to calculate Z,. One then calculates the approximate variance-covariance matrix for the six parameters from equation (11). This basic technique can be extended in an obvious way to chains of coupled rigid bodies with any number of joints.
EXPERIMENTALEVALUATION This portion of the paper describes an experiment utilizing a spatial mechanical linkage to evaluate the accuracy attainable using the methods described above.
540
NORMAN R. MILL~R,RO~~RTSHAPIRO andTIR)m M. MCLAUGHLIN
The experimentul linknge While the techniques described in the first part of this paper were under development, a spatial linkage was constructed. This linkage is shown in Fig. 2. A dc. motor with a speed controller drives the input crank which carries a Hooke “universal” joint. This joint is attached to the coupler link which is attached to the output link with a ball and socket joint. The coupler link has true spatial motion. The coupler was, therefore, chosen as the experimental body, The distance from the center of the ball to the amter of the Hooke joint is 25.4ctn This is approximately the length of the human forearm. While many investigators have used light emitting diodes (LED) as photographic markers, it was decided to use small balls, since this type of marker must be used for activities requiring bright illumination. The balls arc 1.8cm in diameter, which, it turns out, is larger than needed. Some error is introduced by assuming that the center of the ball’s photographic image lies at the same point as the image of a geometrical point at the ball’s center would lie. Calculations show that this error is smaller than the rcaolution error of the Vangard Motion Analyzer (VMA) used to quantify the photograph coordinates. The coupler carries 8 photographic markers on metal rods of various lengths. The 4 longest stand 4.8 cm from the centerline of the coupler. The 2 shorter markers at the Hooke joint end arc 3 cm long and the similar pair at the ball joint end are 4.3 cm long. The shorter markers approximately represent markers affixed directly to the surface of the forearm, while the longer markers represent points attached to extenders. A number of other markers arc attached to other parts of the machine. These markers serve as reference points for the computation of the DLT parameters.
1 sec. A frame rate of 100frames per set was more than sufficient for this experiment. Higher speed motions encountered in actual biomechanics experiments will require higher frame rates. The cameras were placed in accordance with the recommendations of Marzan (1976); that is, the cameras were separated by onethird of the camera to subject distance with the photographic axes very slightly convergent. In this case, the camera-sub&t distance was about 2.3 m and the inter-camera distance, 0.8 m. The important point here is that this is not critical. The first step of the experiment was to photograph the mechanism while in a stationary reference position. This reference position (which was set by a mechanical stop) was such that the Hooke joint was at its topmost position. The cameras were not moved after photographing the mechanism in this position. The frames shot at this reference position were later used to establish the camera space to object space DLT parameters. The mechanism was then put into motion and was filmed simultaneously by the two cameras during severalcycles of motion. A potentiometer was attached to the input crank shaft and simultaneously fired, by means of a voltage comparator circuit., the light emitting diode event markers in each camera, Unfortunately, as will be discussed later, this only served to approximately synchronize the cameras to each other and the linkage. After processing, point data was collected from each f&n&p usinga VMA.The data was input to and stored on files on the campus CYBER 175 computer. The frame rates of the two cameras were similar, but not identical.One of the filmstripswas chosen as a reference and the second filmstrip was artificially synchronized to it by cubic spline interpolation. This process is an obvious source of error.
Experimental techniques
Data were gathered by filming the linkage, using two Locam, pin registered high-speed cameras loaded with Kodak RAR 2498film on an Estar-AH base. Both cameras were set to film at 100 frames per sec. The speed of the motor driving the linkage was adjusted so that the linkage completed one cycle in approximately
Hooke Joint
r4-
-&put Fig. 2(a).
Lir 4
Schematicdiagramof experimentallinkage.
Data reduction The position of all the photographic markers in the reference positi& was calculated using the standard kinematic techniques of Uicker et al. (1964). The DLT parameters were calculated using various combinations of marker points. The routine reported in Marzan (1975, 1976) was used for these computations. The DLT parameters from both cameras were used to back calculate the positions of the marker points allowing errors to be calculated. In addition, as will be described below, several sets of DLT parameters were used to calculate the position of the coupler as a function of time. From these results, it was concluded that the use of the 8 points on the mechanism’s coupler, along with the 2 points on the arm carrying the ball joint, produced the best overall result. Using 10 points to calculate the DLT parameters is a bare minimum. More marker points should have been provided in the region near the coupler. This points out the advisability of providing markers throughout the region in which a body is to move.
Fig. 2(b). Experimental spatial linkage in reference position.
541
segments of biomechanical systems
The root mean square errors in estimating the position of the marker points using the 10 point DLT parameters are, respectively, for X, Y and Z, 0.176 cm, 0.428 cm and 0.468 cm. Clearly, these errors are’large. This is rather surprising in view of the excellent results obtained by Professor Karrera and his students using 35 mm still cameras. (Their results are typically two orders of magnitude more accurate.) Any technique of this type is critically dependent on the equipment used. Problems with cameras, film, comparitor or the skill of the personnel collecting the data will adversely influence the results. It is unclear where the weak link in this chain lies. Obviously, the poor precision afforded by the estimated DLT parameters will adversely influence the balance of the experimental results. We shall see, however, that we will be able to estimate the position of the body to a higher accuracy. The DLT algorithm was used to estimate more than the basic 11 parameters. This procedure was unsuccessful. The program provided by Prof. Karrera does not seem to work for more than 11 DLT parameters. A program was written which used Karrera’s 11 DLT parameters to estimate starting values for the 10 parameters implicit in equations (3). The program then used the modified algorithm of Levenberg-Marquardt (Brown and Dennis, 1972) to obtain true least squares estimates of the 10 parameters implicit in equation (3) and then to estimate the various parameters of equation (4). It was found that this procedure was expensive computationally and only provided a limited reduction in the error of the estimate of the position of the marker points as calculated using the photographic data from the two cameras. As a result, it was concluded that inaccuracy was due to errors in the data taken at the reference position and not due to deficiencies in the DLT algorithm. The 11 DLT parameters originally estimated for each camera were utilized throughout the balance of the computations. A series of Fortran programs were written for the CYBER 175 computer system to perform the computations described in the first part of this paper. The results of these computations follow.
543
ball joint and the Y axis is directed roughly downward when the hook joint is in its topmost position. Figure 3(a) shows the location of the coupler center as a function of time, as computed from kinematic theory. The mean coordinates of the coupler position have been extracted to produce plots centered about zero. Figure 3(b) shows the corresponding computed Eulerian angles. Again, the mean angles have been extracted from the data prior to plotting. These plots serve to indicate the range of motion of the linkage. The techniques described in the first part of this paper were next used to estimate the spatial parameters from the cinematographic data. The errors between these estimated values and the calculated values shown in Fig. 3 were evaluated. Figure 4 shows a typical plot of error vs frame. In addition, the 95% confidence interval about the estimated value of the parameter was computed for each frame for each parameter, as discussed above. The theoretical value of the parameter was subtracted from the upper and lower bound on this interval, and the results plotted in the figure. If we had no knowledge of the true value of the parameter, we would hope that the 95% confidence interval about the parameter estimate would embrace the true parameter value. In our plot, we would hope to see these bounds embrace zero error. We see that in general this condition is met. Figure 5 shows the same error plot as Fig. 4 with the 95% confidence region calculated using the more complete estimate of the variance described above. As can be seen, this confidence region is not notably different from the one obtained using the simpler method. Since this technique more than doubles computation time, the first estimate of the confidence region is to be preferred. The data plotted in both Figs. 4 and 5 were computed using all 8 markers carried by the coupler and data from both cameras. Several additional runs were made in which data from : (1) 8 markers but only the left-hand camera were Used;
(2) 8 markers but only the right-hand camera were Used;
RESULTS
The theoretically correct transformation from the ground coordinate system to a coordinate system moving with the coupler of the linkage was calculated to correspond to each set of photographic data throughout the cycle. Calculations were performed as described in Uicker et al. (1964). The ground coordinate system has its center at the point where the input crank shaft projects from its bearing. The X axis is directed downward, the Z along the input shaft toward the middle of the apparatus and the Y axis in the proper direction as defined by the right-hand convention. The coordinate system attached to the coupler has its origin at the geometric center of the link, the local X axis runs from this center toward the
(3) only the shortest 4 markers on the coupler and from both cameras were used ; (4) all eight markers and both cameras were used, but the position of the input crank was calculated so only the two Hooke joint parameters were estimated. The true root mean square error over one cycle for each parameter under each condition is indicated in Table 1. As can be seen, the error for the left-hand camera seems to be significantly larger than the error for the right-hand camera. Notice that the results with 4 markers on the coupler are not significantly worse than those for 8 markers. This is surprising and probably reflects errors in synchronization between the two cameras and the mechanism. Finally, the results from the estimation of only the two Hooke joint
NORMAN R. MILLER, Roem
544
Frame
SHAPIROand THOMAS M. MCLAUGHLIN
Number
Fig. 3(a). Position of coupler center as computed from kinematic theory.
parameters are extremely good. The error in x is zero, since specifying the input crank position fixes this parameter. This section will be concluded with a few words concerning the synchronization problem alluded to above. Synchronization was achieved by firing the light emitting diode event markers in both cameras using a potentiometer on the input crank as a trigger. Unfortunately, due to the pin registered drive of the cameras and the resultant intermittent motion of the film, there is no way of knowing to within less than one frame, the synchronization of the two cameras to each other and the mechanism. The Locam Camera can be obtained with an option to output a signal indicating shutter opening. With this information, one could largely overcome the synchronization problem. Best of ah, it is now possible to obtain phase-locked shutter systems. This would allow nearly perfect synchronization. COMPUTATlON OF SECOND TIME DERIVATWES OF THE PARAMETERS
pfl.564rad
,,...”
-41 0
_.’
6= 4.697rad I
20
40
60
Correspandii Frame
80
loo
I20
Numtw
Fig. 3(b). Angular position of coupler as computed from kinematic theory.
Table
Usually the object of measuring the spatial parameters of the segments of the musculo-skeletal system is to 6nd their tirst- and second-time derivatives and use them in the calculation of kinetic energy or in the computation of forces and moments using the Euler equations. Many smoothing and differentiation techniques are in common use. Experimental position data were available in this study and the theoretically correct second time derivatives of the spatial parameters could be calculated. This being the case, it was decided to evaluate the accuracy of several smoothing and differentiation schemes. Note that error estimates for each parameter at each frame are available, thus weighted least squares techniques can be used to fit smoothing functions. Figures 6(a) and (b) show the calculated second derivative of each parameter. This computation involves a rather direct extension of the techniques of
1. Root mean square error over one cycle for various types of data reduction
dondition
RMS error over cycle a (rad) z, (cm)
B (rad)
Y(rad)
0.019
0.014
0.010
0.692
0.023
0.017
0.014
0.295
0.442
0.018
0.014
0.012
0.080
0.117
0.160
0.019
0.011
0.008
0.066
0.071
0.044
0.0
0.009
0.004
x0 (cm)
Y0 (cm)
0.073
0.096
0.149
left camera only
0.106
0.947
8 markcrsright camera only
0.046
4 markersboth cameras 8 markersestimated only two Hook joint parameters
8 markers on coupkrboth cameras 8 markers-
545
Segments of biomechanical systems
_ *r
30
_ ‘.
4-
.::
‘. ..
. ‘.
..
;.
.’
..‘..
.., ..IJpper
of 95
‘.
,_
.;
I
r
Bound %
.‘-’ Conhdence
.,
Lower Bound
-6 0
20
40
80
60 Frmne
, 120
loo
-30
1
0
Number
I 20
40
Cornrpondmg
Fig. 4. Error in YOparameter as a function of frame number showing the estimated error band based on the 95% confidence interval about the estimated parameter value.
Upper Bound
40
80
60
Frame
of 95% Confidence
‘.
._..’
IO0
120
Number
Fig. 5. Error in YOparameter as a limction of frame number showing the estimated error band based on the 95% confidence interval about the estimated parameter value (error band based on more exact estimate). 60
I
I
60
80
Frame
I
I
100
120
Number
Fig. 6(b). Calculated angular acceleration of coupler.
Denavit et al. (1965). Basically, one calculates the firstand second-time derivatives of the global coordinates of three points on the body and writes the first total time derivatives of equation (6). He then solves for the first time derivatives of the spatial parameters. Differentiating quation (6) again, one obtains quations for the second time derivatives. Various techniques for smoothing and numerical diffematiation were applied to the #mtographic po sition data. They include the following: (1) simple second-order finite differences; (2) digital filtering wing a Butterworth filter followed by finite di&rences (Winter et aL, 1974); (3) fitting of polynomials by weighted linear least squares, followed by differentiation of the polynomials ; (4) fitting of a Fourier series via weighted linear least squares followed by differentiation of the resulting Fourier series (Cappozzo et al., 1975); (5) fitting of a cubic spline using the weighting technique described in Reinsch (1967) and McLaughlin et ol. (1977). The weights for the spline estimate were obtained as follows : If (s.e.), is the standard error of the parameter estimate for a given’ parameter and frame, D, is the difference between the spline value and data value at a frame and NF is the number of frames in the data set. Then 5
[D&e.),]’
I NF.
I=1
-60
I 0
20
40 Correspondlng
60 Frame
80
too
120.
Number
Fig. 6(a). Calculated acceleration of coupler center.
Table 2 lists the root mean square error over a cycle for each parameter’s second derivative using each smoothing technique. Notice that the digital filter, the Fourier series and the cubic spline all work well. The Fourier series is not a good technique to choose for smoothing a. (Indeed, the best technique would be to fit a straight line.) It would appear that different techniques are applicable
546
NORMAN
R. MILL& Rotwtr SHAPIRO and Txoh~ls M. MCLAUGHLIN
Table 2. Root mean square error in second parameter time derivatives over one cycle for various smoothing and differentiation techniques RMS error over cycle 20 (radiec’) ~cm/sec’)
x0 (cmkc’ ) Second order tinite ditkrences
316
Digital filmring (2 Hz break for linear data; 3 Hz for angular data) Polynomials (3rd dcgrtc)
1258
3.50
6.09
26.7
28.2
1250
32
7.68
.. @ad/se? )
25
1.29
16.6
B (rad/sec’)
0.101
35
1.08
0.96
13.8
11.6
Fourier series (X-2nd harmonic; Y, Z, a-fundamental ; j?,y-5th harmonic )
4.72
1.91
4.07
N.A.
2.90
1.40
Cubic spline
3.66
4.72
5.41
0.95
1.02
1.11
to different parameters This simply emphasixes that there is no one “best” smoothing technique and that plots and common sense are indispensible if one wishes to obtain reasonable results. Figures 7(a) and (b) show the error resulting from use of the cubic spiine as a function of frame number. These error plots should be compared with the calculated accekrations in Fig 6. Notice that the magnitude of the maximum error in X0 is about 8cm/sec2, while the piak magnitude of x,, is about 55 cm/sec2 ; thus, the ratio of peak error to peak magnitude is g/55 - 0.15. As one can see, while the accuracy is not exceptional, one can obtain an order of magnitude evaluation of acceleration using this technique. Further, these results can probably be taken as an indication of the minimum accuracy obtainable by these methods.
2-
-2 . *'
3O
20
40 60 80 Frame Number
NM
I I20
Fig. 7(b). Angular acceleration error - cubic splint
CONCLUSIONS
-I”
0
20
40 60 80 Frame Number
100
I20
Fig. 7(a). Linear acceleration error - cubic spline.
A technique for photographically obtaining the spatial kinematic parameters of a rigid body has been described. A minimum of one camera and three points on the body are required. Accuracy can be improved by the use of multiple cameras. The method provides estimates of the errors on the calculated parameters. The method extends in a natural way to estimating the positions of coupled chains of rigid bodies. In this case, a reduced set of parameters needs to be estimated. A controlled experiment was conducted to study the accuracy of the method. Very good positional accuracy was achieved and a reasonable level of acceleration accuracy was obtained by the use of several smoothing
and differentiation
techniques.
547
Segments of biomechanicai systems
Draper, N. and Smith, H. (1966) Applied Regression Analysis. Wiley, New York. Erdman, A. G., Forstrom, R., Dorman, F., Ah&en, F:, Beer, R., Isaacson, R. J. and Speidel, M. T. (1975) Experimental motion analysis via a stereo-photography technique: ap plication to the human jaw and mechanical linkages, ASME Paper No. 75-DET-58. Greenwood, D. T. (1965) Principles OfDynamics, pp. 333-335. Prentice-Hail, Englewood Cliffs. Marxan, T. and Karrara, H. M. (1975) A computer program for direct linear transformation of the colinarity condition and some applications of it. Symp. on Close Range Photogrammetric Systems. American Society of Photogrammetry, Fails Church. Marxan. T. (1976) Rational design for close-range photogrammetry. Ph.D. Dissertation. University of Illinois at Urbana-Champaign. McLaughlin, T. M., Diiiman, C. J. and Lardner, T. J. (1977) Biomechanicai analysis with cubic spiine functions. Res. Q. 48, 570-582. McLaughlin, T. M. and Miller, N. R. (1980) Techniques for the evaluation of loads on the forearm prior to impact in tennis strokes. J. Mech. Design Trans. ASME, in press. Miller, D. 1. and Petak, K. (1973) Three-dimensional cinematography. Kinesiology III. American Association for Health, Physical Education and Recreation, Washington. Reinsch, D. H. (1967) Smoothing by spiine functions. h’um. Math. 10, 177-183. Shapiro, R. (1978) The direct linear transformation method for threedimensional cinematography. Res. Q. 49, 197-205. Suh, C. H. (1974) The fundamentals of computer aided X-ray analysis of the spine. J. Biomechanics 7, 161-169. Uicker, Denavit, J. and Harm&erg, R. S. (1964) An iterative method for the displacement analysis of spatial mechanisms. J. appl. Mech. 31. Trans. ASME, Se&s E, g6, pp. 309-3 14. Williams, R. and Seireg, A. (1977) Intaractive computer modeling of the muscuio-skeletal system. IEEE Trans. biomed. Engr. BME-24,213-219. Winter, D. A., Sidwail, H. G. and Hobson, D. A. (1974) Measurement and reduction of noise in kinematics of locomotion. J. Biomechanics 17, 157-159.
Acknowledgements - The authors would like to thank Professor Karrera of the Department of Civil Engineering, University of Illinois at Urban-Champaign, for supplying us with a copy of the program used to estimate the DLT parameters (Marxan and Karrera, 1975).
REFERENCES Abdei-A&, Y. I. and Karara, I-l. M. (1971) Direct linear transformation from comparator coordinates into object space coordinates in close-range photogrammetry. ASP Symp. on Close-Range Photogrammetry. American Society of Photogrammetry, Fails Church. Ayoub, M. A., Ayoub, M. M. and Ramsey, J. D. (1970) A steriometric system for measuring human motion. Hum. Factors 12, 523-535.
Brown R. H., Burstein, A. H., Nash, C. L. and Shock, C. C. (1976) Spinal analysis using a thradimensionai radiographic technique. J. Biomechanics 9, 355-365. Brown, K. M. and Dennis, J. E. (1972) Derivative free anaiogues of the Lever&erg-Marquardt and Gauss aigorithms for nonlinear least squares approximations. Num. Math. 18, 289-297. Bullock, M. and Harley, I. A. (1972) The measurement of three-dimensional body movements by the use of photogrammetry. Ergonomics 15, 309-322. Cappozxo, A., Leo, T. and Pedatti, A. (1975) A general computing method for the analysis of human locomotion. J. Biomechanics 18,307-320. Chao, E. Y. and Morrey, B. F. (1978) Three-dimensional rotation of the elbow. J. Biomechanics 11, 57-73. Crowninsbield, R. D, Johnston, R. C., Andrews, J. G. and Brand, R. A. (1978) A biomechanicai investigation of the human hip. J. Biunechanics 11, 75-85. Denavit, J. and Harms&erg, R. S. (1955) A kinematicnotation for lower pair machanisms based on matrices. J. appl. Mech. 22. Trans. ASME T7,215-221. Denavit, J., Harm&erg, R. S., Rti and Uicker. (1965) Velocity, acceleration, and static force analyses of spatial linkages. J. appl. Mech. 32. Trans. ASME. Series E, lt7, 903-Y 10.
APPENDIX: The Tran$ormation
Matrix in Terms of Eulerian Angles
1 cos u,, cos ya,
-
cosaIl m
Y#I
sin a,, sin &, -sin apl cos /I,, sin yI, sm a#, cos Y,~
-sin aII cos &, cos ye, - sm aorsm ysr - cos a,, sin &,
T,, = + cos e,t cos B,r sin ysl sin /J,, sin ys,
+ cos a,1 cos 8,, cos Y#I sin B,t cos Y,I
cos B,,
0021-9290/80m01-0535
$02.00/0
A TECHNIQUE FOR OBTAINING SPATIAL KINEMATIC PARAMETERS OF SEGMENTS OF BIOMECHANICAL SYSTEMS FROM CINEMATOGRAPHIC DATA* NORMANR. MILLER Assistant Professor of Mechanical Engineering and of Bioengineering, Department of Mechanical and Industrial Engineering, University of Illinois at Urbana-Champaign, Urbana, IL 61801, U.S.A.,
ROBERTSHAPIRO Assistant Professor of Physical Education, Department of Physical Education, Northern Illinois University, DeKalb, IL 60115, U.S.A.
and THOMASM. MCLAUGHLIN Assistant Professor of Physical Mucation, Department of Health, Physical Education and Recreation, Auburn University, Auburn AL 36830, U.S.A.
Abatract - Ibis paper describes a technique for establishing the three-dimensional kinematic parameters of a body using only the basic cinematographic
equipment available in most biomechanics
laboratories.
determined using only one camera, but accuracy is improved by utilizing data from two or more aimubaneous, syncbronixed film strips. In general, at least 3 points of known location must be available on each body. More points can be utitizcd to obtain a better solution in a Icaat squares aenae. Fewer than 3 poiatgaaeachbodyuencededwhantbebodiesarecaupladbyjomts~witbkno~cburctaiPtigrsina ma&me or approximately in biomecbanicnl systems. Tbe pammemrr wedad for biomcclrpaial anaIyais via matrix mctbods are obtainable dimctIy. The accuracy of the method is explored by means of a contr&d experiment u&g a tbree-dimcnaional mechanical linkage. In addition, the data was utilized to invcatigate the accuracy of several smoothing and differentiation techniques. Parameters can be
INTRODUCTION
et al. (1978) also used diode triades to define body
In recent years there has been increasing interest among workers in biomechanics in methods of obtaining the three-dimensional kinematic parameters of the various segments of the musculo-skeletal systems of humans or animals. A number of papers have described three-dimensional radiographic techniques (Chao and Morrey, 1978; Suh, 1974; Brown et al., 1976). These techniques are, of course, primarily suitable for determining static geometrical configurations. A number of authors have discussed various three-dimensional cinematographic technique-s (Bullock and Harley, 1972; Miller and Petake, 1973; Ayoub ef al., 1970; Erdman et al., 1975). Most of these authors have limited their discussion to techniques for finding the spatial coordinates of points from two photographs. Erdman et al. (1975) uses triades of light emitting diodes attached to body segments to determine the spatial parameters of the body via the method of the intersection of 3 spheres. Crowninshield
* Received 14 November 1978; in revisedjorm
7 May 1979.
positions. These techniques do not take into account the inevitable experimental error in photographically obtained point coordinates. This paper describes a rather general technique for obtaining the spatial location of segments of biomechanical systems. The emphasis is on simplicity of equipment. The technique also provides an estimate of the errors in positional parameters. THE RELATIONSHIPBETWEENCAMERASPACE AND OBJECTSPACE Every point in the portion of space within a camera’s view may be thought to have a unique mapping into the two-dimensional space represented by the photographic image. Unfortunately, the converse is not true, that is, a point on the image does not represent a unique point in space but rather a line in space. If multiple cameras are used to film the same scene, each camera will have its own unique mapping. If X, Y and 2 are the three-dimensional coordinates of a point in space and ui and vi are the corresponding photographic coordinates of the image of the point in the
NORMAN R. MILLER, ROEJEWSHAPIROand THOMASM. MCLAUGHLIN
536
photograph taken by the ith camera, then for the ith camera, there exists the following relationship:
L -cJ
Lz-zoJ
with C, = C/L,;C, = C/n,,; and the m,, are elements of the rotation matrix M. The systematic error terms, Aq and Ar, can be considered to be due to symmetrical and asymmetrical lens distortion and expressed as follows : Aq = q’(K,R’
-I- K,R4
+ KaR6
+ . ti.)
@a) + P,(R’ + 2q”) + 2P,q’r’ where C denotes camera principal distance; S scale factor; up and vi,, photo coordinates of the principle Ar = r’(K,R’ + K,R4 + KaR6 + . ..) point; [M], 3 x 3 rotation matrix ; and X0. Y,,,and (4b) ZO,three-dimensional coordinates of the camera per+ P,(R’ + 2r”) + 2P,q’r’, spective center. whereq’ = q - qO;r’ = r - to; R2 = q12 + r’l; K,,K2 This relationship is based on the basic theoretical and K s are coefficients of symmetrical lens distortion ; concept of photogrammetry; that is, that the photoand P, and P2 are coefficients of asymmetrical lens graph, being a perfect plane, is a central projection of distortion. the object space. In practice, the film coordinates are Equations (3) and (4) express the relationship not measured directly but instead are determined by between three-dimensional coordinates in space and the use of an optical comparator. If we assume that the the corresponding two-dimensional coordinates as coordinate systems of the photograph and the commeasured by an optical comparator. It is normal parator are parallel, we can make a transformation photogrammetric practice to determine the constants from the comparator coordinates to the photo coorin these equations by careful calibration of the cameras dinates and correct for linear film deformation, linear and comparators. In addition, the photogrammetric lens distortion and comparator errors as follows: camera must be constructed to image fiducial marks on the film to establish the image coordinate system. nt - u, - k(q + & - 40) (2) The specialized equipment required is extremely vi - vp = I& + Ar - ro), expensive. where J,, and I, are scale factors for each ax is ; q and r Abdel-Aziz and Karara (1971) regrouped the paraare the comparator coordinates corresponding to u meters of equations (3) to form the following and v; q. and r. are coordinates of the principal point relationships : in the comparator coordinate system ; and Aq and Ar L,X+L,Y+L,Z+L, _ are systematic errors in the comparator coordinates. q + Lpx + LlOY + LllZ + 1 Figure 1 illustrates the fundamental relationships of equations (1) and (2). Upon eliminating 1 from &X+L,Y+L,Z+Ls (5) equation (1) and substituting equations (2), we obtain r + Ar = b,x + L,oY + L,rZ + 1
4=
4+&-go
They solve for the 11 unknown parameters (L,, i = l-11) in quations (5) and as many of the parameters in equations (4) as desired by utilizing a two-pass m&Y Yo) m&Z-Zo) linear regression technique called the Direct Linear \ (3) Transformation (DLT) technique. The 11 parameters are referred to as the DLT parameters. This technique only provides estimates for the m21W-Xo)+ m22(Y-Yo) + m,,(Z-Z,) parameters of equations (3) and (4). but in practice it ‘mdX-X,) + MY-Y,) + m&Z-Z,)
mllW-X0) + m12(Y-Yo) + m,,(Z-Zo)
=-
cqmslW-Xo) + - +
r+Ar-r.
= _
c
Fig. 1. Relationship between spatial, photo and comparator coordinate systems.
Segments of biomcchanical systems
seems to perform in a very satisfactory fashion. The original application of the technique was for use with still cameras. Shapiro (1978) reports the use of the method with high speed cinematographic cameras to measure the acceleration of a ball in a free fall. Once the parameters in equations (3) or (5) and possibly equations (4) have been determined for two or more cameras, the equations can be used to determine the three-dimensional coordinates of any point which appears in at least two photographs. Each photograph provides data for one set of equations (5) in the three unknown coordinates, X, Y and 2. Thus, two photographs provide four equations in three unknowns which by simple algebraic rearrangement can be solved by least squares for the unknown coordinates.
537
in 4 x 4 (see Denavit and Hartenberg, notation as
1955) matrix
(6)
where 1
0
0
0
Tp,=
LOCATION OF A BODY IN SPACE WHICH CARRIES THREE POINTS WITH KNOWN COORDINATES
We have described a reasonably simple technique for determining the spatial coordinates of points falling in the field of view of two or more cameras. Let us assume that three points are attached to a body. Further, we assume that the position of the points are measured relative to a local coordinate system attached to the body. For simplicity, we assume that point A with local coordinates 0,O and 0 global coordinates (as measured photographically) X, Y, and 2, Jies at the origin of the local system. We assume point B (local coordinates xr,, 0 and 0, global coordinates X, Y, and 2,) lies on the local positive X axis and finally that point C(x, y, and 0 ; X, Y, and Z,) lies in the part of the local x, y plane where both x and y are greater than or equal to zero. Now if:5 G, are the three unit vectors in the global system andi’,j’ and G’are the unit vectors in the local system, i’s
ii =
G’ =i’ T,
I’
ZJG + (Y*- r,y + (z*-z,)y
(X, - x,); + (Y, - Y.h' + (Z, -
[(X*-X,,’
(X,-x,)i+ x
(Y,-Y,)j+
(Z,-ZJG
$1;’ x jq;
G’ XT’.
The position of the body can then be expressed by the origin of the local system (X,, Y., Z,) and the rotation matrix
J&I =
i.1 r.i’ 1 I 1
-_ ;.;’
G .T
cos (X, y)
cos (X, 2)
cos(Y,y)
cos(Y,z)
cos (Z, y)
cos (Z, z)
i.i’
T-G’ i.G’
G-i’
G.G’
.
The complete relationship between coordinates in the ground and local coordinate systems can be expressed
X0, Ye, Z. are the coordinates of the origin of the local coordinate system as measured in the global coordinate system. M,, can he expressed in terms of three Eulerian angles (see, for example, Greenwood, 1965). The set of angular rotations a about the Z axis, followed by fi about the new X axis and finahy y about the new Z axis will be used in this paper. The transformation matrix in terms of &se an&s is shown in the Appendix. Once the trrnrtanartion matrix has been calculated from the photographic data, one can extract the three Eukrian angles from the expanded version of the matrix.
A LEAST SQUARES TECHNIQUE FOR OBTAINING THE SIX SPATIAL PARAMETERS OF A RIGID BODY
The technique described above for obtaining the spatial coordinates of a body from photographic data suffers from two deficiencies. First, it is implicitly assumed that there is no error in the photographically obtained three-dimensional coordinates. Second, if more than three points are attached to the body, there is no way to take advantage of the additional information to improve the estimate of the spatial parameters. (One can, of course, calculate the-distance. between points on the body and in this manner eliminate points which are completely in error.) To overcome these problems, Jet us observe first of all that we are interested in the six spatial parameters of the body in question. Thus, let us combine equations (5) and (6), eliminating the global coordinates X, Y and Z. (We could just as well combine equations 3 and 6.) We write q+Aq=tj=NgfDg
and r+Ar=i=NrIDr
(7)
538
NORMANR. MILLER ROBERT SHAPIRO and THOMASM. MCLAUGHLIN
where the superscript 0 indicates evaluation using the current estimate of the six spatial parameters and AX0 = X0 - X$; AYo = Y. - Yg, etc. We can then write for the entire system of equations
where 1
Ng = L, + [0 Li Ls LJT,
’
of:,,
-
Y
ax0
z
hl% I
ar,
ay
af:,, a7
3ff12 _... ax0
1 Dg = 1 + (0 La Lro
af;,, af:,, -...-
’
kc
Y
AYo
Z
MO
and Nr and Dr are defined in a similar manner. We see that 4 and i are functions of the six spatial parameters, or
ami
6,
-j&f,,
Yo,
20,
a, A
YX
Aa
AB
63)
L
aft, _... ax0
aft,
L
AY
ay
where the subscript i indicates the local point on the
body and the subscript i indicates the camera The power ofthe technique now becomes apparent. ,Let us assume that we have only one camera (i = 1) and three points on the body (all visible to the camera). Then we can write quation (8) three times (i = 1,2 and 3), giving a total of six non-linear equations in the six unknown parameters, X0, Yo, Z,,, a, /I and y. These equations are independent and, given suitable initial estimates for the parameters can be solved iteratively by Newton’s method. Let us now examine the general case where we have more than three points on the body and more than one camera. We have one set of equations for each point visible in each camera’s image. Thus if all the points on the body are visible in all cameras, we will have 2(i x i) equations. (Notice that points need not be visible in all cameras to be used in the solution.) It is now possible to solve for a set of spatial parameters which minimizes the differences between the left- and right-hand sides of these equations in a least squares sense. Because of the highly non-linear nature of these equations, normal linear least squares does not apply. Many non-linear least squares algorithms exist, but it was decided to try the linearization technique described in Chapter 10 of Draper and Smith (1966) because of its close similarity to the iterative solution technique used in mechanism theory (see Uicker et al., 1964). Using a truncated Taylor series expansion, we write
+ aft* &#
+
da
art1 -A/?+-Ay ag
We now apply standard linear least squares techniques to this set of linearized equations; that is & z (gZ,)-iZ;F,,.
AN ESTIMATION OF THE ERROR IN THE SPATIAL PARAMEIEP ESTIMATE An
estimate of the variance-covariance parameters can be obtained from
$+Xo+$AY,,+...+dyAy, 0
0
(10)
We calculate X0 = AX: + X0, Y. = Y8 + A Yo, etc. and repeat the process using these new estimates of the spatial parameters in equation (9). The procedure continues until all elements of Ist are less than some arbitrary convergence criteria. Normally, the initial set of parameters used in quation (9) for the first frame in a set of fimstrips is obtained from the three coordinate technique described above. If this is not possible (say, only one camera was used), one can generally guess the parameters to an accuracy sufficiently close for convergence of the algorithm After the first set of parameters has been calculated, the starting values for the next frame are simply the final parameters from the last calculation. In practice, the linearization technique described above has been quite successful. Some problems with convergence might be experienced if the frame rate were slow relative to the velocity of the object. In such a situation, the initial values of the parameters provided by the preceeding calculation might not be sufficiently close to ensure convergence. Under these conditions, a more powerful algorithm such as described by Brown and Dennis (1972) would probably ensure convergence of quation (8).
aft1 av
and similarly, i&-f;,+
(9)
matrix of the
aft2 V(S) = (Z;Zo)-‘S2,
(11)
Segments of biomcchanical systems where (ZgTZO)-l is the same matrix used to evaluate equation (10) during the last iteration. S2 is an estimate of the variance of the residuals. It is calculated by summing the squares of the elements of vector Fb in equation (9) and dividing by the number of elements in F0 (that is, n) less six (the number of parameters estimated or p). For true linear regression, the joint confidence region for the six parameters would he an ellipsoid in six-dimensional space. This would be approximately true of this non-linear case (see pp. 64 and 273 of Draper and Smith, 1966). A rough estimate of the confidence band associated with each parameter can be calculated as
s.e. (&)J,
where Ci is the estimated value of one of the parameters ; t(u, 1 - (1/2)x,) is the “student-t” statistic with c==n - p degrees of freedom and at a risk level of a, Estimated s.e. ({J is the square root of the ith diagonal term of matrix (ZiZ,)- ’ in equation (11). This error estimate has a very small computational cost since the inverted matrix in equation (11) already exists. It can be used to assess the accuracy of the estimated parameters based solely on the data taken during an experiment Thus, it is most useful. In calculating this error estimate, the tacit assump tion is made that the DLT parameters are error free. Since, of course, this is not the case, these error estimates are changed in an unknown and complex manner. To obtain an idea of theefkt of the estimated DLT parameters, one can imagine 8 non-linear least squares approximation involving simultaneously the DLT parameters and the spatial parameters. The matrix Z, in equation (9) now includes partial derivatives with respect to each DLT parameter. In addition, vector F, will contain the residuals from fitting the DLT parameters for each camera. An equation analogous to equation (11) can be used to obtain an augmented variance-covariance matrix. Naturally, this is computationally expensive. Fortunately, as will be shown, the resulting error estimate is not significantly different from that obtained using equation (11). THE REDUCIION OF THE NUMBER OF PARAMETERS TO BE ESTIMATED BY THE USE OF GEOMETBK CONSTBAlNTS
Frequently in biomechanics, one is concerned with the measurement of the spatial parameters of a system of coupled bodies. It also frequently happens that one knows quite accurately the position of one body and wishes to find the position of a body attached to it by a joint. An example of this type occured in the study of tennis carried out by McLaughlin and Miller (1980). In that study, the position of a body consisting of the tennis racket and subject’s hand were relatively well known (having been measured by the methods de-
539
scribed above). It was desired to find the position of the forearm. Now one may consider the wrist to be approximately a Hooke or “Universal” joint having two degrees of freedom. Thus, to find the location of the forearm, one needs only estimate the two Hooke joint parameters. This can be done by a simple extension of the techniques described above. In this case, T*i = T,,@,(OT,,,
(12)
where T,,i is the 4 x 4 transformation matrix of the type described in equation (6) and used in equation (7). T,j is the known transformation from the fixed coordinate system to the joint in question. It may or may not have been found by the techniques described above. @At) is the 4 x 4 transformation for the joint type in question. For a Hooke joint, t represents two parameters consisting of the angles of the joint. T,! is a fixed known transformation relating the local coordinate system at the joint to the system in which the local coordinates on the body are known. Equation (12) can be substituted into equations (7) and a set of equations (8) can be written where the parameters in A,, andJ;1, sre only those of the joint. The Taylor series expansion will only involve the parameters for the joint. A{ in equation (9) wiii also only have as many elements as joint parameters. An iterative solution using equation (10) is used as above. The variance-covariance matrix of the estimpted joint vsri8bles can be found from an equation of the same form as equation (11). Finally, the confidence interval for the joint parameters can be calculated exactly as described above for the spetiai parameters. If one plans to use analysis techniques based upon 4 x 4 methods, such as the MAN program (see Williams, 1977), the estimates of the joint parameters will be sufficient. If, however, one wants the complete set of spatial parameters for the second body, one may use the estimated joint parameters in equation (12) and extract the six spatial parameters from the definition of the transformation matrix. One may estimate the variance-covariance matrix for the six spatial parameters by first constructing equation (9) using F, from the last stage of the evaluation of the joint parameters and using the six estimates of the spatial parameters to calculate Z,. One then calculates the approximate variance-covariance matrix for the six parameters from equation (11). This basic technique can be extended in an obvious way to chains of coupled rigid bodies with any number of joints.
EXPERIMENTALEVALUATION This portion of the paper describes an experiment utilizing a spatial mechanical linkage to evaluate the accuracy attainable using the methods described above.
540
NORMAN R. MILL~R,RO~~RTSHAPIRO andTIR)m M. MCLAUGHLIN
The experimentul linknge While the techniques described in the first part of this paper were under development, a spatial linkage was constructed. This linkage is shown in Fig. 2. A dc. motor with a speed controller drives the input crank which carries a Hooke “universal” joint. This joint is attached to the coupler link which is attached to the output link with a ball and socket joint. The coupler link has true spatial motion. The coupler was, therefore, chosen as the experimental body, The distance from the center of the ball to the amter of the Hooke joint is 25.4ctn This is approximately the length of the human forearm. While many investigators have used light emitting diodes (LED) as photographic markers, it was decided to use small balls, since this type of marker must be used for activities requiring bright illumination. The balls arc 1.8cm in diameter, which, it turns out, is larger than needed. Some error is introduced by assuming that the center of the ball’s photographic image lies at the same point as the image of a geometrical point at the ball’s center would lie. Calculations show that this error is smaller than the rcaolution error of the Vangard Motion Analyzer (VMA) used to quantify the photograph coordinates. The coupler carries 8 photographic markers on metal rods of various lengths. The 4 longest stand 4.8 cm from the centerline of the coupler. The 2 shorter markers at the Hooke joint end arc 3 cm long and the similar pair at the ball joint end are 4.3 cm long. The shorter markers approximately represent markers affixed directly to the surface of the forearm, while the longer markers represent points attached to extenders. A number of other markers arc attached to other parts of the machine. These markers serve as reference points for the computation of the DLT parameters.
1 sec. A frame rate of 100frames per set was more than sufficient for this experiment. Higher speed motions encountered in actual biomechanics experiments will require higher frame rates. The cameras were placed in accordance with the recommendations of Marzan (1976); that is, the cameras were separated by onethird of the camera to subject distance with the photographic axes very slightly convergent. In this case, the camera-sub&t distance was about 2.3 m and the inter-camera distance, 0.8 m. The important point here is that this is not critical. The first step of the experiment was to photograph the mechanism while in a stationary reference position. This reference position (which was set by a mechanical stop) was such that the Hooke joint was at its topmost position. The cameras were not moved after photographing the mechanism in this position. The frames shot at this reference position were later used to establish the camera space to object space DLT parameters. The mechanism was then put into motion and was filmed simultaneously by the two cameras during severalcycles of motion. A potentiometer was attached to the input crank shaft and simultaneously fired, by means of a voltage comparator circuit., the light emitting diode event markers in each camera, Unfortunately, as will be discussed later, this only served to approximately synchronize the cameras to each other and the linkage. After processing, point data was collected from each f&n&p usinga VMA.The data was input to and stored on files on the campus CYBER 175 computer. The frame rates of the two cameras were similar, but not identical.One of the filmstripswas chosen as a reference and the second filmstrip was artificially synchronized to it by cubic spline interpolation. This process is an obvious source of error.
Experimental techniques
Data were gathered by filming the linkage, using two Locam, pin registered high-speed cameras loaded with Kodak RAR 2498film on an Estar-AH base. Both cameras were set to film at 100 frames per sec. The speed of the motor driving the linkage was adjusted so that the linkage completed one cycle in approximately
Hooke Joint
r4-
-&put Fig. 2(a).
Lir 4
Schematicdiagramof experimentallinkage.
Data reduction The position of all the photographic markers in the reference positi& was calculated using the standard kinematic techniques of Uicker et al. (1964). The DLT parameters were calculated using various combinations of marker points. The routine reported in Marzan (1975, 1976) was used for these computations. The DLT parameters from both cameras were used to back calculate the positions of the marker points allowing errors to be calculated. In addition, as will be described below, several sets of DLT parameters were used to calculate the position of the coupler as a function of time. From these results, it was concluded that the use of the 8 points on the mechanism’s coupler, along with the 2 points on the arm carrying the ball joint, produced the best overall result. Using 10 points to calculate the DLT parameters is a bare minimum. More marker points should have been provided in the region near the coupler. This points out the advisability of providing markers throughout the region in which a body is to move.
Fig. 2(b). Experimental spatial linkage in reference position.
541
segments of biomechanical systems
The root mean square errors in estimating the position of the marker points using the 10 point DLT parameters are, respectively, for X, Y and Z, 0.176 cm, 0.428 cm and 0.468 cm. Clearly, these errors are’large. This is rather surprising in view of the excellent results obtained by Professor Karrera and his students using 35 mm still cameras. (Their results are typically two orders of magnitude more accurate.) Any technique of this type is critically dependent on the equipment used. Problems with cameras, film, comparitor or the skill of the personnel collecting the data will adversely influence the results. It is unclear where the weak link in this chain lies. Obviously, the poor precision afforded by the estimated DLT parameters will adversely influence the balance of the experimental results. We shall see, however, that we will be able to estimate the position of the body to a higher accuracy. The DLT algorithm was used to estimate more than the basic 11 parameters. This procedure was unsuccessful. The program provided by Prof. Karrera does not seem to work for more than 11 DLT parameters. A program was written which used Karrera’s 11 DLT parameters to estimate starting values for the 10 parameters implicit in equations (3). The program then used the modified algorithm of Levenberg-Marquardt (Brown and Dennis, 1972) to obtain true least squares estimates of the 10 parameters implicit in equation (3) and then to estimate the various parameters of equation (4). It was found that this procedure was expensive computationally and only provided a limited reduction in the error of the estimate of the position of the marker points as calculated using the photographic data from the two cameras. As a result, it was concluded that inaccuracy was due to errors in the data taken at the reference position and not due to deficiencies in the DLT algorithm. The 11 DLT parameters originally estimated for each camera were utilized throughout the balance of the computations. A series of Fortran programs were written for the CYBER 175 computer system to perform the computations described in the first part of this paper. The results of these computations follow.
543
ball joint and the Y axis is directed roughly downward when the hook joint is in its topmost position. Figure 3(a) shows the location of the coupler center as a function of time, as computed from kinematic theory. The mean coordinates of the coupler position have been extracted to produce plots centered about zero. Figure 3(b) shows the corresponding computed Eulerian angles. Again, the mean angles have been extracted from the data prior to plotting. These plots serve to indicate the range of motion of the linkage. The techniques described in the first part of this paper were next used to estimate the spatial parameters from the cinematographic data. The errors between these estimated values and the calculated values shown in Fig. 3 were evaluated. Figure 4 shows a typical plot of error vs frame. In addition, the 95% confidence interval about the estimated value of the parameter was computed for each frame for each parameter, as discussed above. The theoretical value of the parameter was subtracted from the upper and lower bound on this interval, and the results plotted in the figure. If we had no knowledge of the true value of the parameter, we would hope that the 95% confidence interval about the parameter estimate would embrace the true parameter value. In our plot, we would hope to see these bounds embrace zero error. We see that in general this condition is met. Figure 5 shows the same error plot as Fig. 4 with the 95% confidence region calculated using the more complete estimate of the variance described above. As can be seen, this confidence region is not notably different from the one obtained using the simpler method. Since this technique more than doubles computation time, the first estimate of the confidence region is to be preferred. The data plotted in both Figs. 4 and 5 were computed using all 8 markers carried by the coupler and data from both cameras. Several additional runs were made in which data from : (1) 8 markers but only the left-hand camera were Used;
(2) 8 markers but only the right-hand camera were Used;
RESULTS
The theoretically correct transformation from the ground coordinate system to a coordinate system moving with the coupler of the linkage was calculated to correspond to each set of photographic data throughout the cycle. Calculations were performed as described in Uicker et al. (1964). The ground coordinate system has its center at the point where the input crank shaft projects from its bearing. The X axis is directed downward, the Z along the input shaft toward the middle of the apparatus and the Y axis in the proper direction as defined by the right-hand convention. The coordinate system attached to the coupler has its origin at the geometric center of the link, the local X axis runs from this center toward the
(3) only the shortest 4 markers on the coupler and from both cameras were used ; (4) all eight markers and both cameras were used, but the position of the input crank was calculated so only the two Hooke joint parameters were estimated. The true root mean square error over one cycle for each parameter under each condition is indicated in Table 1. As can be seen, the error for the left-hand camera seems to be significantly larger than the error for the right-hand camera. Notice that the results with 4 markers on the coupler are not significantly worse than those for 8 markers. This is surprising and probably reflects errors in synchronization between the two cameras and the mechanism. Finally, the results from the estimation of only the two Hooke joint
NORMAN R. MILLER, Roem
544
Frame
SHAPIROand THOMAS M. MCLAUGHLIN
Number
Fig. 3(a). Position of coupler center as computed from kinematic theory.
parameters are extremely good. The error in x is zero, since specifying the input crank position fixes this parameter. This section will be concluded with a few words concerning the synchronization problem alluded to above. Synchronization was achieved by firing the light emitting diode event markers in both cameras using a potentiometer on the input crank as a trigger. Unfortunately, due to the pin registered drive of the cameras and the resultant intermittent motion of the film, there is no way of knowing to within less than one frame, the synchronization of the two cameras to each other and the mechanism. The Locam Camera can be obtained with an option to output a signal indicating shutter opening. With this information, one could largely overcome the synchronization problem. Best of ah, it is now possible to obtain phase-locked shutter systems. This would allow nearly perfect synchronization. COMPUTATlON OF SECOND TIME DERIVATWES OF THE PARAMETERS
pfl.564rad
,,...”
-41 0
_.’
6= 4.697rad I
20
40
60
Correspandii Frame
80
loo
I20
Numtw
Fig. 3(b). Angular position of coupler as computed from kinematic theory.
Table
Usually the object of measuring the spatial parameters of the segments of the musculo-skeletal system is to 6nd their tirst- and second-time derivatives and use them in the calculation of kinetic energy or in the computation of forces and moments using the Euler equations. Many smoothing and differentiation techniques are in common use. Experimental position data were available in this study and the theoretically correct second time derivatives of the spatial parameters could be calculated. This being the case, it was decided to evaluate the accuracy of several smoothing and differentiation schemes. Note that error estimates for each parameter at each frame are available, thus weighted least squares techniques can be used to fit smoothing functions. Figures 6(a) and (b) show the calculated second derivative of each parameter. This computation involves a rather direct extension of the techniques of
1. Root mean square error over one cycle for various types of data reduction
dondition
RMS error over cycle a (rad) z, (cm)
B (rad)
Y(rad)
0.019
0.014
0.010
0.692
0.023
0.017
0.014
0.295
0.442
0.018
0.014
0.012
0.080
0.117
0.160
0.019
0.011
0.008
0.066
0.071
0.044
0.0
0.009
0.004
x0 (cm)
Y0 (cm)
0.073
0.096
0.149
left camera only
0.106
0.947
8 markcrsright camera only
0.046
4 markersboth cameras 8 markersestimated only two Hook joint parameters
8 markers on coupkrboth cameras 8 markers-
545
Segments of biomechanical systems
_ *r
30
_ ‘.
4-
.::
‘. ..
. ‘.
..
;.
.’
..‘..
.., ..IJpper
of 95
‘.
,_
.;
I
r
Bound %
.‘-’ Conhdence
.,
Lower Bound
-6 0
20
40
80
60 Frmne
, 120
loo
-30
1
0
Number
I 20
40
Cornrpondmg
Fig. 4. Error in YOparameter as a function of frame number showing the estimated error band based on the 95% confidence interval about the estimated parameter value.
Upper Bound
40
80
60
Frame
of 95% Confidence
‘.
._..’
IO0
120
Number
Fig. 5. Error in YOparameter as a limction of frame number showing the estimated error band based on the 95% confidence interval about the estimated parameter value (error band based on more exact estimate). 60
I
I
60
80
Frame
I
I
100
120
Number
Fig. 6(b). Calculated angular acceleration of coupler.
Denavit et al. (1965). Basically, one calculates the firstand second-time derivatives of the global coordinates of three points on the body and writes the first total time derivatives of equation (6). He then solves for the first time derivatives of the spatial parameters. Differentiating quation (6) again, one obtains quations for the second time derivatives. Various techniques for smoothing and numerical diffematiation were applied to the #mtographic po sition data. They include the following: (1) simple second-order finite differences; (2) digital filtering wing a Butterworth filter followed by finite di&rences (Winter et aL, 1974); (3) fitting of polynomials by weighted linear least squares, followed by differentiation of the polynomials ; (4) fitting of a Fourier series via weighted linear least squares followed by differentiation of the resulting Fourier series (Cappozzo et al., 1975); (5) fitting of a cubic spline using the weighting technique described in Reinsch (1967) and McLaughlin et ol. (1977). The weights for the spline estimate were obtained as follows : If (s.e.), is the standard error of the parameter estimate for a given’ parameter and frame, D, is the difference between the spline value and data value at a frame and NF is the number of frames in the data set. Then 5
[D&e.),]’
I NF.
I=1
-60
I 0
20
40 Correspondlng
60 Frame
80
too
120.
Number
Fig. 6(a). Calculated acceleration of coupler center.
Table 2 lists the root mean square error over a cycle for each parameter’s second derivative using each smoothing technique. Notice that the digital filter, the Fourier series and the cubic spline all work well. The Fourier series is not a good technique to choose for smoothing a. (Indeed, the best technique would be to fit a straight line.) It would appear that different techniques are applicable
546
NORMAN
R. MILL& Rotwtr SHAPIRO and Txoh~ls M. MCLAUGHLIN
Table 2. Root mean square error in second parameter time derivatives over one cycle for various smoothing and differentiation techniques RMS error over cycle 20 (radiec’) ~cm/sec’)
x0 (cmkc’ ) Second order tinite ditkrences
316
Digital filmring (2 Hz break for linear data; 3 Hz for angular data) Polynomials (3rd dcgrtc)
1258
3.50
6.09
26.7
28.2
1250
32
7.68
.. @ad/se? )
25
1.29
16.6
B (rad/sec’)
0.101
35
1.08
0.96
13.8
11.6
Fourier series (X-2nd harmonic; Y, Z, a-fundamental ; j?,y-5th harmonic )
4.72
1.91
4.07
N.A.
2.90
1.40
Cubic spline
3.66
4.72
5.41
0.95
1.02
1.11
to different parameters This simply emphasixes that there is no one “best” smoothing technique and that plots and common sense are indispensible if one wishes to obtain reasonable results. Figures 7(a) and (b) show the error resulting from use of the cubic spiine as a function of frame number. These error plots should be compared with the calculated accekrations in Fig 6. Notice that the magnitude of the maximum error in X0 is about 8cm/sec2, while the piak magnitude of x,, is about 55 cm/sec2 ; thus, the ratio of peak error to peak magnitude is g/55 - 0.15. As one can see, while the accuracy is not exceptional, one can obtain an order of magnitude evaluation of acceleration using this technique. Further, these results can probably be taken as an indication of the minimum accuracy obtainable by these methods.
2-
-2 . *'
3O
20
40 60 80 Frame Number
NM
I I20
Fig. 7(b). Angular acceleration error - cubic splint
CONCLUSIONS
-I”
0
20
40 60 80 Frame Number
100
I20
Fig. 7(a). Linear acceleration error - cubic spline.
A technique for photographically obtaining the spatial kinematic parameters of a rigid body has been described. A minimum of one camera and three points on the body are required. Accuracy can be improved by the use of multiple cameras. The method provides estimates of the errors on the calculated parameters. The method extends in a natural way to estimating the positions of coupled chains of rigid bodies. In this case, a reduced set of parameters needs to be estimated. A controlled experiment was conducted to study the accuracy of the method. Very good positional accuracy was achieved and a reasonable level of acceleration accuracy was obtained by the use of several smoothing
and differentiation
techniques.
547
Segments of biomechanicai systems
Draper, N. and Smith, H. (1966) Applied Regression Analysis. Wiley, New York. Erdman, A. G., Forstrom, R., Dorman, F., Ah&en, F:, Beer, R., Isaacson, R. J. and Speidel, M. T. (1975) Experimental motion analysis via a stereo-photography technique: ap plication to the human jaw and mechanical linkages, ASME Paper No. 75-DET-58. Greenwood, D. T. (1965) Principles OfDynamics, pp. 333-335. Prentice-Hail, Englewood Cliffs. Marxan, T. and Karrara, H. M. (1975) A computer program for direct linear transformation of the colinarity condition and some applications of it. Symp. on Close Range Photogrammetric Systems. American Society of Photogrammetry, Fails Church. Marxan. T. (1976) Rational design for close-range photogrammetry. Ph.D. Dissertation. University of Illinois at Urbana-Champaign. McLaughlin, T. M., Diiiman, C. J. and Lardner, T. J. (1977) Biomechanicai analysis with cubic spiine functions. Res. Q. 48, 570-582. McLaughlin, T. M. and Miller, N. R. (1980) Techniques for the evaluation of loads on the forearm prior to impact in tennis strokes. J. Mech. Design Trans. ASME, in press. Miller, D. 1. and Petak, K. (1973) Three-dimensional cinematography. Kinesiology III. American Association for Health, Physical Education and Recreation, Washington. Reinsch, D. H. (1967) Smoothing by spiine functions. h’um. Math. 10, 177-183. Shapiro, R. (1978) The direct linear transformation method for threedimensional cinematography. Res. Q. 49, 197-205. Suh, C. H. (1974) The fundamentals of computer aided X-ray analysis of the spine. J. Biomechanics 7, 161-169. Uicker, Denavit, J. and Harm&erg, R. S. (1964) An iterative method for the displacement analysis of spatial mechanisms. J. appl. Mech. 31. Trans. ASME, Se&s E, g6, pp. 309-3 14. Williams, R. and Seireg, A. (1977) Intaractive computer modeling of the muscuio-skeletal system. IEEE Trans. biomed. Engr. BME-24,213-219. Winter, D. A., Sidwail, H. G. and Hobson, D. A. (1974) Measurement and reduction of noise in kinematics of locomotion. J. Biomechanics 17, 157-159.
Acknowledgements - The authors would like to thank Professor Karrera of the Department of Civil Engineering, University of Illinois at Urban-Champaign, for supplying us with a copy of the program used to estimate the DLT parameters (Marxan and Karrera, 1975).
REFERENCES Abdei-A&, Y. I. and Karara, I-l. M. (1971) Direct linear transformation from comparator coordinates into object space coordinates in close-range photogrammetry. ASP Symp. on Close-Range Photogrammetry. American Society of Photogrammetry, Fails Church. Ayoub, M. A., Ayoub, M. M. and Ramsey, J. D. (1970) A steriometric system for measuring human motion. Hum. Factors 12, 523-535.
Brown R. H., Burstein, A. H., Nash, C. L. and Shock, C. C. (1976) Spinal analysis using a thradimensionai radiographic technique. J. Biomechanics 9, 355-365. Brown, K. M. and Dennis, J. E. (1972) Derivative free anaiogues of the Lever&erg-Marquardt and Gauss aigorithms for nonlinear least squares approximations. Num. Math. 18, 289-297. Bullock, M. and Harley, I. A. (1972) The measurement of three-dimensional body movements by the use of photogrammetry. Ergonomics 15, 309-322. Cappozxo, A., Leo, T. and Pedatti, A. (1975) A general computing method for the analysis of human locomotion. J. Biomechanics 18,307-320. Chao, E. Y. and Morrey, B. F. (1978) Three-dimensional rotation of the elbow. J. Biomechanics 11, 57-73. Crowninsbield, R. D, Johnston, R. C., Andrews, J. G. and Brand, R. A. (1978) A biomechanicai investigation of the human hip. J. Biunechanics 11, 75-85. Denavit, J. and Harms&erg, R. S. (1955) A kinematicnotation for lower pair machanisms based on matrices. J. appl. Mech. 22. Trans. ASME T7,215-221. Denavit, J., Harm&erg, R. S., Rti and Uicker. (1965) Velocity, acceleration, and static force analyses of spatial linkages. J. appl. Mech. 32. Trans. ASME. Series E, lt7, 903-Y 10.
APPENDIX: The Tran$ormation
Matrix in Terms of Eulerian Angles
1 cos u,, cos ya,
-
cosaIl m
Y#I
sin a,, sin &, -sin apl cos /I,, sin yI, sm a#, cos Y,~
-sin aII cos &, cos ye, - sm aorsm ysr - cos a,, sin &,
T,, = + cos e,t cos B,r sin ysl sin /J,, sin ys,
+ cos a,1 cos 8,, cos Y#I sin B,t cos Y,I
cos B,,