[Chimera-users] Get angles in degree out of Rotation matrix

[Thomas Goddard]

Hi Michael,

  It is not very common to specify a 3-dimensional rotation in terms
of a x-axis rotation, y-axis rotation and z-axis rotation because this
decomposition depends on the order of the x, y, and z rotations.
If you rotate 90 degrees about the x-axis, then 90 about the y-axis
you get a different orientation than if you do the y-axis rotation
first followed by the x-axis rotation.

  Chimera has a routine to produce the rotation axis (a unit vector in
3-dimensions) and rotation angle about that vector corresponding to a
3x3 rotation matrix.  You could use that in Chimera with some Python code.
If you just want the rotation angle and don't need the axis then the
formula is simple

	Tr(R) = 1 + 2*cos(theta)

where
	Tr(R) = sum of rotation matrix diagonal elements

so the rotation angle

	theta = arccos((Tr(R) - 1)/2)

  Below I give the Chimera C++ code that also computes the axis.  There
are corner cases handled by this code for 0 degree and 180 degree rotations.

	Tom


void Xform::getRotation(Vector *axis, Real *angle) const
{
	Real cosTheta = (rot[0][0] + rot[1][1] + rot[2][2] - 1) / 2;
	if (fuzzyEqual(cosTheta, 1)) {
		*angle = 0;
		(*axis)[0] = (*axis)[1] = 0;
		(*axis)[2] = 1;
		return;
	}
	if (!fuzzyEqual(cosTheta, -1)) {
		*angle = degrees(acos(cosTheta));
		Real sinTheta = sqrt(1 - cosTheta * cosTheta);
		(*axis)[0] = (rot[2][1] - rot[1][2]) / 2 / sinTheta;
		(*axis)[1] = (rot[0][2] - rot[2][0]) / 2 / sinTheta;
		(*axis)[2] = (rot[1][0] - rot[0][1]) / 2 / sinTheta;
		return;
	}
	*angle = 180;
	if (rot[0][0] >= rot[1][1]) {
		if (rot[0][0] >= rot[2][2]) {
			// rot00 is maximal diagonal term
			(*axis)[0] = sqrt(rot[0][0] - rot[1][1] - rot[2][2]
					  			+ 1.0f) / 2;
			Real halfInverse = 1 / (2 * (*axis)[0]);
			(*axis)[1] = halfInverse * rot[0][1];
			(*axis)[2] = halfInverse * rot[0][2];
		} else {
			// rot22 is maximal diagonal term
			(*axis)[2] = sqrt(rot[2][2] - rot[0][0] - rot[1][1]
					  			+ 1.0f) / 2;
			Real halfInverse = 1 / (2 * (*axis)[2]);
			(*axis)[0] = halfInverse * rot[0][2];
			(*axis)[1] = halfInverse * rot[1][2];
		}
	} else {
		if (rot[1][1] >= rot[2][2]) {
			// rot11 is maximal diagonal term
			(*axis)[1] = sqrt(rot[1][1] - rot[0][0] - rot[2][2]
					  			+ 1.0f) / 2;
			Real halfInverse = 1 / (2 * (*axis)[1]);
			(*axis)[0] = halfInverse * rot[0][1];
			(*axis)[2] = halfInverse * rot[1][2];
		} else {
			// rot22 is maximal diagonal term
			(*axis)[2] = sqrt(rot[2][2] - rot[0][0] - rot[1][1]
					  			+ 1.0f) / 2;
			Real halfInverse = 1 / (2 * (*axis)[2]);
			(*axis)[0] = halfInverse * rot[0][2];
			(*axis)[1] = halfInverse * rot[1][2];
		}
	}
}