# Rectangular Bezier Patches

Rectangular Bezier Patches are one way to bring Bezier curves into the 3rd dimension as a Bezier surface. Below is a rendered image of a quadratic Bezier rectangle (degree of (2,2)) and a cubic Bezier rectangle (degree of (3,3)) taken as screenshots from a shadertoy demo I created that renders these in real time. Links at bottom of post!

## Intuition

Imagine that you had a Bezier curve with some number of control points. Now, imagine that you wanted to animate those control points over time instead of having a static curve.

One way to do this would be to just have multiple sets of control points as key frames, and just linearly interpolate between the key frames over time. You’d get something that might look like the image below (lighter red = farther back in time).

That is a simple and intuitive way to animate a Bezier curve, and is probably what you thought of immediately. Interestingly though, since linear interpolation is really a degree 1 Bezier curve, this method is actually using a degree 1 Bezier curve to control each control point!

What if we tried a higher order curve to animate each control point? Well… we could have three sets of control points, so that each control point was controlled over time by a quadratic curve. We could also try having four sets of control points, so that each control point was controlled over time by a cubic curve.

We could have any number of sets of control points, to be able to animate the control points over time using any degree curve.

Now, instead of animating the curve over TIME, what if we controlled it over DISTANCE (like, say, the z-axis, or “depth”). Look at the image above and think of it like you are looking at a surface from the side. If you took a bunch of the time interpolations as slices and set them next to each other so that there were no gaps between them, you’d end up with a smooth surface. TA-DA! This is how a Rectangular Bezier Patch is made.

Note that the degree of the curve on one axis doesn’t have to match the degree of the curve on the other axis. You could have a cubic curve where each control point is controlled by a linear interpolation, or you could have a degree 5 curve where each control point is controlled by degree 7 curves. Since there are two degrees involved in a Bezier rectangle, you describe it’s order with two numbers. The first example is degree (3,1) and the second example is degree (5,7).

Higher Dimensions

While you thinking about this, I wanted to mention that you could animate a bezier rectangle over time, using bezier curves to control those control points. If you then laid that out over distance instead of time, you’d end up with a rectangular box Bezier solid. If you are having trouble visualizing that, don’t feel dumb, it’s actually four dimensional!

You can think of it like a box that has a value stored at every (x,y,z) location, and those values are controlled by Bezier formulas so are smooth and are based on control points. It’s kind of a strange concept but is useful in some situations.

Say you made a 3d hot air baloon game and wanted to model temperature of the air at differently locations to simulate thermals. One way you could do this would be to store a bunch of temperatures in a 3d grid. Another way might involve using a grid of rectangular box Bezier solids perhaps. One benefit to the Bezier solid representation is that the data points are much smoother than a grid would be, and another is that you could make the grid much less dense.

Now, let’s say that you wanted to animate the thermals over time. You could use a fifth dimensional bezier hypercube solid. Let’s move on, my brain hurts 😛

## Math

The equation for a Bezier Rectangle is:

$\mathbf{p}(u, v) = \sum_{i=0}^n \sum_{j=0}^m B_i^n(u) \; B_j^m(v) \; \mathbf{k}_{i,j}$

$\mathbf{p}(u, v)$ is the point on the surface that you get after you plug in the parameters. $u$ and $v$ are the parameters to the surface and should be within the range 0 to 1. These are the same thing as the $t$ you see in Bezier curves, but there are two of them since there are two axes.

There are two Sigmas (summations) which mean that it’s a double for loop.

One of the for loops make $i$ go from 0 to $n$ and the other makes $j$ go from 0 to $m$. $m$ and $n$ are the degree of each axis.

$B_i^n(u)$ and $B_i^n(u)$ are Bernstein polynomials (aka binomial expansion terms) just as you see in Bezier Curves – there is just one per axis.

Lastly comes the control points $\mathbf{k}_{i,j}$. The number of control on one axis are multiplied by the number of control points on the other axis.

A biquadratic Bezier patch has a degree of (2,2) and has 3 control points on one axis, and 3 control points on the other. That means that it has 9 control points total.

A bicubic Bezier patch has a degree of (3,3) with 4 control points on each axis, for a total of 16 control points.

If you had a patch of degree (7,1), it would have 8 control points on one axis and 2 control points on the other axis, and so would also have 16 control points total, but they would be laid out differently than a bicubic Bezier patch.

As far as actually calculating points on a curve, the above only calculates the value for a single axis for the final point on the curve. If you have three dimensional control points (X,Y,Z), you have to do the above math for each one to get the final result. This is the same as how it works for evaluating Bezier curves.

## Code

#include
#include

typedef std::array TFloat3;
typedef std::array TFloat3x3;

const TFloat3x3 c_ControlPointsX =
{
{
{ 0.7f, 0.8f, 0.9f },
{ 0.2f, 0.5f, 0.4f },
{ 0.6f, 0.3f, 0.1f },
}
};

const TFloat3x3 c_ControlPointsY =
{
{
{ 0.2f, 0.8f, 0.5f },
{ 0.6f, 0.9f, 0.3f },
{ 0.7f, 0.1f, 0.4f },
}
};

const TFloat3x3 c_ControlPointsZ =
{
{
{ 0.6f, 0.5f, 0.3f },
{ 0.7f, 0.1f, 0.9f },
{ 0.8f, 0.4f, 0.2f },
}
};

void WaitForEnter ()
{
printf("Press Enter to quit");
fflush(stdin);
getchar();
}

float QuadraticBezier (const TFloat3& p, float t)
{
float s = 1.0f - t;
float s2 = s * s;
float t2 = t * t;

return
p[0] * s2 +
p[1] * 2.0f * s * t +
p[2] * t2;
}

float BiquadraticBezierPatch(const TFloat3x3& p, float u, float v)
{
TFloat3 uValues;
}

int main(int argc, char **argv)
{
// how many values to display on each axis. Limited by console resolution!
const int c_numValues = 4;

printf("Bezier rectangle:n");
for (int i = 0; i < c_numValues; ++i)
{
float iPercent = ((float)i) / ((float)(c_numValues - 1));
for (int j = 0; j < c_numValues; ++j)
{
if (j == 0)
printf("  ");
float jPercent = ((float)j) / ((float)(c_numValues - 1));
float valueX = BiquadraticBezierPatch(c_ControlPointsX, jPercent, iPercent);
float valueY = BiquadraticBezierPatch(c_ControlPointsY, jPercent, iPercent);
float valueZ = BiquadraticBezierPatch(c_ControlPointsZ, jPercent, iPercent);
printf("(%0.2f, %0.2f, %0.2f) ", valueX, valueY, valueZ);
}
printf("n");
}
printf("n");

WaitForEnter();
return 0;
}


And here is the output it gives:

Note that in the program above, I evaluate the surface points by evaluating one axis and then the other. This is basically the same as how I explained it at the top, where I’m effectively animating the control points over distance, then evaluating the curve slice of the surface at that specific distance.

You could also write it another way though, where you literally expand the mathematical formula to get just one expression to evaluate that takes all control points at once. I like the simplicity (of understanding) of the method I used, but the other method works just as well.

## The Rendering

It’s easy enough to calculate values on a Bezier Rectangle, but what if you want to draw one?

One way is to tessellate it, or break it up into triangles and then render the triangles. You can think of it like trying to render a grid, where each point of the grid is moved to be where ever the Bezier rectangle function says it should be.

Raytracing against these objects in the general case is very difficult however, because it basically comes down to solving equations of very high degree.

Raymarching against these objects is also difficult unfortunately because while raymarching only needs to know “am i above the shape, or underneath it?”, knowing what u,v to plug into the equation to get the height most relevant to a random point in space is also very difficult. Not as difficult as the raytracing equations, but probably just as much out of reach.

But never fear, as always, you can cheat!

If you read my post about one dimensional (explicit) Bezier curves (One Dimensional Bezier Curves), you may remember that math gets easier if you use one dimensional control points. The same is actually true with Bezier rectangles!

For the ray marching case, you can march a point through space, and plug the x,z coordinate of the point into the Bezier rectangle function as u,v values and the number that comes out you can treat as a y coordinate.

Now, ray marching a Bezier rectangle is the same as ray marching any old height map (check links section for more info on that).

What I did in my demos, is since i knew that the curve was constrained to 0-1 on the x and z axis, and the y axis min and max was the control point min and maxes, I did a raytrace of that bounding box to get a minimum and maximum distance that the ray was inside that box. From there, I did raymarching from that min time to the max time along the ray, considering the ray as hitting the surface whenever the distance from the ray to the surface on the y axis (rayPos.y – bezierRectangle.y) changed sign.

After I had a hit, I got the height of the curve slightly offset on the x axis, then slightly offset on the z axis to get a triangle that I could calculate a surface normal from, to do lighting and shading with.

There is room for improvement in the ray marching though. I evenly divide the space through the box by a specific amount to control the size of the steps. A better way to do this I think would be to get the gradient of the function and use that to get a distance estimate (check links section below for more information). I could use that value to control the distance the ray marches at each step, and should be able to march through the box much quicker.

Also, as the link on terrain marching explains, you can usually take farther steps when the ray is farther from the camera, because the eye notices less detail. I removed that since the Bezier rectangles are pretty close to the camera, but it probably still would be helpful. Also, it would DEFINITELY be helpful in the case of the “Infinite Bezier Rectangles” scene.

I am pretty sure you could directly raytrace an explicit Bezier rectangle (one who has one dimensional control points) – at least for low degrees. I personally don’t know how you would do that, but I think it might boil down to solving a 4th degree function or something else “reasonable” based on a similar question I had about Bezier triangles on the mathematics stack exchange site (link below).

Another Way To Render

There is another way to render Bezier surfaces using ray based methods that I didn’t use but want to mention.

A property of Bezier curves and surfaces is that they are guaranteed to be completely contained by the convex hull created by their control points.

Another property of Bezier curves and surfaces is that you can use the De Casteljeau algorithm to cut them up. For instance you could cut a Bezier curve into two different Bezier curves, and the same holds for Bezier surfaces.

Using these two properties, there is an interesting way to be able to tell if a ray intersects a bezier curve or not, which is:

1. If the line misses the convex hull, return a miss
2. If the convex hull is smaller than a pixel, return a hit
3. Otherwise, cut the Bezier object into a couple smaller Bezier objects
4. Recurse for each smaller Bezier object

Yes, believe it or not, that is a real technique! It’s called Bezier Clipping and there is a research paper in the links section below that talks about some of the details of using that rendering technique.

Lastly, I wanted to mention that the above is completely about Bezier rectangles, but there is no reason you couldn’t extend these rectangles to use rational Bezier functions, or be based on B-splines or NURBS, or even go a different direction and make hermite surfaces or catmull-rom surfaces, or even make surfaces that used exotic basis functions of your own crafting based on trigonometric functions or whatever else!

# Lagrange Interpolation

Lagrange interpolation is a way of crafting a function from a set of data points..

In the past I’ve seen reference to Lagrange interpolation in relation to audio programming like, for helping make a soft knee for a limiter, but it can be used wherever you need to make a function from some data points.

## What’s It Do?

Lagrange interpolation is a way of crafting a $y=f(x)$ function from a set of $(x,y)$ data pairs. The resulting function passes through all the data points you give it (like a Catmull-Rom spline does), so can be used to find a function to interpolate between data sets.

You can’t give two value pairs that have the same x value, but the data points don’t have to be evenly spaced.

Also, if you give $N$ data points, you’ll get out a function that is a $N-1$ degree polynomial. So, if you interpolate two data points, you’ll get a degree 1 polynomial (a line). If you interpolate three data points, you’ll get a degree 2 polynomial (a quadratic).

The function will be quite messy, but you can use algebra, or wolframalpha.com (or the like) to simplify it for you to a simpler equation.

Lagrange interpolation is subject to Runge’s Phenomenon, so the more data points you have, the more the interpolation tends to get “squirly” near the edges and shoot off up high or down low, instead of smoothly interpolating between data values.

## How’s It Do It?

Well, to make any kind of curve from data points, if we want the curve to pass through those data points, one way would be to come up with a set of functions to multiply each data point by.

Each function must evaluate to 1 when the curve is at that control point, it should be zero when the curve is at any other control point. Between control points, the function can take any value, but if you make it continuous / smooth, the curve will be continuous and smooth, so that’s usually what is desired.

When we have those functions, to get a point on the curve we just multiply each control point by it’s corresponding function (called a basis function), and we sum up the results.

The pseudocode below is how this works and is the basic functionality of most common curve types:

// The basic way to evaluate most any type of curve
float PointOnCurve (float t, float *controlPoints, int numControlPoints)
{
float value = 0.0f;

for (int i = 0; i < numControlPoints; ++i)
value += controlPoints[i] * ControlPointFunction(i, t);

return value;
}

float ControlPointFunction (int i, float t)
{
// return the ith control point function evaluated at time t.
// aka return f(t) for the ith control point.
}


What makes Lagrange interpolation different than other curve types is the basis functions it uses.

## The Math

If you aren’t used to seeing a capital pi, or a laplacian style cursive l in equations, it’s about to get a bit mathy!

If you feel like skipping to the next section, I don’t blame you, but if you are feeling brave, you should try and follow along, because I’m going to slowly walk through each and every symbol to help explain what’s going on and why.

Let’s say that you are want to be able to interpolate between $k+1$ data points:

$(x_0, y_0)\ldots(x_k, y_k)$

The formula for calculating a Lagrange interpolated value is this:

$L(x) := \sum_{j=0}^{k} y_j \ell_j(x)$

The capital sigma ($\sum_{j=0}^{k}$) just means that we are going to loop a variable j from 0 to k (including k), and we are going to sum up the total of everything on the right for all values of j. When you see a capital sigma, think sum (note they both start with an s).

The next thing after the sigma is $y_j$. That is just the y value from our jth control point. That is essentially controlPoints[j].y.

After that comes the last part $\ell_j(x)$. That is just the function for the jth control point that we multiply the control point by (aka the basis function), evaluated for the specific value x.

Since there is no operator between this function and the control point, that means we multiply them together. So yeah… that crazy math just says “multiply each control point by it’s basis function, and sum up the results”, just like our pseudo code above does!

The second equation we need to look at is the definition of the basis functions for each control point. Here is the formula that describes the jth basis function, for the jth control point:

$\ell_j(x) := \prod_{\begin{smallmatrix}0\le m\le k\\ m\neq j\end{smallmatrix}} \frac{x-x_m}{x_j-x_m}$

First is the capital pi $\prod_{\begin{smallmatrix}0\le m\le k\\ m\neq j\end{smallmatrix}}$. This means that we are going to do a loop, but instead of adding the results of the loop, we are going to multiply them together. Where a capital sigma means sum, capital pi means product.

The notation for product is a bit different here than in the sigma though which may be a bit tricky to read at first. Instead of explicitly saying that $m$ should go from 0 to $k$, the notation $latex 0\le m\le k\\$ says that implicitly. That same notation can be used with sigma, or the more explicit style notation could be used with pi.

The pi also has this notation next to it $m\neq j$. That means that the case where $m$ equals $j$ should be skipped.

Finally, on to the second part: $\frac{x-x_m}{x_j-x_m}$. This part is pretty easy to read. $x$ is the parameter to the function of course, $x_m$ is just controlPoints[m].x where $m$ is the index variable of our product loop ($\prod$), and $x_j$ is just controlPoints[j].x where $j$ is the index variable of our summation loop ($\sum$).

Let’s say that $k$ was 2 because we had 3 data pairs. Our three basis functions would be:

$\ell_0(x) := \frac{x-x_1}{x_0-x_1} * \frac{x-x_2}{x_0-x_2}$
$\ell_1(x) := \frac{x-x_0}{x_1-x_0} * \frac{x-x_2}{x_1-x_2}$
$\ell_2(x) := \frac{x-x_0}{x_2-x_0} * \frac{x-x_1}{x_2-x_1}$

Which means that our final Lagrange interpolation function would be:

$L(x) := y_0 * \frac{x-x_1}{x_0-x_1} * \frac{x-x_2}{x_0-x_2} + y_1 * \frac{x-x_0}{x_1-x_0} * \frac{x-x_2}{x_1-x_2} + y_2 * \frac{x-x_0}{x_2-x_0} * \frac{x-x_1}{x_2-x_1}$

That is quite a mouth full, but hopefully you understand how we came up with that!

$x_i$ is just controlPoints[i].x and $y_i$ is just controlPoints[i].y.

## Math Intuition

The intuition here is that we need to come up with a set of functions to multiply each control point by, such that when the function’s x value is at the control point’s x value, the function should evaluate to 1. When the function’s x value is at a different control points x value, the function should evaluate to 0. The rest of the time, the function can evaluate to whatever it wants, although again, having it have smooth values is nice to making a good curve.

So the first problem is, how do we make a function evaluate to 0 when x is at a different control point?

The easy way would be to multiply a bunch of terms together of this form $(x - x_i)$, but make sure and not include the x of the actual control point that we are multiplying against.

That is exactly what it does with the numerator in the product notation of the basis function.

$\ell_j(x) := \prod_{\begin{smallmatrix}0\le m\le k\\ m\neq j\end{smallmatrix}} \frac{x-x_m}{x_j-x_m}$

Note that $j$ is the index of the current control point that we are calculating the basis function for. All values of x, that isn’t the x value of a control point will evaluate to non zero.

The denominator value is there so that when x is the value of the control point that we care about, that the function will evaluate to 1.

It does this by figuring out what the value of the numerator will be when x is at the control point, and then makes that be the value that it divides by, so that it’s 1 at that x value.

Not too much to it. Pretty simple stuff, but powerful as well!

## Extending to 2D and Beyond

Lagrange interpolation is a one dimensional interpolation scheme, meaning that if you have data points of the form (x,y), it can give you an interpolated y value based on an x value you give it. The interpolation it does can never give two different y values for the same x.

If you want to extend this technique to interpolating a curve through two dimensional data points, or even higher, you need to do interpolation independently for each axis and use a “parametric” value for that axis.

For instance, if you needed to interpolate a curve through 3 dimensional points, you would have data points like this:

X Points = $(t_{x,0}, x_0)\ldots(t_{x,k+1}, x_{k+1})$
Y Points = $(t_{y,0}, y_0)\ldots(t_{y,k+1}, y_{k+1})$
Z Points = $(t_{z,0}, z_0)\ldots(t_{z,k+1}, y_{k+1})$

And then you would interpolate on each axis by the t value to get your X, Y and Z axis values. This should look familiar, because this is how higher dimensional Bezier curves work; you evaluate them per axis based on a parametric value per axis (s,t,u,etc).

You could use the same t values on each axis, or they could be completely independent. You don’t even need to have the same number of points for each axis!

You might wonder how this differs from the standard interpolation in the 2D case. Check the demos in the link section below to really get a grasp of the difference, but in essence, with standard (1D) interpolation, you can never have two x values that evaluate to 2 different y values. Extending it like the above into two dimensions by parameterizing each axis lets you get around that limitation and you can make true 2d shapes.

Lastly, it is possible to make Lagrange interpolated surfaces! I won’t go into the details (perhaps a future post!), but if you know how to make a bezier rectangle by doing a tensor product (basically having X axis Bezier curves, multiplied by Y axis Bezier curves), you can accomplish a Lagrange surface in a really similar way.

## Sample Code

This sample code is written for readability, but could easily be optimized for faster execution. Also, from what I hear, the second form of Barycentric Lagrange Interpolation is touted as the fastest form of Lagrange interpolation, since many values can be pre-calculated and re-used for different values of x.

#include <stdio.h>
#include <vector>

struct SPoint
{
float x;
float y;
};

typedef std::vector<SPoint> TPointList;

void WaitForEnter ()
{
printf("Press Enter to quit");
fflush(stdin);
getchar();
}

// calculates the lagrange basis function y value for control point "controlPointIndex" at x value "x"
float LagrangeBasis (const TPointList& pointList, size_t controlPointIndex, float x)
{
// this is the pi "inner loop" multiplication work
float value = 1.0f;
for (size_t i = 0, c = pointList.size(); i < c; ++i) {
if (i != controlPointIndex)
value *= (x - pointList[i].x) / (pointList[controlPointIndex].x - pointList[i].x);
}
return value;
}

// returns a value at x, using lagrange interpolation over the specified list of (x,y) pairs
float LagrangeInterpolate (const TPointList& pointList, float x)
{
// this is the sigma "outer loop" summation work
float sum = 0.0f;
for (size_t controlPointIndex = 0, c = pointList.size(); controlPointIndex < c; ++controlPointIndex)
sum += pointList[controlPointIndex].y * LagrangeBasis(pointList, controlPointIndex, x);
return sum;
}

int main (int argc, char **argv)
{
// show some 1d interpolated values
// note that the points don't need to be sorted on X, but it makes for easier to read examples
{
// (x,y) pairs
const TPointList points =
{
{ 0.0f, 1.1f },
{ 1.6f, 8.3f },
{ 2.3f, 6.5f },
{ 3.5f, 4.7f },
{ 4.3f, 3.1f },
{ 5.9f, 7.5f },
{ 6.8f, 0.0f }
};

// show values interpolated from x = 0, to x = max x
printf("1d interpolated values.  y = L(t)n");
const float c_numPoints = 10;
for (int i = 0; i < c_numPoints; ++i)
{
float percent = ((float)i) / (float(c_numPoints - 1));
float x = points.back().x * percent;
float y = LagrangeInterpolate(points, x);
printf("  (%0.2f, %0.2f)n", x, y);
}
printf("n");
}

// show some 2d interpolated values
// also note that x and y don't have to have matching t values!
{
// (t, x) pairs
const TPointList pointsX =
{
{ 0.0, 0.0f},
{ 1.0, 1.6f},
{ 2.0, 2.3f},
{ 3.0, 3.5f},
{ 4.0, 4.3f},
{ 5.0, 5.9f},
{ 6.0, 6.8f}

};
// (t, y) pairs
const TPointList pointsY =
{
{ 0.0f, 1.1f },
{ 1.0f, 8.3f },
{ 2.0f, 6.5f },
{ 3.0f, 4.7f },
{ 4.0f, 3.1f },
{ 5.0f, 7.5f },
{ 6.0f, 0.0f }
};

// show values interpolated from t = 0, to t = max t, on each axis
printf("2d interpolated values.  x = L(t_x), y = L(t_y)n");
const float c_numPoints = 10;
for (int i = 0; i < c_numPoints; ++i)
{
float percent = ((float)i) / (float(c_numPoints - 1));

// calculate x
float tx = pointsX.back().x * percent;
float x = LagrangeInterpolate(pointsX, tx);

// calculate y
float ty = pointsY.back().x * percent;
float y = LagrangeInterpolate(pointsY, ty);

printf("  (%0.2f, %0.2f)n", x, y);
}
printf("n");
}

WaitForEnter();
return 0;
}


And here’s the programs output:

## Final Notes

