Quaternions: How

Game developers have a love-hate relationship with quaternions. Most understand that they are somehow mathematically superior to Euler-Angles for representing rotations. Most also know that they are more compact, and easier to interpolate than rotation matrices. Other than that, quaternions remain a mystery.

I have a very deep and powerful love relationship with quaternions. For me, knowing why they work is more important than knowing how they work. However, most programmers need to see the how of quaternions before they care to ask about why. I will attempt to bestow my quaternion mojo upon you at this time, if you are willing. For now, I will just talk about HOW.

We will discuss the basic things that any game programmer actually does with quaternions, such as:

• Representing a rotation using an axis and an angle
• Rotating vectors
• Concatenating rotations
• Blending between an initial and final rotation
• Integrating quaternions using an angular velocity

How to Represent a Quaternion:

A quaternion can be represented by 4 numbers. A common struct definition might be

struct Quat
{
     float x;
     float y;
     float z;
     float w;
};

The x, y, and z components of a quaternion represent coordinates in standard 3D space. Really. Some people might try to make it more complicated than that, but they are mistaken. When using a quaternion to represent a standard vector, it is customary, but unnecesary to set the w component to zero. It is convenient to define the quaternion with the w last, so that you can easily typecast it to a 3D vector. A quaternion IS a vector, in every sense of the term vector.

You may ask: if a quaternion IS a vector, why all the fuss? That is an excellent question the answer of which plumbs the very foundations of the universe itself. (You think I'm joking, but I'm not!!) I will give you the answer when I'm ready to discuss the why of the quaternion.

When talking about quaternions it is convenient to use (scalar, vector) notation. The scalar I'm referring to is the w component, and the vector is the x, y, and z components. For instance, I might use the following notation

q = (w, V)

In order to use a quaternion to represent a rotation, you need to know the angle θ of the rotation, and the axis n around which you are rotating. The axis n is a normalized 3D vector, and the angle θ is e.g. a float. The rotation quaternion is defined as:

r = ( cos(θ/2), n sin(θ/2) )

Looking at this, you may think that it almost makes sense; you wonder why θ/2, and not just θ? This is a wonderful question, and the answer is both beautiful and profound. (You think I'm joking, but I'm not!!) But I'm not answering why yet, I'm just discussing how. One thing that is very important about this formula: the angle is in radians, NOT degrees!

To generate a quaternion with a given axis and angle, you may envision creating a function like this

Quat QuatFromAxisAngle(const Vector& axis, float angleInRadians)
{
     Quat result;
     float angle = angleInRadians / 2.0f;
     float sinAngle = sin(angle);
     Vector n = VecNormalize(axis);
     result.w = cos(angle);
     result.x = n.x * sinAngle;
     result.y = n.y * sinAngle;
     result.z = n.z * sinAngle;
     return result;
}

You now have a quaternion that represents a rotation.

How To Rotate a Vector With a Quaternion:

Quaternions aren't just a set of 4 numbers, they are an algebra. This means that there is a procedure defined to add and multiply these quantities. You can perform a vector rotation using a formula involving quaternion multiplication. But I'm not going to talk about that yet, because most of the time you aren't going to use the quaternion rotation formula.

The first rule of rotating vectors with quaternions is: don't rotate vectors with quaternions! Quaternions are great for representing rotations, but when you get to the point in your code where you are doing the actual rotation computation, it's better to convert to a matrix form. Using the quaternion multiplication rule to rotate a vector isn't too expensive, but nothing is more efficient than a matrix for transforming vectors. In fact, it is cheaper to convert the quaternion to a matrix, and use the matrix to rotate the vector, than to use the quaternion formula.

As it turns out, most transforming of vectors is done by the graphics library and not by you. However, transformations are ubiquitously represented by matrices in graphics libraries. If you want to send your rotation to the graphics library, you will need to convert it to a matrix. Therefore it is essential for you to know how to convert a quaternion to a matrix.

To convert a quaternion to a matrix use this function:

Matrix MatrixFromQuaternion(const Quat& q)
{
     Matrix result;
     //helper quantities - calculate these up front
     //to avoid redundancies
     float xSq = q.x * q.x;
     float ySq = q.y * q.y;
     float zSq = q.z * q.z;
     float wSq = q.w * q.w;
     float twoX = 2.0f * q.x;
     float twoY = 2.0f * q.y;
     float twoW = 2.0f * q.w;
     float xy = twoX * q.y;
     float xz = twoX * q.z;
     float yz = twoY * q.z;
     float wx = twoW * q.x;
     float wy = twoW * q.y;
     float wz = twoW * q.z;
     //fill in the first row
     result.m00 = wSq + xSq - ySq - zSq;
     result.m01 = xy - wz;
     result.m02 = xz + wy;
     //fill in the second row
     result.m10 = xy + wz;
     result.m11 = wSq - xSq + ySq - zSq;
     result.m12 = yz - wx;
     //fill in the third row
     result.m20 = xz - wy;
     result.m21 = yz + wx;
     result.m22 = wSq - xSq - ySq + zSq;
     return result;
}

This function does not assume that the input quaternion is normalized. A quaternion only represents a rotation if it is normalized. If it is not normalized, then there is also a uniform scale that accompanies the rotation. Normalizing a quaternion is similar to normalizing a vector. You just have to take into account that the quaternion has 4 components. To be explicit, here is a function to normalize a quaternion.

Quat QuatNormalize(const Quat& q)
{
     Quat result;
     float sq = q.x * q.x;
     sq += q.y * q.y;
     sq += q.z * q.z;
     sq += q.w * q.w;
     //detect badness
     assert(sq > 0.1f);
     float inv = 1.0f / sqrt(sq);
     result.x = q.x * inv;
     result.y = q.y * inv;
     result.z = q.z * inv;
     result.w = q.w * inv;
     return result;
}

Now, you may be wondering how I came up with this magic formula to convert quaternions to matrices. I can do it, because I know how to use the quaternion rotation formula. You may also be wondering why you would ever mess around with quaternions if you just are going to convert it to a matrix anyway. Keep reading... you will see.

How to Concatenate Rotations:

When you multiply two transformation matrices together, the result (aside from any numerical error) is also a transformation matrix. Transformation matrices can rotate, translate, scale, and skew. However, in many cases the only operation being performed by a transformation matrix is a rotation.

If you multiply two 3x3 rotation matrices together, there are 27 total terms that need to be evaluated. If you multiply 2 quaternions together, there are only 16. This number can actually be optimized a bit, but even with a naive approach you can perform about 1.5 quaternion multiplies for every matrix multiply. Thus, if you know your transformations only involve rotations, using a quaternion is a very good thing.

The order of quaternion multiplication is important, so you need to keep track. Just remember, the first rotation is on the right, and the second rotation is on the left.

q_T = q₂q₁

This multiplication is a quaternion multiplication. The way that quaternion multiplication is defined is one of the things that makes quaternions good at representing rotations.

You can express the quaternion multiplication in terms of standard vector operations, such as dot and cross products.

A_qB_q = (ab - A•B, aB + bA + A×B)

A function which multiplies two quaternions together may be defined in terms of the components as follows:

Quat QuatMultiply(const Quat& q1, const Quat& q2)
{
     Quat result;
     result.w = q1.w*q2.w - q1.x*q2.x - q1.y*q2.y - q1.z*q2.z;
     result.x = q1.w*q2.x + q1.x*q2.w + q1.y*q2.z - q1.z*q2.y;
     result.y = q1.w*q2.y + q1.y*q2.w + q1.x*q2.z - q1.z*q2.x;
     result.z = q1.w*q2.z + q1.z*q2.w + q1.x*q2.y - q1.y*q2.x;
     return result;
}

This is just an expanded version of the vector operations given previously.

How to Interpolate a Quaternion:

This is a very huge topic, but I will boil it down for you and give you the goods.

Interpolating a quaternion is useful when smoothly varying between an initial and final rotation. Interpolation is good for finding arbitrary in-between values of rotation. This is employed especially in character animation systems. It is possible to interpolate rotation matrices, but the interpolated matrix may not be a size and shape preserving matrix. Needless to say, interpolating a quaternion is a bajillion times easier than interpolating rotation matrices.

There is one interesting property of quaternions that comes into play when dealing with interpolation. If a rotation can be represented by a quaternion q, then the quaterion -q also represents the same rotation. Why is that? I'm not going to explain it right now, other than to say that it is connected to the very fabric of reality. (You think I'm joking but I'm not!) What you need to worry about is which one of these quaternions you are going to use.

To describe the difference between q and -q, consider that you turn a quarter turn to your left. Esentially this is the same as turning 3/4 turn to your right. One turn is the "short" turn and the other is the "long"one. When representing a static orientation it is irrelevant if a quaternion represents the short, or long path, because it just sits in the final position and you don't get to see the in-between values. However, when you are blending it surely does make a difference.

When blending between an initial and a final quaternion, there is some ambiguity as to if we are taking the "short" way or the "long" way. It seems like the right thing to always blend on the shortest path. Given the two input quaternions, it is possible to determine which way we are going to blend. You can check this by examining the sign of the 4D dot product of the inputs. If the sign is negative, then you know you are going to be blending the long way.

So, what do you do if you find out that you are blending the long way? Simply flip the sign on one of your input quaternions. remember q and -q represent the same rotation. Flipping the sign on one of your inputs will flip the sign of the 4D dot product.

