Maths - Conversion Matrix to Axis Angle

Summary of results

A 3×3 matrix has 9 numbers therefore it contains replicated information, so there are many ways to derive the rotation from the numbers, here is a possible conversion:

angle = acos(( m00 + m11 + m22 - 1)/2)
x = (m21 - m12)/√((m21 - m12)2+(m02 - m20)2+(m10 - m01)2)
y = (m02 - m20)/√((m21 - m12)2+(m02 - m20)2+(m10 - m01)2)
z = (m10 - m01)/√((m21 - m12)2+(m02 - m20)2+(m10 - m01)2)

There are two singularities at angle = 0° and angle = 180°, in these cases the above formula may not work (as pointed out by David) so we have to test for these cases separately. At 0° the axis is arbitrary (any axis will produce the same result), at 180° the axis is still relevant so we have to calculate it.

Derivation of Equations

I can think of two methods to derive the equations:

I have used both of these and so that I don't interrupt the flow I have put them on this page.

Java Code

Here is java code to do conversion, the handling of the singularities at 0° and 180° is explained in the next section below.
/**
This requires a pure rotation matrix 'm' as input.
*/
public axisAngle toAxisAngle(matrix m) {
  double angle,x,y,z; // variables for result
	double epsilon = 0.01; // margin to allow for rounding errors
	double epsilon2 = 0.1; // margin to distinguish between 0 and 180 degrees
	// optional check that input is pure rotation, 'isRotationMatrix' is defined at:
	// http://www.euclideanspace.com/maths/algebra/matrix/orthogonal/rotation/
	assert isRotationMatrix(m) : "not valid rotation matrix" ;// for debugging
	if ((Math.abs(m[0][1]-m[1][0])< epsilon)
	  && (Math.abs(m[0][2]-m[2][0])< epsilon)
	  && (Math.abs(m[1][2]-m[2][1])< epsilon)) {
		// singularity found
		// first check for identity matrix which must have +1 for all terms
		//  in leading diagonaland zero in other terms
		if ((Math.abs(m[0][1]+m[1][0]) < epsilon2)
		  && (Math.abs(m[0][2]+m[2][0]) < epsilon2)
		  && (Math.abs(m[1][2]+m[2][1]) < epsilon2)
		  && (Math.abs(m[0][0]+m[1][1]+m[2][2]-3) < epsilon2)) {
			// this singularity is identity matrix so angle = 0
			return new axisAngle(0,1,0,0); // zero angle, arbitrary axis
		}
		// otherwise this singularity is angle = 180
		angle = Math.PI;
		double xx = (m[0][0]+1)/2;
		double yy = (m[1][1]+1)/2;
		double zz = (m[2][2]+1)/2;
		double xy = (m[0][1]+m[1][0])/4;
		double xz = (m[0][2]+m[2][0])/4;
		double yz = (m[1][2]+m[2][1])/4;
		if ((xx > yy) && (xx > zz)) { // m[0][0] is the largest diagonal term
			if (xx< epsilon) {
				x = 0;
				y = 0.7071;
				z = 0.7071;
			} else {
				x = Math.sqrt(xx);
				y = xy/x;
				z = xz/x;
			}
		} else if (yy > zz) { // m[1][1] is the largest diagonal term
			if (yy< epsilon) {
				x = 0.7071;
				y = 0;
				z = 0.7071;
			} else {
				y = Math.sqrt(yy);
				x = xy/y;
				z = yz/y;
			}	
		} else { // m[2][2] is the largest diagonal term so base result on this
			if (zz< epsilon) {
				x = 0.7071;
				y = 0.7071;
				z = 0;
			} else {
				z = Math.sqrt(zz);
				x = xz/z;
				y = yz/z;
			}
		}
		return new axisAngle(angle,x,y,z); // return 180 deg rotation
	}
	// as we have reached here there are no singularities so we can handle normally
	double s = Math.sqrt((m[2][1] - m[1][2])*(m[2][1] - m[1][2])
		+(m[0][2] - m[2][0])*(m[0][2] - m[2][0])
		+(m[1][0] - m[0][1])*(m[1][0] - m[0][1])); // used to normalise
	if (Math.abs(s) < 0.001) s=1; 
		// prevent divide by zero, should not happen if matrix is orthogonal and should be
		// caught by singularity test above, but I've left it in just in case
	angle = Math.acos(( m[0][0] + m[1][1] + m[2][2] - 1)/2);
	x = (m[2][1] - m[1][2])/s;
	y = (m[0][2] - m[2][0])/s;
	z = (m[1][0] - m[0][1])/s;
   return new axisAngle(angle,x,y,z);
}

