Introduction to Quaternions for 3D Rotation Computations

History and Motivation

Why do we care about quaternions? Short answer: 3D rotations. Longer answer: for centuries, mathematicians have tried to find new, simpler algorithms to perform geometry computations. A common problem is to compute rotations of geometric figures.

Two Dimensions

For example, in two dimensions, suppose you define a triangle by the Cartesian coordinates of its vertices, such as $(0,0)$ , $(0,1)$ , and $(1,0)$ , which gives us an isosceles right triangle whose legs have length $1$ and whose hypotenuse has length $\sqrt{2}$ .

If one rotates it $90$ degrees counterclockwise about the point $(0,0)$ , the transformed triangle now has vertices $(0,0)$ , $(-1,0)$ , and $(0,1)$ . How would we compute that? One method is to use trigonometry; the method can be summarized as: convert the points to polar form, add $90$ degrees ( $\pi/2$ radians) to the angles, and then convert the points back to Cartesian form. With a little algebra, one can collapse these three steps into one formula.

If you are like me, you don’t like to memorize such formulas, and you are looking for something mathematically “natural”. It turns out that one can use complex numbers: the points of the (original, unrotated) triangle are now the three complex numbers $0$ , $i$ , and $1$ . To rotate by $90$ degrees counterclockwise, you only need to multiply the numbers by $i$ , giving $0$ , $-1$ , and $i$ . This is easier to remember. How did we know that $i$ is $90$ degrees counterclockwise? Well, we still need some trigonometry for that: convert $(1,0)$ , or in complex numbers, $1$ , to polar form, rotate by $90$ degrees counterclockwise by adding $90$ to the angle, then convert back to the complex number $i$ . Unlike before, we only need to do this once instead of 3 times. We then multiply the complex representations of the three points by the result $i$ . This is really advantageous for more complex figures, say a sprite with 100 vertices.

Three dimensions

However, we care about 3D rotations. How would we do this in three dimensions? This problem was the lifelong focus of the Irish mathematician William Rowan Hamilton. He tried to find a three-dimensional version of complex numbers so he could rotate figures in three dimensions by multiplying numbers instead of by converting 3D Cartesian coordinates to, for example, spherical coordinates, adjusting the angles, and then converting back to Cartesian coordinates.

His original plan was unsuccessful: mathematicians would later prove that, in technical language, “one cannot find a three-dimensional associative division ring over the real numbers”. In fact, the theorem showed more: “The only associative division rings over the reals have dimensions 1, 2, or 4”. We know the one-dimensional case is the real numbers, the two-dimensional case is the complex numbers. The next case is four-dimensions. These turn out to be the quaternions.

Hamilton could not find his three-dimensional analog of the complex numbers, but he found a four-dimensional analog instead. He was able to use it for computing 3D rotations as he had hoped. Modern day graphics cards still use Hamilton’s quaternions for this purpose, this being the most efficient way of rotating meshes having millions of vertices.

What Is a Quaternion?

Just as a complex number looks like $a+bi$ where $a$ and $b$ are real numbers, a quaternion looks like $w + xi + yj + zk$ where $w, x, y, z$ are real numbers.

Multiplication of quaternions is done algebraically, except you have to remember that it is not commutative. for example, $ij=k$ but $ji=-k$ . Real numbers commute with all other quaternions, so $2i = i2$ and $3i(-1)jk2i = -6ijki$ .

Associativity and distributivity still hold, so with care, ordinary algebraic simplification operations can be done. The simplification rules are: $i^2 = j^2 = k^2 = ijk = -1$ , and all other algebraic operations on quaternions can be derived from these.

Example Computations

