Rodrigues' rotation formula

In many areas of robotics, computer vision and graphic is a need to rotate points or objects by certain angle around a given axis. One way to do it is to use Rodrigue's rotation formula.

R=I+[n]_\times \sin\alpha+[n]_\times^2(1-\cos\alpha)

Where the [n]_\times is a skew symmetric of normalized vector of the rotation axis and \alpha is the rotation angle in radians. The result is a 3x3 rotation matrix R

It's simple to use and implement in every programming language but what's under the bonnet? How does it work? Ok lets work it out. There are actually at least two ways how to derive this formula, one purely geometric and the other more mathematical based on the matrix exponential. I will describe the geometric approach first.

So let's start from the beginning.  The Rodrigues formula is used to perform angle-axis rotation, or simply said it rotates a point in the space around a given axis by given angle. It also outputs standard rotation matrix that I'm sure most people reading this article are familiar with. The most basic way of representing a rotation are Euler Angles.

Euler Angles

Any rotation in the space is represented by three 2D rotations in planes xy,  yz, and xz. The matrix form of these rotation is

    \[A=\begin{pmatrix}\cos\phi & \sin\phi & 0\\-\sin\phi & \cos\phi & 0\\0 & 0 & 1\end{pmatrix} \]

    \[B=\begin{pmatrix}1 & 0 & 0\\0 & \cos\theta & \sin\theta\\0 & -\sin\theta & \cos\theta\end{pmatrix}\]

    \[C=\begin{pmatrix}\cos\theta & 0 & -\sin\theta \\0 & 1 & 0\\\sin\theta & 0 &\cos\theta\\\end{pmatrix}\]

The rotation matrix is a product of the three discrete 2D rotations.

R=ABC

It's clear the result depends on the order in which we perform the rotations. Rotating around X->Y->Z axes results in different final position than X->Z->Y. This is also given by that fact matrix multiplication is not commutative. The Euler angles rotation is simple to understand and visualize however is not really practical in many applications. Among problems associated with this technique is that the rotation is not smooth (a small change in position of the object can result in dramatic change of one or more Euler angles) and there is also risk of gimbal lock. Angle-axis rotation solves these problems and on top of that is also very simple to understand and work with.

Angle-axis

Angle-axis rotation is more common the one may think. For instance everything on the Earth is rotation around the Earth's axis which is not the same as the solar system or our galaxy axes. Other example would be a car. Wheels are rotating around their axes that are different from the axis of the steering wheel or the engine crankshaft. Let's look at some more mathematical description.

Let say we have a point V that we want to rotate by angle \alpha around axis n to the final position U

To do that we first define a vector v from the origin to the point V and a plane in the origin perpendicular to the axis of rotation n We also normalize the vector representing the axis of rotation |n|=1 Second step is to project the vector v on the axis of rotation and also on the plane in the origin. The vector projected on the rotation axis n is called v_p We will calculate it as a vector projection which is basically a dot product of the vector v with the unit vector n in the direction of n

v_p=(n\cdot v)n

To get the vector v_r we simply subtract v_p from v

v_r=v-v_p

The vector v_r is called vector rejection Now we create a vector w perpendicular to the plane v,v_p We will use this vector to rotate the image of the v_r by the angle \alpha on the plane in the origin.

w=n\times v

It turns out that the cross product between unit vector n and the vector v has the same length as the vector rejection v_r To do that we perform a 2D rotation of the vector v_r in the plane v_r,w. We will call this new vector v_{rr}

v_{rr}=v_r\cos\alpha+w\sin\alpha We are almost done. There's one more step. To get from the plane in the origin to the point U we just need to add vectors v_{rr} and v_p

u=v_{rr}+v_p Let's recap all equations and pull them together.

v_p=(n\cdot v)n
v_r=v-v_p
w=n\times v
v_{rr}=v_r\cos\alpha+w\sin\alpha
u=v_{rr}+v_p

If we substitute all equations we get
u=(v - (n\cdot v)n)\cos\alpha + (n\times v) \sin\alpha +(n\cdot v)n u=v\cos\alpha-(n\cdot v)n\cos\alpha + (n\times v) \sin\alpha + (n\cdot v)n u=v\cos\alpha+(n\cdot v)n(1-\cos\alpha) + (n\times v) \sin\alpha

We can simplify this equation by converting the dot and cross products to multiplication. Assuming that all vectors are column vectors the dot product is calculated as:

v_p=(v\cdot n)n = (n^Tv)n = nn^Tv

To calculate the cross product we use a skew-symmetric matrix trick 🙂 We first convert the n vector to the skew-symmetric matrix and then multiply it with the vector v The vector n in the skew-symmetric matrix form is

[n]_\times=\begin{pmatrix}0 & -n_3 & n_2\\n_3 & 0 & -n_1\\-n_2 & n_1 & 0\end{pmatrix}

w=n\times v = [n]_\times v

Substituting the dot and cross product equations we get the Rodrigues formula.

u=v\cos\alpha+nn^Tv(1-\cos\alpha) + [n]_\times v\sin\alpha

