 Dr. Dobb's is part of the Informa Tech Division of Informa PLC

This site is operated by a business or businesses owned by Informa PLC and all copyright resides with them. Informa PLC's registered office is 5 Howick Place, London SW1P 1WG. Registered in England and Wales. Number 8860726.

Forward Difference Calculation of Bezier Curves

November 1997/Forward Difference Calculation of Bezier Curves

Bezier curves are used widely in graphics systems because they are intuitive and flexible. They can also be easy to compute.

Cubic Bezier curves are commonly used elements in graphics systems. They're used in Adobe's Postscript, for example, and in many programs that output primarily PostScript, such as Adobe Illustrator and Macromedia FreeHand. Bezier curves in general have some nice properties that make them ideal for constructing vector-based graphics. This article shows how to render cubic Bezier curves using forward differencing. The Rendering problem is a common one, and it is especially useful in illustrating two concepts that play well together: the mathematical definition of the cubic Bezier curve, and the forward differencing technique as applied to calculation of polynomials.

The biggest obstacle I encountered when I first studied these concepts was the mathematical notation. It had been several years since I had done anything really mathematical, and I had become a bit rusty. However, I assure you that this article involves nothing more advanced than basic high school algebra. I've tried to avoid skipping too many steps when explaining the mathematics. I apologize if I've carried the explanation a bit to the extreme. Looking at a new equation for the first time is a lot like looking at someone else's code and trying to figure out what it does. Since I'm trying to explain some new ideas, not challenge the reader to solve a puzzle, I've decided to err on the verbose side.

Properties of Bezier Curves

A Bezier curve is described in terms of a control polygon, whose vertices are commonly called control points. An example of a cubic Bezier curve is shown in Figure 1. A Bezier curve begins at exactly the first control point (K0 in Figure 1) and terminates at exactly the last control point (K3 in Figure 1) . This makes it trivial to construct more complicated curves, since two curves may be connected simply by making the last point of one curve coincident with the first point of another.

Another property of the Bezier curve is that it leaves its first control point on a tangent to the line segment between the first and second control points and enters the final control point on a tangent to the line segment connecting the last and next-to-last control point. A Bezier curve is also completely bounded by the convex hull of its control points. The convex hull is a complex polygon enclosing all of the control points. This bounding property turns out to be quite useful in practice. Furthermore, a Bezier curve is invariant under affine transformations. An affine transformation is a scale, shear, translation, rotation, or any combination thereof. Since a Bezier curve is invariant under such transformations, transforming the control points of a curve has the same effect as transforming the curve itself.

Bezier curves may be of any degree, but in practice cubic bezier curves seem to be the most common. A cubic Bezier curve is described by four points: the first and fourth describe the endpoints of the curve; the second and third control the tangents of the curve at the endpoints, as well as the magnitude of the curvature. I focus on cubic Bezier curves for this article because they are easy to work with, quite common, and not least of all, I'm familiar with them from practice.

Formulas for Cubic Bezier Curves

The cubic Bezier curve in Figure 1 is described by four control points, K0 through K3. Shortly I will give a simple formula for converting the four control points into a pair of parametric equations. The equations will be in this form:

```x(t) = axt3 + bxt2 + cxt + dx
y(t) = ayt3 + byt2 + cyt + dy
```

Note that instead of expressing y directly in terms of x with a single equation of the form y = f(x), the parametric equations above express x and y separately. x and y are expressed as functions of a third parameter t, where t is a "parameter" that varies over the range 0 to 1. Parametric equations are flexible; they can represent curves that cannot be represented by a single function y = f(x).

At this point I also introduce a form of notational shorthand: since I so frequently discuss pairs of numbers — the x and y values at a point, or the components of a vector in the x and y directions — I represent an x,y pair as a single capital letter. Thus, I can represent the parametric equations above as:

```P(t) = ( x(t), y(t) )
= At3 + Bt2 + Ct + D
```

where A = (ax, ay), B = (bx, by), etc.

Now, given the four control points K0 = (x0, y0), K1 = (x1, y1), K2 = (x2, y2), and K3 = (x3, y3), here is the formula for computing the coefficients of the parametric equation P(t):

```A = -K0 + 3K1 + -3K2 + K3
B = 3K0 + -6K1 + 3K2
C = -3K0 + 3K1
D = K0
```

or in simple scalar form

```ax = -K0x + 3K1x - 3K2x + K3x
ay = -K0y + 3K1y - 3K2y + K3y
```

etc.

The simplest way to render a cubic bezier curve is to evaluate P(t) at several different places in 0..1, computing a sequence of points along the curve. A drawing routine can then approximate the curve by drawing line segments between each of the computed points. (Obviously, the more points there are, the closer the approximation is to the true curve.) This is essentially the approach I follow here, but with suitable optimizations.