Now that you know how to do all this stuff I wanted to share a couple more pieces of info.

Firstly, it’s kind of weird to call this “Lagrange Interpolation”. A better term is to call this the “Lagrange Form of Polynomial Interpolation”. The reason for that is that if you have some number of data points, there exists only one unique minimal order polynomial (lowest degree of x possible) that fits those points. That is due to the “unisolvence theorem” that you can read more about here: Wikipedia: Polynomial interpolation.

What that means is that if you were to use a different type of polynomial interpolation – such as newton interpolation – the result you get out is algebraically equivalent to the one you’d get from this Lagrange form. There are pros and cons to using different forms of polynomials, but that’s out of the scope of this post so go read about them if you are interested!

Speaking of that, even though this sample code is focused on interpolation using the Lagrange form, this technique is really great at being able to just come up with some simpler f(x) function that passes through specific data points. In this way, you can kind of “bake out” a custom f(x) function to do interpolation for specific values, that doesn’t need all the moving parts of the Lagrange form. For example, if you make the formula for lagrange interpolation of 3 specific value pairs and then simplify, will get out a simple quadratic function in the form of $y=Ax^2+Bx+C$!

Here are some interactive demos I made to let you play with Lagrange interpolation to get a feel for how it works, and it’s strengths and weaknesses:
One Dimensional Lagrange Interpolation
Two Dimensional Lagrange Interpolation

I also found these links really helpful in finally understanding this topic:
Lagrange Interpolation
Lagrange’s Interpolation Formula

Want to follow the rabbit hole a little deeper? Check out how sinc interpolation relates to the Lagrange form!
The ryg blog: sinc and Polynomial interpolation

# The De Casteljau Algorithm for Evaluating Bezier Curves

Over the past year or so I’ve been digging fairly deeply into curves, mostly into Bezier curves specifically.

While digging around, I’ve found many mentions of the De Casteljau algorithm for evaluating Bezier curves, but never much in the way of a formal definition of what the algorithm actually is, or practical examples of it working.

Now that I understand the De Casteljau algorithm, I want to share it with you folks, and help there be more useful google search results for it.

The De Casteljau algorithm is more numerically stable than evaluating Bernstein polynomials, but it is slower. Which method of evaluating Bezier curves is more appropriate is based on your specific usage case, so it’s important to know both.

If you are looking for the mathematical equation of a Bezier curve (the Bernstein form which uses Bernstein basis functions), you have come to the right place, but the wrong page! You can find that information here: Easy Binomial Expansion & Bezier Curve Formulas

Onto the algorithm!

## The De Casteljau Algorithm

The De Casteljau algorithm is actually pretty simple. If you know how to do a linear interpolation between two values, you have basically everything you need to be able to do this thing.

In short, the algorithm to evaluate a Bezier curve of any order $N$ is to just linearly interpolate between two curves of degree $N-1$. Below are some examples to help show some details.

The simplest version of a Bezier curve is a linear curve, which has a degree of 1. It is just a linear interpolation between two points $A$ and $B$ at time $t$, where $t$ is a value from 0 to 1. When $t$ has a value of 0, you will get point $A$. When $t$ has a value of 1, you will get point $B$. For values of t between 0 and 1, you will get points along the line between $A$ and $B$.

The equation for this is super simple and something you’ve probably seen before: $P(t) = A*(1-t) + B*t$.

The next simplest version of a Bezier curve is a quadratic curve, which has a degree of 2 and control points $A,B,C$. A quadratic curve is just a linear interpolation between two curves of degree 1 (aka linear curves). Specifically, you take a linear interpolation between $A,B$, and a linear interpolation between $B,C$, and then take a linear interpolation between those two results. That will give you your quadratic curve.

The next version is a cubic curve which has a degree of 3 and control points $A,B,C,D$. A cubic curve is just a linear interpolation between two quadratic curves. Specifically, the first quadratic curve is defined by control points $A,B,C$ and the second quadratic curve is defined by control points $B,C,D$.

The next version is a quartic curve, which has a degree of 4 and control points $A,B,C,D,E$. A quartic curve is just a linear interpolation between two cubic curves. The first cubic curve is defined by control points $A,B,C,D$ and the second cubic curve is defined by control points $B,C,D,E$.

So yeah, an order $N$ Bezier curve is made by linear interpolating between two Bezier curves of order $N-1$.

## Redundancies

While simple, the De Casteljau has some redundancies in it, which is the reason that it is usually slower to calculate than the Bernstein form. The diagram below shows how a quartic curve with control points $A,B,C,D,E$ is calculated via the De Casteljau algorithm.

Compare that to the Bernstein form (where $s$ is just $(1-t)$)

$P(t) = A*s^4 + B*4s^3t + C*6s^2t^2 + D*4st^3 + E*t^4$

The Bernstein form removes the redundancies and gives you the values you want with the least amount of moving parts, but it comes at the cost of math operations that can give you less precision in practice, versus the tree of lerps (linear interpolations).

## Sample Code

Pretty animations and intuitive explanations are all well and good, but here’s some C++ code to help really drive home how simple this is.

#include <stdio.h>

void WaitForEnter()
{
printf("Press Enter to quit");
fflush(stdin);
getchar();
}

float mix(float a, float b, float t)
{
// degree 1
return a * (1.0f - t) + b*t;
}

float BezierQuadratic(float A, float B, float C, float t)
{
// degree 2
float AB = mix(A, B, t);
float BC = mix(B, C, t);
return mix(AB, BC, t);
}

float BezierCubic(float A, float B, float C, float D, float t)
{
// degree 3
float ABC = BezierQuadratic(A, B, C, t);
float BCD = BezierQuadratic(B, C, D, t);
return mix(ABC, BCD, t);
}

float BezierQuartic(float A, float B, float C, float D, float E, float t)
{
// degree 4
float ABCD = BezierCubic(A, B, C, D, t);
float BCDE = BezierCubic(B, C, D, E, t);
return mix(ABCD, BCDE, t);
}

float BezierQuintic(float A, float B, float C, float D, float E, float F, float t)
{
// degree 5
float ABCDE = BezierQuartic(A, B, C, D, E, t);
float BCDEF = BezierQuartic(B, C, D, E, F, t);
return mix(ABCDE, BCDEF, t);
}

float BezierSextic(float A, float B, float C, float D, float E, float F, float G, float t)
{
// degree 6
float ABCDEF = BezierQuintic(A, B, C, D, E, F, t);
float BCDEFG = BezierQuintic(B, C, D, E, F, G, t);
return mix(ABCDEF, BCDEFG, t);
}

int main(int argc, char **argv)
{
struct SPoint
{
float x;
float y;
};

SPoint controlPoints[7] =
{
{ 0.0f, 1.1f },
{ 2.0f, 8.3f },
{ 0.5f, 6.5f },
{ 5.1f, 4.7f },
{ 3.3f, 3.1f },
{ 1.4f, 7.5f },
{ 2.1f, 0.0f },
};

//calculate some points on a sextic curve!
const float c_numPoints = 10;
for (int i = 0; i < c_numPoints; ++i)
{
float t = ((float)i) / (float(c_numPoints - 1));
SPoint p;
p.x = BezierSextic(controlPoints[0].x, controlPoints[1].x, controlPoints[2].x, controlPoints[3].x, controlPoints[4].x, controlPoints[5].x, controlPoints[6].x, t);
p.y = BezierSextic(controlPoints[0].y, controlPoints[1].y, controlPoints[2].y, controlPoints[3].y, controlPoints[4].y, controlPoints[5].y, controlPoints[6].y, t);
printf("point at time %0.2f = (%0.2f, %0.2f)n", t, p.x, p.y);
}

WaitForEnter();
}


Here’s the output of the program:

Thanks to wikipedia for the awesome Bezier animations! Wikipedia: Bézier curve