Quaternions

Gerald Farin

January 27, 2004

1 Motivation

For animation, objects have to be moved through space. If the object’s orientation changes during such a motion, a rotation is invariably involved. A typical problem is that the object is only defined at several key positions, and the intermediate ones have to be computed. In other words, we have to be able to interpolate between rotations. This is not trivial: if R_{1} and R_{2} are two rotation matrices, then the blend

Q = (1 − t)R_{1} + tR_{2}

will, in general, not be a rotation matrix since we cannot guarantee that Q^{T}Q = I as is required for rotation matrices. A different entity, called a quaternion, may also be used to describe rotations, but in addition is capable of handling the interpolation problem.

2 Introduction

Let R be rotation matrix, such that

R^{T} R = I.

First we show that R has an eigenvalue 1, i.e., there is a unit vector u, corresponding to this eigenvalue, such that

Au = u.

Assume λ is the eigenvector corresponding to u. Then

(λu)^{T}(λu) = λ^{2}u^{T}u = u^{T}R^{T}Ru = u^{T}u = 1,

hence λ^{2} = ±1. To show that in fact +1 is an eigenvector requires some more work.

Thus in order to find the axis of a general rotation matrix, all we have to do is to find its eigenvector corresponding to the eigenvalue +1. The angle Φ by which v is rotated around u to yield v′ can be computed from

trace(R) = 1 + 2 cos Φ.

Example Let a rotation matrix be given by

R = [\begin{matrix}τ^{2} & τ^{2} − τ & τ^{2} + τ \\ τ^{2} + τ & τ^{2} & τ^{2} − τ \\ τ^{2} − τ & τ^{2} + τ & τ^{2} \end{matrix} ]

with τ = 1/\sqrt{3}. The matrix has an eigenvector u = [τ, τ, τ]^{T} corresponding to the eigenvalue 1:

[\begin{matrix}τ^{2} & τ^{2} − τ & τ^{2} + τ \\ τ^{2} + τ & τ^{2} & τ^{2} − τ \\ τ^{2} − τ & τ^{2} + τ & τ^{2} \end{matrix} ] [\begin{matrix}τ \\ τ \\ τ \end{matrix}] = [\begin{matrix}τ^{3} + τ^{3} − τ^{2} + τ^{3} + τ^{2} \\ τ^{3} + τ^{2} + τ^{3} + τ^{3} − τ^{2} \\ τ^{3} − τ^{2} + τ^{3} + τ^{2} + τ^{3} \end{matrix}] = [\begin{matrix}τ \\ τ \\ τ \end{matrix}]

We have trace(R) = 3τ^{2} = 1. Thus cos Φ = 0 and thus Φ = 90. If v′ = Rv, then v′ is in the plane which contains v and has u as its normal. To see this, note that (u^{T}v)u is the perpendicular projection of v onto u, marked by a solid circle in Figure 1. If we can show that the projection of v′ onto u is the same vector, our claim is proved. That projection is given by (u^{T}v′)u. But (noting that Ru = u):

(u^{T}v′)u = (u^{T}R^{T}Rv)u = (u^{T}v)u,

thus proving our claim.

Figure 1: The vector v is rotated into v′ around axis u.

Conversely, we may ask how to compute the rotation matrix R if we are given the axis vector u and the rotation angle Φ. This matrix turns out to be (check any linear algebra text)

[\begin{matrix}u_{1}^{2} + c(1 − u_{1}^{2}) & (1 − c)u_{1}u_{2} − su_{3} & (1 − c)u_{1}u_{3} + su_{2} \\ (1 − c)u_{1}u_{2} + su_{3} & u_{2}^{2} + c(1 − u_{2}^{2}) & (1 − c)u_{2}u_{3} − su_{1} \\ (1 − c)u_{1}u_{3} − su_{2} & (1 − c)u_{2}u_{3} + su_{1} & u_{3}^{2} + c(1 − u_{3}^{2}) \end{matrix}] ,  (1)

where c = cos Φ, s = sin Φ. For a sanity check, verify that trace(R) = 1 + 2c!

The essential information defining a rotation is given by u and Φ. This information is collected in the quaternion

û = [\begin{matrix}su \\ c \end{matrix}]  (2)

Clearly we can rebuild the matrix R from this information. But instead of working with matrices, some applications related to animation are formulated much more elegantly using quaternions.

3 Definitions

Since quaternions were originally invented as a generalization of complex numbers, we now turn to a formulation which is more general than the above. Now, a quaternion is defined as

q̂ = [\begin{matrix}x \\ w \end{matrix}]

where

x = ix + jy + kz

Note that on one hand, we write x as a vector, on the other hand as a generalized complex number. It’s not a totally bad idea to think of i, j, k as the unit direction vectors of IE^{3}, but they are better thought of as generalized complex numbers:

