Mac B Spline interpolation

Duke Leto

macrumors regular
Original poster
Mar 17, 2008
166
0
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.

Code:
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

macrumors 6502
Apr 24, 2008
480
7
I know this has nothing to do with your question, but it looks like you're leaking memory in this function...
 

ncl

macrumors member
Aug 16, 2008
58
0
I am not familiar with this but there is something wrong in your initialization. Shouldn't
Code:
P[i] = malloc(sizeof(Point)*i);
be
Code:
P[i] = malloc(sizeof(Point)*[b]d[/b]);
?
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.
 

Duke Leto

macrumors regular
Original poster
Mar 17, 2008
166
0
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

macrumors 6502a
Jun 8, 2005
601
0
West-Europe
...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

macrumors regular
Original poster
Mar 17, 2008
166
0
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:
Code:
	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.
 

Duke Leto

macrumors regular
Original poster
Mar 17, 2008
166
0
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

macrumors 6502a
Jun 8, 2005
601
0
West-Europe
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

macrumors regular
Original poster
Mar 17, 2008
166
0
deBoor:
http://en.wikipedia.org/wiki/De_Boor's_algorithm
http://www.idav.ucdavis.edu/education/CAGDNotes/Deboor-Cox-Calculation.pdf

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

Code:
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

macrumors newbie
Jan 20, 2011
1
0
Germany
Highlighting and dragging points with GLUT

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!
 
Register on MacRumors! This sidebar will go away, and you'll see fewer ads.