Path: news.mathworks.com!not-for-mail
From: "Bruno Luong" <brunoluong@yahoo.com>
Newsgroups: comp.soft-sys.matlab
Subject: Re: angle from rotation matrix
Date: Wed, 19 Dec 2007 22:14:58 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 98
Message-ID: <fkc552$si7$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> <fk6q3s$dns$1@fred.mathworks.com> <fk6vf6$9gt$1@fred.mathworks.com> <1dtem3965vls6097t13ven16m2keemmsqb@4ax.com> <fk7vji$9is$1@fred.mathworks.com> <fkbqtl$at5$1@fred.mathworks.com> <fkbs1d$25q$1@fred.mathworks.com>
Reply-To: "Bruno Luong" <brunoluong@yahoo.com>
NNTP-Posting-Host: webapp-05-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1198102498 29255 172.30.248.35 (19 Dec 2007 22:14:58 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Wed, 19 Dec 2007 22:14:58 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1184112
Xref: news.mathworks.com comp.soft-sys.matlab:443087


OK, I have made a small study to evaluate the methods. A
million random runs are done for three method:
- Direct method (D)
- MATLAB Eigen value method (E)
- James's quaternion (Q).

Next, I compute the median and mean of the L^2 errors
(100000 runs) in the rotation direction for three methods,
and here is the result:

Median results:
- D: 0.1241E-15
- E: 0.2077E-15
- Q: 0.1360E-15

Mean results:
- D: 0.6487088853611E-13
- E: 0.59115297358E-15
- Q: 0.274215765605692E-11

Conclusion:
- (E) is amazingly very stable, even in the difficult
situation where angle~0. It clearly outperforms the other
two methods n this case.
- (D) is less stable in difficult situation, but perform
overall 40% better than (E) for easy situation.
- (Q): Same remark as with (D). However it's (relatively)
poorly performed when angle close to 0.

Bellow is the code if you want to try it:

Bruno

%%%%%%%%%%%%%%%%%%%%%%%%

ntest=100000;

error=zeros(ntest,3);
for n=1:ntest
    theta=2*pi*(rand()-0.5);
    
    u=randn(3,1); u=u/norm(u);

    I=eye(3);
    R=[0 -u(3) u(2);
        u(3) 0 -u(1);
        -u(2) u(1) 0];
    T=u*u';
    A=cos(theta)*I + (1-cos(theta))*T + sin(theta)*R;

    d=diag(A);
    t=sum(d);
    uR=[A(3,2)-A(2,3);...
        A(1,3)-A(3,1); ...
        A(2,1)-A(1,2)];
    uR2=uR.^2;
    s2=sum(uR2);
    a=(3-t)+s2;
    if a>0
        theta=atan2(sqrt(s2),t-1);
        u2=(sum((A+A'-t*eye(3)).^2,2)+((4*d-2*t)+1))+uR2;
        absu=sqrt(u2.*(u2>0))/sqrt(a);
        signu=sign(uR);
        signu(signu==0)=1; % we don't want the nil case
        udir=absu.*signu;
        udir=udir/norm(udir);
    else
        warning('Rotation matrix close to indentity');
        udir=[1 0 0]';
        theta=0;
    end

    [P D]=eig(A);
    [dummy k]=min(diag(D)-1);
    ueig=P(:,k);
    
    [a v] = angleaxis(dc2quat(A));
    uqua=v';
    
    error(n,1)=min(norm(u-udir),norm(u+udir));
    error(n,2)=min(norm(u-ueig),norm(u+ueig));
    error(n,3)=min(norm(u-uqua),norm(u+uqua));
    
%     if (error(n,1)>4e-12)
%         return;
%     end
    
end

median(error,1)
mean(error,1)

semilogy(error);
legend('Direct','Eigen','Quaternion');

end