Why do the eul2quat and quaternion functions differ on the calculated results?

The following discussion is all about external rotation, which corresponds to the point rotation in matlab：

I use the built-in functions "eul2quat" and "quaternion" for quaternion conversion of the same set of Euler angles data, but the results are slightly different and the signs are also different? How can I explain this, thanks in advance!

format long
eulars = [-0.001180313	-0.1101197	-3.122032;
-0.001413606	-0.1086281	-3.119092;
-0.001928661	-0.1062565	-3.116415];
a = eul2quat(eulars,'XYZ') % first X, then Y, final Z ?

a = 3×4

0.009797824155183 0.055023629743008 -0.001127456560182 -0.998436343315603 0.011271864801883 0.054275960421124 -0.001316453379056 -0.998461482549404 0.012621931438076 0.053086905293953 -0.001631383177894 -0.998508791109027

b = quaternion(eulars,'euler','XYZ','point') % first X, then Y, final Z ?

b = 3×1 quaternion array
      0.00973287223990636 -   0.0550351559117731i +  5.1011347229122e-05j -    0.998436978586816k
       0.0111951287419489 -   0.0542918401622062i +  9.4978608694494e-05j -     0.99846234589158k
       0.0125195213825984 -   0.0531111499712373i + 0.000294403625928257j -    0.998510080399586k

But I changed the order of eul2quat to "ZYX" and it seems that the differences in values are equal, but the output quaternions are in the opposite order.

a_OK = eul2quat(eulars,'ZYX') % first X, then Y, final Z , but the output order is a bit counter-intuitive, in the form of [w,z,y,x]

a_OK = 3×4

0.009732872239906 -0.998436978586816 0.000051011347229 -0.055035155911773 0.011195128741949 -0.998462345891580 0.000094978608694 -0.054291840162206 0.012519521382598 -0.998510080399586 0.000294403625928 -0.053111149971237

Answer 1

[Brian Fanous]

Quick takeaways :

• Both eul2quat and the quaternion class use intrinsic rotations.
• eul2quat([a b c], 'XYZ') is equivalent to quaternion([a b c], 'euler', 'XYZ', 'frame')

Some more details about the quaternion class:

Frame Rotations

When you type quaternion([a b c], 'euler', 'XYZ', 'frame') that means find a quaternion q, such that

for a 3-by-1 vector v and rotation matrices R. (Obviously v on the right hand side is a pure quaternion with imaginary parts equal to v). The right hand side is quaternion frame rotation - hence frame in the quaternion() call. Also note that is the rotation matrix that rotates the frame of reference clockwise around the x axis by a. A new reference frame is achieved after each subsequent rotation because (as noted above), rotations are intrinsic.

Point Rotations

When you type quaternion([a b c], 'euler', 'XYZ', 'point') that means find a quaternion q, such that

for a 3-by-1 vector v and rotation matrices T. (Again, v on the right hand side is a pure quaternion with imaginary parts equal to v). The right hand side is quaternion point rotation - hence point in the quaternion() call. Also note that is the rotation matrix that rotates a point clockwise around the x axis by a. Note here that

[Paul]

@Brian Fanous

Does any Matlab function in any toolbox to your knowledge use extrinsic rotations?

[Brian Fanous]

I don't believe so. But you can convert an extrinsic sequence to an intrinsic one.

Answer 2

[Sam Chak]

Hi @cui

I remember that previously we have talked about that Euler angles and quaternion are not unique. You need to check the rotation matrix. See example below:

eul = [-0.001180313	-0.1101197	-3.122032]
qXYZ = angle2quat(eul(1), eul(2), eul(3), 'XYZ')

dcm1 = angle2dcm(eul(1), eul(2), eul(3), 'XYZ')
dcm2 = quat2dcm(qXYZ)

dcm1 - dcm2

[cui,xingxing]

Right, the previous discussion was that the Euler angles are not unique, but the rotation matrix is. This time I think the quaternion should also be unique, thank you for your answer,@Sam Chak， what is the function "angle2dcm" in your code? Angle to quaternion?