Thanks for the following corrections to this code:

Singularities

So how can we find this case? When we look at the formula for AxisAngle to Matrix the asymmetrical component about the leading diagonal is due to the sin(angle) component.

So when sin(0)=sin(180) = 0
Then the matrix [m]t = [m]

where:

Note: the transpose of a normalised matrix represents the inverse transform, so this is saying that rotation by 180° is the same as rotation by -180° and rotation by 0° is the same as rotation by -0°.

So if we test for symmetry:
Abs(m01-m10)<0.001 and Abs(m02-m20)<0.001 and Abs(m12-m21)<0.001
then
Angle = 0° or 180°

The angle = 0 case should be easy to find because the matrix will be:
1 0 0
0 1 0
0 0 1

So we can test for this case with:

Abs(m01+m10)<0.001 and Abs(m02+m20)<0.001 and Abs(m12+m21)<0.001

At 180° we need to calculate the axis, to do this we go back to the formula for AxisAngle to Matrix
c = cos(180) = -1
s = sin (180) = 0
t = 1-c = 2

so the matrix is:

2*x*x-1 2*x*y 2*x*z
2*x*y 2*y*y -1 2*y*z
2*x*z 2*y*z 2*z*z*-1

so we can calculate the axis terms x,y and z using the leading diagonal terms:
x = ±√((m00+1)/2)
y = ±√((m11+1)/2)
z = ±√((m22+1)/2)

There are a couple of problems with this:

To get round these problems we only use one of the diagonal tems and calculate the other terms relative to that. This means that it does not matter if we choose positive or negative for the first term we calculate. To understand this we must remember that we are are working with the specific case of rotation by 180° (π radians) which is the same as -180° (-π), reversing the angle is the same as reversing the axis, so if x,y,z is the axis we want then -x,-y,-z is also valid. So it is only the relative signs of these terms that matter and we can get these from the non diagonal terms of the matrix.

So, for instance, we can calculate the axis from the first column of the matrix:

x = Math.sqrt((m[0][0]+1)/2); // since m00= 2*x*x-1
y = m[0][1]/(2*x); // since m01= 2*x*y
z = m[0][2]/(2*x); // since m02= 2*x*z

This will only work if m[0][0]+1 >0, that is m[0][0] > -1, therefore make sure that we are taking the root of a positive number and to minimise any rounding errors we compare the diagonal terms and use only the column containing the highest diagonal term.

Also to minimise rounding errors we use the average of opposite terms, because: m[0][1]=m[1][0]=m[0][1]+m[1][0])/2
Using the average might help if the matrix has become de-orthogonalised due to rounding errors.

x = Math.sqrt((m[0][0]+1)/2);
y = (m[0][1]+m[1][0])/(4*x);
z = (m[0][2]+m[2][0])/(4*x);

So does this avoid any possible problems? (square root of negative or divide by zero). Because terms in a normalised matrix are between +1 and -1 then I think the worst case would be:

-1 0 0
0 -1 0
0 0 -1

This could cause divide by zero and any rounding errors could cause root of negative. So we need to test for this condition.

Issues

Examples

In order to check if the axis is reversed when we have a minus angle here are 2 examples, example 1 represents + 30° about the z axis and example 2 represents - 30 degrees about the z axis, so these examples should produce either, equal and opposite angles, or equal and opposite axes.

Example 1

In this example heading is + 30 degrees about the z axis.

matrix =
cos(heading) = 0.866 sin(heading) = 0.5 0
-sin(heading) = -0.5 cos(heading) = 0.866 0
0 0 1

angle = acos ( ( m00 + m11 + m22 - 1)/2)

angle = acos ( ( 0.866 + 0.866 + 1 - 1)/2)

angle = acos ( 0.866 )

angle = 30 degrees

x = (m21 - m12) = 0 - 0 =0
y = (m02 - m20) = 0 - 0 =0
z = (m10 - m01) = -0.5 - 0.5 = -1

example 2

now the same example with heading = -30 degrees about the z axis.

matrix =
cos(heading) = 0.866 sin(heading) = -0.5 0
-sin(heading) = 0.5 cos(heading) = 0.866 0
0 0 1

angle = 30 degrees

x = (m21 - m12) = 0 - 0 =0
y = (m02 - m20) = 0 - 0 =0
z = (m10 - m01) = 0.5 + 0.5 = 1

So reversing the angle reverses the axis instead of the angle. However these are equivalent (we can reverse both the angle and axis and still represent the same rotation). So this is equivalent to -30° about z = -1 axis.

This gives me some confidence that whatever quadrant the rotation is in the above formula will convert this correctly to the correct axis angle.

example 3

matrix =
0.96608673169969 -0.25800404198456 -0.01050433974302
0.25673182392846 0.95537412871306 0.14611312318926
-0.02766220194012 -0.14385474794174 0.98921211783846

angle = acos ( ( 0.966 + 0.955 + 0.989 - 1)/2)

angle = acos (0.955336489125605)

angle = 17.3 degrees

x = (m21 - m12) = -0.14385474794174 - 0.14611312318926 = -0.289967871131
y = (m02 - m20) = -0.01050433974302 + 0.02766220194012 = 0.0171578621971
z = (m10 - m01) = 0.25673182392846 + 0.25800404198456 = 0.51473586591302

This axis can then be normalised if required.

Here I am not sure that the matrix represents a pure rotation? The matrix is not skew symmetric?

Specially m02 and m20 as they are both small and have the same sign. This leads me to the question about how sensitive the conversion is to input errors? Perhaps a small error here might put the axis into the wrong quadrant. And this error would be increased when the axis is normalised. If the original matrix may not be orthogonal then, perhaps we should use a different conversion method, depending on the ratio of values in the leading diagonal as we did here: matrix to quaternion page.

Another Example

we take the 90 degree rotation from this: to this:

As shown here the matrix for this rotation is:

[R] =
1 0 0
0 0 -1
0 1 0

So using the above result:

angle = acos ( ( m00 + m11 + m22 - 1)/2) = acos (0) = 90 degrees
x = (m21 - m12)/√((m21 - m12)2+(m02 - m20)2+(m10 - m01)2) = 2/√(4) = 1
y = (m02 - m20)/√((m21 - m12)2+(m02 - m20)2+(m10 - m01)2) = 0
z = (m10 - m01)/√((m21 - m12)2+(m02 - m20)2+(m10 - m01)2) = 0

angle = 90°
axis = 1,0,0

As you can see here, this gives the result that we are looking for.

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° steps here

Here is a simple javascript calculator, input matrix values then press calculate button:

axis x  
axis y  
axis z  
angle radians
type

metadata block
see also:

 

Correspondence about this page

Book Shop - Further reading.

Where I can, I have put links to Amazon for books that are relevant to the subject, click on the appropriate country flag to get more details of the book or to buy it from them.

cover Introductory Techniques for 3-D Computer Vision by Emanuele Trucco, Alessandro Verri

 

This site may have errors. Don't use for critical systems.

Copyright (c) 1998-2014 Martin John Baker - All rights reserved - privacy policy.