Quaternions

Main Reference for this page: Computer Animation: Algorithms and Techniques by Rick Parent (Chapter 2.2.4, Chapter 3.3, Appendix B.3.4, Appendix B.4.9, and Appendix B.4.10)

Other references:

Quaternions

1.     when interpolating intermediate orientations (with respect to rotation), if the object (or joint) has 3-degrees-of-freedom (can rotate on the x, y, and z axis), it is much simpler to use quaternions

2.     a quaternion is a 4x1 matrix of real numbers [s  x  y  z] which can also be represented as [s  v] with s being the scalar and v a 3D vector [x y z]

3.     a quaternion is a substitute for angle and axis representation

4.     quaternions are used to easily interpolate between key-framed rotation

5.     quaternions are also used to eliminate the gimbal lock effect

- gimbal lock occurs when 2 angles of rotation line up on top (or very near to) of each other

- this would cause your object to inadvertently be rotated in certain ways, introducing an error into your rotations

- Example: check out http://www.anticz.com/eularqua.htm.  Theres a pretty picture and everything!! (as well as downloadable openGL code)

- the axis and angle info of a quaternion can be viewed in two ways:

1.     orientation relative to initial object space definition (used for key-frame animation)

2.     orientation to be applied to object space definition (used for compound transforms where many rotation are compounded into a single rotation)

### Basic Quaternion Math

[s1, v1] + [s2, v2] = [s1+ s2,   v1 + v2]

2.      Multiplication (associative but not commutative)

[s1,  v1] ∙ [s2,  v2] = [(s1s2)  (v1  v2),      (s1v2) + (s2v1) + (v1 x v2)]

3.      Multiplicative identity

[1, (0, 0, 0)]

4.      The magnitude (very similar to length; also called norm) of a quaternion q, denoted ||q||, is equal to

(s2 + x2 + y2 + z2)½

5.      The inverse of the quaternion q

q-1 = (1/||q||)2 ∙ [s,  v]                  where q = [s,   v]

6.      A unit-length quaternion can be created by dividing a given quaternion q by its magnitude: quq / (||q||)

The special unit-length quaternion [1, (0, 0, 0)] is obtained by multiplying a quaternion by its inverse:  q·q-1 = [1, (0, 0, 0)]

7.      A point v or a vector from the origin to point v can be represented as
[0 v] or [0, v]                               (Parents book uses [0, v])

8.      Given a quaternion:

q = [sv] = q0 + qV       (where qV = Xi + Yj + Zk)

its conjugate is:

q* = [s,  v] = q0  qV   (just like conjugate for complex numbers)

9.      The dot product of two quaternions [s1,  v1] and [s2,  v2] is

[s1,  v1]  [s2,  v2] = s1 s2 + v1  v2

## Representing Rotations Using Quaternions

A rotation of angle T about axis (x, y, z) can be represented as:

q = RT, (x, y, z) = [cos(T/2),  (sin(T/2) ∙ (x, y, z))]

#### Consider q:

q = R -T, -(x, y, z) = [cos(T/2),  (sin(T/2) ∙ (x, y, z))]

= [cos(T/2),  (sin(T/2) ∙ (x, y, z))]

= [cos(T/2),  (sin(T/2) ∙ (x, y, z))]

note that q = q since they represent the same orientation

# Rotating Vectors Using Quaternions

To apply a rotation operation, R, to a vector, v, using quaternion, q, we do the following:

v = R (v) = qvq-1         (based on a Counterclockwise (CCW) Rotation for a Right-Handed system)

Example: we want to rotate point v (0, 3, 0) around the z-axis by an angle of T = 60° to get point v.

1.      Find q (see Representing Rotations Using Quaternions above):

q = RT, (x, y, z)  =  [cos(T/2),  (sin(T/2) ∙ (x, y, z))]

We know that the vector corresponding to a z-axis rotation is (0, 0, 1) so:

q  =  R60, (0, 0, 1)

=  [cos(60/2),  (sin(60/2) ∙ (0, 0, 1))]

=  [0.8660, (0.5 ∙ (0, 0, 1))]

q  =  [0.8660, (0, 0, 0.5)]

2.  Now we need qs inverse, q-1:

We know that:

q-1 = (1/||q||)2 ∙ [s,  v]

