Home > Mathematica > Around the Twist

Around the Twist

I have become addicted to Mathematica Stackexchange. The unanswered questions are opportunities to think and improve one’s Mathematical and Mathematica skills as well as, hopefully, provide a useful answer to the questioner.

The community is active, vigilant. Comments, though sometimes terse, appear in general well motivated and  targeted and helpful.

On a recent question, I had the extremely embarrassing but ultimately rewarding experience of making conceptual and typographical errors. The question was could the axis and angle of rotation matrix (3D) be  recovered from the matrix.  My obvious starting point was that the eigenvector  associated with the eigenvalue 1 was parallel to the axis. The other eigenvalues provide the angle…then I went somewhat around the twist trying to get the axis orientation and angle to match.

A high reputation and experienced user,Jens, came up with a neat approach that worked for almost every rotation matrix (in practical terms it was universal: failing only at $-\pi, 0, \pi$ and for axis $(a,a,a), a\neq 0$. (I was moderately annoyed when a moderator rebuked my answer for its failure for $\pi$ and arbitrary axis which was also the case for approach of the experienced user).  For practical purposes, esp testing random angles $(-\pi,\pi)$ and random axes the the algorithm was universal (only failing if deliberately inputting an exception case).

My algorithm was still fatally flawed.  After, an embarrassing attempt to resolve this and completely erroneous critique of the excellent well constructed answer referred to above, I returned to the “real world”.

I decided I would use the principles of the solution from Jens, particularly once deriving the axis generating and orthonormal system (axis and two vectors orthogonal to the axis and each other), rotating these vectors and using their projections onto one of them to determine the angle of rotation. I persisted with the eigensystem approach as it dealt with finding the axis of the $\pi$ case as well as the other cases whereas the null space approach  of antisymmetric decomposition  ($R-R^T$) where $R$ is rotation matrix) yields  a zero matrix for $\pm \pi$., I dealt with axis $(a,a,a)$ by mapping this to $(a,-a,a)$ and just rotated the other axis vectors  which deals with the other cases (just another permutation).  At its heart this is Jens approach. The exercise for me was to go away (not look at the code though remembering the elements) and write my version of the idea.

So, what I finally came up with (again this is a learning exercise based on the idea of Jens):

 f[u_] := If[u == IdentityMatrix[3], {{0, 0, 0}, 0}, Module[ {ev, evec, axis, axisn, vec, veco, nvec, angle}, {ev, evec} = Eigensystem[u]; axis = First@Extract[evec, Position[ev, _?Positive]]; vec = If[Equal[axis], {1, -1, 1} axis, RotateLeft[axis]]; {axisn, veco} = Orthogonalize[{axis, vec}]; nvec = Cross[veco, axisn]; angle = N@Arg[Complex @@ (Dot[u.#, veco] & /@ {veco, nvec})]; {angle, axisn}]]; 
Also borrowing the testing approach:
 test[] := Module[ {angle = RandomReal[{-Pi, Pi}], axis = RandomReal[{0, 1}, {3}], rot, a1, a2}, rot = RotationMatrix[angle, axis]; {a1, a2} = f[rot]; Chop[a1 Dot[Normalize@axis, a2] - angle] === 0]

And @@ Table[test[], {1000}]

yielded True on multiple passes…

My errors were instructive and apart from the embarrassment, it was finally satisfying to understand the problem, my misconceptions and work through it myself.

Post script: as it turns out the Presentations package has a function that does this: RotationAngleAndAxes