This notebook is part of https://github.com/AudioSceneDescriptionFormat/splines, see also https://splines.readthedocs.io/.

[back to overview](bezier.ipynb)

# De Casteljau's Algorithm

There are several ways that lead to Bézier curves, one was already shown in
[the notebook about Hermite curves](hermite-uniform.ipynb#Relation-to-Bézier-Splines)
(but only for cubic curves).
In this notebook,
we will derive Bézier curves of arbitrary polynomial degree utilizing
[De Casteljau's algorithm](https://en.wikipedia.org/wiki/De_Casteljau's_algorithm).

### Preparations

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


```python
%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](utility.py)
and
[helper.py](helper.py).


```python
from utility import NamedExpression, NamedMatrix
from helper import plot_basis
```

Let's prepare a few symbols for later use ...


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

... and a helper function for plotting:


```python
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.plot(
        *np.asarray(points).T,
        color='lightgrey',
        linestyle=':',
        marker='x',
        markeredgecolor='black',
    )
    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](casteljau.py):


```python
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 (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 it like this:

\begin{equation*}
\boldsymbol{p}_{0,1}(t) = (1 - t) \boldsymbol{x}_0 + t \boldsymbol{x}_1,
\end{equation*}
where in both cases $t \in [0, 1]$.
Since we will be needing quite a few of those linear interpolations,
let's create a helper function:


```python
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:


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


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

Doesn't look like much, but those are the
[Bernstein bases](https://en.wikipedia.org/wiki/Bernstein_polynomial)
for degree 1.
It doesn't get much more interesting if we plot them:


```python
plot_basis(*b1)
```

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


```python
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
```

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$):


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


```python
_.doit()
```

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


```python
M_B1.I
```

Anywho, let's calculate points on the curve by using the Bernstein basis functions:


```python
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
```


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


```python
plot_curve(linear, points)
```


```python
show_casteljau_animation(points)
```

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


```python
p01
```

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


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

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


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

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


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


```python
plot_basis(*b2)
```


```python
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
```


```python
M_B2.I
```


```python
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
```


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


```python
plot_curve(quadratic, points)
```


```python
show_casteljau_animation(points)
```

#### Quadratic Tangent Vectors

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


```python
v02 = p02.diff(t)
```

... at the beginning and the end of the curve:


```python
v02.evaluated_at(t, 0)
```


```python
v02.evaluated_at(t, 1)
```

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


```python
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$.

Combination of $\boldsymbol{x}_2$ and $\boldsymbol{x}_3$:


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

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


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

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


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

This leads to the cubic Bernstein bases:


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

Those are of course the same Bernstein bases
as we found in
[the notebook about Hermite splines](hermite-uniform.ipynb#Relation-to-Bézier-Splines).


```python
plot_basis(*b3)
```


```python
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
```


```python
M_B3.I
```


```python
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
```


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


```python
plot_curve(cubic, points)
```


```python
show_casteljau_animation(points)
```

#### Cubic Tangent Vectors

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


```python
v03 = p03.diff(t)
```

... at the beginning and the end of the curve:


```python
v03.evaluated_at(t, 0)
```


```python
v03.evaluated_at(t, 1)
```

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](hermite-uniform.ipynb#Relation-to-Bézier-Splines)
and to
[construct uniform Catmull--Rom splines using Bézier segments](catmull-rom-uniform.ipynb#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]$):


```python
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:


```python
p03.evaluated_at(t, 0)
```


```python
p03.evaluated_at(t, 1)
```

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:


```python
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
```

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


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

Of course, those are the same matrices as shown in the
[notebook about uniform cubic Hermite splines](hermite-uniform.ipynb#Relation-to-Bézier-Splines).

### Degree 4 (Quartic)

By now you know the drill,
we 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:


```python
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))
p04
```

Kinda long, but anyway, let's try to extract the Bernstein bases:


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


```python
plot_basis(*b4)
```


```python
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
```


```python
M_B4.I
```


```python
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
```


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


```python
plot_curve(quartic, points)
```


```python
show_casteljau_animation(points)
```

#### Quartic Tangent Vectors

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


```python
v04 = p04.diff(t)
```

... at the beginning and the end of the curve:


```python
v04.evaluated_at(t, 0)
```


```python
v04.evaluated_at(t, 1)
```

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


```python
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](https://en.wikipedia.org/wiki/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](../python-module/splines.rst#splines.Bernstein).


```python
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),
])
```

60cac34219f4cf0b29888630e43a100ff6a399b2
ipynb
Jupyter Notebook
doc/euclidean/bezier-de-casteljau.ipynb
AudioSceneDescriptionFormat/splines
0faf87ffddf3be01087959f80d9c43a4da2ae862
MIT
2021-06-21T11:41:02.000Z
2022-02-15T19:18:47.000Z
doc/euclidean/bezier-de-casteljau.ipynb
AudioSceneDescriptionFormat/splines
0faf87ffddf3be01087959f80d9c43a4da2ae862
MIT
2018-07-09T17:47:29.000Z
2022-01-28T19:55:55.000Z
doc/euclidean/bezier-de-casteljau.ipynb
AudioSceneDescriptionFormat/splines
0faf87ffddf3be01087959f80d9c43a4da2ae862
MIT
2020-05-15T18:34:42.000Z
2021-12-28T06:29:54.000Z
Qwen/Qwen-72B
1. YES
2. YES