1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
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;
}
|