It solved two problems of the trivial conversion method in Matrix and Quaternion FAQ:

- Use arbitrarily ordered axes (like ZXY) in addition to XYZ
- Remove the discontinuity happening at Ax = +/- 90 degrees

- The algorithm always return the minimum set of absolute Euler angles (closest to {0;0;0}, as defined by the sum of squares of the angles) from the two alternatives. Often it is desirable to get the set of angles that are closest not to identity {0;0;0} but to some predefined rotation {Rx;Ry;Rz}.
- When gimbal lock occured, the X rotation (or corresponding rotation if axes were reordered) was set to zero. This may not be the best solution, if we need the Euler angles to be as close to some specific set of values as possible.

Stable conversio of matrix to Euler angles closest to given set

The matrices have row-major order, so they have to be transposed for OpenGL.

**Base**is a set of three angles {Bx;By;Bz} that will be used as reference while determining the best set of resulting angles.

void Matrix3x4f::getEulerAngles (f32 angles[], const f32 base[]) const { float aX, aY, aZ; float aX2, aY2, aZ2; aY = asin (c[2]); aY2 = (aY >= 0) ? M_PI_F - aY : -M_PI_F - aY; if (fabs (c[2]) < 0.99999f) { // Over 0.25 degrees from right angle float C = cos (aY); float C2 = cos (aY2); float trX = c[10] / C; float trY = -c[6] / C; aX = atan2 (trY, trX); aX2 = atan2 (-c[6] / C2, c[10] / C2); trX = c[0] / C; trY = -c[1] / C; aZ = atan2 (trY, trX); aZ2 = atan2 (-c[1] / C2, c[0] / C2); } else { // Less than 0.25 degrees from right angle // Gimbal lock occured aX = aX2 = base[0]; // Keep X rotation the same as base float Sx = sin (aX); float Cx = cos (aX); if (fabs (Cx) > fabs (Sx)) { float Sz = (c[4] + (Sx/Cx) * c[8]) / (Sx*Sx/Cx + Cx); aZ = aZ2 = asin (Sz); } else { float Sz = (c[8] + (Cx/Sx) * c[4]) / (Cx*Cx/Sx + Sx); aZ = aZ2 = asin (Sz); } } float d1X = clamp_pi (aX - base[0]); float d1Y = clamp_pi (aY - base[1]); float d1Z = clamp_pi (aZ - base[2]); float d2X = clamp_pi (aX2 - base[0]); float d2Y = clamp_pi (aY2 - base[1]); float d2Z = clamp_pi (aZ2 - base[2]); if ((d1X * d1X + d1Y * d1Y + d1Z * d1Z) < (d2X * d2X + d2Y * d2Y + d2Z * d2Z)) { angles[X] = aX; angles[Y] = aY; angles[Z] = aZ; } else { angles[X] = aX2; angles[Y] = aY2; angles[Z] = aZ2; } }

If aY is close to +/- 90 degrees (gimbal lock occured) we want to set aX to the base value instead of setting it to zero.

Also we calculate the differences of both alternate sets of angles from base and choose the final set based on the squared sum of those differences.

**Clamp_pi**simply converts angle to -PI...PI range.

Conversion of matrix to arbitrarily ordered Eular angles

The matrices have row-major order, so they have to be transposed for OpenGL.

**Order**is array of three integers, determining the order of rotations ({0;1;2} is XYZ, {2;0;1} is ZXY).

**Base**values should be ordered according to order.

void Matrix3x4f::getEulerAngles (f32 angles[], const int order[], const f32 base[]) const { float sign = ((order[1] == (order[0] + 1)) || (order[2] == (order[1] + 1))) ? 1.0f : -1.0f; // Permute matrix to switched axes space Elea::Matrix3x4f rs; for (int row = 0; row < 3; row++) { for (int col = 0; col < 3; col++) { rs[4 * row + col] = c[4 * order[row] + order[col]]; } } rs[2] *= sign; rs[6] *= sign; rs[8] *= sign; rs[9] *= sign; // New base float nbase[3]; for (int i = 0; i < 3; i++) nbase[i] = base[i]; nbase[2] *= sign; // Decompose to Euler angles rs.getEulerAngles (angles, nbase); angles[2] *= sign; }

The trick is to permute the rotation matrix in such way, that the XYZ angles of new matrix are properly ordered decomposition angles of the original matrix and then do the trivial decomposition.

Hopefully this time the algorithm is final ;-)

Have fun!