[Edit 09/29/05 - explained how to extract negative scale]

[Edit 02/26/07 - fixed a math bug. Oops! Thanks, Karsten Isakovic!]

[Edit 06/24/09 - fixed another math bug. Oops! Thanks, Aaron!] [Edit 06/24/09 - Added highlighting to matrix elements so you can see start/stop] [Edit 08/08/10 - Clarified solving for x and z]

While this might almost sound like the title of a bad Zombie movie, it’s all about taking matrices apart. While matrices are a great vehicle to handle transformations (in our case, in 3D graphics), they are a complete pain for humans to read. At least they are for me – whenever I debug 3D code and I look at a matrix, intense pain sets in. The set of values is essentially meaningless.

So what do we need to do? We need to split it apart, into its’ basic components. These would usually be scale, rotation, and translation. Theoretically, there could be skew, but few 3D render engines use that.

If you’re working with a 3D library, you usually got a call that does exactly what I described above. In D3D, its D3DXMatrixDecompose (if your SDK is newer that Summer ’04). But this call has a couple of drawbacks. First, it ties you to a specific sequence of rotations – Yaw,Pitch,Roll for D3D. If your rotations are applied in any other order, the values don’t mean much. And I prefer my rotations in XYZ order, so that’s not an option.

Even worse, though – most decomposers return you a quaternion. And while quaternions are really nifty beasts, they are even less readable than a matrix. Now, granted, you can then convert the quaternion into angles – but at this time, we’ve had a lot of pain already. So, let’s instead find out how to do that decomposition ourselves.

Theoretically, we could take the code for that decomposition from Shoemake, Euler Angle Conversion, in Graphics Gems IV. But first of all, that code is completely unreadable (and gives a great example why macros and short variable names are bad). And second, it doesn’t handle two exceptional cases that can occur when decompositing. So, let’s do it by hand.

Step one: Extract the translation. That’s trivial, really. Just pick the first three values in the last row (or column, if you’re doing column-major matrices)

Step two: The scale. That’s not much more difficult – take the 3×3 submatrix that contains the scale, and compute the “length” of the three rows. Done. (Well, almost. This only accounts for positive scale – but negative scale is done quite easily, too. The determinant is your friend!)

[Edit - how exactly is the determinant your friend? ]

So how do we get if the scale is negative or positive? After all, the length computation always gives a positive. The answer is in the determinant of the matrix. I really can’t explain here how to compute it, or what it is, exactly – but head over to our friends at gamedev.net to get a decent primer on matrices.

Anyways – we know that none, or up to three of our scales are negative. However, scaling negatively on two axes is the equivalent of rotating 180 degrees. (Insert pi/2 in the three basic transformation matrices, and you will see). That means I cannot detect if two scales are negative. Three negative scales is the same as a rotation by 180 degrees and one negative scale.

So all we can do is detect if we should have a negative scale or not. We can’t detect which axis – so we just need to pick one. And the sign to detect this is that the determinant of the normalized rotation matrix is less than 0.

If that’s the case, simply make one of the scales negative. I usually pick X – but really, it doesn’t matter. Mathematically, it’s all the same.

[End Edit]

Step three is where it’s getting interesting. We need to extract the rotation from that same 3×3 submatrix, but the values seem gibberish at a first glance. So let’s see how we arrived at our rotation matrix. Here are the three basic rotation matrices:

X-Rotation:

1 0 0
0 Cos(x) Sin(x)
0-Sin(x)Cos(x)

Y-Rotation:

Cos(y) 0 Sin(y)
010
Sin(y)0Cos(y)

Z-Rotation:

Cos(z) Sin(z) 0
-Sin(z)Cos(z)0
001

We multiply them in the desired order (XYZ) and get:

Cos(y).Cos(z) Cos(y).Sin(z) -Sin(y)
Cos(z).Sin(x).Sin(y) – Sin(z).Cos(x) Sin(x).Sin(y).Sin(z) + Cos(x).Cos(z) Sin(x).Cos(y)
Cos(x).Cos(z).Sin(y) + Sin(x).Sin(z) Cos(x).Sin(y).Sin(z) – Sin(x).Cos(z) Cos(x).Cos(y)

Slightly less mysterious than the data we see in our debugger, but the question remains: How do we get our angles from that? Let’s start easy – if you look at row 1, column 3, you’ll see that value is only dependent on a single angle – Y rotation.

So, easily enough, we get:

rotation.y=asin( - matrix._13 );

Leaves us with x and z.

Unfortunately, there’s not a single matrix cell that is only defined by x or z. So we need to find two cells that both depend on the same variables and are closely enough related that we can actually solve for one of those variables

Staring a bit harder at the matrix, we find a nice pair of values in row 2&3, column 3: Sin(x).Cos(y) and Cos(x).Cos(y).

