167

I have one triangle in 3d space that I am tracking in a simulation. Between time steps I have the the previous normal of the triangle and the current normal of the triangle along with both the current and previous 3d vertex positions of the triangles.

Using the normals of the triangular plane I would like to determine a rotation matrix that would align the normals of the triangles thereby setting the two triangles parallel to each other. I would then like to use a translation matrix to map the previous onto the current, however this is not my main concern right now.

I have found this website http://forums.cgsociety.org/archive/index.php/t-741227.html that says I must

  • determine the cross product of these two vectors (to determine a rotation axis)
  • determine the dot product ( to find rotation angle)
  • build quaternion (not sure what this means)
  • the transformation matrix is the quaternion as a 3 by 3 ( not sure)

Any help on how I can solve this problem would be appreciated.

user1084113
  • 1,849
  • 4
  • 13
  • 10
  • 2
    This may be helpful: http://gamedev.stackexchange.com/questions/20097/how-to-calculate-a-3x3-rotation-matrix-from-2-direction-vectors – mboratko Aug 08 '12 at 21:21
  • 13
    I know it's not your main concern right now, but I suspect it will become a concern later: There's no reason to expect that after applying an arbitrary rotation aligning the normals the triangles will be related by a translation -- you'd still have to rotate around the normal to align them. Stated differently, it's not clear that aligning the normals is a good first step, since you then have to perform two separate rotations. You can get the axis of the full rotation by taking the cross-product of the changes in the differences between two pairs of vertex positions. – joriki Aug 09 '12 at 05:13
  • 2
    That is a very good point I had not even though of that – user1084113 Aug 09 '12 at 14:23
  • @joriki I am not sure what you mean by "changes in the differences" do you mean to cross the vectors (A' - A) with (B'- B) where A, B are vertex positions of the previous triangle and the primes are the positions of the corresponding vertices on the current triangle, (I am having a hard time picturing as to why this cross product would give an axis that would allow the entire rotation) – user1084113 Aug 10 '12 at 13:43
  • 3
    @user1084113: No, that would be the cross-product of the changes in two vertex positions; I was talking about the cross-product of the changes in the differences between two pairs of vertex positions, which would be $((A-B)-(A'-B'))\times((B-C)\times(B'-C'))$. This gives you the axis of rotation (except if it lies in the plane of the triangle) because the translation drops out due to the differences, so this is purely the change in two different vectors due to rotation. The translation wouldn't drop out in your version. – joriki Aug 10 '12 at 13:51
  • @joriki Sorry if this may sound silly. Given that the resulting vector from the cross product is the axis of rotation, through which point does this axis of rotation pass through? Is it through the origin or one of the points of the triangle, I am trying to physically imagine it. – user1084113 Aug 10 '12 at 14:29
  • 4
    @user1084113: There's no canonical answer to that question. An orientation-preserving isometry of $\mathbb R^3$ can be written as a composition of a rotation and a translation (in either order), but this decomposition is not unique -- you can shift the rotation axis and compensate by performing a different translation. Thus you can only determine the direction of the axis from the data, not its position. You can choose the position arbitrarily, e.g. through the origin, and choose the translation accordingly. – joriki Aug 10 '12 at 14:33
  • Does this answer your question? [$n$-dimensional rotation along a 2D arbitrary plane](https://math.stackexchange.com/questions/501943/n-dimensional-rotation-along-a-2d-arbitrary-plane) – Michael Hoppe Apr 14 '20 at 18:10
  • @joriki should it be `((B-C)-(B'-C'))`? Also what do we do once we have the axis of rotation? – Ken Aug 06 '21 at 21:50

20 Answers20

145

Suppose you want to find a rotation matrix $R$ that rotates unit vector $a$ onto unit vector $b$.

Proceed as follows:

Let $v = a \times b$

Let $s = \|v\|$ (sine of angle)

Let $c = a \cdot b$ (cosine of angle)

Then the rotation matrix R is given by: $$R = I + [v]_{\times} + [v]_{\times}^2\frac{1-c}{s^2},$$

where $[v]_{\times}$ is the skew-symmetric cross-product matrix of $v$, $$[v]_{\times} \stackrel{\rm def}{=} \begin{bmatrix} \,\,0 & \!-v_3 & \,\,\,v_2\\ \,\,\,v_3 & 0 & \!-v_1\\ \!-v_2 & \,\,v_1 &\,\,0 \end{bmatrix}.$$

The last part of the formula can be simplified to $$ \frac{1-c}{s^2} = \frac{1-c}{1-c^2} = \frac{1}{1+c}, $$ revealing that it is not applicable only for $\cos(\angle(a, b)) = -1$, i.e., if $a$ and $b$ point into exactly opposite directions.

Nico Schlömer
  • 1,263
  • 8
  • 20
Jur van den Berg
  • 1,474
  • 1
  • 10
  • 2
  • 5
    I confirm that this works and gives answers identical to those in [my answer](http://math.stackexchange.com/a/897677/76513). – Kuba hasn't forgotten Monica Aug 14 '14 at 22:46
  • I tried to implement this, however I am getting an issue when the cross products is 0. – user2970916 Jan 14 '15 at 14:40
  • 1
    cross product can't be 0 because a and b are unit vectors – Julien__ May 23 '15 at 22:45
  • 7
    @Julien__ what if a == b? – Jerem May 30 '15 at 12:22
  • 4
    Indeed... I forgot that case. In that case v = 0 and the formula gives R = I. So here is a better answer : when the cross product is 0, R = I – Julien__ May 31 '15 at 13:08
  • @Julien__ Or -I if a==-b – user877329 Jun 01 '15 at 11:14
  • yes, the formula "doesn't work" if a == -b – Julien__ Jun 04 '15 at 08:48
  • So what to do in this case. -I turns the vector, but it flips all my normals. – user877329 Jun 15 '15 at 16:21
  • -I is not even a rotation matrix. If a == -b you can pick a rotation of pi about any axis perpendicular to a – Sean Dec 09 '15 at 16:42
  • Correct me if Im wrong, but do you assume we know the angle of rotation ? – Blue_Elephant Jul 01 '16 at 10:52
  • Can someone clarify the exceptions where this formula fails? Is it when c==0, and c ==Pi? Are there other instances? – jwarton Jul 29 '16 at 20:18
  • 1
    @jwarton c is the cosine of the angle between the vectors, and the function fails for parallel vectors, meaning when c == 1 and c == -1. – PfhorSlayer Jul 12 '17 at 18:50
  • 1
    @jwarton The reason the formula fails for anti-parallel vectors is that the rotation is not unique. Remember, the columns of R represent the 3 axes of the rotated coordinate system (expressed in terms of the original coordinates). If there is not a unique set of axes to represent the rotated coordinate system, that is when the formula fails. – Steven B. Segletes Jun 22 '18 at 13:28
  • 2
    I get the wrong answer when trying this: import numpy as np a = np.array([1, 0, 0], dtype=np.float64) b = np.array([0, 0, 1], dtype=np.float64) v = np.cross(a, b) s = np.linalg.norm(v) c = np.dot(a, b) vx = np.array([[0, -v[2], v[1]], [v[2], 0, -v[0]], [-v[1], v[0], 0]]) r = np.eye(3) + vx + vx*vx* (1-c)/(s**2) print(r) – user972014 Sep 12 '18 at 12:38
  • 8
    @user972014 note that you should use `np.dot` when squaring `vx` matrix, not a regular multiplication. – taketwo Feb 07 '19 at 10:49
  • 2
    Note, the formula above follows directly from Rodrigues' rotation formula. – Patwie Jun 14 '19 at 13:43
  • 1
    Do you apply this on the right or left? – user3180 Aug 25 '19 at 09:30
  • @Blue_Elephant No, he doesn't assume the rotation angle is known – Nathan Sep 14 '19 at 09:52
  • 2
    The Rodrigues formula states: $R = I + s[v]_\times + (1-c)[v]_\times^2$. Why is that the same as your solution? – kafman Aug 27 '20 at 17:51
  • 2
    Similar to @kafman I think the formula written here is wrong. The sentiment of the answer is correct but I think the author mixed up the angle with s. I have not found any reformulation of the rodriguez formula as stated here. – Marc HPunkt Oct 23 '20 at 13:34
  • Is the rotation matrix defined here in the active or passive sense? https://en.wikipedia.org/wiki/Active_and_passive_transformation – Sterling Dec 01 '20 at 16:17
  • I believe this is in the active sense. Ran a simple test using the example 2D matrix given in the wiki link above, and the matrix that is output seems to be consistent with the active interpretation, that is if a rotation by "theta" is defined as CCW in the active sense and we are looking at a right-handed coordinate system. See also [Rotation matrix](https://en.wikipedia.org/wiki/Rotation_matrix) – Sterling Dec 05 '20 at 06:23
  • For a more in depth explanation: https://www.theochem.ru.nl/%7Epwormer/Knowino/knowino.org/wiki/Rotation_matrix.html#Vector_rotation – Ian Zurutuza Apr 21 '21 at 07:21
66

Using Kjetil's answer answer, with process91's comment, we arrive at the following procedure.

Derivation

We are given two unit column vectors, $A$ and $B$ ($\|A\|=1$ and $\|B\|=1$). The $\|\circ\|$ denotes the L-2 norm of $\circ$.

First, note that the rotation from $A$ to $B$ is just a 2D rotation on a plane with the normal $A \times B$. A 2D rotation by an angle $\theta$ is given by the following augmented matrix: $$G=\begin{pmatrix} \cos\theta & -\sin\theta & 0 \\ \sin\theta & \cos\theta & 0 \\ 0 & 0 & 1 \end{pmatrix}.$$

Of course we don't want to actually compute any trig functions. Given our unit vectors, we note that $\cos\theta=A\cdot B$, and $\sin\theta=||A\times B||$. Thus $$G=\begin{pmatrix} A\cdot B & -\|A\times B\| & 0 \\ \|A\times B\| & A\cdot B & 0 \\ 0 & 0 & 1\end{pmatrix}.$$

This matrix represents the rotation from $A$ to $B$ in the base consisting of the following column vectors:

  1. normalized vector projection of $B$ onto $A$: $$u={(A\cdot B)A \over \|(A\cdot B)A\|}=A$$

  2. normalized vector rejection of $B$ onto $A$: $$v={B-(A\cdot B)A \over \|B- (A\cdot B)A\|}$$

  3. the cross product of $B$ and $A$: $$w=B \times A$$

Those vectors are all orthogonal, and form an orthogonal basis. This is the detail that Kjetil had missed in his answer. You could also normalize $w$ and get an orthonormal basis, if you needed one, but it doesn't seem necessary.

The basis change matrix for this basis is: $$F=\begin{pmatrix}u & v & w \end{pmatrix}^{-1}=\begin{pmatrix} A & {B-(A\cdot B)A \over \|B- (A\cdot B)A\|} & B \times A\end{pmatrix}^{-1}$$

Thus, in the original base, the rotation from $A$ to $B$ can be expressed as right-multiplication of a vector by the following matrix: $$U=F^{-1}G F.$$

One can easily show that $U A = B$, and that $\|U\|_2=1$. Also, $U$ is the same as the $R$ matrix from Rik's answer.

2D Case

For the 2D case, given $A=\left(x_1,y_1,0\right)$ and $B=\left(x_2,y_2,0\right)$, the matrix $G$ is the forward transformation matrix itself, and we can simplify it further. We note $$\begin{aligned} \cos\theta &= A\cdot B = x_1x_2+y_1y_2 \\ \sin\theta &= \| A\times B\| = x_1y_2-x_2y_1 \end{aligned}$$

Finally, $$U\equiv G=\begin{pmatrix} x_1x_2+y_1y_2 & -(x_1y_2-x_2y_1) \\ x_1y_2-x_2y_1 & x_1x_2+y_1y_2 \end{pmatrix}$$ and $$U^{-1}\equiv G^{-1}=\begin{pmatrix} x_1x_2+y_1y_2 & x_1y_2-x_2y_1 \\ -(x_1y_2-x_2y_1) & x_1x_2+y_1y_2 \end{pmatrix}$$

Octave/Matlab Implementation

The basic implementation is very simple. You could improve it by factoring out the common expressions of dot(A,B) and cross(B,A). Also note that $||A\times B||=||B\times A||$.

GG = @(A,B) [ dot(A,B) -norm(cross(A,B)) 0;\
              norm(cross(A,B)) dot(A,B)  0;\
              0              0           1];

FFi = @(A,B) [ A (B-dot(A,B)*A)/norm(B-dot(A,B)*A) cross(B,A) ];

UU = @(Fi,G) Fi*G*inv(Fi);

Testing:

> a=[1 0 0]'; b=[0 1 0]';
> U = UU(FFi(a,b), GG(a,b));
> norm(U) % is it length-preserving?
ans = 1
> norm(b-U*a) % does it rotate a onto b?
ans = 0
> U
U =

   0  -1   0
   1   0   0
   0   0   1

Now with random vectors:

> vu = @(v) v/norm(v);
> ru = @() vu(rand(3,1));
> a = ru()
a =

   0.043477
   0.036412
   0.998391
> b = ru()
b =

   0.60958
   0.73540
   0.29597
> U = UU(FFi(a,b), GG(a,b));
> norm(U)
ans =  1
> norm(b-U*a)
ans =    2.2888e-16
> U
U =

   0.73680  -0.32931   0.59049
  -0.30976   0.61190   0.72776
  -0.60098  -0.71912   0.34884

Implementation of Rik's Answer

It is computationally a bit more efficient to use Rik's answer. This is also an Octave/MatLab implementation.

ssc = @(v) [0 -v(3) v(2); v(3) 0 -v(1); -v(2) v(1) 0]
RU = @(A,B) eye(3) + ssc(cross(A,B)) + \
     ssc(cross(A,B))^2*(1-dot(A,B))/(norm(cross(A,B))^2)

The results produced are same as above, with slightly smaller numerical errors since there are less operations being done.

  • 1
    Matlab notation for Rik's Implementation: `function R=fcn_RotationFromTwoVectors(A, B) v = cross(A,B); ssc = [0 -v(3) v(2); v(3) 0 -v(1); -v(2) v(1) 0]; R = eye(3) + ssc + ssc^2*(1-dot(A,B))/(norm(v))^2;` – phyatt Aug 20 '15 at 21:55
  • @phyatt Doesn't the lambda syntax work in Matlab, too? – Kuba hasn't forgotten Monica Aug 20 '15 at 22:16
  • 1
    I've used lamdas in other languages, and I haven't seen it in use in the Matlab projects... and sure enough, its there. http://www.mathworks.com/help/matlab/matlab_prog/anonymous-functions.html Thanks @KuberOber. – phyatt Aug 20 '15 at 22:21
  • Rik's answer doesn't work if a == b, does your answer work in this case? – Ruts Nov 10 '15 at 08:53
  • 1
    @Ruts It can't work, since you can't normalize a zero-length vector `a-dot(a,a)*a`. You have to check for `a==b` first :) – Kuba hasn't forgotten Monica Nov 10 '15 at 17:06
  • A tiny point, I don't think it makes a difference but shouldn't $w$ be defined as $w = A\times B$ due to the fact that we are dealing with a right handed coordinate system and the way you defined $G$? – Dipole Sep 23 '16 at 16:58
  • The tests pass, so I think it's right as it is, but if you disagree I'm all ears. – Kuba hasn't forgotten Monica Sep 23 '16 at 17:11
  • who is "Rik" referring to? – Krupip Apr 18 '18 at 20:02
  • @snb Look for "Rik" within the answer. It links to an answer provided by someone who at some point called themselves "Rik" :) – Kuba hasn't forgotten Monica Apr 19 '18 at 04:10
  • Why w is normal? it's not always the case – user972014 Sep 12 '18 at 12:57
  • @user972014 `w` is normal to `B` and `A` both - as I understand, that's a property of the cross product in non-degenerate cases. If you can provide an example where `w` wouldn't be "normal" - I'll certainly reconsider! A vector is normal to some other object, such as another vector or a plane, so you'd have to mention that - otherwise what you say has no useful meaning. And just to be explicit: a plane is uniquely normal to either of a pair of opposite vectors, and any vector on that plane will also be normal to both. An infinite set of unit vectors is normal to any given unit vector. – Kuba hasn't forgotten Monica Sep 12 '18 at 13:13
  • The sentence said "Those vectors are all orthogonal and normal, and form an orthonormal basis". I thought that it meant that w is of unit length. – user972014 Sep 12 '18 at 15:25
  • I guess you could normalize it, but it doesn't seem to matter. It is an orthogonal and mostly normal basis then :) I probably used "normal" in an ambiguous fashion. – Kuba hasn't forgotten Monica Sep 12 '18 at 16:24
24

Rodrigues's rotation formula gives the result of a rotation of a vector $a$ about an axis of rotation $k$ through the angle $\theta$. We can make use of this by realizing that, in order to bring a normalized vector $a$ into coincidence with another normalized vector $b$, we simply need to rotate $a$ about $k=(a+b)/2$ by the angle $\pi$. With this, one gets the beautiful $$ R = 2 \frac{(a+b)(a+b)^T}{(a+b)^T(a+b)} - I. $$ This works in any dimension.

For the case $d=3$, one can derive a different matrix from Rodrigues's rotation formula, namely the rotation around the orthogonal $(a\times b) / \|a\times b\|$: $$ R_3 = \frac{1}{\|a\|\|b\|} \left(\langle a, b\rangle I + [a \times b]_\times + \frac{\|a\| \|b\| - \langle a, b\rangle}{\|a \times b\|^2} (a \times b) (a \times b)^T\right). $$ This matrix is the identity for $a = b$.

Nico Schlömer
  • 1,263
  • 8
  • 20
  • 2
    I can confirm that this works. Checked using Wolfram Mathematica that $R\cdot a=b$ with the assumptions that $||a||=||b||=1$. – Ruslan Apr 28 '18 at 19:30
  • It might be instructive to note that another way of thinking (that of course leads to a very similar result) are Householder reflections – Bort May 10 '18 at 19:24
  • 4
    While this will map vectors $a$ to $b$, it can add a lot of unnecessary "twist". For example, if $a=b=e_z$ this formula will produce a 180-degree rotation about the $z$-axis rather than the identity one might expect. – Alec Jacobson Mar 24 '19 at 19:54
  • I assume I is the identity matrix, but what is T? – omikes Apr 01 '19 at 17:46
  • @omikes Transpose. – Nico Schlömer Apr 01 '19 at 21:33
  • Could it be that the formula on Wikipedia at https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula#Matrix_notation lacks a $\cos \theta$ before the identity matrix? – XavierStuvw Mar 23 '20 at 10:17
  • Octave code snippet `k = (a / norm(a) + b / norm(b)) / 2; R = 2 * k * k' / (k' * k) - eye(3); det(R)` assuming that `a` and `b` are arbitrary 3D column vectors. This implements the top formula. – XavierStuvw Mar 23 '20 at 13:57
  • I might be late on this, but the first formula does not work in any dimension. Rotation around an axis only makes sense in 3 dimensions - so does Rodrigues's formula. – Matrefeytontias Aug 13 '20 at 07:56
  • 2
    @Matrefeytontias No, it really does work in any dimension. – Nico Schlömer Aug 13 '20 at 09:38
  • The first formula seems to produce an R that has the first two columns sign-reversed compared to other answers here. This still rotates a to b, but for other vectors, results in sign-reversed x and y elements of the result. Any ideas why? – jspencer Jul 15 '21 at 06:33
9

The MATLAB code for any dimension greater than one is

u = a/norm(a);                      % a and b must be column vectors
v = b/norm(b);                      % of equal length
N = length(u);
S = reflection( eye(N), v+u );      % S*u = -v, S*v = -u 
R = reflection( S, v );             % v = R*u

where

function v = reflection( u, n )     % Reflection of u on hyperplane n.
%
% u can be a matrix. u and v must have the same number of rows.

v = u - 2 * n * (n'*u) / (n'*n);
return

See this for background on how this works.

T L Davis
  • 378
  • 2
  • 7
8

From the top of my head (do the checking yourself)

  1. Let the given vectors in $R^3$ be $A$ and $B$. For simplicity assume they have norm 1 and are not identical.

  2. Define $C$ as the cross product of $A$ and $B$. We want an orthogonal matrix $U$ such that $UA=B$ and $UC=C$.

  3. First change bases into the new base $(U_1,u_2,u_3)=(A,B,C)$. In this new basis the matrix doing the job is simply $G=\left(\begin{smallmatrix} 0&1&0\\1&0&0\\0&0&1\end{smallmatrix}\right)$.

  4. Then we need the basis shift matrix to the new basis. Write the coordinates of the vectors in the old base as simply $A=(a_1,a_2,a_3), B=(b_1,b_2,b_3), C=(c_1,c_2,c_3)$. Then the basis shift matrix can be seen to be $\left( \begin{smallmatrix} a_1&b_1&c_1\\a_2&b_2&c_2\\a_3&b_3&c_3 \end{smallmatrix}\right)^{-1}$.

  5. The result is now simply $U=F^{-1} G F$, which is an orthogonal matrix rotating $A$ into $B$.

XavierStuvw
  • 298
  • 1
  • 12
kjetil b halvorsen
  • 6,117
  • 5
  • 35
  • 45
  • 1
    Actually, I think there may be further trouble - your matrix $U$ is only guaranteed to be orthogonal if $F$ is, which is only the case if $A$ and $B$ are orthogonal. In other instances, the angles of various other vectors to which $U$ is applied may be stretched or compressed. – mboratko Aug 09 '12 at 19:21
  • Your thought process is a good one, though - don't give up on it yet! I think there is a solution this way if you consider different original basis vectors which still have an easily represented form of $A$ and $B$. You could, for instance, choose the projection of $A$ on $B$, the rejection of $A$ on $B$, and $A \times B$. – mboratko Aug 09 '12 at 19:50
  • 5
    Just so that people aren't derailed: this answer doesn't work as-is, one has to take process91's comments into account. +1 since it provided valuable inspiration. – Kuba hasn't forgotten Monica Aug 14 '14 at 22:26
4

As you read from the thread: http://forums.cgsociety.org/archive/index.php/t-741227.html

Note that this will not align all 3 axes of the teapots, only the Z axis. We have no knowledge of the full orientation of the teapots in space, only know where the Z are pointing at. theTM will be the value you are looking for.

There is NO unique Matrix that could rotate one unit vector to another. Simply because the solution to 3 equations with 9 arguments does not unique.

Since you have the plane (not only the normal vector), a way to find a unique rotation matrix between two coordinate system would be: do the non-unique rotation twice!

That is

  1. Find a orthogonal vector in the same plane of interest with A and B respectively. Say $A_o \bot A$ in the original plane and $B_o \bot B$ in the target plane.
  2. use previous answers: This or this to find the rotation matrix within the $A\times B$ plane. Assume matrix is $U_{AB}$
  3. Apply the matrix $U_{AB}$ to $A_o$ result in a new $B'_o=U_{AB}A_o$.
  4. This $B'_o$ can be used as one of the new inputs to the same rotation funtion from the previous answers together with $B_o$. Another rotation matrix $U_{B'B}$ is calculated in $B'\times B$ plane.
  5. Since $B'\times B$ plane is perpendicular to $A$, $U_{B'B}$ would not influence the transform of $A\to B$. i.e. $B=U_{AB}A=U_{B'B}U_{AB}A$. and the final coordinate transform matrix should be $U=U_{B'B}U_{AB}$
  6. consider yourself for corner cases. :)

