Before playing with 3D space and plain quaternions, we should first focus on a simpler problem: 2D rotations. To do this, I will have to introduce a few notions about complex numbers. This may look odd at first sight, but it is necessary to understand quaternions.

Complex numbers are the set of numbers described by two reals a and b and the imaginary number, i, through the notation z=a+bi. The imaginary number has the property that its square is equal to -1 (i²=-1). If you are not very comfortable with this, just think of it as the writing of a concept but don't expect i to have a deep meaning about the universe; it's just some fancy writing. Still not convinced? Well, look at the following matrices:

If you write a bit of math, you will see that J²=-U and I bet this time you will not complain about it! Well, our imaginary number is not that different. Actually, you will always be able to replace any complex number z=a+bi by the matrix system Z=aU+bJ.

Let's focus back on complex number through the usual notation.

Multiplying two complex numbers is just the development of the product taking care of the imaginary number property i²=-1:

Still, this does not tell us why complex number would have anything to do with rotations. Well, actually, it does.

A complex number has two degree of freedom (its two real number a and b). For example, 3+i and 2+2i are two different complex numbers and they do not mix. At any time, we can find back the real (a in our definition) and imaginary (b in our definition) parts of these numbers. So we may decide, for some silly but actually extremely clever reasons, to write 2D vectors as complex number. A vector v (v_x,v_y) could then be written either as v_x+v_yi or v_y+v_xi. I will keep the former notation (v_x+v_yi).

Adding and subtracting vectors in their complex number notation is the same as in the vector notation. But an interesting thing happens when we multiply two complex numbers together which usually do not have sense in the vector world. Let's study the product of v=v_x+v_yi with z=cos(θ)+sin(θ)i:

which is nothing but the vector v rotated by an angle θ around its origin. If you don't believe me, just write z using the matrix notation used before:

You will surely recognize the matrix of a 2D rotation out there! So the product of Z with V=v_xU+v_yJ is the rotation of v by an angle θ. And this is the key to understand quaternion-based rotations.

But before jumping to the more complex (sorry for the joke) quaternions, there are still a few things to introduce.

Let's focus on the inverse rotation z^-1; that is, the rotation of an angle -θ. We may write it directly as:

