Is there an algorithm for converting quaternion rotations to Euler angle rotations?

Question

Is there an existing algorithm for converting a quaternion representation of a rotation to an Euler angle representation? The rotation order for the Euler representation is known and can be any of the six permutations (i.e. xyz, xzy, yxz, yzx, zxy, zyx). I've seen algorithms for a fixed rotation order (usually the NASA heading, bank, roll convention) but not for arbitrary rotation order.

Furthermore, because there are multiple Euler angle representations of a single orientation, this result is going to be ambiguous. This is acceptable (because the orientation is still valid, it just may not be the one the user is expecting to see), however it would be even better if there was an algorithm which took rotation limits (i.e. the number of degrees of freedom and the limits on each degree of freedom) into account and yielded the 'most sensible' Euler representation given those constraints.

I have a feeling this problem (or something similar) may exist in the IK or rigid body dynamics domains.

---

Solved: I just realised that it might not be clear that I solved this problem by following Ken Shoemake's algorithms from Graphics Gems. I did answer my own question at the time, but it occurs to me it may not be clear that I did so. See the answer, below, for more detail.

---

Just to clarify - I know how to convert from a quaternion to the so-called 'Tait-Bryan' representation - what I was calling the 'NASA' convention. This is a rotation order (assuming the convention that the 'Z' axis is up) of zxy. I need an algorithm for all rotation orders.

Possibly the solution, then, is to take the zxy order conversion and derive from it five other conversions for the other rotation orders. I guess I was hoping there was a more 'overarching' solution. In any case, I am surprised that I haven't been able to find existing solutions out there.

In addition, and this perhaps should be a separate question altogether, any conversion (assuming a known rotation order, of course) is going to select one Euler representation, but there are in fact many. For example, given a rotation order of yxz, the two representations (0,0,180) and (180,180,0) are equivalent (and would yield the same quaternion). Is there a way to constrain the solution using limits on the degrees of freedom? Like you do in IK and rigid body dynamics? i.e. in the example above if there were only one degree of freedom about the Z axis then the second representation can be disregarded.

---

I have tracked down one paper which could be an algorithm in this pdf but I must confess I find the logic and math a little hard to follow. Surely there are other solutions out there? Is arbitrary rotation order really so rare? Surely every major 3D package that allows skeletal animation together with quaternion interpolation (i.e. Maya, Max, Blender, etc) must have solved exactly this problem?

Answer 1

This looks like a classic case of old technology being overlooked - I managed to dig out a copy of Graphics Gems IV from the garage and it looks like Ken Shoemake has not only an algorithm for converting from Euler angles of arbitrary rotation order, but also answers most of my other questions on the subject. Hooray for books. If only I could vote up Mr. Shoemake's answer and reward him with reputation points.

I guess a recommendation that anybody working with Euler angles should get a copy of Graphics Gems IV from their local library and read the section starting page 222 will have to do. It has to be the clearest and most concise explanation of the problem I have read yet.

---

Here's a useful link I have found since - http://www.cgafaq.info/wiki/Euler_angles_from_matrix - This follows the same system as Shoemake; the 24 different permutations of rotation order are encoded as four separate parameters - inner axis, parity, repetition and frame - which then allows you to reduce the algorithm from 24 cases to 2. Could be a useful wiki in general - I hadn't come across it before.

To old link provided seems to be broken here is another copy of "Computing Euler angles from a rotation matrix ".

Answer 2

In a right-handed Cartesian coordinate system with Z axis pointing up, do this:

struct Quaternion
{
    double w, x, y, z;
};

void GetEulerAngles(Quaternion q, double& yaw, double& pitch, double& roll)
{
    const double w2 = q.w*q.w;
    const double x2 = q.x*q.x;
    const double y2 = q.y*q.y;
    const double z2 = q.z*q.z;
    const double unitLength = w2 + x2 + y2 + z2;    // Normalised == 1, otherwise correction divisor.
    const double abcd = q.w*q.x + q.y*q.z;
    const double eps = 1e-7;    // TODO: pick from your math lib instead of hardcoding.
    const double pi = 3.14159265358979323846;   // TODO: pick from your math lib instead of hardcoding.
    if (abcd > (0.5-eps)*unitLength)
    {
        yaw = 2 * atan2(q.y, q.w);
        pitch = pi;
        roll = 0;
    }
    else if (abcd < (-0.5+eps)*unitLength)
    {
        yaw = -2 * ::atan2(q.y, q.w);
        pitch = -pi;
        roll = 0;
    }
    else
    {
        const double adbc = q.w*q.z - q.x*q.y;
        const double acbd = q.w*q.y - q.x*q.z;
        yaw = ::atan2(2*adbc, 1 - 2*(z2+x2));
        pitch = ::asin(2*abcd/unitLength);
        roll = ::atan2(2*acbd, 1 - 2*(y2+x2));
    }
}

Answer 3

I've been looking for several days for a similar solution, and I finally ran across this website that has an algorithm for converting quaternions to arbitrary Euler and Tait-Bryan rotations!

Here's the link: http://bediyap.com/programming/convert-quaternion-to-euler-rotations/

And here's the code:

///////////////////////////////
// Quaternion to Euler
///////////////////////////////
enum RotSeq{zyx, zyz, zxy, zxz, yxz, yxy, yzx, yzy, xyz, xyx, xzy,xzx};

void twoaxisrot(double r11, double r12, double r21, double r31, double r32, double res[]){
  res[0] = atan2( r11, r12 );
  res[1] = acos ( r21 );
  res[2] = atan2( r31, r32 );
}

void threeaxisrot(double r11, double r12, double r21, double r31, double r32, double res[]){
  res[0] = atan2( r31, r32 );
  res[1] = asin ( r21 );
  res[2] = atan2( r11, r12 );
}

void quaternion2Euler(const Quaternion& q, double res[], RotSeq rotSeq)
{
    switch(rotSeq){
    case zyx:
      threeaxisrot( 2*(q.x*q.y + q.w*q.z),
                     q.w*q.w + q.x*q.x - q.y*q.y - q.z*q.z,
                    -2*(q.x*q.z - q.w*q.y),
                     2*(q.y*q.z + q.w*q.x),
                     q.w*q.w - q.x*q.x - q.y*q.y + q.z*q.z,
                     res);
      break;

    case zyz:
      twoaxisrot( 2*(q.y*q.z - q.w*q.x),
                   2*(q.x*q.z + q.w*q.y),
                   q.w*q.w - q.x*q.x - q.y*q.y + q.z*q.z,
                   2*(q.y*q.z + q.w*q.x),
                  -2*(q.x*q.z - q.w*q.y),
                  res);
      break;

    case zxy:
      threeaxisrot( -2*(q.x*q.y - q.w*q.z),
                      q.w*q.w - q.x*q.x + q.y*q.y - q.z*q.z,
                      2*(q.y*q.z + q.w*q.x),
                     -2*(q.x*q.z - q.w*q.y),
                      q.w*q.w - q.x*q.x - q.y*q.y + q.z*q.z,
                      res);
      break;

    case zxz:
      twoaxisrot( 2*(q.x*q.z + q.w*q.y),
                  -2*(q.y*q.z - q.w*q.x),
                   q.w*q.w - q.x*q.x - q.y*q.y + q.z*q.z,
                   2*(q.x*q.z - q.w*q.y),
                   2*(q.y*q.z + q.w*q.x),
                   res);
      break;

    case yxz:
      threeaxisrot( 2*(q.x*q.z + q.w*q.y),
                     q.w*q.w - q.x*q.x - q.y*q.y + q.z*q.z,
                    -2*(q.y*q.z - q.w*q.x),
                     2*(q.x*q.y + q.w*q.z),
                     q.w*q.w - q.x*q.x + q.y*q.y - q.z*q.z,
                     res);
      break;

    case yxy:
      twoaxisrot( 2*(q.x*q.y - q.w*q.z),
                   2*(q.y*q.z + q.w*q.x),
                   q.w*q.w - q.x*q.x + q.y*q.y - q.z*q.z,
                   2*(q.x*q.y + q.w*q.z),
                  -2*(q.y*q.z - q.w*q.x),
                  res);
      break;

    case yzx:
      threeaxisrot( -2*(q.x*q.z - q.w*q.y),
                      q.w*q.w + q.x*q.x - q.y*q.y - q.z*q.z,
                      2*(q.x*q.y + q.w*q.z),
                     -2*(q.y*q.z - q.w*q.x),
                      q.w*q.w - q.x*q.x + q.y*q.y - q.z*q.z,
                      res);
      break;

    case yzy:
      twoaxisrot( 2*(q.y*q.z + q.w*q.x),
                  -2*(q.x*q.y - q.w*q.z),
                   q.w*q.w - q.x*q.x + q.y*q.y - q.z*q.z,
                   2*(q.y*q.z - q.w*q.x),
                   2*(q.x*q.y + q.w*q.z),
                   res);
      break;

    case xyz:
      threeaxisrot( -2*(q.y*q.z - q.w*q.x),
                    q.w*q.w - q.x*q.x - q.y*q.y + q.z*q.z,
                    2*(q.x*q.z + q.w*q.y),
                   -2*(q.x*q.y - q.w*q.z),
                    q.w*q.w + q.x*q.x - q.y*q.y - q.z*q.z,
                    res);
      break;

    case xyx:
      twoaxisrot( 2*(q.x*q.y + q.w*q.z),
                  -2*(q.x*q.z - q.w*q.y),
                   q.w*q.w + q.x*q.x - q.y*q.y - q.z*q.z,
                   2*(q.x*q.y - q.w*q.z),
                   2*(q.x*q.z + q.w*q.y),
                   res);
      break;

    case xzy:
      threeaxisrot( 2*(q.y*q.z + q.w*q.x),
                     q.w*q.w - q.x*q.x + q.y*q.y - q.z*q.z,
                    -2*(q.x*q.y - q.w*q.z),
                     2*(q.x*q.z + q.w*q.y),
                     q.w*q.w + q.x*q.x - q.y*q.y - q.z*q.z,
                     res);
      break;

    case xzx:
      twoaxisrot( 2*(q.x*q.z - q.w*q.y),
                   2*(q.x*q.y + q.w*q.z),
                   q.w*q.w + q.x*q.x - q.y*q.y - q.z*q.z,
                   2*(q.x*q.z + q.w*q.y),
                  -2*(q.x*q.y - q.w*q.z),
                  res);
      break;
    default:
      std::cout << "Unknown rotation sequence" << std::endl;
      break;
   }
}

