Curve Rendering (part 1): B-Spline Curves && de de Boor Algorithm

While taking a graphics class last fall, we touched upon curve rendering (one of the most interesting graphics topics IMO). We learned of two different types of curves - Bezier and B-Spline, and two algorithm to generate these curves - de Casteljau and de Boor. Because I'm not a mathematician, these algorithms were pretty difficult to understand. But after spending a long time studying them, I got a fairly good understanding of how each work. I don't want to bore you with an explanation of how the algorithm works because it can be found everywhere on the web. Instead, I want to share with you my implementation of the two algorithms. In this particular post, however, I will be sharing my implementation of the de Boor algorithm, which is used to generate B-Spline curves:

struct Point {
 Point() : x(0.0), y(0.0) {}
 Point(const double xIn, const double yIn) : x(xIn), y(yIn) {}

 Point operator=(const Point rhs) {
   x = rhs.x;
   y = rhs.y;
   return *this;
 }
 Point operator+(const Point rhs) const {
   return Point(x + rhs.x, y + rhs.y);
 }
 Point operator*(const double m) const {
   return Point(x * m, y * m);
 }
 Point operator/(const double m) const {
   return Point(x / m, y / m);
 }

 double x, y;
};

void deBoor(const Point* ctrlPoints, const int numCtrlPts,
                       const float* knots, const int numKnots,
                       const int k, const int resolution,
                       Point* bsplinePoints) {
 /*
  ctrlPoints -- our control points for this curve
  numCtrlPts -- number of control points
  knots -- list of knots for this curve
  numKnots -- number of knots we currently have, this is equal to (numCtrlPts - 1 + k)
  k -- order of this spline curve
  resolution -- number of spline points
  */

   int initialX = k - 1, finalX = numCtrlPts, numIterations = 0;
   float initial = knots[initialX], final = knots[finalX];
   float du = (final - initial) / (double) resolution;
  
   for(float u = initial; u < final; u += du) {
     /* find where u lies */
     int I = initialX;
     for(; I <= finalX - 1; I++) {
       if(u >= knots[I] && u < knots[I + 1]) {
         break;
       }
     }

   /* initialize set of points used in this layer of the calculation */
   Point myD[k][numCtrlPts];
   for(int i = 0; i < numCtrlPts; i++) {
     myD[0][i] = ctrlPoints[i];
   }
  
   /* calculate spine points this "layer" */  
   for(int r = 1; r <= (k - 1); r++) {
     for(int i = I - (k - 1); i <= (I - r); i++) {
       float alpha = (knots[i + k] - u), alpha2 = (u - knots[i + r]);
       float bottem = (knots[i + k] - knots[i + r]);
       if(bottem == 0) {
         alpha = alpha2 = 0;
         bottem = 1;
       }
       myD[r][i] = myD[r - 1][i] * (alpha / bottem) + myD[r - 1][i + 1] * (alpha2 / bottem);
     }
   }

   /* set output B-Spline point */
   bsplinePoints[numIterations++] = myD[k - 1][I - (k - 1)];
  }
}

Check out the source code.