Matlab Script for Rodrigues' rotation formula

Rotation matrix is a 3x3 unitary matrix which rotates one 3D vector to another. Assuming two unit 3D vectors k and v and their angle \theta, the rotation matrix is calculated by using Rodrigues' formula (refer to wiki for more details http://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula):


Here is an example of Matlab script of how to calculate the matrix:

% This script is used to calculate rotation matrix based on Rodrigues formula

clear all; close all;

% acc measurement

in0 = [-1.213062233 0.472455803 9.666190333];

in1 = [0 0 1];

% Reference

in0Norm = in0 ./ sqrt(in0*in0.');

in1Norm = in1 ./ sqrt(in1*in1.');

axisRaw = cross(in0Norm, in1Norm);

axis = axisRaw./sqrt(axisRaw*axisRaw.');

angle = acos(dot(in0Norm, in1Norm));

M = [axis(1)^2+(1-axis(1)^2)*cos(angle) axis(1)*axis(2)*(1-cos(angle))-axis(3)*sin(angle) axis(1)*axis(3)*(1-cos(angle))+axis(2)*sin(angle);

axis(1)*axis(2)*(1-cos(angle))+axis(3)*sin(angle) axis(2)^2+(1-axis(2)^2)*cos(angle) axis(2)*axis(3)*(1-cos(angle))-axis(1)*sin(angle);

axis(1)*axis(3)*(1-cos(angle))-axis(2)*sin(angle) axis(2)*axis(3)*(1-cos(angle))+axis(1)*sin(angle) axis(3)^2+(1-axis(3)^2)*cos(angle)];

In this script, in0 and in1 are two 3D vectors. One example is that in0 and in1 are readings of accelerometer for a smart phone at different orientations. The output, matrix M, is the rotation matrix.

[Added on 07/30/2016] Rodrigues' rotation formula can be also understand in this way. The rotation matrix R is:
R = I + sin(angle)*K + (1-cos(angle))*K*K
where K is the matrix for cross product as K = [0 -axis(3) axis(2); axis(3) 0 -axis(1); -axis(2) axis(1) 0]. How to interpret this formula? For the input vector of in0, K*K*in0 becomes -in0, and K*in0 is a vector perpendicular to in0. Thus, R*in0 become sin(angle)*in0_perpendicular+cos(angle)*in0_parallel.