Answer 4

I have posted my paper titled "Quaternion to Euler Angle Conversion for Arbitrary Rotation Sequence Using Geometric Methods" on my website at noelhughes.net. I also have algorithms for converting any set of Euler angles to a quaternion and quaternion to/from direction cosine matrix which I will post this weekend. These are also on Martin Bakers website, though a little difficult to find. Google my name, Noel Hughes, and quaternions and you should find it.

Answer 5

I solve it this way:

step 1: Make sure which convention for Euler rotation you want, say, zyx.

step 2: Compute the analytical rotation matrix for the rotation. For example, if you want R(zyx),

**R***zyx* = **R***x*( phi ) * **R***y*( theta ) * **R***z*( psi ), where the elements become

R11 =  cos(theta)*cos(psi)
R12 = -cos(theta)*sin(psi)
R13 =  sin(theta)
R21 =  sin(psi)*cos(phi) + sin(theta)*cos(psi)*sin(phi)
R22 =  cos(psi)*cos(phi) - sin(theta)*sin(psi)*sin(phi)
R23 = -cos(theta)*sin(phi)
R31 =  sin(psi)*sin(phi) - sin(theta)*cos(psi)*cos(phi)
R32 =  cos(psi)sin(phi) + sin(theta)*sin(psi)*cos(phi)
R33 =  cos(theta)*cos(phi) 

step 3: By inspection, you can find the sin or tan for the three angles using the elements above. In this example,

tan(phi) = -R23/R33

sin(theta) = -R13

tan(psi) = -R12/R11

step 4: Compute the rotation matrix from your quaternion (see wikipedia), for the elements you need to compute the angles as in 3) above.

Other conventions can be computed using the same procedure.

Answer 6

Here is a paper I wrote on converting a quaternion to Euler angles.

Link 1

I have also put a number of documents at this location discussing various aspects of quaternions, Euler angles and rotation matrices (DCM).

Link 2

Answer 7

For those that stumble upon this page while Googling, I recently found derivations for these conversions for all 12 intrinsic Tait-Bryan (1-2-3, 3-2-1, etc.) and Proper Euler (1-2-1, 3-1-3, etc.) rotation sequences in the following two references:

• NASA Paper, July 1977: Euler > DCM, Euler > Quaternion, DCM > Euler
• Birdy's Notebook, Oct 2013: Quaternion > Euler

Thanks to frodo2975 for the second link.

Answer 8

Wikipedia shows how you can use the parts of the quaternion and calculate the euler angles.