Path: news.mathworks.com!not-for-mail
From: "Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid>
Newsgroups: comp.soft-sys.matlab
Subject: Re: angle from rotation matrix
Date: Mon, 17 Dec 2007 07:47:11 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 42
Message-ID: <fk59hv$n5r$1@fred.mathworks.com>
References: <fjvbqk$gbt$1@fred.mathworks.com> <fk0cpd$g0v$1@fred.mathworks.com> <fk0iq9$k0j$1@fred.mathworks.com>
Reply-To: "Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid>
NNTP-Posting-Host: webapp-02-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1197877631 23739 172.30.248.37 (17 Dec 2007 07:47:11 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Mon, 17 Dec 2007 07:47:11 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1187260
Xref: news.mathworks.com comp.soft-sys.matlab:442709


"Bruno Luong" <brunoluong@yahoo.com> wrote in message <fk0iq9$k0j
$1@fred.mathworks.com>...
> One can use either eig() or direct calculation from A:
> 
> %%%% First method 
> 
> D=eig(A);
> [dummy k]=max(abs(D-1));
> theta=angle(D(k))
> 
> %%%%% Second method, returned angle in [0,pi]
> 
> c=(trace(A)-1);
> s=norm([A(3,2)-A(2,3);...
>         A(1,3)-A(3,1); ...
>         A(2,1)-A(1,2)]);
> theta=atan2(s,c) % better than using acos or asin
> 
> Bruno
---------
Hello again, Bruno.

  It has occurred to me that the second method of yours would not be robust 
with respect to accuracy of the rotation axis when the rotation angle theta is 
near pi.  In that case all the elements of your quantity u, computed from 
matrix A, would necessarily be very small differences of typically large 
quantities and subsequent division by norm(u) would yield excessive roundoff 
errors in the normalized u.  If norm(u) were zero in this case, the answer 
[1,0,0] would be completely wrong, in general.

  I have not been able to think of an alternative method of doing the problem 
in a robust manner, except by performing eigenvector analysis on it, as in 
your first method, and in mine.  In your first method it looks as though you 
would still need to compute u from the A matrix and take its dot product with 
the eigenvector whose eigenvalue is 1 in order to determine whether or not 
that dot product is equal to, and therefore has the same sign as, 2*sin(theta), 
where theta has been determined from one of the complex eigenvalues.  If 
not, the eigenvector should be reversed in direction, or the sign of theta 
reversed, so that the rotation axis unit vector would match the rotation angle.

Roger Stafford