Bézier curve evaluation

From CGAFaq

An efficient and elegant way of evaluating a Bézier curve is using de Casteljau's algorithm. This approach turns out to be a special case of blossoming. de Casteljau's algorithm says that given a Bézier curve defined by a series of control points the curve can be evaluated at a point by the following process:

• Compute new points where , and where .
• With these new points compute points
• Continue this process until only one point remains, .

It can be shown that this point is equal to the point on the curve. (One way to show this is to expand out the linear interpolations and show that the result gives the same weighted sum as the Bézier curve defined in terms of the Bernstein polynomials.) The derivative of the curve is equal to the degree of the Bézier times the difference between the penultimate pair of points, . (Second derivatives and beyond can also be computed in a similar manner along the way; see a reference like Farin--cited below--for details.)

For example, the following code evaluates a cubic Bézier curve with given control points. If the deriv pointer is not NULL, the derivative of the curve is stored there. For the case of control points in three dimensions, this routine would need to be called three times, once for the x, y, and z control point coordinates. (Or it could be modified to handle that case directly.)

#define LERP(a, b, t) ((1.-(t))*(a) + ((t)*(b)))
 
double Bezier1D(double cp[4], double t, double *deriv)
{
    double work[3];

    work[0] = LERP(cp[0], cp[1], t);
    work[1] = LERP(cp[1], cp[2], t);
    work[2] = LERP(cp[2], cp[3], t);

    work[0] = LERP(work[0], work[1], t);
    work[1] = LERP(work[1], work[2], t);

    if (deriv)
        *deriv = 3. * (work[1] - work[0]);

    return LERP(work[0], work[1], t);
}

For a Bézier curve of arbitrary degree, the following routine evaluates it at the given position.

// cp should be of length degree+1
double Bezier1D(double *cp, double t, int degree, double *deriv)
{
    double *work = new double[degree];

    // first one is special to do lerp from cp[] into work[]
    for (int i = 0; i < degree; ++i)
        work[i] = LERP(cp[i], cp[i+1], t);

    // don't do last one in this loop so we can compute derivative if needed
    for (int i = 1; i < degree-1; ++i)
        for (int j = 0; j < degree-i; ++j)
            work[j] = LERP(work[j], work[j+1], t);

    if (deriv)
        *deriv = degree * (work[1] - work[0]);

    double val = LERP(work[0], work[1], t);
    delete[] work;
    return val;
}

Further References

• Gerald Farin, Curves and Surfaces for Computer Aided Geometric Design, Academic Press, 1997, 4th Edition (Chapter 3)