# Calculate angular velocity from rotation matrix difference

I am working on something for a game. I need to calculate the angular velocity, however in my situation I only have access to the previous rotation matrix and the current rotation matrix (and some other variables like dt or age of the object).

My angular velocity is hereby constant between t0 and t1.

The normal formula for the angular velocity W(t) = dR(t)/dt * R(t)^-1 can't be used here, since I have discrete values. Is there a possibility to calculate the angular velocity between timestep t0 and t1 with the given R0 and R1 rotation matrix of an object?

I am not very advanced with matrix math, so please go easy on me.

• I don't have complete control over the whole pipeline and the system wasn't designed with this in mind so I it is very hard for me to go to the code pieces which control the rotation and calculate it from there, especially because I have a parenting object system too. Mar 18 '21 at 7:39

You can get the difference between the two transformations by multiplying the new one by the inverse of the old. (If these matrices are pure rotations, with no translation/scale/shear, then this inverse is just the transpose)

$$R_\text{diffference} = R_1 \times R_0^{-1}$$

Now you can extract the quaternion representing the rotation from this difference matrix. The axis of rotation then points along the imaginary part of the quaternion, and the angle of travel is twice the arc cosine of the real component:

q = MatrixToQuaternion(R_difference)
axis = Normalize(q.xyz)
angle = 2 * Acos(q.w)

Or you can extract the angle and axis directly from the matrix as described here

Now to convert this to an angular velocity, you can simply multiply the axis unit vector by the angle travelled, divided by the time step over which the travel occurred.

• multiply the axis unit vector by the angle travelled this confused me, where does the axis come into play when looking for angular velocity $ω$ ? Mar 18 '21 at 12:22
• Because angular velocity is a vector quantity. It has both a magnitude (speed of spin) and a direction (axis of rotation). Mar 18 '21 at 12:23
• I see what you mean. That actually determines if the velocity is positive or negative, right? Always depending on the choice between "right-handed" or "left-handed" coordinate system. Mar 18 '21 at 12:28
• Sort of. Negating the speed or reversing the axis of rotation both have the effect of reversing the angular velocity vector (as we'd hope they would, since both should have exactly the same effect on the motion of the object). But because we're in 3 dimensions, there's not an unambiguous notion of the "sign" of a velocity as there is in 2D. We can basically treat all speeds as positive, with the direction of the vector giving us a richer representation of direction than just "positive or negative". Now we can represent "signs" like x+, y+, z+, x-, y-, z-, and every mixture in between. Mar 18 '21 at 12:52
• Using the arc cosine is a bad idea for small angles, essentially you are doing sqrt(1-(1-x*x)) in floating point, losing most digits for x<1e-4 and all for x<1e-8. One would be better served with the size of the imaginary part, its arc sine, or combining both, Atan2(Norm(q.xyz),q.w). Then reuse the norm computation for the normalization. Mar 20 '21 at 9:14

Assuming the matrices you have are pure rotation matrices (i.e. they do not contain scale information) and describe a rotation around an axis which is the same as the one you want to calculate the angular velocity for, then you can extract their angle $$\θ\$$ directly.

A rotation of angle $$\θ\$$ around axis $$\z\$$ in 3D space for example is described by the augmented matrix $$R=\begin{pmatrix} \cos\theta & -\sin\theta & 0 & 0\\ \sin\theta & \cos\theta & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{pmatrix}.$$ It is possible that you do not have the augmented matrix at hand in which case it will be a 3x3 matrix with the last row and last column trimmed off.

Hopefully the matrices you have are of a similar form (i.e. they describe a rotation around one of the basis vectors). Once you have extracted $$\θ_0\$$ and $$\θ_1\$$ you can find their difference $$\Δθ\$$ and from there you can use the angular velocity formula $$ω = \frac{Δθ}{Δt}$$ If the matrices you have are not around one of the basis vectors then this becomes more complicated but the logic stands. Check this for rotation matrices around an arbitrary axis

• yeah, the rotation matrix is around all three axis. It is the complete global rotation. Mar 18 '21 at 10:25
• A rotation is always around 1 axis, it's just that if it is an arbitrary axis then it can be broken down into 3 distinct rotations around each of the basis vectors. I think the url I linked could help you retrieve the angle $θ$ within the rotation matrix and use that. Mar 18 '21 at 10:28

The step rotation is $$\Q=R_1R_0^{-1}\$$ and the goal is to find its "logarithm" $$\X=α(\hat n\times)\$$ for the identity $$\Q=e^{X}\$$ where the unit vector $$\\hat n\$$ is the rotation axis, $$\(\hat n×)\$$ the skew-symmetric matrix realizing the cross product and $$\α\$$ the rotation angle.

If $$\Q\$$ is close to the identity, $$\X\$$ can be selected close to the zero matrix. Then it might be sufficient to approximate the matrix exponential by a Padé approximation and solve $$(2I+X)(2I-X)^{-1}=Q\iff X=2(Q-I)(Q+I)^{-1}=2I-4(Q+I)^{-1}$$ This will give a skew-symmetric result within the boundaries of the floating point accuracy. Read off $$X=α(\hat n\times)=\pmatrix{0&-z&y\\z&0&-x\\-y&x&0},~~ α=\|(x,y,z)\|, \hat n=\frac{(x,y,z)}{α}$$