Blossoms and Bézier curves

An degree $$n$$ Bézier curve is given by the following formula $$C(t) = \sum_{i=0}^n a_i B_i^n(t).$$ 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: $$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$$ 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: $$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$$ 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 $$\frac{dC}{dt} = n a^n_0[t,\dots, t,1],$$ 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: $$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],$$ 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$$: $$\tilde{a}_i^0= a_0^{n+1}[0,0, \dots, 1,1]$$ 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: $$C(t) = c_0 + c_1 t +c_2 t^2$$ for some coefficients $$c_i$$. If we define $$t_1 = t_2 = t$$, we might be able to convince ourselves that $$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$$ 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);
}
}
}