This page was generated from doc/euclidean/bezier-de-casteljau.ipynb. Interactive online version: .

De Casteljau’s Algorithm

There are several ways that lead to Bézier curves, one (but only for cubic curves) was already shown in the notebook about Hermite curves. In this notebook, we will derive Bézier curves of arbitrary polynomial degree utilizing De Casteljau’s algorithm.

Preparations

%config InlineBackend.print_figure_kwargs = {'bbox_inches': None}
import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
sp.init_printing()

We import a few utilities and helpers from the files utility.py and helper.py.

from utility import NamedExpression, NamedMatrix
from helper import plot_basis, plot_vertices_2d

Let’s prepare a few symbols for later use …

t, x0, x1, x2, x3, x4 = sp.symbols('t, xbm:5')

… and a helper function for plotting:

def plot_curve(func, points, dots=30, ax=None):
    if ax is None:
        ax = plt.gca()
    times = np.linspace(0, 1, dots)
    ax.plot(*func(points, times).T, '.')
    plot_vertices_2d(points)
    ax.set_title(func.__name__ + ' Bézier curve')
    ax.axis('equal')

We also need to prepare for the animations we will see below. This is using code from the file casteljau.py:

from casteljau import create_animation

def show_casteljau_animation(points, frames=30, interval=200):
    ani = create_animation(points, frames=frames, interval=interval)
    display({
        'text/html': ani.to_jshtml(default_mode='reflect'),
        'text/plain': 'Animations can only be shown in HTML output, sorry!',
    }, raw=True)
    plt.close()  # avoid spurious figure display

Degree 1 (Linear)

After all those preparations, let’s start with the trivial case: A Bézier spline of degree 1 is just a piecewise linear curve connecting all the control points. There are no “off-curve” control points that could bend the curve segments.

Assuming that we have two control points, \(\boldsymbol{x}_0\) and \(\boldsymbol{x}_1\), we can set up a linear equation:

\begin{equation*} \boldsymbol{p}_{0,1}(t) = \boldsymbol{x}_0 + t (\boldsymbol{x}_1 - \boldsymbol{x}_0). \end{equation*}

Another way to write the same thing is like this:

\begin{equation*} \boldsymbol{p}_{0,1}(t) = (1 - t) \boldsymbol{x}_0 + t \boldsymbol{x}_1, \end{equation*} where in both cases \(0 \le t \le 1\). These linear interpolations are sometimes also called affine combinations. Since we will be needing quite a few of those linear interpolations, let’s create a helper function:

def lerp(one, two):
    """Linear interpolation.

    The parameter *t* is expected to be between 0 and 1.

    """
    return (1 - t) * one + t * two

Now we can define the equation in SymPy:

p01 = NamedExpression('pbm_0,1', lerp(x0, x1))
p01

$\displaystyle \boldsymbol{p}_{0,1} = t \boldsymbol{x}_{1} + \boldsymbol{x}_{0} \cdot \left(1 - t\right)$

b1 = [p01.expr.expand().coeff(x.name).factor() for x in (x0, x1)]
b1

$\displaystyle \left[ 1 - t, \ t\right]$

Doesn’t look like much, but those are the Bernstein bases for degree 1. It doesn’t get much more interesting if we plot them:

plot_basis(*b1)

If you want to convert this to coefficients for the monomial basis \([t, 1]\) instead of the Bernstein basis functions, you can use this matrix:

M_B1 = NamedMatrix(
    r'{M_\text{B}^{(1)}}',
    sp.Matrix([[c.coeff(x) for x in (x0, x1)]
               for c in p01.expr.as_poly(t).all_coeffs()]))
M_B1

$\displaystyle {M_\text{B}^{(1)}} = \left[\begin{matrix}-1 & 1\\1 & 0\end{matrix}\right]$

Applying this matrix leads to the coefficients of the linear equation mentioned in the beginning of this section (\(\boldsymbol{p}_{0,1}(t) = t (\boldsymbol{x}_1 - \boldsymbol{x}_0) + \boldsymbol{x}_0\)):

sp.MatMul(M_B1.expr, sp.Matrix([x0, x1]))

$\displaystyle \left[\begin{matrix}-1 & 1\\1 & 0\end{matrix}\right] \left[\begin{matrix}\boldsymbol{x}_{0}\\\boldsymbol{x}_{1}\end{matrix}\right]$

_.doit()

