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.

DIFFERENTIATION

I 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 VELOCITY

1. 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,

coming soon on the Journal of Biomechanics; DOI:

10.1016/j.jbiomech.2007.09.030; when published on line, 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).

With kind regards,

Paolo de Leva

-----Messaggio originale-----

Da: * Biomechanics and Movement Science listserver

[mailto:BIOMCH-L@NIC.SURFNET.NL] Per conto di Rettig, Oliver

Inviato: venerdì 9 novembre 2007 15.48

A: BIOMCH-L@NIC.SURFNET.NL

Oggetto: Estimation of angular velocity/acceleration from 3x3 matrices

Hallo,

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. 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 regards

Oliver Rettig

---

Oliver Rettig

Stiftung Orthopädische Universitätsklinik Heidelberg

Ganglabor

Schlierbacher Landstr. 200a

69118 Heidelberg

Germany

Tel: +49 6221-96 6720

Fax: +49 6221-96 6725

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.

DIFFERENTIATION

I 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 VELOCITY

1. 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,

coming soon on the Journal of Biomechanics; DOI:

10.1016/j.jbiomech.2007.09.030; when published on line, 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).

With kind regards,

Paolo de Leva

-----Messaggio originale-----

Da: * Biomechanics and Movement Science listserver

[mailto:BIOMCH-L@NIC.SURFNET.NL] Per conto di Rettig, Oliver

Inviato: venerdì 9 novembre 2007 15.48

A: BIOMCH-L@NIC.SURFNET.NL

Oggetto: Estimation of angular velocity/acceleration from 3x3 matrices

Hallo,

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. 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 regards

Oliver Rettig

---

Oliver Rettig

Stiftung Orthopädische Universitätsklinik Heidelberg

Ganglabor

Schlierbacher Landstr. 200a

69118 Heidelberg

Germany

Tel: +49 6221-96 6720

Fax: +49 6221-96 6725