angle and axis of rotation 2: the two correct answers

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

N = \left(\begin{array}{lll} 0 & c & -b \\ -c & 0 & a \\ b & -a & 0\end{array}\right)

and then we construct

A = I + N sin \theta + N^2 \left(1-cos \theta \right)

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 \theta\ .

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

trace[A] = 1 + 2\ cos \theta

to get the angle \theta\ .

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 A = a_{ij}\ , the axis of rotation was given by

i \left(a_{23} - a_{32}\right) +  j \left(a_{31} - a_{13}\right) +  k \left(a_{12} - a_{21}\right)

where i, j, k are the usual unit coordinate-axis vectors. That is, the components are

\left(\left(a_{23} - a_{32}\right),\   \left(a_{31} - a_{13}\right),\ \left(a_{12} - a_{21}\right)\right)\ .

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 A - A^T\ . Let

A = I + N sin \theta + N^2 \left(1-cos \theta \right)

and we get

A = \left(\begin{array}{lll} \left(-b^2-c^2\right) (1-\cos (\theta ))+1 & a b (1-\cos (\theta ))+c \sin (\theta ) & a c   (1-\cos (\theta ))-b \sin (\theta ) \\ a b (1-\cos (\theta ))-c \sin (\theta ) & \left(-a^2-c^2\right) (1-\cos (\theta ))+1 & b c   (1-\cos (\theta ))+a \sin (\theta ) \\ a c (1-\cos (\theta ))+b \sin (\theta ) & b c (1-\cos (\theta ))-a \sin (\theta ) &   \left(-a^2-b^2\right) (1-\cos (\theta ))+1\end{array}\right)

and then A - A^T\ is

A - A^T = \left(\begin{array}{lll} 0 & 2 c \sin (\theta ) & -2 b \sin (\theta ) \\ -2 c \sin (\theta ) & 0 & 2 a \sin (\theta ) \\ 2 b \sin (\theta ) & -2 a \sin (\theta ) & 0\end{array}\right)

Just look at those off-diagonal terms! They each involve one of a, b, c, and a common factor of sin \theta\ . Your first reaction might be that we don’t know sin \theta\ , but that doesn’t matter. We can write the vector 2\ sin \theta\ (a,\ b,\ c) 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 sin \theta\ .

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 cos \theta\ from the trace of the rotation matrix,
  • Get sin \theta\ from a (or b or c) and the appropriate off-diagonal term in A - A^T\ .

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 sin \theta\ 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 2\ sin \theta\ (a, b, c) from the off-diagonal terms,
  • Get the vector (a, b, c) as a unit vector,
  • Get cos \theta\ from the trace of the rotation matrix.
  • Get sin \theta\ from a (or b or c) and the appropriate off-diagonal term in A - A^T\ ,

In either case, having both the cosine and the sine of an angle determines it (mod 2\ \pi\ , 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 sin \theta\ ; pick the other sign for the vector, get -sin \theta\ . These are the two correct pairs.

Let me emphasize that we work with two matrices: the rotation matrix A and the skew matrix A - A^T\ .

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

A = \left(\begin{array}{lll} 0.90956 & -0.414415 & -0.0310051 \\ 0.414851 & 0.909845 & 0.00899314 \\ 0.0244829 & -0.0210423 & 0.999479\end{array}\right)

I ask Mathematica® for the eigendecomposition, and select the eigenvector with eigenvalue 1. It is:

\{-0.0361149,\ -0.0667194,\ 0.997118\}

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 1 + 2\ cos \theta\ , so we get the cosine of the angle:

\cos \left(\theta \right)\to 0.909442

Next we compute A - A^T\ :

\left(\begin{array}{lll} 0. & -0.829266 & -0.055488 \\ 0.829266 & 0. & 0.0300354 \\ 0.055488 & -0.0300354 & 0.\end{array}\right)

We know that the (2,3) location is 2\ a\ sin \theta\ , so we have…

0.0300354=2 a \sin (\theta )

Plug in a…

0.0300354=-0.0722298 \sin (\theta )

which gives us

\{\sin (\theta )\to -0.415831\}

Now, we want the angle whose tangent is \frac{\sin \theta}{\cos \theta}\ , 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)

\theta = -0.428857

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

\left(\begin{array}{lll} 0 & 2 c \sin (\theta ) & -2 b \sin (\theta ) \\ -2 c \sin (\theta ) & 0 & 2 a \sin (\theta ) \\ 2 b \sin (\theta ) & -2 a \sin (\theta ) & 0\end{array}\right)

