From: Pertti Lounesto <lounesto@torstai.hit.fi>
Newsgroups: sci.math
Subject: Re: rotation matrix
Date: 28 Aug 1997 13:39:37 +0300

poseyj@bek.com (Jeff Posey) writes:
> Given a 3x3 rotation matrix, how can you extract a corresponding 3D
> rotation axis vector and the angle of rotation about this vector?

Let U be the rotation matrix.  Then A=(U-I)/(U+I) is antisymmetric,
and x=(A23,-A13,A12) is the rotation axis.  The angle of rotation
is arccos((trace(U)-1)/2).
-- 
Pertti Lounesto          http://www.hit.fi/~lounesto

==============================================================================

From: glynnec@ix.netcom.com(Glynne Casteel)
Newsgroups: sci.math
Subject: Re: rotation matrix
Date: 30 Aug 1997 03:04:47 GMT

In <7ld8my4mmu.fsf@torstai.hit.fi> Pertti Lounesto
<lounesto@torstai.hit.fi> writes: 
>
>poseyj@bek.com (Jeff Posey) writes:
>> Given a 3x3 rotation matrix, how can you extract a corresponding 3D
>> rotation axis vector and the angle of rotation about this vector?
>
>Let U be the rotation matrix.  Then A=(U-I)/(U+I) is antisymmetric,
>and x=(A23,-A13,A12) is the rotation axis.  The angle of rotation
>is arccos((trace(U)-1)/2).
>

I like the expression for the angle of rotation, but you don't need to
do a matrix inversion to obtain the axial vector.  Here is a much
cheaper way to obtain the vector.
-----------------------------------------------------

In terms of a unit vector, n, and an angle of rotation, b, a rotation
tensor is given by:

    R = I + N sin(b) + N^2 (1-cos(b))
        where N = skew(n) = n x I
        and I = identity tensor
    
A little vector algebra will show that   N^2 = nn - I.

The trace of R yields the following:
    tr(R) = tr(I) + tr(N) sin(b) + tr(N^2) (1-cos(b))
          = 3   +   0    +  tr(nn-I)(1-cos(b))
          = 3  +  tr(nn)(1-cos(b))  -  tr(I)(1-cos(b))
          = 3  +  (1-cos(b))   - 3(1-cos(b))
          = 1 + 2 cos(b)

        where I have made liberal use of the fact that
            tr(I) = 3
            tr(N) = 0  (since it's skew-symmetric)
            tr(nn) = n.n = 1  (since it's a unit vector)

We can rearrange this last result to yield an expression for cos(b) in
terms of tr(R):
    cos(b) = (1/2) (tr(R)-1)

Since R is a rotation matrix, the transpose (or inverse) is the result
of rotating about the same n-vector but in the opposite direction (-b).
This reverses the sign of sin(b), but nothing else:

    R^t = I - N sin(b) + N^2 (1-cos(b))

We can use this to eliminate the N-term:

    (1/2) (R+R^t) - I = N^2 (1-cos(b))
                      = nn(1-cos(b)) - I(1-cos(b))

    (1/2) (R+R^t) - I cos(b) = nn(1-cos(b))

We're almost there.  Just rearrange this last equation for nn:

         (1/2)(R+R^t) - (1/2)(tr(R)-1)I
    nn = ------------------------------
             (1/2)(3-tr(R))


         (R+R^t) + (1-tr(R))I
    nn = ---------------------
             (3-tr(R))


If we call this last matrix M, then we can pick off the components of
the n-vector from the diagonal of M:

    n_1 = sqrt(M_11)
    n_2 = sqrt(M_22)
    n_3 = sqrt(M_33)



Sure, it was a lot of work to derive, but calculating the trace and
transpose of R is trivial compared to calculating the inverse of (R+I).

~Glynne