and:

||q|| = (s2 + x2 + y2 + z2)½

With q  = [0.8660, (0, 0, 0.5)],

||q|| = (s2 + x2 + y2 + z2)½

= (0.86602 + 02 + 02 + 0.52)½

= (0.75 + 0.25)½

= (1)½

||q|| = 1

Now we can find q-1:

q-1 = (1/||q||)2 ∙ [s,  v]

= (1/1)2 ∙ [0.8660,  (0, 0, 0.5)]

q-1 = [0.8660, (0, 0, 0.5)]

3.  We can now rotate point v using the quaternion q.

v = qvq-1

What we should do:

A.     Perform the calculations for qv

B.     Then do calculations for (qv) ∙ q-1

The result should be:  v = [0,  (2.5981, 1.5, 0)]

Recall that for multiplication of quaternions q1 and q2:

q1q2 = [s1,  v1] ∙ [s2,  v2] = [(s1s2 )  (v1  v2) ,      (s1v2) + (s2v1) + (v1 x v2)]

Step A: Perform the calculations for qv:

Before we start, we should find (v1 x v2), where q = [s1, v1], v = v2 and s2 = 0:

(v1 x v2) =  (0, 0, 0.5) x (0, 3, 0)

=  (  (v1y v2z  v1zv2y),   (v1zv2x  v1x v2z),   (v1xv2y  v1yv2x)  )

=  (  (0 ∙  0  0.5 ∙ 3),   (0.5 ∙ 0  0 ∙  0),   (0 ∙ 3  0 ∙ 0)  )

=  (  (0  1.5),    (0  0),    (0  0)  )

(v1 x v2) =  (1.5, 0, 0)

Now compute qv:

qv  = [0.8660, (0, 0, 0.5)] ∙ [0, (0, 3, 0)]

= [(s1s2 )  (v1  v2),      (s1v2) + (s2v1) + (v1 x v2)]

= [(0.8660 ∙ 0)  ((0, 0, 0.5)  (0, 3, 0)),      (0.8660 ∙ (0, 3, 0)) + (0 ∙ (0, 0, 0.5)) + ((0, 0, 0.5) x (0, 3, 0))]

= [0  (0 + 0 + 0),      (0, 2.5981, 0) + (0, 0, 0) + ((0, 0, 0.5) x (0, 3, 0))]

= [0,  (0, 2.5981, 0) + (0, 0, 0) + ((0, 0, 0.5) x (0, 3, 0))]

= [0,  (0, 2.5981, 0) + (0, 0, 0) + ( 1.5, 0, 0)]

qv  = [0,  (1.5, 2.5981, 0)]

Step B: Perform the calculations for (qv) ∙ q-1:

Before we start, we should find (v1 x v2), where qv = [s1, v1] and q-1 = [s2, v2]

(v1 x v2) =  (1.5, 2.5981, 0) x (0, 0, 0.5)

=  (  (v1y v2z  v1zv2y),   (v1zv2x  v1x v2z),   (v1xv2y  v1yv2x)  )

=  (  (2.5981 ∙ (0.5)  0 ∙ 0),   (0 ∙ 0  (1.5) ∙ (0.5)),   ((1.5) ∙ 0  2.5981 ∙ 0)  )

=  (  (1.2991  0),   (0  0.75),   (0  0)  )

(v1 x v2) =  (1.2991, 0.75, 0)

v = R (v)

= R ([0, (0, 3, 0)])

= qvq-1

= [0,  (1.5, 2.5981, 0)] ∙ [0.8660, (0, 0, 0.5)]

= [(s1s2 )  (v1  v2),      (s1v2) + (s2v1) + (v1 x v2)]

= [(0 ∙ 0.8660)  ((1.5, 2.5981, 0)  (0, 0, 0.5)),      (0 ∙ (0, 0, 0.5)) + (0.8660 ∙ (1.5, 2.5981, 0)) + ((1.5, 2.5981, 0) x (0, 0, 0.5))]

= [0  (0 + 0 + 0),      (0, 0, 0) + (1.2990, 2.25, 0) + ((1.5, 2.5981, 0) x (0, 0, 0.5))]

= [0,  (0, 0, 0) + (1.2990, 2.25, 0) + (1.2991, 0.75, 0)]

v = [0,  (2.5981, 1.5, 0)]

Note: Rotations represented by quaternions p and q can be compounded and represented as:

Rq (Rp (v)) = q ∙ (pvp-1) ∙ q-1

= qpvqp-1

= Rqp (v)

Mathematical Background on Quaternions

See: http://mathworld.wolfram.com/Quaternion.html

A quaternion is somewhat like a complex number.  Recall that for a complex number, such as 6.3+4.5i, there is a real part (here 6.3) and an imaginary part (here 4.5i, where i denotes the square root of 1, i.e., (1)0.5).  We know that i2 = 1.

A key formula for quaternion algebra states:

i2 = j2 = k2 = ijk = 1

If we try to interpret these values with i as denoting the square root of 1, then so do j and k, but then ijk = 1 makes no sense.  Any easy interpretation to understand is that i, j, k, and 1 are 4x4 matrices called the basis matrices, as described below.

As well, a quaternion can be written as q = a + bi + cj + dk, which appears to be a real part and three imaginary parts.

A quaternion can also be treated as a 2x2 matrix of complex numbers.

Here, w and z are complex numbers, a, b, c, and d are real numbers, and w* and z* are the complex conjugates of w and z.

Quaternion basis matrices (these were not the original formulation, since quaternions were invented before matrices):

(William) Hamiltons Rules

i2 = j2 = k2 = 1

ij = -ji = k

jk = -kj = i

ki = -ik = j

More facts:

Recall that the inverse of the quaternion q

q-1 = (1/||q||)2 ∙ [s,  v]                  where q = [s,   v]

Another formula is:

q-1 = q* / (qq*)

Proof:  check using q = [s, v] = [s, (x, y, z)]

q-1 = q* / (qq*)

= [s, -v] / [s, v][s, -v] = [s, -v] /  [(s s)  (v  -v),      (s)(-v) + (s v) + (v x -v)]

= [s, -v] / [s2  (-x2 y2  z2),    (0, 0, 0) + (0, 0, 0)]

= [s, -v] / [s2 + x2 + y2 +z2, (0, 0, 0)]

= [s, -v] / [||q||2, (0, 0, 0)]                                 // by definition of ||q||

= (1/||q||2) [s, -v] / [1, (0, 0, 0)]                         // move constant to the front

= (1/||q||2) [s, -v]                                               // since [1, (0, 0, 0)] is the multiplicative identity element

For unit-length quaternions, where ||q|| = 1

q-1 = (1/||q||)2 ∙ [s,  v] = (1/1)2 ∙ [s,  v] =  [s,  v] = q*

Interpolation of Rotations Represented by Quaternions

·        linear interpolation of rotation using quaternions generates non-uniform motion around the circle of rotation

·        the magnitude (length) of q is removed from q resulting in a unit quaternion

·        a point of the unit sphere in 4D space is calculated by: q / (||q||)

·        the interpolated motion is non-uniform because, given two key-framed positions, linearly interpolated intermediate orientations (rotations) form a straight-line, equal-interval path.  When this path is mapped onto the unit sphere, the distance between the points is no longer equidistant.

·        thus, to get uniform motion, one must interpolate directly on the surface of the unit sphere (along the arc between key-frame 1 and key-frame 2)

·        first, find the shortest path between the two quaternion points (p and q)

o       we recall that q and q represent the same orientation

o       therefore, a rotation from quaternion q to quaternion p is the same as a rotation from q to p

o       one of the paths will be shorter than the other (the one with the smaller angle)

o       recall that one definition of a dot product is a  b  = |a|| ||b|| cos T, where a and b are vectors, ||a|| and ||b|| are the magnitudes of the vectors, and T is the angle between the two vectors (such that 0 <= T <= π, in radians).

o       This could be rewritten as (a  b) / (||a|| ||b||) =  cos T

o       The same relationship holds for quaternions.  For unit quaternions, it can be simplified to a  b  =  cos T.

o       To find the shortest path (either to p or p), we first find the cosine of the angle between q and p:

cos T    = q  p

= (sq ∙ sp ) + (vq   vp)

o       if cos T is positive, the path from q to p is shorter.  If cos T is less than or equal to zero, the path from q to p is shorter.

Lets see if that makes sense intuitively.  We are interested in looking at the angle on either side of q.  Consider the cosine for each of these angles.  Remember that for any angle T, the cosine is positive from 0 to π/2 = 90°.  For angles of π/2 to π, it is negative.  See the graph of the cosine function below.

For example in the circle picture shown above, T is greater than π/2, so cos T will be negative, and thus q is not closer to p.  On the other hand, the other angle, 180  T is less than π/2, so q is closer to p.  In the diagram below, the pink shaded area indicates the area around q where cos is positive.

·        Spherical linear interpolation gives intermediate points on the sphere between two points on the sphere

·        once the shortest path has been selected, the spherical linear interpolation (slerp) of the path may be calculated (given parameter u: 0 ≤ u ≤ 1)

slerp (qpu) = {((sin((1  u) ∙ T))/(sin T)) ∙ q} + {(sin(uT)) / (sin T) ∙ p}

·        For example, think of T = 90° = π/2.  We get:

= {((sin((1  u) ∙ π/2))/(sin π/2)) ∙ q} + {(sin(u ∙ π/2)) / (sin π/2) ∙ p}

When u = 0, we get:

= {((sin((1) π/2))/(sin π/2)) ∙ q} + {(sin(0 ∙ π/2)) / (sin π/2) ∙ p}

= {((sin(π/2))/(sin π/2)) ∙ q} + {(sin(0 ∙ π/2)) / (sin π/2) ∙ p}

= {1 ∙ q} + {0 ∙ p}

= q

When u = 0.5, we get:

= {((sin((1  0.5) π/2))/(sin π/2)) ∙ q} + {(sin(0.5 ∙ π/2)) / (sin π/2) ∙ p}

= {((sin((0.5) π/2))/(sin π/2)) ∙ q} + {(sin(0.5 ∙ π/2)) / (sin π/2) ∙ p}

= {((sin(π/4))/(sin π/2)) ∙ q} + {(sin(π/4)) / (sin π/2) ∙ p}

Since sin(π/2) = 1, we can reduce this to:

= {sin(π/4) ∙ q} + {sin(π/4) ∙ p}

We are half way through rotating by 90° and the influences of q and p are equal, as they should be.

·        Many graphics-oriented languages have a built-in slerp function

·

 VB Public Shared Function Slerp( _ ByValq1AsQuaternion, _ ByValq2AsQuaternion, _ ByValtAsSingle _ )AsQuaternion C# (C++) public static Quaternion Slerp( Quaternion q1, Quaternion q2, float t ); JScript public static function Slerp( q1:Quaternion, q2:Quaternion, t:float ) : Quaternion;

·        the interpolated points must now be normalized (made into unit quaternions)

·        once the points have been normalized, in order to avoid first-order (tangential) discontinuity, the points must go through cubic Bezier smoothing

o       to smooth the curve, we take two consecutive points on the curve, find their control points, fit all four points to a cubic Bezier curve, and then find the midpoint (or point corresponding to u) of the curve and average it with the initial point

o       the initial curve should fit a cubic equation

o       as a walkthrough, lets first try it with 2D points

P(u) = au3 + bu2 + cu + d

o       assume each point gets two control points (aka internal control points), one before (bn) the point and one after (an) the point.  To define control points given a curve:

§         to define the control points for point pn, add vector pn -1 to pn to point pn

§         once the vector has been added, find the average of the new point with point pn +1

§         we now have control point an

§         to calculate control point bn, draw a vector from an to point pn

§         add that vector to point pn

§         we now have both control points

o       once all the control points for all interpolated points have been defined, define a cubic Bezier curve segment with points (pn, an, bn+1, pn+1)

o       a way of finding a midpoint (or point corresponding to u) on the curve, is De Casteljaus method of geometrically constructing a cubic Bezier curve:

§         given a parameter u, lines are drawn between given points (2 points - pn and pn+1 - plus 2 control points - an and bn+1) and a new point is drawn on the line at u distance from the initial point

§         thus a point is found on the curve representing u

o       the new point on the curve (represented by a star) can now be averaged with the old point to smooth the curve

·        the newly smoothed 4D points can now be used as inbetween frames containing linear motion

·        to transfer to a 4D space, we use quaternions for all points, concatenated rotations for vector additions, and slerping to find the average orientation