Edit 19 Sept 2010. I strongly recommend using the eigenvector rather than trying to make a unit vector out of the off-diagonal terms. Find “edit” below.
Back in what has turned out to be my most popular post, axis and angle of rotation, I showed how to switch between the angle and axis of rotation and the rotation matrix.
Given the axis of rotation as a unit vector (a, b, c), and the angle of rotation, we construct the matrix N
and then we construct
and that’s a rotation. To be specific, it is an attitude matrix for a rotation of coordinate axes, about the axis (a, b, c) through the angle .
Conversely, given a rotation matrix, we find the axis of rotation as the eigenvector whose eigenvalue is 1, and then we use the formula for the trace of a rotation matrix
to get the angle .
The problem is that the eigenvector and the angle are each determined only to within a sign. There are 4 possible pairings, only 2 of which are correct. Given a candidate pair, we construct the rotation matrix from it, and compare it to the given rotation matrix.
I now know how to uniquely determine the two correct pairs. We need no longer test the possibilities.
The key was a sidebar comment on p. 74 in Kuipers’ “Quaternions and Rotation Sequences”. What he said was that given an arbitrary orthogonal 3×3 matrix , the axis of rotation was given by
where i, j, k are the usual unit coordinate-axis vectors. That is, the components are
I struggled with that for a while, silly me, before I realized that instead of trying to prove it for an arbitrary orthogonal 3×3, I should just use the magic formula for the rotation matrix in terms of the axis and angle.
Once I realized that, the obvious calculation is . Let
and we get
and then is
Just look at those off-diagonal terms! They each involve one of a, b, c, and a common factor of . Your first reaction might be that we don’t know , but that doesn’t matter. We can write the vector from the off-diagonal terms, and then we can get (a, b, c) just by making a unit vector from it. And once we have any one of a, b, c, we can get .
So Kuipers was right, and we can get the axis of rotation from the skew part of the rotation matrix. What he never said was that we could get the sine of the angle, too. And that, combined with the cosine, determines the quadrant.
Given a rotation matrix (representing a rotation of the coordinate axes), there are two ways to proceed. One is to enhance the existing algorithm:
- Get the rotation axis (a, b, c) as the unit eigenvector with eigenvalue 1,
- Get from the trace of the rotation matrix,
- Get from a (or b or c) and the appropriate off-diagonal term in .
The other way to go is to forget about the eigenvector calculation. Edit: No. Get the eigenvector. The following alternative can be numerically unsound. If is small, then the off-diagonal elements will be small, and making a unit vector out of them may not be accurate.
- Get the vector (a, b, c) from the off-diagonal terms,
- Get the vector (a, b, c) as a unit vector,
- Get from the trace of the rotation matrix.
- Get from a (or b or c) and the appropriate off-diagonal term in ,
In either case, having both the cosine and the sine of an angle determines it (mod , of course, but that’s not the issue). We know which quadrant we’re in.
Oh, we get one or the other of the two correct pairs of axis and angle by the two possible choices of the sign of (a, b, c). Pick one sign for the vector, get ; pick the other sign for the vector, get . These are the two correct pairs.
Let me emphasize that we work with two matrices: the rotation matrix A and the skew matrix .
Let’s do it both ways. Here’s the attitude matrix to get to Mars’ coordinate frame. (It rotates the solar system coordinate axes to Mars coordinate axes.)
I ask Mathematica® for the eigendecomposition, and select the eigenvector with eigenvalue 1. It is:
and that is a unit vector, whose components we call a, b, c, respectively.
We compute that the trace of the rotation matrix is 2.81888 and that’s , so we get the cosine of the angle:
Next we compute :
We know that the (2,3) location is , so we have…
Plug in a…
which gives us
Now, we want the angle whose tangent is , with the correct signs of the sine and cosine. I think that most software, nowadays, has a two-argument arctangent function which does exactly that, and Mathematica® certainly does. I get that the angle is (in radians)
And that is what we found originally: a valid answer was the eigenvector returned by Mathematica® and the negative angle.
Let’s try it the other way. From the general skew matrix
and the specific skew matrix
we construct the vector
which happens to be the negative of the eigenvector we had found. Make a unit vector out of it, which gives us values for (a, b, c):
That’s pedagogically very nice: we just happen to have made the other sign choice for the axis.
As before, we know that the (2,3) location is , so we have…
Plug in the value of a…
which gives us
which, as it must be, is the negative of the previous value for the sine.
The is still the same, so we find the new two-argument arctangent:
This time we got the positive angle. Of course, because we changed the sign on the rotation axis.
We have both of the two correct answers, and only the two correct answers. Having one, of course, we could write the other. I didn’t work the problem twice in order to get two solutions; I worked it twice to illustrate the two possible algorithms. I was very glad that they just happened to give the two different correct answers.
We no longer have to pick one of two unit vectors from column A and one of two angles from column B, and test pairs to see if they reproduce the rotation matrix.
As for the two algorithms, take your pick. My preference is to get the eigenvector, but Mathematica® makes that almost trivial for me to do. If it’s easier for you to get (a, b, c) from , go for it.