[Sam Chak]

angle2dcm converts the rotation angles to direction cosine matrix (aka rotation matrix). https://www.mathworks.com/help/aerotbx/ug/angle2dcm.html

angle2quat converts rotation angles to quaternion, and it works similarly like eul2quat.

Answer 3

[Paul]

Hi cui,

I think that both eul2quat and quaternion('frame') are using, right-handed, intrinsic rotations and the sequence of rotations is in the order of the sequence argument to those functions. Furthermore, eul2quat is doing a frame rotation.

Define three Euler angles

rng(100);
theta = rand(1,3);

Quaternion using eul2quat, intrinsic angles are: first around Z, second around Y', third around X''.

q1 = eul2quat(theta,'ZYX')

q1 = 1×4

0.9404 0.1646 0.1866 0.2317

Quaternion based on frame rotation

q2 = quaternion(theta,'euler','ZYX','frame');
compact(q2)

ans = 1×4

0.9404 0.1646 0.1866 0.2317

So those two are the same.

The confusion, at least for me, comes in when generating the associatied direction cosine matrices (DCMs). I'm going to generate a DCM that I understand as a frame rotation for the intrinsic ZYX sequence

Zrot = [cos(theta(1)) sin(theta(1)) 0; -sin(theta(1)) cos(theta(1)) 0;0 0 1];
Yrot = [cos(theta(2)) 0 -sin(theta(2)); 0 1 0; sin(theta(2)) 0 cos(theta(2))];
Xrot = [1 0 0;0 cos(theta(3)) sin(theta(3)); 0 -sin(theta(3)) cos(theta(3))];
DCMp = Xrot*Yrot*Zrot;

Given a vector v resolved in coordinate frame xyz expressed in a column, DCMp resolves that same vector into frame x''y''z'' via post-multiplication by v

v = rand(3,1);
vp = DCMp*v;

That works the same using q2

DCM2 = rotmat(q2,'frame');
v2 = DCM2*v

v2 = 3×1

0.6642 -0.2641 0.4664

However, the doc page for quat2rotm says that DCM1 is used by pre-multiplying by transpose(v), which is actually called "post-multiply format" here! (side note: I very much dislike this idea of pre-multiplication by a row vector). So we'll do that, and then transpose the result to get it back to column vector for comparison

