De Boor’s Algorithm in C++

Here is the c++ code for the De Boor algorithm. It calculates a point C(x) on a B-Spline curve of any degree.

Point deBoor(int k,int degree, int i, double x, double* knots, Point *ctrlPoints)
{   // Please see wikipedia page for detail
	// note that the algorithm here kind of traverses in reverse order
	// comapred to that in the wikipedia page
	if( k == 0)
		return ctrlPoints[i];
	else
	{   
		double alpha = (x-knots[i])/(knots[i+degree+1-k]-knots[i]);
		return (deBoor(k-1,degree, i-1, x, knots, ctrlPoints)*(1-alpha )+deBoor(k-1,degree, i, x, knots, ctrlPoints)*alpha );
	}
}

The Point class is defined as the follow:

class Point
{
public:
	Point(){x=0.;y=0.;z=0.;};
	// copy operator
	Point operator=(const Point pt) ;
	Point operator+(const Point pt) const;
	//Point operator-(const Point pt) const;
	Point operator*(double m) const;
	Point operator/(double m) const;
	double x,y,z;
};

Point Point::operator=(const Point pt)
{
	x = pt.x;
	y = pt.y;
	z = pt.z;
	return *this;
}
Point Point::operator+(const Point pt) const
{
	Point temp;
	temp.x = x + pt.x;
	temp.y = y + pt.y;
	temp.z = z + pt.z;
	return temp;
}
Point Point::operator*(double m) const
{
	Point temp;
	temp.x = x*m;
	temp.y = y*m;
	temp.z = z*m;
	return temp;
}
Point Point::operator/(double m) const
{
	Point temp;
	temp.x = x/m;
	temp.y = y/m;
	temp.z = z/m;
	return temp;
}
Advertisements

22 Responses to De Boor’s Algorithm in C++

  1. Robert says:

    Great algortithm. It looks nice and clean, however: I can’t get it to run.
    Do you refer to the wikipedia article at “http://en.wikipedia.org/wiki/De_Boor_algorithm”?
    And, what are the parameters k and i.

    I am running it with 11 knots and 7 control points, then use degree = 3 (11-7-1), k=3 and i=6 (7-1). I got a result, but it looks not right to me.

    Sorry, but it looks like I am missing a point. Thx a lot for your help.

    • chi3x10 says:

      Hi, Robert,
      That’s the part I didn’t implement well.

      “i” should be the interval “x” value resides in the knot vectors. For example, given knots vector [0, 1, 2, 3, 4, 5, 6, 7, 8], if “x” is 2.5 then “i” is 2. Here is the code I used to determine the value of i (given x and the knot vectors). knot is the knot vectors and ti is the total number of knots.

      int WhichInterval(double x, double *knot, int ti)
      {
      for(int i=1;i<ti-1;i++)
      {
      if(x<knot[i])
      return(i-1);
      else if(x == knot[ti-1])
      return(ti-1);
      }
      return -1;
      }

      • mjn says:

        Hi chi3x10.
        I think there are some problems in your function WhichInterval().
        In the for loop, variable i can be ti – 1. Here is my modification:
        int WhichInterval(double x, double *knot, int ti)
        {
        int index = -1;

        for(int i = 1; i <= ti – 1; i++)
        {
        if(x < knot[i]) {
        index = i – 1;
        break;
        }
        }
        if(x == knot[ti – 1]) {
        index = ti – 1;
        }
        return index;
        }
        Am I right?
        mjn

      • chi3x10 says:

        I don’t think there is a problem. I actually used this code in my other projects. Maybe you use different parameterization?

      • AD says:

        Hi
        I think there is a problem :
        If “x” is 0.5 with [0,1,2,3,…] as knot vector,
        then “i” is 0.
        But in the first loop recursion of deBoor function, i-1 is called.
        However knots[-1] doesn’t exit.
        Am I right?
        A.D

      • AD says:

        Oh sorry, I’m wrong.
        In fact, I work with open curves. Therefore the domain of my curve is [u_p,u_(n-p)] where u is the knot vector, p the degree of my B-Spline and n+1 the number of control points.
        I apologize.
        A.D

  2. Robert says:

    Hi chi3x10,

    thank you very much for your quick reply and this explanation. It’s looking much better now and I see the match with the wikipedia page.

    Cheers,

    Robert

  3. Robert says:

    Hi again,
    after I got this working I found out that it was not what I am looking for 😦
    I need a approximation method which has only a few “supporting points”, e.g. measured values at varying intervals. I wanted to calculate the values in between. My idea was to put a b-spline “near” the existing points and calculate the approx. values at the missing points. However, using the “x” parameter I am counting along the curve, but not necessarily along the x-axis. So I can’t fill in a x from my x-axis and get the y of the spline.
    Maybe I just did not found the right word to google for the right algorithm.

  4. Никита says:

    Hi Robert,
    Can you give me all algorithm De Boor and send on this mail adres Konstantine92@inbox.ru? Please!!!

  5. wsdb says:

    hi,
    can you tell me whats parameter “int k” in your algorithm

    thanks
    wsdb

  6. Psycho_Dad says:

    Hi chi3x10!

    I’m in the middle of implementing something similar to your algorithm in C, using OpenGL to draw a curve.

    I managed to draw a curve but it won’t start and end at my first and last control point. I know this is related to the multiplicity of the knot values. However if I increment the knot values multiplicity at the first and last knot value in my knot vector the curve would start from the first control point but still won’t end in the last one. It seems as if half of my curve is good, and the other half is
    “weird”, because it will end in “0,0” coordinates.

    So then after my struggling I found your algorithm, copied it bit-by-bit (I’m using ANSI C btw., so I had to modify a bit), but I face the exact same problem with your solution.

    I’m not a mathematician or anything, I’ve had a hard time understanding how the de Boor algorithm works (I roughly followed this exlantion: http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/de-Boor.html).

    Could help me with some information regarding knot values and control points? Let’s suppose I have the following informations at hand:
    “P[]” – control points array with “p” elements.
    “k” – the degree of the B-spline curve approximating the polygonial chain defined by the “p” control points.
    “U[]” – knot values array with “u” elements. I don’t know for sure, but for me the minimum number of knot values is “u” = “p” + “k” + 1. As I increment the multiplicity at “U[0]” and “U[u-1]” this is increasing.
    “X[]” – parameter values array with “x” elements. (I’m using uniform paramter values; I calculate the parameter values with a for loop:
    for(double x = 0.0; x <= 1.0; x += 0.1) { ..

    Could clarify how the elements of my knot vector should be defined to have a curve which start and ends at the first and last control point? If you have any advice I'm all ears. 🙂

    Thanks in advance!

  7. Peter Fong says:

    Hi, if I was to plot the curve out, where should I put my makepixel fuction at?

  8. Ihor says:

    Hello, thanks for the algorithm. I would like to know is your code suitable if there are knots with multiplicities ? Thank you in advance for the answer!

  9. ir says:

    This is not an efficient implementation of de boor algorithm. There are many excessive computations in your code.

    • SonicMouse says:

      why not post an update then?

      What do you mean by “excessive” computations? Are you referring to the Point class’s operator overrides not passing const parameters by reference? Or are you talking about the DeBoor computation only?

  10. SonicMouse says:

    Why not post something on how to call this. Right now your post is in 1st place on Google for “Cox Deboor C++”, but I guarantee people are glazing over it because there is no example on how to use it.

  11. JR says:

    Due to inaccurate precision i suggest using an EPSILON (eg. 1e-10) to determin if two doubles x1, x2 are equal.
    so changing x1 == x2 into abs(x1 – x2) < EPSILON

  12. MarcusHe says:

    chi3x10
    I don’t know whether you can see my question . Pls help me.
    As WikiPedia says , first we should locate where the X is.
    However knotsCount=ControlPointsCount+degrees+1 , for example 6 control points ,10 knots as [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0],3 degrees . I set the inserted u as [0.15,0.25,0.35,0.45,0.55,0.65,0.75,0.85,0.95].
    When X=0.75, i here is 6, could you return ctrlPoints[6]???
    6 is out of index, Now way to go on.
    I don’t know whether i create the right all of the inserted u ??

  13. seek h says:

    Please tell me which header files are needed to be added in this code and can you post the whole source code.

    – Some one new to programming

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: