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

Non-Uniform Catmull–Rom-Like Rotation Splines#

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

Shoemake [Sho85], section 6: “Questions”

In the uniform case we have used De Casteljau’s algorithm with Slerp to create a “cubic” rotation spline. To extend this to the non-uniform case, we can transform the parameter \(t \to \frac{t - t_i}{t_{i+1} - t_i}\) for each spline segment – as shown in the notebook about non-uniform Euclidean Bézier splines. This is implemented in the class splines.quaternion.DeCasteljau.

Assuming the control points at the start and the end of each segment are given (from a sequence of quaternions to be interpolated), we’ll also need a way to calculate the missing two control points. For inspiration, we can have a look at the notebook about non-uniform (Euclidean) Catmull–Rom splines which provides these equations:

\begin{align*} \boldsymbol{v}_i &= \frac{ \boldsymbol{x}_{i+1} - \boldsymbol{x}_i }{ t_{i+1} - t_i } \\ \boldsymbol{\dot{x}}_i &= \frac{ (t_{i+1} - t_i) \, \boldsymbol{v}_{i-1} + (t_i - t_{i-1}) \, \boldsymbol{v}_i }{ t_{i+1} - t_{i-1} } \\ \boldsymbol{\tilde{x}}_i^{(+)} &= \boldsymbol{x}_i + \frac{(t_{i+1} - t_i) \, \boldsymbol{\dot{x}}_i}{3} \\ \boldsymbol{\tilde{x}}_i^{(-)} &= \boldsymbol{x}_i - \frac{(t_i - t_{i-1}) \, \boldsymbol{\dot{x}}_i}{3} \end{align*}

With the relative rotation \(\delta_i = q_{i+1} {q_i}^{-1}\) we can try to “translate” this to quaternions (using some vector operations in the tangent space):

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

where \(\vec{\rho}_{i}\) is the angular velocity along the great arc from \(q_i\) to \(q_{i+1}\) within the parameter interval from \(t_i\) to \(t_{i+1}\) and \(\vec{\omega}_i\) is the angular velocity of the Catmull–Rom-like quaternion curve at the control point \(q_i\) (which is reached at parameter value \(t_i\)). Finally, \(\tilde{q}_i^{(-)}\) and \(\tilde{q}_i^{(+)}\) are the Bézier-like control quaternions before and after \(q_i\), respectively.

[1]:
from splines.quaternion import UnitQuaternion

def cr_control_quaternions(qs, ts):
    q_1, q0, q1 = qs
    t_1, t0, t1 = ts
    rho_in = q_1.rotation_to(q0).log_map() / (t0 - t_1)
    rho_out = q0.rotation_to(q1).log_map() / (t1 - t0)
    w0 = ((t1 - t0) * rho_in + (t0 - t_1) * rho_out) / (t1 - t_1)
    return [
        UnitQuaternion.exp_map(-w0 * (t0 - t_1) / 3) * q0,
        UnitQuaternion.exp_map(w0 * (t1 - t0) / 3) * q0,
    ]

This approach is also implemented in the class splines.quaternion.CatmullRom.

To illustrate this, let’s load NumPy, a few helpers from helper.py and splines.quaternion.canonicalized().

[2]:
import numpy as np
np.set_printoptions(precision=4)
from helper import angles2quat, animate_rotations, display_animation
from splines.quaternion import canonicalized

The following function can create a closed spline using the above method to calculate control quaternions.

[3]:
from splines.quaternion import DeCasteljau

def catmull_rom_curve(rotations, grid):
    """Create a closed Catmull-Rom-like quaternion curve."""
    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, *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 = cr_control_quaternions(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)

To try this out, we need a few example quaternions and time instances:

[4]:
rotations1 = [
    angles2quat(0, 0, 180),
    angles2quat(0, 45, 90),
    angles2quat(90, 45, 0),
    angles2quat(90, 90, -90),
    angles2quat(180, 0, -180),
    angles2quat(-90, -45, 180),
]
[5]:
grid1 = 0, 0.5, 2, 5, 6, 7, 9
[6]:
cr = catmull_rom_curve(rotations1, grid1)
[7]:
def evaluate(spline, frames=200):
    times = np.linspace(
        spline.grid[0], spline.grid[-1], frames, endpoint=False)
    return spline.evaluate(times)
[8]:
ani = animate_rotations(evaluate(cr))
[9]:
display_animation(ani, default_mode='loop')

Parameterization#

Instead of choosing arbitrary time intervals between control quaternions (via the grid argument), we can calculate time intervals based on the control quaternions themselves.

[10]:
rotations2 = [
    angles2quat(90, 0, -45),
    angles2quat(179, 0, 0),
    angles2quat(181, 0, 0),
    angles2quat(270, 0, -45),
    angles2quat(0, 90, 90),
]

We have seen uniform parameterization already in the previous notebook, where each parameter interval is set to 1:

[11]:
uniform = catmull_rom_curve(rotations2, grid=range(len(rotations2) + 1))

For chordal parameterization of Euclidean splines, we used the Euclidean distance as basis for calculating the time intervals. For rotation splines, it makes more sense to use rotation angles, which are proportional to the lengths of the great arcs between control quaternions:

[12]:
angles = np.array([
    a.rotation_to(b).angle
    for a, b in zip(rotations2, rotations2[1:] + rotations2[:1])])
angles
[12]:
array([1.7027, 0.0349, 1.7027, 2.5936, 1.7178])

The values are probably easier to understand when we show them in degrees:

[13]:
np.degrees(angles)
[13]:
array([ 97.5592,   2.    ,  97.5592, 148.6003,  98.4211])
[14]:
chordal_grid = np.concatenate([[0], np.cumsum(angles)])
[15]:
chordal = catmull_rom_curve(rotations2, grid=chordal_grid)

For centripetal parameterization of Euclidean splines, we used the square root of the Euclidean distances, here we use the square root of the rotation angles:

[16]:
centripetal_grid = np.concatenate([[0], np.cumsum(np.sqrt(angles))])
[17]:
centripetal = catmull_rom_curve(rotations2, grid=centripetal_grid)
[18]:
ani = animate_rotations({
    'uniform': evaluate(uniform),
    'centripetal': evaluate(centripetal),
    'chordal': evaluate(chordal),
})
[19]:
display_animation(ani, default_mode='loop')

The class splines.quaternion.CatmullRom provides a parameter alpha that allows arbitrary parameterization between uniform and chordal – see also parameterized parameterization of Euclidean splines.