[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.
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:
We multiply them in the desired order (XYZ) and get:
|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.