20
$\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$
0

7 Answers 7

15
$\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$
9
  • 1
    $\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, 2016 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$ Aug 8, 2016 at 9:57
  • $\begingroup$ That was dumb on many levels...bad me. $\endgroup$ Aug 8, 2016 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, 2018 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$ Oct 23, 2018 at 21:47
11
$\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$
3
  • 1
    $\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, 2013 at 21:53
  • 1
    $\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, 2014 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, 2014 at 23:40
4
$\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$
4
  • 1
    $\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, 2015 at 12:12
  • 1
    $\begingroup$ Excellent reference. :) $\endgroup$ Jul 18, 2016 at 3:26
  • 1
    $\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$ Dec 25, 2018 at 11:24
  • 1
    $\begingroup$ I guess I should have pointed this out myself, eh? Thanks. :) $\endgroup$ Dec 25, 2018 at 16:25
2
$\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$
1
  • 1
    $\begingroup$ Note that for $k$ a unit vector, $A^2 = kk^T - I$. $\endgroup$
    – user7530
    Feb 2, 2013 at 21:40
1
$\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$
2
  • $\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, 2014 at 14:14
  • 1
    $\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$ Jul 19, 2016 at 17:21
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$
1
  • $\begingroup$ What does the link have to do with this? $\endgroup$
    – Vladhagen
    Nov 6, 2014 at 5:01
0
$\begingroup$

For a solution that works in any dimensions and leaves vectors perpendicular to $a$ and $b$ unchanged, see Does this orthogonal matrix parametrization have a name?

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .