Getting angular velocities

[opherdonchin]

Following examples in the documentation, I had no difficulty taking my marker data and turning it into frames, and from frames into Eueler joint angles. I'm sure that I can use numpy to get a derivative if I want, but it would be great to know if there is a best practice way to do this within the toolkit and whether there are any gotchas I should be aware of.

Thanks,
Opher

[felixchenier]

Hi @opherdonchin

I'm glad to know that you got your joint angles easily. Now, for joint angular velocities, I have a small good news, but a more serious bad news for you.

Good news: It is very easy to derivate a TimeSeries directly in Kinetics Toolkit using ktk.filters.deriv.

Bad news: In almost every case, we cannot simply derivate the Euler angles to obtain the angular velocity. The three Euler angles are highly dependent one to the other. Therefore, as changing one Euler angle usually affects at least one another, then the angular speed on a given axis is also dependent of the angular speed on the other axes. I don't know of any method to obtain joint velocities that is simple, intuitive, and general. You can take a look at this paper where they calculated these velocities first by calculating a velocity matrix, then mapped this matrix to the axes corresponding to their Euler angle convention:

Slawinski, J., Bonnefoy, A., Ontanon, G., Leveque, J.M., Miller, C., Riquet, A., Chèze, L., Dumas, R., 2010. Segment-interaction in sprint start: Analysis of 3D angular velocity and kinetic energy in elite sprinters. Journal of Biomechanics 43, 1494–1502. https://doi.org/10.1016/j.jbiomech.2010.01.044

To my understanding, this method needs the three rotation axes to be orthogonal, which means it works for cardan angles (XYZ, XZY, etc.) but would not make sense for proper Euler angles (ZYZ, XZX, etc.) in which case it would be very difficult to interpret anyway.

Good luck! 😁

[opherdonchin]

Hi Felix,

After a little perusal of the paper you sent, it looks like the the theory and some applications for deriving homogeneuos velocity or acceleration matrices are worked out in two papers I reference below. It doesn't seem like figuring out the required matrix algebra should be too hard, but what isn't clear to me is whether I can use your function get_euler_angles on the velocity and acceleration matrixes and get the angular velocity.

Do you have insight on that?

Thanks,
Opher

Giovanni Legnani, Federico Casolo, Paolo Righettini, Bruno Zappa,
A homogeneous matrix approach to 3D kinematics and dynamics — I. Theory, Mechanism and Machine Theory. Mechanism and Machine Theory 31(5):573-587, 1996 https://doi.org/10.1016/0094-114X(95)00100-D.

Giovanni Legnani, Federico Casalo, Paolo Righettini, Bruno Zappa, A homogeneous matrix approach to 3D kinematics and dynamics—II. Applications to chains of rigid bodies and serial manipulators, Mechanism and Machine Theory 31(5):589-605, 1996.
https://doi.org/10.1016/0094-114X(95)00101-4.

[felixchenier]

@opherdonchin I love these kinds of discussion. Everybody learns! By the way, thank you a lot for these papers.

OK let's first rule out the get_euler_angles function. This function extract angles from a homogeneous transform, using trigonometric equations. In your case, you don't have a homogeneous transform but a velocity matrix, which don't use sines and cosines but arcs instead. Therefore get_euler_angles will definitely not work.

Be sure to double-check what I say here because I'm quite new to 3D angular velocity calculation. What these papers give you are a velocity matrix W that expresses the angular velocity (wx, wy, wz) and the linear velocity (vx, vy, vz) of a body relative to a given frame. If you calculate W_proximal for your proximal segment in the global reference frame, then W_distal for your distal segment, also in the global reference frame, then you can calculate

W_joint = W_distal - W_proximal

, which will also be expressed in global coordinates.