Now that we have discussed that little tid-bit, let's move on to interpolation formulas.

There are a few different interpolation formulas, but two main ones: NLerp is a linear interpolation of the components that is followed by a normalization of the interpolated quaternion, to ensure that it represents a rotation. Slerp is a spherical interpolation which interpolates in a spherical space, rather than in the cartesian space of the coordinates. The interpolant of the slerp function moves at a constant motion, while the interpolant of the NLerp has some non-linear acceleration.

Heres the quick and dirty: Don't mess around with Slerp, even though you think it might be the more "correct" thing to do. It is too expensive, and has too many special cases that need to be considered. There are some complicated schemes that try to closely approximate the Slerp function, but it just isn't worth it. Just use NLerp. Especially for computationally strapped code.

In fact, I'm not even going to show you how to SLerp. You can consult google if you really want to know.

Here is a blending function that uses NLerp

Quat QuatBlend(const Quat& i, const Quat& f, float blend)
{
     Quat result;
     float dot = i.w*f.w + i.x*f.x + i.y*f.y + i.z*f.z;
     float blendI = 1.0f - blend;
     if(dot < 0.0f)
     {
          Quat tmpF;
          tmpF.w = -f.w;
          tmpF.x = -f.x;
          tmpF.y = -f.y;
          tmpF.z = -f.z;
          result.w = blendI*i.w + blend*tmpF.w;
          result.x = blendI*i.x + blend*tmpF.x;
          result.y = blendI*i.y + blend*tmpF.y;
          result.z = blendI*i.z + blend*tmpF.z;
     }
     else
     {
          result.w = blendI*i.w + blend*f.w;
          result.x = blendI*i.x + blend*f.x;
          result.y = blendI*i.y + blend*f.y;
          result.z = blendI*i.z + blend*f.z;
     }
     result = QuatNormalize(result);
     return result;
}

This function has a singularity when the difference between the initial and final quaternions is a 180 degree rotation. This is due to the fact that the axis of rotation for the blend becomes ambiguous. You could potentially detect for this case, and decide on using the "up" vector for the axis of the blend. Or you could break up the blend into a few steps. This singularity is something that shows up in any interpolation scheme, not just NLerp.

How to Integrate a Quaternion:

Updating the dynamical state of a rigid body is referred to as integration. If you represent the orientation of this body with a quaternion, you will need to know how to update it. This is done with the following quaternion formula.

q' = Δq q

We calculate Δq using a 3D vector ω whose magnitude represents the angular velocity, and whose direction represents the axis of the angular velocity. We also use the time step Δt over which the velocity should be applied. Δq is still a rotation quaternion, and has the same form involving sines and cosines of a half angle. We use the angular velocity and time step to construct a vector θ, whose magnitude is the half angle, and whose direction is the axis.

θ = ωΔt/2

Note: I've included the factor of 1/2, which shows up inside the trig functions of the rotation quaternion. Expressing the rotation quaternion in terms of this vector you have

Δq = ( cos(θ), (θ/|θ|) sin(θ) )

This works well, however this formula becomes numerically unstable as |θ| approaches zero. If we can detect that |θ| is small, we can safely use the Taylor series expansion of the sin and cos functions. The "low angle" version of this formula is

Δq = (1 - |θ|²/2, θ - θ|θ|²/6)

We use the first 3 terms of the Taylor series expansion, so we should ensure that the fourth term is less than machine precision before we use the "low angle" version. The fourth term of the expansion is

|θ|⁴/24 < ε

Here is a sample function for integrating a quaternion with a given angular velocity and time step

Quat QuatIntegrate(const Quat& q, const Vector& omega, float deltaT)
{
     Quat deltaQ;
     Vector theta = VecScale(omega, deltaT * 0.5f);
     float thetaMagSq = VecMagnitudeSq(theta);
     float s;
     if(thetaMagSq * thetaMagSq / 24.0f < MACHINE_SMALL_FLOAT)
     {
          deltaQ.w = 1.0f - thetaMagSq / 2.0f;
          s = 1.0f - thetaMagSq / 6.0f;
     }
     else
     {
          float thetaMag = sqrt(thetaMagSq);
          deltaQ.w = cos(thetaMag);
          s = sin(thetaMag) / thetaMag;
     }
     deltaQ.x = theta.x * s;
     deltaQ.y = theta.y * s;
     deltaQ.z = theta.z * s;
     return QuatMultiply(deltaQ, q);
}

This is basically it! You now know how to accomplish all of the main tasks that any game programmer will usually bump up against relating to quaternions. If you are brave, you can move on to my next post, which covers a lot of details concerning the WHY of quaternions.