Reversing RotateAxisAngle back to angles

Question

I'm trying to figure out how to reverse RotateAxisAngle to get back rotations around these arbitrary axes (or equivalent rotations that yield same net rotation, doesn't have to be identical). Does anyone know how to do it? I'm using MathGeoLib, but I don't see an opposite way, to return back the angles about axes, when all you have is the matrix.

Here's the forward direction code (RotateAxisAngle is from MathGeoLib):

float4x4 matrixRotation = float4x4::RotateAxisAngle(axisX, ToRadian(rotation.x));
matrixRotation = matrixRotation * float4x4::RotateAxisAngle(axisY, ToRadian(rotation.y));
matrixRotation = matrixRotation * float4x4::RotateAxisAngle(axisZ, ToRadian(rotation.z));

Now I want to get back to the degrees, about these arbitrary axes, in the same order (well, pull off Z, then Y, then X), so if I did it again, forward direction, would yield the same net rotations.

Here's the sample/matrix corresponding to that set of rotations I posted above, if it helps, on reversing back to it:

axisX:
x   0.80878228  float
y   -0.58810818 float
z   0.00000000  float
Rot about that axis:
30.000000   float

axisY:
x   0.58811820  float
y   0.80877501  float
z   0.00000000  float
Rot about that axis:
60.000000   float

axisZ:
x   0.00000000  float
y   0.00000000  float
z   1.0000000   float
Rot about that axis:
40.000000   float

That forms this matrix, which is stored to a file, and need to retrieve the rotations about above axes (without any info about rotations originally used)

[4][4]
[0x0]   0.65342271  float
[0x1]   -0.51652151 float
[0x2]   0.55339342  float
[0x3]   0.00000000  float
[0x0]   0.69324547  float
[0x1]   0.11467478  float
[0x2]   -0.71151978 float
[0x3]   0.00000000  float
[0x0]   0.30405501  float
[0x1]   0.84856069  float
[0x2]   0.43300733  float
[0x3]   0.00000000  float
[0x0]   0.00000000  float
[0x1]   0.00000000  float
[0x2]   0.00000000  float
[0x3]   1.0000000   float

Answer 1

OK, I'm going to take another stab at this. My first answer was for XYZ order of rotations. This answer is for ZYX order, now that I know more about how MathGeoLib works.

MathGeoLib represents position vectors as column vectors v = [x y z 1]^T where ^T is the transpose operator which flips rows to columns (and vice versa). Rotation matrices pre-multiply column vectors. So if we have a matrix Rx(s) representing a rotation around the x-axis by s degrees, then a rotation Ry(t) representing a rotation about the y-axis by t degrees, then rotation Rz(u) representing a rotation about the z-axis by u degrees, and we combine them and multiply with v as Rx(s) Ry(t) Rz(u) v, we're actually applying the z rotation first. But we can still work out the angles from the combined matrix, it's just that the formulas are going to be different from the more common XYZ order.

We have the upper left blocks of the rotation matrices as follows. (The fourth row and column are all 0s except for the diagonal element which is 1; that never changes in the calculations that follow so we can safely ignore.) MathGeoLib appears to use left-handed coordinates, so the rotation matrices are:

        [1      0      0]          [ cos t  0  sin t]          [ cos u -sin u  0]
Rx(s) = [0  cos s -sin s], Ry(t) = [     0  1      0], Rz(u) = [ sin u  cos u  0]
        [0  sin s  cos s]          [-sin t  0  cos t]          [     0      0  1]