Brute-Force Computation

Here's a simple example of a rendering routine:

```numSteps = 20;
stepSize = 1.0 / (double)numSteps;
MoveTo((int)dx, (int)dy);
for (i = 1; i < numSteps; i++) {
t = stepSize * (double) i;
x = ax * (t * t * t) + bx * (t * t) + cx * t + dx;
y = ay * (t * t * t) + by * (t * t) + cy * t + dy;
LineTo((int)x, (int)y);
}
```

An obvious way to optimize this is to compute the tn terms as follows:

```x = dx;
y = dy;
term = t;
x += term * cx;
y += term * cy;
term *= t;
x += term * bx;
y += term * by;
term *= t;
x += term * ax;
y += term * ay;
```

This is a common optimization which is known as Horner's rule. However, for the cubic 2-dimensional case, this algorithm requires eight multiplies and six additions per step. Fortunately there's something better. Using an algorithm called forward differencing, each step can be computed with six additions and no multiplies at all.

Introduction to Forward Differencing

Forward differencing is an efficient scheme for evaluating a polynomial at many uniform steps. The idea is simple. Suppose you want to evaluate a cubic polynomial of t at several uniform steps. Derive a difference function that expresses the difference between p(t), the value of the polynomial evaluated at t, and p(t+h), the polynomial evaluated at t + h, where h is the uniform step size. For a cubic polynomial and a constant step size, this difference is a quadratic polynomial of t.

Clearly, you can evaluate a quadratic at every step more efficiently than a cubic polynomial. However, if this scheme will work for a cubic equation, it ought to work for a quadratic as well. So, applying the above technique to a quadratic yields a difference function which is a linear function of t. Finally, what works for the quadratic ought to work for the linear, so applying the technique to a linear function yields a difference function which is a constant.

Now it appears that instead of evaluating a cubic polynomial at every step, we can compute it at the initial value of t and then compute each subsequent step by just doing some additions. This is the forward differencing technique.

This idea is pretty simple, but an illustration in code might be helpful. Listing 1 shows part of a program that evaluates a polynomial at multiple uniform steps in t, again, using the brute-force method. After evaluating the polynomial, the program computes the polynomial's forward differences by subtracting the next value of the polynomial from the current value. The program stores these forward differences in the array col3. Next, the program applies the same process to the forward differences themselves. This produces the "second" forward differences, which are stored in the array col4. Finally, the program computes the forward differences of the second forward differences to produce the third forward differences, which are stored in the array col5. The results are shown below:

```t       p(t)     col. 3   col 4.   col 5.
------     ------     ------   ------   ------
0.0000  4.0000   0.0830  -0.0220   0.0180
0.1000  4.0830   0.0610  -0.0040   0.0180
0.2000  4.1440   0.0570   0.0140   0.0180
0.3000  4.2010   0.0710   0.0320   0.0180
0.4000  4.2720   0.1030   0.0500   0.0180
0.5000  4.3750   0.1530   0.0680   0.0180
0.6000  4.5280   0.2210   0.0860   0.0180
```

Since the polynomial is cubic, the third forward differences are all the same. Note that given the initial value of p(t), col3, col4, and col5, you can calculate subsequent values of the polynomial by doing a series of additions at each step. You need no explicit knowledge of the polynomial. The snippet below shows how it is done:

```p         = points;
firstFD   = col3;
secondFD  = col4;
thirdFD   = col5;
for (i = 0; i < NUM_STEPS - 3; i++){
p         += firstFD;
firstFD   += secondFD;
secondFD  += thirdFD;
printf("%3.4f\n",p);
}
```

Finding Initial Values

You can use the above scheme to compute any polynomial at uniform steps, but it would be nice to have a more direct way of finding the initial forward differences. In fact that's easy enough to do.

Consider the case of a simple line:

```p(t) = at + b
```

which is to be evaluated at a sequence of steps of size h:

```p(t + h) - p(t) = a(t + h) + b - (at + b)
p(t + h) - p(t) = at + ah + b - at - b
p(t + h) - p (t) = ah        (Eq. 1)
```

This tells us that if we know p(t0) for some t0, then we know that p(t0 + h) = p(t0) + ah, where ah is always constant. If we know any point on the line, we can compute points for as many fixed-size steps from that point as we want, by simply adding ah on each iteration.

```p(t + h) - p(t) =
a(t + h)2 + b(t + h) + c - (at2 + bt + c)
p(t + h) - p(t) =
a(t2 + 2ht + h2) + bt + bh + c - at2 - bt - c
p(t + h) - p(t) =
at2 + 2aht + ah2 + bt + bh + c - at2 - bt -c
p(t + h) - p(t) =
2aht + ah2 + bh           (Eq. 2) ```

The difference isn't constant this time but instead is a linear function of t. However, a program can use equation (1) to compute the difference function at each step, so it can evaluate p(t) at many steps using a system of two forward differences. We can do the same thing for the cubic case, but since you've probably got the idea now, I omit the intermediate steps in the derivation of the cubic's difference function.

```P(t + h) - p (t) =
a(t + h)3 + b(t + h)2 + c(t + h) + d -
(a3 + b2 + ct + d)
...
P(t + h) - p (t) =
(3ah)t2 + (3ah2 + 2bh)t + ah3 + bh2 + ch    (Eq. 3)
```

Equation (3) is a quadratic equation to which a program can apply equation (2). The program can evaluate p(t) at many different steps using a system of three forward differences. Now I show how to derive the above formula.

Let P(t) = At2 + Bt2 + Ct + D. Plugging a=A, b=B, c=C, and d=D into Equation (3), gives:

```P(t + h) - P(t) = (3Ah)t2 + (3Ah2 + 2Bh)t + (Ah3 + Bh2 + Ch)
```

We can use this formula to compute an initial first forward difference at t = t0:

```F1(t0) = (3Ah)t02 + (3Ah2 + 2Bh)t0 +
(Ah3 + Bh2 + Ch)    (Eq. 4)
```

Now setting a=(3Ah), b=(3Ah2 + 2Bh), and applying Equation (2) yields the initial "second" forward difference at t = t0:

```F2(t0) = 2(3Ah)ht + (3Ah)h2 + (3Ah2 + 2Bh)h
F2(t0) = 6Ah2t + 6Ah3 + 2Bh2
F2(t0) = 6Ah2t0 + 6Ah3 + 2Bh2        (Eq. 5)
```

Finally, setting a=6Ah2 and applying Equation (1) yields an initial third forward difference at t = t0:

```p(t + h) - p (t) = ah
F3(t0) = 6Ah3         (Eq. 6)
```

Rendering with Forward Differences

Recall that a cubic Bezier curve is represented by four control points: K0, K1, K2, K3, which produce the more familiar power series form: P(t) = At3 + Bt2 + Ct + D according to the following formulas:

```A = -K0 + 3K1 + -3K2 + K3
B = 3K0 + -6K1 + 3K2
C = -3K0 + 3K1
D = K0
```

It is now possible to apply equations (4), (5), and (6) to the above to compute the forward differences. Since Bezier curves are defined over the interval 0..1 we want to use t0 = 0, which simplifies the equations by making all the t terms drop out. This yields:

```P = D            // initial point
F1 =  Ah3 + Bh2 + Ch
F2 =  6Ah3 + 2Bh2
F3 =  6Ah3
```

Coding the above in C is straightforward, keeping in mind that each of A, B, C, and D is a 2-vector rather than a simple scalar. Thus, the code in Listing 2 evaluates two different polynomials simultaneously.

The code in Listing 2 generates a sequence of points along the curve defined by the control points (10.0, 70.0), (50.0, 10.0), (150.0, 10.0), and (200.0, 180.0). Another function may then render an approximation of the curve by drawing lines between adjacent points. Such a rendering appears in Figure 2.

The code in Listing 2 suffers from the deficiency of a one-size-fits-all approach. In practice it's desirable to compute just enough points along the curve to maintain some minimum standard of accuracy to the true curve. The step size needed to achieve this goal will vary from curve to curve. In fact, the best compromise between accuracy and efficiency for many curves might best be accomplished by varying the step size when computing points along a single curve. A large step size is most appropriate for areas of relatively low curvature. A smaller step size is appropriate for areas of high curvature.

Forward differencing would at first glance appear to be a poor choice for such a scheme. However, there is a variant of forward differencing, known as adaptive forward differencing, that works in exactly this fashion. If you're interested, you may want to read the paper by Lien et al, referenced below. o

Acknowledgement

I was introduced to forward differencing by Charles Loop.

References

Bob Wallis. "Tutorial on Forward Differencing," Graphics Gems (AP Professional, 1990).

Donald Hearn and M. Pauline Baker. Computer Graphics, C Version (Prentice-Hall, 1997).

S. Lien, M. Shantz, and V. Pratt. "Adaptive Forward Differencing for Rendering Curves and Surfaces," Computer Graphics (SIGGRAPH), 1987.

Curtis Bartley is a software engineer for Microsoft, where he works in the Internet Multimedia Tools group. He can be reached by email at [email protected]

More Insights To upload an avatar photo, first complete your Disqus profile. | View the list of supported HTML tags you can use to style comments. | Please read our commenting policy.

Featured Reports Featured Whitepapers 