To learn a little about how NURBS curves work I decided to implement a NURBS class in C#. Following several different resources online I have created a basic class. Much of its functionality I have tested and verified the behaviours of. However, when evaluating points on the curves I find myself getting into situations where a singularity occurs within the NURBS basis function. I have implemented the basis function as per the algorithm below.
My implementation of this function can be seen here in C#.
// i is the knot index, filled by iterating over control points which are less than the number of knots
// p is the degree
// x is the position to evaluate at
private double N (int i, int p, double x) {
double ti0 = knotVector[i]; //u_{i}
double ti1 = knotVector[i+1]; //u_{i+1}
double tip = knotVector[i+p]; // u_{i+p}
double tip1 = knotVector[i+p+1];// u_{i+p+1}
if (p == 0) {
if (ti0 <= x && x < ti1) {
return 1;
} else {
return 0;
}
}
double fi0 = (x - ti0) / (tip - ti0);
double gi1 = (tip1 - x) / (tip1 - ti1);
return fi0 * N(i, p-1, x)
+ gi1 * N(i+1, p-1, x);
}
Where this basis function is used in the property below to evaluate the coordinates of the point in 3d space.
public Vec3 this[double t] {
get {
double x = 0, y = 0, z = 0;
double rationalWeight = 0;
double[] ns = new double[this.ControlPoints.Length];
for (int i = 0; i < this.ControlPoints.Length; i++) {
double n = N(i, this.Degree, t);
ns[i] = n;
double temp = n * this.ControlPoints[i].Weight;
rationalWeight += temp;
}
for (int i = 0; i < this.ControlPoints.Length; i++) {
double temp = ns[i];
x += this.ControlPoints[i].X * this.ControlPoints[i].Weight * temp / rationalWeight;
y += this.ControlPoints[i].Y * this.ControlPoints[i].Weight * temp / rationalWeight;
z += this.ControlPoints[i].Z * this.ControlPoints[i].Weight * temp / rationalWeight;
}
return new Vec3(x,y,z);
}
}
For instance, for this curve
- degree 3
- points [(-4, -4, 0), (-2, 4, 0), (2, -4, 0), (4, 4, 0)]
- weights [1, 1, 1, 1]
- knots [0,0,0,0,1,1,1,1]
I end up getting a cases, such as at x = 0.5, where (tip - ti0) and (tip1 - ti1) are 0 resulting in divide by 0 and NaN issues. Instead of getting a correct value (as verified from the curve link earlier), I am simply getting "Not a Number" from the function. From my reading so far I have found nothing to describing dealing singularities. This indicates to me that I am missing something in my implementation or have misunderstood how parts of it work.
Does anyone have an idea where I have gone wrong with this implementation of the basis function?
