Code covered by the BSD License  

Highlights from
Function to Convert between DCM, Euler angles, Quaternions, and Euler vectors

4.61538

4.6 | 13 ratings Rate this file 146 Downloads (last 30 days) File Size: 5.71 KB File ID: #20696

Function to Convert between DCM, Euler angles, Quaternions, and Euler vectors

by

 

14 Jul 2008 (Updated )

Function to convert rotation data between 4 types: DCM, Euler Angles, Quaternions, and Euler Param.

| Watch this File

File Information
Description

SpinCalc is a consolidated matlab function that will convert any rotation data between the 4 types included. Will also convert between 2 different Euler angle set types.

Multiple orientations can be input. For N orientations:
DCM ===> 3x3xN multidimensional array
EA### ===> Nx3 matrix
Euler Vectors ===> Nx4 matrix
Quaternions ===> Nx4 matrix

Input includes error check flag that will warn when Euler angles approach singularity or when appropriate values deviate from unity. Fatal errors issued for improper DCM's etc.

*****NOTE TO USERS*****
I have gotten many questions regarding translation to Euler angle sets. When converting data to Euler angles, you MUST make sure the orientation you are translating is not near a singularity. Singular Euler sets are orientations which cannot be uniquely converted to 3 variables in that particular rotation order. The singular sets are as follows:

Type 1 Rotations: 123 - 132 - 213 - 231 - 321 - 312
Singular if second rotation angle is -90 or 90 degrees.

Type 2 Rotations: 121 - 131 - 212 - 232 - 313 - 323
Singular if second rotation angle is 0 or 180 degrees.

SpinCalc should now detect when input DCM, EV, or Q correspond to a singular Euler set output. It will prohibit output in such an event.

Naturally when converting from these singular Euler angle sets to other data types, you will receive a correct answer. Unfortunately you cannot convert that output back to the correct Euler angle set. This is why singular Euler input is prohibited.

Some users may need SpinCalc for the very reason of converting singular Euler sets to the other types. In those cases, removing the prohibiting error check shouldn't be too difficult. Contact me for help if need be.

Acknowledgements

This file inspired Euler Angle, Dcm, Quaternion, And Euler Vector Conversion/Teaching Gui, Spin Conv, and Quaternion.M.

MATLAB release MATLAB 7.5 (R2007b)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (43)
16 Oct 2014 Daniel Lucena

What am I supposed to add instead of tol and ichk? I replaced with 0.0001 and 0 respectvely but it doesn't work. Thanks.

31 May 2014 Jing

it is very helpful

09 Apr 2014 spencer

Are the euler angles defined as "body-fixed" or "world-fixed"?

27 Mar 2014 John Fuller

Steve,
Honestly I just use "eps" for most purposes. 9 times out of 10 it doesn't constrain the problem too much. When it does fail, you can ratchet it higher for compliance.

Thanks,
John

27 Mar 2014 Steve

Hey, thanks so much for sharing this code! It helps me a lot.
but I have one single question:
Can you give a recommendation for a tolerance value, for a conversion from rotationmatrix to Eulerangles?
Thanks!

29 Jan 2013 Paolo de Leva

That's bad. My second last posting was totally deleted by MATLAB Central, possibly because it contained special characters...

I wrote that in your latest release of SpinCalc, you introduced a new bug in line 238.

In that line (which I cannot copy here because it contains special characters), you should use the element-by-element OR operator (single vertical bar), rather than the short-circuit OR operator (double bar) which does not work when N is greater than 1.

In short, you should go back to previous version of line 238.

Please also check my 4 advices posted yesterday.

29 Jan 2013 Paolo de Leva

MATLAB Central keeps truncating everything I write after special characters. In my latest posting, the logical OR operators and whatever was written after them were deleted. Anyway, you can check previous version of line 238. You should use the element-by-element OR operator (single vertical bar), rather than the short-circuit OR operator (double bar).

28 Jan 2013 Paolo de Leva

I have been trying to post this text, but MATLAB Central keeps truncating it. Here is the complete version, with special characters removed:

1) Please update the release version number.

2) Line 388:

OUTPUT=mod([psi,theta,phi]*180/pi,360);

Proposed change:

OUTPUT=[psi,theta,phi]*180/pi;

THETA is computed in your code either using ACOS or ASIN. In MATLAB, both ACOS(X) and ASIN(X) return values within a limited range. For -1 =< X =< 1,

ACOS_X = ACOS(X)*(180/pi)

always returns values ranging from 0 to 180 degrees. This almost coincides with the range (from 0.1 to 179.9 degrees) accepted by SpinCalc as input for THETA (with EA type 2 rotations). Similarly, for -1 =< X =< 1,

ASIN_X = ASIN(X)*(180/pi)

always returns angles ranging from -90 to 90 degrees. Again, this almost coincides with the range (from -89.9 to 89.9 degrees) accepted by SpinCalc as input for THETA (with EA type 1 rotations). If X < -1 or X > 1, then the result ACOS_X or ASIN_X is complex, and SpinCalc returns an error message.
Computing MOD(ACOS_X,360) is useless, because

THETA = MOD(ACOS_X, 360) = ACOS_X

