How do I implement a Bézier curve in C++?

I'd like to implement a Bézier curve. I've done this in C# before, but I'm totally unfamiliar with the C++ libraries. How should I go about creating a quadratic curve?

void printQuadCurve(float delta, Vector2f p0, Vector2f p1, Vector2f p2);

Clearly we'd need to use linear interpolation, but does this exist in the standard math library? If not, where can I find it?

I'm using Linux.

Recently I ran across the same question and wanted to implemented it on my own. This image from Wikipedia helped me:

The following code is written in C++ and shows how to compute a quadratic bezier.

int getPt( int n1 , int n2 , float perc )
{
    int diff = n2 - n1;

    return n1 + ( diff * perc );
}    

for( float i = 0 ; i < 1 ; i += 0.01 )
{
    // The Green Line
    xa = getPt( x1 , x2 , i );
    ya = getPt( y1 , y2 , i );
    xb = getPt( x2 , x3 , i );
    yb = getPt( y2 , y3 , i );

    // The Black Dot
    x = getPt( xa , xb , i );
    y = getPt( ya , yb , i );

    drawPixel( x , y , COLOR_RED );
}

With (x1|y1), (x2|y2) and (x3|y3) being P0, P1 and P2 in the image. Just for showing the basic idea...

For the ones who ask for the cubic bezier, it just works analogue (also from Wikipedia):

This answer provides Code for it.

Here is a general implementation for a curve with any number of points.

vec2 getBezierPoint( vec2* points, int numPoints, float t ) {
    vec2* tmp = new vec2[numPoints];
    memcpy(tmp, points, numPoints * sizeof(vec2));
    int i = numPoints - 1;
    while (i > 0) {
        for (int k = 0; k < i; k++)
            tmp[k] = tmp[k] + t * ( tmp[k+1] - tmp[k] );
        i--;
    }
    vec2 answer = tmp[0];
    delete[] tmp;
    return answer;
}

Note that it uses heap memory for a temporary array which is not all that efficient. If you only need to deal with a fixed number of points you could hard-code the numPoints value and use stack memory instead.

Of course, the above assumes you have a vec2 structure and operators for it like this:

struct vec2 {
    float x, y;
    vec2(float x, float y) : x(x), y(y) {}
};

vec2 operator + (vec2 a, vec2 b) {
    return vec2(a.x + b.x, a.y + b.y);
}

vec2 operator - (vec2 a, vec2 b) {
    return vec2(a.x - b.x, a.y - b.y);
}

vec2 operator * (float s, vec2 a) {
    return vec2(s * a.x, s * a.y);
}

You have a choice between de Casteljau's method, which is to recursively split the control path until you arrive at the point using a linear interpolation, as explained above, or Bezier's method which is to blend the control points.

Bezier's method is

 p = (1-t)^3 *P0 + 3*t*(1-t)^2*P1 + 3*t^2*(1-t)*P2 + t^3*P3 

for cubics and

 p = (1-t)^2 *P0 + 2*(1-t)*t*P1 + t*t*P2

for quadratics.

t is usually on 0-1 but that's not an essential - in fact the curves extend to infinity. P0, P1, etc are the control points. The curve goes through the two end points but not usually through the other points.

Did you use a C# library earlier?

In C++, no standard library function for Bezier curves is available (yet). You can of course roll your own (CodeProject sample) or look for a math library.

This blogpost explains the idea nicely but in Actionscript. Translation should not be much of a problem.

To get an individual point (x, y) along a cubic curve at a given percent of travel (t), with given control points (x1, y1), (x2, y2), (x3, y3), and (x4, y4) I expanded De Casteljau’s algorithm and rearranged the equation to minimize exponents:

//Generalized code, not C++
variables passed to function:   t,  x1, y1, x2, y2, x3, y3, x4, y4
variables declared in function: t2, t3, x,  y
t2 = t * t
t3 = t * t * t
x = t3*x4 + (3*t2 - 3*t3)*x3 + (3*t3 - 6*t2 + 3*t)*x2 + (3*t2 - t3 - 3*t + 1)*x1
y = t3*y4 + (3*t2 - 3*t3)*y3 + (3*t3 - 6*t2 + 3*t)*y2 + (3*t2 - t3 - 3*t + 1)*y1

(t) is a decimal value between 0 and 1 (0 <= t <= 1) that represents percent of travel along the curve.

The formula is the same for x and y, and you can write a function that takes a generic set of 4 control points or group the coefficients together:

t2 = t * t
t3 = t * t * t

A = (3*t2 - 3*t3)
B = (3*t3 - 6*t2 + 3*t)
C = (3*t2 - t3 - 3*t + 1)

x = t3*x4 + A*x3 + B*x2 + C*x1
y = t3*y4 + A*y3 + B*y2 + C*y1

For quadratic functions, a similar approach yields:

t2 = t * t

A = (2*t - 2*t2)
B = (t2 - 2*t + 1)

x = t2*x3 + A*x2 + B*x1
y = t2*y3 + A*y2 + B*y1