B Spline interpolation

[Duke Leto]

I'm trying to implement de Boor's algorithm for B splines, but I'm having some issues.

The rest of the code contains the control points, and it creates a uniform integer knot vector {0, 1, 2, 3, .... } It goes through the knots and calculates the point on the spline for that point. After that, it creates a .png file using the results from the de Boor algorithm.

Point coxDeBoor(float* u, Point* pts, const float x, const int n, int d)
{
	// u is the knot vector
	// pts contains the control points
	// x is the parameter we're calculating for
	// n is the number of control points
	// d is the degree
	
	int i, l, j = 0;
	
	// find the knot span that x is in
	for(i=d-1; i<=n+1; i++)
	{
		if (x>=u[i] && x<u[i+1]) {
			l = i;
			break;
		}
	}
	
	float tau, X, Y;
	
	// initialize our points
	Point** P = (Point**)malloc(sizeof(Point*)*d);
	Point testPt;
	for(i=0; i<d; i++)
	{
		P[i] = malloc(sizeof(Point)*i);
		P[i][0] = testPt;
	}
	
	// these variables are used to spread out the equations (for debugging purposes)
	Point old, new;
	float u1, u2;
	
	// loop through the levels
	for(i=1; i<d; i++)
	{
		// loop through the affected points
		for(j=i; j<d; j++)
		{
			// the calculations
			u1 = u[j];
			u2 = u[j+d-i];
			tau = (x-u1)/(u2-u1);
			
			old = P[j-1][i-1];
			new = P[j][i-1];
			
			X = (1-tau)*old.x + tau*new.x;
			Y = (1-tau)*old.y + tau*new.y;
			
			P[j][i].x = X;
			P[j][i].y = Y;
		}
	}
	
	// return the last value
	return P[d-1][d-1];
}

After a couple of days of debugging, I'm pretty sure the problem is in here. Here are some .png files of the output.

The control points:


The output:


As you can see, there are skips in the spline for some reason, and it goes past the control points, instead of staying well inside. I've also determined that the jumps happen when the x value passes knot values. This leads me to believe that the problem is in the calculation of tau, but I'm not sure. Any help is appreciated!!

[Sander]

I know this has nothing to do with your question, but it looks like you're leaking memory in this function...

[ncl]

I am not familiar with this but there is something wrong in your initialization. Shouldn't

P[i] = malloc(sizeof(Point)*i);

be

P[i] = malloc(sizeof(Point)*d);

?
And testPt is not initialized. Also, it looks like pts and l aren't used in this function. And finally, P is never freed.

If I understand correctly the wikipedia article, P[*][0] should be initialized with the control points in l-d…l and the "for j" loop should go from l-d+i to l.

[lloyddean]

And the parameter 'pts' is never referenced within the function.

[Duke Leto]

Yes, I am leaking memory, but I am not too worried about that at the moment.

In the original code, I used the pts parameter to set testPoint so that I could see it. That is where the initialization occurs. I accidentally changed it in this code. Similarly, the second array did go to d, not i. That was a change I made to test something.

[MrFusion]

Originally Posted by Duke Leto

...it goes past the control points, instead of staying well inside.

That doesn't surprise me. At least not if it is a third order, which I am most familiar with. If you calculate the third order (by hand if you are so inclined), you will see that there are 4 additional points required to complete the recursion. Depending on the definition of your control point, it is either 4 points after your last knot point, or two before the first and two after the last one. These four points have to be defined on the knot axis, but don't require a value for the control point.

Can you also post the recursive formula in math notation? Your implementation differs from the formula I have.
Also, x can be any value in the range you want to plot. The spline curve is a function of x, how can you then calculate it inside the loop?

[Duke Leto]

I have the extra knot values. The control points are supposed to be a convex hull for the spline, so they must stay inside. There are n+d knot values, and x ranges from d to n.

I have fixed a good bit by changing the tau calculation to:

	ii = j+l-d+1;
	u1 = u[ii];
	u2 = u[ii+d-i+1];
	tau = (x-u1)/(u2-u1);

However, there are gaps. Whenever x passes a knot value, the spline point skips a bit.

[MrFusion]

I see that you calculate only 1 tau, and x as well.
The B-splines as I am familiar with them look like in the attachments.
That is why I was asking you to post your formula in math notation.

[Duke Leto]

I'm using deBoor's algorithm.

Anyway, its all fixed now, so if you want to see a sample (I made this today from boredom), then here's a
download.

Its a GLUT application.
n to add a point (in the location of the mouse)
d t delete a point
click/drag to move them
-/= to decrease/increase resolution
2-9 to change the degree
p to toggle viewing the control polygon

I'll post the code if its requested.

[MrFusion]

Originally Posted by Duke Leto

I'm using deBoor's algorithm.

I am not a mathematician nor a computer scientist. I don't know deBoor's algorithm. I am sure there are other people here who also don't know deBoor's algorithm. That doesn't mean we are not capable of finding an error in your code if you had also provided the algorithm.

I'll post the code if its requested.

That would be nice for anyone stumbling on to this thread who has the same problem.

[Duke Leto]

deBoor:
http://en.wikipedia.org/wiki/De_Boor's_algorithm
http://www.idav.ucdavis.edu/educatio...alculation.pdf

^^
The second link is better at explaining, the first gives you a general idea.

Code:

struct pt {
	float x;
	float y;
};
typedef struct pt Point;

void generateRegularBSpline(const int n, Point* controlPts, const int d, Point* splinePts, int numSplinePts)
{
	float*u = malloc(sizeof(float)*(n+d));
	int i;
	for(i=0; i<n+d; i++)
	{
		u[i]=i;
	}
	float begin = u[d];
	float end = u[n];
	float x;
	for(i=0; i<numSplinePts; i++)
	{
		x = (begin + i*(end-begin)/numSplinePts);
		splinePts[i] = coxDeBoor(u, controlPts, x, n, d);
	}
	free(u);
}

Point coxDeBoor(float* u, Point* pts, const float x, const int n, int d)
{
	// u is the knot vector
	// pts contains the control points
	// x is the parameter we're calculating for
	// n is the number of control points
	// d is the degree
	
	int i, l, j = 0;
	// find the knot span that x is in
	for(i=d-1; i<=n+1; i++)
	{
		if (x>=u[i] && x<u[i+1]) {
			l = i;
			break;
		}
	}
	
	float tau, X, Y;
	
	// initialize our points
	Point** P = (Point**)malloc(sizeof(Point*)*(d+1));
	for(i=0; i<=d; i++)
	{
		P[i] = (Point*)malloc(sizeof(Point)*(i+1));
		P[i][0] = pts[l-d+i];
		for(int k = 1; k < i+1; k++)
		{
			P[i][k].x = 0;
			P[i][k].y = 0;
		}
	}
	
	// these variables are used to spread out the equations (for debugging purposes)
	Point old, new;
	float u1, u2;
	//int ii;
	// loop through the levels
	for(i=1; i<=d; i++)
	{
		// loop through the affected points
		for(j=i; j<=d; j++)
		{
			// the calculations
			u1 = u[j+l-d];
			u2 = u[j+l-i+1];
			tau = (x-u1)/(u2-u1);
			
			old = P[j-1][i-1];
			new = P[j][i-1];
			
			X = (1-tau)*old.x + tau*new.x;
			Y = (1-tau)*old.y + tau*new.y;
			
			P[j][i].x = X;
			P[j][i].y = Y;
		}
	}
	
	Point returnPoint = P[d][d];
	for(i=0; i<=d; i++)
	{
		free(P[i]); 
	}
	free(P);
	// return the last value
	return returnPoint;
}

The function you call is generateRegularBSpline. It can be modified it you want to change the knot vector.

(const int n, Point* controlPts, const int d, Point* splinePts, int numSplinePts)
n: number of control points
controlPts: pointer to the first control point
d: degree of the spline
splinePts: pointer to where you want the spline points to be output (its awkward, I know, but you catch the drift)
numSplinePts: number of spline points you want to calculate (and have room for in splinePts)

[little gauss]

Originally Posted by Duke Leto

I'm using deBoor's algorithm.

Anyway, its all fixed now, so if you want to see a sample (I made this today from boredom), then here's a
download.

Its a GLUT application.
n to add a point (in the location of the mouse)
d t delete a point
click/drag to move them
-/= to decrease/increase resolution
2-9 to change the degree
p to toggle viewing the control polygon

I'll post the code if its requested.

I'd be very keen and grateful to learn from you how to highlight and drag rectangular points or similar in OpenGL-GLUT code as you have very nicely demonstrated in your Spline Creator application. I'd be most grateful if you would send me your GLUT code or post it somewhere accessible for my personal study and experimentation. Many thanks and best regards and wishes in advance!