What is $ji$ ? We reason as follows:
$ji = ji(1)$
$= ji(-1)(-1)$
$= ji(ijk)(k^2)$
$= j(ii)j(kk)k$
$= j(i^2)j(k^2)k$
$= j(-1)j(-1)k$
$= (-1)(-1)(jj)k$
$= (1)(-1)k$
$= -k$ . Thus we have proved $ji=-k$ .
Using a similar procedure, we can prove all of the following:
$ij=k$
$ji=-k$
$jk=i$
$kj=-i$
$ki=j$
$ik=-j$ . Using these as theorems, we can now simplify arbitrary quaternion sums, differences, and products.
Let us simplify $(1+2i+3j+4k)(1-2i-3j-4k)$ :
Applying the distributive property repeatedly to $(1+2i+3j+4k)(1-2i-3j-4k)$ gives:
$1-2i-3j-4k+2i-2i(2i)-2i(3j)-2i(4k+3j-3j(2i)-3j(3j)-3j(4k)+4k-4k(2i)-4k(3j)-4k(4k)$
$=1-2i-3j-4k+2i-4i^2-6ij- 8ik+3j-6ji-9j^2-12jk+4k-8ki-12kj-16k^2$
$=1-2i-3j-4k+2i+4-6k+8j+3j+6k+9-12i+4k-8j+12i+16$ . Cancelling terms gives:
$1+4+9+16 = 30$ .
The previous example can be generalized by using variables instead of numbers for coefficients. If w, x, y, and z are real numbers, following the same procedure gives us:
$(w+xi+yj+zk)(w-xi-yj-zk) = w^2+x^2+y^2+z^2$ .

Dividing Quaternions

Using example 4 as a theorem, we can now divide quaternions. Since

$(w+xi+yj+zk)(w-xi-yj-zk) = w^2+x^2+y^2+z^2$ , and $w^2+x^2+y^2+z^2$ is a real number, we can do algebra to get the following formula: $1/(w+xi+yj+zk) = (w-xi-yj-zk)/(w^2+x^2+y^2+z^2)$ . The last division can be done because the denominator is an ordinary real number.

As an example, let us compute $1/i$ . $i$ is just $0+1i+0j+0k$ , so using the formula, $1/i = (0-1i-0j-0k)/(0^2+1^2+0^2+0^2) = -i$ . Thus, $1/i = -i$ . This of course could have been shown by noting $i^2=-1$ and dividing both sides of this equation by $-i$ . We do one more example.

$k/(i+j)$
$=k(1/(0+1i+1j+0k))$
$=k((0-1i-1j-0k)/(0^2+1^2+1^2+0^2))$
$=k(-i-j)/2$
$=-ki/2 - kj/2$
$=-j/2 + i/2$
$=i/2-j/2$ .

Quaternion Formulas

Rather than giving the algebraic procedures for manipulating quaternions by hand, we can apply those procedures to general quaternions (with variables as coefficients, like $w+xi+yj+zk$ ) to produce formulas that can be put in our code. Here are the basic formulas. Let $w+xi+yj+zk$ and $w'+x'i+y'j+z'k$ be quaternions, where the $w,w',x,x',y,y',z,z'$ are real numbers. Then (and these can be proved through algebra just as with the above calculations):

Addition: $(w+xi+yj+zk)+(w'+x'i+y'j+z'k) = (w+w') + (x+x')i + (y+y')j + (z+z')k$
Subtraction: $(w+xi+yj+zk)-(w'+x'i+y'j+z'k) = (w-w') + (x-x')i + (y-y')j + (z-z')k$
Multiplication: $(w+xi+yj+zk) (w'+x'i+y'j+z'k) = ww' - xx' - yy' - zz' + (wx'+xw' + yz' -zy')i + (wy'+yw'+zx'-xz')j + (wz'+zw'+xy'-x'y)k$
Division: $(w+xi+yj+zk) /(w'+x'i+y'j+z'k) = (w+xi+yj+zk)(w'-x'i-y'j-z'k)/({w'}^2+{x'}^2+{y'}^2+{z'}^2)$

Rotations in 3D

Suppose we want to rotate some figure $F$ which may be a mesh defined by millions of vertices. One way to specify the rotation is by specifying its rotation in three different directions, with an axis parallel to the x axis, an axis parallel to the y axis, and an axis parallel to the z axis. Order matters! Rotating 30 degrees about the x axis then 20 degrees about the y axis is not going to give the same result as rotating 20 degrees about the y axis and then 30 degrees about the x axis.

In addition to specifying the axis of rotation (a vector, so the x axis is $(1,0,0)$ , the y axis is $(0,1,0)$ , and the z axis is $(0,0,1)$ ), we need a pivot point $(a,b,c)$ to rotate around. Thus, a rotation can be specified with three quantities: a vector axis $(t,u,v)$ which has to be scaled so that it is a unit vector, that is, so that $t^2+u^2+v^2=1$ , a pivot point $(a,b,c)$ , and the angle of rotation $\theta$ . Note any nonzero vector can be scaled to a unit vector by replacing $(t,u,v)$ with $(t/l,u/l,v/l)$ where $l$ is the length of the vector, which is equal to $l = \sqrt{t^2+u^2+v^2}$ .

Using Quatnerions

We wish to convert this data to a quaternion that can be used to easily compute the rotations of all the vertices of the mesh $F$ . This quaternion is then:

$q = \cos(\theta/2) + \sin(\theta/2)(ti+uj+vk)$ .

Then, if the point of $F$ that we want to rotate has coordinates $(x,y,z)$ , we first make this into a quaternion $p = (x-a)i+(y-b)j+(z-c)k$ . Then, compute the product $qpq^{-1}$ , then express that as a quaternion $qpq^{-1}=w'+x'i+y'j+z'k$ , and the new rotated point is $(x'+a,y'+b,z'+c)$ .

Note, computing $q^{-1}$ is easy for this specific value of $q$ : $q^{-1} = \cos(\theta/2) - \sin(\theta/2)(ti+uj+vk)$ (this works out because $(t,u,v)$ is a unit vector, that is, $t^2+u^2+v^2=1$ , and also because the square of the sine of something plus the square of the cosine of the same thing is 1; after that, it’s just algebra).

So in summary, the steps are: translate the vertices of $F$ so that its pivot point is moved to the origin; Rotate about the specified axis by sending each vertex $p$ to $qpq^{-1}$ ; then translate the result back to the original pivot point.

Three Axes of Rotation

But we want to rotate, in general, about three axes, a different angle for each axis. How do we do that? Well, the information is there if you read between the lines! Do it like this:

Let $q_1$ be the value of $q$ used to perform the rotation about the $X$ axis, which is the vector $(1,0,0)$ , with the angle we want. Similarly, let $q_2$ be the quaternion for the $Y$ axis and the angle for that, and $q_3$ be for the $Z$ axis and the angle for that. Then, rotating about the $X$ axis first, then the $Y$ axis, and then the $Z$ axis, is the same as using the quaternion $q_3q_2q_1$ for rotating. Thus each point $p$ is moved to the point $q_3q_2q_1p(q_3q_2q_1)^{-1}$ which can also be written as $q_3q_2q_1pq_1^{-1}q_2^{-1}q_3^{-1}$ . Remember, order matters with quaternions: they do not generally commute! As before, we do this after translating $F$ from the pivot point to the origin, and after that, we translate the changed $F$ back to the pivot point.

Example

As an example, let us work out by hand (part of) the rotation of a cube. A former boss told me, you don’t understand a computation unless you can do it (a small case or a representative part of it) by hand.

So, suppose we have a cube whose eight vertices are at $(1,0,0)$ , $(1,0,1)$ , $(1,1,0)$ , $(1,1,1)$ , $(2,0,0)$ , $(2,0,1)$ , $(2,1,0)$ , and $(2,1,1)$ . This is a 1x1x1 cube whose faces are oriented perpendicular to the $X$ , $Y$ , and $Z$ axes and whose center is at $(1.5,0.5,0.5)$ . We shall use this center as the pivot point, and rotate the cube 45 degrees counterclockwise (when looking downward) about the $Z$ axis (or rather, the line through the pivot point that is parallel to the $Z$ axis). Here is a screenshot from Blender showing the cube:

Note that, as is somewhat standard in many graphics applications, the $X$ axis is red, the $Y$ axis is green, and the $Z$ axis is blue.

First let us compute q, the rotation quaternion. $q = \cos(\theta/2) + \sin(\theta/2)(ti+uj+vk)$ , where $\theta=45$ degrees (or $\pi/4$ radians), $t=u=0$ and $v=1$ . Thus,

$q = \cos(\pi/8) + \sin(\pi/8)k \approx 0.92387953251 + 0.38268343236k$ .

and

$q^{-1} = \cos(\pi/8) - \sin(\pi/8)k \approx 0.92387953251 - 0.38268343236k$ .

Now, for each vertex, we translate by $(-1.5,-0.5,-0.5)$ , make a quaternion $p$ out of it, compute $p'=qpq^{-1}$ , then translate $p'$ by $(1.5,0.5,0.5)$ . Here we go:

$(1,0,0) \longrightarrow (-0.5,-0.5,-0.5)$
$\longrightarrow -0.5i-0.5j-0.5k$
$\longrightarrow(0.92387953251 + 0.38268343236k) (-0.5i-0.5j-0.5k) (0.92387953251 - 0.38268343236)$
$\longrightarrow -0.707107j-0.5k$
$\longrightarrow (0,-0.707107,-0.5)$
$\longrightarrow (1.5,0.207107,0)$
so $(1,0,0)$ gets rotated to $(1.5, 0.207107,0)$ .

Similarly, we have the following rotations:

$(1,0,1)$ rotates to $(1.5, 0.207107, 1)$
$(1,1,0)$ rotates to $(0.792893 ,0.5, 0)$
$(1,1,1)$ rotates to $(0.792893 ,0.5, 1)$
$(2,0,0)$ rotates to $(2.20711 ,0.5, 0)$
$(2,0,1)$ rotates to $(2.20711 ,0.5, 1)$
$(2,1,0)$ rotates to $(1.5, 1.20711, 0)$
and $(2,1,1)$ rotates to $(1.5, 1.20711, 1)$ .

Pseudocode

Let F be a mesh with any number of vertices.

Let (a_x, a_y, a_z) be the axis of rotation, and (o_x, o_y, o_z) be the pivot.
Suppose we want to rotate by an angle t, in radians, counterclockwise
when looking down the axis (from positive to negative).

Let q = [q_w, q_x, q_y, q_z] = [cos(t/2), a_x sin(t/2), a_y sin(t/2), a_z sin(t/2)] and its
inverse q’ = [q’_w, q’_x, q’_y, q’_z] = [cos(t/2), -a_x sin(t/2), -a_y sin(t/2), -a_z sin(t/2)]

For each vertex (v_x, v_y, v_z) of the mesh F:

Let p = [p_w, p_x, p_y, p_z] = [0, v_x-o_x, v_y-o_y, v_z-o_z]

Compute p’ = [p’_w, p’_x, p’_y, p’_z] = qpq’ as follows:

p’_w = 0 (since it is equal to p_w (q_w^2 + q_x^2 + q_y^2 + q_z^2), and p_w=0)

p’_x = 2 (p_z q_w q_y + p_y q_x q_y – p_y q_w q_z + p_z q_x q_z) +
p_x (q_w^2 + q_x^2 – q_y^2 – q_z^2)

p’_y = 2 (-p_z q_w q_x + p_x q_x q_y + p_x q_w q_z + p_z q_y q_z) +
p_y (q_w^2 – q_x^2 + q_y^2 – q_z^2)

p’_z = 2 (p_y q_w q_x – p_x q_w q_y + p_x q_x q_z + p_y q_y q_z) +
p_z (q_w^2 – q_x^2 – q_y^2 + q_z^2)

Then, compute (v’_x, v’_y, v’_z) = (p’_x+o_x, p’_y+o_y, p’_z+o_z).

Replace the vertex (v_x, v_y, v_z) in F with the rotated vertex (v’_x, v’_y, v’_z)

End For

S	M	T	W	T	F	S
1	2	3	4	5	6	7
8	9	10	11	12	13	14
15	16	17	18	19	20	21
22	23	24	25	26	27	28
29	30	31