Bezier curves always sounded like something super-advanced, something that only mathematically-minded, skilled professionals who wrote Path Routines for operating systems wrote. After much searching and reading, I found out that one of the reasons Bezier curves are so popular is that they are so easy. Unfortunately, a lot of the articles out there are heavy on the math but light on the implementation.

A search on Bezier curves will quickly inform you that two French guys named De Casteljau and Bezier independently invented the Bezier curve in order to design cars, which have all sorts of curves which are not easily represented by normal functions. You can find plenty of discussion on De Castlejau’s algorithm, along with a nice animation cribbed from Wikipedia, but they inevitably end with “but don’t use this in real life.” And you can also find the equation for Bezier’s algorithm. But pseudocode? Nowhere.

Eventually I realized that nobody gives pseudocode because the the equation is the algorithm:

def bezierCurveAt(u, points):
    p = [0.0, 0.0]  # assuming 2D control points
    # Note that n in the summation is len(points) - 1.
    # Also note that in Python, xrange(0, m) goes from 0, 1, ... m - 1.
    for i in xrange(0, len(points)):
        B = bernstein(i, len(points) - 1, u)
        p[0] += B * points[i][0]
        p[1] += B * points[i][1]
    return p

def calcCurve(controlPoints, nSegments):
    points = []
    for i in xrange(0, nSegments + 1):
        u = float(i) / float(nSegments)
        points.append(bezierCurveAt(u, controlPoints))

So what are these Bernstein coefficients? Well, from a coding standpoint, they are really easy:

def bernstein(i, n, u):
    binomial = float(math.factorial(n)) / (float(math.factorial(i)) * float(math.factorial(n - i)))
    return binomial * math.pow((1.0 - u), n - i) * math.pow(u, i)

The Theory

That is not very satisfying from an understanding standpoint, however. Over the years I have observed that there are a bunch of magic functions that show up in math. Can’t integrate a Gaussian function? No problem, use the error function, problem magically solved! When I studied Physics, we were happily deriving the equations of the spherical model of a hydrogen atom (well, the professor was happily deriving it, I was copying Greek letters as fast as possible, and reading the book after the lecture to figure out what was going on). The equations cannot be solved with the methods of solving differential equations familiar to undergraduate students. No problem, Legendre polynomials to the rescue! We spent a week deriving the hydrogen atom, but nobody bothered explaining Legendre polynomials. The book, which was quite good, said something to the effect of “conveniently, it happens that this equation can be solved with Legendre polynomials.” We, however, are not going to take that easy road.

Bernstein coefficients belong to the Bernstein polynomials, which are named after Sergei Bernstein, who created them (or discovered, depending on your philosophical perspective on mathematics) to prove that any curve can be approximated by polynomials. Basically it is the same idea as a Fourier series, which approximates any periodic (cyclical) function as a sum of sines or cosines. A Fourier series approximates a periodic function as

f(x) ≈ A + B sin(αx) + C sin(2αx) + D sin(3αx) + ... X_n sin(nαx)

(where A, B, etc. are appropriate constants). Bernstein polynomials do the same sort of thing, but for any function, not just periodic ones:

f(x) ≈ Ax⁰ + Bx¹ + Cx² + ... X_n xⁿ

It is only valid over the range [0, 1], but that is not a huge problem, just scale and translate your function appropriately. Approximating a function is pretty much exactly the use case that Bezier was trying to solve—a digital approximation to the smooth curves of cars—so it is not surprising that Bernstein polynomials show up in Bezier curves. Indeed, one would almost expect it.

Bezier Surfaces

Bezier curves turn out to be so easy that we still have time to look at Bezier surfaces, which are basically a 3D mesh, instead of a 2D line. I recently wanted a model of the famous Utah Teapot for computer graphics (the so-called “teapotahedron”), but it turns out it is defined as a series of Bezier surfaces. I was completely unable to find a teapot model that used vertices that I could just send to glDrawArrays. Since I was intimidated by the fancy name “bezier curve,” I was going to compile Mesa and modify glEvalMesh2 to spit out the vertices that glutSolidTeapot generated. Hah! Too much effort. Let’s just write a short Python script.

def bezierSurfaceAt(u, v, points, nUPoints, nVPoints):
    p = [0.0, 0.0, 0.0]  # assuming 3D control points
    for j in xrange(0, nVPoints):
        for i in xrange(0, nUPoints):
            Bi = bernstein(i, nUPoints - 1, u)
            Bj = bernstein(j, nVPoints - 1, v)
            cp = points[j * nUPoints + i]
            p[0] += Bi * Bj * cp[0]
            p[1] += Bi * Bj * cp[1]
            p[2] += Bi * Bj * cp[2]
    return p

This gives us the vertices, but we still need the normals. The normal is simply the tangent in the u direction crossed with the tangent in the v direction. Since the tangent of a curve is simply its derivative,

It turns out that the derivative of a Bernstein polynomial is a sum of lower order Bernstein polynomials:

Thus our normal is:

So we can calculate our normals with a straightforward implementation of the equation:

def bezierNormalAt(u, v, points, nUPoints, nVPoints):
    uTangent = [0.0, 0.0, 0.0]
    for j in xrange(0, nVPoints):
        for i in xrange(0, nUPoints):
            Bi1 = bernstein(i - 1, nUPoints - 2, u)
            Bi2 = bernstein(i, nUPoints - 2, u)
            Bj = bernstein(j, nVPoints - 1, v)
            a = nUPoints * (Bi1 - Bi2) * Bj
            cp = points[j * nUPoints + i]
            uTangent[0] += a * cp[0]
            uTangent[1] += a * cp[1]
            uTangent[2] += a * cp[2]

    vTangent = [0.0, 0.0, 0.0]
    for j in xrange(0, nVPoints):
        for i in xrange(0, nUPoints):
            Bi = bernstein(i, nUPoints - 1, u)
            Bj1 = bernstein(j - 1, nVPoints - 2, v)
            Bj2 = bernstein(j, nVPoints - 2, v)
            a = nUPoints * Bi * (Bj1 - Bj2)
            cp = points[j * nUPoints + i]
            vTangent[0] += a * cp[0]
            vTangent[1] += a * cp[1]
            vTangent[2] += a * cp[2]

    return normalized(crossProduct(uTangent, vTangent))

I have packaged this up in a short Python script (teapot.py) which will generate the vertices, normals, and vertex ids for the Utah teapot. These are suitable for passing directly to glDrawElements() with GL_VERTICES.