validate:

function:

% Implementation of Rik's Answer:
ssc = @(v) [0 -v(3) v(2); v(3) 0 -v(1); -v(2) v(1) 0];
RU = @(A,B) eye(3) + ssc(cross(A,B)) + ...
     ssc(cross(A,B))^2*(1-dot(A,B))/(norm(cross(A,B))^2);

Test:

> a=[1 0 0]'; b=[0 1 0]';
> ao=[0 1 0]'; bo=[-sqrt(2)/2,0,sqrt(2)/2]';
> Uab = RU(a,b);
> Ucd = RU(Uab*ao,bo);
> U = Ucd*Uab;
> norm(b-U*a) % does Ucd influence a to b?
ans = 0
> norm(bo-U*ao) % does Ucd influence a-orthogonal and b-orthogonal?
ans = 0
> U
U =

     0   -0.7071    0.7071
1.0000         0         0
     0    0.7071    0.7071
  • Please, use Latex. – ً ً Jul 15 '16 at 09:21
  • 1
    @MithleshUpadhyay I don't get it? I feel somewhat difficult when editing answers but I did use some LaTeX for equations. They looks nice on my screen. Is there anything I can do better? If there are tricks and tips I'm willing to hear from you. – Kangqiao Zhao Jul 19 '16 at 11:47
4

Use Rodrigues' rotation formula (See the section "Conversion to rotation matrix"). $\cos\theta$ is the dot product of the normalised initial vectors and $\sin\theta$ can be determined from $\sin^2\theta + \cos^2\theta =1$

Stefan Hansen
  • 24,191
  • 7
  • 51
  • 79
glennr
  • 141
  • 2
  • I don't think $\sin \theta$ can be determined in that way. There is sign ambiguity. (If it weren't so, then we'd get the same rotation matrix for going from $a$ to $b$ and viceversa.) – leonbloy Aug 08 '17 at 23:59
  • Could it be that the formula on Wikipedia at https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula#Matrix_notation lacks a cos before the identity matrix? I do not catch the explanation why this has been dropped. – XavierStuvw Mar 23 '20 at 10:30
3

The quaternion is a 4-dimensional complex number: http://en.wikipedia.org/wiki/Quaternion used to describe rotations in space. A quaternion (like a complex number) has a polar representation involving the exponential of the arguments (rotations), and a magnetude multiplier. Building the quaternion comes from the cross product (the product of the complex components), which will give you the argument in those 3 dimensions, you'll then get a number from that in the form A+Bi+Cj+Dk, and write it out in the matrix form described in the article there.

An easier way would be to simply fingure out what your original vectors are in the 4-space, and take the appropriate inverse operations to get your resultant quaternion (without going through the dot/cross product steps) but that requires a good foundation in hypercomplex algebra.