$\displaystyle \left[\begin{matrix}- \boldsymbol{x}_{0} + \boldsymbol{x}_{1}\\\boldsymbol{x}_{0}\end{matrix}\right]$

In case you ever need that, here’s the inverse:

M_B1.I

$\displaystyle \left({M_\text{B}^{(1)}}\right)^{-1} = \left[\begin{matrix}0 & 1\\1 & 1\end{matrix}\right]$

Anyhow, let’s calculate points on the curve by using the Bernstein basis functions:

def linear(points, times):
    """Evaluate linear Bézier curve (given by two points) at given times."""
    return np.column_stack(sp.lambdify(t, b1)(times)) @ points

points = [
    (0, 0),
    (1, 0.5),
]

plot_curve(linear, points)

show_casteljau_animation(points)

Once Loop Reflect

I know, not very exciting. But it gets better!

Degree 2 (Quadratic)

Now we consider three control points, \(\boldsymbol{x}_0\), \(\boldsymbol{x}_1\) and \(\boldsymbol{x}_2\). We use the linear interpolation of the first two points from above …

p01

$\displaystyle \boldsymbol{p}_{0,1} = t \boldsymbol{x}_{1} + \boldsymbol{x}_{0} \cdot \left(1 - t\right)$

… and we do the same thing for the second and third point:

p12 = NamedExpression('pbm_1,2', lerp(x1, x2))
p12

$\displaystyle \boldsymbol{p}_{1,2} = t \boldsymbol{x}_{2} + \boldsymbol{x}_{1} \cdot \left(1 - t\right)$

Finally, we make another linear interpolation between those two results:

p02 = NamedExpression('pbm_0,2', lerp(p01.expr, p12.expr))
p02

$\displaystyle \boldsymbol{p}_{0,2} = t \left(t \boldsymbol{x}_{2} + \boldsymbol{x}_{1} \cdot \left(1 - t\right)\right) + \left(1 - t\right) \left(t \boldsymbol{x}_{1} + \boldsymbol{x}_{0} \cdot \left(1 - t\right)\right)$

From this, we can get the Bernstein basis functions of degree 2:

b2 = [p02.expr.expand().coeff(x.name).factor() for x in (x0, x1, x2)]
b2

$\displaystyle \left[ \left(t - 1\right)^{2}, \ - 2 t \left(t - 1\right), \ t^{2}\right]$

plot_basis(*b2)

M_B2 = NamedMatrix(
    r'{M_\text{B}^{(2)}}',
    sp.Matrix([[c.coeff(x) for x in (x0, x1, x2)]
               for c in p02.expr.as_poly(t).all_coeffs()]))
M_B2

$\displaystyle {M_\text{B}^{(2)}} = \left[\begin{matrix}1 & -2 & 1\\-2 & 2 & 0\\1 & 0 & 0\end{matrix}\right]$

M_B2.I

$\displaystyle \left({M_\text{B}^{(2)}}\right)^{-1} = \left[\begin{matrix}0 & 0 & 1\\0 & \frac{1}{2} & 1\\1 & 1 & 1\end{matrix}\right]$

def quadratic(points, times):
    """Evaluate quadratic Bézier curve (given by three points) at given times."""
    return np.column_stack(sp.lambdify(t, b2)(times)) @ points

points = [
    (0, 0),
    (0.2, 0.5),
    (1, -0.3),
]

plot_curve(quadratic, points)

show_casteljau_animation(points)

Once Loop Reflect

Quadratic Tangent Vectors

For some more insight, let’s look at the first derivative of the curve (i.e. the tangent vector) …

v02 = p02.diff(t)

… at the beginning and the end of the curve:

v02.evaluated_at(t, 0)

$\displaystyle \left.{\frac{d}{d t} \boldsymbol{p}_{0,2}}\right\rvert_{t=0} = - 2 \boldsymbol{x}_{0} + 2 \boldsymbol{x}_{1}$

v02.evaluated_at(t, 1)

$\displaystyle \left.{\frac{d}{d t} \boldsymbol{p}_{0,2}}\right\rvert_{t=1} = - 2 \boldsymbol{x}_{1} + 2 \boldsymbol{x}_{2}$

This shows that the tangent vector at the beginning and end of the curve is parallel to the line from \(\boldsymbol{x}_0\) to \(\boldsymbol{x}_1\) and from \(\boldsymbol{x}_1\) to \(\boldsymbol{x}_2\), respectively. The length of the tangent vectors is twice the length of those lines.

