Blossoms and Bézier curves

Table of contents

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:

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.

· 5 min read
../