15
$\begingroup$

I have written an algorithm for solving the following problem: Given two 3d-vectors, say: $a,b$, find rotation of $a$ so that its orientation matches $b$.

However, I am not sure if the following algorithm works in all cases:

1) Find axis and angle using cross product and dot product:

$$\mathbf{x}=\frac{a\times b}{||a\times b||}\\ \theta=\cos^{-1}(\frac{a\cdot b}{||a||\cdot ||b||})$$

3) Find rotation matrix using exponential map:

$$\mathbf{R}=e^{\mathbf{A}\theta} =\mathbf{I}+\sin(\theta)\cdot \mathbf{A}+\left(1-\cos(\theta)\right) \cdot \mathbf{A}^{2}$$

where $\mathbf{A}$ is a skew-symmetric matrix corresponding to $\mathbf{x}$:

$$\mathbf{A}=[\mathbf{x}]_{\times}=\begin{bmatrix}0 & -\mathbf{x}_{3} & \mathbf{x}_{2} \\ \mathbf{x}_{3} & 0 & -\mathbf{x}_{1} \\ -\mathbf{x}_{2} & \mathbf{x}_{1} & 0\end{bmatrix}$$

Notes:

The axis is computed using cross product as this gives vector perpendicular to both $a$ and $b$. Only direction of the axis is important, hence it is divided by its magnitude. However, I am not sure if $\mathbf{x}$ will always have the proper direction (the result can be $-\mathbf{x}$ instead of $\mathbf{x}$?).

The rotation matrix is computed using Rodrigues' rotation formula.

Finally, the vector $\mathbf{R}a$ should have same direction as $b$.

I have tested this numerically and it seems working, but I would like to be sure the formulas work for any two $a,b$.

$\endgroup$
7
$\begingroup$

This is the right general approach, but the corner case $\|a\times b\| \approx 0$ must be handled.

If $\theta < \epsilon,$ $R=I$.

If $\pi-\theta < \epsilon$, you can choose for $\mathbf{x}$ any vector orthogonal to $\mathbf{a}$, for instance $\mathbf{x} = \frac{\mathbf{a} \times e_i}{\|\mathbf{a}\times e_i\|}$, where $i$ is the index of the component of $\mathbf{a}$ with least magnitude.

$\endgroup$
  • $\begingroup$ I was aware of the first case (very small angle), but you have shown the angle approaching 180° is also important to handle. Thanks. $\endgroup$ – Libor Feb 2 '13 at 21:53
  • $\begingroup$ Libor, user7530: Can anyone explain the origin of rotation matrix formula? I've no idea about rotation matrix and skew matrix. I'd also like to know what happens at 0 and 180 to handle the cases separately? I also have no idea about orthogonal vector. I can understand till finding the angle between vectors and the axis of rotation. $\endgroup$ – cegprakash Jun 9 '14 at 23:17
  • $\begingroup$ @cegprakash Have you read the article on Wikipedia about the Rodrigues Rotation Formula? What you ask is a rather broad topic, you may want to try asking specific, new questions. $\endgroup$ – user7530 Jun 9 '14 at 23:40
9
$\begingroup$

I know this is a long-dead and well-answered question, but I woke this morning thinking "All of the answers involve an "if" statement, i.e., all are discontinuous functions of the inputs. Is there a continuous answer, i.e., a continuous function $$ R: S^2 \times S^2 \to SO(3) $$ that takes a pair of vectors $u, v$ to a rotation matrix $R(u, v)$ with the property that $R(u, v) u = v$?"

The paper that Tomas Moller and I wrote back in 1999, for instance, uses "the coordinate vector corresponding to the smallest entry of $w$" for some vector $w$, which doesn't vary continuously as a function of $w$. And I wondered, "Did we really do as well as possible, or might there have been a continuous solution?"

The answer is "no." But the proof uses a bit of topology.

Fix the vector $u$ (set it to be $e_1$, for instance), and look at the map $$ K : S^2 \to SO(3) : v \mapsto R(u, v). $$ Then compose this with the map $$ H : SO(3) \to S^2 : M \mapsto Mu. $$ The composite map $$ H\circ K: S^2 \to S^2 $$ is $v \mapsto R(u, v)u = v$, i.e., the identity map on $S^2$. But that means that $$ (H\circ K)_{*} : H_2(S^2) \to H_2(S^2), $$ the induced map on second homology, must be the identify from $\mathbb Z $ to $\mathbb Z$. But since $$ (H\circ K)_{*} = H_{*} \circ K_{*}$$ this map must factor through $H_2(SO(3)) = 0$, which is impossible.

Thus: There's no continuous solution to the "rotate one vector to another" problem, a fact that I should have mentioned back in our original paper. Sigh. Hindshight is 20-20.

$\endgroup$
  • $\begingroup$ Thanks for elaborate answer. I needed this for a software project, so it's "good enough" for me to just calculate the rotation. When it works in practice, I am happy. I believe your contribution is very helpful for someone more mathematically oriented. $\endgroup$ – Libor Jul 19 '16 at 14:30
  • $\begingroup$ Assuming can do a branch free select you can form a branch free version by selecting between $\frac{1-c}{1-c^2}$ and $\frac{1}{1+c}$ if I'm thinking correctly. $\endgroup$ – MB Reynolds Aug 8 '16 at 9:57
  • $\begingroup$ That was dumb on many levels...bad me. $\endgroup$ – MB Reynolds Aug 8 '16 at 15:28
  • $\begingroup$ I am glad this answer exists. Could someone post a reference to which topological concepts are being used? What is the second homology, and how does it relate between $S^2$ and $\mathbb{Z}$? $\endgroup$ – Gus Oct 23 '18 at 21:30
  • $\begingroup$ The key concepts are (1) The notion of homology groups, (2) The fact that continuous maps induce homomorphisms on homology groups in a nice way, (3) the computation of $H_2(S^2)$, which is easy once you know about homology, and (4) the computation of $H_2(SO(3))$, for which you need to know that $SO(3)$ has $S^3$ as its universal cover (the "quaternion representation of rotations"), and the computation of $H_2(S^3)$, which is even easier than #3. The first three are covered in pretty much any intro to algebraic topology. I like Marvin Greenberg's book, but there's surely something more modern. $\endgroup$ – John Hughes Oct 23 '18 at 21:47
