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

Uniform Catmull–Rom Splines#

Catmull and Rom [CR74] presented a class of splines which can be described mathematically, in its most generic form, with what is referred to as equation (1):

\begin{equation*} F(s) = \frac{ \sum x_i(s) w_i(s) }{ \sum w_i(s) }, \end{equation*} where the part \(w_i(s)/\sum w_i(s)\) is called blending functions.

Since the blending functions presented above are, as of now, completely arbitrary we impose some constraints in order to make them easier to use. We shall deal only with blending functions that are zero outside of some given interval. Also we require that \(\sum w_i(s)\) does not vanish for any \(s\). We shall normalize \(w_i(s)\) so that \(\sum w_i(s) = 1\) for all \(s\).

Catmull and Rom [CR74], section 3, “Blending Functions”

The components of the equation are further constrained to produce an interpolating function:

Consider the following case: Let \(x_i(s)\) be any function interpolating the points \(p_i\) through \(p_{i+k}\) and let \(w_i(s)\) be zero outside \((s_{i-1}, s_{i+k+1})\). The function \(F(s)\) defined in equation (1) will thus be an interpolating function. Intuitively, this says that if all of the functions that have an effect at a point, pass through the point, then the average of the functions will pass through the point.

Catmull and Rom [CR74], section 2: “The Model”

Typo Alert

The typo “\(p_i\) through \(s_{i+k}\)” has been fixed in the quote above.

A polynomial of degree \(k\) that pass[e]s through \(k+1\) points will be used as \(x(s)\). In general it will not pass through the other points. If the width of the interval in which \(w_i(s)\) is non zero is less than or equal to \(k+2\) then \(x_i(s)\) will not affect \(F(s)\) outside the interpolation interval. This means that \(F(s)\) will be an interpolating function. On the other hand if the width of \(w_i(s)\) is greater than \(k+2\) then \(x_i(s)\) will have an effect on the curve outside the interpolation interval. \(F(s)\) will then be an approximating function.

Catmull and Rom [CR74], section 2: “The Model”

After limiting the scope of the paper to interpolating splines, it is further reduced to uniform splines:

[…] in the parametric space we can, without loss of generality, place \(s_j=j\).

Catmull and Rom [CR74], section 2: “The Model”

Whether or not generality is lost, this means that the rest of the paper doesn’t give any hints on how to construct non-uniform splines. For those who are interested nevertheless, we show how to do that in the notebook about non-uniform Catmull–Rom splines and once again in the notebook about the Barry–Goldman algorithm.

After the aforementioned constraints and the definition of the term cardinal function

Cardinal function: a function that is \(1\) at some knot, \(0\) at all other knots and can be anything in between the other knots. It satisfies \(F_i(s_j) = \delta_{ij}\).

Catmull and Rom [CR74], section 1: “Introduction”

… the gratuitously generic equation (1) is made a bit more concrete:

If in equation (1) we assume \(x_i(s)\) to be polynomials of degree \(k\) then this equation can be reduced to a much simpler form:

\begin{equation*} F(s) = \sum_j p_j C_{jk}(s) \end{equation*} where the \(C_{jk}(s)\) are cardinal blending functions and \(j\) is the knot to which the cardinal function and the point belong and each \(C_{jk}(s)\) is a shifted version of \(C_{0,k}(s)\). \(C_{0,k}(s)\) is a function of both the degree \(k\) of the polynomials and the blending functions \(w(s)\):

\begin{equation*} C_{0,k}(s) = \sum_{i=0}^k \Big[ \prod_{\substack{j=i-k\\j \ne 0}}^i \left(\frac{s}{j}+1\right) \Big] w(s+i) \end{equation*}

In essence we see that for a polynomial case our cardinal functions are a blend of Lagrange polynomials. When calculating \(C_{0,k}(s)\), \(w(s)\) should be centered about \(\frac{k}{2}\).

Catmull and Rom [CR74], section 4: “Calculating Cardinal Functions”

This looks like something we can work with, even though the blending function \(w(s)\) is still not defined.

[1]:
import sympy as sp

We use \(t\) instead of \(s\):

[2]:
t = sp.symbols('t')
[3]:
i, j, k = sp.symbols('i j k', integer=True)
[4]:
w = sp.Function('w')
[5]:
C0k = sp.Sum(
        sp.Product(
            sp.Piecewise((1, sp.Eq(j, 0)), ((t / j) + 1, True)),
            (j, i - k, i)) * w(t + i),
        (i, 0, k))
