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:
normalized vector projection of $B$ onto $A$: $u={(A\cdot B)A \over \|(A\cdot B)A\|}=A$
normalized vector rejection of $B$ onto $A$: $v={B-(A\cdot B)A \over \|B- (A\cdot B)A\|}$
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.