summaryrefslogtreecommitdiff
path: root/dviware/quicspool/src/spline.c
diff options
context:
space:
mode:
Diffstat (limited to 'dviware/quicspool/src/spline.c')
-rw-r--r--dviware/quicspool/src/spline.c75
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;
+}