Something wrong with quat2angle() function?

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

[James Tursa]

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
------------------------------

[James Tursa]

Maybe this is related to another bug report:

https://www.mathworks.com/support/bugreports/1551143?s_tid=answers_rc2-1_p4_Topic

[Paul]

Any idea how to interpret quat2angle when the norm of the input quaternion is very different from unity, as in the examples on the quat2angle doc page?

[James Tursa]

This function, along with many other quaternion related functions in the various toolboxes, normalizes the input quaternions when a unit quaternion is expected. So you will get the same answer with any multiple of the unit quaternion. E.g.,

>> [yaw, pitch, roll] = quat2angle([1 0 1 0])
yaw =
     0
pitch =
    1.5708
roll =
     0
>> [yaw, pitch, roll] = quat2angle([5 0 5 0])
yaw =
     0
pitch =
    1.5708
roll =
     0
>> [yaw, pitch, roll] = quat2angle([sqrt(2)/2 0 sqrt(2)/2 0])
yaw =
     0
pitch =
    1.5708
roll =
     0

[Paul]

@James Tursa

That's about what I expected. Too bad that doc page doesn't explicitly state that the input will be normalized regardless of what its norm actually is. Maybe there's an overarching statement about that somehwere else in the toolbox doc. Still, it looks a bit odd to show examples with an input quaternion that doesn't even come close to having a unit norm. I guess one example of such might be useful if the intent is to highlight that feature of the function. I always though that the normalization was to smooth out small rounding errors that can accumulate in forming a quaternion, not to just brush aside gross errors. But that's their choice I guess.

@Nico de Korte

Did you find any other similar error cases? I'm pretty sure there are more. Did you file a bug report? If so, what was the response? If it actually turns into a formal bug report, can you post the link here?

[James Tursa]

@Paul

Yes I agree that the doc is lacking in the quaternion normalization assumptions, and many examples are misleading and sometimes just plain wrong.

[Paul]

@James Tursa

Can you post a link to an example (or two or ...) that's just plain wrong? I'm very curious about that, particularly if in the Aerospace Toolbox.