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

Spherical Quadrangle Interpolation (Squad)

Supposedly, Shoemake [Sho87] introduced the Squad method, but sadly, his paper doesn’t seem to be available anywhere online. We have to rely on third-party accounts, for example by Watt and Watt [WW92], who state in section 15.3.8, “Parametrization of Orientation”:

Boehm [Boe82] shows how, given a Bézier curve segment \((b_0, b_1, b_2, b_3)\) one can derive the quadrangle points \((b_0, S_1, S_2, b_3)\) […]

The mathematical significance of this construction is that it shows how to construct a cubic as a series of three linear interpolations of the quadrangle points. Shoemake [Sho87] takes this construction onto the surface of the four-dimensional hypersphere by constructing a spherical curve, using three spherical linear interpolations of a quadrangle of unit quaternions. This he defines as \(\operatorname{squad}()\), where: \begin{equation*} \operatorname{squad}(b_0, S_1, S_2, b_3, u) = \operatorname{slerp}( \operatorname{slerp}(b_0, b_3, u), \operatorname{slerp}(S_1, S_2, u), 2u(1-u) ) \end{equation*}

—Watt and Watt [WW92], p. 366

Note

The original text of the above quote uses \(ut\) instead of \(u\) in the first three instances, which is most likely a typo.

Given a series of quaternion keys one can construct a cubic segment across keys \(q_i\) and \(q_{i+1}\) by constructing a quadrangle of quaternions \((q_i, a_i, b_{i+1}, q_{i+1})\) where \(a_i\), \(b_{i+1}\) have to be determined. These inner quadrangle points are chosen in such a way to ensure that continuity of tangents across adjacent cubic segments is guaranteed. The derivation for the inner quadrangle points is difficult, involving as it does the calculus and exponentiation of quaternions and we will just quote the results, referring the interested reader to Shoemake [Sho87]: \begin{equation*} a_i = b_i = q_i \exp\left( - \frac{\ln\left(q_i^{-1} q_{i+1}\right) + \ln\left(q_i^{-1} q_{i-1}\right)}{4} \right) \end{equation*}

—Watt and Watt [WW92], p. 366

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

\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*}

Even though the quote above claimed that “the derivation for the inner quadrangle points is difficult”, we can still 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

Presumably, Shoemake [Sho87] uses uniform parameter intervals and doesn’t talk about the non-uniform case. At least Watt and Watt [WW92] don’t. 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

Now we see more differences, but the two are still quite close.

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

0.06068087545750728

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

times = 0, 0.75, 1.2, 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

With more uneven time values, the behavior of the Squad curve becomes more and more erratic.