File Exchange

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

version 1.11 (5.71 KB) by

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

4.52381
20 Ratings

Updated

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.

Stefan Grandl

Yun Peng

### Yun Peng (view profile)

Dear John,

Thank you so much in advance :)

John Fuller

### John Fuller (view profile)

Devendra,
It sounds like your conversion string command doesn't match what the program needs to instruct it which operation to do. The conversion string should be something like:

'EA321toDCM'

Where there is always a "to" between the two rotation format types. The file has a header that lists all of the possible conversions.

Thanks,
John

Hello John,
Hope u can assist me please,

??? Input argument "CONVERSION" is undefined.

Error in ==> SpinCalc at 62
i_type=strfind(lower(CONVERSION),'to');

Devendra

Geoffrey Vaquette

John Fuller

### John Fuller (view profile)

Ali,
If you are looking for a DCM that pre-multiplies a vector to move that vector, it would appear like the DCM is transposed. This function is intended produce output that rotates the frame about the vector, not the vector within the frame. If you only require DCM's which rotate the vector within the frame, you will always need to transpose the output.

John

Ali Aghaeifar

### Ali Aghaeifar (view profile)

Conversion from quaternion to DCM generates transpose of desired matrix.

Thomas Corie

### Thomas Corie (view profile)

Very useful. I was attempting something similar, but this is what I had in mind and more. Thank you!

HAN SUN

Felix Goldberg

### Felix Goldberg (view profile)

Great work.

shangcheng "sam" wang

### shangcheng "sam" wang (view profile)

Very practical tool set. if literature reference could be provided, that would be awesome.

John Fuller

### John Fuller (view profile)

Tim, apologies for the late reply, I lost my login email for quite a while and stopped receiving updates regarding SpinCalc.

The value is indeed 2 degrees, and the commenting is likely outdated. I'll be submitting an updated version shortly that will tweak this. You may modify the error-check value as you see fit though (some people remove it altogether).

Tim Marks

### Tim Marks (view profile)

Thank you for this very useful tool.
I think there is a small error in the comments at the top of the .m file (what you see if you type "help SpinCalc"). The comments state "Prohibits output if middle angle is within 0.1 degree of singularity value." But in fact, output is prohibited if middle angle is within 2 degrees of singularity value.

Daniel Lucena

### Daniel Lucena (view profile)

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.

Jing

spencer

### spencer (view profile)

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

John Fuller

### John Fuller (view profile)

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

Steve

### Steve (view profile)

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!

Paolo de Leva

### Paolo de Leva (view profile)

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.

Paolo de Leva

### Paolo de Leva (view profile)

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

Paolo de Leva

### Paolo de Leva (view profile)

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

Paolo de Leva

### Paolo de Leva (view profile)

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

Paolo de Leva

### Paolo de Leva (view profile)

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

Paolo de Leva

### Paolo de Leva (view profile)

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

Paolo de Leva

### Paolo de Leva (view profile)

--------

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

Paolo de Leva

### Paolo de Leva (view profile)

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

Paolo de Leva

### Paolo de Leva (view profile)

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

John Fuller

### John Fuller (view profile)

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

Paolo de Leva

### Paolo de Leva (view profile)

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

J. Hu

### J. Hu (view profile)

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

Georg Wiora

James Coburn

zhaopeng QIU

Paul

### Paul (view profile)

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.

Paolo de Leva

### Paolo de Leva (view profile)

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

Paolo de Leva

### Paolo de Leva (view profile)

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

John Fuller

### John Fuller (view profile)

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

Paolo de Leva

### Paolo de Leva (view profile)

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

John Fuller

### John Fuller (view profile)

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

John Fuller

### John Fuller (view profile)

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.

Paolo de Leva

### Paolo de Leva (view profile)

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?

John Fuller

### John Fuller (view profile)

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

Gaetano Giulio Violante

### Gaetano Giulio Violante (view profile)

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]

chandrakala Gowda

Oliver Woodford

Oliver Woodford

### Oliver Woodford (view profile)

Very useful. Thanks.

bryan

### bryan (view profile)

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

James Tursa

### James Tursa (view profile)

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

Joel

Jørgen

### Jørgen (view profile)

Okay, thank you for your feedback.

John Fuller

### John Fuller (view profile)

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.

Jørgen

### Jørgen (view profile)

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.

John Fuller

### John Fuller (view profile)

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.

Stanley

### Stanley (view profile)

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?

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.

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