i^{2} = j^{2} = k^{2} = −1; ij = k, jk = i, ki = j.

We also have ij = −ji etc.

The product q̂_{1}q̂_{2} of two quaternions may again expressed as a quaternion. It is found by multiplying

(ix_{1} + jy_{1} + kz_{1} + w_{1})(ix_{2} + jy_{2} + kz_{2} + w_{2})

which gives

q̂_{1}q̂_{2} = [\begin{matrix}w_{1}x_{2} + x_{1}w_{2} + y_{1}z_{2} − z_{1}y_{2} \\ w_{1}y_{2} − x_{1}z_{2} + y_{1}w_{2} + z_{1}x_{2} \\ w_{1}z_{2} + x_{1}y_{2} − y_{1}x_{2} + z_{1}w_{2} \\ w_{1}w_{2} − x_{1}x_{2} − y_{1}y_{2} − z_{1}z_{2} \end{matrix}] .

It will be helpful to rewrite this in matrix form:

q̂_{1}q̂_{2} = [\begin{matrix}w_{1} & −z_{1} & y_{1} & x_{1} \\ z_{1} & w_{1} & −x_{1} & y_{1} \\ −y_{1} & x_{1} & w_{1} & z_{1} \\ −x_{1} & −y_{1} & −z_{1} & w_{1} \end{matrix}] [\begin{matrix}x_{2} \\ y_{2} \\ z_{2} \\ w_{2} \end{matrix}] ,

which we abbreviate as

q̂_{1}q̂_{2} = M_{1}q̂_{2}.  (3)

There is a second possibility:

q̂_{1}q̂_{2} = [\begin{matrix}w_{2} & z_{2} & −y_{2} & x_{2} \\ −z_{2} & w_{2} & x_{2} & y_{2} \\ y_{2} & −x_{2} & w_{2} & z_{2} \\ −x_{2} & −y_{2} & −z_{2} & w_{2} \end{matrix}] [\begin{matrix}x_{1} \\ y_{1} \\ z_{1} \\ w_{1} \end{matrix}]

which we abbreviate to

q̂_{1}q̂_{2} = M_{2}q̂_{1}.  (4)

This multiplication is clearly not commutative!

Yet a third possibility for the computation of q̂_{1}q̂_{2} is given by

q̂_{1}q̂_{2} = [\begin{matrix}x_{1} ∧ x_{2} + w_{2}x_{1} + w_{1}x_{2} \\ w_{1}w_{2} − x_{1}x_{2} \end{matrix}]

Other operations on quaternions are most easily understood by knowing they are equivalent to matrices. For example, addition is defined by

q̂ + p̂ = [\begin{matrix}q_{v} + p_{v} \\ q_{w} + p_{w} \end{matrix}] .

Another property is linearity:

q̂(αp̂_{1} + βp_{2}) = αq̂p̂_{1} + βq̂p̂_{2}.

Since the i, j, k parts are imaginary numbers, we may define the complex conjugate q̂^{∗} of a quaternion q̂. First, recall the standard complex conjugate of a complex number:

(a + ib)^{∗} = a − ib.

For quaternions, the same pattern holds:

q̂^{∗} = [\begin{matrix}−x \\ w \end{matrix}] .

The norm (length) of a complex number is given by |a + ib|^{2} = (a + ib)(a + ib)^{∗} = (a + ib)(a − ib) = a^{2} − aib + aib − i^{2}b = a^{2} + b^{2}. The same holds for quaternions:

∥q̂∥^{2} = q̂q̂^{∗}.

Note that q̂q̂^{∗} yields a real number because of the properties of the imaginary units in quaternions.

A unit quaternion is one with ∥q̂∥^{2} = 1 and ties in with the quaternions introduced in Section 3. It is thus of the form of (5):

q̂ = [\begin{matrix}sin Φu \\ cos Φ \end{matrix}]  (5)

We see that q̂q̂^{∗} = s^{2}(u_{1}^{2} + u_{2}^{2} + u_{3}^{2}) + c^{2} = 1 since we assumed that u is a unit vector. In what follows, we will only deal with unit quaternions.

Finally, the identity quaternion:

î = [\begin{matrix}0 \\ 1 \end{matrix}] .

It is often convenient to simply treat this quaternion as the scalar 1.

We may also define the inverse q̂^{−1} of q̂, defined by

q̂^{−1}q̂ = q̂q̂^{−1} = 1.

It is given by

q̂^{−1} = \frac{q̂^{∗}}{∥q̂∥^{2}} ,

and is verified by

\frac{q̂^{−1}q̂^{∗}}{∥q̂∥^{2}} = 1.

4 Quaternions and Matrices

Rotating a vector v around an axis u by an angle 2Φ may be expressed by v′ = Rv, where R is given by (1). In terms of quaternions, the corresponding equation is

