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: Tue, 18 Dec 2007 04:38:39 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 74
Message-ID: <fk7isf$cgo$1@fred.mathworks.com>
References: <fjvbqk$gbt$1@fred.mathworks.com> <fk0cpd$g0v$1@fred.mathworks.com> <fk0iq9$k0j$1@fred.mathworks.com> <fk59hv$n5r$1@fred.mathworks.com> <2shcm3pkforks50k3jjkjo0bdkioir9jq2@4ax.com> <fk6eis$ra6$1@fred.mathworks.com> <fk6ief$3kq$1@fred.mathworks.com>
Reply-To: "Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid>
NNTP-Posting-Host: webapp-03-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1197952719 12824 172.30.248.38 (18 Dec 2007 04:38:39 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Tue, 18 Dec 2007 04:38:39 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1187260
Xref: news.mathworks.com comp.soft-sys.matlab:442859


"Bruno Luong" <b.luong@fogale.fr> wrote in message <fk6ief$3kq
$1@fred.mathworks.com>...
> Roger,
> 
> Thank you or pointing out the accuracy problem when
> theta~pi. Here is the change I made that work better. Though
> I don't think the issue when theta~0 could be overcome.
> 
> Bruno
> 
> c=(trace(A)-1);
> u=[A(3,2)-A(2,3);...
>    A(1,3)-A(3,1); ...
>    A(2,1)-A(1,2)];
> s=norm(u);
> theta=atan2(s,c)
> 
> if abs(2-c)>s % <- Use diagonal terms to estimate u
>     signu=sign(u);
>     signu(signu==0)=1; % we don't want the nil case 
>     u=sqrt((2*diag(A)-c)/(2-c)).*signu
> elseif s~=0 % <- use off-diagonal terms to estimate u
>     u=u/s
> else
>     warning('A close to identity');
>     u=[1 0 0]
> end
----------
  That's certainly an improvement over the former method, Bruno.  However, 
there are still some circumstances where it makes errors in the rotation axis 
direction much larger than with the eigenvector method.  This happens when 
the rotation angle is greater than pi/2 and when the rotation axis lies almost 
in the plane of a pair of coordinate axes.  The rotation angle, on the other 
hand, is determined just as accurately as that found with complex eigenvalue 
angles.

  For example, set

 p = .00000001234567;
 q = .65432101234567;
 r = sqrt(1-p^2-q^2); % = 0.75621690856720

for the three direction cosines of a unit vector along the rotation axis and

 w = .7*pi;

for the rotation angle.  The result with your new method is:

 u =
   0.00000000836198
   0.65432101234567
   0.75621690856720

where the first component is off in the ninth decimal place.  Doing it with 
matlab's 'eig' function gives a fully accurate value.

  This latter fact shows that the information is actually present in matrix A but 
that the new method is still not managing to fully extract it.  I suspect that a 
method which is as robust as the eigenvector method would have to make 
use of the sums of the matrix elements, A(3,2)+A(2,3), A(1,3)+A(3,1), and A
(2,1)+A(1,2), as well as their differences.

  If we go to the extreme case of a rotation angle exactly equal to pi, then A 
would be a symmetric matrix, and it is clear that we need to use more than 
just the three diagonal elements of A to determine the signs of the elements 
of u.  Of course, it is true in this case that the rotation axis can only be 
determined as a line, not the direction along the line, since rotating by pi in 
either the clockwise or counterclockwise directions about a line achieve the 
same result.  Nevertheless, of the eight possible sign combinations of 
elements of u, we ought to be able to reduce it to only two possibilities, but 
off-diagonal matrix elements would surely need to be taken into account.

Roger Stafford