Extract a robust rotation matrix from an arbitray 3x3 matrix

extractRotation.c++

void extractRotation(const Matrix3d &A, Quaterniond &q,
	const unsigned int maxIter)
{
	for (unsigned int iter = 0; iter < maxIter; iter++)
	{
		Matrix3d R = q.matrix();
		Vector3d omega = (R.col(0).cross(A.col(0)) + R.col
			(1).cross(A.col(1)) + R.col(2).cross(A.col(2))
			) * (1.0 / fabs(R.col(0).dot(A.col(0)) + R.col
				(1).dot(A.col(1)) + R.col(2).dot(A.col(2))) +
			1.0e-9);
			double w = omega.norm();
			if (w < 1.0e-9)
				break;
			q = Quaterniond(AngleAxisd(w, (1.0 / w) * omega)) *
			q;
			q.normalize();
	}
}

This is the paper's code of A Robust Method to Extract the Rotational Part of Deformations by Muller etc.

I think that the epsilon value should be inside 1/(x+epsilon):
(1.0 / (fabs(R.col(0).dot(A.col(0)) +
R.col(1).dot(A.col(1)) +
R.col(2).dot(A.col(2))) +
1.0e-9));

Certainly seems like it, in order to avoid div0. The original paper has it outside the divide as well in the code, but it is right in the text.