This page was generated from doc/rotation/catmull-rom-non-uniform.ipynb. Interactive online version:

Non-Uniform Catmull–Rom-Like Rotation Splines§

What is the best way to allow varying intervals between sequence points in parameter space?

—[Sho85], section 6: “Questions”

notebook about non-uniform (Euclidean) Catmull–Rom splines

\begin{align*} \vec{v}_i &= \frac{ \vec{v}_{i,\text{in}} (t_{i+1} - t_i) + \vec{v}_{i,\text{out}} (t_i - t_{i-1}) }{ t_{i+1} - t_{i-1} } \\ \\ \vec{v}_{i,\text{in}} &= \frac{\vec{x}_i - \vec{x}_{i-1}}{t_i - t_{i-1}} \\ \vec{v}_{i,\text{out}} &= \frac{\vec{x}_{i+1} - \vec{x}_i}{t_{i+1} - t_i} \end{align*}

“translated” to quaternions …

\begin{align*} \tilde{q}_i^{(-)} &= \exp\left(\vec{\omega}_i \frac{t_i - t_{i-1}}{3}\right)^{-1} q_i \\ \tilde{q}_i^{(+)} &= \exp\left(\vec{\omega}_i \frac{t_{i+1} - t_i}{3}\right) q_i \\ \\ \vec{\omega}_i &= \frac{ \vec{\omega}_{i,\text{in}} (t_{i+1} - t_i) + \vec{\omega}_{i,\text{out}} (t_i - t_{i-1}) }{ t_{i+1} - t_{i-1} } \\ \vec{\omega}_{i,\text{in}} &= \frac{\ln(q_{i,\text{in}})}{t_i - t_{i-1}} \\ \vec{\omega}_{i,\text{out}} &= \frac{\ln(q_{i,\text{out}})}{t_{i+1} - t_i} \\ q_{i,\text{in}} &= q_i {q_{i-1}}^{-1} \\ q_{i,\text{out}} &= q_{i+1} {q_i}^{-1} \end{align*}

factor of \(\frac{1}{3}\) because we are dealing with cubic splines

[1]:

from splines.quaternion import UnitQuaternion

[2]:

def control_quaternions1(qs, ts):
    q_1, q0, q1 = qs
    t_1, t0, t1 = ts
    q_in = q0 * q_1.inverse()
    q_out = q1 * q0.inverse()
    w_in = q_in.log_map() / (t0 - t_1)
    w_out = q_out.log_map() / (t1 - t0)
    w0 = ((t1 - t0) * w_in + (t0 - t_1) * w_out) / (t1 - t_1)
    return [
        UnitQuaternion.exp_map(-w0 * (t0 - t_1) / 3) * q0,
        UnitQuaternion.exp_map(w0 * (t1 - t0) / 3) * q0,
    ]

similar, but not quite identical:

[3]:

def control_quaternions2(qs, ts):
    q_1, q0, q1 = qs
    t_1, t0, t1 = ts
    q_in = (q0 * q_1.inverse())**((t1 - t0) / (3 * (t0 - t_1)))
    q_out = (q1 * q0.inverse())**((t0 - t_1) / (3 * (t1 - t0)))
    q_tangent = ((q_out * q_in.inverse())**(1 / 2) * q_in)**2
    return [
        (q_tangent**((t0 - t_1) / (t1 - t_1))).inverse() * q0,
        q_tangent**((t1 - t0) / (t1 - t_1)) * q0,
    ]

[4]:

import numpy as np

helper.py

[5]:

from helper import angles2quat, animate_rotations, display_animation

splines.quaternion.DeCasteljau class

[6]:

from splines.quaternion import DeCasteljau, canonicalized

[7]:

def create_closed_curve(rotations, grid, control_quaternion_func):
    assert len(rotations) + 1 == len(grid)
    rotations = rotations[-1:] + rotations + rotations[:2]
    # Avoid angles of more than 180 degrees (including the added rotations):
    rotations = list(canonicalized(rotations))
    first_interval = grid[1] - grid[0]
    last_interval = grid[-1] - grid[-2]
    extended_grid = [grid[0] - last_interval] + list(grid) + [grid[-1] + first_interval]
    control_points = []
    for qs, ts in zip(
            zip(rotations, rotations[1:], rotations[2:]),
            zip(extended_grid, extended_grid[1:], extended_grid[2:])):
        q_before, q_after = control_quaternion_func(qs, ts)
        control_points.extend([q_before, qs[1], qs[1], q_after])
    control_points = control_points[2:-2]
    segments = list(zip(*[iter(control_points)] * 4))
    return DeCasteljau(segments, grid)

[8]:

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

[9]:

grid = np.array([0, 0.5, 2, 5, 6, 7, 10])

[10]:

s1 = create_closed_curve(rotations, grid, control_quaternions1)
s2 = create_closed_curve(rotations, grid, control_quaternions2)

[11]:

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

for comparison, Barry–Goldman

[12]:

from splines.quaternion import BarryGoldman

[13]:

bg = BarryGoldman(rotations, grid)

[14]:

ani = animate_rotations({
    '1': evaluate(s1),
    '2': evaluate(s2),
    'Barry–Goldman': evaluate(bg),
}, figsize=(6, 2))

[15]:

display_animation(ani, default_mode='loop')

[16]:

max(max(map(abs, q.xyzw)) for q in (evaluate(s1) - evaluate(s2)))

[16]:

0.01051294090598398

Parameterization§

[17]:

rotations = [
    angles2quat(90, 0, -45),
    angles2quat(179, 0, 0),
    angles2quat(181, 0, 0),
    angles2quat(270, 0, -45),
    angles2quat(0, 90, 90),
]

[18]:

uniform2 = create_closed_curve(rotations, range(len(rotations) + 1), control_quaternions1)

chordal parameterization

[19]:

angles = np.array([
    np.arccos(a.dot(b))
    #np.arccos(a.dot(b)) * 2
    #a.dot(b)
    #(b * a.inverse()).angle
    #(b - a).norm
    for a, b in zip(rotations, rotations[1:] + rotations[:1])])
angles

[19]:

array([0.85136439, 0.01745329, 0.85136439, 1.29678212, 0.85888576])

[20]:

chordal_grid = np.concatenate([[0], np.cumsum(angles)])

[21]:

chordal2 = create_closed_curve(rotations, chordal_grid, control_quaternions1)

centripetal parameterization

[22]:

centripetal_grid = np.concatenate([[0], np.cumsum(np.sqrt(angles))])

[23]:

centripetal2 = create_closed_curve(rotations, centripetal_grid, control_quaternions1)

[24]:

ani = animate_rotations({
    'uniform': evaluate(uniform2),
    'chordal': evaluate(chordal2),
    'centripetal': evaluate(centripetal2),
}, figsize=(6, 2))

[25]:

display_animation(ani, default_mode='loop')

The other method is very similar:

[26]:

uniform1 = create_closed_curve(rotations, range(len(rotations) + 1), control_quaternions2)
chordal1 = create_closed_curve(rotations, chordal_grid, control_quaternions2)
centripetal1 = create_closed_curve(rotations, centripetal_grid, control_quaternions2)

[27]:

max(max(map(abs, q.xyzw)) for q in (evaluate(uniform1) - evaluate(uniform2)))

[27]:

0.0026683373798292997

[28]:

max(max(map(abs, q.xyzw)) for q in (evaluate(chordal1) - evaluate(chordal2)))

[28]:

0.002259641902391918

[29]:

max(max(map(abs, q.xyzw)) for q in (evaluate(centripetal1) - evaluate(centripetal2)))

[29]:

0.002486960572741559