So now we have a way to rotate point v around the rotation axis n by angle \alpha. To get the rotation matrix R we need to split this equation in the R and v components.

u=Rv

thus
u=(I\cos\alpha+nn^T(1-\cos\alpha) + [n]_\times \sin\alpha)v
R=I\cos\alpha+nn^T(1-\cos\alpha) + [n]_\times \sin\alpha

We can still simplify it a bit. We use the identity of the outer product

nn^T=[n]_\times^2+I

Where I is the identity matrix.

This identity can be proved by following  a few steps. First of all lets write down nn^T in the term of element-wise operations.

{nn}^T=\begin{bmatrix}n_1 \\ n_2 \\ n_3\end{bmatrix}\begin{bmatrix}n_1 & n_2 & n_3\end{bmatrix}=\begin{bmatrix}{n_1}^2 & n_1n_2 & n_1n_3\\n_1n_2 & {n_2}^2 &n_2n_3\\n_1n_3 & n_2n_3 &{n_3}^2\end{bmatrix}

With a little bit of insight we may realize this is very close to the skew symmetric matrix of the n squared [n]_\times^2

[n]_\times=\begin{bmatrix}0 & -n_3 & n_2\\n_3 & 0 &-n_1\\-n_2 & n_1 & 0\end{bmatrix}

If we square it we get

[n]_\times^2=\begin{bmatrix}-{n_3}^2-{n_2}^2 & n_1n_2 & n_1n_3\\n_1n_2 & -{n_3}^2-{n_1}^2&n_2n_3\\n_1n_3 & n_2n_3 &-{n_2}^2-{n_1}^2\end{bmatrix}

Now we can see that it really looks like  nn^T apart from the diagonal. The last step is to somehow turn the diagonal of the [n]_\times^2 to the diagonal of nn^T This is quite easy. Since the vector n has unit length we can just add one to the [n]_\times^2 diagonal to get nn^T For example

    \[{n_1}^2+{n_2}^2+{n_3}^2=1 \quad \Rightarrow \quad {n_1}^2=-{n_2}^2-{n_3}^2+1\]

So now we can write

\begin{bmatrix}{n_1}^2 & n_1n_2 & n_1n_3\\n_1n_2 & {n_2}^2 &n_2n_3\\n_1n_3 & n_2n_3 &{n_3}^2\end{bmatrix}=\begin{bmatrix}-{n_3}^2-{n_2}^2&n_1n_2&n_1n_3\\n_1n_2 &-{n_3}^2-{n_1}^2&n_2n_3\\n_1n_3&n_2n_3 &-{n_2}^2-{n_1}^2\end{bmatrix}+I

which is

{nn}^T =[n]_\times^2+I

Using this identity in our equation we get

R=I\cos\alpha+([n]_\times^2+I)(1-\cos\alpha) + [n]_\times\sin\alpha
R=I\cos\alpha+[n]_\times^2+I-[n]_\times^2\cos\alpha-I\cos\alpha+[n]_\times\sin\alpha
R=I+[n]_\times^2-[n]_\times^2\cos\alpha+[n]_\times\sin\alpha

And finally we get the Rodriques formula
R=I+[n]_\times\sin\alpha+[n]_\times^2(1-\cos\alpha)

Next time I will show how the same formula is derived from angular velocity of  a rigid object.

This entry was posted in Computer Vision. Bookmark the permalink.

4 Responses to Rodrigues' rotation formula

  1. machinegunmax says:

    Fantastic proof. Took me a while to get through, but I got there in the end. Thanks for doing this!!!

  2. skysoup says:

    Thanks for this! The pictures really help with visualizing things.
    I understand almost all of it except for 2 things in the simplification part .

    1.) For this part, vp=(v⋅n)n=(nTv)n=nTnv (converting the dot and cross products to multiplication)
    How were you able to swap the order of (nTv)n to nTnv?
    I thought matrix multiplication was non-commutative? (Since vectors are just matrices with one of its dimensions being 1)
    It also appears differently in the next step:
    u=vcosα+nnTv(1−cosα)+[n]×vsinα (nnTv != nTnv)

    2.) For this part, nnT=[n]2×+I (Taking the identity of the outer product)
    I'm not too clear on how this identity holds true.
    I tried giving concrete values to n1=1, n2=2, and n3=3 of the vector n and I'm getting different results on the 2 sides of the equation.

    Once again, thanks for this awesome explanation!

    • admin says:

      Hi skysoup,
      Sorry for this late reply and BIG thanks for your comment and spotting my mistake. The (v⋅n)n is indeed nnTv and not nTnv.

      Regarding the outer product identity I have update the article with the proof. The identity holds true only for normalized vector n which is probably why it didn't work for you. In the Matlab/Octave you can try this.
      a=[1;2;3];
      b=a/norm(a);
      b*b'
      b_skew = [ 0 -b(3) b(2) ; b(3) 0 -b(1) ; -b(2) b(1) 0];
      b_skew^2 + eye(3)

Leave a Reply