2
$\begingroup$

Let $F_n(x) = x - {2 \over n . n} n (n \ . \ x)$ be the transformation that reflects $x$ through the plane that is perpendicular to $n$. Composing two reflections gives a rotation: if the angle from $a$ to $b$ is $\phi$ then $F_b(F_a(x))$ rotates $x$ on the plane spanned by $a$ and $b$ by $2 \phi$.

Given the normalized vectors $\hat a=a/|a|$, $\hat b=b/|b|$, $\hat c=c/|c|$, where $c = \hat a + \hat b$, the angle from $\hat a$ to $\hat c$ is $\phi$, half the angle from $\hat a$ to $\hat b$.

In fact, $F_{\hat c}(F_{\hat a}(x))$ rotates $x$ by $\phi$ on the plane spanned by $a$ and $b$. This is valid even if $a$ and $b$ are parallel. This avoids computing a cross product, inverse cosine, sine and cosine, or division by a magnitude that can be arbitrarily close to zero.

$\endgroup$
  • $\begingroup$ On second thought, this still isn't valid if the two vectors point in opposite directions. Other than expecting the user to explicitly provide a bisecting vector, reflect on (t-a) and then (t-b) where t is a coordinate vector corresponding to the smallest component of the vector a. This is detailed in this paper: cs.brown.edu/~jfh/papers/Moller-EBA-1999/paper.pdf $\endgroup$ – Centrinia Mar 18 '15 at 12:12
  • $\begingroup$ Excellent reference. :) $\endgroup$ – John Hughes Jul 18 '16 at 3:26
  • $\begingroup$ @JohnHughes Cute:) I was just googling this stuff myself, and, just to update Centrinia's now-dead link, cs.brown.edu/research/pubs/pdfs/1999/Moller-1999-EBA.pdf is where your paper's at now. $\endgroup$ – John Forkosh Dec 25 '18 at 11:24
  • $\begingroup$ I guess I should have pointed this out myself, eh? Thanks. :) $\endgroup$ – John Hughes Dec 25 '18 at 16:25
0
$\begingroup$

I'm not clear on why you have a factor of $A^2$ in your expression for $R$. In particular, wikipedia lists the matrix form for the Rodrigues formula as

$$R = I \cos \theta + A \sin \theta + (1-\cos \theta) x x^T$$

$\endgroup$
  • $\begingroup$ Note that for $k$ a unit vector, $A^2 = kk^T - I$. $\endgroup$ – user7530 Feb 2 '13 at 21:40
0
$\begingroup$
function [ R ] = RotAtoB( a,b )
    x = [a(2)*b(3) - b(2)*a(3);a(3)*b(1) - b(3)*a(1);a(1)*b(2) - b(1)*a(2)];
    x = x/norm(x);
    theta = acos(a'*b/(norm(a)*norm(b)));
    A = [0    -x(3)  x(2)
         x(3)   0   -x(1)
        -x(2)  x(1)   0  ];
    R = eye(3) + sin(theta)*A + (1-cos(theta))*A^2;
end

Calculating the powers of matrix : Nth power of a square matrix and the Binet Formula for Fibonacci sequence

$\endgroup$
  • $\begingroup$ What does the link have to do with this? $\endgroup$ – Vladhagen Nov 6 '14 at 5:01
0
$\begingroup$

I have a simpler method comes from Erigen's "Mechanics of Continua". Here R is rotational matrix that rotate vector "a" align with vector "b"

Matlab Code:

%%%%%% Rotate vector a align with vector b%%%%%%%%%%
syms ax ay az bx by bz k real

a=[ax ay az]'

au=a./sqrt(ax^2+ay^2+az^2)

b=[bx by bz]'

bu=b./sqrt(bx^2+by^2+bz^2)

R=[bu(1)*au(1) bu(1)*au(2) bu(1)*au(3);

    bu(2)*au(1) bu(2)*au(2) bu(2)*au(3)

    bu(3)*au(1) bu(3)*au(2) bu(3)*au(3)]

% You can verify it by type

c=R*a

cu=c./sqrt(c(1)^2+c(2)^2+c(3)^2)

simple(bu-cu)

%the result is zero means c(a after rotation) and b are aligned with each other.  

simple(sqrt(c(1)^2+c(2)^2+c(3)^2)-sqrt(c(1)^2+c(2)^2+c(3)^2))

%the result is zero means c(a after rotation) and a are of the same length


%%%%%%%End%%%%%%%%%%%%
$\endgroup$
  • $\begingroup$ Thanks. Could you please point me to the theory? The matrix seems to be simply $R=ab^{T}$ where $a$ and $b$ are normalized vectors. It is not clear to me why $R^{T}=R^{-1}$ (i.e. $R^{T}R=I$) in that case. The book is unfortunately hard to obtain and too much hassle just to get this information... $\endgroup$ – Libor Jul 26 '14 at 14:14
  • $\begingroup$ The matrix $R$ is just $b a^t$ (if $a$ and $b$ are unit vectors, hence it has rank 1, and therefore is not a rotation matrix. So this answer, although simple, is wrong. $\endgroup$ – John Hughes Jul 19 '16 at 17:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.