# Axis/Angle from rotation matrix

With r = RotationMatrix[a, {x, y, z}] I can compute a 3D rotation matrix from its axis/angle representation.

Given a 3D rotation matrix r, how can I compute a and {x, y, z}?

Example:

r = {{0.966496, -0.214612, 0.14081},
{0.241415, 0.946393, -0.214612},
{-0.0872034, 0.241415, 0.966496}}


The result should be a = 20. Degree and {x, y, z} = {2, 1, 2}/3 (or equivalent).

Edit: I am fine with any answer that gives the same r when applied to RotationMatrix.

• I'm not sure this question is well posed in the sense that the answer is non-unique. A simple example would be r=RotationMatrix[1,{1,0,0}]=RotationMatrix[-1,{-1,0,0}] are the same. So which one are you intending to find? Commented Aug 6, 2013 at 13:08
• The one where the angle is in [0,2 pi]. I think it is well posed as angle and axis are directly linked and you can provide a unique answer with a perfectly reasonable additional assumption. Commented Aug 6, 2013 at 13:43
• Actually there was an answer earlier which was just fine but it was deleted? Commented Aug 6, 2013 at 13:46
• @Danvil it was deleted because it was incorrect, so he deleted it. It was incorrect, as it worked only for the example you provided and could not be generalized. Commented Aug 6, 2013 at 14:42
• @Danvil I'm thankful I noticed the wink. It would have been a bad day, otherwise. :P Commented Aug 6, 2013 at 16:38

There is no need to use Eigensystem or Eigenvectors to find the axis of a rotation matrix. Instead, you can read the axis vector components off directly from the skew-symmetric matrix

$$a \equiv R^T-R$$

In three dimensions (which is assumed in the question), applying this matrix to a vector is equivalent to applying a cross product with a vector made up of the three independent components of $a$:

{1, -1, 1}Extract[a, {{3, 2}, {3, 1}, {2, 1}}]


This one-line method of finding the axis is applied in the following function. To get the angle of rotation, I construct two vectors ovec, nvec perpendicular to the axis and to each other, to find the cosine and sine of the angle using the Dot product (could equally have used Projection). To get a first vector ovec that is not parallel to the axis, I permute the components of the axis vector using the fact that

Solve[{x, -y, z} == {y, z, x}, {x, y, z}]
(* ==> {{x -> 0, y -> 0, z -> 0}} *)


which means the above permutation with sign change of a nonzero axis vector is always different from the axis. This is sufficient to use Orthogonalize and Cross to get the desired orthogonal vectors.

axisAngle[m_] := Module[
{axis, ovec, nvec
},
{axis, ovec} =
Orthogonalize[{{1, -1, 1} #, Permute[#, Cycles[{{1, 3, 2}}]]}] &@
Extract[m - Transpose[m], {{3, 2}, {3, 1}, {2, 1}}];
(* nvec is orthogonal to axis and ovec: *)
nvec = Cross[axis, ovec];
{axis,
Arg[Complex @@ (((m.ovec).# &) /@ {ovec, nvec})]}
]


The angle is calculated with Arg instead of ArcTan[x, y] here because the latter throws an error for x = y = 0.

Here I test the results of the function for 100 random rotation matrices:

testRotation[] := Module[
{m, a, axis, ovec, nvec,
v = Normalize[RandomReal[{0, 1}, {3}]],
α = RandomReal[{-Pi, Pi}], angle
},
m = RotationMatrix[α, v];
{axis, angle} = axisAngle[m];
Chop[
angle Dot[v, axis] - α
] === 0
]

And @@ Table[testRotation[], {100}]

(* ==> True *)


In the test, I have to account for the fact that if the function axisAngle defined the axis vector with the opposite sign as the random test vector, I have to reverse the sign of the rotation angle. This is what the factor Dot[v, axis] does.

Explanation of how the axis results from a skew-symmetric matrix

If $\vec{v}$ is the axis of rotation matrix $R$, then we have both $R\vec{v} = \vec{v}$ and $R^T\vec{v} = \vec{v}$ because $R^T$ is just the inverse rotation. Therefore, with $a \equiv R^T-R$ as above, we get

$$a \vec{v} = \vec{0}$$

Now the skew-symmetric property $a^T = -a$, which can be seen from its definition, means there are exactly three independent matrix element in $a$. They can be arranged in the form of a 3D vector $\vec{w}$ which must have the property $a \vec{w} = 0$. This vector is obtained in the Extract line above. In fact, $a \vec{x} = \vec{w}\times \vec{x}$ for all $\vec{x}$, and hence if $a \vec{x} = 0$ then $\vec{x}\parallel\vec{w}$. Therefore, the vector $\vec{v}$ is also parallel to $\vec{w}$, and the latter is a valid representation of the rotation axis.

Edit 2: speed considerations

Since the algorithm above involves only elementary operations that can be compiled, it makes sense that a practical application of this approach would use Compile. Then the function could be defined as follows (keeping the return values arranged as above):

Clear[axisAngle1, axisAngle]

axisAngle1 =
Compile[{{m, _Real, 2}},
Module[{axis, ovec, nvec, tvec, w, w1},
tvec = {m[[3, 2]] - m[[2, 3]], m[[3, 1]] - m[[1, 3]],
m[[2, 1]] - m[[1, 2]]};
If[tvec == {0., 0., 0.},
{#/Sqrt[#.#] &[#[[Last@Ordering[N[Abs[#]]]]] &[
1/2 (m + {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}})]],
If[Sum[m[[i, i]], {i, 3}] == 3, 0, Pi] {1, 1, 1}},
axis = {1, -1, 1} tvec;
axis = axis/Sqrt[axis.axis];
w = {tvec[[2]], tvec[[3]], tvec[[1]]};
ovec = w - axis Dot[w, axis];
nvec = Cross[axis, ovec];
w1 = m.ovec;
{axis, {1, 1, 1} ArcTan[w1.ovec, w1.nvec]}
]
]
];

axisAngle[m_] := {#1, Last[#2]} & @@ axisAngle1[m]


The results are the same as for the previous definition of axisAngle, but I now get a much faster execution as can be seen in this test:

tab = RotationMatrix @@ # & /@
Table[{RandomReal[{-Pi, Pi}],
Normalize[RandomReal[{0, 1}, {3}]]}, {100}];
timeAvg =
Function[func,
Do[If[# > 0.3, Return[#/5^i]] & @@ Timing@Do[func, {5^i}], {i, 0,
15}], HoldFirst];

timeAvg[axisAngle /@ tab]

(* ==> 0.000801259 *)


This is more than an order of magnitude faster than the un-compiled version. I removed Orthogonalize from the code because I didn't find it in the list of compilable functions.

Note that Eigensystem is not in that list, either.

Edit 3

The first version of axisAngle demonstrated the basic math, but the compiled version axisAngle1 (together with the re-defined axisAngle as a wrapper) is faster. One thing that was missing was the correct treatment of the edge case where the rotation is by exactly $\pi$ in angle.

I added that fix only to the compiled version (axisAngle1) because I think that's the more practical version anyway. The trivial case of zero rotation angle was already included in the earlier version.

To explain the added code, first note that for angle $\pi$ you can't read off the axis from $R^T - R$ because the resulting matrix vanishes. To get around this singular case, we can use the geometric fact that a rotation by $\pi$ is equivalent to an inversion in the plane perpendicular to the rotation axis given by the unit vector $\vec{n}$. Therefore, if we form the sum of a vector $\vec{v}$ and its $\pi$-rotated counterpart, the components transverse to the rotation axis cancel and the result is always parallel to the axis. In matrix form,

$$(R+1)\vec{v} = 2\vec{n}(\vec{n}\cdot\vec{v}) = 2\left(\vec{n}\vec{n}^T\right)\vec{v}$$

Since this holds for all vectors, it is a matrix identity. The right-hand side contains a matrix $\vec{n}\vec{n}^T$ which must have at least one row that's nonzero. This row is proportional to $\vec{n}^T$, so you can read of the axis vector directly from $(R+1)$, again without any eigenvalue computations.

• Beautiful! (+1) My hamfisted approach to this has been an opportunity to learn. Commented Aug 6, 2013 at 22:28
• Note (as you have outlined) you could also obtain the axis from NullSpace[m-Transpose[m]] Commented Aug 7, 2013 at 4:17
• @ubpdqn True, but that takes more computation... as I said, you can read the axis off literally from the matrix a itself.
– Jens
Commented Aug 7, 2013 at 5:15
• @J.M. Regarding the paper, it seems that the referee observed what I would also say: that the basic idea behind the paper is well-known (and this also applies to my proof).
– Jens
Commented Jul 6, 2015 at 2:59
• @David That matrix corresponds to the edge case I mention in edit 3: It's a rotation by $\pi$. The result of the compiled function in edit 3 give the correct result. Note that such a matrix is degenerate and therefore the eigenvector with eigenvalue $-1$ is not uniquely determined.
– Jens
Commented Oct 2, 2022 at 16:44

Here is a slightly shorter re-implementation of Jens's method, which uses the "Pixar method" (itself a modification of Frisvad's method) to generate an orthogonal basis:

axisAngleC = Compile[{{m, _Real, 2}},
Module[{ang, axis, ovec, nn, nvec, rm, s, w, w1, wm, xx, yy, zz},
axis = {m[[3, 2]] - m[[2, 3]], m[[1, 3]] - m[[3, 1]], m[[2, 1]] - m[[1, 2]]};
nn = Norm[axis];
If[nn == 0., (* singular case *)
ang = If[Total[Transpose[m, {1, 1}]] == 3., 0., π];
rm = 0.5 (m + {{1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}});
axis = rm[[Ordering[Max /@ Abs[rm], -1][[1]]]];
axis /= Norm[axis],
(* non-singular case *)
axis /= nn; {xx, yy, zz} = axis;
s = 2 UnitStep[zz] - 1; w = -1/(s + zz); w1 = xx yy w;
ovec = {1 + s w xx xx, s w1, -s xx};
nvec = {w1, s + w yy yy, -yy};
wm = m.ovec; ang = Arg[wm.ovec + I wm.nvec]];
Append[axis, ang]]];


Check with some randomly generated axis-angle pairs:

(* https://mathematica.stackexchange.com/a/111693 *)
rodrigues[th_, axis_?VectorQ] :=
First[LinearAlgebraMatrixPolynomial[{{1, Sin[th], 2 Sin[th/2]^2}},
-LeviCivitaTensor[3, List].Normalize[axis]]]

Max[Table[With[{v = Normalize[RandomVariate[NormalDistribution[], 3]],
th = RandomReal[π]},
Norm[axisAngleC[rodrigues[th, v]] - Append[v, th], ∞]], {10^3}]] // Chop
0


Check singular cases:

axisAngleC[IdentityMatrix[3]]
{0., 0., 1., 0.}

Max[Table[With[{v = Normalize[RandomVariate[NormalDistribution[], 3]]},
Norm[Cross[Most[axisAngleC[RotationMatrix[π, v]]], v], ∞]], {10^3}]] // Chop
0


Here is a version for exact input that returns an inactivated RotationMatrix[]:

axisAngle[m_?OrthogonalMatrixQ] /; Det[m] == 1 :=
Module[{ang, axis, ovec, nn, nvec, rm, s, w, w1, wm, xx, yy, zz},
axis = {m[[3, 2]] - m[[2, 3]], m[[1, 3]] - m[[3, 1]], m[[2, 1]] - m[[1, 2]]};
nn = Simplify[Norm[axis]];
If[nn == 0, (* singular case *)
ang = π Boole[Total[Diagonal[m]] < 3];
rm = (m + IdentityMatrix[3])/2;
axis = Normalize[Extract[rm, Ordering[Max /@ Abs[rm], -1]]],
(* non-singular case *)
{xx, yy, zz} = Simplify[axis/nn];
s = 2 UnitStep[zz] - 1; w = -1/(s + zz); w1 = xx yy w;
ovec = {1 + s w xx xx, s w1, -s xx};
nvec = {w1, s + w yy yy, -yy};
wm = m.ovec; ang = Arg[Simplify[wm.ovec + I wm.nvec]]];
Inactive[RotationMatrix][ang, axis]]


Since this question still seems to be alive after some time, I'm giving a solution from Presentations, which I sell. It has a routine RotationAngleAndAxis and here it is used on the example.

<< Presentations

r = {{0.966496, -0.214612, 0.14081}, {0.241415,
0.946393, -0.214612}, {-0.0872034, 0.241415, 0.966496}};
RotationAngleAndAxis[r]

(* {0.349066, {0.666667, 0.333333, 0.666667}} *)


Here is a short example of its use that checks the result. The axis vector is normalized. Notice that reversing the angle of rotation and the axis gives an equivalent rotation.

rotations =
Table[{RandomReal[{-π, π}],
Normalize[Array[RandomReal[{-1, 1}] &, {3}]]}, {2}]

(* {{2.90598, {-0.596373, -0.74938,
0.287697}}, {-2.44158, {0.331763, -0.943343, 0.00605196}}} *)

rotationMatrices = RotationMatrix @@@ rotations;
anglesAndAxes = RotationAngleAndAxis /@ rotationMatrices

(* {{2.90598, {-0.596373, -0.74938, 0.287697}}, {2.44158, {-0.331763,
0.943343, -0.00605196}}} *)

resultingMatrices = RotationMatrix @@@ anglesAndAxes;

rotationMatrices - resultingMatrices // Chop
Total[Flatten[%]]

(* {{{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}, {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}} *)
(* 0 *)


The following is a test on 1000 initial rotation specifications. Notice also that I am comparing the resulting rotation matrices by a slightly risky procedure.

rotations =
Table[{RandomReal[{-π, π}],
Normalize[Array[RandomReal[{-1, 1}] &, {3}]]}, {1000}];
rotationMatrices = RotationMatrix @@@ rotations;
Timing[anglesAndAxes = RotationAngleAndAxis /@ rotationMatrices;]
resultingMatrices = RotationMatrix @@@ anglesAndAxes;
Total[Flatten[rotationMatrices - resultingMatrices // Chop]]

(* {0.873606, Null} *)
(* 0 *)


Presentations also has an EulerAngles routine that will return the two sets of Euler angles for any of the possible sequence of rotations specified by strings. For example, here are the two sets of Euler angles of the example for two different rotation sequences.

EulerAngles[r, "ZXZ"]
EulerAngles[r, "XYZ"]

(* {{-0.346633, 0.259587, 0.580661}, {2.79496, -0.259587, -2.56093}} *)
(* {{0.244775, 0.0873143, 0.244775}, {-2.89682, 3.05428, -2.89682}} *)


The axis of rotation $n$ can be obtained as the eigenvector corresponding to the eigenvalue $1$.

r = {{0.966496, -0.214612, 0.14081}, {0.241415,
0.946393, -0.214612}, {-0.0872034, 0.241415, 0.966496}};

es = Eigensystem[r];
pos = Nearest[Thread[es[[1]] -> Range[3]], 1][[1]];
n = es[[2, pos]]
(*  {0.666667, 0.333333, 0.666666} *)


The angle of rotation can be obtained directly from the trace of $r$.

ArcCos[(Tr[r] - 1)/2] 180/π
(* 20. *)


(This was my earlier and longer way to obtain the angle of rotation. As rcollyer pointed out this will have issues if $\pi k,k\in \mathcal{Z}$)

The rotation matrix in terms of the angle of rotation $\theta$ and axis of rotation $\{n1, n2, n3\}$.

R[t_, {n1_, n2_, n3_}] :=
With[{n = {{0, -n3, n2}, {n3, 0, -n1}, {-n2, n1, 0}}, id = IdentityMatrix[3]},
Inverse[id - t n].(id + t n)]


Here $t =\tan \left(\frac{\theta }{2}\right)$.

We get three solutions from the element in the first position.

eqns = Flatten[R[t, es[[2, pos]]] - r] // Simplify;
sols = NSolve[eqns[[1]] == 0, t]
(* {{t -> 3.88248*10^16}, {t -> 0.176327}, {t -> -0.176327}} *)


Choose the one which satisfies all the entries.

Chop[eqns /. sols, 10^-5];
Position[%, {0 ..}];
solst = Extract[sols, %]
(*  {{t -> 0.176327}}  *)


Obtain it in terms of degrees.

2 ArcTan[t] 180/π /. solst
(* {20.} *)

• It doesn't seem to work with a rotation of $\pi$ about the {0,0,1}. Commented Aug 6, 2013 at 14:39
• @rcollyer, that is to be expected because the rotation group $SO(3)$ cannot be covered by a single chart. Each chart will have its cases that are problematic. I used this parametrization and it worked for the OP's question. Still using this and going with sols = NSolve[eqns[[2]] == 0, t] we get $t=0$ for your case. The technique I used for the axis of rotation should be fine for all cases, however. Commented Aug 6, 2013 at 14:55
• True, but I think it is amenable to decomposition into the generators of rotation, i.e. the rotations should have eigenvalues of the form $\pmatrix{e^{i \theta} & 1 & e^{-i \theta}}$, or for real rotations $\pmatrix{\cos\theta + i \sin\theta & 1 & \cos\theta - i \sin\theta}$. Thus, the angle is easily recoverable. As to the direction, I have not yet looked through your code. Commented Aug 6, 2013 at 15:05
• Thanks for pointing that out. With your approach the angle of rotation is ambiguous for the OP's question, i.e we don't know if it is +20 or -20, whereas with the parametrization I used only + 20 satisfied all conditions. Probably quaternions will give a unique correct solution, but I do not have the material before me now. Commented Aug 6, 2013 at 15:20
• Well you can choose the angle to be between $0$ and $2 \pi$, to fix the branch. Commented Aug 6, 2013 at 15:25

My +1 to Jens.

For known Euler angles:

   {Φ, Θ, Ψ} = {0.1, 0.2, 0.3};
m = RotationMatrix[Φ, {0, 0, 1}].RotationMatrix[Θ, {1, 0,
0}].RotationMatrix[Ψ, {0, 0, 1}];
axisAngle[m_] :=Module[{axis, ovec, nvec}, {axis, ovec} =
Orthogonalize[{{1, -1, 1} #, Permute[#, Cycles[{{1, 3, 2}}]]}] &@
Extract[m - Transpose[m], {{3, 2}, {3, 1}, {2, 1}}];
nvec = Cross[axis, ovec];
{Arg[Complex @@ (((m.ovec).# &) /@ {ovec, nvec})], axis}];
N[axisAngle[m]]
(*  {0.446615, {0.448552, -0.0450053, 0.892623}}  *)


I am afraid that all of the above answers lack Mathematica's one line simplicity.

One can simply use HodgeDual@*Chop@*MatrixLog to get (-) the angle multiplied by the unit vector, where MatrixLog inverts the matrix exponent, and HodgeDual collects the skew-symmetric elements (the Choping is so that the HodgeDual will recognize the dimension correctly, alternatively, one can use LeviCivitaTensor[3]).

In a one line function:

R |-> With[{v = HodgeDual@Chop@MatrixLog@R}, Inactive[RotationMatrix][Norm@v, -Normal@Normalize@v]]