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

Non-Uniform Catmull–Rom Splines

[CR74] describes only the uniform case, but it is straightforward to extend the method to non-uniform splines.

The method comprises using three linear interpolations (and extrapolations) between neighboring pairs of the four relevant control points and then blending 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 reversed. This means that equivalently, two (overlapping) quadratic Lagrange interpolations can be used, followed by linearly blending the two resulting points.

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

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

import sympy as sp
sp.init_printing()

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 interpolated 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*.

    """
    assert len(xs) == len(ts)
    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.

p4.simplify()

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

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} \left(- \boldsymbol{x}_{4} + \boldsymbol{x}_{5}\right) + \Delta_{4}^{2} \left(- \boldsymbol{x}_{3} + \boldsymbol{x}_{4}\right)}{\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} \boldsymbol{x}_{5} + \Delta_{4}^{2} \boldsymbol{x}_{6} - \Delta_{5}^{2} \boldsymbol{x}_{4} + \Delta_{5}^{2} \boldsymbol{x}_{5}}{\Delta_{4} \Delta_{5} \left(\Delta_{4} + \Delta_{5}\right)}$

Both results lead to the same general expression:

\begin{equation*} \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}) } \end{equation*}

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}\).

[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 is not describing the tangents of Overhauser splines (as presented in [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 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

Basis Matrix

We already have the correct result, but if we want to derive our basis matrix, we have to re-scale this a bit. The parameter is supposed to go from \(0\) to \(1\) instead of from \(t_4\) to \(t_5\):

p4_normalized = p4.expr.subs(t, t * (t5 - t4) + t4)

M_CR = NamedMatrix(
    r'{M_{\text{CR},4}}',
    sp.Matrix([[c.expand().coeff(x).factor() for x in (x3, x4, x5, x6)]
               for c in p4_normalized.as_poly(t).all_coeffs()]))
M_CR.subs(deltas).simplify()

$\displaystyle {M_{\text{CR},4}} = \left[\begin{matrix}- \frac{\Delta_{4}^{2}}{\Delta_{3} \left(\Delta_{3} + \Delta_{4}\right)} & \frac{\Delta_{4} \left(\Delta_{3} + \Delta_{4} + \Delta_{5}\right)}{\Delta_{3} \left(\Delta_{4} + \Delta_{5}\right)} & - \frac{\Delta_{4} \left(\Delta_{3} + \Delta_{4} + \Delta_{5}\right)}{\Delta_{5} \left(\Delta_{3} + \Delta_{4}\right)} & \frac{\Delta_{4}^{2}}{\Delta_{5} \left(\Delta_{4} + \Delta_{5}\right)}\\\frac{2 \Delta_{4}^{2}}{\Delta_{3} \left(\Delta_{3} + \Delta_{4}\right)} & - \frac{\Delta_{4} \left(\Delta_{3} + 2 \Delta_{4} + 2 \Delta_{5}\right)}{\Delta_{3} \left(\Delta_{4} + \Delta_{5}\right)} & \frac{\Delta_{4} \left(\Delta_{3} + \Delta_{4} + 2 \Delta_{5}\right)}{\Delta_{5} \left(\Delta_{3} + \Delta_{4}\right)} & - \frac{\Delta_{4}^{2}}{\Delta_{5} \left(\Delta_{4} + \Delta_{5}\right)}\\- \frac{\Delta_{4}^{2}}{\Delta_{3} \left(\Delta_{3} + \Delta_{4}\right)} & \frac{- \Delta_{3} + \Delta_{4}}{\Delta_{3}} & \frac{\Delta_{3}}{\Delta_{3} + \Delta_{4}} & 0\\0 & 1 & 0 & 0\end{matrix}\right]$

We can try to manually simplify this a bit more:

M_CR.subs(deltas).simplify().subs([[e.factor(), e] for e in [
    delta4 / (delta4 + delta5) + delta4 / delta3,
    -delta4 / (delta3 + delta4) - delta4 / delta5,
    -delta4 / (delta4 + delta5) - 2 * delta4 / delta3,
    2 * delta4 / (delta3 + delta4) + delta4 / delta5,
]])

$\displaystyle {M_{\text{CR},4}} = \left[\begin{matrix}- \frac{\Delta_{4}^{2}}{\Delta_{3} \left(\Delta_{3} + \Delta_{4}\right)} & \frac{\Delta_{4}}{\Delta_{4} + \Delta_{5}} + \frac{\Delta_{4}}{\Delta_{3}} & - \frac{\Delta_{4}}{\Delta_{3} + \Delta_{4}} - \frac{\Delta_{4}}{\Delta_{5}} & \frac{\Delta_{4}^{2}}{\Delta_{5} \left(\Delta_{4} + \Delta_{5}\right)}\\\frac{2 \Delta_{4}^{2}}{\Delta_{3} \left(\Delta_{3} + \Delta_{4}\right)} & - \frac{\Delta_{4}}{\Delta_{4} + \Delta_{5}} - \frac{2 \Delta_{4}}{\Delta_{3}} & \frac{2 \Delta_{4}}{\Delta_{3} + \Delta_{4}} + \frac{\Delta_{4}}{\Delta_{5}} & - \frac{\Delta_{4}^{2}}{\Delta_{5} \left(\Delta_{4} + \Delta_{5}\right)}\\- \frac{\Delta_{4}^{2}}{\Delta_{3} \left(\Delta_{3} + \Delta_{4}\right)} & \frac{- \Delta_{3} + \Delta_{4}}{\Delta_{3}} & \frac{\Delta_{3}}{\Delta_{3} + \Delta_{4}} & 0\\0 & 1 & 0 & 0\end{matrix}\right]$

We can even introduce two new symbols in order to simplify it yet a bit more:

phi, psi = sp.symbols('Phi4 Psi4')
phi_psi_subs = {
    phi: delta4 / (delta3 + delta4),
    psi: delta4 / (delta4 + delta5),
}
phi_psi_subs

$\displaystyle \left\{ \Phi_{4} : \frac{\Delta_{4}}{\Delta_{3} + \Delta_{4}}, \ \Psi_{4} : \frac{\Delta_{4}}{\Delta_{4} + \Delta_{5}}\right\}$

sp.Matrix([
    [
        -phi * delta4 / delta3,
        psi + delta4 / delta3,
        -phi - (delta4 / delta5),
        psi * delta4 / delta5,
    ], [
        phi * 2 * delta4 / delta3,
        -psi - 2 * delta4 / delta3,
        2 * phi + delta4 / delta5,
        -psi * delta4 / delta5,
    ], [
        -phi * delta4 / delta3,
        (delta4 - delta3) / delta3,
        delta3 / (delta3 + delta4),
        0
    ], [0, 1, 0, 0]
])

$\displaystyle \left[\begin{matrix}- \frac{\Delta_{4} \Phi_{4}}{\Delta_{3}} & \Psi_{4} + \frac{\Delta_{4}}{\Delta_{3}} & - \frac{\Delta_{4}}{\Delta_{5}} - \Phi_{4} & \frac{\Delta_{4} \Psi_{4}}{\Delta_{5}}\\\frac{2 \Delta_{4} \Phi_{4}}{\Delta_{3}} & - \Psi_{4} - \frac{2 \Delta_{4}}{\Delta_{3}} & \frac{\Delta_{4}}{\Delta_{5}} + 2 \Phi_{4} & - \frac{\Delta_{4} \Psi_{4}}{\Delta_{5}}\\- \frac{\Delta_{4} \Phi_{4}}{\Delta_{3}} & \frac{- \Delta_{3} + \Delta_{4}}{\Delta_{3}} & \frac{\Delta_{3}}{\Delta_{3} + \Delta_{4}} & 0\\0 & 1 & 0 & 0\end{matrix}\right]$

assert sp.simplify(
    _.subs(phi_psi_subs) - M_CR.expr.subs(deltas)) == sp.Matrix.zeros(4, 4)

Just to make sure that \(M_{\text{CR},i}\) is consistent with the result from uniform Catmull–Rom splines, let’s set all \(\Delta_i\) to \(1\):

uniform = {
    t3: 3,
    t4: 4,
    t5: 5,
    t6: 6,
    M_CR.name: sp.Symbol(r'{M_\text{CR,uniform}}'),
}

M_CR.subs(uniform).pull_out(sp.S.Half)

$\displaystyle {M_\text{CR,uniform}} = \frac{1}{2} \left[\begin{matrix}-1 & 3 & -3 & 1\\2 & -5 & 4 & -1\\-1 & 0 & 1 & 0\\0 & 2 & 0 & 0\end{matrix}\right]$