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

Non-Uniform Catmull–Rom Splines

Catmull and Rom [CR74] describe only the uniform case, but it is straightforward to extend their method to non-uniform splines.

The method creates three linear interpolations (and extrapolations) between neighboring pairs of the four relevant control points and then blends the three resulting points with a quadratic B-spline basis function.

As we have seen in the notebook about uniform Catmull–Rom splines and as we will again see in the notebook about the Barry–Goldman algorithm, the respective degrees can be swapped. This means that equivalently, two (overlapping) quadratic Lagrange interpolations can be used, followed by linearly blending the two resulting points.

Since the latter is both easier to implement and easier to wrap one’s head around, we’ll use it in the following derivations.

We will derive the tangent vectors at the segment boundaries, which will later serve as a starting point for deriving non-uniform Kochanek–Bartels splines. See the notebook about the Barry–Goldman algorithm for an alternative (but closely related) derivation.

import sympy as sp
sp.init_printing()

As usual, we look at the fifth polynomial segment \(\boldsymbol{p}_4(t)\) from \(\boldsymbol{x}_4\) to \(\boldsymbol{x}_5\), where \(t_4 \le t \le t_5\). Later, we will generalize this to an arbitrary polynomial segment \(\boldsymbol{p}_i(t)\) from \(\boldsymbol{x}_i\) to \(\boldsymbol{x}_{i+1}\), where \(t_i \le t \le t_{i+1}\).

x3, x4, x5, x6 = sp.symbols('xbm3:7')

t, t3, t4, t5, t6 = sp.symbols('t t3:7')

We use some tools from utility.py:

from utility import NamedExpression, NamedMatrix

As shown in the notebook about Lagrange interpolation, it can be implemented using Neville’s algorithm:

def lerp(xs, ts, t):
    """Linear interpolation.

    Returns the interpolated value at time *t*,
    given the two values *xs* at times *ts*.

    """
    x_begin, x_end = xs
    t_begin, t_end = ts
    return (x_begin * (t_end - t) + x_end * (t - t_begin)) / (t_end - t_begin)

def neville(xs, ts, t):
    """Lagrange interpolation using Neville's algorithm.

    Returns the interpolated value at time *t*,
    given the values *xs* at times *ts*.

    """
    if len(xs) != len(ts):
        raise ValueError('xs and ts must have the same length')
    while len(xs) > 1:
        step = len(ts) - len(xs) + 1
        xs = [
            lerp(*args, t)
            for args in zip(zip(xs, xs[1:]), zip(ts, ts[step:]))]
    return xs[0]

Alternatively, sympy.interpolate() could be used.

We use two overlapping quadratic Lagrange interpolations followed by linear blending:

p4 = NamedExpression(
    'pbm4',
    lerp([
        neville([x3, x4, x5], [t3, t4, t5], t),
        neville([x4, x5, x6], [t4, t5, t6], t),
    ], [t4, t5], t))

Note

Since the two invocations of Neville’s algorithm overlap, some values that are used by both are unnecessarily computed by both. It would be more efficient to calculate each of these values only once.

The Barry–Goldman algorithm avoids this repeated computation.

But here, since we are using symbolic expressions, this doesn’t really matter because the redundant expressions should be simplified away by SymPy.

The following expressions can be simplified by introducing a few new symbols \(\Delta_i\):

delta3, delta4, delta5 = sp.symbols('Delta3:6')
deltas = {
    t4 - t3: delta3,
    t5 - t4: delta4,
    t6 - t5: delta5,
    t5 - t3: delta3 + delta4,
    t6 - t4: delta4 + delta5,
    t6 - t3: delta3 + delta4 + delta5,
    # A few special cases that SymPy has a hard time resolving:
    t4 + t4 - t3: t4 + delta3,
    t6 + t6 - t3: t6 + delta3 + delta4 + delta5,
}

Tangent Vectors

To get the tangent vectors at the control points, we just have to take the first derivative …

pd4 = p4.diff(t)

… and evaluate it at \(t_4\) and \(t_5\):

start_tangent = pd4.evaluated_at(t, t4)
start_tangent.subs(deltas).simplify()

$\displaystyle \left.{\frac{d}{d t} \boldsymbol{p}_{4}}\right\rvert_{t=t_{4}} = \frac{- \Delta_{3}^{2} \boldsymbol{x}_{4} + \Delta_{3}^{2} \boldsymbol{x}_{5} - \Delta_{4}^{2} \boldsymbol{x}_{3} + \Delta_{4}^{2} \boldsymbol{x}_{4}}{\Delta_{3} \Delta_{4} \left(\Delta_{3} + \Delta_{4}\right)}$

end_tangent = pd4.evaluated_at(t, t5)
end_tangent.subs(deltas).simplify()

$\displaystyle \left.{\frac{d}{d t} \boldsymbol{p}_{4}}\right\rvert_{t=t_{5}} = \frac{\Delta_{4}^{2} \left(- \boldsymbol{x}_{5} + \boldsymbol{x}_{6}\right) + \Delta_{5}^{2} \left(- \boldsymbol{x}_{4} + \boldsymbol{x}_{5}\right)}{\Delta_{4} \Delta_{5} \left(\Delta_{4} + \Delta_{5}\right)}$

Both results lead to the same general expression (which is expected, since the incoming and outgoing tangents are supposed to be equal):

\begin{align*} \boldsymbol{\dot{x}}_i &= \frac{ (t_{i+1} - t_i)^2 (\boldsymbol{x}_i - \boldsymbol{x}_{i-1}) + (t_i - t_{i-1})^2 (\boldsymbol{x}_{i+1} - \boldsymbol{x}_i) }{ (t_{i+1} - t_i)(t_i - t_{i-1})(t_{i+1} - t_{i-1}) } \\ &= \frac{ {\Delta_i}^2 (\boldsymbol{x}_i - \boldsymbol{x}_{i-1}) + {\Delta_{i-1}}^2 (\boldsymbol{x}_{i+1} - \boldsymbol{x}_i) }{ \Delta_i \Delta_{i-1} (\Delta_i + \Delta_{i-1}) } \end{align*}

Equivalently, this can be written as:

\begin{align*} \boldsymbol{\dot{x}}_i &= \frac{ (t_{i+1} - t_i) (\boldsymbol{x}_i - \boldsymbol{x}_{i-1}) }{ (t_i - t_{i-1})(t_{i+1} - t_{i-1}) } + \frac{ (t_i - t_{i-1}) (\boldsymbol{x}_{i+1} - \boldsymbol{x}_i) }{ (t_{i+1} - t_i)(t_{i+1} - t_{i-1}) } \\ &= \frac{ \Delta_i (\boldsymbol{x}_i - \boldsymbol{x}_{i-1}) }{ \Delta_{i-1} (\Delta_i + \Delta_{i-1}) } + \frac{ \Delta_{i-1} (\boldsymbol{x}_{i+1} - \boldsymbol{x}_i) }{ \Delta_i (\Delta_i + \Delta_{i-1}) } \end{align*}

An alternative (but very similar) way to derive these tangent vectors is shown in the notebook about the Barry–Goldman algorithm.

And there is yet another way to calculate the tangents, without even needing to obtain a cubic polynomial and its derivative: Since we are using a linear blend of two quadratic polynomials, we know that at the beginning (\(t = t_4\)) only the first quadratic polynomial has an influence and at the end (\(t = t_5\)) only the second quadratic polynomial is relevant. Therefore, to determine the tangent vector at the beginning of the segment, it is sufficient to get the derivative of the first quadratic polynomial.

first_quadratic = neville([x3, x4, x5], [t3, t4, t5], t)

sp.degree(first_quadratic, t)

$\displaystyle 2$

first_quadratic.diff(t).subs(t, t4)

