This page was generated from doc/rotation/de-casteljau.ipynb. Interactive online version: .

De Casteljau’s Algorithm

In [Sho85], which famously introduces quaternions to the field of computer graphics, Shoemake suggests to apply a variant of the De Casteljau’s Algorithm to a quaternion control polygon, using Slerp instead of linear interpolations.

def slerp(one, two, t):
    return (two * one.inverse())**t * one

import numpy as np

helper.py

from helper import angles2quat, animate_rotations, display_animation

“Cubic”

[Sho85] only talks about the “cubic” case, consisting of three nested applications of Slerp.

The resulting curve is of course not simply a polynomial of degree 3, but something quite a bit more involved. Therefore, we use the term “cubic” in quotes.

Shoemake doesn’t talk about the “degree” of the curves at all, they are only called “spherical Bézier curves”.

def de_casteljau(q0, q1, q2, q3, t):
    slerp_0_1 = slerp(q0, q1, t)
    slerp_1_2 = slerp(q1, q2, t)
    slerp_2_3 = slerp(q2, q3, t)
    return slerp(
        slerp(slerp_0_1, slerp_1_2, t),
        slerp(slerp_1_2, slerp_2_3, t),
        t,
    )

q0 = angles2quat(45, 0, 0)
q1 = angles2quat(0, 0, 0)
q2 = angles2quat(-90, 90, -90)
q3 = angles2quat(-90, 0, 0)

times = np.linspace(0, 1, 100)

ani = animate_rotations(
    [de_casteljau(q0, q1, q2, q3, t) for t in times],
    figsize=(3, 2),
)

display_animation(ani, default_mode='once')

Once Loop Reflect

Arbitrary “Degree”

splines.quaternion.DeCasteljau class

from splines.quaternion import DeCasteljau

s = DeCasteljau([
    [
        angles2quat(0, 0, 0),
        angles2quat(90, 0, 0),
    ],
    [
        angles2quat(90, 0, 0),
        angles2quat(0, 0, 0),
        angles2quat(0, 90, 0),
    ],
    [
        angles2quat(0, 90, 0),
        angles2quat(0, 0, 0),
        angles2quat(-90, 0, 0),
        angles2quat(-90, 90, 0),
    ],
], grid=[0, 1, 3, 6])

times = np.linspace(s.grid[0], s.grid[-1], 100)

ani = animate_rotations(s.evaluate(times), figsize=(3, 2))

display_animation(ani, default_mode='once')

Once Loop Reflect

Constant Angular Speed

Is there a way to construct a curve parameterized by arc length? This would be very useful.

—[Sho85], section 6: “Questions”

from splines import ConstantSpeedAdapter

s1 = DeCasteljau([[
    angles2quat(90, 0, 0),
    angles2quat(0, -45, 90),
    angles2quat(0, 0, 0),
    angles2quat(180, 0, 180),
]])

s2 = ConstantSpeedAdapter(s1)

ani = animate_rotations({
    'non-constant speed': s1.evaluate(np.linspace(s1.grid[0], s1.grid[-1], 100)),
    'constant speed': s2.evaluate(np.linspace(s2.grid[0], s2.grid[-1], 100)),
}, figsize=(5, 2))

display_animation(ani, default_mode='reflect')

Once Loop Reflect

Joining Curves

In section 4.2, [Sho85] provides two function definitions:

\begin{align*} \operatorname{Double}(p, q) &= 2 (p \cdot q) q - p\\ \operatorname{Bisect}(p, q) &= \frac{p + q}{\|p + q\|} \end{align*}

def double(p, q):
    return 2 * p.dot(q) * q - p

def bisect(p, q):
    return (p + q).normalize()

Given three successive key quaternions \(q_{n-1}\), \(q_n\) and \(q_{n+1}\), these functions are used to compute control quaternions \(b_n\) (controlling the incoming tangent of \(q_n\)) and \(a_n\) (controlling the outgoing tangent of \(q_n\)):

\begin{align*} a_n &= \operatorname{Bisect}(\operatorname{Double}(q_{n-1}, q_n), q_{n+1})\\ b_n &= \operatorname{Double}(a_n, q_n) \end{align*}

It is unclear where these equations come from, we only get a little hint:

For the numerically knowledgeable, this construction approximates the derivative at points of a sampled function by averaging the central differences of the sample sequence.

—[Sho85], footnote on page 249

def shoemake_control_quaternions(q_1, q0, q1):
    """Shoemake's control quaternions (part 1).

    Given three key quaternions, return the control quaternions
    preceding and following the middle one.

    """
    a = bisect(double(q_1, q0), q1)
    b = double(a, q0).normalize()
    return b, a

Normalization of \(b_n\) is not explicitly mentioned in the paper, but even though the results have a length very close to 1.0, we still have to call normalize() to turn the Quaternion result into a UnitQuaternion.

But wait, we are not finished yet!

In a later section, the last missing piece of the puzzle is unveiled:

A simple check proves the curve touches \(q_n\) and \(q_{n+1}\) at its ends. A rather challenging differentiation shows it is tangent there to the segments determined by \(a_n\) and \(b_{n+1}\). However, as with Bézier’s original curve, the magnitude of the tangent is three times that of the segment itself. That is, we are spinning three times faster than spherical interpolation along the arc. Fortunately we can correct the speed by merely truncating the end segments to one third their original length, so that \(a_n\) is closer to \(q_n\) and \(b_{n+1}\) closer to \(q_{n+1}\).

—[Sho85], section 4.4: “Tangents revisited”

No further hints are given about how to actually implement this, but it should work something like this:

def shoemake_corrected_control_quaternions(q_1, q0, q1):
    """Shoemake_s control quaternions, corrected."""
    b, a = shoemake_control_quaternions(q_1, q0, q1)
    return q0.rotation_to(b)**(1/3) * q0, q0.rotation_to(a)**(1/3) * q0

We create a helper function for building a spline, because we will re-use the same thing further below:

def create_closed_shoemake_curve(rotations):
    rotations = rotations + rotations[:2]
    control_points = []
    for q_1, q0, q1 in zip(rotations, rotations[1:], rotations[2:]):
        b, a = shoemake_corrected_control_quaternions(q_1, q0, q1)
        control_points.extend([b, q0, q0, a])
    control_points = control_points[-2:] + control_points[:-2]
    segments = list(zip(*[iter(control_points)] * 4))
    return DeCasteljau(segments)

We don’t want to worry about end conditions here, therefore we create a closed curve.

Let’s come up with some example rotations and see how Shoemake’s curve-joining method works.

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

s = create_closed_shoemake_curve(rotations)

times = np.linspace(s.grid[0], s.grid[-1], 200)

ani = animate_rotations(s.evaluate(times), figsize=(3, 2))

display_animation(ani, default_mode='loop')

Once Loop Reflect

Joining Curves, Another Approach

uniform Catmull–Rom-like quaternion splines

splines.quaternion.CatmullRom

from splines.quaternion import CatmullRom

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

ani = animate_rotations({
    "Shoemake's method": s.evaluate(times),
    'Other approach': cr.evaluate(times),
}, figsize=(5, 2))

display_animation(ani, default_mode='loop')

Once Loop Reflect

There is no visible difference between the two approaches. They are not quite identical, though:

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

0.0082592514906695