An degree \(n\) Bézier curve is given by the following formula \begin{equation} C(t) = \sum_{i=0}^n a_i B_i^n(t). \end{equation} The coefficients \(a_i\) are called control points are \(n\)-dimensional vectors to produce an \(n\)-dimensional curve. The function \(B_i^n(t) = \binom{n}{i}t^i(1-t)^{n-i}\) is known as the \(i\)th Bernstein basis function of degree \(n\) (with \(\binom{n}{i} = \frac{n!}{i!(n-i)!})\). As with most polynomials, you can evaluate \(C(t)\) at a given value of \(t\) by evaluating the basis functions at \(t\) and using the above formula.
But Bézier curves have a cool property: you can compute \(C(t)\) directly from the control points without evaluating the basis functions. The algorithm to achieve this is called de Casteljau’s algorithm and can be expressed compactly in the following recursive formula: \begin{equation} a_i^r = (1-t)a^{r-1}_i + ta_{i+1}^{r-1},\quad r = 1,\dots, n, \quad i=0,\dots, n-r \end{equation} with \(a_i^0 = a_i\). Evaluating \(a_0^n\) for a given \(t\) will produce the value \(C(t)\), i.e., \(a_0^n = C(t)\). Each iteration of \(r\) performs a linear interpolation between the control points (or the intermediate control points from the previous recursive evaluation). In pseudocode form, we can write this recursively as:
def de_casteljau(r, a):
if r == 0:
return a[i]
else:
# a^{r-1} and a^r are arrays of coefficients
a^{r-1} = de_casteljau(r-1, a)
for i in range(n):
a^r[i] = (1-t) * a^{r-1}[i] + t * a^{r-1}[i+1]
return a^r
Why are we talking about flowers?
This algorithm becomes much more interesting when we ask: what happens if we vary \(t\) with each recursive step? This means that we now have a vector \(\mathbf{t} = (t_1, t_2, \dots, t_n)\) and our equation becomes: \begin{equation} a_i^r = (1-t_r)a^{r-1}_i + t_ra_{i+1}^{r-1},\quad r = 1,\dots, n, \quad i=0,\dots, n-r \end{equation} again with \(a_i^0 = a_i\). This is what we call the blossom of the control points \(a_i\) above over the values \(\mathbf{t}\), which we’ll write as \(a_i^r[\mathbf{t}]\) or \(a_i^r[t_1,\dots, t_n]\) for the intermediate blossom levels.
To evaluate the full blossom, we write \(a_0^n[\mathbf{t}]\), as with de Casteljau. The pseudocode looks very similar to de Casteljau:
def blossom(r, a):
if r == 0:
return a[i]
else:
a^{r-1} = blossom(r-1, a)
for i in range(r):
a^r[i] = (1-t) * a^{r-1}[i] + t * a^{r-1}[i+1]
return a^r
When \(t_i = t\) for each \(i\), we recover de Casteljau’s algorithm, but we’re now free to vary \(t\) throughout the algorithm.
Nice properties of blossoms
This might seem like a trivial generalization, but blossoms have some interesting uses:
- Blossoms are symmetric: \(a_0^n[\mathbf{t}] = a_0^n[\pi(\mathbf{t})]\), where \(\pi\) is a permutation (i.e., reordering) of the entries of \(\mathbf{t}\).
- Blossoms are multiaffine: \(a_0^n[bt_1 + ct_2, \dots] = ba_0^n[t_1, \dots] + ca_0^n[t_2,\dots]\)
- We can express the \(i\)th control point as \(a_i^0 = a_0^n[\mathbf{\tau}_i]\) with the vector \(\tau_i = (0, 0,\dots, 1,1)\) containing \(n-i\) zeros and \(i\) ones.
- Similarly, we can write down the Bézier form of a subcurve on the domain \([c,d]\) with in terms of blossoms. We can compute the \(i\)th control point of the subcurve by evaluating \(a_0^n[\mathbf{\eta}_i]\) with \(\eta_i= (c,c,\dots, d,d)\) being a vector of \(n-i\) copies of \(c\) and \(i\) copies of \(d\).
- We can differentiate Bézier curves trivially with blossoms, using the following expression \begin{equation} \frac{dC}{dt} = n a^n_0[t,\dots, t,1], \end{equation} evaluting the blossom with \(n-1\) copies of \(t\) followed by a single 1.
- We can elevate the degree of a Bézier curve by summing over various blossoms: \begin{equation} a_0^{n+1}[t_1, \dots, t_{n+1}] = \frac{1}{n+1}\sum_{i=0}^{n+1} a_0^n[t_1,\dots,t_{n+1}\mid t_i], \end{equation} where the notation \(t_1,\dots,t_{n+1}\mid t_i\) means that the entry \(t_i\) is omitted from the sequence. This isn’t the most efficient way to evaluate an elevated degree Bézier curve, but it does lead to a compact formula for the elevated curve’s control points \(\tilde{a}_i^0\): \begin{equation} \tilde{a}_i^0= a_0^{n+1}[0,0, \dots, 1,1] \end{equation} with \(n+1-i\) zeroes and \(i\) ones as arguments in the blossom.
Every polynomial in one variable has a unique form as a blossom, but blossoms are mostly used in the context of Bézier curves. For some intuition to see why this might be true, lets look at the quadratic case. We know that \(a_0^n[t,t] = C(t)\), and since \(C(t)\) is a polynomial, we can write it as: \begin{equation} C(t) = c_0 + c_1 t +c_2 t^2 \end{equation} for some coefficients \(c_i\). If we define \(t_1 = t_2 = t\), we might be able to convince ourselves that \begin{equation} a_0^n[t_1,t_2] = c_0^\prime + c_1^\prime t_1 + c_2^\prime t_2 + c_3^\prime t_1t_2 \end{equation} for some other coefficients \(c^\prime_i\). By equating these two equations, we can see that \(c_0 = c_0^\prime, c_1 = \frac{c_1^\prime + c_2^\prime}{2},\) and \(c_3 = c_3^\prime\). A cubic example is worked out here on page 4.
de Castejau implementation
Here’s a simple C++ implementation of de Casteljau’s algorithm, using Eigen. It’s fairly simple to implement without many surprises. The full implementation can be found in this commit of nanospline.
// _dim and _degree are dimension and degree of the Bezier curve.
using Scalar = double;
using ControlPoints = Eigen::Matrix<Scalar, _dim, _degree+1>
...
ControlPoints de_casteljau(Scalar t, int num_recursions) const {
const auto degree = Base::get_degree();
if (num_recursions < 0 || num_recursions > degree) {
throw invalid_setting_error( "Number of de Casteljau recursion
cannot exceeds degree");
}
if (num_recursions == 0) {
// get original control points at the bottom of the recursion
return Base::m_control_points;
} else {
ControlPoints ctrl_pts = de_casteljau(t, num_recursions-1);
assert(ctrl_pts.rows() >= degree+1-num_recursions);
for (int i=0; i<degree+1-num_recursions; i++) {
// ctrl_pts.row(i) gets the i-th control point.
ctrl_pts.row(i) = (1.0-t) * ctrl_pts.row(i) +
t * ctrl_pts.row(i+1);
}
return ctrl_pts;
}
}
Blossom implementation
For our blossom implementation, we can essentially just add a vector input instead of a single value of \(t\).
For a slightly cleaner implementation, we use two for loops instead of an explicit recursion.
The full implementation is available in the Bezier<...>::evaluate function in nanospline.
using Scalar = double;
using ControlPoints = Eigen::Matrix<Scalar, _dim, _degree+1>
using BlossomVector = Eigen::Matrix<Scalar, _degree, 1>
...
void blossom(const BlossomVector& blossom_vector, int degree,
ControlPoints& control_pts) const {
for (int r = 1; r <= degree; r++) {
for (int j = degree; j >= r; j--) {
Scalar t = blossom_vector(r- 1);
control_pts.row(j) =(1. - t) * control_pts.row(j - 1) +
t * control_pts.row(j);
}
}
}
For more information
Blossoms are explained in the most detail in Gerald Farin’s Curves and Surfaces for Computer-Aided Graphics and Design. For more Bézier algorithms and useful visualizations, the Primer on Bézier curves is a great resource. These lecture notes provide a bit more detail about Bézier and B-Spline blossoms.