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:
 Invert AxisAngle to Matrix equations
 Convert AxisAngle to quaternion then convert to Matrix
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: // https://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:
 [m]^{t} = transpose of matrix (exchange rows with columns)
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(m01m10)<0.001 and Abs(m02m20)<0.001 and Abs(m12m21)<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 = 1c = 2
so the matrix is:
2*x*x1  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:
 The result of the square root may be positive or negative so we need more information to know which is the case.
 If the values are very close to zero then rounding errors may flip the sign.
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*x1 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 deorthogonalised 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
 For issues concerning acos see here.
 Most trig functions (except openGL) use radians not degrees.
ExamplesIn 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 1In this example heading is + 30 degrees about the z axis.
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 example 2now the same example with heading = 30 degrees about the z axis.
angle = 30 degrees x = (m21  m12) = 0  0 =0 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
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 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] = 

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: