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

Check this out.

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

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". –  muntoo 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. –  Ｊ. Ｍ. 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. –  Ｊ. Ｍ. 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;
}
};

-

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.

-