#ifndef lint static char *rcs = "\$Header: spline.c,v 1.1 88/01/15 13:05:27 simpson Rel \$"; #endif /* * \$Log: spline.c,v \$ * Revision 1.1 88/01/15 13:05:27 simpson * initial release * * Revision 0.1 87/12/18 11:20:40 simpson * beta test * */ #include static double U[50]; /* Returns a set of (m+1) spline curve points given a set of (n+1) control * points. (t-1) is the degree of the polynomial used for the B-splines. * Cubic polynomials (i.e., t=4) are usually sufficient. [3] designates * (x,y,z). For a discussion of the spline curve algorithm, see the book * ``Computer Graphics'', Donald Hearn and M. Pauline Baker, Prentice-Hall, * 1986, pp. 200-2. */ void spline(controlpoints, n, curvepoints, m, t) double controlpoints[][3]; int n; double curvepoints[][3]; int m; int t; { int j, k; double N(), temp, u; for (j = 0; j <= n + t; j++) { if (j < t) U[j] = 0.0; else if (t <= j && j <= n) U[j] = j - t + 1; else U[j] = n - t + 2; } for (j = 0; j <= m; j++) { u = ((double)j / m)*(n-t+2-.00000001); curvepoints[j][0] = curvepoints[j][1] = curvepoints[j][2] = 0.0; for (k = 0; k <= n; k++) { temp = N(k, t, u); curvepoints[j][0] += controlpoints[k][0] * temp; curvepoints[j][1] += controlpoints[k][1] * temp; curvepoints[j][2] += controlpoints[k][2] * temp; } } } static double N(k, t, u) int k; int t; double u; { double firstterm, secondterm; if (t == 1) if (U[k] <= u && u < U[k+1]) return 1.0; else return 0.0; if (U[k+t-1]-U[k] < 1.0e-10) /* Test for zero with real #s */ firstterm = 0.0; else firstterm = ((u-U[k])/(U[k+t-1]-U[k]))*N(k,t-1,u); if (U[k+t]-U[k+1] < 1.0e-10) secondterm = 0.0; else secondterm = ((U[k+t]-u)/(U[k+t]-U[k+1]))*N(k+1,t-1,u); return firstterm + secondterm; }