Vilid
  • 401
  • 3
  • 7
  • 1
    Just to be clear: from the article, the dot product will be your scalar part, and the cross product the vector part, so with dot product: x and scalar product iy + jz + ka the quaternion q would be q = x +iy +jz +ka. – Vilid Aug 08 '12 at 21:51
  • Not exactly. If the dot product is a, and the cross product is (x, y, z), it would be cos(a/2) + (x * sin(a/2))i + (y * sin(a/2))j + ( z * sin(a/2))k. Rotation quaternions are unitary. – Jeff Jul 19 '13 at 02:33
2

Here is a Matlab function that can be used to calculated the rotation from one vector to another.

Example 1:

>> v1=[1 2 3]';
>> v2=[4 5 6]';
>> fcn_RotationFromTwoVectors(v1, v2)

ans =

0.9789 0.0829 0.1870
 -0.0998 0.9915 0.0829
 -0.1785 -0.0998 0.9789

Example 2:

>> v1=[1 2 0]';
>> v2=[3 4 0]';
>> fcn_RotationFromTwoVectors(v1, v2)

ans =

0.9839 0.1789 0
 -0.1789 0.9839 0
 0   0     1.0000

Function:

function R=fcn_RotationFromTwoVectors(v1, v2)
% R*v1=v2
% v1 and v2 should be column vectors and 3x1