C0k
[5]:
$\displaystyle \sum_{i=0}^{k} w{\left(i + t \right)} \prod_{j=i - k}^{i} \begin{cases} 1 & \text{for}\: j = 0 \\1 + \frac{t}{j} & \text{otherwise} \end{cases}$

Blending Functions#

Catmull and Rom [CR74] leave the choice of blending function to the reader. They show two plots (figure 1 and figure 3) for a custom blending function stitched together from two Bézier curves, but they don’t show the cardinal function nor an actual spline created from it.

The only other concrete suggestion is to use B-spline basis functions as blending functions. A quadratic B-spline basis function is shown in figure 2 and both cardinal functions and example curves are shown that utilize both quadratic and cubic B-spline basis functions (figures 4 through 7). No mathematical description of B-spline basis functions is given, instead they refer to Gordon and Riesenfeld [GR74]. That paper provides a pair of equations (3.1 and 3.2) that can be used to recursively construct B-spline basis functions. Simplified to the uniform case, this leads to the base case (i.e. degree zero) …

[6]:
B0 = sp.Piecewise((0, t < i), (1, t < i + 1), (0, True))
B0
[6]:
$\displaystyle \begin{cases} 0 & \text{for}\: i > t \\1 & \text{for}\: t < i + 1 \\0 & \text{otherwise} \end{cases}$

… which can be used to obtain the linear (i.e. degree one) basis functions:

[7]:
B1 = (t - i) * B0 + (i + 2 - t) * B0.subs(i, i + 1)

We can use one of them (where \(i = 0\)) as blending function:

[8]:
w1 = B1.subs(i, 0)

With some helper functions from helper.py we can plot this.

[9]:
from helper import plot_sympy, grid_lines
[10]:
plot_sympy(w1, (t, -0.2, 2.2))
grid_lines([0, 1, 2], [0, 1])
../_images/euclidean_catmull-rom-uniform_30_0.svg

The quadratic (i.e. degree two) basis functions can be obtained like this:

[11]:
B2 = (t - i) / 2 * B1 + (i + 3 - t) / 2 * B1.subs(i, i + 1)

For our further calculations, we use the function with \(i=-1\) as blending function:

[12]:
w2 = B2.subs(i, -1)
[13]:
plot_sympy(w2, (t, -1.2, 2.2))
grid_lines([-1, 0, 1, 2], [0, 1])
../_images/euclidean_catmull-rom-uniform_35_0.svg

This should be the same function as shown by Catmull and Rom [CR74] in figure 2.

Cardinal Functions#

The first example curve in the paper (figure 5) is a cubic curve, constructed using a cardinal function with \(k=1\) (i.e. using linear Lagrange interpolation) and a quadratic B-spline basis function (as shown above) as blending function.

With the information so far, we can construct the cardinal function \(C_{0,1}(t)\), using our quadratic B-spline blending function w2 (which is, as required, centered about \(\frac{k}{2}\)):

[14]:
C01 = C0k.subs(k, 1).replace(w, lambda x: w2.subs(t, x)).doit().simplify()
C01
[14]:
$\displaystyle \begin{cases} 0 & \text{for}\: t < -2 \\\frac{\left(t + 1\right) \left(t + 2\right)^{2}}{2} & \text{for}\: t < -1 \\- \frac{3 t^{3}}{2} - \frac{5 t^{2}}{2} + 1 & \text{for}\: t < 0 \\\frac{3 t^{3}}{2} - \frac{5 t^{2}}{2} + 1 & \text{for}\: t < 1 \\\frac{\left(1 - t\right) \left(t - 2\right)^{2}}{2} & \text{for}\: t < 2 \\0 & \text{otherwise} \end{cases}$
[15]:
plot_sympy(C01, (t, -2.2, 2.2))
grid_lines(range(-2, 3), [0, 1])
../_images/euclidean_catmull-rom-uniform_40_0.svg

This should be the same function as shown by Catmull and Rom [CR74] in figure 4.

The paper does not show that, but we can also try to flip the respective degrees of Lagrange interpolation and B-spline blending. In other words, we can set \(k=2\) to construct the cardinal function \(C_{0,2}(t)\), this time using the linear B-spline blending function w1 (which is also centered about \(\frac{k}{2}\)) leading to a total degree of 3:

