Skip to main content

linear algebra - 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.




Answer



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.


Comments

Popular posts from this blog

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...

What is and isn't a valid variable specification for Manipulate?

I have an expression whose terms have arguments (representing subscripts), like this: myExpr = A[0] + V[1,T] I would like to put it inside a Manipulate to see its value as I move around the parameters. (The goal is eventually to plot it wrt one of the variables inside.) However, Mathematica complains when I set V[1,T] as a manipulated variable: Manipulate[Evaluate[myExpr], {A[0], 0, 1}, {V[1, T], 0, 1}] (*Manipulate::vsform: Manipulate argument {V[1,T],0,1} does not have the correct form for a variable specification. >> *) As a workaround, if I get rid of the symbol T inside the argument, it works fine: Manipulate[ Evaluate[myExpr /. T -> 15], {A[0], 0, 1}, {V[1, 15], 0, 1}] Why this behavior? Can anyone point me to the documentation that says what counts as a valid variable? And is there a way to get Manpiulate to accept an expression with a symbolic argument as a variable? Investigations I've done so far: I tried using variableQ from this answer , but it says V[1...