% 1. rotation vector
w=cross(v1,v2);
w=w/norm(w);
w_hat=fcn_GetSkew(w);
% 2. rotation angle
cos_tht=v1'*v2/norm(v1)/norm(v2);
tht=acos(cos_tht);
% 3. rotation matrix, using Rodrigues' formula
R=eye(size(v1,1))+w_hat*sin(tht)+w_hat^2*(1-cos(tht));

function x_skew=fcn_GetSkew(x)
x_skew=[0 -x(3) x(2);
 x(3) 0 -x(1);
 -x(2) x(1) 0];
Shiyu
  • 4,740
  • 3
  • 29
  • 41
2

Kuba Ober and Leyu Wang's answer works great. Here, is a python implementation of the same algorithm.

import numpy as np
import math


def rotation_matrix(A,B):
# a and b are in the form of numpy array

   ax = A[0]
   ay = A[1]
   az = A[2]

   bx = B[0]
   by = B[1]
   bz = B[2]

   au = A/(np.sqrt(ax*ax + ay*ay + az*az))
   bu = B/(np.sqrt(bx*bx + by*by + bz*bz))

   R=np.array([[bu[0]*au[0], bu[0]*au[1], bu[0]*au[2]], [bu[1]*au[0], bu[1]*au[1], bu[1]*au[2]], [bu[2]*au[0], bu[2]*au[1], bu[2]*au[2]] ])


   return(R)