You might have already seen this coming, but it turns out that the last line in De Casteljau’s algorithm (\(\boldsymbol{p}_{1,2}(t) - \boldsymbol{p}_{0,1}(t)\) in our case) is exactly half of the tangent vector (at any given \(t \in [0, 1]\)).

assert (v02.expr - 2 * (p12.expr - p01.expr)).simplify() == 0

In case you are wondering, the factor 2 comes from the degree 2 of our quadratic curve.

Degree 3 (Cubic)

Let’s now consider four control points, \(\boldsymbol{x}_0\), \(\boldsymbol{x}_1\), \(\boldsymbol{x}_2\) and \(\boldsymbol{x}_3\).

By now, the pattern should be clear: We take the result from the first three points from above and linearly interpolate it with the result for the three points \(\boldsymbol{x}_1\), \(\boldsymbol{x}_2\) and \(\boldsymbol{x}_3\), which we will derive in the following.

We still need the combination of \(\boldsymbol{x}_2\) and \(\boldsymbol{x}_3\) …

p23 = NamedExpression('pbm_2,3', lerp(x2, x3))
p23

$\displaystyle \boldsymbol{p}_{2,3} = t \boldsymbol{x}_{3} + \boldsymbol{x}_{2} \cdot \left(1 - t\right)$

… which we are using to calculate the combination of \(\boldsymbol{x}_1\), \(\boldsymbol{x}_2\) and \(\boldsymbol{x}_3\) …

p13 = NamedExpression('pbm_1,3', lerp(p12.expr, p23.expr))
p13

$\displaystyle \boldsymbol{p}_{1,3} = t \left(t \boldsymbol{x}_{3} + \boldsymbol{x}_{2} \cdot \left(1 - t\right)\right) + \left(1 - t\right) \left(t \boldsymbol{x}_{2} + \boldsymbol{x}_{1} \cdot \left(1 - t\right)\right)$

… which we need for the combination of \(\boldsymbol{x}_0\), \(\boldsymbol{x}_1\), \(\boldsymbol{x}_2\) and \(\boldsymbol{x}_3\):

p03 = NamedExpression('pbm_0,3', lerp(p02.expr, p13.expr))
p03

$\displaystyle \boldsymbol{p}_{0,3} = t \left(t \left(t \boldsymbol{x}_{3} + \boldsymbol{x}_{2} \cdot \left(1 - t\right)\right) + \left(1 - t\right) \left(t \boldsymbol{x}_{2} + \boldsymbol{x}_{1} \cdot \left(1 - t\right)\right)\right) + \left(1 - t\right) \left(t \left(t \boldsymbol{x}_{2} + \boldsymbol{x}_{1} \cdot \left(1 - t\right)\right) + \left(1 - t\right) \left(t \boldsymbol{x}_{1} + \boldsymbol{x}_{0} \cdot \left(1 - t\right)\right)\right)$

This leads to the cubic Bernstein bases:

b3 = [p03.expr.expand().coeff(x.name).factor() for x in (x0, x1, x2, x3)]
b3

$\displaystyle \left[ - \left(t - 1\right)^{3}, \ 3 t \left(t - 1\right)^{2}, \ - 3 t^{2} \left(t - 1\right), \ t^{3}\right]$

Those are of course the same Bernstein bases as we found in the notebook about Hermite splines.

plot_basis(*b3)

M_B3 = NamedMatrix(
    r'{M_\text{B}^{(3)}}',
    sp.Matrix([[c.coeff(x) for x in (x0, x1, x2, x3)]
               for c in p03.expr.as_poly(t).all_coeffs()]))
M_B3

$\displaystyle {M_\text{B}^{(3)}} = \left[\begin{matrix}-1 & 3 & -3 & 1\\3 & -6 & 3 & 0\\-3 & 3 & 0 & 0\\1 & 0 & 0 & 0\end{matrix}\right]$

M_B3.I

$\displaystyle \left({M_\text{B}^{(3)}}\right)^{-1} = \left[\begin{matrix}0 & 0 & 0 & 1\\0 & 0 & \frac{1}{3} & 1\\0 & \frac{1}{3} & \frac{2}{3} & 1\\1 & 1 & 1 & 1\end{matrix}\right]$

def cubic(points, times):
    """Evaluate cubic Bézier curve (given by four points) at given times."""
    return np.column_stack(sp.lambdify(t, b3)(times)) @ points

