Addition-only Bezier cubic splines

It’s not so often you need to code your own bezier curves nowadays. If you do, here’s a fast algorithm which uses fewer arithmetic operations — and additions only. I felt very clever when I discovered this a dozen years ago, but I wasn’t the first.

Following the usual formula, to plot a bezier cubic spline you must perform three multiplications and three additions for each dimension, at every iteration. For a two-dimensional curve, that’s twelve floating-point operations each iteration.

The standard method is as follows (p0..p3 are the control points; t is the progress variable):

void drawBezierByTraditionalMethod(
  const Vector &p0,
  const Vector &p1,
  const Vector &p2,
  const Vector &p3,
  const double tIncrement)
{
    register double t;
    register double ax, ay, bx, by, cx, cy, dx, dy;

    // Find a, b, c, d for x part of equation
    ax = (p1.x - p2.x) / 2 + (p3.x - p0.x) / 6;
    bx = (p2.x + p0.x) / 2 - p1.x - 3 * ax;
    cx = p1.x - p0.x - ax - bx;
    dx = p0.x;
    // Find a, b, c, d for y part of equation
    ay = (p1.y - p2.y) / 2 + (p3.y - p0.y) / 6;
    by = (p2.y + p0.y) / 2 - p1.y - 3 * ay;
    cy = p1.y - p0.y - ay - by;
    dy = p0.y;
    // Draw curve
    for (t = 0.0; t <= 3.0000001; t += tIncrement)
    {
        // Add vertex at point:
        // x = ax*t^3 + bx*t^2 + cx*t + dx
        // y = ay*t^3 + by*t^2 + cy*t + dy
        drawPointToScreenBuffer((int) (((ax * t + bx) * t + cx) * t + dx),
          (int) (((ay * t + by) * t + cy) * t + dy));
    }
}

This works fine, but if you’re keen on drawing lots of curves fast, it can be sped up. The following version uses only three floating-point additions per dimension per iteration — i.e., 6 additions per iteration for a two-dimensional curve (and 6 assignments).

This works through making small increments to running x and y values, but each cycle the increments themselves change by an increment, which itself changes as well each cycle. That’s three ‘levels’ of increments, which are determined from the first, second and third derivatives of the cubic equation.

void drawBezierByFastMethod(
  const Vector &p0,
  const Vector &p1,
  const Vector &p2,
  const Vector &p3,
  const double tIncrement)
{
    register double x, y, incr1x, incr1y, incr2x, incr2y, incr3x, incr3y;
    const double tISquared = tIncrement * tIncrement;
    const double tICubed = tISquared * tIncrement;
    register int count;

    incr3x = (3 * (p1.x - p2.x) + p3.x - p0.x) * tICubed;
    incr3y = (3 * (p1.y - p2.y) + p3.y - p0.y) * tICubed;
    incr2x = (2*p0.x + 4*p2.x - 5*p1.x - p3.x) * tISquared - incr3x;
    incr2y = (2*p0.y + 4*p2.y - 5*p1.y - p3.y) * tISquared - incr3y;
    incr1x = ((18*p1.x - 11*p0.x - 9*p2.x + 2*p3.x) * tIncrement
      - 3 * incr2x - 2 * incr3x) / 6;
    incr1y = ((18*p1.y - 11*p0.y - 9*p2.y + 2*p3.y) * tIncrement
      - 3 * incr2y - 2 * incr3y) / 6;
    x = p0.x;
    y = p0.y;
    const int numSteps = (int) (3 / tIncrement);

    // Draw first vertex
    drawPointToScreenBuffer((int) x, (int) y);
    // Draw rest of curve
    for (count = 1; count <= numSteps; count ++)
    {
        incr2x += incr3x;
        incr2y += incr3y;
        incr1x += incr2x;
        incr1y += incr2y;
        x += incr1x;
        y += incr1y;
        drawPointToScreenBuffer((int) x, (int) y);
    }
}