This yields us with two equations:

(Eq.1) matrix._23 = Sin(x).Cos(y)
 (Eq.2) matrix._33 = Cos(x).Cos(y)

If we divide equation (1) by equation (2), the Cos(y) term cancel each others out. That leaves us with

(Eq.3) matrix._23/matrix._33 = Sin(x)/Cos(x)

Simplified, this gives:

(Eq.3.1) matrix._23/matrix._33 = Tan(x)

So finally, we solve for x:

(Eq.4) atan(matrix._23/matrix/_33) = x

Dividing those two values could of course get an overflow error, but nicely enough there’s a library function that takes care of that, atan2.

Now apply the same idea to row 1, column 1&2, and you’ll get Tan(z), too. So, turning that into code, we get:

rotation.x = atan2( matrix._23, matrix._33 );
rotation.z = atan2( matrix._12, matrix._11 );

But wait – we’re not quite done. We haven’t removed the scale from our matrix before doing this decompostion. We need to divide _11 by scale.x, and _33 by scale.z. Then again, any one of those two scales could be zero. So, no division here – we simply multiply the dividend with the scale instead. I.e.:

rotation.x = atan2( matrix._23scale.x, matrix._33 );
rotation.z = atan2( matrix._12scale.y, matrix._11 );

Are we there yet? Not quite! What if Cos(y) is zero? Then our above equations don’t work out – the division is undefined. Those are the two exceptions where the Gems IV code breaks down. Now, what we need to do in those cases is: We pin one of the two values arbitrarily to 0. (It doesn’t matter what value we pick there, but that makes our life easier. This is the mathematical equivalent of gimbal lock). Let’s set X-rotation to 0. We now know that Cos(x) = +1 and Sin(y) = +/-1. (And of course, Sin(x)=0,Cos(y)=0)

That lets us simplify the nasty equations in the matrix a bit. matrix._22 becomes Cos(z), and matrix._21 becomes -Sin(z). Which, again, gives us Tan(z), but without a division by zero.

All we need to do to detect this special case is see if Sin(y) is less than -0.999 or greater than 0.999 (or whatever margin of error we want to leave).

And finally, we have our rotation. Our matrix is readable by humans.

Commentary

  1. vince wrote on 14. Feb 2005

    Both Bloglines and Firefox freak out with your RSS feed link and can’t parse the link.

  2. NeARAZ wrote on 22. Feb 2005

    I usually tend to think of matrices in terms of coordinate frame. First three rows are your X, Y, Z axes, last one is “position”. It helps a lot, like: need to move “forward” – add something*ZAxis to the “position”. Actually my matrix class has things like Vector3& getAxisX() – note the reference, so you can do stuff like matrix.getOrigin() += matrix.getAxisZ() * velocity; I think over last several years, I didn’t need the “rotation angles” any single time. Oh well, maybe that’s just me :)

  3. Robert Blum wrote on 22. Feb 2005

    I’m working on tools to edit game content – and that means oftentimes somebody wants to see how an object is rotated – hence the need to see the angles.

  4. 3dlamer wrote on 29. Sep 2005

    Cool article except for this part

    > Done. (Well, almost. This only accounts for positive scale – > but negative scale is done quite easily, too. The determinant > is your friend!)

    Everywhere else in your article you actually give the answer but in that one case you assume we know what you’re talking about. I think most readers that know how to apply a determinant to figure out negative scale would not need this article in the first place. Any chance you could add that little piece of info?

  5. Robert Blum wrote on 29. Sep 2005

    Your wish is my command ;) Hope that helps!

  6. Karsten Isakovic wrote on 08. May 2006

    There is a minor bug in your formlar for the matrix.32 – in the first term there should be a Sin(z) instead of Cos(z). So matrix.32 should read Cos(x).Sin(y).Sin(z) – Sin(x).Cos(z) .

  7. Rob wrote on 05. Jul 2006

    Just wonderful to find an xyz rotation matrix for reference.

  8. Bill wrote on 29. May 2008

    Hi, You wrote: “Staring a bit harder at the matrix, we find a nice pair of values in row 2&3, column 3: Sin(x).Cos(y) and Cos(x).Cos(y). If we divide those two, we get Sin(x)/Cos(x) = Tan(x).”

    But the first pair of values is negative and so Tan(x) doesnt work this way….have you simply forgotten the minus sign? Thanks Bill

  9. Administrator wrote on 02. Jun 2008

    Thanks, Bill. Updated the article.

  10. Administrator wrote on 03. Jun 2008

    @Bill: Indeed. Great catch. I’ve updated it

  11. cb wrote on 16. Apr 2009

    Are we there yet? Not quite! What if Cos(y) is zero?

    Is there any way you could add code for this case, as well? Thanks!

  12. Andy Geers wrote on 09. May 2009

    Assuming the description of how to apply the scaling is correct, then I think there’s an error in your code: rotation.x = atan2( matrix._23scale.x, matrix._33 ); rotation.z = atan2( matrix._12scale.y, matrix._11 );

    should really read as follows: rotation.x = atan2( matrix._23scale.z, matrix._33 ); rotation.z = atan2( matrix._12scale.x, matrix._11 );

    … since you said you want to divide matrix._11 by scale.x and matrix._33 by scale.z

  13. Aaron wrote on 19. Jun 2009

    I think there’s a mistake: rotation.y=asin( matrix._13 ) I think should be rotation.y=asin( -matrix._13 )

  14. Rachel 'Groby' Blum wrote on 24. Jun 2009

    Hey Aaron – you are perfectly right. It’s amazing how many typos one can add transcribing math :)

  15. Tim wrote on 08. Aug 2010

    Hey Rachel – this is great. Thank you. The one think I’m stumped on is your statement:

    “We find a nice pair of values in row 2&3, column 3: Sin(x).Cos(y) and Cos(x).Cos(y). If we divide those two, we get Sin(x)/Cos(x) = Tan(x)”

    I don’t understand what allows us to divide m23 by m33. They are just two different values in two cells. How did you determine you could divide?

    I’m sure I’m missing something obvious.

    Thanks

  16. Administrator wrote on 08. Aug 2010

    Hi Tim!

    I hope I’ve clarified a bit how I arrived at this – basically, I’m solving a system with two unknowns for one of them, taking advantage of common terms between them

  17. Tim wrote on 09. Aug 2010

    Thank you. That helps. I kept trying to somehow make _23 and _33 equal each other. What you’re saying makes perfect sense.

    Again, great article!

  18. Vitaly wrote on 12. Apr 2011

    The tricks are great, but math is all wrong from the very beginning. 1) Correct rotation matrices (common rotation direction is assumed): 1 0 0 0 cosx -sinx 0 sinx cosx

    cosy 0 siny 0 1 0 -siny 0 cosy

    cosz -sinz 0 sinz cosz 0 0 0 1

    2) If you rotate around X, then around Y and then around Z, you should multiply Rz by Ry and then by Rx: R = (RzRy)Rx. Correct final rotation matrix is given here http://en.wikipedia.org/wiki/Rotation_matrix#General_rotations phi is x, theta is y and psi is z

  19. Kapil wrote on 19. Apr 2012

    Is it just me but I feel we could have done all the calculations with just the following elements of the above Ratation matrix: _12, _13, _23. step 1) _13 gives Y (the way you described) step 2) _23 gives X (find cos(Y) now that you have Y from step 1, substitute in _23 and get Sin(x) hence X step 3) same as step2 use _12 instead

    This will also avoid all the scaling business.

  20. Huy wrote on 07. Mar 2013

    I agree with Vitaly. The math is all wrong from the beginning. Assuming the scale factor is 1 then it should be:

    y = asin(matrix._31)); z = atan2(-matrix._21, matrix._11); x = atan2(-matrix._32, matrix._33);

  21. Poh wrote on 31. Jul 2013

    Hi, first of all thanks for putting up such a great article. It helped me a lot!

    I am confused with this line for the negative scaling part.

    “If that’s the case, simply make one of the scales negative. I usually pick X – but really, it doesn’t matter. Mathematically, it’s all the same.”

    If the there is one negative scale say on the Y axis like this:

    [ 1 0 0 ] [ 0 -1 0 ] [ 0 0 1 ]

    Using the way you described it would end up as: [-1 0 0 ] [ 0 1 0 ] [ 0 0 1 ]

    How would it be mathematically the same? Is there something that I’m seeing wrongly?

    Thanks a lot for your help!

  22. Bruce McElwee wrote on 20. Sep 2014

    Orthogonal Matrix Angle Decomposition I worked on this problem extensively to come up with the correct input angles for a given matrix. I created the matrix with a Z . Y . X . Y . Z composition from points created using a Z . Y . X matrix. I decomposed the matrix into X, Y, and Z. The problem in the math comes from the inability to determine a precise value for the SIN since -SIN(x) = SIN(-x) and affects the other values of the matrix. The solution I came up with works in all cases. This solution will always accurately recreate the matrix. The solution will be the exact original angles if the original input for Y was such that -90<Y<90. The solution creates the correct matrix for other values of Y from -PI to -90 and 90 to PI, but the angles will not be the original input angles. The solution for Y=-90 or Y=90 does not allow for any derivation of the original X and Z angles from the matrix. In this case, an "if statement" is used to 0 the X value and apply a special value to Z. In a dynamic application, tracking of the X and Z as one of these Y values is approached could allow interpretation of the X and Z quantities.

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

    Hope this helps, it works very well for my application. Bruce

Leave a reply