All you have to do next (I believe) would be to choose which local coordinate system to use to express your angular velocities. For instance, if you want to express velocities in the distal segment coordinates (M_distal, using the notation M for the homogeneous transform as in Legnani's papers), you would use:

W_joint_local = ktk.geometry.get_local_coordinates(W_joint, M_distal)

At this point, you can extract w_x, w_y and w_z from W_joint_local, knowing that x, y and z are now the lateral, longitudinal and anterior axes of the distal bone according to your convention:

w_x = W_joint_local[:, 2, 1]
w_y = W_joint_local[:, 0, 2]
w_z = W_joint_local[:, 1, 0]

I may be plain wrong, so please be sure to check by yourself and continue to investigate, but this is the starting point I would use. This also sparks lots of ideas for new features in Kinetics Toolkit, once I have time to sit and master this important topic.

[felixchenier]

@opherdonchin You know what? I've also been reading around, because you tickled my curiosity, and I think we may be overthinking this problem.

Let's say we work on the shoulder, and we use the tilt-and-torsion method (shameless plug for a paper of mine on a variation of Euler angles to correct for humeral rotation bias 😋: Chénier, F., Alberca, I., Faupin, A., Gagnon, D.H., 2022. Interpreting the tilt-and-torsion method to express shoulder joint kinematics. Clinical Biomechanics 92, 105573. https://doi.org/10.1016/j.clinbiomech.2022.105573

Following this convention, the plane of elevation, elevation, and humeral rotation of the right arm are calculated using Euler's YXY angles:

• $\theta_\text{plane of elevation} = \theta_1(t)$
• $\theta_\text{elevation} = - \theta_2(t)$
• $\theta_\text{internal rotation} = \theta_3(t) - \theta_1(t)$

Then by definition:

• $\omega_\text{plane of elevation} = d(\theta_1)/dt$
• $\omega_\text{elevation} = d(- \theta_2)/dt$
• $\omega_\text{internal rotation} = d(\theta_3 - \theta_1)/dt$

In other words (and other conventions), if you have a TimeSeries that you defined as being the joint flexion, then numerically deriving this TimeSeries necessarily gives the flexion speed. There's no way it would give something else.

It's all about of how you define your outcome measures. In my previous post, the velocity was represented differently, as "the angular joint velocity expressed in distal segment coordinates". Both methods simply have different meanings.

• Deriving the Euler angles may be more intuitive in terms of clinical interpretation (ex.: flexion speed, rotation speed, etc.). However, it is sensitive to gimbal lock and it is less obvious which rotations are around which axis (for proper Euler angles, the three axes are not even orthogonal).

• Finding the velocity matrix and expressing it in one coordinate system of choice (say the distal segment) has a more obvious geometrical interpretation as we know explicitly what are the rotation axes, it is much more practical for calculating joint power as the three velocities are orthogonal and can be easily expressed in the same reference frame that the joint moments, and it's not sensitive to gimbal lock. But it is less related to clinical terms such as flexion, rotation, abduction, etc.

Therefore what I think is that it really depends on what results you want to report and in what aim.

Still, this is an interesting topic!

[opherdonchin]

Hi Felix,

What you say makes a lot of sense.

That is, if my fundamental interest is in the Euler angles than the derivative of the Euler angles seems like the right thing to calculate. If I want the trajectory of my segment in space, then the rotation matrix seems like the right thing to look at.

Here are a few more thoughts and questions:

1. I posted a similar question on Biomech-L and although nobody responded, I did some work looking through the archives and posted my own response. You might be interested in what I found.
2. Off hand, it seems like the derivative of the rotation matrix ought to give you the derivatives of the Euler angles somehow, but I'm still not clear on how you would go between them. This feels like a sign I don't really understand things yet.
3. The general consensus from my reading is that the best way (numerically) to get the derivative of the rotation matrix is actually from markers themselves. That is, to use linear regression to estimate the translation and rotation that would get your from the markers at time t-1 to the markers at time t. It seems to me like it would be kind of fun, and would have a clear general solution to do this, but see #2.

Thanks for the enthusiasm!
Opher

I wonder if there is a general rule for how these equations will work or if I need to go through i