This page was generated from doc/rotation/barry-goldman.ipynb. Interactive online version: Binder badge.

Barry–Goldman Algorithm§

We can try to use the Barry–Goldman algorithm for non-uniform Euclidean Catmull–Rom splines using Slerp instead of linear interpolations, just as we have done with De Casteljau’s algorithm.

[1]:
def slerp(one, two, t):
    return (two * one.inverse())**t * one
[2]:
def barry_goldman(rotations, times, t):
    q0, q1, q2, q3 = rotations
    t0, t1, t2, t3 = times
    return slerp(
        slerp(
            slerp(q0, q1, (t - t0) / (t1 - t0)),
            slerp(q1, q2, (t - t1) / (t2 - t1)),
            (t - t0) / (t2 - t0)),
        slerp(
            slerp(q1, q2, (t - t1) / (t2 - t1)),
            slerp(q2, q3, (t - t2) / (t3 - t2)),
            (t - t1) / (t3 - t1)),
        (t - t1) / (t2 - t1))

Example:

[3]:
import numpy as np

helper.py

[4]:
from helper import angles2quat, animate_rotations, display_animation
[5]:
q0 = angles2quat(0, 0, 0)
q1 = angles2quat(90, 0, 0)
q2 = angles2quat(90, 90, 0)
q3 = angles2quat(90, 90, 90)
[6]:
t0 = 0
t1 = 1
t2 = 3
t3 = 3.5
[7]:
frames = 50
[8]:
ani = animate_rotations({
    'Barry–Goldman (q0, q1, q2, q3)': [
        barry_goldman([q0, q1, q2, q3], [t0, t1, t2, t3], t)
        for t in np.linspace(t1, t2, frames)
    ],
    'Slerp (q1, q2)': slerp(q1, q2, np.linspace(0, 1, frames)),
}, figsize=(6, 2))
[9]:
display_animation(ani, default_mode='once')

splines.quaternion.BarryGoldman class

[10]:
from splines.quaternion import BarryGoldman
[11]:
import numpy as np

helper.py

[12]:
from helper import angles2quat, animate_rotations, display_animation
[13]:
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),
]
[14]:
grid = np.array([0, 0.5, 2, 5, 6, 7, 9])
[15]:
bg = BarryGoldman(rotations, grid)

For comparison … Catmull–Rom-like quaternion spline

splines.quaternion.CatmullRom class

[16]:
from splines.quaternion import CatmullRom
[17]:
cr = CatmullRom(rotations, grid, endconditions='closed')
[18]:
def evaluate(spline, samples=200):
    times = np.linspace(spline.grid[0], spline.grid[-1], samples, endpoint=False)
    return spline.evaluate(times)
[19]:
ani = animate_rotations({
    'Barry–Goldman': evaluate(bg),
    'Catmull–Rom-like': evaluate(cr),
}, figsize=(4, 2))
[20]:
display_animation(ani, default_mode='loop')
[21]:
rotations = [
    angles2quat(90, 0, -45),
    angles2quat(179, 0, 0),
    angles2quat(181, 0, 0),
    angles2quat(270, 0, -45),
    angles2quat(0, 90, 90),
]
[22]:
s_uniform = BarryGoldman(rotations)
s_chordal = BarryGoldman(rotations, alpha=1)
s_centripetal = BarryGoldman(rotations, alpha=0.5)
[23]:
ani = animate_rotations({
    'uniform': evaluate(s_uniform, samples=300),
    'chordal': evaluate(s_chordal, samples=300),
    'centripetal': evaluate(s_centripetal, samples=300),
}, figsize=(6, 2))
[24]:
display_animation(ani, default_mode='loop')

Constant Angular Speed§

Not very efficient, De Casteljau’s algorithm is faster because it directly provides the tangent.

[25]:
from splines import ConstantSpeedAdapter
[26]:
class BarryGoldmanWithDerivative(BarryGoldman):

    delta_t = 0.000001

    def evaluate(self, t, n=0):
        """Evaluate quaternion or angular velocity."""
        if not np.isscalar(t):
            return np.array([self.evaluate(t, n) for t in t])
        if n == 0:
            return super().evaluate(t)
        elif n == 1:
            # NB: We move the interval around because
            #     we cannot access times before and after
            #     the first and last time, respectively.
            fraction = (t - self.grid[0]) / (self.grid[-1] - self.grid[0])
            before = super().evaluate(t - fraction * self.delta_t)
            after = super().evaluate(t + (1 - fraction) * self.delta_t)
            # NB: Double angle
            return (after * before.inverse()).log_map() * 2 / self.delta_t
        else:
            raise ValueError('Unsupported n: {!r}'.format(n))
[27]:
s = ConstantSpeedAdapter(BarryGoldmanWithDerivative(rotations, alpha=0.5))

Takes a long time!

[28]:
ani = animate_rotations({
    'non-constant speed': evaluate(s_centripetal),
    'constant speed': evaluate(s),
}, figsize=(4, 2))
[29]:
display_animation(ani, default_mode='loop')