Soumya
  • 21
  • 1
2

You could say you are looking for a transformation between two orthonormal bases:

$$ M*[\vec{i},\vec{j},\vec{k}]=[\vec{i}',\vec{j}',\vec{k}'] $$ where

  • $\vec{i}$ and $\vec{i}'$ are the two vectors in question ("from" and "to")

and the missing parts can be picked in a way which suits the purpose:

  • $\vec{j}=\vec{j}'=\frac{\vec{i}\times\vec{i}'}{||\vec{i}\times\vec{i}'||}$, the axis of rotation, which is left in place and is orthogonal to the $\vec{i},\vec{i}'$ plane (and has to be normalized)
  • $\vec{k}=\vec{i}\times\vec{j}, \vec{k}'=\vec{i}'\times\vec{j}'$, because you need a pair of 3rd vectors for completing the bases.

Since $[\vec{i},\vec{j},\vec{k}]$ is an orthogonal matrix, there is no need for 'real' inversion and the transformation is

$$ M=[\vec{i}',\vec{j}',\vec{k}']*[\vec{i},\vec{j},\vec{k}]^T $$ Except when $\vec{i},\vec{i}'$ do not stretch a plane (and thus $||\vec{i}\times\vec{i}'||=0$). This calculation will not produce $I$, e.g. dies on both $\vec{i}=\vec{i}'$ and $\vec{i}=-\vec{i}'$

Disclaimer: Sorry about the necro, and especially for the partial repeat: I see Kjetil's answer, but I simply do not understand what and why that skewed matrix is doing there, and while Kuba's answer says it builds on Kjetil's, it introduces trigonometry on top of that, slightly defeating the idea (of course I understand that the trigonometric part is expressed with dot/cross products at the end)

Plus this was my first time with LaTeX, I feel $\vec{i}~\vec{i}'$ look a bit too similar, but $\vec{i}~\vec{i'}$ is just plain ugly. And writing that fraction properly is way above me.

tevemadar
  • 216
  • 1
  • 8
1

This is an interesting problem. It guarantees that you rotate the unit vector $a$ to coincide with $b$. Still it would not be enough to rotate a rigid body to make it coincide with another. I post an example of this below but first let me point a few important references:

  1. @glennr correctly pointed out the Rodrigues rotation formula . Note that this formula is the same as the one shown by @Jur var den vec here. The scaling factors are due to normalization of the cross and dot products since the vector $v$ in the Wikipedia website (here $a$) is arbitrary and here $a$ is unit. The Wikipedia page only requires vector manipulations (need to be good at taking cross products).
  2. Here is an elegant proof that does not depend too much on geometrical constructions but it requires simple linear ordinary differential equations for matrices, in addition to the cross product representation as a matrix operator.

  3. The final reference is the most complete document that I have found on this topic. It is titled ROTATION: A review of useful theorems involving proper orthogonal matrices referenced to three dimensional physical space.

Now, for an example. This stack exchange link is about a question I posted on the rotation of a tetrahedron. I found and answer to my question by two consecutive rotations of matrices. The apex of the tetrahedron was at $(0,0,\sqrt(3))$ and I wanted it to be at $(-1,-1,1)$. Of course I wanted to verify that all vertices were rotated to their correct place and the product of matrices: \begin{eqnarray} Y = \left ( \begin{array}{ccc} \cos \alpha & -\sin \alpha & 0 \\ \sin \alpha & \cos \alpha & 0 \\ 0 & 0 & 1 \end{array} \right ) = \left ( \begin{array}{ccc} \frac{\sqrt{2}}{2} & -\frac{\sqrt{2}}{2} & 0 \\ \frac{\sqrt{2}}{2} & \frac{\sqrt{2}}{2} & 0 \\ 0 & 0 & 1 \end{array} \right ) \end{eqnarray}

and \begin{eqnarray} P = \left ( \begin{array}{ccc} \cos \theta & 0 & -\sin \theta \\ 0 & 1 & 0 \\ \sin \theta & 0 & \cos \theta \end{array} \right ) = \left ( \begin{array}{ccc} \frac{\sqrt{3}}{3} & 0 & -\frac{\sqrt{2}}{\sqrt{3}} \\ 0 & 1 & 0 \\ \frac{\sqrt{2}}{\sqrt{3}} & 0 & \frac{\sqrt{3}}{3} \\ \end{array} \right ) \end{eqnarray}

That is the matrix: \begin{eqnarray*} \frac{1}{6} \left ( \begin{array}{ccc} \sqrt{6} & -3\sqrt{2} & - 2 \sqrt{3} \\ \sqrt{6} & 3 \sqrt{2} & -2 \sqrt{3} \\ 2 \sqrt{6} & 0 & 2 \sqrt{3} \end{array} \right ) \end{eqnarray*}

did the job.

I thought that the Rodrigues rotation (the one being considered here) would do the job but it did not. I computed the Rodrigues rotation for this problem. Here is the result:

\begin{eqnarray*} R = \left ( \begin{array}{ccc} \frac{\sqrt{3}+1}{2 \sqrt{3}} & -\frac{\sqrt{3}-1}{2 \sqrt{3}} & -\frac{1}{\sqrt{3}} \\ -\frac{\sqrt{3}-1}{2 \sqrt{3}} & \frac{\sqrt{3}+1}{2 \sqrt{3}}& -\frac{1}{\sqrt{3}} \\ \frac{1}{\sqrt{3}} & \frac{1}{\sqrt{3}} & \frac{1}{\sqrt{3} } \end{array} \right ) \end{eqnarray*}

This matrix correctly maps the vector $(0,0,1)$ into the vector $\frac{1}{\sqrt{3}} (-1,-1,1)$ as promised, but the other three vectors on the base of the pyramid (tetrahedron) are not correctly mapped into their positions.

The meaning of this is that, to map rigid bodies through rotations, there is not a magic formula where only one matrix could be found and instead a multiplication of elementary (or non elmentary) matrices is required to do the job, unless we are lucky. Right?

XavierStuvw
  • 298
  • 1
  • 12
Herman Jaramillo
  • 2,568
  • 21
  • 25
1

Given two unit vectors $\hat a$ and $\hat c$, reflecting a vector $x$ across the orthogonal complement of $\hat a$ and then for $\hat c$ will rotate the part of $x$ in the span of $\hat a$ and $\hat c$ by twice the angle from $\hat a$ to $\hat c$. Letting $\hat c = \frac{\hat a + \hat b}{|\hat a+\hat b|}$ be the unit vector that bisects $\hat a$ and $\hat b$, the composition of the two reflections $R_{\hat c} \circ R_{\hat a}$ will rotate $\hat a$ to $\hat b$, and any vector $x$ by the angle from $\hat a$ to $\hat b$. Recall that the reflection transformation is $R_n(x) = x - 2 \frac{x \cdot n}{n \cdot n} n$.

Centrinia
  • 51
  • 1
1

General solution for n dimensions in matlab / octave:

%% Build input data n = 4; a = randn(n,1); b = randn(n,1);

%% Compute Q = rotation matrix A = a*b'; [V,D] = eig(A'+A); [~,idx] = min(diag(D)); v = V(:,idx); Q = eye(n) - 2*(v*v');

%% Validate Q is correct b_hat = Q'*a*norm(b)/norm(a); disp(['norm of error = ' num2str(norm(b_hat-b))]) disp(['eigenvalues of Q = ' num2str(eig(Q)')])