More importantly, computing MOD(ASIN_X,360) is inconsistent with other parts of your code, because for negative values of ASIN_X (i.e. -90 =< ASIN_X < 0 degrees)

THETA = MOD(ASIN_X, 360) == ASIN_X + 360

which means that, for negative values of ASIN_X, THETA will range from 270 =< ASIN_X < 360 degrees, which is considered to be out range when used as input for SpinCalc (with EA type 1 rotations).

3) Line 396:

sing_chk=...
[find(abs(theta*180/pi) < 0.1); ...
find(abs(theta*180/pi-180) < 0.1);...
find(abs(theta*180/pi-360)) < 0.1];

Proposed change:

sing_chk=...
[find(abs(theta*180/pi) < 0.1); ...
find(abs(theta*180/pi-180) < 0.1);...
find(abs(theta*180/pi-360) < 0.1)];

4) Every time you convert a Nx1 or Nx3 array containing angles in radians into an array containing angles in degrees, or vice versa (and you do it a lot of times in your code, including the above mentioned lines 388 and 396), you use this syntax:

Y = X*180/pi
or
Y = X*pi/180

This way, you repeat the division (by PI or by 180) N times. I suggest you to use this sintax, which is much more efficient as it computes the division only once (see also functions deg2rad and rad2deg):

Y = X*(180/pi)
or
Y = X*(pi/180).

28 Jan 2013 Paolo de Leva

I am trying to post this text, containing a few notes about SpinCalc, but the web site keeps truncating it. Here is the complete version, with a few special characters removed:

1) Please update the release version number.

2) Line 388:

OUTPUT=mod([psi,theta,phi]*180/pi,360);

Proposed change:

OUTPUT=[psi,theta,phi]*180/pi;

THETA is computed in your code either using ACOS or ASIN. In MATLAB, both ACOS(X) and ASIN(X) return values within a limited range. For -1 <= X <= 1,

ACOS_X = ACOS(X)*(180/pi)

always returns values ranging from 0 to 180 degrees. This almost coincides with the range (from 0.1 to 179.9 degrees) accepted by SpinCalc as input for THETA (with "input type" EA### and "Euler type" 2). Similarly, for -1 <= X <= 1,

ASIN_X = ASIN(X)*(180/pi)

always returns angles ranging from -90 to 90 degrees. Again, this almost coincides with the range (from -89.9 to 89.9 degrees) accepted by SpinCalc as input for THETA (with "input type" EA### and "Euler type" 1). If X > 1 or X < -1, then the result ACOS_X or ASIN_X is complex, and SpinCalc returns an error message.
Computing MOD(ACOS_X,360) is useless, because

THETA = MOD(ACOS_X, 360) = ACOS_X

More importantly, computing MOD(ASIN_X,360) is inconsistent with other parts of your code, because for negative values of ASIN_X (i.e. -90 <= ASIN_X < 0 degrees)

THETA = MOD(ASIN_X, 360) == ASIN_X + 360

which means that, for negative values of ASIN_X, THETA will range from 270 <= ASIN_X < 360 degrees, which is considered to be out range when used as input for SpinCalc (with "input type" EA### and "Euler type" 1).

3) Line 396:

sing_chk=...
[find(abs(theta*180/pi) < 0.1); ...
find(abs(theta*180/pi-180) < 0.1);...
find(abs(theta*180/pi-360)) < 0.1];

Proposed change:

sing_chk=...
[find(abs(theta*180/pi) < 0.1); ...
find(abs(theta*180/pi-180) < 0.1);...
find(abs(theta*180/pi-360) < 0.1)];

4) Every time you convert a Nx1 or Nx3 array containing angles in radians into an array containing angles in degrees, or vice versa (and you do it a lot of times in your code, including the above mentioned lines 388 and 396), you use this syntax:

Y = X*180/pi
or
Y = X*pi/180

This way, you repeat the division (by PI or by 180) N times. I suggest you to use this sintax, which is much more efficient as it computes the division only once (see also functions deg2rad and rad2deg):

Y = X*(180/pi)
or
Y = X*(pi/180).

28 Jan 2013 Paolo de Leva

I am trying to post this text, containing a few notes about SpinCalc, but the web site keeps truncating it. Here is the complete version, with a few special characters removed:

1) Please update the release version number.

2) Line 388:

OUTPUT=mod([psi,theta,phi]*180/pi,360);

Proposed change:

OUTPUT=[psi,theta,phi]*180/pi;

THETA is computed in your code either using ACOS or ASIN. In MATLAB, both ACOS(X) and ASIN(X) return values within a limited range. For -1<=X<=1,

ACOS_X = ACOS(X)*(180/pi)