(Note the position of the - sign in Ry(t); it's there because we think of the coordinates in cyclic order. Rx(s) rotates y and z; Ry(t) rotates z and x; Rz(u) rotates x and y. Since Ry(t) rotates z and x not in alphabetical order but in cyclic order, the rotation is opposite in direction from what you expect for alphabetical order.

Now we multiply the matrices in the correct order. Rx(s) Ry(t) is

[1      0      0][ cos t  0  sin t]   [       cos t      0        sin t] 
[0  cos s -sin s][     0  1      0] = [ sin s*sin t  cos s -sin s*cos t]
[0  sin s  cos s][-sin t  0  cos t]   [-cos s*sin t  sin s  cos s*cos t]

The product of that with Rz(u) is

[       cos t      0        sin t][ cos u -sin u  0] 
[ sin s*sin t  cos s -sin s*cos t][ sin u  cos u  0] =
[-cos s*sin t  sin s  cos s*cos t][     0      0  1]

[                   cos t*cos u                   -cos t*sin u        sin t]
[ sin s*sin t*cos u+cos s*sin u -sin s*sin t*sin u+cos s*cos u -sin s*cos t]
[-cos s*sin t*cos u+sin s*sin u  cos s*sin t*sin u+sin s*cos u  cos s*cos t]

So we can figure out the angles as follows:

tan s = -(-sin s * cos t)/(cos s * cos t) = M23/M33 => s = -arctan2(M23,M33)
sin t = M13 => t = arcsin(M13)
tan u = -(-cos t * sin u)/(cos t * cos u) = M12/M11 => u = -arctan2(M12,M11)

If we are to implement those calculations, we need understand how the matrix is indexed in MathGeoLib. The indexing is row major, just like the mathematical notation, but indexing starts at 0 (computer style) rather than at 1 (math style) so the C++ formulas you want are

s = -atan2(M[1][2],M[2][2]);
t = asin(M[0][2]);
u = -atan2(M[0][1],M[0][0]);

The angles are returned in radians so will need to be converted to degrees if desired. You should test that result in the case when axes for the rotations Z, Y, and X are in standard position (001), (010), and (100).

If we are to reverse a rotation about non-standard axes, as in your example, the problem becomes more difficult. However, I think it can be done by a "change of coordinates". So if our rotation mystery matrix is matrixRotation, I believe you can just form the "conjugate" matrix

M = coordinateChangeMatrix*matrixRotation*coordinateChangeMatrix^{-1}

and then use the above formulas. Here coordinateChangeMatrix would be the matrix

[Xaxis0 Xaxis1 Xaxis2 0]
[Yaxis0 Yaxis1 Yaxis2 0]
[Zaxis0 Zaxis1 Zaxis2 0]
[     0      0      0 1]

where the rotation X-axis is (Xaxis0,Xaxis1,Xaxis2). In your example those numbers would be (0.808...,-0.588...,0). You should make sure that the rotation matrix is orthonormal, i.e., the dot product of Xaxis with itself is 1, the dot product of Xaxis with another axis is 0, and the same for any other axis. If the coordinate change matrix is not orthonormal, the calculation might still work, but I don't know for sure.

The inverse of the coordinate change matrix can be calculated using float4x4::inverseOrthonormal or if it's not orthonormal you can use float4x4::inverse but as I mentioned I don't know how well that will work.

Answer 2

If you just want a rotation that reverses the rotation you've got in one step, you can invert the rotation matrix. float4x4::InverseOrthonormal should work, and it's fast and accurate. float4x4::Inverse will also work but it's slower and less accurate.

If you really want to recover the angles, it goes something like this. (There are numerous different conventions, even for X-Y-Z; I think this one matches, but you may have to take the transpose of the matrix or make some other modification. If this doesn't work I can suggest alternatives.) First we follow the Wikipedia article for a description of the Euler Angles to Matrix conversion. In the resulting matrix, we have

A11 = cos theta cos psi
A21 = -cos theta sin psi
A31 = sin theta
A32 = -sin phi cos theta
A33 = cos phi cos theta

where phi is the rotation around the x-axis, theta is the rotation around the y-axis, and psi is rotation around the z-axis. To recover the angles, we do

phi = -arctan2(A32,A33)
theta = arcsin(A31)
psi = -arctan2(A21,A11)

The angles may not match the original angles exactly, but the rotations should match. arctan2 is the two-argument form of the arctan function, which takes into account the quadrant of the point represented by the argument, and deals correctly with 90 degree angles.

Given the way your rotations are represented, I think you may have to use the transpose instead. That's easy: you just swap the indices in the above formulas:

phi = -arctan2(A23,A33)
theta = arcsin(A13)
psi = -arctan2(A12,A11)

If neither of those work, I could take a closer look at the MathGeoLib library and figure out what they're doing.

Update

I neglected to take into account information about the axes of rotation in my previous answer. Now I think I have a plan for dealing with them.

The idea is to "change coordinates" then do the operations as above in the new coordinates. I'm a little hazy on the details, so the process is a little "alchemical" at the moment. The OP should try various combinations of my suggestions and see if any of them work ... there aren't too many (just 4 ... for the time being).

The idea is to form a coordinate change matrix using the coordinates of the rotation axes. We do it like this:

axisX: 0.80878228 -0.58810818 0.00000000 0.00000000
axisY: 0.58811820  0.80877501 0.00000000 0.00000000
axisZ: 0.00000000  0.00000000 1.0000000  0.00000000
and..: 0.00000000  0.00000000 0.00000000 1.0000000

I have just taken the three 3-vectors axisX, axisY, axisZ, padded them with 0 at the end, and added the row [0 0 0 1] at the bottom.

I also need the inverse of that matrix. Since the coordinate system is an orthonormal frame, the inverse is the transpose. You can use the InverseOrthonormal function in the library; all it does is form the transpose.

Now take your mystery matrix, and pre-multiply it by the coordinate change matrix, and post-multiply it by the inverse of the coordinate change matrix. Then apply one of the two calculations above using inverse trig functions. Crossing my fingers, I think that's it ...

If that doesn't work, then pre-multiply the mystery matrix by the inverse of the coordinate change matrix, and post-multiply by the coordinate change matrix. Then apply one or the other of the sets of trig formulas.

Does it work?