MATLAB: Something wrong with quat2angle() function

Question

Hi MATLAB community,

I was wondering why the following is happening:

As you can see, conversion for the first case is right (quaternion indeed represents this rotation) but the second rotation is wrong. One would expect angleX2 to be 90 degrees (pi/2) as well, next to angleY2 being 90 degrees.

Any ideas?

Best,

Nico

Answer 1

I agree that this looks like a bug, and you should submit a bug report to TMW for this. I ran it in R2017a and R2020a and both versions have the bug. I first started with this post:

https://www.mathworks.com/matlabcentral/answers/244488-quat2angle-function?s_tid=srchtitle

I modified that code to use your 'ZYX' and 'XYZ' rotation sequences for q = [0.5 0.5 0.5 0.5], and I compared that result to the results of a conversion to direction cosine matrix and a quaternion rotation multiply. The 'ZYX' matched but the 'XYZ' didn't. Hence I conclude that it is a bug in the quat2angle( ) function. If you perturb the quaternion slightly (use a small value for f, e.g. 1e-10) it gets the correct result. But it can't handle that exact quaternion properly. The code I used for this is:

disp('------------------------------')
disp('quat2angle( ) test')
% sample data
f = 1e-10; % Use 0 here to demonstrate quat2angle bug
q = [0.5 0.5 0.5 0.5] + (rand(1,4)*2-1)*f;
q = q / norm(q)
vi = [1 2 3]';
% quat2angle ZYX
disp('------------------------------')
disp('dcm = quat2dcm(q)')
dcm = quat2dcm(q)
[a3,a2,a1] = quat2angle(q,'ZYX');
c = cos(a3);
s = sin(a3);
r3 = [ c  s  0;
      -s  c  0;
       0  0  1];
c = cos(a2);
s = sin(a2);
r2 = [ c  0 -s;
       0  1  0;
       s  0  c];
c = cos(a1);
s = sin(a1);
r1 = [ 1  0  0;
       0  c  s;
       0 -s  c];
disp('dce = result of individual rotation matrix multiplies, ZYX')
dce = r1 * r2 * r3
% compare rotation methods

vb = dcm * vi;
ve = dce * vi;
vq = quatmultiply(quatconj(q),[0 vi']);
vq = quatmultiply(vq,q);
vr = quatrotate(q,vi');
disp('------------------------------')
disp('vi = input inertial vector')
disp(vi')
disp('vb = dcm * vi')
disp(vb')
disp('ve = dce * vi')
disp(ve')
disp('vq = conj(q) * vi * q')
disp(vq(2:4))
disp('vr = quatrotate(q,vi'')')
disp(vr)
disp('q')
disp(q);
disp('angle2quat(angles) for ZYX')
disp(angle2quat(a3,a2,a1,'ZYX'))
disp('------------------------------')
% quat2angle XYZ
disp('------------------------------')
[a3,a2,a1] = quat2angle(q,'XYZ');
c = cos(a3);
s = sin(a3);
r3 = [ 1  0  0;
       0  c  s;
       0 -s  c];
c = cos(a2);
s = sin(a2);
r2 = [ c  0 -s;
       0  1  0;
       s  0  c];
c = cos(a1);
s = sin(a1);
r1 = [ c  s  0;
      -s  c  0;
       0  0  1];
disp('dce = result of individual rotation matrix multiplies, XYZ')
dce = r1 * r2 * r3
% compare rotation methods
vb = dcm * vi;
ve = dce * vi;
vq = quatmultiply(quatconj(q),[0 vi']);
vq = quatmultiply(vq,q);
vr = quatrotate(q,vi');
disp('------------------------------')
disp('vi = input inertial vector')
disp(vi')
disp('vb = dcm * vi')
disp(vb')
disp('ve = dce * vi')
disp(ve')
disp('vq = conj(q) * vi * q')
disp(vq(2:4))
disp('vr = quatrotate(q,vi'')')
disp(vr)
disp('q')
disp(q);
disp('angle2quat(angles) for XYZ')
disp(angle2quat(a3,a2,a1,'XYZ'))
disp('------------------------------')

And the results I get for f=0 are:

>> quat2angle_test
------------------------------
quat2angle( ) test
f =
     0
q =
    0.5000    0.5000    0.5000    0.5000
------------------------------
dcm = quat2dcm(q)
dcm =
     0     1     0
     0     0     1
     1     0     0
dce = result of individual rotation matrix multiplies, ZYX
dce =
    0.0000    1.0000         0
   -0.0000    0.0000    1.0000
    1.0000   -0.0000    0.0000
------------------------------
vi = input inertial vector
     1     2     3
vb = dcm * vi
     2     3     1
ve = dce * vi
     2     3     1
vq = conj(q) * vi * q
     2     3     1
vr = quatrotate(q,vi')
     2     3     1
q
    0.5000    0.5000    0.5000    0.5000
angle2quat(angles) for ZYX
    0.5000    0.5000    0.5000    0.5000
------------------------------
------------------------------
dce = result of individual rotation matrix multiplies, XYZ
dce =
    0.0000         0   -1.0000
         0    1.0000         0
    1.0000         0    0.0000
------------------------------
vi = input inertial vector
     1     2     3
vb = dcm * vi
     2     3     1
ve = dce * vi
   -3.0000    2.0000    1.0000
vq = conj(q) * vi * q
     2     3     1
vr = quatrotate(q,vi')
     2     3     1
q
    0.5000    0.5000    0.5000    0.5000
angle2quat(angles) for XYZ
    0.7071         0    0.7071         0
------------------------------

You can see the discrepancy in the dcm produced and the different quaternion produced by angle2quat when using the 'XYZ' rotation sequence. There may be other quaternion and sequence combinations that produce erroneous errors for edge cases also.

When you perturb the q slightly (small non-zero f = 1e-10), things work correctly:

>> quat2angle_test
------------------------------
quat2angle( ) test
f =
   1.0000e-10
q =
    0.5000    0.5000    0.5000    0.5000
------------------------------
dcm = quat2dcm(q)
dcm =
    0.0000    1.0000   -0.0000
   -0.0000    0.0000    1.0000
    1.0000   -0.0000    0.0000
dce = result of individual rotation matrix multiplies, ZYX
dce =
    0.0000    1.0000   -0.0000
   -0.0000    0.0000    1.0000
    1.0000   -0.0000    0.0000
------------------------------
vi = input inertial vector
     1     2     3
vb = dcm * vi
    2.0000    3.0000    1.0000
ve = dce * vi
    2.0000    3.0000    1.0000
vq = conj(q) * vi * q
    2.0000    3.0000    1.0000
vr = quatrotate(q,vi')
    2.0000    3.0000    1.0000
q
    0.5000    0.5000    0.5000    0.5000
angle2quat(angles) for ZYX
    0.5000    0.5000    0.5000    0.5000
------------------------------
------------------------------
dce = result of individual rotation matrix multiplies, XYZ
dce =
    0.0000    1.0000         0
   -0.0000         0    1.0000
    1.0000   -0.0000    0.0000
------------------------------
vi = input inertial vector
     1     2     3
vb = dcm * vi
    2.0000    3.0000    1.0000
ve = dce * vi
    2.0000    3.0000    1.0000
vq = conj(q) * vi * q
    2.0000    3.0000    1.0000
vr = quatrotate(q,vi')
    2.0000    3.0000    1.0000
q
    0.5000    0.5000    0.5000    0.5000
angle2quat(angles) for XYZ
    0.5000    0.5000    0.5000    0.5000
------------------------------