which is the conjugate of z (the complex conjugate of z=a+bi is z'=a-bi). You can check that zz^-1=z^-1z=1 which is the property of inverse function (that is, rotation of an angle θ followed or preceded by a rotation of a angle -θ has no overall effect).

From all of this, we may say that any number z=a+bi represents a rotation if a²+b²=1. Furthermore, the inverse rotation will be equal to the complex conjugate.

You may wonder why on earth would someone use complex numbers to write rotations since matrices are working just fine for this. Well, actually there are a few advantages:

1/ The composition of two rotations is equal to the product of two complex numbers which requires four multiplications, one subtraction and one addition while the matrix notation uses eight multiplications and four additions. It is then faster to compose rotations using complex numbers.

2/ The composition of many rotations is more stable in the complex number notation. We have seen that the norm of the complex number (a²+b²) has to be equal to unity to represent a rotation. After many compositions, the norm may, for some reasons, become slightly different from unity which results in a degenerated rotation (a rotation that also has a scaling). The same thing happens with matrices too except that complex number are much easier to regenerate than matrices. All we have to do to regenerate a pure rotation from a degenerate rotation is to divide the rotation complex number by its norm (normalization process, just like for vectors).

3/ Interpolation, which is often used in animation, is much easier to do with complex number. Just use linear interpolation between two rotation complex numbers and normalize the result to have an interpolated rotation. This is in no-way trivial to do with matrices! I will talk later on interpolation; but you can already remember that complex number performs better for that too.

Finally, we should note that converting from complex number system to a matrix system is easy because by definition:

This is often useful to integrate complex number rotation into traditional matrix-based pipelines.

So why not using the same approach for 3D rotations? The problem is not new and was first raised by Sir William Rowan Hamilton (1805-1865). He first tried to solve it by adding a second imaginary number, j (j²=-1), to increase the dimensionality to 3 and match the 3D vectors degree of freedom. He got:

with no clue to what ij or ji would represent.

He finally found an answer in 1843 by adding yet another imaginary number, k, with the rules i²=j²=k²=ijk=-1. By developing these rules, a few sub-rules emerge: ij=k, jk=i, ki=j, ji=-k, kj=-i and ik=-j. It is these hyper-complex numbers with these exact sets of rules that allows representation of 3D rotations. We call them quaternions. Again I would like to stress, and to demystify them, that the ability to perform 3D rotations has nothing to do with the deep nature of imaginary numbers but only on the multiplication rules we have imposed. Actually, we could have done the same things using a matrix notation just like in the 2D case. If you are interested in the matrix notation of quaternions, have a look on Internet, it's easy to find them.

So we now have our quaternion thing that we write q=w+xi+yj+zk for which the product of two quaternions gives:

leading to an horribly long equation that is painful to read and decipher. To make it a bit easier, people often write quaternions as the sum of a scalar (w) and a vector (v): q=w+v. That way, the previous equation simplifies to:

because by definition of vector mathematics,

Here again, it's only some fancy way of writing; there are nothing magical and adding scalar to vectors has no other meaning than representing a quaternion. Remember that mathematics is only a way to manipulate ideas, not to cast spells on numbers. So you can use any notation you like as long as you keep the underlying sense of the objects you are manipulating.

Anyway, let's go back to our quaternion beast.

One of the thing we should note with the previous equations is that they are not commutative, i.e. q₁q₂≠q₂q₁ unless their vector parts are on the same line (in which case v₁xv₂=0 and q₁q₂=q₂q₁).

Hamilton also noticed (without demonstrating) that the conjugate product of a unit quaternion q=cos(θ/2)+sin(θ/2)n with a vector v, noted qvq^*, represents the rotation of v around the axis n by an angle θ. Indeed, and considering v as a pure-imaginary quaternion (a quaternion with no scalar term), we may write:

which you may recognize as the Euler-Rodrigues formula describing the rotation of the vector v around the axis n by an angle θ.

The reason why the transformation of 3D vectors takes such a weird shape compared to the 2D case is not obvious. It has something to do with the fact that using only qv do not conserve the length of v while qvq^* does. Also, neither qv nor vq^* is a vector because their scalar part is not null. The θ/2 can be explained by the fact that we are rotating in two steps (and hence twice the given angle θ/2). This is one of the things about quaternions that are more difficult to explain. Let's assume it since we have proven it to work.

One of the beautiful properties of that conjugate product is that the composition of two quaternions p and q (of angle θ₁ and θ₂ respectively) represents a rotation of θ₂ followed by the rotation θ₁ around the q and p respective axis. Indeed:

This is because the conjugate of a quaternion product is equal to the commuted product of their conjugate:

Also, the conjugate of a unit quaternion represents its rotation inverse (q^-1=q^*):

Finally, and just like for 2D rotations, quaternions will be useless if we can't transform them to rotation matrices since all the CGI pipelines are designed to cope with these. We can do this by expanding the conjugate product:

which gives after some ordering:

that we can write as a matrix-vector product v'=Mv with the following matrix:

Because the conjugate product of a unit quaternion represents a rotation, then the matrix M also represents a 3D rotation if a²+b²+c²+w²=1.

Interpolation consists of finding a smooth path between discrete key events. This is common in animation systems were the orientation of a body is only given at intervals such as to save memory when storing the animation file on disk. There are several ways to do interpolation but we will focus here on linear interpolations only.

Let us consider first a 2D rotation scenario where we would like to go from a rotation described by z₁ to a rotation described by z₂. We may write the linear interpolation z(χ) with 0≤χ≤1:

but z(χ) no longer represents a rotation because its norm is not guaranteed to be unit. Check with z₁=1 and z₂=i at χ=0.5. The solution is to normalize z(χ) such that it now represents a valid rotation. The procedure is known as the LERP method (for linear interpolation):

The LERP operation is known to be of minimal torque because it follows the shortest circular path from z₁ to z₂. On the other hand, it does not have a constant speed along this path because of the normalization operation. Finally, it commutes following the law:

When the speed-constness becomes an issue and the commutation property is no longer required, a popular interpolation formula is the SLERP method (for spherical linear interpolation). To describe the SLERP operation, we first need to build an orthonormal basis through the vector z₁ and z₂:

where α is the angle between the vectors z₁ and z₂ and k a normalization factor. We may compute k as:

we may then express the SLERP method by creating a vector z(χ) going from z₁ to z₂ using a rotation expressed in the basis z₁/z₂':

It is however not the most known form of the SLERP formula. We may re-arrange it as:

which is the well-known Shoemake formula (1985) written with vectors instead of quaternions:

An interesting property is that:

Finally, I should mention that the SLERP interpolation is computationally heavier than the LERP one. So if you don't require the constness of speed or if the commuting property is more important to you, it's better to use LERP and not SLERP. Don't fall for SLERP just because it is popular.

I hope you know better understand quaternion-based rotations and that you will no longer look at them as black magic :)