This notebook is part of https://github.com/AudioSceneDescriptionFormat/splines, see also https://splines.readthedocs.io/.
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
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
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)
display(HTML(ani.to_jshtml(default_mode='reflect')))
plt.close() # avoid spurious figure display
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
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 for degree 1 (https://en.wikipedia.org/wiki/Bernstein_polynomial).
It doesn't get much more interesting if we plot them:
sp.plot(*b1, (t, 0, 1));
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
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]))
_.doit()
If you ever need that, here's the inverse:
M_B1.I
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)
I know, not very exciting. But it gets better!
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
... and we do the same thing for the second and third point:
p12 = NamedExpression('pbm_1,2', affine_combination(x1, x2))
p12
Finally, we make another affine combination of those two results:
p02 = NamedExpression('pbm_0,2', affine_combination(p01.expr, p12.expr))
p02
Bernstein basis functions:
b2 = [p02.expr.expand().coeff(x.name).factor() for x in (x0, x1, x2)]
b2
sp.plot(*b2, (t, 0, 1));
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
M_B2.I
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)
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)
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]$).
(v02.expr - 2 * (p12.expr - p01.expr)).simplify()
In case you are wondering, the factor 2 comes from the degree 2 of our quadratic curve.
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
Combination of $\boldsymbol{x}_1$, $\boldsymbol{x}_2$ and $\boldsymbol{x}_3$:
p13 = NamedExpression('pbm_1,3', affine_combination(p12.expr, p23.expr))
p13
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
Bernstein bases:
b3 = [p03.expr.expand().coeff(x.name).factor() for x in (x0, x1, x2, x3)]
b3
TODO: show that those are the same Bernstein bases as in the notebook about Hermite splines
sp.plot(*b3, (t, 0, 1));
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
M_B3.I
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)
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)
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.
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()
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)
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 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
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)
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
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
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
sp.plot(*b4, (t, 0, 1));
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
M_B4.I
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)
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)
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]$):
(v04.expr - 4 * (p14.expr - p03.expr)).simplify()
Again, the factor 4 comes from the degree 4 of our curve.
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),
])