NoseKnowsAll
  • 2,001
  • 14
  • 22
Orion
  • 11
  • 2
1

Sadly I don't have enough points to comment on the accepted answer but as others have noted, the formula doesn't work when a == -b.

To solve this edge case you have to create a normal vector of a by for example using the formula found here (a,b and c being the components of the vector):

function (a,b,c)  
{
    return  c<a  ? (b,-a,0) : (0,-c,b)  
}

then make the rotation matrix by rotating vector a around this normal by Pi.

jaga
  • 111
  • 1
1

You can easily do all this operation using the Vector3 library.

The following four steps worked for me.

Vector3D axis = Vector3D.CrossProduct(v1, v2);

if (axis.Magnitude != 0)
{
    axis.Normalize();
    double AngleFromZaxis = Vector3D.AngleBetween(new Vector3D(0, 0, 1), vAxis);
    Vector3D vAxis = new Vector3D(axis.X, axis.Y, axis.Z);
    Matrix3D m = Matrix3D.Identity;
    Quaternion q = new Quaternion(vAxis, AngleFromZaxis);

    m.RotateAt(q, centerPoint);
    MatrixTransform3D mT = new MatrixTransform3D(m);

    group.Children.Add(mT);

    myModel.Transform = group;
}
Sen Jacob
  • 663
  • 9
  • 11