and the specific skew matrix

\left(\begin{array}{lll} 0. & -0.829266 & -0.055488 \\ 0.829266 & 0. & 0.0300354 \\ 0.055488 & -0.0300354 & 0.\end{array}\right)

we construct the vector

2\ sin \theta\ (a,\ b,\ c) = \{0.0300354,\ 0.055488,\ -0.829266\}

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):

(a,\ b,\ c) = \{0.0361149,\ 0.0667194,\ -0.997118\}

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 2\ a\ sin \theta\ , so we have…

0.0300354=2 a\ \sin (\theta )

Plug in the value of a…

0.0300354=0.0722298 \sin (\theta )

which gives us

\{\sin (\theta )\to 0.415831\}

which, as it must be, is the negative of the previous value for the sine.

The cos \theta is still the same, so we find the new two-argument arctangent:

\theta = 0.428857

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 A - A^T\ , go for it.

Advertisements

8 Responses to “angle and axis of rotation 2: the two correct answers”

  1. Dave Cromley Says:

    I can’t tell you, Rip, how thankful I am for these articles. I worked for a couple of weeks at getting the eigenvector and the angle from the (rigid) rotation matrix, and was unable to do it. I even got into quaternions and thought that would do it, but maybe I didn’t use them correctly.

    What I want it for is to do rotations in thinBasic TBGL (Thin Basic Graphics Language). TBGL is based on openGL — I don’t know any more than that. At any rate, the command to rotate is

    TBGL_Rotate angle, evx, evy, evz

    where evx and evy and evz are any multiple of the eigenvector, and the angle is in degrees.

    My methods were getting into the bizarre “gimbal lock” and was unstable. So your articles were just what I wanted. My code from them work great.

    Here is my (thinBasic) code to get from the rotation matrix to the angle and evx, evy, and evz:

    EVX = m23-m32: EVY = m31-m13: EVZ = m12-m21
    EVLen = sqr(EVX^2 + EVY^2 + EVZ^2)
    Cosa = (m11 + m22 + m33 – 1.) / 2. ‘ (trace-1)/2
    Sina = EVLen / 2.
    TBGL_Rotate atan2(Sina, Cosa), EVX, EVY, EVZ

    Two points:
    1) The Eigenvector is not made to be a unit vector.
    2) The equation for Sina = EVX /2a is inadequate because “a” could be 0. A good equation for Sina is miraculously Sina = EVLen / 2.

    I’m going to follow up with a URL to my postings on community.thinbasic.com. I think I can get a clean executable of my demo program. By “clean”, I mean it stands alone and doesn’t leave anything anywhere. If so, I’ll have a URL for that also.

    Thanks again. I’m so excited, I can hardly sleep. Dave C

  2. rip Says:

    Hi Dave,
    I am delighted to hear that some of my posts were a significant help. I want to take a closer look at your code, since you were so kind as to post it.

    Please let me know if you have any questions or difficulties.

    Thanks again.

  3. Chris Says:

    Awesome set of posts. You helped me solve a major problem when (20-years out of college) I have forgotten more math than I remember. Great way of explaining the fundamentals of Axis/Angle rotation.

  4. rip Says:

    Thanks, Chris.

    I’m glad it made sense, and that it helped with something. I appreciate the feedback.

  5. Daniel C Says:

    Hello, I’m trying to solve the following problem: I have an orthonormal coordinate system with axes known
    u = ((x1, y1, z1), (x2, y2, z2), (x3, y3, z3))
    and another orthonormal with known axes
    u’ = ((x’1, y’1, z’1), (x’2, y’2, Z’2), (x’3, y’3, z’3))
    and need finding the axis of rotation and put in the coordinate system
    (1,0,0), (0,1,0), (0,0,1)
    and find the angle of rotation, still do not understand how, appreciate your help

    (sorry for my bad English)

    • rip Says:

      Hi Daniel,

      If I understand correctly, you want to to find the axis-angle representation of the rotation matrix which relates u to u’. Well, if you write

      u = u’ R

      then

      R = u u'^{-1} = u u'{^T}

      and this post showed how to get the axis-angle of any R. Note that I often use ‘ to denote a transpose, but in your notation, u and u’ are distinct rotation matrices.

  6. Daniel C Says:

    I found something interesting on page 4 of this pdf

    http://www.seas.upenn.edu/~meam520/notes02/EulerChasles4.pdf

    • rip Says:

      Hi Daniel,

      It appears to me that the PDF matches my description of how to find the sine of the angle: compute A – A’ and solve for one of the sines.

      rip


Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: