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

Barry–Goldman Algorithm

The Barry–Goldman algorithm – named after Barry and Goldman [BG88] – can be used to calculate values of non-uniform Catmull–Rom splines. We have also applied this algorithm to rotation splines.

Catmull and Rom [CR74] describe “a class of local interpolating splines” and Barry and Goldman [BG88] describe “a recursive evaluation algorithm for a class of Catmull–Rom splines”, by which they mean a sub-class of the original class, which only contains splines generated from a combination of Lagrange interpolation and B-spline blending:

In particular, they observed that certain choices led to interpolatory curves. Although Catmull and Rom discussed a more general case, we will restrict our attention to an important class of Catmull–Rom splines obtained by combining B-spline basis functions and Lagrange interpolating polynomials. […] They are piecewise polynomial, have local support, are invariant under affine transformations, and have certain differentiability and interpolatory properties.

—Barry and Goldman [BG88], section 1: “Introduction”

The algorithm can be set up to construct curves of arbitrary degree (given enough vertices and their parameter values), but here we only take a look at the cubic case (using four vertices), which seems to be what most people mean by the term Catmull–Rom splines.

The algorithm is a combination of two sub-algorithms:

The Catmull–Rom evaluation algorithm is constructed by combining the de Boor algorithm for evaluating B-spline curves with Neville’s algorithm for evaluating Lagrange polynomials.

—Barry and Goldman [BG88], abstract

Combining the two will lead to a multi-stage algorithm, where each stage consists of only linear interpolations (and extrapolations).

We will use the algorithm here to derive an expression for the tangent vectors, which will show that the algorithm indeed generates non-uniform Catmull–Rom splines.

Triangular Schemes

Barry and Goldman [BG88] illustrate the presented algorithms using triangular evaluation patterns, which we will use here in a very similar form.

As an example, let’s look at the most basic building block: linear interpolation between two given points (in this case \(\boldsymbol{x}_4\) and \(\boldsymbol{x}_5\) with corresponding parameter values \(t_4\) and \(t_5\), respectively):

\begin{equation*} \begin{array}{ccccc} && \boldsymbol{p}_{4,5} && \\ & \frac{t_5 - t}{t_5 - t_4} && \frac{t - t_4}{t_5 - t_4} & \\ \boldsymbol{x}_4 &&&& \boldsymbol{x}_5 \end{array} \end{equation*}

The values at the base of the triangle are known, and the triangular scheme shows how the value at the apex can be calculated from them.

In this example, to obtain the linear polynomial \(\boldsymbol{p}_{4,5}\) one has to add \(\boldsymbol{x}_4\), weighted by the factor shown next to it (\(\frac{t_5 - t}{t_5 - t_4}\)), and \(\boldsymbol{x}_5\), weighted by the factor next to it (\(\frac{t - t_4}{t_5 - t_4}\)).

The parameter \(t\) can be chosen arbitrarily, but in this example we are mostly interested in the range \(t_4 \le t \le t_5\). If the parameter value is outside this range, the process is more appropriately called extrapolation instead of interpolation. Since we will need linear interpolation (and extrapolation) quite a few times, let’s define a helper function:

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)

Neville’s Algorithm

We have already seen this algorithm in our notebook about Lagrange interpolation, where we have shown the triangular scheme for the cubic case – which is also shown by Barry and Goldman [BG88] in figure 2. In the quadratic case, it looks like this:

\begin{equation*} \begin{array}{ccccccccc} &&&& \boldsymbol{p}_{3,4,5} &&&& \\ &&& \frac{t_5 - t}{t_5 - t_3} && \frac{t - t_3}{t_5 - t_3} &&& \\ && \boldsymbol{p}_{3,4} &&&& \boldsymbol{p}_{4,5} && \\ & \frac{t_4 - t}{t_4 - t_3} && \frac{t - t_3}{t_4 - t_3} & & \frac{t_5 - t}{t_5 - t_4} && \frac{t - t_4}{t_5 - t_4} & \\ \boldsymbol{x}_{3} &&&& \boldsymbol{x}_{4} &&&& \boldsymbol{x}_{5} \end{array} \end{equation*}

import matplotlib.pyplot as plt
import numpy as np

Let’s try to plot this for three points:

points = np.array([
    (0, 0),
    (0.5, 2),
    (3, 0),
])

In the following example plots we show the uniform case (with \(t_3=3\), \(t_4=4\) and \(t_5=5\)), but don’t worry, the algorithm works just as well for arbitrary non-uniform time values.

plot_times = np.linspace(4, 5, 30)

plt.scatter(*np.array([
    lerp(
        [lerp(points[:2], [3, 4], t), lerp(points[1:], [4, 5], t)],
        [3, 5], t)
    for t in plot_times]).T)
plt.plot(*points.T, 'x:g')
plt.axis('equal');

Note that the quadratic curve is defined by three points but we are only evaluating it between two of them (for \(4 \le t \le 5\)).

De Boor’s Algorithm

This algorithm [dB72] can be used to calculate B-spline basis functions.

The quadratic case looks like this:

\begin{equation*} \begin{array}{ccccccccc} &&&& \boldsymbol{p}_{3,4,5} &&&& \\ &&& \frac{t_5 - t}{t_5 - t_4} && \frac{t - t_4}{t_5 - t_4} &&& \\ && \boldsymbol{p}_{3,4} &&&& \boldsymbol{p}_{4,5} && \\ & \frac{t_5 - t}{t_5 - t_3} && \frac{t - t_3}{t_5 - t_3} & & \frac{t_6 - t}{t_6 - t_4} && \frac{t - t_4}{t_6 - t_4} & \\ \boldsymbol{x}_{3} &&&& \boldsymbol{x}_{4} &&&& \boldsymbol{x}_{5} \end{array} \end{equation*}

The cubic case is shown by Barry and Goldman [BG88] in figure 1.

plt.scatter(*np.array([
    lerp(
        [lerp(points[:2], [3, 5], t), lerp(points[1:], [4, 6], t)],
        [4, 5], t)
    for t in plot_times]).T)
plt.plot(*points.T, 'x:g')
plt.axis('equal');

Combining Both Algorithms

Catmull and Rom [CR74] show (in figure 5) an example where linear interpolation is followed by quadratic B-spline blending to create a cubic curve.

We can re-create this example with the building blocks from above:

• At the base of the triangle, we put four known vertices.

• Consecutive pairs of these vertices form three linear interpolations (and extrapolations), resulting in three interpolated (and extrapolated) values.

• On top of these three values, we arrange a quadratic instance of de Boor’s algorithm (as shown above).

This culminates in the final value of the spline (given an appropriate parameter value \(t\)) at the apex of the triangle, which looks like this:

\begin{equation*} \def\negspace{\!\!\!\!\!\!} \begin{array}{ccccccccccccc} &&&&&& \boldsymbol{p}_{3,4,5,6} &&&&&& \\ &&&&& \negspace \frac{t_5 - t}{t_5 - t_4} \negspace && \negspace \frac{t - t_4}{t_5 - t_4} \negspace &&&&& \\ &&&& \boldsymbol{p}_{3,4,5} &&&& \boldsymbol{p}_{4,5,6} &&&& \\ && & \negspace \frac{t_5 - t}{t_5 - t_3} \negspace && \negspace \frac{t - t_3}{t_5 - t_3} \negspace & & \negspace \frac{t_6 - t}{t_6 - t_4} \negspace && \negspace \frac{t - t_4}{t_6 - t_4} \negspace & && \\ && \boldsymbol{p}_{3,4} &&&& \boldsymbol{p}_{4,5} &&&& \boldsymbol{p}_{5,6} && \\ & \negspace \frac{t_4 - t}{t_4 - t_3} \negspace && \negspace \frac{t - t_3}{t_4 - t_3} \negspace & & \negspace \frac{t_5 - t}{t_5 - t_4} \negspace && \negspace \frac{t - t_4}{t_5 - t_4} \negspace & & \negspace \frac{t_6 - t}{t_6 - t_5} \negspace && \negspace \frac{t - t_5}{t_6 - t_5} \negspace & \\ \boldsymbol{x}_3 &&&& \boldsymbol{x}_4 &&&& \boldsymbol{x}_5 &&&& \boldsymbol{x}_6 \end{array} \end{equation*}

Here we are considering the fifth spline segment \(\boldsymbol{p}_{3,4,5,6}(t)\) (represented at the apex of the triangle) from \(\boldsymbol{x}_4\) to \(\boldsymbol{x}_5\) (to be found at the base of the triangle) which corresponds to the parameter range \(t_4 \le t \le t_5\). To calculate the values in this segment, we also need to know the preceding control point \(\boldsymbol{x}_3\) (at the bottom left) and the following control point \(\boldsymbol{x}_6\) (at the bottom right). But not only their positions are relevant, we also need the corresponding parameter values \(t_3\) and \(t_6\), respectively.

This same triangular scheme is also shown by Yuksel et al. [YSK11] in figure 3, except that here we shifted the indices by \(+3\).

Another way to construct a cubic curve with this algorithm would be to swap the degrees of interpolation and blending, in other words:

• Instead of three linear interpolations (and extrapolations), apply two overlapping quadratic Lagrange interpolations using Neville’s algorithm (as shown above) to \(\boldsymbol{x}_3\), \(\boldsymbol{x}_4\), \(\boldsymbol{x}_5\) and \(\boldsymbol{x}_4\), \(\boldsymbol{x}_5\), \(\boldsymbol{x}_6\), respectively. Note that the interpolation of \(\boldsymbol{x}_4\) and \(\boldsymbol{x}_5\) appears in both triangles but has to be calculated only once – see also figures 3 and 4 by Barry and Goldman [BG88].

• This will occupy the lower two stages of the triangle, yielding two interpolated values.

• Those two values are then linearly blended in the final stage.

Readers of the notebook about uniform Catmull–Rom splines may already suspect that, for others it might be a revelation: both ways lead to exactly the same triangular scheme and therefore they are equivalent!

The same scheme, but only for the uniform case, is also shown by Barry and Goldman [BG88] in figure 7, and they casually mention the equivalent cases (with \(m\) being the degree of Lagrange interpolation and \(n\) being the degree of the B-spline basis functions):

Note too from Figure 7 that the case \(n=1\), \(m=2\) […] is identical to the case \(n=2\), \(m=1\) […]

—Barry and Goldman [BG88], section 3: “Examples”

Not an Overhauser Spline

Equally casually, they mention:

Finally, the particular case here is also an Overhauser spline [Ove68].

—Barry and Goldman [BG88], section 3: “Examples”

This is not true. Overhauser splines – as described by Overhauser [Ove68] – don’t provide a choice of parameter values. The parameter values are determined by the Euclidean distances between control points, similar, but not quite identical to chordal parameterization. Calculating a value of a Catmull–Rom spline doesn’t involve calculating any distances.

For completeness’ sake, there are two more combinations that lead to cubic splines, but they have their limitations:

• Cubic Lagrange interpolation, followed by no blending at all, which leads to a cubic spline that’s not \(C^1\) continuous (only \(C^0\)), as shown by Barry and Goldman [BG88] in figure 8.

• No interpolation at all, followed by cubic B-spline blending, which leads to an approximating spline (instead of an interpolating spline), as shown by Barry and Goldman [BG88] in figure 5.

Note

Here we are using the time instances of the Lagrange interpolation also as B-spline knots. Barry and Goldman [BG88] show a more generic formulation of the algorithm with separate parameters \(s_i\) and \(t_i\) in equation (9).

Step by Step

The triangular figure above looks more complicated than it really is. It’s just a bunch of linear interpolations and extrapolations.

Let’s go through the figure above, piece by piece.

import sympy as sp

t = sp.symbols('t')

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

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

We use some custom SymPy-based tools from utility.py:

from utility import NamedExpression, NamedMatrix

First Stage

In the center of the bottom row, there is a straightforward linear interpolation from \(\boldsymbol{x}_4\) to \(\boldsymbol{x}_5\) within the interval from \(t_4\) to \(t_5\).

p45 = NamedExpression('pbm_4,5', lerp([x4, x5], [t4, t5], t))
p45

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

Obviously, this starts at:

p45.evaluated_at(t, t4)

$\displaystyle \left.{\boldsymbol{p}_{4,5}}\right\rvert_{t=t_{4}} = \boldsymbol{x}_{4}$

… and ends at:

p45.evaluated_at(t, t5)

$\displaystyle \left.{\boldsymbol{p}_{4,5}}\right\rvert_{t=t_{5}} = \boldsymbol{x}_{5}$

The bottom left of the triangle looks very similar, with a linear interpolation from \(\boldsymbol{x}_3\) to \(\boldsymbol{x}_4\) within the interval from \(t_3\) to \(t_4\).

p34 = NamedExpression('pbm_3,4', lerp([x3, x4], [t3, t4], t))
p34

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

However, that’s not the parameter range we are interested in! We are interested in the range from \(t_4\) to \(t_5\). Therefore, this is not actually an interpolation between \(\boldsymbol{x}_3\) and \(\boldsymbol{x}_4\), but rather a linear extrapolation starting at \(\boldsymbol{x}_4\) …

p34.evaluated_at(t, t4)

$\displaystyle \left.{\boldsymbol{p}_{3,4}}\right\rvert_{t=t_{4}} = \boldsymbol{x}_{4}$

… and ending at some extrapolated point beyond \(\boldsymbol{x}_4\):

p34.evaluated_at(t, t5)

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

Similarly, at the bottom right of the triangle there isn’t a linear interpolation from \(\boldsymbol{x}_5\) to \(\boldsymbol{x}_6\), but rather a linear extrapolation that just reaches \(\boldsymbol{x}_5\) at the end of the parameter interval (i.e. at \(t=t_5\)).

p56 = NamedExpression('pbm_5,6', lerp([x5, x6], [t5, t6], t))
p56

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

p56.evaluated_at(t, t4)

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

p56.evaluated_at(t, t5)

$\displaystyle \left.{\boldsymbol{p}_{5,6}}\right\rvert_{t=t_{5}} = \boldsymbol{x}_{5}$