points = [
    (0, 0.3),
    (0.2, 0.5),
    (0.1, 0),
    (1, 0.2),
]

plot_curve(cubic, points)

show_casteljau_animation(points)

Once Loop Reflect

Cubic Tangent Vectors

As before, let’s look at the derivative (i.e. the tangent vector) of the curve …

v03 = p03.diff(t)

… at the beginning and the end of the curve:

v03.evaluated_at(t, 0)

$\displaystyle \left.{\frac{d}{d t} \boldsymbol{p}_{0,3}}\right\rvert_{t=0} = - 3 \boldsymbol{x}_{0} + 3 \boldsymbol{x}_{1}$

v03.evaluated_at(t, 1)

$\displaystyle \left.{\frac{d}{d t} \boldsymbol{p}_{0,3}}\right\rvert_{t=1} = - 3 \boldsymbol{x}_{2} + 3 \boldsymbol{x}_{3}$

This shows that the tangent vector at the beginning and end of the curve is parallel to the line from \(\boldsymbol{x}_0\) to \(\boldsymbol{x}_1\) and from \(\boldsymbol{x}_2\) to \(\boldsymbol{x}_3\), respectively. The length of the tangent vectors is three times the length of those lines. This also means that if the begin and end positions \(\boldsymbol{x}_0\) and \(\boldsymbol{x}_3\) as well as the corresponding tangent vectors \(\boldsymbol{\dot{x}}_0\) and \(\boldsymbol{\dot{x}}_3\) are given, it’s easy to calculate the two missing control points:

\begin{align*} \boldsymbol{x}_1 &= \boldsymbol{x}_0 + \frac{\boldsymbol{\dot{x}}_0}{3}\\ \boldsymbol{x}_2 &= \boldsymbol{x}_3 - \frac{\boldsymbol{\dot{x}}_3}{3} \end{align*}

This can be used to turn uniform Hermite splines into Bézier splines and to construct uniform Catmull–Rom splines using Bézier segments.

We can now also see that the last linear segment in De Casteljau’s algorithm (\(\boldsymbol{p}_{1,3}(t) - \boldsymbol{p}_{0,2}(t)\) in this case) is exactly a third of the tangent vector (at any given \(t \in [0, 1]\)):

assert (v03.expr - 3 * (p13.expr - p02.expr)).simplify() == 0

Again, the factor 3 comes from the degree 3 of our curve.

Cubic Bézier to Hermite Segments

We now know the tangent vectors at the beginning and the end of the curve, and obviously we know the values of the curve at the beginning and the end:

p03.evaluated_at(t, 0)

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

p03.evaluated_at(t, 1)

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

With these four pieces of information, we can find a transformation from the four Bézier control points to the two control points and two tangent vectors of a Hermite spline segment:

M_BtoH = NamedMatrix(
    r'{M_\text{B$\to$H}}',
    sp.Matrix([[expr.coeff(cv) for cv in [x0, x1, x2, x3]]
               for expr in [
                   x0,
                   x3,
                   v03.evaluated_at(t, 0).expr,
                   v03.evaluated_at(t, 1).expr]]))
M_BtoH

$\displaystyle {M_\text{B$\to$H}} = \left[\begin{matrix}1 & 0 & 0 & 0\\0 & 0 & 0 & 1\\-3 & 3 & 0 & 0\\0 & 0 & -3 & 3\end{matrix}\right]$

And we can simply invert this if we want to go in the other direction, from Hermite to Bézier:

M_BtoH.I.pull_out(sp.S.One / 3)

$\displaystyle {M_\text{B$\to$H}}^{-1} = \frac{1}{3} \left[\begin{matrix}3 & 0 & 0 & 0\\3 & 0 & 1 & 0\\0 & 3 & 0 & -1\\0 & 3 & 0 & 0\end{matrix}\right]$

Of course, those are the same matrices as shown in the notebook about uniform cubic Hermite splines.

Degree 4 (Quartic)

By now you know the drill, let’s consider five control points, \(\boldsymbol{x}_0\), \(\boldsymbol{x}_1\), \(\boldsymbol{x}_2\), \(\boldsymbol{x}_3\) and \(\boldsymbol{x}_4\), which lead to more linear interpolations:

p34 = NamedExpression('pbm_3,4', lerp(x3, x4))
p24 = NamedExpression('pbm_2,4', lerp(p23.expr, p34.expr))
p14 = NamedExpression('pbm_1,4', lerp(p13.expr, p24.expr))
p04 = NamedExpression('pbm_0,4', lerp(p03.expr, p14.expr))

