This page was generated from doc/rotation/squad.ipynb. Interactive online version: .

Spherical Quadrangle Interpolation (Squad)

The Squad method was introduced by Shoemake [Sho87]. For a long time, his paper was not available online, but thanks to the nice folks at the Computer History Museum (who only suggested a completely voluntary donation), it is now available as PDF file on their website.

The main argument for using Squad over De Casteljaus’s algorithm with Slerp is computational efficiency:

Boehm [Boe82], in comparing different geometric controls for cubic polynomial segments, describes an evaluation method using “quadrangle points” which requires only 3 Lerps, half the number needed for the Bézier method adapted in Shoemake [Sho85].

—Shoemake [Sho87]

Given the start and end points \(p\) and \(q\) of a curve segment and the so-called quadrangle points \(a\) and \(b\), Shoemake provides an equation for Squad:

The interpretation of this algorithm is simple: \(p\) and \(q\) form one side of a quadrilateral, \(a\) and \(b\) the opposite side; the sides may be non-parallel and non-coplanar. The two inner Lerps find points on those sides, then the outer Lerp finds a point in between. Essentially, a simple parabola drawn on a square is subjected to an arbitrary bi-linear warp, which converts it to a cubic. Transliterated into Slerps, Boehm’s algorithm gives a spherical curve, \begin{equation*} \operatorname{Squad}(p, a, b, q; \alpha) = \operatorname{Slerp}( \operatorname{Slerp}(p, q; \alpha), \operatorname{slerp}(a, b; \alpha); 2(1-\alpha)\alpha ) \end{equation*}

—Shoemake [Sho87]

Shoemake also derives equations for the quadrangle points, which involves differentiation of Squad and assuming tangent vectors similar to uniform Euclidean Catmull–Rom splines.

Given a series of quaternions \(q_n\), use of Squad requires filling in values \(a_n\) and \(b_n\) on both sides of the interpolation points, so that each “cubic” segment is traced out by \(\operatorname{Squad}(q_n, a_n, b_{n+1}, q_{n+1}; \alpha)\) […] the values for \(a_n\) and \(b_n\) are given by \begin{equation*} a_n = b_n = q_n \exp\left( - \frac{\ln\left(q_n^{-1} q_{n+1}\right) + \ln\left(q_n^{-1} q_{n-1}\right)}{4} \right) \end{equation*}

—Shoemake [Sho87]

Note

Allegedly, the proof of continuity of tangents by Shoemake [Sho87] is flawed. Kim et al. [KKS96] and Dam et al. [DKL98] provide new proofs, in case somebody wants to look that up.

The equation for the inner quadrangle points uses relative rotations in the local frame of reference defined by \(q_i\). Since we have mainly used rotations in the global frame of reference so far, we can also rewrite this equation to the equivalent form (changing the index \(n\) to \(i\) while we are at it)

\begin{equation*} a_i = b_i = \exp\left( -\frac{\ln\left(q_{i+1} q_i^{-1}\right) + \ln\left(q_{i-1} q_i^{-1}\right)}{4} \right) \, q_i. \end{equation*}

We can try to get some intuition by looking at the Euclidean case. Euclidean quadrangle interpolation is shown in a separate notebook and we know how to calculate outgoing and incoming quadrangle points for uniform Euclidean Catmull–Rom splines:

\begin{equation*} \boldsymbol{\bar{x}}_i^{(+)} = \boldsymbol{\bar{x}}_i^{(-)} = \boldsymbol{x}_i - \frac{ (\boldsymbol{x}_{i+1} - \boldsymbol{x}_i) + (\boldsymbol{x}_{i-1} - \boldsymbol{x}_i) }{4}. \end{equation*}

With a bit of squinting, we can see that this is analogous to the quaternion equation shown above.

To show an example, we import splines.quaternion.Squad and a few helper functions from helper.py …

from splines.quaternion import Squad
from helper import angles2quat, animate_rotations, display_animation

… we define a sequence of rotations …

rotations = [
    angles2quat(0, 0, 0),
    angles2quat(90, 0, -45),
    angles2quat(-45, 45, -90),
    angles2quat(135, -35, 90),
    angles2quat(90, 0, 0),
]

… and create a Squad object:

sq = Squad(rotations)

For comparison, we use splines.quaternion.CatmullRom with the same sequence of rotations:

from splines.quaternion import CatmullRom

cr = CatmullRom(rotations, endconditions='closed')

import numpy as np

def evaluate(spline, frames=200):
    times = np.linspace(
        spline.grid[0], spline.grid[-1], frames, endpoint=False)
    return spline.evaluate(times)

ani = animate_rotations({
    'Squad': evaluate(sq),
    'Catmull–Rom-like': evaluate(cr),
})
display_animation(ani, default_mode='loop')

Once Loop Reflect

