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

De Casteljau’s Algorithm With Slerp

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

def slerp(one, two, t):
    """Spherical Linear intERPolation."""
    return (two * one.inverse())**t * one

We’ll also need NumPy and a few helpers from helper.py:

import numpy as np
from helper import angles2quat, plot_rotation, plot_rotations
from helper import animate_rotations, display_animation

“Cubic”

Shoemake [Sho85] only talks about the “cubic” case, consisting of three nested applications of Slerp. Since this is done in a curved space, 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 cubic_de_casteljau(q0, q1, q2, q3, t):
    """De Casteljau's algorithm of "degree" 3 using Slerp."""
    if not np.isscalar(t):
        # If t is a list, return a list of unit quaternions
        return [cubic_de_casteljau(q0, q1, q2, q3, t) for t in 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,
    )

To illustrate this, let’s define 4 unit quaternions that we can use as control points:

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

plot_rotation({'q0': q0, 'q1': q1, 'q2': q2, 'q3': q3});

plot_rotations(
    cubic_de_casteljau(q0, q1, q2, q3, np.linspace(0, 1, 9)),
    figsize=(8, 1))

We can see that the curve starts with the first rotation and ends with the last one. The two middle control quaternions q1 and q2 influence the shape of the rotation curve but they are not part of the interpolant themselves.

ani = animate_rotations(
    cubic_de_casteljau(q0, q1, q2, q3, np.linspace(0, 1, 100)))

display_animation(ani, default_mode='reflect')

Once Loop Reflect

Arbitrary “Degree”

The class splines.quaternion.DeCasteljau allows arbitrary numbers of unit quaternions per segment and therefore arbitrary “degrees”:

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])

ani = animate_rotations(s.evaluate(np.linspace(s.grid[0], s.grid[-1], 100)))

display_animation(ani, default_mode='reflect')

Once Loop Reflect

Constant Angular Speed

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

—Shoemake [Sho85], section 6: “Questions”

Remember arc-length parameterization of Euclidean splines? We used the class splines.UnitSpeedAdapter which happens to be implemented in a way that it is also usable for rotation splines, how convenient! The only requirement is that the second derivative of the wrapped spline yields an angular velocity vector, which is nothing else than the instantaneous rotation axis scaled by the angular speed.

from splines import UnitSpeedAdapter

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

s2 = UnitSpeedAdapter(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)),
})

display_animation(ani, default_mode='reflect')

Once Loop Reflect

Joining Curves

Until now, we have assumed that four control quaternions are given for each “cubic” segment.

If a list of quaternions is given, which is supposed to be interpolated, the intermediate control quaternions can be computed from neighboring quaternions as shown in the notebook about uniform Catmull–Rom-like quaternion splines.