Vtyjkyuk

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=beginpmatrix
costheta & -sintheta & 0 \
sintheta & costheta & 0 \
0 & 0 & 1
endpmatrix.$$

Of course we don't want to actually compute any trig functions. Given our unit vectors, we note that $costheta=Acdot B$, and $sintheta=||Atimes B||$. Thus
$$G=beginpmatrix
Acdot B & -|Atimes B| & 0 \
|Atimes B| & Acdot B & 0 \
0 & 0 & 1endpmatrix.$$

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$$




2. normalized vector rejection of $B$ onto $A$: $$v=$$




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

Those vectors are all orthogonal and normal, and form an orthonormal basis. This is the detail that Kjetil had missed in his answer.

The basis change matrix for this basis is:
$$F=beginpmatrixu & v & w endpmatrix^-1=beginpmatrix A & & B times Aendpmatrix^-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^-1G 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,0right)$ and $B=left(x_2,y_2,0right)$, the matrix $G$ is the forward transformation matrix itself, and we can simplify it further. We note
$$beginaligned
costheta &= Acdot B = x_1x_2+y_1y_2 \
sintheta &= | Atimes B| = x_1y_2-x_2y_1
endaligned$$

Finally,
$$Uequiv G=beginpmatrix
x_1x_2+y_1y_2 & -(x_1y_2-x_2y_1) \
x_1y_2-x_2y_1 & x_1x_2+y_1y_2
endpmatrix$$
and
$$U^-1equiv G^-1=beginpmatrix
x_1x_2+y_1y_2 & x_1y_2-x_2y_1 \
-(x_1y_2-x_2y_1) & x_1x_2+y_1y_2
endpmatrix$$

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 $||Atimes B||=||Btimes 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.

At the top of my head (do the checking yourself)
Let the given vectors in $R^3$ be $A$ and $B$, for simplicity assume they have norm 1, and assume they are not identical. Define $C$ as the cross product of $A$ and $B$. We want an orthogonal matrix $U$ such that $UA=B$ and $UC=C$.
First change bases, to 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(beginsmallmatrix 0&1&0\1&0&0\0&0&1endsmallmatrixright)$. 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( beginsmallmatrix a_1&b_1&c_1\a_2&b_2&c_2\a_3&b_3&c_3 endsmallmatrixright)^-1$. The result is now simply $U=F^-1 G F$, which is an orthogonal matrix rotating $A$ into $B$.

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.

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.

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

Rodrigues's rotation formula gives the result of a rotation of a vector $a$ around a rotation axis $k$ by the angle $theta$. We can make use of this by realizing that, to rotate a unit vector $a$ into $b$, we simply need to rotate $a$ by $pi$ around $(a+b)/2$. With this, one gets the beautiful
$$
R = 2 frac(a+b)(a+b)^T(a+b)^T(a+b) - I.
$$

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];

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 $Atimes B$ plane. Assume matrix is $U_AB$


3. Apply the matrix $U_AB$ to $A_o$ result in a new $B'_o=U_ABA_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 $Ato B$. i.e. $B=U_ABA=U_B'BU_ABA$. and the final coordinate transform matrix should be $U=U_B'BU_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

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

The following four steps worked for me.

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

if (axis.Magnitude != 0)

 axis.Normalize();
 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;

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)')])

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

$$
M*[veci,vecj,veck]=[veci',vecj',veck']
$$
where

• $veci$ and $veci'$ are the two vectors in question ("from" and "to")

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

• $vecj=vecj'=fracvecitimesveci'$, the axis of rotation, which is left in place and is orthogonal to the $veci,veci'$ plane (and has to be normalized)


• $veck=vecitimesvecj, veck'=veci'timesvecj'$, because you need a pair of 3rd vectors for completing the bases.

Since $[veci,vecj,veck]$ is an orthogonal matrix, there is no need for 'real' inversion and the transformation is

$$
M=[veci',vecj',veck']*[veci,vecj,veck]^T
$$
Except when $veci,veci'$ do not stretch a plane (and thus $||vecitimesveci'||=0$). This calculation will not produce $I$, e.g. dies on both $veci=veci'$ and $veci=-veci'$

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 $veci~veci'$ look a bit too similar, but $veci~veci'$ is just plain ugly. And writing that fraction properly is way above me.

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^-1t + 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.

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 will post an example of this below but first let me point a few important references. @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). 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. 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 proer orthogonal matries 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:
begineqnarray Y = left ( beginarrayccc cos alpha &
-sin alpha & 0 \ sin alpha & cos alpha & 0 \
0 & 0 & 1 endarray right )
= left ( beginarrayccc fracsqrt22 & -fracsqrt22 & 0 \ fracsqrt22 & fracsqrt22 & 0 \
0 & 0 & 1 endarray right )
endeqnarray

and
begineqnarray P = left (
beginarrayccc cos theta & 0 & -sin theta \ 0 & 1 & 0
\ sin theta & 0 & cos theta endarray right )
= left ( beginarrayccc fracsqrt33 & 0 & -fracsqrt2sqrt3 \
0 & 1 & 0 \ fracsqrt2sqrt3 & 0 & fracsqrt33 \ endarray right )
endeqnarray

That is the matrix:
begineqnarray*
frac16 left (
beginarrayccc
sqrt6 & -3sqrt2 & - 2 sqrt3 \
sqrt6 & 3 sqrt2 & -2 sqrt3 \
2 sqrt6 & 0 & 2 sqrt3
endarray
right )
endeqnarray*

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:

begineqnarray*
R =
left (
beginarrayccc
fracsqrt3+12 sqrt3 & -fracsqrt3-12 sqrt3 & -frac1sqrt3 \
-fracsqrt3-12 sqrt3 & fracsqrt3+12 sqrt3& -frac1sqrt3 \
frac1sqrt3 & frac1sqrt3 & frac1sqrt3
endarray
right )
endeqnarray*

This matrix correctly maps the vector $(0,0,1)$ into the vector
$frac1sqrt3 (-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?

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)

One way to proceed is as following:

Start by constructing one orthonormal basis for each of the vectors $vecn_1$ and $vecn_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=beginbmatrixvecu_1 & vecv_1 & vecw_1endbmatrix$$

and
$$R_2=beginbmatrixvecu_2 & vecv_2 & vecw_2endbmatrix$$

Then the rotation matrix for aligning $vecn_1$ onto $vecn_2$ becomes

$$R=R_2R_1^Tvecn_1$$

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

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 = frachat 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 fracx cdot nn cdot n n$.

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.

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.