[16]:
C02 = C0k.subs(k, 2).replace(w, lambda x: w1.subs(t, x)).doit().simplify()

And as it turns out, this is exactly the same thing!

[17]:
assert C01 == C02

By the way, we come to the same conclusion in our notebook about the Barry–Goldman algorithm, which means that this is also true in the non-uniform case.

Many authors nowadays, when using the term Catmull–Rom spline, mean the cubic spline created using exactly this cardinal function.

As we have seen, this can be equivalently understood either as three linear interpolations (more exactly: one interpolation and two extrapolations) followed by quadratic B-spline blending or as two overlapping quadratic Lagrange interpolations followed by linear blending. The two equivalent approaches are illustrated by means of animations in the notebook about non-uniform Catmull–Rom splines.

Example Plot#

[18]:
import matplotlib.pyplot as plt
import numpy as np

To quickly check how a spline segment would look like when using the cardinal function we just derived, let’s define a few points …

[19]:
vertices = np.array([
    (-0.1, -0.5),
    (0, 0),
    (1, 0),
    (0.5, 1),
])

… and plot \(F(t)\) (or \(F(s)\), as it has been called originally):

[20]:
plt.scatter(*np.array([
    sum([vertices[i] * C01.subs(t, s - i + 1) for i in range(4)])
    for s in np.linspace(0, 1, 20)]).T)
plt.plot(*vertices.T, 'x:g');
../_images/euclidean_catmull-rom-uniform_53_0.svg

For calculating more than one segment, and also for creating non-uniform Catmull–Rom splines, the class splines.CatmullRom can be used. For more plots, see the notebook about properties of Catmull–Rom splines.

Basis Polynomials#

The piecewise expression for the cardinal function is a bit unwieldy to work with, so let’s bring it into a form we already know how to deal with.

We are splitting the piecewise expression into four separate pieces, each one to be evaluated at \(0 \le t \le 1\). We are also reversing the order of the pieces, to match our intended control point order:

[21]:
b_CR = sp.Matrix([
    expr.subs(t, t + cond.args[1] - 1)
    for expr, cond in C01.args[1:-1][::-1]]).T
b_CR.T
[21]:
$\displaystyle \left[\begin{matrix}- \frac{t \left(t - 1\right)^{2}}{2}\\\frac{3 t^{3}}{2} - \frac{5 t^{2}}{2} + 1\\- \frac{3 \left(t - 1\right)^{3}}{2} - \frac{5 \left(t - 1\right)^{2}}{2} + 1\\\frac{t^{2} \left(t - 1\right)}{2}\end{matrix}\right]$
[22]:
from helper import plot_basis
[23]:
plot_basis(*b_CR, labels=sp.symbols('xbm_i-1 xbm_i xbm_i+1 xbm_i+2'))
../_images/euclidean_catmull-rom-uniform_58_0.svg

For the following sections, we are using a few tools from utility.py:

[24]:
from utility import NamedExpression, NamedMatrix

Basis Matrix#

[25]:
b_monomial = sp.Matrix([t**3, t**2, t, 1]).T
M_CR = NamedMatrix(r'{M_\text{CR}}', 4, 4)
control_points = sp.Matrix(sp.symbols('xbm3:7'))

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

[26]:
p4 = NamedExpression('pbm4', b_monomial * M_CR.name * control_points)
p4
[26]:
$\displaystyle \boldsymbol{p}_{4} = \left[\begin{matrix}t^{3} & t^{2} & t & 1\end{matrix}\right] {M_\text{CR}} \left[\begin{matrix}\boldsymbol{x}_{3}\\\boldsymbol{x}_{4}\\\boldsymbol{x}_{5}\\\boldsymbol{x}_{6}\end{matrix}\right]$

From the basis polynomials and the control points, we can already calculate \(\boldsymbol{p}_4(t)\)

[27]:
p4.expr = b_CR.dot(control_points).expand().collect(t)
p4
[27]:
$\displaystyle \boldsymbol{p}_{4} = t^{3} \left(- \frac{\boldsymbol{x}_{3}}{2} + \frac{3 \boldsymbol{x}_{4}}{2} - \frac{3 \boldsymbol{x}_{5}}{2} + \frac{\boldsymbol{x}_{6}}{2}\right) + t^{2} \left(\boldsymbol{x}_{3} - \frac{5 \boldsymbol{x}_{4}}{2} + 2 \boldsymbol{x}_{5} - \frac{\boldsymbol{x}_{6}}{2}\right) + t \left(- \frac{\boldsymbol{x}_{3}}{2} + \frac{\boldsymbol{x}_{5}}{2}\right) + \boldsymbol{x}_{4}$

… and with a little bit of squinting, we can directly read off the coefficients of the basis matrix:

[28]:
M_CR.expr = sp.Matrix([
    [b.get(m, 0) for b in [
        p4.expr.expand().coeff(cv).collect(t, evaluate=False)
        for cv in control_points]]
    for m in b_monomial])
M_CR.pull_out(sp.S.Half)
[28]:
$\displaystyle {M_\text{CR}} = \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]$

Catmull and Rom [CR74] show this matrix in section 6.

In case you want to copy&paste it, here’s a plain text version:

[29]:
print(_.expr)
(1/2)*Matrix([
[-1,  3, -3,  1],
[ 2, -5,  4, -1],
[-1,  0,  1,  0],
[ 0,  2,  0,  0]])

And, in case somebody needs it, its inverse looks like this:

[30]:
M_CR.I
[30]:
$\displaystyle {M_\text{CR}}^{-1} = \left[\begin{matrix}1 & 1 & -1 & 1\\0 & 0 & 0 & 1\\1 & 1 & 1 & 1\\6 & 4 & 2 & 1\end{matrix}\right]$
[31]:
print(_.expr)
Matrix([[1, 1, -1, 1], [0, 0, 0, 1], [1, 1, 1, 1], [6, 4, 2, 1]])

Tangent Vectors#

To get the tangent vectors, we simply have to take the first derivative …

[32]:
pd4 = p4.diff(t)

… and evaluate it at the beginning and the end of the segment:

[33]:
start_tangent = pd4.evaluated_at(t, 0)
start_tangent
[33]:
$\displaystyle \left.{\frac{d}{d t} \boldsymbol{p}_{4}}\right\rvert_{t=0} = - \frac{\boldsymbol{x}_{3}}{2} + \frac{\boldsymbol{x}_{5}}{2}$
[34]:
end_tangent = pd4.evaluated_at(t, 1)
end_tangent
[34]:
$\displaystyle \left.{\frac{d}{d t} \boldsymbol{p}_{4}}\right\rvert_{t=1} = - \frac{\boldsymbol{x}_{4}}{2} + \frac{\boldsymbol{x}_{6}}{2}$

These two expressions can be generalized to – as already shown in the notebook about Catmull–Rom properties:

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

Using Bézier Segments#

The above equation for the tangent vectors can be used to construct Hermite splines or, after dividing them by 3, to obtain the control points for cubic Bézier spline segments:

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

[35]:
x4, x5 = control_points[1:3]
[36]:
x4tilde = x4 + start_tangent.expr / 3
x4tilde
[36]:
$\displaystyle - \frac{\boldsymbol{x}_{3}}{6} + \boldsymbol{x}_{4} + \frac{\boldsymbol{x}_{5}}{6}$
[37]:
x5tilde = x5 - end_tangent.expr / 3
x5tilde
[37]:
$\displaystyle \frac{\boldsymbol{x}_{4}}{6} + \boldsymbol{x}_{5} - \frac{\boldsymbol{x}_{6}}{6}$

Using Quadrangle Interpolation#

Remember the notebook about quadrangle interpolation? It showed us how to calculate the quadrangle points given the Bézier control points:

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

[38]:
x4bar = 3 * x4tilde / 2 - x5 / 2
x4bar
[38]:
$\displaystyle - \frac{\boldsymbol{x}_{3}}{4} + \frac{3 \boldsymbol{x}_{4}}{2} - \frac{\boldsymbol{x}_{5}}{4}$
[39]:
x5bar = 3 * x5tilde / 2 - x4 / 2
x5bar
[39]:
$\displaystyle - \frac{\boldsymbol{x}_{4}}{4} + \frac{3 \boldsymbol{x}_{5}}{2} - \frac{\boldsymbol{x}_{6}}{4}$

Generalizing these expressions and juggling the terms around a bit, we get

\begin{equation*} \boldsymbol{\bar{x}}_i^{(+)} = \boldsymbol{\bar{x}}_i^{(-)} = \boldsymbol{x}_i - \frac{ (\boldsymbol{x}_{i+1} - \boldsymbol{x}_i) + (\boldsymbol{x}_{i-1} - \boldsymbol{x}_i) }{4}. \end{equation*}