always returns values ranging from 0 to 180 degrees. This almost coincides with the range (from 0.1 to 179.9 degrees) accepted by SpinCalc as input for THETA (with "input type" EA### and "Euler type" 2). Similarly, for -1<=X<=1,

ASIN_X = ASIN(X)*(180/pi)

always returns angles ranging from -90 to 90 degrees. Again, this almost coincides with the range (from -89.9 to 89.9 degrees) accepted by SpinCalc as input for THETA (with "input type" EA### and "Euler type" 1). If X>1 or X<-1, then the result ACOS_X or ASIN_X is complex, and SpinCalc returns an error message.
Computing MOD(ACOS_X,360) is useless, because

THETA = MOD(ACOS_X, 360) = ACOS_X

More importantly, computing MOD(ASIN_X,360) is inconsistent with other parts of your code, because for negative values of ASIN_X (i.e. -90<=ASIN_X<0 degrees)

THETA = MOD(ASIN_X, 360) == ASIN_X + 360

which means that, for negative values of ASIN_X, THETA will range from 270<=ASIN_X<360 degrees, which is considered to be out range when used as input for SpinCalc (with "input type" EA### and "Euler type" 1).

3) Line 396:

sing_chk=...
[find(abs(theta*180/pi)<0.1); ...
find(abs(theta*180/pi-180)<0.1);...
find(abs(theta*180/pi-360))<0.1];

Proposed change:

sing_chk=...
[find(abs(theta*180/pi)<0.1); ...
find(abs(theta*180/pi-180)<0.1);...
find(abs(theta*180/pi-360)<0.1)];

4) Every time you convert a Nx1 or Nx3 array containing angles in radians into an array containing angles in degrees, or vice versa (and you do it a lot of times in your code, including the above mentioned lines 388 and 396), you use this syntax:

Y = X*180/pi
or
Y = X*pi/180

This way, you repeat the division (by pi or by 180) N times. I suggest you to use this sintax, which is much more efficient as it computes the division only once:

Y = X*(180/pi)
or
Y = X*(pi/180).

27 Jan 2013 Paolo de Leva

2) Line 388:

OUTPUT=mod([psi,theta,phi]*180/pi,360);

Proposed change:

OUTPUT=[psi,theta,phi]*180/pi;

THETA is computed in your code either using ACOS or ASIN. In MATLAB, both ACOS(X) and ASIN(X) return values within a limited range. For -1<=X<=1,

ACOS_X = ACOS(X)*(180/pi)

