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

Properties of Catmull–Rom Splines

Catmull and Rom [CR74] present a whole class of splines with a whole range of properties. Here we only consider one member of this class which is a cubic polynomial interpolating spline with \(C^1\) continuity and local support. Nowadays, this specific case is typically simply referred to as Catmull–Rom spline.

This type of spline is very popular because they are very easy to use. Only a sequence of control points has to be specified, the corresponding tangents are calculated automatically from the given points. Using those tangents, the spline can be implemented using cubic Hermite splines. Alternatively, spline values can be directly calculated with the Barry–Goldman algorithm.

To calculate the spline values between two control points, the preceding and the following control points are needed as well. The tangent vector at any given control point can be calculated from this control point, its predecessor and its successor. Since Catmull–Rom splines are \(C^1\) continuous, incoming and outgoing tangent vectors are equal.

The following examples use the Python class splines.CatmullRom to create both uniform and non-uniform splines. Only closed splines are shown, other end conditions can also be used, but they are not specific to this type of spline.

import matplotlib.pyplot as plt
import numpy as np
np.set_printoptions(precision=4)

Apart from the splines module …

import splines

… we also import a few helper functions from helper.py:

from helper import plot_spline_2d, plot_tangent_2d

Let’s choose a few points for an example:

points1 = [
    (-1, -0.5),
    (0, 2.3),
    (1, 1),
    (4, 1.3),
    (3.8, -0.2),
    (2.5, 0.1),
]

Without specifying any time values, we get a uniform spline:

s1 = splines.CatmullRom(points1, endconditions='closed')

fig, ax = plt.subplots()
plot_spline_2d(s1, ax=ax)

Tangent Vectors

In the uniform case, the tangent vectors at any given control point are parallel to the line connecting the preceding point and the following point. The tangent vector has the same orientation as that line but only half its length. In other (more mathematical) words:

\begin{equation*} \boldsymbol{\dot{x}}_i = \frac{\boldsymbol{x}_{i+1} - \boldsymbol{x}_{i-1}}{2} \end{equation*}

This is illustrated for two control points in the following plot:

for idx, color in zip([2, 5], ['purple', 'hotpink']):
    plot_tangent_2d(
        s1.evaluate(s1.grid[idx], 1),
        s1.evaluate(s1.grid[idx]), color=color, ax=ax)
    ax.plot(
        *s1.evaluate([s1.grid[idx - 1], s1.grid[idx + 1]]).T,
        '--', color=color, linewidth=2)
fig

We can see here that each tangent vector is parallel to and has half the length of the line connecting the preceding and the following vertex, just as promised.

However, this will not be true anymore if we are using non-uniform time instances:

times2 = 0, 1, 2.2, 3, 4, 4.5, 6

s2 = splines.CatmullRom(points1, grid=times2, endconditions='closed')

plot_spline_2d(s2, ax=ax)
for idx, color in zip([2, 5], ['green', 'crimson']):
    plot_tangent_2d(
        s2.evaluate(s2.grid[idx], 1),
        s2.evaluate(s2.grid[idx]), color=color, ax=ax)
fig

In the non-uniform case, the equation for the tangent vector gets quite a bit more complicated:

\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*}

The derivation of this equation is shown in a separate notebook.

Equivalently, this can be written as:

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

Also equivalently, with \(\boldsymbol{v}_i = \frac{\boldsymbol{x}_{i+1} - \boldsymbol{x}_i}{t_{i+1} - t_i}\), it can be written as:

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

Wrong Tangent Vectors

Some sources provide a simpler equation which is different from the tangent vector of a Catmull–Rom spline (except in the uniform case):

\begin{equation*} \boldsymbol{\dot{x}}_i \overset{?}{=} \frac{\boldsymbol{v}_{i-1} + \boldsymbol{v}_i}{2} = \frac{1}{2} \left( \frac{\boldsymbol{x}_i - \boldsymbol{x}_{i-1}}{t_i - t_{i-1}} + \frac{\boldsymbol{x}_{i + 1} - \boldsymbol{x}_i}{t_{i + 1} - t_i} \right) \end{equation*}

class MeanVelocity(splines.CatmullRom):

    @staticmethod
    def _calculate_tangent(points, times):
        x_1, x0, x1 = np.asarray(points)
        t_1, t0, t1 = times
        v_1 = (x0 - x_1) / (t0 - t_1)
        v0 = (x1 - x0) / (t1 - t0)
        return (v_1 + v0) / 2

Until April 2023, Wikipedia showed yet a simpler equation. They mentioned that “this assumes uniform parameter spacing”, but since \(t_{i - 1}\) and \(t_{i + 1}\) appeared in the equation, it might be tempting to use it for the non-uniform case as well. We’ll see below how that turns out.

The authors of the page don’t seem to have been quite sure about this equation, because it has changed over time. This was shown until mid-2021:

\begin{equation*} \boldsymbol{\dot{x}}_i \overset{?}{=} \frac{\boldsymbol{x}_{i + 1} - \boldsymbol{x}_{i - 1}}{t_{i + 1} - t_{i - 1}} \end{equation*}

class Wikipedia1(splines.CatmullRom):

    @staticmethod
    def _calculate_tangent(points, times):
        x_1, _, x1 = np.asarray(points)
        t_1, _, t1 = times
        return (x1 - x_1) / (t1 - t_1)

And this slight variation was shown since then until April 2023:

\begin{equation*} \boldsymbol{\dot{x}}_i \overset{?}{=} \frac{1}{2}\frac{\boldsymbol{x}_{i + 1} - \boldsymbol{x}_{i - 1}}{t_{i + 1} - t_{i - 1}} \end{equation*}

class Wikipedia2(splines.CatmullRom):

    @staticmethod
    def _calculate_tangent(points, times):
        x_1, _, x1 = np.asarray(points)
        t_1, _, t1 = times
        return (1/2) * (x1 - x_1) / (t1 - t_1)

The first one is correct in the uniform case (which the Wikipedia page assumes), but not in the general non-uniform case, as we’ll see in a moment.

The second one is obviously wrong in the case where all intervals are of length \(1\) (i.e. \(t_{i+1} - t_i = t_i - t_{i-1} = 1\)):

\begin{equation*} \frac{\boldsymbol{x}_{i + 1} - \boldsymbol{x}_{i - 1}}{4} \ne \frac{\boldsymbol{x}_{i + 1} - \boldsymbol{x}_{i - 1}}{2} = \boldsymbol{\dot{x}}_i \end{equation*}

Since April 2023, the page is showing the correct equation for the uniform case.

The X3D standard (version 3.3) even suggests to use different incoming and outgoing tangents, which destroys \(C^1\) continuity!

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

class X3D(splines.KochanekBartels):
    # We derive from KochanekBartels because the
    # incoming and outgoing tangents are different:
    @staticmethod
    def _calculate_tangents(points, times, _ignored):
        x_1, _, x1 = np.asarray(points)
        t_1, t0, t1 = times
        incoming = (t1 - t0) * (x1 - x_1) / (t1 - t_1)
        outgoing = (t0 - t_1) * (x1 - x_1) / (t1 - t_1)
        return incoming, outgoing

To illustrate the different choices of tangent vectors, we use the vertex data from Lee [Lee89], figure 6:

points3 = [
    (0, 0),
    (10, 25),
    (10, 24),
    (11, 24.5),
    (33, 25),
]

Deciding between “right” and “wrong” tangent vectors is surprisingly hard, because most of the options look somewhat reasonable in most cases. However, we can try to use quite extreme vertex positions and we can use centripetal parameterization (see below) and check if its guaranteed properties hold for different choices of tangent vectors.

def plot_spline(cls, linestyle='-', **args):
    # alpha=0.5 => centripetal parameterization
    spline = cls(points3, alpha=0.5)
    plot_spline_2d(
        spline, label=cls.__name__, chords=False,
        marker=None, linestyle=linestyle, **args)

plot_spline(MeanVelocity, linestyle=':')
plot_spline(X3D, linestyle='-.')
plot_spline(Wikipedia1)
plot_spline(Wikipedia2, linestyle='--')
plot_spline(splines.CatmullRom, linewidth=3)
plt.axis([9, 13, 23.9, 25.6])
plt.legend();

As we can immediately see, the tangents from X3D are utterly wrong and the first one from Wikipedia is also quite obviously broken. The other two don’t look too bad, but they slightly overshoot, and according to Yuksel et al. [YSK11] that is something that centripetal Catmull–Rom splines are guaranteed not to do.

Again, to be fair to the Wikipedia article’s authors, they mentioned that uniform parameter spacing is assumed, so their equation is not supposed to be used in this non-uniform context. The equation has been changed in the meantime to avoid confusion.

Cusps and Self-Intersections

Uniform parametrization typically works very well if the (Euclidean) distances between consecutive vertices are all similar. However, if the distances are very different, the shape of the spline often turns out to be unexpected. Most notably, in extreme cases there might be even cusps or self-intersections within a spline segment.

def plot_catmull_rom(*args, **kwargs):
    plot_spline_2d(splines.CatmullRom(*args, endconditions='closed', **kwargs))

points4 = [
    (0, 0),
    (0, 0.5),
    (1.5, 1.5),
    (1.6, 1.5),
    (3, 0.2),
    (3, 0),
]

plot_catmull_rom(points4)

We can try to compensate this by manually selecting some non-uniform time instances:

times4 = 0, 0.2, 0.9, 1, 3, 3.3, 4.5

plot_catmull_rom(points4, times4)

Time values can be chosen by trial and error, but there are also ways to choose the time values automatically, as shown in the following sections.

Chordal Parameterization

One way to go about this is to measure the (Euclidean) distances between consecutive vertices (i.e. the chordal lengths) and simply use those distances as time intervals:

distances = np.linalg.norm(np.diff(points4 + points4[:1], axis=0), axis=1)
distances

array([0.5   , 1.8028, 0.1   , 1.9105, 0.2   , 3.    ])

times5 = np.concatenate([[0], np.cumsum(distances)])
times5

array([0.    , 0.5   , 2.3028, 2.4028, 4.3133, 4.5133, 7.5133])

plot_catmull_rom(points4, times5)

This makes the speed along the spline nearly constant, but the distance between the curve and its longer chords can become quite huge.

Centripetal Parameterization

As a variation of the previous method, the square roots of the chordal lengths can be used to define the time intervals [Lee89].

times6 = np.concatenate([[0], np.cumsum(np.sqrt(distances))])
times6

array([0.    , 0.7071, 2.0498, 2.366 , 3.7482, 4.1954, 5.9275])

plot_catmull_rom(points4, times6)

The curve takes its course much closer to the chords, but its speed is obviously far from constant.

Centripetal parameterization has the very nice property that it guarantees no cusps and no self-intersections, as shown by Yuksel et al. [YSK11]. The curve is also guaranteed to never “move away” from the successive vertex:

When centripetal parameterization is used with Catmull–Rom splines to define a path curve, the direction of motion for the object following this path will always be towards the next key-frame position.

—Yuksel et al. [YSK11], Section 7.2: “Path Curves”

Parameterized Parameterization

It turns out that the previous two parameterization schemes are just two special cases of a more general scheme for obtaining time intervals between control points:

\begin{equation*} t_{i+1} = t_i + |\boldsymbol{x}_{i+1} - \boldsymbol{x}_i|^\alpha, \text{ with } 0 \le \alpha \le 1. \end{equation*}

In the Python class splines.CatmullRom, the parameter alpha can be specified.

def plot_alpha(alpha, label):
    s = splines.CatmullRom(points4, alpha=alpha, endconditions='closed')
    plot_spline_2d(s, label=label)

plot_alpha(0, r'$\alpha = 0$ (uniform)')
plot_alpha(0.5, r'$\alpha = 0.5$ (centripetal)')
plot_alpha(0.75, r'$\alpha = 0.75$')
plot_alpha(1, r'$\alpha = 1$ (chordal)')
plt.legend(loc='center', numpoints=3);

As can be seen here – and as Yuksel et al. [YSK11] demonstrate to be generally true – the uniform curve is farthest away from short chords and closest to long chords. The chordal curve behaves contrarily: closest to short chords and awkwardly far from long chords. The centripetal curve is closer to the uniform curve for long chords and closer to the chordal curve for short chords, providing a very good compromise.

Any value between \(0\) and \(1\) can be chosen for \(\alpha\), but \(\alpha = \frac{1}{2}\) (i.e. centripetal parameterization) stands out because it is the only one of them that guarantees no cusps and self-intersections:

In this paper we prove that, for cubic Catmull–Rom curves, centripetal parameterization is the only parameterization in this family that guarantees that the curves do not form cusps or self-intersections within curve segments.

—Yuksel et al. [YSK11], abstract

[…] we mathematically prove that centripetal parameterization of Catmull–Rom curves guarantees that the curve segments cannot form cusps or local self-intersections, while such undesired features can be formed with all other possible parameterizations within this class.

—Yuksel et al. [YSK11], Section 1: “Introduction”

Cusps and self-intersections are very common with Catmull–Rom curves for most parameterization choices. In fact, as we will show here, the only parameterization choice that guarantees no cusps and self-intersections within curve segments is centripetal parameterization.

—Yuksel et al. [YSK11], Section 3: “Cusps and Self-Intersections”