Maths - Conversion Quaternion to Euler

/** assumes q1 is a normalised quaternion */

public void set(Quat4d q1) {
	double test = q1.x*q1.y + q1.z*q1.w;
	if (test > 0.499) { // singularity at north pole
		heading = 2 * atan2(q1.x,q1.w);
		attitude = Math.PI/2;
		bank = 0;
		return;
	}
	if (test < -0.499) { // singularity at south pole
		heading = -2 * atan2(q1.x,q1.w);
		attitude = - Math.PI/2;
		bank = 0;
		return;
	}
    double sqx = q1.x*q1.x;
    double sqy = q1.y*q1.y;
    double sqz = q1.z*q1.z;
    heading = atan2(2*q1.y*q1.w-2*q1.x*q1.z , 1 - 2*sqy - 2*sqz);
	attitude = asin(2*test);
	bank = atan2(2*q1.x*q1.w-2*q1.y*q1.z , 1 - 2*sqx - 2*sqz)
}

Note1: I have replaced Math.atan(a/b) with Math.atan2(a,b) as suggested by Mark Elson

Note2: The cutoff point for singulaties is set to 0.499 and -0.499 which corresponds to 86.3 degrees, this can be set closer to 90 if required but we have to accept some inaccuracy around the singulaties. (Thanks to Jack Morrison for corrections)

Note3: as suggested by minorlogic and Andy the above expressions can be made to work with non-normalised quaternion (this saves a CPU heavy square root operation). This uses 1=sqx + sqy + sqz + sqw so:

/** q1 can be non-normalised quaternion */

public void set(Quat4d q1) {
    double sqw = q1.w*q1.w;
    double sqx = q1.x*q1.x;
    double sqy = q1.y*q1.y;
    double sqz = q1.z*q1.z;
	double unit = sqx + sqy + sqz + sqw; // if normalised is one, otherwise is correction factor
	double test = q1.x*q1.y + q1.z*q1.w;
	if (test > 0.499*unit) { // singularity at north pole
		heading = 2 * atan2(q1.x,q1.w);
		attitude = Math.PI/2;
		bank = 0;
		return;
	}
	if (test < -0.499*unit) { // singularity at south pole
		heading = -2 * atan2(q1.x,q1.w);
		attitude = -Math.PI/2;
		bank = 0;
		return;
	}
    heading = atan2(2*q1.y*q1.w-2*q1.x*q1.z , sqx - sqy - sqz + sqw);
	attitude = asin(2*test/unit);
	bank = atan2(2*q1.x*q1.w-2*q1.y*q1.z , -sqx + sqy - sqz + sqw)
}

In order to try to get an intuative understanding of the singularities involved in converting other representations of 3D rotations to Euler angles it may help to look at the way we project the surface of a sphere onto a 2 dimensional map. This is a different case and it only involves two angles, latitude and longitude, but this simpler example may show the principles. This analogy is explained in more detail on this page.

At the north and south poles longitude does not matter so if we are converting some other description of north pole to latitude, longitude it might give any value for longitude, for example, infinity. So we have to be very careful at these points.

Similarly we can map Euler angles to quaternions (4 dimensional hypersphere). This maps a one dimensional space (rotations around 0,1,0 axis) to a two dimensional plane in Euler terms. This is where attitude = 90° and heading, bank vary:

On this plane lines of common orientation are diagonal lines, that is rotation around 0,1,0 axis are given by angle = heading+bank.

Similarly for the south pole.

Example

we take the 90 degree rotation from this:

to this:

As shown here the quaternion for this rotation is: (0.7071+ i 0.7071)

So using the above result:

sqw = q1.w*q1.w = 0.5
sqx = q1.x*q1.x = 0.5
sqy = q1.y*q1.y = 0
sqz = q1.z*q1.z =0
heading = atan2(2.0 * (q1.x*q1.y + q1.z*q1.w),(sqx - sqy - sqz + sqw)) = atan2(0,2) = 0
bank = atan2(2.0 * (q1.y*q1.z + q1.x*q1.w),(-sqx - sqy + sqz + sqw)) = atan2(0.5,0) = 90 degrees
attitude = asin(-2.0 * (q1.x*q1.z - q1.y*q1.w)/sqx + sqy + sqz + sqw) = asin(0/2) = 0

So this gives the correct result shown here, it is banking by 90 degrees, but we have to be very careful about the following issues:

• If we had not used the atan2 function then one of the intermediate steps would involve infinity (division by 0) therefore most maths libraries have an atan2 function to avoid this step.
• Most maths libraries use radians instead of degrees (apart from OpenGL).

Alternative Euler Angle Sequences

I am trying to make this site consistant and therefore I have standarised on one Euler angle sequence as defined on this page. However there are often reasons to use different Euler angle sequences:

• You may be working with a program which uses a different Euler angle sequence.
• You may need the singularities to be at different angles so they are less likely to cause problems.
• The objects that you are modelling may just favor a different set of angles.

Amy de Buitléir has written this document and kindly allowed me to publish it here. This shows how to convert a quaternion to any Euler angle sequence. Instead of defining the quaternion in terms of rotations about the absolute coordinates i, j and k the document defines 3 mutually perpendicular axes e1, e2 and e3. So, to generate the mapping for a given set of Euler angles the user needs to map e1, e2 and e3 to i, j or k in the appropriate order. This involves the value e which seems to define a sort of left or right handedness, but on its own this is not enough to define the Euler angle sequence, will still need to define the e1, e2, e3 to i, j, k mapping.

Angle Calculator and Further examples

I have put a java applet here which allows the values to be entered and the converted values shown along with a graphical representation of the orientation.

Also further examples in 90 degree steps here

Or enter the quaternion below, then click calculate to get the heading, attitude and bank in degrees.