$\displaystyle \frac{\frac{\left(- t_{3} + t_{4}\right) \left(- \boldsymbol{x}_{4} + \boldsymbol{x}_{5}\right)}{- t_{4} + t_{5}} + \frac{\left(- t_{4} + t_{5}\right) \left(- \boldsymbol{x}_{3} + \boldsymbol{x}_{4}\right)}{- t_{3} + t_{4}}}{- t_{3} + t_{5}}$

This can be written as (which is sometimes called the standard three-point difference formula):

\begin{equation*} \boldsymbol{\dot{x}}_i = \frac{ \Delta_i \boldsymbol{v}_{i-1} + \Delta_{i-1} \boldsymbol{v}_i }{ \Delta_{i-1} + \Delta_i }, \end{equation*}

with \(\Delta_i = t_{i+1} - t_i\) and \(\boldsymbol{v}_i = \frac{\boldsymbol{x}_{i+1} - \boldsymbol{x}_i}{\Delta_i}\).

de Boor [dB78] calls this piecewise cubic Bessel interpolation, and it has also been called Bessel tangent method, Overhauser method and Bessel–Overhauser splines.

Note

Even though this formula is commonly associated with the name Overhauser, it does not describe the tangents of Overhauser splines as presented by Overhauser [Ove68].

Long story short, it’s the same as we had above:

assert sp.simplify(_ - start_tangent.expr) == 0

The first derivative of the second quadratic polynomial can be used to get the tangent vector at the end of the segment.

second_quadratic = neville([x4, x5, x6], [t4, t5, t6], t)
second_quadratic.diff(t).subs(t, t5)

$\displaystyle \frac{\frac{\left(- t_{4} + t_{5}\right) \left(- \boldsymbol{x}_{5} + \boldsymbol{x}_{6}\right)}{- t_{5} + t_{6}} + \frac{\left(- t_{5} + t_{6}\right) \left(- \boldsymbol{x}_{4} + \boldsymbol{x}_{5}\right)}{- t_{4} + t_{5}}}{- t_{4} + t_{6}}$

assert sp.simplify(_ - end_tangent.expr) == 0

You might encounter yet another way to write the equation for \(\boldsymbol{\dot{x}}_4\) (e.g. at https://stackoverflow.com/a/23980479/) …

(x4 - x3) / (t4 - t3) - (x5 - x3) / (t5 - t3) + (x5 - x4) / (t5 - t4)

$\displaystyle \frac{- \boldsymbol{x}_{4} + \boldsymbol{x}_{5}}{- t_{4} + t_{5}} - \frac{- \boldsymbol{x}_{3} + \boldsymbol{x}_{5}}{- t_{3} + t_{5}} + \frac{- \boldsymbol{x}_{3} + \boldsymbol{x}_{4}}{- t_{3} + t_{4}}$

… but again, this is equivalent to the equation shown above:

assert sp.simplify(_ - start_tangent.expr) == 0

Using Non-Uniform Bézier Segments

Similar to the uniform case, the above equation for the tangent vectors can be used to construct non-uniform Hermite splines or, after multiplying them with the appropriate parameter interval and dividing them by 3, to obtain the two additional control points for non-uniform cubic Bézier spline segments:

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

This is again using \(\Delta_i = t_{i+1} - t_i\) and \(\boldsymbol{v}_i = \frac{\boldsymbol{x}_{i+1} - \boldsymbol{x}_i}{\Delta_i}\).

x4tilde = x4 + (t5 - t4) * start_tangent.expr / 3

x5tilde = x5 - (t5 - t4) * end_tangent.expr / 3

Using Non-Uniform Quadrangle Interpolation

Just like in the uniform case, we calculate the quadrangle points from the Bézier control points, as shown in the notebook about quadrangle interpolation:

\begin{align*} \boldsymbol{\bar{x}}_i^{(+)} &= \frac{3}{2} \boldsymbol{\tilde{x}}_i^{(+)} - \frac{1}{2} \boldsymbol{x}_{i+1}\\ \boldsymbol{\bar{x}}_i^{(-)} &= \frac{3}{2} \boldsymbol{\tilde{x}}_i^{(-)} - \frac{1}{2} \boldsymbol{x}_{i-1} \end{align*}

x4bar = 3 * x4tilde / 2 - x5 / 2

terms4 = sp.collect(x4bar.expand(), [x3, x4, x5], evaluate=False)

Some manual rewriting leads to this expression:

sp.factor(terms4[x4] + terms4[x5] + terms4[x3]) * x4 - (
    sp.factor(-terms4[x5]) * (x5 - x4) +
    sp.factor(-terms4[x3]) * (x3 - x4))

$\displaystyle \boldsymbol{x}_{4} - \frac{\left(t_{4} - t_{5}\right) \left(- \boldsymbol{x}_{4} + \boldsymbol{x}_{5}\right)}{2 \left(t_{3} - t_{5}\right)} - \frac{\left(t_{4} - t_{5}\right)^{2} \left(\boldsymbol{x}_{3} - \boldsymbol{x}_{4}\right)}{2 \left(t_{3} - t_{4}\right) \left(t_{3} - t_{5}\right)}$

We should make sure that our re-written expression is actually the same as the one we started from:

assert sp.simplify(_ - x4bar) == 0

Now the same for the incoming quadrangle point:

x5bar = 3 * x5tilde / 2 - x4 / 2

terms5 = sp.collect(x5bar.expand(), [x4, x5, x6], evaluate=False)

sp.factor(terms5[x5] + terms5[x6] + terms5[x4]) * x5 - (
    sp.factor(-terms5[x6]) * (x6 - x5) +
    sp.factor(-terms5[x4]) * (x4 - x5))

$\displaystyle \boldsymbol{x}_{5} - \frac{\left(t_{4} - t_{5}\right)^{2} \left(- \boldsymbol{x}_{5} + \boldsymbol{x}_{6}\right)}{2 \left(t_{4} - t_{6}\right) \left(t_{5} - t_{6}\right)} - \frac{\left(t_{4} - t_{5}\right) \left(\boldsymbol{x}_{4} - \boldsymbol{x}_{5}\right)}{2 \left(t_{4} - t_{6}\right)}$

assert sp.simplify(_ - x5bar) == 0

The above expressions can be generalized to (as always with \(\Delta_i = t_{i+1} - t_i\)):

\begin{align*} \boldsymbol{\bar{x}}_i^{(+)} &= \boldsymbol{x}_i - \frac{\Delta_i}{2 (\Delta_{i-1} + \Delta_i)} \left( (\boldsymbol{x}_{i+1} - \boldsymbol{x}_i) + \frac{\Delta_i}{\Delta_{i-1}} (\boldsymbol{x}_{i-1} - \boldsymbol{x}_i) \right)\\ \boldsymbol{\bar{x}}_i^{(-)} &= \boldsymbol{x}_i - \frac{\Delta_{i-1}}{2 (\Delta_{i-1} + \Delta_i)} \left( \frac{\Delta_{i-1}}{\Delta_i} (\boldsymbol{x}_{i+1} - \boldsymbol{x}_i) + (\boldsymbol{x}_{i-1} - \boldsymbol{x}_i) \right) \end{align*}

Animation

To illustrate what two quadratic Lagrange interpolations followed by linear blending might look like, we can generate an animation by means of the file catmull_rom.py, with some help from helper.py:

from catmull_rom import animation_2_1, animation_1_2
from helper import show_animation

vertices = [
    (1, 0),
    (0.5, 1),
    (6, 2),
    (5, 0),
]

times = [
    0,
    1,
    6,
    8,
]

show_animation(animation_2_1(vertices, times))

Once Loop Reflect

In the beginning of this notebook, we claimed that two quadratic interpolations followed by linear blending are easier to understand. To prove this, let’s have a look at what three linear interpolations (and extrapolations) followed by quadratic B-spline blending would look like:

show_animation(animation_1_2(vertices, times))

Once Loop Reflect

Would you agree that this is less straightforward?

If you would rather replace the quadratic B-spline basis function with a bunch of linear interpolations (using De Boor’s algorithm), take a look at the notebook about the Barry–Goldman algorithm.