# 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).

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

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.

## Further reading

The algorithm I presented here may not be the best one, but it should be plenty fast for many applications. If you need more, read these (of which I've read none):

- Adaptive forward differencing for rendering curves and surfaces; Sheue-Ling Lien, Michael Shantz and Vaughan Pratt; Proceedings of the 14th annual conference on Computer graphics, 1987, Pages 111 - 118
- Rendering trimmed NURBS with adaptive forward differencing; M. Shantz and Sheue-Ling Chang; Proceedings of the 15th annual conference on Computer graphics, 1988, Pages 189 - 198
- Rendering cubic curves and surfaces with integer adaptive forward differencing; S.-L. Chang and M. S. R. Rocchetti; Conference proceedings on Computer graphics, 1989, Pages 157 - 166