F*cking Quaternions: How Do They Work?

Quaternions are one of those interesting beasts in that most people (and I'm certainly included) don't know what the heck they are, but manage to use them just fine anyway.

Now you can point me to the Wikipedia page and show me all the lovely formulas and, yeah, three (differently) imaginary numbers and a real, with the crazy \(i^2=j^2=k^2=ijk=-1\) multiplication rules and all, and I can memorize those formulas and apply them, but really I've still got no idea what the damned thing really is...

But that's perfectly alright.

They Just Do! Deal With It

This isn't going to be an explanation of exactly what a quaternion is, and it's especially not an explanation of why it can be used to represent a rotation in 3D. There's no way I'm going to do the underlying mathematics (such as I've studied) justice on this page. This is a collection of some useful operations that allow us to use quaternions to encode rotations for the purpose of making games.

Here's (most of) what you need to know to use quaternions in game development.

First off, we're using them to represent rotations in 3D. We use them for this purpose because they're extremely flexible. Unlike Euler angles, quaternions don't suffer from ambiguous configurations (gimbal lock), at the cost of only one extra float to store (and there are techniques to compress it away in many cases). They're easy to interpolate. They're also a little bit faster for some things: multiplying two quaternions takes 16 multiplications and 12 additions (or subtractions), whereas multiplying a pair of $3\times3$ matrices takes 27 multiplications and 18 additions. And they can be freely converted to matrices for when we need to combine them with other transformations.

Now, because we're only dealing with rotations, we can restrict ourselves to a small subset of quaternions and the operations on them. So first off, we only care about unit quaternions, because those are the ones which represent rotations in \(R^3\).

On to the Code!

First things first, we need some storage space. Quaternions have four components, 3 imaginary numbers and one real. Some math libraries name the imaginary parts x, y, and z after the vector notation used in some math texts, but I'll be calling them i, j, and k, after the imaginary coefficients they represent. The real component is just w, which is fairly standard. Now some math libraries store the real part before the imaginary components, but most which I've used store the real component last, which is why I order my fields as you'll shortly see them. (Remember: it's useful to have math data types laid out the same as those of any middleware you'll use - pointer casting tricks are awesome.)

struct quat
{
    float i, j, k, w;
};

There. Now we need to actually work with this thing...

Initialization

Unit Quaternions (as we're using them here) represent a rotation of some size \(\theta\) about a line which is parallel to some vector \(\vec{a}\) and which passes through the origin. We should therefore be able to construct a quat given those values. And indeed we can:

quat quat_from_aa( const vec3 &a, angle theta )
{
    float hr = theta.rad_val * 0.5F;
    float st = sin( hr );

    quat ret =
    {
        a.x * st,
        a.y * st,
        a.z * st,

        cos( hr )
    };

    return ret;
}

And that's it. One important note: \(\vec{a}\) must be a unit vector.

This process is also easily reversed:

void quat_to_aa( vec3 &a, angle &theta, const quat &q )
{
    theta = rads( 2.0F * acosf( q.w ) );

    if( q.w >= 1.0F - 1e-5F )
    {
        //very close to zero angle, throw back any vector

        a = vec3c( 0, 0, 1 );
    }
    else
    {
        float ist = 1.0F / sqrtf( 1.0F - q.w * q.w );

        a.x = q.i * ist;
        a.y = q.j * ist;
        a.z = q.k * ist;
    }
}

The computation of ist is just one of the trig identities in place of having to compute the \(\arcsin\) (which is usually slower than the square root).

Combining Rotations

Quaternions can be combined in the same way as matrices, by using the appropriate multiplication operator. The Wikipedia article I linked above has the mathematical description of that operator. This is what it looks like as code:

quat quat::operator * ( const quat &q ) const
{
    quat ret;

    ret.i = i * q.w + w * q.i + j * q.k - k * q.j;
    ret.j = j * q.w + w * q.j + k * q.i - i * q.k;
    ret.k = k * q.w + w * q.k + i * q.j - j * q.i;
    ret.w = w * q.w - i * q.i - j * q.j - k * q.k;

    return ret;
}

That's written as an overloaded operator which is a member of the quat type, so this is on the left hand side and q is on the right hand side of the operator (order matters, as with matrices, and as you'd expect).

Interpolating Quaternions

This is my favorite quaternion-related operation, but only because I'm easily amused and it's called slerp. That stands for Spherical Linear Interpolation, and it's especially useful in animation. It's designed to interpolate the rotation such that it varies smoothly from the start to the end orientation along the shortest possible path.

It's also kind of a beast. Here goes...

quat slerp( const quat &a, const quat &b, float t )
{
    float cosTheta;
    float ta, tb;
    bool neg, norm;

    cosTheta = a.i * b.i + a.j * b.j + a.k * b.k + a.w * b.w;

    if( cosTheta < 0.0F )
    {
        //need to negate the result, or we'll go the long way around the sphere
        cosTheta = -cosTheta;
        neg = true;
    }
    else
    {
        neg = false;
    }

    if( cosTheta >= 1.0F - 1e-5F )
    {
        //quats are very close, just lerp to avoid the division by zero
        ta = 1.0F - t;
        tb = t;

        norm = false;
    }
    else if( cosTheta >= 0.2F )
    {
        //quats are somewhat close - normalized lerp will do
        ta = 1.0F - t;
        tb = t;

        norm = true;
    }
    else
    {
        //do full slerp - quickly
        float sinThetaSqr = 1.0F - (cosTheta * cosTheta);
        float sinThetaInv = rsqrt( sinThetaSqr );
        float theta = atan2_pos( sinThetaSqr * sinThetaInv, cosTheta );

        //we're good to SLERP
        ta = sin_0_HalfPi( (1.0F - t) * theta ) * sinThetaInv;
        tb = sin_0_HalfPi( t * theta ) * sinThetaInv;

        norm = false;
    }

    if( neg )
        tb = -tb;

    quat ret;

    ret.i = a.i * ta + b.i * tb;
    ret.j = a.j * ta + b.j * tb;
    ret.k = a.k * ta + b.k * tb;
    ret.w = a.w * ta + b.w * tb;

    if( norm )
    {
        float l = ret.i * ret.i + ret.j * ret.j + ret.k * ret.k + ret.w * ret.w;
        l = rsqrt( l );

        ret.i *= l;
        ret.j *= l;
        ret.k *= l;
        ret.w *= l;
    }

    return ret;
}

Now, some notes are in order. This function is slow, so we're doing a few things to speed it up at the cost of accuracy. What does this mean? It means that you're going to want to rewrite it for full accuracy in certain cases (animations where you've got ridiculously long bones, for instance). Second, it means that some of its central operations can be optimized:

• rsqrt stands for reciprocal square root. It (roughly) computes \(\sqrt{x}^{-1}\).
• atan2_pos computes what the standard atan2 computes, only it's optimized as it only needs to work in the first quadrant (that is, both of its parameters will always be positive).
• sin_0_HalfPi computes the \(\sin\) of its parameter (which is in radians). Again, it's optimized as it only has to handle input angles in the range \([0,2\pi]\).

If you don't care to write fast/approximate trig and reciprocal square root functions, then the following substitutions will work just fine:

#define rsqrt( x )            (1.0F / sqrtf( x ))
#define atan2_pos( y, x )     (atan2f( y, x ))
#define sin_0_HalfPi( a )     (sinf( a ))

Rotating a Single Point

If you're going to rotate a large number of points, then it's more efficient to convert the quaternion into a matrix and use the matrix. But I'm working my way up to that... Anyway, if you take the standard "rotate a point using a quaternion" formula, you get a simple expression which rapidly expands out into a giant mess. Here's what you're left with once you simplify that mess:

vec3 quat::rotate_point( const vec3 &p ) const
{
    vec3 ret;

    float i2 = i + i;
    float j2 = j + j;
    float k2 = k + k;

    float wi2 = w * i2;
    float wj2 = w * j2;
    float wk2 = w * k2;
    float ii2 = i * i2;
    float ij2 = i * j2;
    float ik2 = i * k2;
    float jj2 = j * j2;
    float jk2 = j * k2;
    float kk2 = k * k2;

    ret.x = p.x * (1 - jj2 - kk2) + p.y * (ij2 - wk2) + p.z * (ik2 + wj2);
    ret.y = p.x * (ij2 + wk2) + p.y * (1 - ii2 - kk2) + p.z * (jk2 - wi2);
    ret.z = p.x * (ik2 - wj2) + p.y * (jk2 + wi2) + p.z * (1 - ii2 - jj2);

    return ret;
}

Getting a Matrix

Notice something? This is exactly the same as the single-point formula, except that instead of multiplying by the point's coordinate values, we store the coefficients in a matrix. If you multiplied a single point by that matrix and inlined and simplified everything, you'd get the above function.

m3x3 m3x3c( const quat &q )
{
    float i2 = q.i + q.i;
    float j2 = q.j + q.j;
    float k2 = q.k + q.k;

    float wi2 = q.w * i2;
    float wj2 = q.w * j2;
    float wk2 = q.w * k2;
    float ii2 = q.i * i2;
    float ij2 = q.i * j2;
    float ik2 = q.i * k2;
    float jj2 = q.j * j2;
    float jk2 = q.j * k2;
    float kk2 = q.k * k2;

    m3x3 r;

    r.c[0][0] = 1 - jj2 - kk2;  r.c[1][0] = ij2 - wk2;      r.c[2][0] = ik2 + wj2;
    r.c[0][1] = ij2 + wk2;      r.c[1][1] = 1 - ii2 - kk2;  r.c[2][1] = jk2 - wi2;
    r.c[0][2] = ik2 - wj2;      r.c[1][2] = jk2 + wi2;      r.c[2][2] = 1 - ii2 - jj2;

    return r;
}

Some notes: m3x3 is a $3\times3$ matrix, the member c stands for columns, so r.c[0][1] reads "row one of column zero of the matrix r". Expanding the result to a larger matrix is as easy as padding it out with the requisite row and column of zeroes ending in a one.

Ungetting a Matrix

This is the above, only in reverse.

quat quatc( const m3x3 &a )
{
    quat ret;

    ret.w = sqrtf( max( 0.0F, 1 + a.c[0][0] + a.c[1][1] + a.c[2][2] ) ) * 0.5F; 
    ret.i = sqrtf( max( 0.0F, 1 + a.c[0][0] - a.c[1][1] - a.c[2][2] ) ) * 0.5F; 
    ret.j = sqrtf( max( 0.0F, 1 - a.c[0][0] + a.c[1][1] - a.c[2][2] ) ) * 0.5F; 
    ret.k = sqrtf( max( 0.0F, 1 - a.c[0][0] - a.c[1][1] + a.c[2][2] ) ) * 0.5F; 

    ret.i = _copysignf( ret.i, a.c[1][2] - a.c[2][1] );
    ret.j = _copysignf( ret.j, a.c[2][0] - a.c[0][2] );
    ret.k = _copysignf( ret.k, a.c[0][1] - a.c[1][0] );

    return ret;
}
 
//Returns a float with the absolute value of f and the sign of s.
inline float _copysignf( float f, float s )
{
    unsigned int uf = *(unsigned int*)&f;
    unsigned int us = *(unsigned int*)&s;

    uf = (uf & 0x7FFFFFFF) | (us & 0x80000000);

    return *(float*)&uf;
}

Some notes on this function: the matrix passed in must represent a rotation and only a rotation. That is, its columns must be mutually orthonormal.