0

One way to proceed is as following:

Start by constructing one orthonormal basis for each of the vectors $\vec{n_1}$ and $\vec{n_2}$. This can be done by the trick given in an answer to another question [1]. This will result in two transformation matrices

$$R_1=\begin{bmatrix}\vec{u_1} & \vec{v_1} & \vec{w_1}\end{bmatrix}$$

and $$R_2=\begin{bmatrix}\vec{u_2} & \vec{v_2} & \vec{w_2}\end{bmatrix}$$

Then the rotation matrix for aligning $\vec{n_1}$ onto $\vec{n_2}$ becomes

$$R=R_2{R_1}^T\vec{n_1}$$

[1] https://math.stackexchange.com/q/712065

user877329
  • 187
  • 11
0

Here's how to find the transformation from one triangle to another.

First triangle has vertices $a,b,c$ and normal $n$ and second triangle has vertices $a',b',c'$ and normal $n'$.

First we will find transformation $f(x)=Mx+t$ from reference triangle to the first triangle and another transformation $g(x)=M'x+t'$ from reference triangle to the second triangle. Then the transformation mapping first triangle to the second is $$T(x) = g(f^{-1}(x)) = M'(M^{-1}(x-t)) + t' = Rx + s$$ where $R = M'M^{-1}$ and $s = -M'M^{-1}t + t'$.

As the reference triangle we will use triangle with vertices $(0,0,0),(1,0,0),(0,1,0)$ and normal $(0,0,1)$

Then: $$M = [ b-a, c-a, n]$$ $$t = a $$ $$M' = [ b'-a', c'-a', n' ] $$ $$t' = a' $$

You have to be cautious and give vertices $a,b,c$ and $a',b',c'$ in right order. This has to satisfy $((b-a)\times (c-a))\cdot n > 0$ and the same for second triangle.


If those two triangles are not the same then matrix $R$ will not be orthogonal. But we can find isometry which maps one triangle to the another as close as possible in some sense. For this you can use Kabsch algorithm which is well explained here.

meetar
  • 117
  • 5
tom
  • 4,488
  • 1
  • 17
  • 37
-1

This paper shows detail how to solve this problem. Assume we want to rotate vector f to vector t,

  1. Let v = f x t, u = v/||v||, c = f . t, h = 1-c/1-c^2
  2. The formula of the rotation is \begin{bmatrix} c + hv_x^2 & hv_xv_y-v_z & hv_xv_z+v_y \\ hv_xv_y+v_z & c+hv_y^2 & hv_yv_z-v_x \\ hv_xv_z-v_y & hv_yv_z+v_x & c+hv_z^2 \end{bmatrix}

Python code:

import numpy as np
v = np.cross(f, t)
u = v/np.linalg.norm(v)
c = np.dot(f, t)
h = (1 - c)/(1 - c**2)

vx, vy, vz = v
rot =[[c + h*vx**2, h*vx*vy - vz, h*vx*vz + vy],
      [h*vx*vy+vz, c+h*vy**2, h*vy*vz-vx],
      [h*vx*vz - vy, h*vy*vz + vx, c+h*vz**2]]
  • Note that the definition ```u = v/||v||``` in bullet 1 and the corresponding statement in the Python code ```u = v/np.linalg.norm(v)``` are redundant as $u$ is never used. It is defined in the original paper as preparation for some special-case handling later on. – Pedro Nov 07 '21 at 13:23
-1

I have a simpler method comes from Erigen's "Mechanics of Continua". 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)]

To verify:

c=R*a
cu=c./sqrt(c(1)^2+c(2)^2+c(3)^2)
simple(bu-cu)

A zero result means that $c$ (rotated $a$) 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))

A zero result means that $c$ (rotated $a$) and $a$ are of the same length.

Leyu Wang
  • 17
  • 1