B Spline interpolation

Discussion in 'Mac Programming' started by Duke Leto, Aug 2, 2010.

  1. macrumors regular

    Joined:
    Mar 17, 2008
    #1
    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:
    [​IMG]

    The output:
    [​IMG]

    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!!
     
  2. macrumors 6502

    Joined:
    Apr 24, 2008
    #2
    I know this has nothing to do with your question, but it looks like you're leaking memory in this function...
     
  3. ncl
    macrumors member

    Joined:
    Aug 16, 2008
    #3
    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.
     
  4. macrumors 6502a

    Joined:
    May 10, 2009
    Location:
    Des Moines, WA
    #4
    And the parameter 'pts' is never referenced within the function.
     
  5. thread starter macrumors regular

    Joined:
    Mar 17, 2008
    #5
    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.
     
  6. macrumors 6502a

    Joined:
    Jun 8, 2005
    Location:
    West-Europe
    #6
    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?
     
  7. thread starter macrumors regular

    Joined:
    Mar 17, 2008
    #7
    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.
     
  8. macrumors 6502a

    Joined:
    Jun 8, 2005
    Location:
    West-Europe
    #8
    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.
     

    Attached Files:

  9. thread starter macrumors regular

    Joined:
    Mar 17, 2008
    #9
    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.
     
  10. macrumors 6502a

    Joined:
    Jun 8, 2005
    Location:
    West-Europe
    #10
    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.

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

    Joined:
    Mar 17, 2008
    #11
    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)
     
  12. macrumors newbie

    Joined:
    Jan 20, 2011
    Location:
    Germany
    #12
    Highlighting and dragging points with GLUT

    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!
     

Share This Page