v′ = q̂vq̂^{−1} = q̂vq̂^{∗}

This makes sense since we may write v in quaternion form with a zero real part. This operation corresponds to a rotation by 2Φ if we set w = cos Φ.

Invoking (3) and (4), we obtain

q̂vq̂^{∗} = q̂(M_{2}^{∗}v) = M_{1}M_{2}^{∗}v.

Thus M_{1}M_{2}^{∗} is the desired rotation matrix R. It is explicitly given by

R = [\begin{matrix}w^{2} + x^{2} − y^{2} − z^{2} & 2xy − 2wz & 2xz + 2wy & 0 \\ 2xy + 2wz & w^{2} − x^{2} + y^{2} − z^{2} & 2yz − 2wx & 0 \\ 2xz − 2wy & 2yz + 2wx & w^{2} − x^{2} − y^{2} + z^{2} & 0 \\ & 0 & 0 & x^{2} + y^{2} + z^{2} + w^{2} \end{matrix}]

Example A rotation around the z− axis by 90 degrees is given by the unit quaternion

q̂ = [0, 0, \frac{1}{2}\sqrt{2}, \frac{1}{2}\sqrt{2}]^{T}.

If we use the above to calculate the corresponding rotation matrix, we find

R = [\begin{matrix}0 & −1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{matrix}] .

We may concatenate quaternions as follows. Transform a vector v using quaternion q̂_{1}, then transform the result using a quaternion q̂_{2}. We have

v′ = q̂_{2}(q̂_{1}vq̂_{1}^{∗})q̂_{2}^{∗}

This may easily be shown to yield

v′ = (q̂_{2}q̂_{1})v(q̂_{2}q̂_{1})^{∗}.

Another application is this: suppose you have two unit vectors v_{1} and v_{2}, and you are interested in the rotation which takes v_{1} to v_{2}. That rotation will have the axis

u = \frac{v_{1} ∧ v_{2}}{∥v_{1} ∧ v_{2}∥};

its angle Φ will be determined by cos 2Φ = v_{1}v_{2}. Note also that sin 2Φ = ∥v_{1} ∧ v_{2}∥. The quaternion q̂ which realizes this rotation is then given by

q̂ = [\begin{matrix}sin Φu \\ cos Φ \end{matrix}] = \frac{1}{\sqrt{2}} [\begin{matrix}\sqrt{1 − v_{1}v_{2}}u \\ √ & \overline{1 + v_{1}v_{2}} \end{matrix}]  (6)

where we had to use these trig identities:

cos Φ = \frac{1}{\sqrt{2}}\sqrt{1 + cos 2Φ}, sin Φ = \frac{1}{\sqrt{2}}\sqrt{1 − sin 2Φ}.

Example: Let

v_{1} = [\begin{matrix}1 \\ 0 \\ 0 \end{matrix}] , v_{2} = [\begin{matrix}0 \\ 1 \\ 0 \end{matrix}] , u = [\begin{matrix}0 \\ 0 \\ 1 \end{matrix}] .

Thus

q̂ = \frac{1}{\sqrt{2}} [\begin{matrix}0 \\ 0 \\ 1 \\ 1 \end{matrix}]

which is the quaternion from the previous example in this section.

5 Interpolation using quaternions

Before we discuss quaternions, let us first solve a simple 3D problem. Given are two unit vectors v_{0} and v_{1}, forming an angle α. We would like to rotate v_{0} into v_{1}, creating a sequence of unit vectors v(t) which, for t = 0, starts out at v_{0} and for t = α ends at v_{1}. We write v(t) as

v(t) = f(t)v_{0} + g(t)v_{1}.  (7)

The functions f and g do not sum to unity – v(t) would not be a unit vector then. We have

∥v_{0} ∧ v_{1}∥ = sin α and ∥v_{0} ∧ v(t)∥ = sin t.

Substituting (7) into this yields

v(t) = \frac{sin(α − t)}{sin(alpha)} + \frac{sin t}{sin α}.

For an illustration, see Fig. ??.

When we interpolate between unit quaternions, the same principle applies. While v(t) above traced out a unit circle, quaternions, being 4D vectors, trace out the 4D unit sphere.

Suppose q̂_{0} and q̂_{1} are two unit quaternions, such that at time t = 0 we assume orientation q̂_{0} and at time t = 1, we assume q̂_{1}. What is a meaningful orientation at an arbitrary time t? The answer is given by spherical linear interpolation or slerp:

q̂(t) = \frac{sin((1 − t)Φ)}{sin Φ}q̂_{0} + \frac{sin(tΦ)}{sin Φ}q̂_{1}.

where cos Φ = q̂_{0} · q̂_{1} (dot product, not quaternion multiplication).

Figure 2: Rotation interpolation in 3D. 8