#include #include // Copyright YON - Jan C. Hardenbergh. Permission to copy is // granted as long as this copyright and this URL are maintained: // http://www.jch.com/NURBS/ // Primary reference Les Peigl - On NURBS: A Survey published // in IEEE Computer Graphics and Applications (CG&A) January 1991. // http://www.cs.wpi.edu/~matt/courses/cs563/talks/nurbs.html // this is intended to explore NURBS - how they are evaluated. // It is explicitly unoptimized. // the NURBCurve used as an example here is a circle inscribed // in an equalateral triangle. I found it on page 374, // Mathematical Elements for Computer Graphics, // 2nd ed. David F. Rogers and J Alan Adams. McGraw Hill, 1990. class Point { public: Point() {m_x = 0; m_y = 0;} Point(double x, double y) {m_x = x; m_y = y;} double m_x, m_y; Point *scale(double d) { return (new Point(m_x*d,m_y*d));} void add(Point &other) { m_x += other.m_x; m_y += other.m_y;} void set(double x, double y) {m_x = x; m_y = y;} }; double B(int i, int k, double t, double *knots) { double ret, a, b; if (k > 0) { double n1 = ((t - knots[i]) * B(i,k-1,t,knots)); double d1 = knots[i+k] - knots[i]; double n2 = ((knots[i+k+1] - t) * B(i+1,k-1,t,knots)); double d2 = knots[i+k+1] - knots[i+1]; if (fabs(d1) > 0.0001) a = n1 / d1; else a = 0; if (fabs(d2) > 0.0001) b = n2 / d2; else b = 0; ret = a + b; } //printf(" i = %d, k = %d, ret = %g, a = %g, b = %g\n",i,k,ret,a,b); else // j == 1 { if ((knots[i] <= t) && (t <= knots[i+1])) ret = 1; else ret = 0; } return (ret); } Point *C(double t, int order, int nPoints, Point *points, double *weights, double *knots) { Point *c = new Point(); double rational = 0; for (int i = 0; i < nPoints; i++) { double b = B(i, order, t, knots); Point *p = points[i].scale(b*weights[i]); c->add(*p); rational += b*weights[i]; } // printf("t = %g i = %d, b = %g, t.x = %g, t.y = %g, c.x = %g, c.y = %g, div = %g\n",t,i,b,p->m_x, p->m_y, c->m_x, c->m_y, div); return (c->scale(1.0/rational)); } #define RAD3 1.732 void main(int argc, char **argv) { Point controlPoints[7];// = {{0,0,1},{1,0,.5},{.5,.5*RAD3,1},{0,RAD3,.5},{-.5,.5*RAD3,1},{-1,0,.5},{0,0,1}} double knots[10] = {0,0,0,1,1,2,2,3,3,3}; double weights[7] = {1,.5,1,.5,1,.5,1}; int order = 2; int i; controlPoints[0].set(0,0); controlPoints[1].set(1,0); controlPoints[2].set(.5,.5*RAD3); controlPoints[3].set(0,RAD3); controlPoints[4].set(-.5,.5*RAD3); controlPoints[5].set(-1,0); controlPoints[6].set(0,0); // interesting point 1.5 is smack in the middle. // Point *c15 = C(1.5, order, 7, controlPoints, weights, knots); #define STEPS 15 for (i=0;i<=STEPS;i++) { double t = i*(3.0/STEPS); Point *c = C(t, order, 7, controlPoints, weights, knots); printf("C(%g) x = %g, y = %g\n", t, c->m_x, c->m_y); } }