Computer Graphics · Animation

Slerp, Squad, and the Cost of Interpolating Rotations

Why component-wise quaternion interpolation fails, what slerp computes and what it costs, and how squad and exponential maps extend the same idea to a whole curve.

A rotation channel is not four independent scalars, but it is stored as four, and the storage shape invites the mistake of treating it like a position channel with one extra component. Component-wise interpolation of unit quaternions produces a path that leaves the unit sphere and a speed that varies with the endpoints; for rotations far apart the error is not subtle. The correction is to interpolate along the sphere, and the cost of doing so correctly is a trigonometric evaluation per sample that most of the surrounding system is designed to avoid paying per vertex.

Why the four components cannot be treated independently

A unit quaternion encodes an orientation as a rotation of angle $\theta$ about a unit axis $\hat{n}$: $q = (\cos(\theta/2), \sin(\theta/2)\hat{n})$, scalar part first. The set of unit quaternions is the 3-sphere in four dimensions, and each orientation corresponds to two points on it, $q$ and $-q$, since negating both halves of the half-angle representation produces the same rotation. That double cover is the source of the first practical rule: before interpolating from $q_0$ to $q_1$, check the sign of their dot product, and negate one if it is negative. Two keys that describe nearly the same orientation can be stored with opposite signs by different exporters, and without the flip the interpolation travels the long way around, roughly $2\pi - \theta$ instead of $\theta$. The visible result is a limb that spins most of a full turn between two nearly identical poses.

The linear interpolation $q_0 + u(q_1 - q_0)$ follows a chord through the interior of the sphere. Normalizing it projects back onto the sphere, which turns the chord into a path along a great circle, and this normalized linear interpolation is a decent approximation: the maximum angular error is small when the endpoints are close, and it grows as they separate. For a segment spanning 90 degrees, the midpoint of the chord sits noticeably inside the sphere, and after normalization the sampled angle at the midpoint lags the true half-angle. The path is a great circle either way, so the geometry is right and only the timing is wrong; the rotation arrives at each intermediate orientation at a different time than uniform angular velocity would dictate.

Uniform angular velocity is the property callers actually require. If a joint has to sweep between two orientations over a fixed interval, and holds the joint under a constant-speed constraint, the parameterization must advance the angle linearly in time. Chord interpolation does not: its angular rate is lowest at the segment midpoint. For a hard cut between two very different poses this shows as a rotation that hesitates halfway and then accelerates, which reads as a mechanical snap rather than a turn.

What slerp actually computes

Spherical linear interpolation takes the angle between the two unit quaternions, $\omega = \arccos(q_0 \cdot q_1)$, and advances along the great circle by a fraction of that angle:

$$\mathrm{slerp}(q_0, q_1, u) = \frac{\sin((1-u)\omega)}{\sin \omega}, q_0 + \frac{\sin(u\omega)}{\sin \omega}, q_1.$$

The two coefficients are not arbitrary weights. They are the barycentric coordinates of the point on the great circle, derived from the requirement that the result be a unit quaternion whose angle from $q_0$ is $u\omega$. Dividing by $\sin\omega$ normalizes them so the sum of the sines reconstructs a unit vector, and that is why the expression needs no explicit renormalization. The angle advanced per unit of $u$ is constant, $\omega$ radians across the segment, so the angular velocity is $\omega$ divided by the segment duration.

/* Unit-quaternion slerp with shortest-arc selection.
   Falls back to normalized lerp when the angle is small, where
   sin(w) approaches w and the ratio is numerically unstable. */
static void slerp(const float q0[4], float q1[4], float u, float out[4]) {
    float d = q0[0]*q1[0] + q0[1]*q1[1] + q0[2]*q1[2] + q0[3]*q1[3];
    float s = 1.0f;
    if (d < 0.0f) { d = -d; s = -1.0f; }        /* take the short arc */
    if (d > 0.9995f) {                           /* angle < ~3.6 degrees */
        for (int i = 0; i < 4; i++)
            out[i] = q0[i] + u * (s*q1[i] - q0[i]);
        float n = sqrtf(out[0]*out[0] + out[1]*out[1] + out[2]*out[2] + out[3]*out[3]);
        for (int i = 0; i < 4; i++) out[i] /= n;
        return;
    }
    float w = acosf(d);
    float a = sinf((1.0f - u) * w) / sinf(w);
    float b = sinf(u * w) / sinf(w);
    for (int i = 0; i < 4; i++) out[i] = a * q0[i] + b * (s * q1[i]);
}

The threshold in that routine is a numerical necessity, not a stylistic choice. As $\omega$ goes to zero, $\sin\omega$ goes to zero and both coefficients become ratios of vanishing quantities; the limit is the normalized lerp, which is why implementations switch to it below a dot product near one. The exact threshold is arbitrary; what matters is that the branch exists and that its output is renormalized, because the lerp path accumulates error from the subtraction.

The cost is one inverse cosine and two sines per evaluated sample. Inverse cosine of a value near one loses precision in single precision, which is a second reason for the branch. In a runtime that evaluates one rotation per joint per frame, with a skeleton of sixty joints, the trigonometric cost is payable. In a runtime that must evaluate a rotation per vertex, or per particle, the same cost is not, and the usual architecture is to evaluate the sparse rotation keys on the CPU once per frame and push the resulting orientation into the skinning palette as a matrix. The distinction is the number of consumers: slerp is a per-pose operation, not a per-vertex one.

Making a curve out of segments

Slerp gives a smooth constant-speed arc between two keys and nothing more. At a key, the incoming segment arrives with angular velocity $\omega_i / \Delta t_i$ along one plane, and the outgoing segment leaves with a different rate in a different plane. The value is continuous, the angular velocity is not, and a viewer sees a change of rotational speed exactly at the key. This is the rotational analogue of the $C^1$ break in cubic Hermite segments, with an extra wrinkle: the interpolation happens on a curved manifold, so the smoothing must too.

Squad, spherical quadrangle interpolation, is the standard fix. It treats three consecutive rotations and builds a control point that plays the role of a tangent quaternion:

$$s_i = q_i \exp!\left(-\frac{\log(q_i^{-1} q_{i+1}) + \log(q_i^{-1} q_{i-1})}{4}\right),$$

and then interpolates from $q_i$ toward $q_{i+1}$ along a spherical Bezier with control points $s_i$ and $s_{i+1}$, evaluated as two nested slerps with a shaped parameter. The exponential and logarithm are the quaternion versions of the tangent map: $\log$ takes a unit quaternion to a pure quaternion encoding half the rotation vector, and $\exp$ inverts it. Averaging the two logarithms gives a tangent that points along the bisector of the incoming and outgoing directions, and dividing by four matches the scale of a cubic Hermite’s control points. The result is $C^1$ continuous in angular velocity across keys.

Two costs come with squad. The control points depend on the neighbouring keys, so evaluating one segment requires three keys instead of two, and the control points must be recomputed whenever a key is edited. And the exponential and logarithm each require an inverse cosine plus a sine, so a squad evaluation triples the trigonometric work of slerp. Exporters that support cubic rotation curves therefore usually store the control points explicitly, which folds the cost into the file and keeps the runtime to nested slerps.

Exponential maps and the log of a rotation

The logarithm of a rotation is the axis-angle vector: direction is the axis, length is the angle in radians. Interpolating in that space means interpolating rotation vectors and exponentiating the result, which is cheap when the angle is small because the exponential series truncates well. Its attraction is that angular velocity becomes a plain vector, so a curve can be differentiated by the same rules as a position channel and a physics integrator can integrate it directly.

The attraction is bounded by a real limitation. The map from rotation vectors to rotations is many-to-one: every rotation vector of length $\theta + 2\pi k$ describes the same rotation, so the space is not a single chart but an infinite stack of balls of radius $\pi$. Two orientations whose rotation vectors differ by a full turn are the same pose but far apart in the chart, and interpolating between them spins through a full revolution. The quaternion form hides this by construction, because the double cover collapses the ambiguity to a sign, and the sign flip handles it. Systems that interpolate rotation vectors must therefore wrap angle differences into $(-\pi, \pi]$ before blending, and the wrap creates a discontinuity of its own at the boundary, at exactly the angle where the shortest path flips direction.

At angles near $\pi$ both representations degrade. For quaternions, the scalar component approaches zero, the dot product of nearby keys approaches zero, and slerp’s shortest-arc test becomes decisive on rounding noise: two keys that straddle half a turn can flip sign between samples and produce a visible jump. For rotation vectors, the same neighbourhood is the edge of the chart. The mitigation in both cases is the same, and it is a data problem rather than an arithmetic one: ensure the authored keys have enough intermediate density near half-turn rotations that no single segment spans a large angle.

Where the spherical machinery is skipped

Not every rotation path needs constant angular velocity. A runtime that blends a pose toward a target over a few frames, or interpolates between two clips, is producing a weighted average rather than a curve, and normalized component-wise blending is what the blend equations are written in. It is faster, and the timing error over a handful of frames is below the threshold of visibility. The mistake is applying that reasoning to a long segment; a channel holding a value for thirty frames and then jumping to another uses the segment as a real curve, and the chord-timing artifact becomes a mid-turn hesitation that no amount of downstream smoothing repairs.

Caching is the other place the spherical path is dropped. A runtime that queries the same curve repeatedly at the same parameter can evaluate the trig once and reuse the two coefficients. Coefficient caching per segment, recomputed only when the key pair changes, keeps slerp’s cost out of the inner loop while preserving uniform angular velocity. That is the same trade as precomputing squad control points: move the transcendental work to the point where the data changes, so the sampling path stays arithmetic.

The palette this produces is consumed downstream by joint hierarchies and by blend layers, where a quaternion’s sign convention stops being a curiosity. An additive delta computed against a base pose whose sign differs from the delta’s reference frame composes as a large rotation rather than a small correction, and the mesh deforms through the antipodal pose on its way to the intended one.

References

  1. Animating Rotation with Quaternion Curves — Shoemake, SIGGRAPH 1985
  2. Quaternions and Rotation Sequences — Kuipers, Princeton University Press, 1999
  3. Rotation Representations and Their Conversions — Geometric Tools documentation