As you can see, the two splines are nearly identical, but not quite:

max(max(map(abs, q.xyzw)) for q in (evaluate(sq) - evaluate(cr)))

0.04640377605179979

Non-Uniform Parameterization

Shoemake [Sho87] uses uniform parameter intervals and doesn’t talk about the non-uniform case at all. But we can try! In the notebook about non-uniform Euclidean Catmull–Rom splines we have seen the equations for the Euclidean quadrangle points (with \(\Delta_i = t_{i+1} - t_i\)):

\begin{align*} \boldsymbol{\bar{x}}_i^{(+)} &= \boldsymbol{x}_i - \frac{\Delta_i}{2 (\Delta_{i-1} + \Delta_i)} \left( (\boldsymbol{x}_{i+1} - \boldsymbol{x}_i) + \frac{\Delta_i}{\Delta_{i-1}} (\boldsymbol{x}_{i-1} - \boldsymbol{x}_i) \right)\\ \boldsymbol{\bar{x}}_i^{(-)} &= \boldsymbol{x}_i - \frac{\Delta_{i-1}}{2 (\Delta_{i-1} + \Delta_i)} \left( \frac{\Delta_{i-1}}{\Delta_i} (\boldsymbol{x}_{i+1} - \boldsymbol{x}_i) + (\boldsymbol{x}_{i-1} - \boldsymbol{x}_i) \right) \end{align*}

This can be “translated” to unit quaternions:

\begin{align*} \bar{q}_i^{(+)} &= \exp\left( -\frac{\Delta_i}{2 (\Delta_{i-1} + \Delta_i)} \left( \ln\left(q_{i+1} q_i^{-1}\right) + \frac{\Delta_i}{\Delta_{i-1}} \ln\left(q_{i-1} q_i^{-1}\right) \right) \right)\, q_i\\ \bar{q}_i^{(-)} &= \exp\left( -\frac{\Delta_{i-1}}{2 (\Delta_{i-1} + \Delta_i)} \left( \frac{\Delta_{i-1}}{\Delta_i} \ln\left(q_{i+1} q_i^{-1}\right) + \ln\left(q_{i-1} q_i^{-1}\right) \right) \right)\, q_i \end{align*}

These two equations are implemented in splines.quaternion.Squad.

Being able to use non-uniform time values means that we can create a centripetal Squad spline:

sq2 = Squad(rotations, alpha=0.5)

cr2 = CatmullRom(rotations, alpha=0.5, endconditions='closed')

ani = animate_rotations({
    'Squad': evaluate(sq2),
    'Catmull–Rom-like': evaluate(cr2),
})
display_animation(ani, default_mode='loop')

Once Loop Reflect

The two movements are still very close.

max(max(map(abs, q.xyzw)) for q in (evaluate(sq2) - evaluate(cr2)))

0.05019343803811403

Let’s try some random non-uniform parameter values:

times = 0, 0.75, 1.6, 2, 3.5, 4

sq3 = Squad(rotations, times)

cr3 = CatmullRom(rotations, times, endconditions='closed')

ani = animate_rotations({
    'Squad': evaluate(sq3),
    'Catmull–Rom-like': evaluate(cr3),
})
display_animation(ani, default_mode='loop')

Once Loop Reflect

Now the two movements have some obvious differences.

max(max(map(abs, q.xyzw)) for q in (evaluate(sq3) - evaluate(cr3)))

0.42509139916677563

With more uneven time values, the behavior of the Squad curve becomes more and more erratic. The reason for this might be the fact that the quadrangle control points are in general much further away from the curve than the Bézier control points. To check this, let’s show the angle between adjacent control points in each segment, starting with the Bézier control points of our Catmull–Rom-like spline:

%precision 1
[[np.degrees(q1.rotation_to(q2).angle) for q1, q2 in zip(s, s[1:])]
 for s in cr3.segments]

[[17.0, 106.3, 19.9],
 [22.5, 107.3, 91.0],
 [42.8, 83.2, 44.9],
 [168.4, 209.5, 68.6],
 [22.9, 57.2, 11.3]]

An angle of 180 degree would mean a quarter of a great circle around the unit hypersphere.

Let’s now compare that to the quadrangle control points:

[[np.degrees(q1.rotation_to(q2).angle) for q1, q2 in zip(s, s[1:])]
 for s in sq3.segments]

[[67.6, 206.7, 48.6],
 [62.4, 205.0, 108.4],
 [24.0, 209.4, 17.9],
 [251.1, 259.1, 103.9],
 [11.5, 130.6, 30.1]]

The angles are clearly much larger here.

With even more extreme time values, the control quaternions might even “wrap around” the unit hypersphere, leading to completely wrong movement between the given sequence of rotations. This will at some point also happen with the CatmullRom class, but with Squad it will happen much earlier.