Spline evaluation

gistfile1.cpp

float safeDiv(float a, float b) {
  if (abs(a) < 1e-10 && abs(b) < 1e-10) return 0;
  else return a / b;
}

float b(vector<float> kn, int i, int n, float t) {
  if (n == 0) {
    return (kn[i] <= t < kn[i+1]) ? 1 : 0;
  }
  return 
    safeDiv(t - kn[i], kn[i+n] - kn[i]) * b(kn, i, n-1, t) +
    safeDiv(kn[i+n+1] - t, kn[i+n+1] - kn[i+1]) * b(kn, i+1, n-1, t);
}

point evaluateCurve(curveOrSurface geometry, float t) {
  point result(0, 0, 0);
  int L = geometry.knotsU.size() - 1;

  for (int j = 0; j <= L; j++) {
    float b_j = b(geometry.knotsU, j, geometry.degreeU, t);
    point P_j = geometry.controlPoints[j];

    result = addpts(result, scalarmul(b_j, P_j));
  }

  return result;
}

void drawCurve(curveOrSurface &geometry) {
  assert(geometry.isCurve());

  const int sample_count = geometry.controlPoints.size() * 4;
  assert(sample_count > 1);

  if (geometry.evaluatedPoints.size() != sample_count) {
    geometry.evaluatedPoints.resize(sample_count);

    float max_knot = 0;
    for (int i = 0; i < geometry.knotsU.size(); i++) {
      if (max_knot < geometry.knotsU[i]) {
        max_knot = geometry.knotsU[i];
      }
    }

    for (int i = 0; i < sample_count; i++) {
      float t = max_knot * (float) i / (float) (sample_count - 1.0);
      geometry.evaluatedPoints[i] = evaluateCurve(geometry, t);
    }
  }

  glBegin(GL_LINE_STRIP); for (int i = 0; i < sample_count; i++) {
    printPt(geometry.evaluatedPoints[i]); cout << endl;
    glVertex3f(
        geometry.evaluatedPoints[i].x(),
        geometry.evaluatedPoints[i].y(),
        geometry.evaluatedPoints[i].z());
  } glEnd();
}