DCM1 = quat2rotm(q1);
v1 = (v.' * DCM1).';

Compare the results, they are very, very close.

format long e
[vp - v1, vp - v2]

ans = 3×2

1.0e+00 * 1.110223024625157e-16 0 5.551115123125783e-17 -1.110223024625157e-16 -1.110223024625157e-16 -5.551115123125783e-17

So, we see that angle inputs to eul2quat define a right-handed, frame rotation, but we have to be careful about the pre- and post-multiplication convention.

[cui,xingxing]

Hi, @Paul

Thank you for doing some experiments, but you expressed I personally think only part of the easy to make people accept, the other part should still be the document is not uniformly clear presentation caused. I will show you the following example.

The premise is still a right-handed Cartesian coordinate system, using the common geometric transformation convention that the point coordinates are in the right-hand side of the rotation matrix in the form of column vectors. According to your first experiment, q1 and q2 agree that indeed intrinsic rotations and the sequence of rotations is in the order of the sequence argument to those functions. but from the user's point of view The results for eul2quat are indeed confusing from the user's point of view, because the eul2* family of functions should have a uniform conversion convention according to a general understanding, such as the following to continue your experiment.

rng(100);
theta = rand(1,3); % your origin data
q1 = eul2quat(theta,'ZYX')

q1 = 1×4

0.9404 0.1646 0.1866 0.2317

%% my experiment
R1 = eul2rotm(theta,'ZYX');
R2 = quat2rotm(q1)

R2 = 3×3

0.8230 -0.3743 0.4273 0.4971 0.8385 -0.2231 -0.2748 0.3960 0.8762

We know that the common rotation matrix convention has the following form.

thetaR = rad2deg(theta);
myR = rotz(thetaR(1))*roty(thetaR(2))*rotx(thetaR(3));% common geometric transformation convention
isEqu1 = all((R1-R2)< 1e-10,'all')

isEqu1 = logical
   1

isEqu2 = all((R1-myR)<1e-10,'all')

isEqu2 = logical
   1

This experiment is enough to show that the rotation matrix obtained from the conversion of the eul2rotm function is the result of our conventional understanding of extrinsic rotation/point rotation, which is fully consistent with our generic conventions, but the eul2quat function your experimental conclusion is intrinsic rotation/frame rotation, which is not uniform and not easily accepted by common people.

Zrot = [cos(theta(1)) sin(theta(1)) 0; -sin(theta(1)) cos(theta(1)) 0;0 0 1];
Yrot = [cos(theta(2)) 0 -sin(theta(2)); 0 1 0; sin(theta(2)) 0 cos(theta(2))];
Xrot = [1 0 0;0 cos(theta(3)) sin(theta(3)); 0 -sin(theta(3)) cos(theta(3))];
DCMp = Xrot*Yrot*Zrot;

isEqu3 = all((R1-DCMp)<1e-10,'all')

isEqu3 = logical
   0

It is bizarre that the quaternion q2 rotated around frame/intrinsic is numerically equal to the rotation matrix rotated by point/extrinsic. Look at the following example.

q2 = quaternion(theta,'euler','ZYX','frame') % is same as q1

q2 = quaternion
     0.94043 + 0.16458i + 0.18664j + 0.23166k

isEqu4 = all(rotmat(q2,'point')-myR<1e-10,'all')

isEqu4 = logical
   1

My conclusion is that the use of quaternions in matlab is not clear, although quat2rotm, rotm2quat can be normal recovery interconversion, pre-multiplication and post-multiplication in different toolboxes behave in different ways

[Paul]

@cui

Cotinuing the discussion ...

rng(100);
theta = rand(1,3); % your origin data
q1 = eul2quat(theta,'ZYX');
R1 = eul2rotm(theta,'ZYX');
R2 = quat2rotm(q1);
R1 - R2

ans = 3×3

1.0e-15 * 0.1110 0.0555 -0.1665 -0.0555 0 0.0278 0.1110 -0.1110 0

So R1 and R2 are equal. According to the doc page for quat2rotm R2 is to be used by premultiplying by a row vector, which I thought I illustrated above. The doc page for eul2rotm is silent on this issue, but R1 == R2, so I'll assume it's the same.

Your comment above is "The premise is ... that the point coordinates are in the right-hand side of the rotation matrix in the form of column vectors. " But R1 and R2 use coordinates on the left hand side as row vectors. So to be consistent with the premise, they both need to be transposed

R1p = transpose(R1);
R2p = transpose(R2);

Next, we have my diefintions of Zrot, etc.

Zrot = [cos(theta(1)) sin(theta(1)) 0; -sin(theta(1)) cos(theta(1)) 0;0 0 1];
Yrot = [cos(theta(2)) 0 -sin(theta(2)); 0 1 0; sin(theta(2)) 0 cos(theta(2))];
Xrot = [1 0 0;0 cos(theta(3)) sin(theta(3)); 0 -sin(theta(3)) cos(theta(3))];

As we can see, these are elemental frame rotations that multiply a column vector. For example:

format long e
Zrot - rotmat(quaternion([theta(1) 0 0],'euler','ZYX','frame'),'frame')

ans = 3×3

1.0e+00 * -2.220446049250313e-16 -1.110223024625157e-16 0 1.110223024625157e-16 -2.220446049250313e-16 0 0 0 -2.220446049250313e-16

I'm shocked that the (3,3) element is identically zero. Regardless, DCMp is a frame rotation applied to column vector

DCMp = Xrot*Yrot*Zrot;

Now DCMp, R1p, and R2p are all the same

DCMp - R1p

ans = 3×3

0 0 0 0 0 0 0 0 0

DCMp - R2p

ans = 3×3

1.0e+00 * 1.110223024625157e-16 -5.551115123125783e-17 1.110223024625157e-16 5.551115123125783e-17 0 -1.110223024625157e-16 -1.665334536937735e-16 2.775557561562891e-17 0

Finally, that leave us with

thetaR = rad2deg(theta);
myR = rotz(thetaR(1))*roty(thetaR(2))*rotx(thetaR(3));% common geometric transformation convention

Now it should be apparent that myR is the transpose of the others

DCMp - transpose(myR)

ans = 3×3

1.0e+00 * 0 0 0 0 0 -5.551115123125783e-17 -5.551115123125783e-17 5.551115123125783e-17 0

Whether or not myR is an extrinsic/point/muliply-by-column-on-right rotation I can't say. But it's clear from the formation of myR why it's the transpose of DCMp (and R1p and R2p). But, we can see that myR can be constructed from theta as a point rotation, but as @Brian Fanous state quaternion() uses intrinsic rotations.

myR - rotmat(quaternion(fliplr(theta),'euler','XYZ','point'),'point')

ans = 3×3

1.0e+00 * 0 -5.551115123125783e-17 0 5.551115123125783e-17 0 -5.551115123125783e-17 0 0 0

I don't disagree with this statement: "the document[ation] is not uniformly clear"

[Brian Fanous]

@Paul, @cui,

Is the answer I supplied above sufficient to include in the documentation? Are there other details you would like included? We'd like to make this as clear as possible.

[cui,xingxing]

The expression through the formula is clear enough, thank you for the discussion@Paul@Brian Fanous@Sam Chak

BTW，this similar issue “why point correspondence to rotation matrix calculation is not my expected?” need official clarification. @Brian Fanous

[Paul]

@Brian Fanous,

Since you aksed, and thanks for doing so ....

Please keep in mind that I don't claim to be familiar with all of the relevant toolboxes and doc pages. But from just poking around in response to this Question ....

Each relevant toolbox (Navigation, Aerospace, UAV, etc.) should have a standalone section that desribes the the conventions and definitions used for that toolbox.

The function doc pages in a toolbox that have anything to do with roation matrices need to clearly state if the convention is post-multiply by column or pre-multiply by row, unless the toolbox only ever uses one convention, which can be stated clearly in the standalone section. Similarly, for point/frame. Each doc page related to rotations should have a link back to the standalone page.

In at least one instance, the meaning of pre-mutiply by row was called a "post-multiply format," which was confusing.

It would be good for the standalone section to explicitly write out the single axis rotations (e.g., the R_x, R_y, R_z in your answer above) that the toolbox uses, and perhaps adopt different notaton for the elemental frame (R ) vs. point (T, as in your answer) if appropriate.

To answer your specific question, yes, have the doc define frame and point rotation as you've done in your Answer would be extremely helpful.

I do realize that documenting everything consistently and clearly across all relevant toolboxes could be a daunting task ...

Answer 4

[cui,xingxing]

"eul2*" series of built-in functions such as the order and expression of customary way ambiguous, for the time being, do not discuss the use of the function, waiting for the official documentation update to clarify.

Typical problems such as rigid3d object rotation matrix for postmultiply form, rotx/y/z, eul2* series functions for premultiply form.The form is not uniform, easy to cause misunderstanding, the best way to write their own rotation matrix

issue1 : https://ww2.mathworks.cn/matlabcentral/answers/1720045-how-to-get-the-relative-camera-pose-to-another-camera-pose#answer_964925

issue2: https://ww2.mathworks.cn/matlabcentral/answers/1716235-why-point-correspondence-to-rotation-matrix-calculation-is-not-my-expected#answer_961530

issue3: https://ww2.mathworks.cn/matlabcentral/answers/1726500-why-do-matlab-and-python-convert-euler-angles-to-quaternions-with-different-results?s_tid=srchtitle

https://ww2.mathworks.cn/help/driving/ug/rotations-using-quaternions-in-automated-driving.html