This post shows how to decompose a 3×3 rotation matrix into the 3 elementary Euler angles, sometimes referred to as yaw/pitch/roll, and going the other way around. The technique I’m presenting is based off http://planning.cs.uiuc.edu/node102.html.
If you have ever seen Wikipedia’s entry on Rotation matrix or Euler angles, you will have no doubt been swamped to your neck with maths equations all over the place, depending how tall your neck is. It turns out there is no single correct answer on defining a rotation matrix in terms of Euler angles. There are a few ways to accomplish it and all of them are valid. But that’s okay, I’ll just show one way, which should be adequate for most applications. I won’t go into any maths derivation, aiming to keep this post implementation friendly.
Decomposing a rotation matrix
Given a 3×3 rotation matrix
![R = \left[ \begin{array}{ccc} r_{11} & r_{12} & r_{13} \\ r_{21} & r_{22} & r_{23} \\ r_{31} & r_{32} & r_{33} \end{array} \right]](https://s0.wp.com/latex.php?latex=R+%3D+%5Cleft%5B+%5Cbegin%7Barray%7D%7Bccc%7D+r_%7B11%7D+%26+r_%7B12%7D+%26+r_%7B13%7D+%5C%5C+r_%7B21%7D+%26+r_%7B22%7D+%26+r_%7B23%7D+%5C%5C+r_%7B31%7D+%26+r_%7B32%7D+%26+r_%7B33%7D+%5Cend%7Barray%7D+%5Cright%5D+&bg=ffffff&fg=000000&s=0 "R = \left[ \begin{array}{ccc} r_{11} & r_{12} & r_{13} \\ r_{21} & r_{22} & r_{23} \\ r_{31} & r_{32} & r_{33} \end{array} \right]")
The 3 Euler angles are
Here atan2 is the same arc tangent function, with quadrant checking, you typically find in C or Matlab.
Composing a rotation matrix
Given 3 Euler angles 
![X = \left[ \begin{array}{ccc} 1 & 0 & 0 \\ 0 & \cos\left(\theta_{x}\right) & -\sin\left(\theta_{x}\right) \\ 0 & \sin\left(\theta_x\right) & \cos\left(\theta_{x}\right) \end{array} \right]](https://s0.wp.com/latex.php?latex=X+%3D+%5Cleft%5B+%5Cbegin%7Barray%7D%7Bccc%7D+1+%26+0+%26+0+%5C%5C+0+%26+%5Ccos%5Cleft%28%5Ctheta_%7Bx%7D%5Cright%29+%26+-%5Csin%5Cleft%28%5Ctheta_%7Bx%7D%5Cright%29+%5C%5C+0+%26+%5Csin%5Cleft%28%5Ctheta_x%5Cright%29+%26+%5Ccos%5Cleft%28%5Ctheta_%7Bx%7D%5Cright%29+%5Cend%7Barray%7D+%5Cright%5D+&bg=ffffff&fg=000000&s=0 "X = \left[ \begin{array}{ccc} 1 & 0 & 0 \\ 0 & \cos\left(\theta_{x}\right) & -\sin\left(\theta_{x}\right) \\ 0 & \sin\left(\theta_x\right) & \cos\left(\theta_{x}\right) \end{array} \right]")
![Y = \left[ \begin{array}{ccc} \cos\left(\theta_{y}\right) & 0 & \sin\left(\theta_{y}\right) \\ 0 & 1 & 0 \\ -\sin\left(\theta_{y}\right) & 0 &\cos\left(\theta_{y}\right) \end{array} \right]](https://s0.wp.com/latex.php?latex=Y+%3D+%5Cleft%5B+%5Cbegin%7Barray%7D%7Bccc%7D+%5Ccos%5Cleft%28%5Ctheta_%7By%7D%5Cright%29+%26+0+%26+%5Csin%5Cleft%28%5Ctheta_%7By%7D%5Cright%29+%5C%5C+0+%26+1+%26+0+%5C%5C+-%5Csin%5Cleft%28%5Ctheta_%7By%7D%5Cright%29+%26+0+%26%5Ccos%5Cleft%28%5Ctheta_%7By%7D%5Cright%29+%5Cend%7Barray%7D+%5Cright%5D+&bg=ffffff&fg=000000&s=0 "Y = \left[ \begin{array}{ccc} \cos\left(\theta_{y}\right) & 0 & \sin\left(\theta_{y}\right) \\ 0 & 1 & 0 \\ -\sin\left(\theta_{y}\right) & 0 &\cos\left(\theta_{y}\right) \end{array} \right]")
![Z = \left[ \begin{array}{ccc} \cos\left(\theta_{z}\right) & -\sin\left(\theta_{z}\right) & 0 \\ \sin\left(\theta_{z}\right) & \cos\left(\theta_{z}\right) & 0 \\ 0 & 0 & 1 \end{array} \right]](https://s0.wp.com/latex.php?latex=Z+%3D+%5Cleft%5B+%5Cbegin%7Barray%7D%7Bccc%7D+%5Ccos%5Cleft%28%5Ctheta_%7Bz%7D%5Cright%29+%26+-%5Csin%5Cleft%28%5Ctheta_%7Bz%7D%5Cright%29+%26+0+%5C%5C+%5Csin%5Cleft%28%5Ctheta_%7Bz%7D%5Cright%29+%26+%5Ccos%5Cleft%28%5Ctheta_%7Bz%7D%5Cright%29+%26+0+%5C%5C+0+%26+0+%26+1+%5Cend%7Barray%7D+%5Cright%5D+&bg=ffffff&fg=000000&s=0 "Z = \left[ \begin{array}{ccc} \cos\left(\theta_{z}\right) & -\sin\left(\theta_{z}\right) & 0 \\ \sin\left(\theta_{z}\right) & \cos\left(\theta_{z}\right) & 0 \\ 0 & 0 & 1 \end{array} \right]")
Note on angle ranges
The Euler angles returned when doing a decomposition will be in the following ranges:
If you keep your angles within these ranges, then you will get the same angles on decomposition. Conversely, if your angles are outside these ranges you will still get the correct rotation matrix, but the decomposed values will be different to your original angles.
Code
Download rotation_matrix_demo.m
The Octave/Matlab script contains the decompose/compose function and a demo on using it. It picks random Euler angles, makes a rotation matrix, decomposes it and verifies the results are the same. An example output
octave:1> rotation_matrix_demo Picking random Euler angles (radians) x = -2.6337 y = -0.47158 z = -1.2795 Rotation matrix is: R = 0.25581 -0.77351 0.57986 -0.85333 -0.46255 -0.24057 0.45429 -0.43327 -0.77839 Decomposing R x2 = -2.6337 y2 = -0.47158 z2 = -1.2795 err = 0 Results are correct!
Hi Nghia,
I believe you misplaced number 1 in the Y rotation matrix. Number 1 should be at the centre of Y. As a result, your derivation of the angles is not generally correct. From the RHS of equation 3.42 of the link you cited (http://planning.cs.uiuc.edu/node102.html), one can easily obtain yaw, pitch and roll angles as follows:
\theta_x = arctan(r_{3,2}/r_{3,3})
\theta_y = -arcsin(r_{3,1})
\theta_z = arctan(r_{2,1}/r_{1,1})
By the way, you other posts and source codes on computer vision are amazing. They are really useful to me.
Cheers,
Chuong
Hi,
Well spotted!
I came to the same result for theta_y, using arcsin instead of atan2. I wasn’t sure if there was any ‘gotchas’ at the time, so I left it using the original equation I cited.
Glad you found something useful on the site 🙂
Hi Nghia,
I think this approach has one singularity. If theta_y is +-90 degrees, it will come that r_{1,1}, r_{2,1}, r_{3,2} and r_{3,3} will evaluate to zero due to the cos(theta_y) term. So you end up to atan(0) which, will return always 0, independent of theta_x and theta_z.
If you try this in your matlab code, make sure that you do not have numerical errors on r_{1,1}, r_{2,1}, r_{3,2} and r_{3,3}, after generating the rotation matrix.
Do you have the solution for that case also?
Cheers,
Andre
I’ve tried putting 90 degrees in my demo program, at line 6
y = pi;
x,z are still random between -180,180 degrees. The matrix changes accordingly. The returned values however in the x,z axis are out by 180 degrees.
Try editing the demo code and see if you get the same results as me.
I’ve tried to put in line 6:
y = pi*0.5;
and after the line “R = compose_rotation(x,y,z)”, the lines:
R(1,1) = 0;
R(2,1) = 0;
R(3,2) = 0;
R(3,3) = 0
to kill numerical errors of the composition process. If you read these values before, you will get values very close to zero (e-17 in my machine). They actually should be zero, due to the “cos(pi*0.5)” term in that specific case.
You can let x and z random.
Ooops my bad, wasn’t thinking straight. You’re very right.
I guess the only one around it is to explicitly check for the +- 90 degree situation and use a different rotation composition. Probably one of the many reasons why people hate using matrices to represent rotation.
Indeed. By the way, thank you for answering and congrats for the blog.
Hi Nghia!
Thank you very much for your article. It helped me very much. For about 2 days I couldn’t solve my problem, I’ve messed up in these rotation-conventions. But when I’ve used your formulas, it worked fine!
It’s important to note that the Euler angles described above are applied in the order theta_X -> theta_Y -> theta_Z relative to the world axis, due to the non-commutativity of 3D rotation.
thank you nghiaho, you excelent, i see much your paper. your’s result very good.
Hello: I have a visualization program which expects 3 values for rotation (x, y, z) between 0 and 360 degrees. Now, my fusion library (reading data from sensors) provides either euler angles or quaternions. I didin’t know how to use quaternions to feed the visualization application, so I ‘m working with euler angles.
The problem is I see euler angles values range from (0 to 360), (-180 to 180), and (-90 to 90). For the first case is Ok. For the second, I can add 180 and I’m done (range 0 to 360). Now for the last one, I see if doing a 360 movement with the board (sensors), the values are changing from -90 to +90 and +90 to -90 again, so in this case I’m not sure what to do. Adding 90 will end in (0 to 180) range.
How can I convert the euler angles ranges to the right (0 to 360) value range the visualization application requires?
Many thanks in advance.
Gus
If you want to convert all your angles to [0,360] correctly then I would simply do:
if(angle < 0) angle += 360;
You are aware I hope, that what you’re describing here are actually the Tait-Bryan angles (or the xyz convention) not the Euler angles (cf. Goldstein p. 154 and appendix)? It really doen’t matter which one you use as long as you are precise about which definition you use (which you are in a sense since you give a link to the maths calculations).
Thanks for the clarification! I didn’t actually know that, my maths background is a bit limited 🙂
Hey Nghia,
i also wanna thank you for your helpful tutorial.
I am working with opencv and visp.
I tried to take an angle alpha in radians, compose the matrix for pitch roll and yaw
then multiplied them like Z*Y*X and then tried to decompose the matrix again to
get the angles for x,y and z. After that i compared the angles as degrees.
It seems that the new calculated angles always are 180°-alpha.
What am i doing wrong?
Here is my sourcecode:
#include
#include
#include
#include
#include
#include
#include
using namespace std;
using namespace cv;
int main() {
for (double angleOld = -M_PI; angleOld <= M_PI; angleOld+=(2*M_PI)/360) {
double a = angleOld;
// ——————————————————————-
// composing rotation matrix
// ——————————————————————-
Mat xRotationOld = *(Mat_(3, 3) << 1, 0, 0,
0, cos(a), -sin(a),
0, sin(a), cos(a) );
Mat yRotationOld = *(Mat_(3, 3) << cos(a), 0, sin(a),
0, 1, 0,
-sin(a),0, cos(a) );
Mat zRotationOld = *(Mat_(3, 3) << cos(a), -sin(a),0,
sin(a), cos(a), 0,
0, 0, 1 );
Mat rotationMat = zRotationOld * yRotationOld * xRotationOld;
// ——————————————————————-
// decomposing rotation matrix
// ——————————————————————-
/*
rotationMat:
(m_00 m_01 m_02)
(m_10 m_11 m_12)
(m_20 m_21 m_22)
*/
double m_00 = rotationMat.at(0, 0);
double m_10 = rotationMat.at(1, 0);
double m_20 = rotationMat.at(2, 0);
double m_21 = rotationMat.at(2, 1);
double m_22 = rotationMat.at(2, 2);
// calculate the euler angles
double xRotation = atan2(m_21, m_22); // range: (-pi, pi)
double yRotation = atan2(- m_20, sqrt((m_21*m_21) + (m_22*m_22))); // range: (-(pi/2), pi/2)
double zRotation = atan2(m_10, m_00); // range: (-pi, pi)
// convert radians to degree
double oldAngleDegree = (angleOld * 180) / M_PI;
double xDegree = (xRotation * 180) / M_PI;
double yDegree = (yRotation * 180) / M_PI;
double zDegree = (zRotation * 180) / M_PI;
// check wether angles are the same as before
if (xRotation != angleOld || yRotation != angleOld || zRotation != angleOld) {
cout << oldAngleDegree << "° wrong calculated " << endl;
cout << "x angle = ";
cout << xDegree << "° :: ";
if (xRotation != angleOld) {
cout << "wrong" << endl;
}else {
cout << "correct" << endl;
}
cout << "y angle = ";
cout << yDegree << "° :: ";
if (yRotation != angleOld) {
cout << "wrong" << endl;
}else {
cout << "correct" << endl;
}
cout << "z angle = ";
cout << zDegree << "° :: ";
if (zRotation != angleOld) {
cout << "wrong" << endl;
}else {
cout << "correct" << endl;
}
}
}
}
I hope anyone can help me with that. i guess it is some simple thing that i am overseeing.
Much thx in advance 🙂
Dear Nghiaho,
I tested your function to convert to-and-fro a rotation matrix, however the result did not match.
R =
0.9918 -0.0025 -0.1281
0.0016 1.0000 -0.0070
-0.1281 -0.0067 -0.9917
[x2,y2,z2] = decompose_rotation(R)
x2 =
-3.1348
y2 =
0.1285
z2 =
0.0016
R2= compose_rotation(x2,y2,z2)
R2 =
0.9918 0.0008 -0.1281
0.0016 -1.0000 0.0066
-0.1281 -0.0067 -0.9917
Any idea how I can reconstruct the exact rotation matrix after I transform it to Euler angles?
Thanks and best regards,
Raymond
Your R is not a numerically correct rotation matrix. If you do an SVD on it your singular values are slightly off 1, which means inv(R) = R’ does not hold either. Your determinant is -1, so there’s a reflection which I believe you have to correct as we;;/
I found this formula to duplicate every matrix I could randomly generate. When Y = +-90 degrees, the difference beween the X and Z angles is all that can be determined and this is used as a Z angle setting X to 0.
X = if(or(m31=1,m31=-1),0,atan2(m33+1e-24,m32))
Y= -asin(m31)
Z= if(or(m31=1,m31=-1),-atan2(m22,m12),atan2(m11+1e-24,m21))
Hope this works for you!
I noticed that the addition of the 1e-24 in the two terms is probably not needed. My original thinking was to keep the term from being 0! These equations will decompose to the exact angles for -90<Y<90. In the quadrant 3 and 4 ranges for Y, X= PI-X and Z=PI-Z. However, due to sin(-a) = -sin(a), this condition cannot be determined from the matrix. Therefore there are always two possible solutions except for Y=+-PI, which is impossible to determine the X and Z components seperately from the matrix. In a flowing transition, the limit as Y=+-PI is approached could be used to preserve the X and Z angles!
Sorry for the rambling, but I need to get this right! The second solution from the above equations are
X = -sign(X)*PI-X
Y= -sign(Y)*PI-Y
Z= -sign(Z)*PI-Z
Again however, there is no wy to determine which solution was initially nput to produce the matrix. Nonetheless, both are valid solutions to produce the matrix.
Have a great day
Bruce
I will get this right soon :o)
solution 1
X = if(or(m31=1,m31=-1),0,atan2(m33,m32))
Y= -asin(m31)
Z= if(or(m31=1,m31=-1),-atan2(m22,m12),atan2(m11,m21))
solution 2
X’ = -sign(X)*PI+X
Y’ = sign(Y)*PI-Y
Z’ = -sign(Z)*PI+Z
I am interpreting this from examples running on an excel program. I assure you that it does work !!!!
i am using accelerometer and magnetometer readings to get the heading.
I am using raw readings from mobile sensors, body coordinates and world coordinates are aligned properly, using ENU convention, but i am not able to get the heading using these rotation matrices.
Please help me in writing matlab code to get the heading, pitch and roll are given by accelerometer readings and yaw will be given by using both.
will be highly thankful .
I haven’t done any work with IMU so can’t help you there. I’m sure there’s heaps of sample code out there though.
what about scaling? its should be considered as well, am i wrong?
ICP was intended for 3D points of the same scale eg. from laser scanner
You are right – this was way more reader friendly than Wikipedia. Very well organized. Thank you.
Hi,
I haven’t been getting the expected rotation angles by decomposing the R matrix.
In the process of finding an independent calculation, I tried MS Excel ATAN2() function.
Using the example above, ThetaX = atan2(-0.43327,-0.77839) results in -2.07871, in comparison to -2.633679, which you get and I also get with my .NET math.atan2() code. Have you seen this before with Excel? Am I missing something ?
Anyway, I am trying to validate my code which is based on your registration approach of two 3D data sets. I take two identical sets of data, rotate one set by a known rotation (about one axis only) and I should see that reflected in the decomposition of the R matrix. It works for rot about X, rot about Y, but when I try for a rot about Z, the decomposition shows rotations about all axis and none sensible and a translation. With the X and Y examples I get the correct rotation and no translation.
Sorry for the long post !
The Z rotation data is:
Ref Points:
P1, 230.5033, 112.0004, 111.0795
P2, 174.0455, 133.189, 91.36805
P3, 204.2584, 151.8481, 80.76611
Rot. Points (=Ref.Points rotated by 15deg about Z axis)
P1, 193.6612539, 167.8427229, 111.0795
P2, 133.6431936, 173.696985, 91.36805
P3, 157.9972836, 199.5399655, 80.76611
The Reg matrices are:
Rotation Matrix
0.936, 0.333, 0.112
-0.011, 0.348, -0.937
0.351, -0.876, -0.33
Translation Matrix
-19.172
159.884
226.769
Registration Error
0.007432663
When decomposing R I get
ThetaX = -1.931067191(radians)
ThetaY = -0.358734974(radians)
ThetaZ = -0.011746066(radians)
Which is nothing like the 15deg Z rotation I used…
Any help/suggestions would be greatly appreciated.
(Frustrated PhD student here!)
The determinant of your rotation matrix is -1. This is generally undesirable, and means there is some unwanted reflection going on.
Fantastic explanation. Thank you.
Great collection of useful little functions! It’s the second time i come find your code+explanations.
You should share on MatlabCentral!
Hi,
I’ve tried using the code to split up a rotation and deformation (like stretch) from the H matrix. The code works perfectly when I test it using a triangle which just rotates, but as soon as it stretches and rotates then the rotation matrix is wrong.
Also, for any rotation in the z axis, even without stretching, I receive wrong values.
e.g. initial=
[ 2 0 2
2 2 2
0 0 2]
final=
[2 0 2
2 2 2
2 0 0]
Sorry, the above comment was meant for the ‘FINDING OPTIMAL ROTATION AND TRANSLATION BETWEEN CORRESPONDING 3D POINTS’ post.
This code only works for a pure rotation matrix. It won’t work if there is scaling or shearing involved. You got an example of rotation around the zaxis going wrong?
I am a master student in Incheon, Korea, studying biomedical engineering and I’ve been using your ARToolKit for a long time.
What I want to ask you is the detail info about ‘arTransMat3.c’,
especially the matrix to make a 3-angles. I thought this is a typical euler matrix but it did not match to any sort of euler matrix formula. May I get a comment with this issue? Below is the matrix generating code in ‘arTransMat3.c’.
I did not write ARToolKit or have anything to do with it.
Hello,
In your section titled “Composing a rotation matrix” you show the different rotation matrixs for X Y Z.
I understand that if I change the signs of the sine functions the rotations will go in the opposite direction.
Do you know if I did this how it would effect your decomposition formulas in section “Decomposing a rotation matrix” “The 3 Euler angles are”? if at all?
I have a situation where I am given numbers for a matrix but I know it rotates in opposite direction and I dont know if your formulas for decomposition would work for me. Any suggestion for me on how to decompose a matrix that rotates in opposite direction? meaning the signs with the sine are opposite of what you have. thanks
This could be a different ordering of the Euler angles. They have something like up to 24 different combinations.
Thank you for posting this formula. Here it is coded in Python for use in Blender Game Engine.
from math import atan2, pow, sqrt
orientation = gameObject.worldOrientation
xyzRotations = [
atan2(orientation[2][1], orientation[2][2]),
atan2(orientation[2][0] * -1.0
, sqrt(pow(orientation[2][1], 2)
+ pow(orientation[2][2], 2))),
atan2(orientation[1][0], orientation[0][0])
]
Cheers, Jim
Hello, how can I get the Euler Angles from two vectors, A (initial) and B (final), without the rotation matrix.
I’m guessing you first find a rotation matrix that rotates A to B then decompose that.
Hi, thanks for this article, the reference I found on wikipedia was either in another coordinate system or simply wrong… while your article fixed my problem