# Casteljau's algorithm - practical example

I have a dataset with about 50 points (x,y) and I would like to draw a smooth curve that can pass as closer as possible on those points.

I have heard about Casteljau's algorithm for splines but after hours searching on google I was not able to find a single piece of code I can use.

As far as I understood, to use this algorithm, I have to divide my dataset in groups of 4 points, right? 1234 5678 etc.. and as far as I noticed, my only problem is to find points in the middle. I mean, if I am calculating a curve for points 1234, I already have points 1 and 4 and I need to calculate 2 and 3, right? But it is a mystery to me how to do that.

I would like to ask you guys if you know some code in C, C++ or Objective-C that computes the curves based on datasets with any amount of number.

What I need is: I send the code an array with the dataset and I receive back an array with the points to draw.

My math is crap. So, please give me practical examples. No need to send me to pages with math theory. Looking at these pages makes my brain hurt...

thanks.

-
Since it's a programming example you're after, this is a good candidate to migrate to StackOverflow. – hardmath Jun 7 '11 at 20:54
there they tell me to post here, here you tell me to post there... – Mike Jun 7 '11 at 21:23
As I understand it, Casteljau's algorithm does not tell you how to compute a Bezier curve approximating a set of points. Rather, once you know the equation for the Bezier curve, Casteljau's algorithm gives you a numerically stable way of computing the points on the curve. Also, it seems quite unlikely that someone is going to write code from scratch for you... – Chris Taylor Jun 7 '11 at 21:35
You need to specify the problem a bit more. When you say that you want a smooth curve passing as closely as possible to the points, what do you mean by "smooth", and what is the maximal polynomial order permitted for the curve? Overhauser's cubic will give an interpolation through the points which IIRC has a continuous second derivative. – Peter Taylor Jun 7 '11 at 22:29
I want it to be as much as closer to the real curve (=points) as possible, just smoother. It can be any degree, but it has to be fast and work from any number of points (x,y). – Mike Jun 7 '11 at 22:36

Here's an implementation of the Casteljau algorithm that I just wrote (in Java, though you should be able to convert to C/C++/Ob-C with little effort - I didn't use any high-level language features).

I haven't tried it so I don't know if it's correct - the best I can say is that it compiles. Here your list of $x_i$ is given in the array x, your $y_i$ are given in y and $n$ is the number of points. It takes $O(n^2)$ time to compute the location of one point along the Bezier curve fitting your points, and if you want to plot the curve with $k$ points, it will take $O(kn^2)$ time to run.

public class Casteljau {

private double[] x;
private double[] y;
private int n;

private double[][] b;

public Casteljau(double[] x, double[] y, int n) {
//require x.length = y.length = n
this.x = x;
this.y = y;
this.n = n;
this.b = new double[n][n];
}

private void init(double[] initialValues) {
for (int i = 0; i < n; i++) {
b[0][i] = initialValues[i];
}
}

private double evaluate(double t, double[] initialValues) {
init(initialValues);
for (int j = 1; j < n; j++) {
for (int i = 0; i < n - j; i++) {
b[j][i] = b[j-1][i] * (1-t) + b[j-1][i+1] * t;
}
}
return(b[n-1][0]);
}

public double[] getXYvalues(double t) {
double xVal = evaluate(t, x);
double yVal = evaluate(t, y);
return new double[] {xVal, yVal};
}

}

-
Don't classes count as "high-level"? Or at least, they should be "middle-level". – Mateen Ulhaq Jun 8 '11 at 1:06
thanks!!!!!!!!! – Mike Jun 8 '11 at 2:49
The algorithm as written here can be streamlined a bit. In particular, a triangular recursion such as de Casteljau's can be implemented to use a one-dimensional array instead of the triangular array that is normally used to display the algorithm. This requires careful overwriting, though. – J. M. Nov 24 '11 at 8:24
If I understand correctly, that would reduce the memory requirements to $O(n)$ but keep the time requirements at $O(n^2)$? – Chris Taylor Nov 24 '11 at 8:53
Yep, the routine is still $O(n^2)$ in effort. At least the storage is reduced. – J. M. Dec 5 '11 at 15:28

@Chris Taylor's implementation in C++:

struct XY
{
double x;
double y;
};

class Casteljau
{
private:
double x[];
double y[];
int n;

double b[][];

void init(double initialValues[])
{
for(int i = 0; i < n; i++)
{
b[0][i] = initialValues[i];
}
}

double evaluate(double t, double initialValues[])
{
init(initialValues);
for(int j = 1; j < n; j++)
{
for(int i = 0; i < n - j; i++)
{
b[j][i] = b[j-1][i] * (1-t) + b[j-1][i+1] * t;
}
}
return(b[n-1][0]);
}

public:
Casteljau(double in_x[], double in_y[], int in_n)
{
x = in_x;
y = in_y;
n = in_n;
b = new double[n][n];
}

XY getXYvalues(XY &xy, double t)
{
double xVal = evaluate(t, x);
double yVal = evaluate(t, y);
xy.X = xVal;
xy.Y = yVal;
}
};

-

Check this out.

-
my brain hurts. I need a programming example. – Mike Jun 7 '11 at 20:35

The premise of this question is wrong. Always has been. If you have 50 points and you want to smoothly curve to all of them. Do a smooth curve. It's even built into things like SVG. Basically you set one control point and interpolate control points for all the other points as the reflection of the point vs the previous point. Then the curve goes through every given point and has a nice curvy bit, that can be set based on the position of the one original control point. Bezier curves at n(50) do not actually go that close to all 50 of those points.

-
Nice, I didn't know of that reflection trick. – Yves Daoust Jul 10 at 8:01
It's how one implements the s and t commands for smooth cubic and smooth quad bezier curves as expressed in the SVG specifications. What's more basically everything lets you use SVG so you could just use an SVG rather than implement anything. Just path d = q0,0,?,?,1,1,t x y x y x y x y x y (for all 50 points) ... w3.org/TR/SVG/paths.html#PathDataQuadraticBezierCommands – Tatarize Jul 10 at 8:31
In fact, here's something that just randomly makes up 50 points and draws it as an svg. jsfiddle.net/Tatarize/3n0stp46 – Tatarize Jul 10 at 8:58

Splines and Bezier curves are two different things. The former are usually interpolating, i.e. they pass through the given control points, while the latter pass "close to" them, in a certain sense.

As said by others, a Bezier of degree 50 is certainly not the right approach, it will strongly filter out the variations and won't be that close.

For interpolation, on prefers "piecewise curves", i.e. curves decomposed in several arcs with different equations, ensuring a smooth connection at the control points.

You can refere here: http://www.codeproject.com/Articles/560163/Csharp-Cubic-Spline-Interpolation. Beware to avoid the confusion between function interpolation (you want to reconstruct a function $Y = f(X)$) and curve interpolation (parametric equations $X = f(t)$, $Y = g(t)$). You are in the second case.

-

To interpolate a set of points smoothly I wouldn't use Bezier curves in the way you suggest, I would use something like a Catmull-Rom spline. See this link for details and source code, albeit using DirectX. I'm sure you could find other public libraries available for constructing Catmull-Rom splines if that doesn't suit your requirements.

-