Interpolating Quaternions

There are two methods of interpolation commonly used in games, Slerp (Spherical Linear Interpolation) and Nlerp (Normalized Linear Interpolation). Slerp is torque minimal and has a constant velocity while interpolating, but it is not commutative. Nlerp on the other hand is commutative and torque minimal, but it does not have constant velocity.

Canvas support required WebGL support required

Which method to use largley depends on context. For animations where accuracy does not matter much, Nlerp is cheaper and probably the right choice. Physics simulations that need to be accurate should probably use Slerp. For more info on which one to use, check out Understanding Slerp, then Not Using It by Jon Blow.

Nlerp

The formula for linearly interpolating a quaternion is the same as the formula for linearly interpolating anything else: \(a + t(b - a)\), where \(a\) is the starting quaternion, \(b\) is the ending quaternion and \(t\) is a floating point value from \(0\) to \(1\). The quaternions \(a\) and \(b\) are assumed to be normalized. This can be easily implemented using component wise addition / multiplication and quaternion scaling. However, there are two problems with this approach.

Why would the result not take the shortest path? Both a quaternion and it's negative represent the same rotation, but they are on different hemispheres of the 4D hypersphere. If this happens, the quaternion will interpolate on the longest path, not the shortest. Testing if two quaternions are on the same hemisphere is simple, if the dot product between the two quaternions is greater than 0 they are on the same hemisphere. Moving a quaternion to the other hemisphere is equally easy, negate it.

To fix the normalization issue, just normalize the result before returning it. Putting all of this together we can write a normalized lerp function:

Quaternion Nlerp(Quaternion start, Quaternion end, float t) {
    if (Dot(start, end) < 0.0f) {
        end = Negate(end);
    }

    // start + (end - start) * t
    Quaternion result = Add(start, Scale(Sub(end, start), t));

    return Normalize(result);
}

Mixing

Quaternion addition can be blended using weights. For example, mixing two quaternions (both at 50% strength) could be done like so: Add(Scale(a, 0.5), Scale(b, 0.5)). This is especially usefull when blending animations. Nlerp can also be implemented in terms of mixing two quaternions. When interpolating this way, \(a\) should have a weight of \(1\) when \(t\) is \(0\) and similarly it shuold have a weight of \(0\) when \(t\) is \(1\). On the other hand \(b\) should have a weight of \(0\) when \(t\) is \(0\) and \(1\) when \(t\) is \(1\). This means \(a\) needs to be scaled by \(1 - t\) and \(b\) needs to be scaled by \(t\), like so:

Quaternion Nlerp_Mix(Quaternion start, Quaternion end, float t) {
    float dot = Dot(start, end);

    // start * (1 - t) + end * t
    Quaternion result = Add(Scale(start, 1.0f - t), Scale(end, dot < 0.0f? -t: t));

    return Normalize(result);
}

Power

In order to slerp a quaternion, the formula start + (end - start) * t still applies, but it must be interpreted differently. The multiplication in the formula can no longer be component wise since it should not affect the magnitude of the quaternion. If (end - start) * t can't scale the quaternion, how should the oepration be performed?

Instead of scaling, the only part of the quaternion that needs to change is it's angle. Adjusting the angle of a quaternion is commonly written as an exponent operation: start + (end - start) ^ t. To raise a quaternion to some power, break the quaternion down into it's angle axis representation, scale the angle by whatever the power is and create a new quaternion from this scaled angle axis representation.

Given a quaternion which rotates by \(\theta\) about axis \(\hat{v}\), raising it to some power \(t\) would look like this: \((cos(\frac{t\theta}{2}), \hat{v}sin(\frac{t\theta}{2}))\). The interesting thing about this operation is that only the input to \(sin\) and \(cos\) are changed. That means the result of this operation is a unit quaternion.

Quaternion Pow(Quaternion q, float power) {
    Vector3 angle = 2.0f * acos(q.w);
    float axis = Normalize(Vector3(q.x, q.y, q.z))

    float halfCos = Math.cos((power * angle) * 0.5);
    float halfSin = Math.sin((power * angle) * 0.5);

    return Quaternion(
        axis.x * halfSin, // x
        axis.y * halfSin, // y
        axis.z * halfSin, // z
        halfCos           // w
    );
}

Slerp

Given two quaternions \(a\) and \(b\) there exists a quaternion that rotates from \(a\) to \(b\). This quaternion, let's call it \(q\) can be found by multiplying the inverse of \(a\) by \(b\). Doing this inverse multiplication results in a quaternion that rotates from \(a\) to \(b\), like so \(q = b a^{-1}\). This quaternion needs to be scaled and combined with \(a\).

The interpolation formula start + (end - start) * t can be adjusted to work with quaternions by making the following changes:

The adjusted formula becomes (end * Inverse(start))^t * start. Expressing this in code is trivial. The same double cover issue that was discussed with Nlerp is present here as well. If the dot product of the start and end quaternions is negative, negate the end quaternion.

Quaternion Slerp(Quaternion a, Quaternion b, float t) {
    if (Dot(a, b) < 0) {
        b = Nagate(b);
    }

    // Only needs to be normalized if inputs are not unit length
    // ((b * Inverse(a))^t) * a
    return Mul(Pow(Mul(b, Inv(a)), t), a);
}

The math for the above function works out to be roughly this: \(Slerp(a, b, t) = (ba^{-1})^{t}a\). Reading trough engines like UE4 or Ogre, they implement slerp as: \(Slerp(a, b, t) = a\frac{sin((1 - t)\theta)}{sin(\theta)} + b\frac{sin(t\theta)}{sin(\theta)}\). This implementation is discussed in my Vector Interpolation blog, the method described there works for touples in any dimension 2 or larger.