diff options
Diffstat (limited to 'dviware/quicspool/src/spline.c')
-rw-r--r-- | dviware/quicspool/src/spline.c | 75 |
1 files changed, 75 insertions, 0 deletions
diff --git a/dviware/quicspool/src/spline.c b/dviware/quicspool/src/spline.c new file mode 100644 index 0000000000..73c319463c --- /dev/null +++ b/dviware/quicspool/src/spline.c @@ -0,0 +1,75 @@ +#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 <math.h> + +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; +} |