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

De Casteljau’s Algorithm

See also https://pomax.github.io/bezierinfo/.

There are several ways to get to Bézier curves, one was already shown in the notebook about Hermite curves (but only for cubic curves).

TODO: first explain control polylines and then link to Hermite splines?

Another one is the so-called De Casteljau’s algorithm. (TODO: link to De Casteljau)

One nice aspect of this is that the algorithm can be used for arbitrary polynomial degrees.

A Bézier spline is defined by a so-called control polyline (or control polygon), which comprises a sequence of control points. Some of those control points are part of the final spline curve, others lie outside of it. The degree of a spline segment determines how many “off-curve” control points are between two “on-curve” control points.

For example, in a cubic (degree = 3) Bézier spline there are two (= degree - 1) “off-curve” control points.

Two equally valid viewpoints for what a Bézier spline is:

• A sequence of curve segments, each defined by degree + 1 control points. The first control point of a segment is the same as the last control point of the previous one.

• A sequence of control points that can be used to shape the resulting curve. Every degree’th control point lies on the curve and the others define the shape of the curve segments.

TODO: most well-known: cubic Bézier splines (show screenshot from drawing program, e.g. Inkscape). The two “off-curve” control points are shown as “handles”.

TODO: typical set of constraints on continuity in drawing programs: C0, C1, G1

Preparations

Before we continue, here are are few preparations for the following calculations:

%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 stuff from the file utility.py:

from utility import NamedExpression, NamedMatrix

helper.py

from helper import plot_basis

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, '.')
    ax.scatter(*np.asarray(points).T, marker='x', c='black')
    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
from IPython.display import display, HTML

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, a.k.a. linear

But 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.

Assume that we have two control points, \(\boldsymbol{x}_0\) and \(\boldsymbol{x}_1\) …

… linear equation …:

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

… in other words … this is called affine combination, but we don’t really have to worry about it …

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

… with \(t \in [0, 1]\) (which is called uniform)

TODO: show change of variables for non-uniform curve?

Since we will be needing quite a bunch of those affine combinations, let’s create a helper function:

def affine_combination(one, two):
    return (1 - t) * one + t * two

Now we can define the equation in SymPy:

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

$\displaystyle \boldsymbol{p}_{0,1} = t \boldsymbol{x}_{1} + \boldsymbol{x}_{0} \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 (https://en.wikipedia.org/wiki/Bernstein_polynomial).

It doesn’t get much more interesting if we plot them:

plot_basis(*b1, labels=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]$

If you ever need that, here’s the inverse:

M_B1.I

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

Anywho, 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, a.k.a. quadratic

Consider three control points, \(\boldsymbol{x}_0\), \(\boldsymbol{x}_1\) and \(\boldsymbol{x}_2\) …

We use the affine combinations of the first two points from above …

p01

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

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

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

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

Finally, we make another affine combination of those two results:

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

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

Bernstein basis functions:

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, labels=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 {M_\text{B}^{(2)}}^{-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

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 that 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]\)).

(v02.expr - 2 * (p12.expr - p01.expr)).simplify()

$\displaystyle 0$

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

Degree 3, a.k.a. cubic

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 affine-combine it with the result for the three points \(\boldsymbol{x}_1\), \(\boldsymbol{x}_2\) and \(\boldsymbol{x}_3\).

Combination of \(\boldsymbol{x}_2\) and \(\boldsymbol{x}_3\):

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

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

Combination of \(\boldsymbol{x}_1\), \(\boldsymbol{x}_2\) and \(\boldsymbol{x}_3\):

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

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

Combination of \(\boldsymbol{x}_0\), \(\boldsymbol{x}_1\), \(\boldsymbol{x}_2\) and \(\boldsymbol{x}_3\):

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

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

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]$

TODO: show that those are the same Bernstein bases as in the notebook about Hermite splines

plot_basis(*b3, labels=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 {M_\text{B}^{(3)}}^{-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

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.

We can now see that the last line 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]\)):

(v03.expr - 3 * (p13.expr - p02.expr)).simplify()

$\displaystyle 0$

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

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 Hermite splines:

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.

TODO: show tangent vectors for non-uniform case

Degree 4, a.k.a. quartic

Consider five control points, \(\boldsymbol{x}_0\), \(\boldsymbol{x}_1\), \(\boldsymbol{x}_2\), \(\boldsymbol{x}_3\) and \(\boldsymbol{x}_4\) …

More combinations!

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

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

Kinda long, but anyway, let’s try 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, labels=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 {M_B^{(4)}}^{-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

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]\)):

(v04.expr - 4 * (p14.expr - p03.expr)).simplify()

$\displaystyle 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\)!

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

with the binomial coefficient \({n \choose i} = \frac{n!}{i!(n - i)!}\).

TODO: link to proof?

TODO: show Bernstein polynomials for “quintic” etc.?

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