Second Stage

The second stage of the algorithm involves linear interpolations of the results of the previous stage.

p345 = NamedExpression('pbm_3,4,5', lerp([p34.name, p45.name], [t3, t5], t))
p345

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

p456 = NamedExpression('pbm_4,5,6', lerp([p45.name, p56.name], [t4, t6], t))
p456

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

Those interpolations are defined over a parameter range from \(t_3\) to \(t_5\) and from \(t_4\) to \(t_6\), respectively. In each case, we are only interested in a sub-range, namely from \(t_4\) to \(t_5\).

These are the start and end points at \(t_4\) and \(t_5\):

p345.evaluated_at(t, t4, symbols=[p34, p45])

$\displaystyle \left.{\boldsymbol{p}_{3,4,5}}\right\rvert_{t=t_{4}} = \frac{\left.{\boldsymbol{p}_{3,4}}\right\rvert_{t=t_{4}} \left(- t_{4} + t_{5}\right) + \left.{\boldsymbol{p}_{4,5}}\right\rvert_{t=t_{4}} \left(- t_{3} + t_{4}\right)}{- t_{3} + t_{5}}$

p345.evaluated_at(t, t5, symbols=[p34, p45])

$\displaystyle \left.{\boldsymbol{p}_{3,4,5}}\right\rvert_{t=t_{5}} = \left.{\boldsymbol{p}_{4,5}}\right\rvert_{t=t_{5}}$

p456.evaluated_at(t, t4, symbols=[p45, p56])

$\displaystyle \left.{\boldsymbol{p}_{4,5,6}}\right\rvert_{t=t_{4}} = \left.{\boldsymbol{p}_{4,5}}\right\rvert_{t=t_{4}}$

p456.evaluated_at(t, t5, symbols=[p45, p56])

$\displaystyle \left.{\boldsymbol{p}_{4,5,6}}\right\rvert_{t=t_{5}} = \frac{\left.{\boldsymbol{p}_{4,5}}\right\rvert_{t=t_{5}} \left(- t_{5} + t_{6}\right) + \left.{\boldsymbol{p}_{5,6}}\right\rvert_{t=t_{5}} \left(- t_{4} + t_{5}\right)}{- t_{4} + t_{6}}$

Third Stage

The last step is quite simple:

p3456 = NamedExpression(
    'pbm_3,4,5,6',
    lerp([p345.name, p456.name], [t4, t5], t))
p3456

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

This time, the interpolation interval is exactly the one we are interested in.

To get the final result, we just have to combine all the above expressions:

p3456 = p3456.subs_symbols(p345, p456, p34, p45, p56).simplify()

This expression is quite unwieldy, so let’s not even look at it.

#p3456

Apart from checking whether it’s really cubic …

sp.degree(p3456.expr, t)

$\displaystyle 3$

… and whether it’s really interpolating …

p3456.evaluated_at(t, t4).simplify()

$\displaystyle \left.{\boldsymbol{p}_{3,4,5,6}}\right\rvert_{t=t_{4}} = \boldsymbol{x}_{4}$

p3456.evaluated_at(t, t5).simplify()

$\displaystyle \left.{\boldsymbol{p}_{3,4,5,6}}\right\rvert_{t=t_{5}} = \boldsymbol{x}_{5}$

… the only thing left to do is to check its …

Tangent Vectors

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

pd3456 = p3456.diff(t)

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

pd3456.evaluated_at(t, t4).simplify().simplify()

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

pd3456.evaluated_at(t, t5).simplify()

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

If all went well, this should be identical to the result in the notebook about non-uniform Catmull–Rom splines. As we have mentioned there, it isn’t even necessary to calculate the last interpolation to get the tangent vectors. At the beginning of the interval (\(t = t_4\)), only the first quadratic polynomial \(\boldsymbol{p}_{3,4,5}(t)\) contributes to the final result, while the other one has a weight of zero. At the end of the interval (\(t = t_5\)), only \(\boldsymbol{p}_{4,5,6}(t)\) is relevant. Therefore, we can simply take their tangent vectors at \(t_4\) and \(t_5\), respectively, and we get the same result:

p345.subs_symbols(p34, p45).diff(t).evaluated_at(t, t4).simplify()

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

p456.subs_symbols(p45, p56).diff(t).evaluated_at(t, t5).simplify()

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

Animation

The linear interpolations (and extrapolations) of this algorithm can be shown graphically. By means of the file barry_goldman.py – and with the help of helper.py – we can show an animation of the algorithm:

from barry_goldman import animation
from helper import show_animation

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

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

show_animation(animation(vertices, times))

Once Loop Reflect

If this doesn’t look very intuitive to you, you are not alone. For a different (and probably more straightforward) point of view, have a look at the notebook about non-uniform Catmull–Rom splines.