Decomposing and composing a 3×3 rotation matrix

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]

The 3 Euler angles are

\theta_{x} = atan2\left(r_{32}, r_{33}\right)

\theta_{y} = atan2\left(-r_{31}, \sqrt{r_{32}^2 + r_{33}^2}\right)

\theta_{z} = atan2\left(r_{21}, r_{11}\right)

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 \theta_{x}, \theta_{y}, \theta_{z} , the rotation matrix is calculated as follows:

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]

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]

R = ZYX

Note on angle ranges

The Euler angles returned when doing a decomposition will be in the following ranges:

\theta_{x} \rightarrow \left(-\pi, \pi \right)

\theta_{y} \rightarrow \left(-\frac{pi}{2}, \frac{pi}{2} \right)

\theta_{z} \rightarrow \left(-\pi, \pi \right)

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!

23 thoughts on “Decomposing and composing a 3×3 rotation matrix

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

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

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

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

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

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

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

  8. 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;;/

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

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

Leave a Reply

Your email address will not be published. Required fields are marked *

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>