University of Delaware



Differentiation of biomechanical dataWCR 2007-11-14This was stimulated by Oliver Rettig’s 2007-11-09 posting to Biomech-L.I am interested in how to calculate angular velocity/acceleration from a time series of 3x3 rotation matrices which I get from marker positions measured by a Vicon motion analysis system. Different ways are possible and I want to learn which of the methods are used in practice and which numerical advantages or disadvantages have the different methods. The following possiblities are my startpoint:1. First calculate time series of cardan angles with an arbitray rotation order. Then you can estimate the derivations of these three time series of double values by a) simply calculating differences or by b) low order polynom derivations (a kind of filtering and differentiating in one step, different orders and techniques are possible). One disadvantage of this method can be gimbal lock.2. Use the formula [~omega]=[M'(t)][M(t)]-1 to calculate the Tensor which includes the components of the angular velocity by time derivatives of each component of the rotation Matrix M(t) multiplicated with the inverse of the rotations matrix. For the estimation of the time derivative of the component the same methods from 1. a),b) can be used. In Comparison with 1. gimbal lock should be no problem. 3. You can simply estimate the velocity by calculating dot products of the columns between two sequent matrices. The implementation of Vicon PiG seems to do this by using matrices of frames with a time distance of 0.05s by 120Hz frame rate.4. Because my rotations matrices are typically based on cross products of vectors between markers I can analytical calculate formulars of the angular velocity as function from the derivations of the marker positions.Which of these methods do you use? Are there further methods? Which filtering details you use? Are there further practical details I have to look at? I am also interested in literature.best regardsOliver RettigThe segment angular velocity and acceleration vectors (ω and α) are expressed in terms of the segment’s own anatomical axes, x,y,z (as opposed to the laboratory frame, X,Y,Z). This is useful because the equations of motion for a segment are simpler when expressed with respect to the segment’s anatomical axes, originating at the segment center of mass.A widely used technique is described by Winter (3rd ed, 2005), among others: Filter the displacement data (x(t), y(t), z(t) for each marker) digitally to remove high frequency noise. Use the smoothed marker data to calculate 3 rotation angles at each instant (θ1(t), θ2 (t), θ3 (t)). These angles describe the orientation of the segment relative to the laboratory frame, in terms of a specified rotation sequence, such as Cardan, Euler, or other. [In Rettig’s case, use the rotation matrices produced by Vicon, to back-calculate θ1(t), θ2 (t), θ3 (t). Are these rotation matrices based on marker data that has already been filtered? That’s a ? for Vicon.] Calculate angular velocities d θ1/dt, dθ2/dt, dθ3/dt) by simple finite differencing.? Then compute the segment angular velocities (ωx, ωy, ωz) using a matrix equation (Eq. 7.6b in Winter 3rd ed.). Compute the segment angular accelerations (αx, αy, αz) by simple finite differencing of the segment angular velocities.The appropriate cutoff frequency for the initial filtering will depend on the system being analyzed and the laboratory equipment used. A 2nd order Butterworth applied in both directions is standard in biomechanics. The appropriate cutoff frequency depends on the system being analyzed and the laboratory equipment used. Winter suggests using residual analysis to determine cutoff frequency; with this method he finds cutoff frequencies of 3-6 Hz for lower limb markers in walking. My advice: when choosing low pass filter cutoff frequencies, it is better to err on the high side (i.e. choose a high cutoff frequency, such as 10 Hz for walking, 20 Hz or more for running). If the resulting angular velocities or accelerations look too noisy, go to lower corner frequencies (i.e. "stronger" filtering) if absolutely necessary. (The flaw in this advice is the subjective determination of what is “too noisy”.)Other investigators suggest different approaches to obtaining linear velocity and acceleration signals, including: local polynomial fitting followed by differentiation; filtering of displacement data followed by differentiation to get velocity, then filtering again followed by differentiation to get acceleration; differentiating first, then filtering, and differentiating twice, then filtering... See Hatze (1981) and Giakas & Baltzopoulos (1997) for a good comparison of methods. Hatze compute “optimal” filter cutoff frequencies for displacement, velocity, and acceleration. With this method, the raw displacement data are filtered with one cutoff to get smoothed displacement. The raw displacements are filtered with a different cutoff, then differentiated (by simple finite differencing) to get velocity. The raw data are filtered with a 3rd cutoff, then second order finite difference is computed to get acceleration. G&B conclude Hatze’s method is the best; Winter’s “residual” approach is not the best for velocity and acceleration. In these papers linear, not angular, velocities and accelerations are computed and compared. This difference is significant, as explained in the next paragraph.Several respondents to Rettig’s posting (incl. van den Bogert and Sommer) state that the calculation of rotation angles is ill conditioned near gimbal lock (i.e. when second rotation angle is near 90°). This is a disadvantage of Winter’s recommended method. An alternative is to calculate segment angular velocities and accelerations by a matrix method that does not use Euler or Cardan angles (Rettig’s methods 2 and 3; see postings by van den Bogert, de Leva, Sommer, Vrongistinos).Young Hoo Kwon’s reply to Rettig’s question (2007-11-12):I am familiar with methods 1 and 2 but I actually use method 1 to compute angular velocity. It actually has additional steps to carry out. See the details at angles may have gimbal lock issue but one definite advantage of the orientation angle approach is its intuitiveness. We had a discussion on Biomch-L before but if sine of the second rotation angle approaches to 1 (or cosine -> 0), I treat it as missing and generate the orientation angles later through interpolation for the missing interval. Some people use different rotation sequences at different joints to avoid gimbal lock as well but I personally use XY'Z" or ZX'Y" sequence only, where X is the ML axis, Y is the AP axis, and Z is the longitudinal axis.I filter the orientation angles by using a Butterworth filter before computing the angular velocity. Angular acceleration is simply a time derivative of the angular velocity.Young-Hoo Kwon, Ph.D., Director, Biomechanics Laboratory, Texas Woman's UnivPaolo de Leva’s reply (2007-11-12):Let's call R = [i j k] the 3x3 orientation matrix, where i, j, k are column vectors. Notice that, according to the most widely used convention, i, j, k are the positions of the tips of the versors of the local (technical or anatomical) coordinate system L (versor = unit vector codirectional with a directed axis), represented in the global coordinate system G. DIFFERENTIATIONI would use smoothing quintic spline functions to compute the derivative of each element of R. This is because spline functions can be differentiated very easily (they are piecewise polynomials). With MATLAB, for instance, you can use SPAPS to build the nine spline functions, and FNDER to obtain their first and second derivatives. These derivatives are the linear velocities and accelerations of the tips of the versors. In my opinion, quintic is preferable to cubic spline (this is related to the fact that, in short, quintic spline smooths by reducing jerk, while cubic spline does it by reducing acceleration).COMPUTING ANGULAR VELOCITY1. Why should you first convert R into Cardan angles, if you can use directly method 2?2. [omega]=[M'(t)][M(t)]-1 is called Poisson's equation, as far as I know. And it does not yield an approximation of omega ([~omega]). The approximation comes from the values of M' and M. Obviously, [~omega]=[~M'(t)][~M(t)]-1. This is the correct method, together with the equivalent method which uses a simple linear combination of three cross divisions:2b. [omega]= 0.5 * (i'/i + j'/j + k'/k) = 0.5 * (i x i' + j x j' + k x k'), where "/" stands for orthogonal cross division. But my article on "Anticrossproducts and cross division" has not been published yet (in press, J. Biomech; DOI: 10.1016/j.jbiomech.2007.09.030; I invite you to read sections see sections 5.3. and 5.4, and figures 3 and 8, which in my opinion very nicely show the geometrical interpretation of the angular velocity vector, and its relationship with the linear velocity of the tips of i, j, k).3. This PiG method can only give the cosines of three angles that are not Euler or Cardan angles. How can you compute the angular velocity vector from that? For instance, if by computing the arccosine of the thee dot products you discover that i and j rotated by 0.05 rad in 0.05 s, and k rotated by 0 deg, then you can say that the angular velocity was 1 rad/s about the z axis. But if all the versors rotated by a nonnull angle, then I guess this method just doesn't work.4. The cross products (and the ensuing normalization needed to obtain versors i, j and k) amplify or attenuate the stereophotogrammetric error. In the past, a lot of papers dealt with the problem of smoothing before or after differentiation, and before or after center of mass and "joint center" position estimation. Estimating R means estimating the position of three points, at a distance of 1 mm or 1 m from the origin of L (depending on the unit you are using for length). These points are not real markers, they are "virtual" markers. Whatever you decide, you should consider that the noise in the position of virtual markers is different from the stereophotogrammetric error (error in the reconstructed position of real markers).Paolo de LevaTon van den Bogert’s reply (2007-11-12):> 1. First calculate time series of cardan angles with an arbitray rotation order. [...]. > One disadvantage of this method can be gimbal lock.Young-Hoo Kwon already cited his page which gives the equations to obtain the angular velocity vector from cardan angle derivatives. Near gimbal lock, these derivatives will be very noisy. I suspect (without proof) that this noise will disappear again after going through the angular velocity equations, but still... You have to set a threshold to define when you are so close to gimbal lock to treat this as missing data etc.This method is actually what I generally use in a linked multibody model, you differentiate the generalized coordinates and use forward kinematics equations to get the angular velocities of the body segments.When using this method, it is best to define kinematic variables (generalized coordinates) in a way that gimbal lock is avoided in your particular application.For any differentiation, I would not recommend polynomial fitting but use a proper low pass filter such as Butterworth digital filter, or splines as Paolo de Leva suggested. > 2. Use the formula [~omega]=[M'(t)][M(t)]-1 to calculate the > Tensor which includes the components of the angular velocity In my experience, this works well (Bellchamber & van den Bogert, J Biomech 2000). Still it seems less elegant to differentiate 9 signals when there are only 3 rotational degrees of freedom.> 3. You can simply estimate the velocity by calculating dot > products of the columns between two sequent matrices. The If that is indeed mathematically correct, it must be done after appropriate smoothing as in method 2. Then the results become independent of the frame rate and resolution.> 4. Because my rotations matrices are typically based on cross > products of vectors between markers I can analytical > calculate formulars of the angular velocity as function from > the derivations of the marker positions.Yes, intuitively it should be optimal to estimate angular velocity from marker velocities directly. In the past, I have suggested the following (van den Bogert, Exerc Sports Sci Rev 1994):(1) Obtain marker velocities from raw data via smoothing and differentiation(2) Rigid body model for marker velocity v_i as a function of segment velocity v and angular velocity omega:v_i = v + omega x (p_i - p)where p us the position of the segment origin and p_i is measured marker position.(3) With N markers, these are 3N linear equations with 6 unknowns: v and omega.(4) Write the equations as A*x = b, where A is a 3N x 6 matrix and solve x = (v, omega) using linear least-squares (Matlab: x = A\b)I believe (without proof) that this method would give the least error in angular velocity. And there are no singularities.If you used your analytical expressions, you would be essentially using only 6 of the 9 equations (when you have three markers). It is better to use all information in the marker data to minimize the effect of non-rigidity of the marker set. Also the least-squares numerical method extends nicely to using more than 3 markers on a segment. These are the same reasons as for using least squares to estimate segment position p and rotation matrix R (e.g. ).I suspect that, in practice, there is little difference between the methods you listed, but it would be nice to have that confirmed. Especially in certain "pathological" situations (e.g. near gimbal lock, near-colinear marker placement, high noise, ...Ton van den Bogert, Department of Biomedical Engineering, Cleveland ClinicH.J. (Joe) Sommer’s reply (2007-11-13)> -----Rettig’s original Message----- > I am interested in how to calculate angular > velocity/acceleration from a time serie of 3x3 rotation > matrices which I get from marker positions measured by a > Vicon motion analysis system. > 1. First calculate time series of cardan angles with an > arbitray rotation order. Then you can estimate the > derivations of these three time series of double values by a) > simply calculating differences or by b) low order polynom > derivations (a kind of filtering and differentiating in one > step, different orders and techniques are possible). One > disadvantage of this method can be gimbal lock. Cardan/Bryant/Euler angles are non-linear functions of rotation matrices.? Taking derivatives of non-linear functions particularly near ill-conditioned values (gimbal lock) will magnify noise.> 2. Use the formula [~omega]=3D[M'(t)][M(t)]-1 to calculate the > Tensor which includes the components of the angular velocity > by time derivatives of each component of the rotation Matrix > M(t) multiplicated with the inverse of the rotations matrix. > For the estimation of the time derivative of the component > the same methods from 1. a),b) can be used. In Comparison > with 1. gimbal lock should be no problem.This provides a first-order finite-difference first-derivative estimate of angular velocity with no inherent filtering.> 3. You can simply estimate the velocity by calculating dot > products of the columns between two sequent matrices. The > implementation of Vicon PiG seems to do this by using > matrices of frames with a time distance of 0.05s by 120Hz frame rate. This is essentially the same as method 2 above.> 4. Because my rotations matrices are typically based on cross > products of vectors between markers I can analytical > calculate formulars of the angular velocity as function from > the derivations of the marker positions.Please see Sommer, H. J.? 1992, Determination of First and Second Order Instant Screw Parameters from Landmark Trajectories. ASME J. Mechanical Design 114(2):274-282. Similar methods have been discussed by Angeles and Brodeur/Soutas-Little.This method uses time derivatives of marker motion and provides linear least-squares equations for angular velocity and angular acceleration.? Because you can use higher order derivative estimates or splines applied directly to raw data and not to ill-conditioned non-linear intermediate functions, this method provides improved linear estimates. We have recently extended the method to include sensitivity equations, angular jerk and experimental measurement of axode invariants. Joe H.J. Sommer III, Ph.D., Professor of Mechanical EngineeringDino Vrongistinos’ reply (2007-11-13):The three most direct methods are --- as in Wittenburg (1977)(1) As you and others mentioned either use tensors (solved here) If world base e(1) (lab coordinates) and a rotating base e(2) attached to segment r is attached to e(2) measured as r(2) -- most times the segment rdot = first derivative of r A12dot = first derivative of A12 A12 = transformation matrix that you multiply r(2) to transform as r(1) in Lab coordinates OMEGA = skew symmetric matrix r(1)=A12 * r(2)? then rdot(1)= A12dot*r(2) or rdot(1)=A12*(dr/dt)=A12 * (w x r(2)) = A12 * OMEGA(2) * r(2) A12dot*r(2)=A12*OMEGA(2)*r(2) => A12dot=A12 * OMEGA(2) => OMEGA(2)=A21*A12dot So if you multiply transformation matrix A21 (transformation matrix to convert a vector from e(1) to e(2)) with the dotA12 (first derivative of transformation matrix to convert a vector from e(1) to e(2) your matrix is then 0-w3w2w30-w1-w2w10=A21? A12dotNo need for specific angles -- Cardan etc (2) If you do not need to calculate any other angular kinematics:Using three markers (or the three less noisy markers and after smoothing with your preferred method)r1 r2 r3 markers in 3D (r1 r2 and r3 should be measured once in static condition otherwise the r errors will give wrong omegas) u1, u2, u3 the linear velocities of the three markers in 3D (x is cross product, . is dot product) w = 2*(u1xu2+u2xu3+u3xu1)/(u1.(r2-r3)+u2.(r3-r1)+u3.(r1-r2)) No need for any angular quantities -- even better (3) If quaternion Q=[q0,q1,q2,q3] represents the transformation matrix equivalent q0 = real part or magnitude (or in practice the angle you are rotating) (q1,q2,q3) = imaginary part (or in practice the axis around which you are rotating)Then angular velocities are (q dot is the first derivative of each parameter in terms of time)0w1w2w3=2q0-q1q1q0q2q3q2q3-q2-q3q0 q1-q3q2-q1 q0q0dotq1dotq2dotq3dotKonstantinos "Dino" Vrongistinos, Ph.D.Next reply (2007-11-): ................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download