Calculating Bezier curves fast via forward differencing

Calculating bezier curves is quite processor intensive task, but this algorithm will make calculating at least 10 times faster. Especially with bezier surfaces of bezier spaces, you will need the speed.

All polynomical functions can be calculated via forward differencing. The key idea is to start at some point of the function, move forwards at constant step and use Taylor series to calculate the next value. With polynomical functions, Nth derivate equals 0. Thus the Taylor series can be cut off there and get rid of the error term. Here's the Tailor series (Rn is the error term).

$f(x) = f(x) + f'(x) (t-x) + f''(x) (t-x)^2 / 2! + ... + f^{(n)}(x) (t-x)^n / n! + Rn$

Now, what we want is to calculate the next point f(x+t), which simplifies the formula a little:

$f(x+t) = f(x) + f'(x) t + f''(x) t^2 / 2! + ... + f^{(n)}(x) t^n / n! + Rn$

How does this help us? For a practical example, lets derive everything step by step for the cubic Bezier formula, for which

Since the function is 3rd order, it's 4th derivate equals 0.

Thus, we can rewrite the Taylor series as:

(1)

It calculates the exact value of f(x+t) (since all the other terms in the Taylor series equal 0). Now all we need to know now is f(x), f '(x), f ''(x) and f '''(x). Since bezier curves are drawn from x=0 to x=1, we should calculate the values of these functions at x=0

Now, with these precalculations and formula (1) above (where x = 0 and t = the point we want), we could easily calculate any point in the bezier curve and we'd only need a couple of multiplications to get the desired point. But lets not stop here, because we can get rid of all the slow multiplications.

We need to have a constant step by which we move through the bezier curve. This is t. Now look at the formula to calculate next step:

To get rid of multiplications, we need a variable for f(x), a variable for f '(x)*t, a variable for f ''(x)*t^2 / 2 and a variable for f '''(x)*t^3 / 6. Let's call these f, fd, fdd_per_2 and fddd_per_6, respectively.

The assignment to evaluate the next point f(x+t) becomes:

f = f + fd + fdd_per_2 + fddd_per_6

Now we can calculate the next point without multiplications in the 'loop'. Hooray! But the problem we now have is that to calculate the next point with this same formula, we need to update fd, fdd_per_2 and fddd_per_6 too, since they're supposed to represent f '(x)*t, f ''(x)*t^2 / 2 and f '''(x)*t^3 / 6, and x just increased by t. Back to Taylor then. First we need to fiqure out how to update fd. Directly from the series follows this:

Now lets just multiply both sides by t

That can be read as an assignment:
fd = fd + ? + ?

Time to introduce new variables.

fd = fd + fdd + fddd_per_2

To calculate fdd's next value, we need:

That can be read as an assignment:
fdd = fdd + fddd

We still need to calculate fdd_per_2. Here's the non-surprising result:

That can be read as an assignment:
fdd_per_2 = fdd_per_2 + fddd_per_2

Now, there's still variables fddd, fddd_per_2 and fddd_per_6 that we haven't updated. Well, the next derivate equals zero from there, so these variables are in fact constants. All done, then.

Here's a C code that wraps it up pretty nicely, and you can see how all the variables are first initialized and then updated in the loop.

//the amount of steps in the bezier curve
unsigned int steps=?;
//the 4 bezier points
double p[4]={?,?,?,?};

//== The algorithm itself begins here ==
double f, fd, fdd, fddd, fdd_per_2, fddd_per_2, fddd_per_6;
double t = 1.0 / double(steps);
double temp = t * t;

//I've tried to optimize the amount of
//multiplications here, but these are exactly
//the same formulas that were derived earlier
//for f(0), f'(0)*t etc.
f = p[0];
fd = 3 * (p[1] - p[0]) * t;
fdd_per_2 = 3 * (p[0] - 2 * p[1] + p[2]) * temp;
fddd_per_2 = 3 * (3 * (p[1] - p[2]) + p[3] - p[0]) * temp * t;

fddd = fddd_per_2 + fddd_per_2;
fdd = fdd_per_2 + fdd_per_2;
fddd_per_6 = fddd_per_2 * (1.0 / 3);

for (unsigned int loop=0; loop < steps; loop++) {
  drawBezierpoint(f);

  f = f + fd + fdd_per_2 + fddd_per_6;
  fd = fd + fdd + fddd_per_2;
  fdd = fdd + fddd;
  fdd_per_2 = fdd_per_2 + fddd_per_2;
}

drawBezierpoint(x[3]);

The initialization needs 11 multiplications, but after that you can get all the points you need with very low cost.

The results

To produce a 9x9 bicubic bezier surface, you'll need 143 multiplications and 923 adds. Compared to this Gamasutra article which introduced an algorithm with 1488 muls (!) and 1506 adds. Thanks to that article I could do this comparsion table for 9x9 bicubic bezier surfaces:

Algorithm name	Muls	Adds
This algorithm	143	923
Central Differencing	1488	1506
De Casteljau's algorithm	3672	3672
Direct evaluation of original formula	8586	3978

And note that this new algo produces the exact same results as the original formula, so it is not an approximation.

Calculating Bezier curves fast via forward differencing

The results

Further reading