The resulting expression for \(\boldsymbol{p}_{0,4}(t)\) is quite long and unwieldy (and frankly, quite boring as well), so we are not showing it here.

#p04

Instead, we are using it immediately to extract the Bernstein bases:

b4 = [p04.expr.expand().coeff(x.name).factor() for x in (x0, x1, x2, x3, x4)]
b4

$\displaystyle \left[ \left(t - 1\right)^{4}, \ - 4 t \left(t - 1\right)^{3}, \ 6 t^{2} \left(t - 1\right)^{2}, \ - 4 t^{3} \left(t - 1\right), \ t^{4}\right]$

plot_basis(*b4)

M_B4 = NamedMatrix(
    '{M_B^{(4)}}',
    sp.Matrix([[c.coeff(x) for x in (x0, x1, x2, x3, x4)]
               for c in p04.expr.as_poly(t).all_coeffs()]))
M_B4

$\displaystyle {M_B^{(4)}} = \left[\begin{matrix}1 & -4 & 6 & -4 & 1\\-4 & 12 & -12 & 4 & 0\\6 & -12 & 6 & 0 & 0\\-4 & 4 & 0 & 0 & 0\\1 & 0 & 0 & 0 & 0\end{matrix}\right]$

M_B4.I

$\displaystyle \left({M_B^{(4)}}\right)^{-1} = \left[\begin{matrix}0 & 0 & 0 & 0 & 1\\0 & 0 & 0 & \frac{1}{4} & 1\\0 & 0 & \frac{1}{6} & \frac{1}{2} & 1\\0 & \frac{1}{4} & \frac{1}{2} & \frac{3}{4} & 1\\1 & 1 & 1 & 1 & 1\end{matrix}\right]$

def quartic(points, times):
    """Evaluate quartic Bézier curve (given by five points) at given times."""
    return np.column_stack(sp.lambdify(t, b4)(times)) @ points

points = [
    (0, 0),
    (0.5, 0),
    (0.7, 1),
    (1, 1.5),
    (-1, 1),
]

plot_curve(quartic, points)

show_casteljau_animation(points)

Once Loop Reflect

Quartic Tangent Vectors

For completeness’ sake, let’s look at the derivative (i.e. the tangent vector) of the curve …

v04 = p04.diff(t)

… at the beginning and the end of the curve:

v04.evaluated_at(t, 0)

$\displaystyle \left.{\frac{d}{d t} \boldsymbol{p}_{0,4}}\right\rvert_{t=0} = - 4 \boldsymbol{x}_{0} + 4 \boldsymbol{x}_{1}$

v04.evaluated_at(t, 1)

$\displaystyle \left.{\frac{d}{d t} \boldsymbol{p}_{0,4}}\right\rvert_{t=1} = - 4 \boldsymbol{x}_{3} + 4 \boldsymbol{x}_{4}$

By now it shouldn’t be surprising that the tangent vector at the beginning and end of the curve is parallel to the line from \(\boldsymbol{x}_0\) to \(\boldsymbol{x}_1\) and from \(\boldsymbol{x}_3\) to \(\boldsymbol{x}_4\), respectively. The length of the tangent vectors is four times the length of those lines. The last line in De Casteljau’s algorithm (\(\boldsymbol{p}_{1,4}(t) - \boldsymbol{p}_{0,3}(t)\) in this case) is exactly a fourth of the tangent vector (at any given \(t \in [0, 1]\)):

assert (v04.expr - 4 * (p14.expr - p03.expr)).simplify() == 0

Again, the factor 4 comes from the degree 4 of our curve.

Arbitrary Degree

We could go on doing this for higher and higher degrees, but this would get more and more annoying. Luckily, there is a closed formula available to calculate Bernstein polynomials for an arbitrary degree \(n\) (using the binomial coefficient \({n \choose i} = \frac{n!}{i!(n - i)!}\)):

\begin{equation*} b_{i,n}(x) = {n \choose i} x^i \left( 1 - x \right)^{n - i}, \quad i = 0, \ldots, n. \end{equation*}

This is used in the Python class splines.Bernstein.

show_casteljau_animation([
    (0, 0),
    (-1, 1),
    (-0.5, 2),
    (1, 2.5),
    (2, 2),
    (2, 1.5),
    (0.5, 0.5),
    (1, -0.5),
])

Once Loop Reflect