always returns values ranging from 0 to 180 degrees. This almost coincides with the range (from 0.1 to 179.9 degrees) accepted by SpinCalc as input for THETA (with "input type" EA### and "Euler type" 2). Similarly, for -1<=X<=1,

ASIN_X = ASIN(X)*(180/pi)

always returns angles ranging from -90 to 90 degrees. Again, this almost coincides with the range (from -89.9 to 89.9 degrees) accepted by SpinCalc as input for THETA (with "input type" EA### and "Euler type" 1).
If X>1 or X<-1, then the result ACOS_X or ASIN_X is complex, and SpinCalc returns an error message.
Computing MOD(ACOS_X,360) is useless, because

THETA = MOD(ACOS_X, 360) = ACOS_X

More importantly, computing MOD(ASIN_X,360) is inconsistent with other parts of your code, because for negative values of ASIN_X (i.e. -90<=ASIN_X<0 degrees)

THETA = MOD(ASIN_X, 360) == ASIN_X + 360

which means that, for negative values of ASIN_X, THETA will range from 270<=ASIN_X<360 degrees, which is considered to be out range when used as input for SpinCalc (with "input type" EA### and "Euler type" 1).

27 Jan 2013 Paolo de Leva

A few other advices about your code:
--------

Please update the release version number.

--------

Line 388:

OUTPUT=mod([psi,theta,phi]*180/pi,360);

Proposed change:

OUTPUT=[psi,theta,phi]*180/pi;

THETA is computed in your code either using ACOS or ASIN. In MATLAB, both ACOS(X) and ASIN(X) return values within a limited range. For -1<=X<=1,

ACOS_X = ACOS(X)*(180/pi)

always returns values ranging from 0 to 180 degrees. This almost coincides with the range (from 0.1 to 179.9 degrees) accepted by SpinCalc as input for THETA (with "input type" EA### and "Euler type" 2). Similarly, for -1<=X<=1,

ASIN_X = ASIN(X)*(180/pi)

always returns angles ranging from -90 to 90 degrees. Again, this almost coincides with the range (from -89.9 to 89.9 degrees) accepted by SpinCalc as input for THETA (with "input type" EA### and "Euler type" 1).
If X>1 or X<-1, then the result ACOS_X or ASIN_X is complex, and SpinCalc returns an error message.
Computing MOD(ACOS_X,360) is useless, because

THETA = MOD(ACOS_X, 360) = ACOS_X

More importantly, computing MOD(ASIN_X,360) is inconsistent with other parts of your code, because for negative values of ASIN_X (i.e. -90<=ASIN_X<0 degrees)

THETA = MOD(ASIN_X, 360) == ASIN_X + 360

which means that, for negative values of ASIN_X, THETA will range from 270<=ASIN_X<360 degrees, which is considered to be out range when used as input for SpinCalc (with "input type" EA### and "Euler type" 1).

---------

Line 396:

sing_chk=...
[find(abs(theta*180/pi)<0.1); ...
find(abs(theta*180/pi-180)<0.1);...
find(abs(theta*180/pi-360))<0.1];

Proposed change:

sing_chk=...
[find(abs(theta*180/pi)<0.1); ...
find(abs(theta*180/pi-180)<0.1);...
find(abs(theta*180/pi-360)<0.1)];

--------------------------------

Every time you convert a Nx1 or Nx3 array containing angles in radians into an array containing angles in degrees, or vice versa (and you do it a lot of times in your code, including the above mentioned lines 388 and 396), you use this syntax:

Y = X*180/pi
or
Y = X*pi/180

This way, you repeat the division (by pi or by 180) N times. I suggest you to use this sintax, which is much more efficient as it computes the division only once:

Y = X*(180/pi)
or
Y = X*(pi/180).

27 Jan 2013 Paolo de Leva

Line 388:

OUTPUT=mod([psi,theta,phi]*180/pi,360);

Proposed change:

OUTPUT=[psi,theta,phi]*180/pi;

THETA is computed in your code either using ACOS or ASIN. In MATLAB, both ACOS(X) and ASIN(X) return values within a limited range. For -1<=X<=1,

ACOS_X = ACOS(X)*(180/pi)

always returns values ranging from 0 to 180 degrees. This almost coincides with the range (from 0.1 to 179.9 degrees) accepted by SpinCalc as input for THETA (with "input type" EA### and "Euler type" 2). Similarly, for -1<=X<=1,

ASIN_X = ASIN(X)*(180/pi)

always returns angles ranging from -90 to 90 degrees. Again, this almost coincides with the range (from -89.9 to 89.9 degrees) accepted by SpinCalc as input for THETA (with "input type" EA### and "Euler type" 1).
If X>1 or X<-1, then the result ACOS_X or ASIN_X is complex, and SpinCalc returns an error message.
Computing MOD(ACOS_X,360) is useless, because

THETA = MOD(ACOS_X, 360) = ACOS_X

More importantly, computing MOD(ASIN_X,360) is inconsistent with other parts of your code, because for negative values of ASIN_X (i.e. -90<=ASIN_X<0 degrees)

THETA = MOD(ASIN_X, 360) == ASIN_X + 360

which means that, for negative values of ASIN_X, THETA will range from 270<=ASIN_X<360 degrees, which is considered to be out range when used as input for SpinCalc (with "input type" EA### and "Euler type" 1).

27 Jan 2013 Paolo de Leva

A few other advices about your code:
------------------------------------
Please update the release version number.
------------------------------------
Line 388:

OUTPUT=mod([psi,theta,phi]*180/pi,360);

Proposed change:

OUTPUT=[psi,theta,phi]*180/pi;

THETA is computed in your code either using ACOS or ASIN. In MATLAB, both ACOS(X) and ASIN(X) return values within a limited range. For -1<=X<=1,

ACOS_X = ACOS(X)*(180/pi)

always returns values ranging from 0 to 180 degrees. This almost coincides with the range (from 0.1 to 179.9 degrees) accepted by SpinCalc as input for THETA (with "input type" EA### and "Euler type" 2). Similarly, for -1<=X<=1,

ASIN_X = ASIN(X)*(180/pi)

always returns angles ranging from -90 to 90 degrees. Again, this almost coincides with the range (from -89.9 to 89.9 degrees) accepted by SpinCalc as input for THETA (with "input type" EA### and "Euler type" 1).
If X>1 or X<-1, then the result ACOS_X or ASIN_X is complex, and SpinCalc returns an error message.
Computing MOD(ACOS_X,360) is useless, because

THETA = MOD(ACOS_X, 360) = ACOS_X

More importantly, computing MOD(ASIN_X,360) is inconsistent with other parts of your code, because for negative values of ASIN_X (i.e. -90<=ASIN_X<0 degrees)

THETA = MOD(ASIN_X, 360) == ASIN_X + 360

which means that, for negative values of ASIN_X, THETA will range from 270<=ASIN_X<360 degrees, which is considered to be out range when used as input for SpinCalc (with "input type" EA### and "Euler type" 1).

--------------------------------

Line 396:

sing_chk=...
[find(abs(theta*180/pi)<0.1); ...
find(abs(theta*180/pi-180)<0.1);...
find(abs(theta*180/pi-360))<0.1];

Proposed change:

sing_chk=...
[find(abs(theta*180/pi)<0.1); ...
find(abs(theta*180/pi-180)<0.1);...
find(abs(theta*180/pi-360)<0.1)];

--------------------------------

Every time you convert a Nx1 or Nx3 array containing angles in radians into an array containing angles in degrees, or vice versa (and you do it a lot of times in your code, including the above mentioned lines 388 and 396), you use this syntax:

Y = X*180/pi
or
Y = X*pi/180

This way, you repeat the division (by pi or by 180) N times. I suggest you to use this sintax, which is much more efficient as it computes the division only once:

Y = X*(180/pi)
or
Y = X*(pi/180).

22 Jan 2013 John Fuller

Paolo,
I have incorporated a fix to address that bug, which is largely similar to your snippet. Thanks for catching this!

11 Jan 2013 Paolo de Leva

There is a bug in the IF statements used to check for singularities in the input values for Euler angles. These statements are only working correctly if INPUT size is 1x3. They do not work as expected when INPUT size is Nx3 (with N>1).

FIRST EXAMPLE

Suppose that rotation order type for Euler angles is 2 (first and third rotations about same axis). Let

INPUT = [10 -10 10] % (1x3) row vector

The following IF statement, used in SpinCalc (source code lines 107,108,109) correctly spots the singularity (second Euler angle is -10, hence out of range) and generates an error message:

if INPUT(:,2)<=zeros(size(INPUT,1),1) |...
INPUT(:,2)>=180*ones(size(INPUT,1),1)
error('Error: Second input Euler angle(s) outside 0 to 180 degree range')
end

SECOND EXAMPLE

Let INPUT be a Nx3 list of row vectors (N>1), such as

INPUT = [10 -10 10
10 00 10
10 10 10
10 180 10
10 190 10] % (5x3 row vector list)

The above mentioned IF statement does not spot the singularities in the first, second, fourth and fifth row (where second Euler angle is out of range), and does not generate an error message.

Spincalc contains a group of similar IF statements which contain the same bug (source code lines from 105 to 121).

Here is my solution:

input_size = size(INPUT);
N = input_size(1);

% Identify singularities (second Euler angle out of range)
EA2 = INPUT(:,2); % (Nx1) 2nd Euler angle(s)
ZEROS = zeros(N,1); % (Nx1)
ONES = ones(N,1); % (Nx1)
if rot_1_in==rot_3_in % Type 2 rotation (1st and 3rd rotations about same axis)
if any(EA2<=ZEROS) || any(EA2>=180*ONES)
error('Second input Euler angle(s) outside 0 to 180 degree range')
elseif any(abs(EA2)<2*ONES) || any(abs(EA2)>178*ONES)
if ichk==1
errordlg('Warning: Second input Euler angle(s) near a singularity (0 or 180 degrees).')
end
end
else % Type 1 rotation (rotations about three distinct axes)
if any(abs(EA2)>=90*ONES)
error('Second input Euler angle(s) outside -90 to 90 degree range')
elseif any(abs(EA2)>88*ONES)
if ichk==1
errordlg('Warning: Second input Euler angle(s) near a singularity (-90 or 90 degrees).')
end
end
end

21 Nov 2012 J. Hu

For some rotation matrix INPUT, certain CONVERSION will lead error messages due to singularity. And it has to use other CONVERSION type. It would be helpful if the code can automatically select a error-free CONVERSION for any INPUT.

For example, when rotation matrix INPUT =
[0 0 1;
1 0 0;
0 1 0]
The conversion has to be 'DCMtoEA2xx' in which 2 needs to be the first and xx can be 13 or 31.

When the code is called for repetitive computation of many data sets, an auto selected CONVERSION type would be desired.

Thanks...

25 Oct 2012 Georg Wiora  
31 May 2012 James Coburn  
17 Nov 2011 zhaopeng QIU  
10 Oct 2011 Paul

Thanks for the code, it works great.

Funny fact: the rotation matrix returned by the procrustes function in the Statistics Toolbox is quite often not orthogonal to eps.

05 Jan 2011 Paolo de Leva

This is even better:

"DCM - 3x3xN multidimensional direction cosine matrix which performs coordinate system rotations by pre-multiplying column vectors (V_rot = DCM * V, where V_rot is V resolved in the rotated coordinate system)."

Paolo de Leva

05 Jan 2011 Paolo de Leva

I suggest this help text:

"All rotation representations used by SpinCalc represent the rotation of the coordinate system."

"DCM - 3x3xN multidimensional matrix which pre-multiplies column vectors to rotate the reference system in which they are represented."

Paolo de Leva

Paolo de Leva

05 Jan 2011 John Fuller

Alright, I don't really have the time to go over semantics in the help text of my code. I can agree that the following block

% DCM - 3x3xN multidimensional matrix which pre-multiplies a coordinate
% frame column vector to rotate it to the desired new frame.

might be better as

% DCM - 3x3xN multidimensional matrix which pre-multiplies a coordinate
% frame column vector to calculate its coordinates in the desired
% new frame.

But the fact still remains that the input DCM data remains the same. I wouldn't say it is highly misleading, but that it also can't be expected to teach every user the nuances of rotation linear algebra as you have illustrated.

The code uses input and output in an accepted convention that 'most' serious users in my industry follow and understand, and it stands simply as that. As with all Matlab FEX code it is very important that the user vets code cautiously before use and knows what goes in and out by external testing.

I will update the commenting above, but I'm not pursuing further modifications to SpinCalc because I believe it is well-accepted, heavily tested, and straightforward as it is. If you don't like that, rate it badly and don't use it.

Thanks,
John

05 Jan 2011 Paolo de Leva

OK, John,

I GUESS THIS IS OUR COMMON BACKGROUND:

We agree that the rotation from a reference frame A (original) to another reference frame B (rotated) produces an apparent rotation in the opposite direction of a vector v which is fixed in A. This leads to two conventions to describe rotation. Let’s call them:
1) "FRAME ROTATION"
2) “VECTOR ROTATION" (or “vector space rotation”)
The most important point is that, when we describe a rotation from A to B, we must assume that v is fixed in A (see below). Also, we agree that SpinCalc uses the “FRAME ROTATION” convention, and of course we agree on Rodriguez's rotation formula and all other conversion formulas which are correctly implemented in SpinCalc.

Moreover, as discussed in the previous message, we agree that there are two conventions used to represent rotation with a DCM (pre- or post-multiplication). Your recently updated help text makes clear you are using pre-multiplication of a column vector by the DCM. IF we call R the DCM, and v the original column vector, then R is a matrix such that:

v_rot = R*v.

NOW, HERE'S WHAT WE DISAGREE ABOUT:

We do not agree about your help text. In your second last message, you explained that you implemented the “FRAME ROTATION” convention, but this is not what you wrote in your help text, which defines the DCM as:

"a multidimensional matrix which pre-multiplies a coordinate frame (column) vector to rotate it to the desired new frame."

Let’s call A the initial frame, and B “the desired new frame“.
It unequivocally follows from your text that, in the equation v_rot = R*v, you describe v as one of the versors of A (i.e., a standard basis vector; i.e., either [1,0,0], [0,1,0], or [0,0,1], in 3D). Actually, in this singular case R*v is the same as selecting one of the columns of R, and this description of v is quite questionable, in my opinion, as a (pre- or post-) multiplication by R is typically used to rotate any vector, not only the versors of A. But this is not the most important fault in your definition.

Since your help text describes the rotation of a "column vector", then you unequivocally defined R as a matrix which produces (i.e. represents) a “VECTOR ROTATION”, not a “FRAME ROTATION”. Moreover, the “FRAME ROTATION” convention assumes that v “is fixed in A” (see above). On the contrary, according to your definition, v unequivocally appears fixed in B (not in A).

Most importantly and unfortunately, R does not represent a “FRAME ROTATION” from A to B, but a “FRAME ROTATION” of equal magnitude in the opposite sense, from A to C. Namely, C does not coincide with B, as B rotates together with the vector. For instance, if B is rotated by 45° about z, then C is rotated by -45° about the same axis. Thus, we can rewrite the equation as follows:

vC = R*vA.
(vB = R*vA would be incorrect).

In short, your help text describes a “FRAME ROTATION” from A to B, which is the same as (not the opposite of) a the “VECTOR ROTATION” rotation produced by R (from v to v_rot, or from vA to vC, if you prefer). But unfortunately R neither produces nor represents the frame rotation from A to B.

However, your axis-angle representation describes the opposite rotation, the “FRAME ROTATION” from A to C, in which v is “fixed in A” and observed from C. This is the frame rotation produced and represented by R, but C is not even mentioned in your help text! In other words, your definition of rotation is inconsistent with your implementation, and this is higly misleading.

Paolo de Leva

Paolo

02 Jan 2011 John Fuller

Another note, if you take a look at my other submission called SpinCalcVis it will give a good idea of how input/output is processed by SpinCalc as it is all based on the same parent code.

http://www.mathworks.com/matlabcentral/fileexchange/27653-euler-angle-dcm-quaternion-and-euler-vector-conversionteaching-gui

Thanks,
John

02 Jan 2011 John Fuller

Paolo,
The convention used by SpinCalc is the standard convention where matrix R pre-multiplies the column vector, or R*v. The "by" in the help section could be a bit ambiguous. I can change that when I have the master file back in front of me.

The reason for the inconsistency you are seeing deals with the fact that the R matrix format used in SpinCalc is for the DCM which rotates the coordinate frame about the vector, not the vector about the coordinate frame. Thus, the Euler vector solution corresponds to the vector about which to rotate the frame and by how much.

Your matrix to rotate [1;0;0] 90 degrees about Z is equivalent to rotating the frame -90 about Z, or +90 about -Z. All data input and output for SpinCalc corresponds to rotation of the frame.

02 Jan 2011 Paolo de Leva

The help text of SpinCalc defines a DCM (rotation matrix) as follows:

"3x3xN multidimensional matrix which pre-multiplies by a coordinate frame vector to rotate it to the desired new frame."

The expression “pre-multiplies by” seems ambiguous to me. In my opinion, a rotation matrix R either "pre-multiplies" a column vector v1 (R*v1), or "is pre-multiplied by" a row vector v2 (v2*R).

What did you mean? R*v or v*R?
Would you mind to you make it clear in your help text?

It is quite important to give a crystal clear definition of the DCM. In linear algebra, the most commonly used convention to perform a rotation is by means of a pre-multiplication of a column vector by a DCM. Here's what I get when I follow this convention:

If i, j, k are the versors of a coordinate system xyz, then a rotation matrix R that rotates i [1;0;0] by 90° about the z axis produces a new versor i_new which coincides with j [0;1;0]:

R = [0 -1 0
1 0 0
0 0 1];
i = [1; 0; 0]; % Column vector
i_new = R * i;
disp(i_new)
0
1
0

Thus, if I use SpinCalc to compute the axis-angle vector corresponding to R, I expect it to be [0 0 1 90], however:

SpinCalc('DCMtoEV', R, 1e-5, 0)

ans =

0 0 -1 90

How do you explain this apparent inconsistency? Is this due to the fact that you are not using the above mentioned convention? Again, can you make this clear in your help text?

20 Dec 2010 John Fuller

Gaetano,
I reviewed this line and don't see why it would produce an error as it is written:

if abs(det(INPUT(:,:,ii))+1)<0.05, %if determinant is near -1, DCM is improper

The line you have proposed would actually cause errors for proper DCM matrices, as their determinant should be equal to +1. In your line, the subtraction would cause proper DCMs to be thrown out, and the additional abs function is redundant.

Thanks,
John

19 Dec 2010 Gaetano Giulio Violante

Recently i downloaded this MatLab code, I used it to get euler angles from rotation matrix.

In my opinion there's an error in the code at line 165, and precisely i would modify it in the following way:

" abs(abs(det(INPUT(:,:,ii)))-1)<0.00005"

Anyone can tell me if this observation is correct?

Thank you,
Gaetano.

[Student in Telecommunications Engineering , Polythecnic of Bari - Bari, Italy]

30 Aug 2010 chandrakala Gowda  
02 Feb 2010 Oliver Woodford  
02 Feb 2010 Oliver Woodford

Very useful. Thanks.

20 Oct 2009 bryan

I found this code to be very useful and it seems to work exactly as described. Thank you for your excellent contribution.

15 Jul 2009 James Tursa

After comparing this submission with others in the FEX, I feel the need to bump my rating up.

30 Jun 2009 Joel  
28 Jun 2009 Jørgen

Okay, thank you for your feedback.

26 Jun 2009 John Fuller

Jørgen,
The 323 rotation you specified is a singular Euler set, because the middle angle is 0. Any type 2 rotation where the middle angle is 0 or 180 cannot be uniquely resolved by trying to translate Q, DCM, or EV back to Euler angles. If you plan on translating amongst orientations that are singular when expressed as Euler angles, I would advise you use strictly DCM, Q, or EV because they can uniquely define all orientations.

26 Jun 2009 Jørgen

When using the DCM to EA323 I get some unexpected results.

This is my code:

>> R = [0 -1 0; 1 0 0; 0 0 1]
R =

0 -1 0
1 0 0
0 0 1

>> SpinCalc2('DCMtoEA323',R,eps,1)

ans =

0 0 0

------

The matrix R represents a simple 90 degrees rotation about the z-axis, but the Euler Angles are all zero.

19 Jun 2009 John Fuller

I am uploading an updated version with the boundaries revised on the first and third Euler angles. They must now be between -360 and 360. The bounds on the second Euler angles are going to stay since they prevent ambiguous Euler sets.

Thanks for the comments.

20 May 2009 Stanley

Works well. Agreed with others that some bounds on inputs, while perhaps obscuring internal code, make implementation into other code more cumbersome.
for instance, in this snippet, ea(1) must be in 0->360.

SpinCalc('EA321toDCM',ea,eps,0)

Alowing ea(1) = -10 degrees is easier.
Perhaps a switch to turn of bounds checking?

29 Jul 2008 John Fuller

Thanks for the comments. I am editing the code to accommodate notes 2-5. I'd like to leave the second euler angle constraint between -90 and 90 since reaching these bounds for Type 1 rotations causes singularities. Allowing angles past those bounds, while correct, makes detecting singularities more obscure.

26 Jul 2008 James Tursa

I spot tested all combinations of input/output type against independent code I happened to have and SpinCalc produced the same results, so I would conclude based on this spot checking that the code generally works as advertised. Thus I gave it a fairly high rating. I did have a few problems, but I think these can easily be improved:

------------------------------------------------------------------------------------
1) Sometimes the results were not invertible. For example:

>> q = [1 2 3 4]
q =
1 2 3 4
>> q = q / norm(q)
q =
0.18257418583506 0.36514837167011 0.54772255750517 0.73029674334022
>> SpinCalc('QtoEA213',q,eps,1)
ans =
1.0e+002 *
0.47726310993906 3.52337744339234 0.70346175941947
>> SpinCalc('EA213toQ',ans,eps,1)
??? Error using ==> SpinCalc
Error: Second input Euler angle(s) outside -90 to 90 degree range

In this case, and in other similar cases, the checks for valid Euler Angle input are too tight. Having the second input outside -90 to +90 may not make sense for some physical applications, but it is not an invalid input and should not produce an error.

------------------------------------------------------------------------------------
2) Calculating the euler rotation angle for small angle rotation quaternions is not accurate. The code in question is this line:

MU = 2*acos(Q(:,4));

For example:

>> q = [1 2 3 4e12]
q =
1.0e+012 *
0.00000000000100 0.00000000000200 0.00000000000300 4.00000000000000
>> q = q / norm(q)
q =
0.00000000000025 0.00000000000050 0.00000000000075 1.00000000000000
>> SpinCalc('QtoEV',q,eps,1)
ans =
1 0 0 0
>> SpinCalc('EVtoQ',ans,eps,1)
ans =
0 0 0 1

You can see in this case that the small angle was completely lost in the q(4) component, so using only that component to calculate the rotation angle is not a good method for small angle cases. One fix is to use the atan2 function instead. For example, suppose the following code was used instead of the code above:

QV = sqrt(sum(Q(:,1:3).*Q(:,1:3),2));
MU = 2*atan2(QV,Q(:,4));

Now do the same calculations as above:

>> SpinCalc('QtoEV',q,eps,1)
ans =
0.26726124191242 0.53452248382485 0.80178372573727 0.00000000010719
>> SpinCalc('EVtoQ',ans,eps,1)
ans =
0.00000000000025 0.00000000000050 0.00000000000075 1.00000000000000

You can see that the small angle is calculated properly. Caveat: The atan2 code above is just shown as an example for improvement. It may not be the best overall way to do this.

------------------------------------------------------------------------------------
3) Some of the algorithms used are not robust for inputs that deviate from unity, etc. For example, suppose we start with a proper direction cosine matrix and perturb it a bit:

>> q = [1 2 3 4]
q =
1 2 3 4
>> q = q / norm(q)
q =
0.18257418583506 0.36514837167011 0.54772255750517 0.73029674334022
>> dcm=SpinCalc('QtoDCM',q,eps,1)
dcm =
0.13333333333333 0.93333333333333 -0.33333333333333
-0.66666666666667 0.33333333333333 0.66666666666667
0.73333333333333 0.13333333333333 0.66666666666667
>> dcm = dcm + 1e-7*(rand(3,3)-0.5)
dcm =
0.13333337834626 0.93333333193158 -0.33333333768657
-0.66666669355282 0.33333337246323 0.66666661851703
0.73333334401759 0.13333335954302 0.66666669880738
>> SpinCalc('DCMtoQ',dcm,eps,1)
Warning: Input DCM(s) matrix not orthogonal to precision tolerance.
(etc.)
ans =
0.18257415540401 0.36514836686611 0.54772255130155 0.73029676324370
>> norm(ans)
ans =
1.00000000382748

The algorithm used for the conversion picks up on the fact that the dcm was not quite proper (good), but does not correct the output. So maybe normalizing the returned quaternion might be better in this case. One way to accomplish this efficiently would be to use Markley's algorithm that was recently published in J. Guidance Vol 31 No 2 pp 440-442. This algorithm produces a normalized quaternion result using only one square root.

------------------------------------------------------------------------------------
4) Some of the code could be written more efficiently. For example in this section:

m1=INPUT(:,1); m2=INPUT(:,2); m3=INPUT(:,3); MU=INPUT(:,4)*pi/180;
if sqrt(m1.^2+m2.^2+m3.^2)-ones(N,1)>tol*ones(N,1),
error('Input euler vector(s) components do not constitute a unit vector')
end
if MU<zeros(N,1) | MU>2*pi*ones(N,1),
error('Input euler rotation angle(s) not between 0 and 360 degrees')
end
Q=[m1.*sin(MU/2),m2.*sin(MU/2),m3.*sin(MU/2),cos(MU/2)];

Why pull out the columns individually in m1, m2, and m3? Why not just use INPUT(:,1:3) in the follow-on calculations? And why calculate sin(MU/2) multiple times? Could recode this to calculate sin(MU/2) only once. For instance:

Q=[INPUT(:,1:3).*sin(MU/2),cos(MU/2)];

------------------------------------------------------------------------------------
5) Multiple warning messages appear for systematic problems when one message would suffice. For example:

>> q = [1 2 3 4]
q =
1 2 3 4
>> q = q / norm(q)
q =
0.18257418583506 0.36514837167011 0.54772255750517 0.73029674334022
>> q = q + 1.e-7*(rand(1,4)-.5)
q =
0.18257418030539 0.36514838321335 0.54772258669887 0.73029678552152
>> qq = [q;q;q;q;q;q;q;q;q;q]
qq =
0.18257418030539 0.36514838321335 0.54772258669887 0.73029678552152
0.18257418030539 0.36514838321335 0.54772258669887 0.73029678552152
(etc.)
>> SpinCalc('QtoDCM',qq,eps,1)
ans(:,:,1) =
0.13333335251383 0.93333342235746 -0.33333337639558
-0.66666675533740 0.33333337341210 0.66666670795756
0.73333338560076 0.13333335997254 0.66666675384566
ans(:,:,2) =
0.13333335251383 0.93333342235746 -0.33333337639558
-0.66666675533740 0.33333337341210 0.66666670795756
0.73333338560076 0.13333335997254 0.66666675384566
(etc.)
>> SpinCalc('DCMtoQ',ans,eps,1)
Warning: Input DCM(s) matrix not orthogonal to precision tolerance.Warning: Input DCM(s) matrix not orthogonal to precision tolerance.
(etc.)
ans =
0.18257418458450 0.36514839177157 0.54772259953620 0.73029676840507
0.18257418458450 0.36514839177157 0.54772259953620 0.73029676840507
(etc.)

So now I have 10 error message boxes to clear. If the array was 1000 then I would have 1000 error message boxes to clear. Only one message per function call is needed for these types of cases.

My comments are admittedly a bit picky, but keep in mind that I think these can easily be addressed to improve the code and I gave it a fairly high rating.

James Tursa

Updates
22 Jul 2008

Forgot to include additional capability in description.

23 Jul 2008

Added clarifying comments to help section of file.

19 Jun 2009

Removed bounding constraints on first and third input Euler angles. Bounds on second angles remain intact due to possibilities of ambiguous input Euler sets.

26 Jun 2009

Revised file description to explain Euler angle singularities

30 Jun 2009

Version 1.3 Now detects when input DCM, Q, or EV are too close to Euler singularity. Prohibits output to Euler angles when second angle is within 0.1 degree of singular value.

30 Jun 2009

Modified v1.3 description

22 Apr 2010

Updated tags

22 Apr 2010

Updated tags

04 May 2010

Updated description.

20 May 2010

Fixed summary.

03 Jan 2011

Updated help text to clarify DCM convention.

05 Jan 2011

Updated help text. Again.

25 Jan 2013

Incorporated fix to 2nd Euler Angle checks.

Contact us