src/gen/gensurf.c

#ifndef lint
static char SCCSid[] = "$SunId$ LBL";
#endif

/* Copyright (c) 1989 Regents of the University of California */

/*
 *  gensurf.c - program to generate functional surfaces
 *
 *      Parametric functions x(s,t), y(s,t) and z(s,t)
 *  specify the surface, which is tesselated into an m by n
 *  array of paired triangles.
 *      The surface normal is defined by the right hand
 *  rule applied to (s,t).
 *
 *      4/3/87
 */

#include  "standard.h"

#define  XNAME          "X_"                    /* x function name */
#define  YNAME          "Y_"                    /* y function name */
#define  ZNAME          "Z_"                    /* z function name */

#define  ABS(x)         ((x)>=0 ? (x) : -(x))

#define  pvect(p)       printf(vformat, (p)[0], (p)[1], (p)[2])

char  vformat[] = "%15.9g %15.9g %15.9g\n";
char  tsargs[] = "4 surf_dx surf_dy surf_dz surf.cal\n";
char  texname[] = "Phong";

int  smooth = 0;                /* apply smoothing? */

char  *modname, *surfname;

double  funvalue(), l_hermite(), l_bezier(), l_bspline(), argument();

typedef struct {
        FVECT  p;       /* vertex position */
        FVECT  n;       /* average normal */
} POINT;


main(argc, argv)
int  argc;
char  *argv[];
{
        POINT  *row0, *row1, *row2, *rp;
        int  i, j, m, n;
        char  stmp[256];

        varset("PI", PI);
        funset("hermite", 5, l_hermite);
        funset("bezier", 5, l_bezier);
        funset("bspline", 5, l_bspline);

        if (argc < 8)
                goto userror;

        for (i = 8; i < argc; i++)
                if (!strcmp(argv[i], "-e"))
                        scompile(NULL, argv[++i]);
                else if (!strcmp(argv[i], "-f"))
                        fcompile(argv[++i]);
                else if (!strcmp(argv[i], "-s"))
                        smooth++;
                else
                        goto userror;

        modname = argv[1];
        surfname = argv[2];
        sprintf(stmp, "%s(s,t)=%s;", XNAME, argv[3]);
        scompile(NULL, stmp);
        sprintf(stmp, "%s(s,t)=%s;", YNAME, argv[4]);
        scompile(NULL, stmp);
        sprintf(stmp, "%s(s,t)=%s;", ZNAME, argv[5]);
        scompile(NULL, stmp);
        m = atoi(argv[6]);
        n = atoi(argv[7]);
        if (m <= 0 || n <= 0)
                goto userror;

        row0 = (POINT *)malloc((n+3)*sizeof(POINT));
        row1 = (POINT *)malloc((n+3)*sizeof(POINT));
        row2 = (POINT *)malloc((n+3)*sizeof(POINT));
        if (row0 == NULL || row1 == NULL || row2 == NULL) {
                fprintf(stderr, "%s: out of memory\n", argv[0]);
                quit(1);
        }
        row0++; row1++; row2++;
                                                /* print header */
        printhead(argc, argv);
                                                /* initialize */
        comprow(-1.0/m, row0, n);
        comprow(0.0, row1, n);
        comprow(1.0/m, row2, n);
        compnorms(row0, row1, row2, n);
                                                /* for each row */
        for (i = 0; i < m; i++) {
                                                /* compute next row */
                rp = row0;
                row0 = row1;
                row1 = row2;
                row2 = rp;
                comprow((double)(i+2)/m, row2, n);
                compnorms(row0, row1, row2, n);

                for (j = 0; j < n; j++) {
                                                        /* put polygons */
                        if ((i+j) & 1)
                                putsquare(&row0[j], &row1[j],
                                                &row0[j+1], &row1[j+1]);
                        else
                                putsquare(&row1[j], &row1[j+1],
                                                &row0[j], &row0[j+1]);
                }
        }

        quit(0);

userror:
        fprintf(stderr, "Usage: %s material name ", argv[0]);
        fprintf(stderr, "x(s,t) y(s,t) z(s,t) m n [-s][-e expr][-f file]\n");
        quit(1);
}


putsquare(p0, p1, p2, p3)               /* put out a square */
POINT  *p0, *p1, *p2, *p3;
{
        static int  nout = 0;
        FVECT  norm[4];
        int  axis;
        FVECT  v1, v2, vc1, vc2;
        int  ok1, ok2;
                                        /* compute exact normals */
        fvsum(v1, p1->p, p0->p, -1.0);
        fvsum(v2, p2->p, p0->p, -1.0);
        fcross(vc1, v1, v2);
        ok1 = normalize(vc1) != 0.0;
        fvsum(v1, p2->p, p3->p, -1.0);
        fvsum(v2, p1->p, p3->p, -1.0);
        fcross(vc2, v1, v2);
        ok2 = normalize(vc2) != 0.0;
        if (!(ok1 | ok2))
                return;
                                        /* compute normal interpolation */
        axis = norminterp(norm, p0, p1, p2, p3);

                                        /* put out quadrilateral? */
        if (ok1 & ok2 && fdot(vc1,vc2) >= 1.0-FTINY*FTINY) {
                printf("\n%s ", modname);
                if (axis != -1) {
                        printf("texfunc %s\n", texname);
                        printf(tsargs);
                        printf("0\n13\t%d\n", axis);
                        pvect(norm[0]);
                        pvect(norm[1]);
                        pvect(norm[2]);
                        fvsum(v1, norm[3], vc1, -0.5);
                        fvsum(v1, v1, vc2, -0.5);
                        pvect(v1);
                        printf("\n%s ", texname);
                }
                printf("polygon %s.%d\n", surfname, ++nout);
                printf("0\n0\n12\n");
                pvect(p0->p);
                pvect(p1->p);
                pvect(p3->p);
                pvect(p2->p);
                return;
        }
                                        /* put out triangles? */
        if (ok1) {
                printf("\n%s ", modname);
                if (axis != -1) {
                        printf("texfunc %s\n", texname);
                        printf(tsargs);
                        printf("0\n13\t%d\n", axis);
                        pvect(norm[0]);
                        pvect(norm[1]);
                        pvect(norm[2]);
                        fvsum(v1, norm[3], vc1, -1.0);
                        pvect(v1);
                        printf("\n%s ", texname);
                }
                printf("polygon %s.%d\n", surfname, ++nout);
                printf("0\n0\n9\n");
                pvect(p0->p);
                pvect(p1->p);
                pvect(p2->p);
        }
        if (ok2) {
                printf("\n%s ", modname);
                if (axis != -1) {
                        printf("texfunc %s\n", texname);
                        printf(tsargs);
                        printf("0\n13\t%d\n", axis);
                        pvect(norm[0]);
                        pvect(norm[1]);
                        pvect(norm[2]);
                        fvsum(v2, norm[3], vc2, -1.0);
                        pvect(v2);
                        printf("\n%s ", texname);
                }
                printf("polygon %s.%d\n", surfname, ++nout);
                printf("0\n0\n9\n");
                pvect(p2->p);
                pvect(p1->p);
                pvect(p3->p);
        }
}


comprow(s, row, siz)                    /* compute row of values */
double  s;
register POINT  *row;
int  siz;
{
        double  st[2];
        register int  i;
                                        /* compute one past each end */
        st[0] = s;
        for (i = -1; i <= siz+1; i++) {
                st[1] = (double)i/siz;
                row[i].p[0] = funvalue(XNAME, 2, st);
                row[i].p[1] = funvalue(YNAME, 2, st);
                row[i].p[2] = funvalue(ZNAME, 2, st);
        }
}


compnorms(r0, r1, r2, siz)              /* compute row of averaged normals */
register POINT  *r0, *r1, *r2;
int  siz;
{
        FVECT  v1, v2, vc;
        register int  i;

        if (!smooth)                    /* not needed if no smoothing */
                return;
                                        /* compute middle points */
        while (siz-- >= 0) {
                fvsum(v1, r2[0].p, r1[0].p, -1.0);
                fvsum(v2, r1[1].p, r1[0].p, -1.0);
                fcross(r1[0].n, v1, v2);
                fvsum(v1, r0[0].p, r1[0].p, -1.0);
                fcross(vc, v2, v1);
                fvsum(r1[0].n, r1[0].n, vc, 1.0);
                fvsum(v2, r1[-1].p, r1[0].p, -1.0);
                fcross(vc, v1, v2);
                fvsum(r1[0].n, r1[0].n, vc, 1.0);
                fvsum(v1, r2[0].p, r1[0].p, -1.0);
                fcross(vc, v2, v1);
                fvsum(r1[0].n, r1[0].n, vc, 1.0);
                normalize(r1[0].n);
                r0++; r1++; r2++;
        }
}


int
norminterp(resmat, p0, p1, p2, p3)      /* compute normal interpolation */
register FVECT  resmat[4];
POINT  *p0, *p1, *p2, *p3;
{
#define u  ((ax+1)%3)
#define v  ((ax+2)%3)

        register int  ax;
        double  eqnmat[4][4];
        FVECT  v1;
        register int  i, j;

        if (!smooth)                    /* no interpolation if no smoothing */
                return(-1);
                                        /* find dominant axis */
        VCOPY(v1, p0->n);
        fvsum(v1, v1, p1->n, 1.0);
        fvsum(v1, v1, p2->n, 1.0);
        fvsum(v1, v1, p3->n, 1.0);
        ax = ABS(v1[0]) > ABS(v1[1]) ? 0 : 1;
        ax = ABS(v1[ax]) > ABS(v1[2]) ? ax : 2;
                                        /* assign equation matrix */
        eqnmat[0][0] = p0->p[u]*p0->p[v];
        eqnmat[0][1] = p0->p[u];
        eqnmat[0][2] = p0->p[v];
        eqnmat[0][3] = 1.0;
        eqnmat[1][0] = p1->p[u]*p1->p[v];
        eqnmat[1][1] = p1->p[u];
        eqnmat[1][2] = p1->p[v];
        eqnmat[1][3] = 1.0;
        eqnmat[2][0] = p2->p[u]*p2->p[v];
        eqnmat[2][1] = p2->p[u];
        eqnmat[2][2] = p2->p[v];
        eqnmat[2][3] = 1.0;
        eqnmat[3][0] = p3->p[u]*p3->p[v];
        eqnmat[3][1] = p3->p[u];
        eqnmat[3][2] = p3->p[v];
        eqnmat[3][3] = 1.0;
                                        /* invert matrix (solve system) */
        if (!invmat(eqnmat, eqnmat))
                return(-1);                     /* no solution */
                                        /* compute result matrix */
        for (j = 0; j < 4; j++)
                for (i = 0; i < 3; i++)
                        resmat[j][i] =  eqnmat[j][0]*p0->n[i] +
                                        eqnmat[j][1]*p1->n[i] +
                                        eqnmat[j][2]*p2->n[i] +
                                        eqnmat[j][3]*p3->n[i];
        return(ax);

#undef u
#undef v
}


/*
 * invmat - computes the inverse of mat into inverse.  Returns 1
 * if there exists an inverse, 0 otherwise.  It uses Gaussian Elimination
 * method.
 */

invmat(inverse,mat)
double mat[4][4],inverse[4][4];
{
#define SWAP(a,b,t) (t=a,a=b,b=t)

        double  m4tmp[4][4];
        register int i,j,k;
        register double temp;

        bcopy((char *)mat, (char *)m4tmp, sizeof(m4tmp));
                                        /* set inverse to identity */
        for (i = 0; i < 4; i++)
                for (j = 0; j < 4; j++)
                        inverse[i][j] = i==j ? 1.0 : 0.0;

        for(i = 0; i < 4; i++) {
                /* Look for raw with largest pivot and swap raws */
                temp = FTINY; j = -1;
                for(k = i; k < 4; k++)
                        if(ABS(m4tmp[k][i]) > temp) {
                                temp = ABS(m4tmp[k][i]);
                                j = k;
                                }
                if(j == -1)     /* No replacing raw -> no inverse */
                        return(0);
                if (j != i)
                        for(k = 0; k < 4; k++) {
                                SWAP(m4tmp[i][k],m4tmp[j][k],temp);
                                SWAP(inverse[i][k],inverse[j][k],temp);
                                }

                temp = m4tmp[i][i];
                for(k = 0; k < 4; k++) {
                        m4tmp[i][k] /= temp;
                        inverse[i][k] /= temp;
                        }
                for(j = 0; j < 4; j++) {
                        if(j != i) {
                                temp = m4tmp[j][i];
                                for(k = 0; k < 4; k++) {
                                        m4tmp[j][k] -= m4tmp[i][k]*temp;
                                        inverse[j][k] -= inverse[i][k]*temp;
                                        }
                                }
                        }
                }
        return(1);

#undef SWAP
}


eputs(msg)
char  *msg;
{
        fputs(msg, stderr);
}


wputs(msg)
char  *msg;
{
        eputs(msg);
}


quit(code)
{
        exit(code);
}


printhead(ac, av)               /* print command header */
register int  ac;
register char  **av;
{
        putchar('#');
        while (ac--) {
                putchar(' ');
                fputs(*av++, stdout);
        }
        putchar('\n');
}


double
l_hermite()                     
{
        double  t;
        
        t = argument(5);
        return( argument(1)*((2.0*t-3.0)*t*t+1.0) +
                argument(2)*(-2.0*t+3.0)*t*t +
                argument(3)*((t-2.0)*t+1.0)*t +
                argument(4)*(t-1.0)*t*t );
}


double
l_bezier()
{
        double  t;

        t = argument(5);
        return( argument(1) * (1.+t*(-3.+t*(3.-t))) +
                argument(2) * 3.*t*(1.+t*(-2.+t)) +
                argument(3) * 3.*t*t*(1.-t) +
                argument(4) * t*t*t );
}


double
l_bspline()
{
        double  t;

        t = argument(5);
        return( argument(1) * (1./6.+t*(-1./2.+t*(1./2.-1./6.*t))) +
                argument(2) * (2./3.+t*t*(-1.+1./2.*t)) +
                argument(3) * (1./6.+t*(1./2.+t*(1./2.-1./2.*t))) +
                argument(4) * (1./6.*t*t*t) );
}
Revision:	1.7
Committed:	Wed Mar 7 11:14:36 1990 UTC (34 years, 1 month ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 1.6:	+17 -3 lines
Log Message:	added bezier and bspline functions
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2			static char SCCSid[] = "$SunId$ LBL";
3			#endif
4	greg	1.2
5	greg	1.7	/* Copyright (c) 1989 Regents of the University of California */
6
7	greg	1.2	/*
8	greg	1.1	* gensurf.c - program to generate functional surfaces
9			*
10			* Parametric functions x(s,t), y(s,t) and z(s,t)
11			* specify the surface, which is tesselated into an m by n
12			* array of paired triangles.
13			* The surface normal is defined by the right hand
14			* rule applied to (s,t).
15			*
16			* 4/3/87
17			*/
18
19	greg	1.5	#include "standard.h"
20	greg	1.1
21			#define XNAME "X_" /* x function name */
22			#define YNAME "Y_" /* y function name */
23			#define ZNAME "Z_" /* z function name */
24
25	greg	1.4	#define ABS(x) ((x)>=0 ? (x) : -(x))
26
27	greg	1.3	#define pvect(p) printf(vformat, (p)[0], (p)[1], (p)[2])
28	greg	1.1
29			char vformat[] = "%15.9g %15.9g %15.9g\n";
30	greg	1.3	char tsargs[] = "4 surf_dx surf_dy surf_dz surf.cal\n";
31			char texname[] = "Phong";
32	greg	1.1
33	greg	1.3	int smooth = 0; /* apply smoothing? */
34	greg	1.1
35	greg	1.3	char modname, surfname;
36	greg	1.1
37	greg	1.7	double funvalue(), l_hermite(), l_bezier(), l_bspline(), argument();
38	greg	1.3
39			typedef struct {
40			FVECT p; /* vertex position */
41			FVECT n; /* average normal */
42			} POINT;
43
44
45	greg	1.1	main(argc, argv)
46			int argc;
47			char *argv[];
48			{
49	greg	1.3	POINT row0, row1, row2, rp;
50	greg	1.1	int i, j, m, n;
51			char stmp[256];
52
53			varset("PI", PI);
54			funset("hermite", 5, l_hermite);
55	greg	1.6	funset("bezier", 5, l_bezier);
56	greg	1.7	funset("bspline", 5, l_bspline);
57	greg	1.1
58			if (argc < 8)
59			goto userror;
60
61			for (i = 8; i < argc; i++)
62			if (!strcmp(argv[i], "-e"))
63			scompile(NULL, argv[++i]);
64			else if (!strcmp(argv[i], "-f"))
65			fcompile(argv[++i]);
66	greg	1.3	else if (!strcmp(argv[i], "-s"))
67			smooth++;
68	greg	1.1	else
69			goto userror;
70
71	greg	1.3	modname = argv[1];
72			surfname = argv[2];
73	greg	1.1	sprintf(stmp, "%s(s,t)=%s;", XNAME, argv[3]);
74			scompile(NULL, stmp);
75			sprintf(stmp, "%s(s,t)=%s;", YNAME, argv[4]);
76			scompile(NULL, stmp);
77			sprintf(stmp, "%s(s,t)=%s;", ZNAME, argv[5]);
78			scompile(NULL, stmp);
79			m = atoi(argv[6]);
80			n = atoi(argv[7]);
81			if (m <= 0 \|\| n <= 0)
82			goto userror;
83
84	greg	1.4	row0 = (POINT )malloc((n+3)sizeof(POINT));
85			row1 = (POINT )malloc((n+3)sizeof(POINT));
86			row2 = (POINT )malloc((n+3)sizeof(POINT));
87	greg	1.3	if (row0 == NULL \|\| row1 == NULL \|\| row2 == NULL) {
88	greg	1.1	fprintf(stderr, "%s: out of memory\n", argv[0]);
89			quit(1);
90			}
91	greg	1.4	row0++; row1++; row2++;
92	greg	1.3	/* print header */
93	greg	1.1	printhead(argc, argv);
94	greg	1.4	/* initialize */
95			comprow(-1.0/m, row0, n);
96	greg	1.3	comprow(0.0, row1, n);
97			comprow(1.0/m, row2, n);
98	greg	1.4	compnorms(row0, row1, row2, n);
99	greg	1.3	/* for each row */
100	greg	1.1	for (i = 0; i < m; i++) {
101			/* compute next row */
102	greg	1.3	rp = row0;
103	greg	1.1	row0 = row1;
104	greg	1.3	row1 = row2;
105			row2 = rp;
106	greg	1.4	comprow((double)(i+2)/m, row2, n);
107			compnorms(row0, row1, row2, n);
108	greg	1.1
109			for (j = 0; j < n; j++) {
110	greg	1.3	/* put polygons */
111			if ((i+j) & 1)
112			putsquare(&row0[j], &row1[j],
113			&row0[j+1], &row1[j+1]);
114			else
115			putsquare(&row1[j], &row1[j+1],
116			&row0[j], &row0[j+1]);
117	greg	1.1	}
118			}
119
120			quit(0);
121
122			userror:
123			fprintf(stderr, "Usage: %s material name ", argv[0]);
124	greg	1.3	fprintf(stderr, "x(s,t) y(s,t) z(s,t) m n [-s][-e expr][-f file]\n");
125	greg	1.1	quit(1);
126			}
127
128
129	greg	1.3	putsquare(p0, p1, p2, p3) /* put out a square */
130			POINT p0, p1, p2, p3;
131			{
132			static int nout = 0;
133			FVECT norm[4];
134			int axis;
135			FVECT v1, v2, vc1, vc2;
136			int ok1, ok2;
137			/* compute exact normals */
138			fvsum(v1, p1->p, p0->p, -1.0);
139			fvsum(v2, p2->p, p0->p, -1.0);
140			fcross(vc1, v1, v2);
141			ok1 = normalize(vc1) != 0.0;
142			fvsum(v1, p2->p, p3->p, -1.0);
143			fvsum(v2, p1->p, p3->p, -1.0);
144			fcross(vc2, v1, v2);
145			ok2 = normalize(vc2) != 0.0;
146			if (!(ok1 \| ok2))
147			return;
148			/* compute normal interpolation */
149			axis = norminterp(norm, p0, p1, p2, p3);
150
151			/* put out quadrilateral? */
152			if (ok1 & ok2 && fdot(vc1,vc2) >= 1.0-FTINY*FTINY) {
153			printf("\n%s ", modname);
154			if (axis != -1) {
155			printf("texfunc %s\n", texname);
156			printf(tsargs);
157			printf("0\n13\t%d\n", axis);
158			pvect(norm[0]);
159			pvect(norm[1]);
160			pvect(norm[2]);
161			fvsum(v1, norm[3], vc1, -0.5);
162			fvsum(v1, v1, vc2, -0.5);
163			pvect(v1);
164			printf("\n%s ", texname);
165			}
166			printf("polygon %s.%d\n", surfname, ++nout);
167			printf("0\n0\n12\n");
168			pvect(p0->p);
169			pvect(p1->p);
170			pvect(p3->p);
171			pvect(p2->p);
172			return;
173			}
174			/* put out triangles? */
175			if (ok1) {
176			printf("\n%s ", modname);
177			if (axis != -1) {
178			printf("texfunc %s\n", texname);
179			printf(tsargs);
180			printf("0\n13\t%d\n", axis);
181			pvect(norm[0]);
182			pvect(norm[1]);
183			pvect(norm[2]);
184			fvsum(v1, norm[3], vc1, -1.0);
185			pvect(v1);
186			printf("\n%s ", texname);
187			}
188			printf("polygon %s.%d\n", surfname, ++nout);
189			printf("0\n0\n9\n");
190			pvect(p0->p);
191			pvect(p1->p);
192			pvect(p2->p);
193			}
194			if (ok2) {
195			printf("\n%s ", modname);
196			if (axis != -1) {
197			printf("texfunc %s\n", texname);
198			printf(tsargs);
199			printf("0\n13\t%d\n", axis);
200			pvect(norm[0]);
201			pvect(norm[1]);
202			pvect(norm[2]);
203			fvsum(v2, norm[3], vc2, -1.0);
204			pvect(v2);
205			printf("\n%s ", texname);
206			}
207			printf("polygon %s.%d\n", surfname, ++nout);
208			printf("0\n0\n9\n");
209			pvect(p2->p);
210			pvect(p1->p);
211			pvect(p3->p);
212			}
213			}
214
215
216	greg	1.1	comprow(s, row, siz) /* compute row of values */
217			double s;
218	greg	1.3	register POINT *row;
219	greg	1.1	int siz;
220			{
221	greg	1.4	double st[2];
222			register int i;
223			/* compute one past each end */
224	greg	1.1	st[0] = s;
225	greg	1.4	for (i = -1; i <= siz+1; i++) {
226			st[1] = (double)i/siz;
227			row[i].p[0] = funvalue(XNAME, 2, st);
228			row[i].p[1] = funvalue(YNAME, 2, st);
229			row[i].p[2] = funvalue(ZNAME, 2, st);
230	greg	1.1	}
231	greg	1.3	}
232
233
234			compnorms(r0, r1, r2, siz) /* compute row of averaged normals */
235			register POINT r0, r1, *r2;
236			int siz;
237			{
238			FVECT v1, v2, vc;
239	greg	1.4	register int i;
240	greg	1.3
241			if (!smooth) /* not needed if no smoothing */
242			return;
243			/* compute middle points */
244	greg	1.4	while (siz-- >= 0) {
245	greg	1.3	fvsum(v1, r2[0].p, r1[0].p, -1.0);
246			fvsum(v2, r1[1].p, r1[0].p, -1.0);
247			fcross(r1[0].n, v1, v2);
248			fvsum(v1, r0[0].p, r1[0].p, -1.0);
249			fcross(vc, v2, v1);
250			fvsum(r1[0].n, r1[0].n, vc, 1.0);
251			fvsum(v2, r1[-1].p, r1[0].p, -1.0);
252			fcross(vc, v1, v2);
253			fvsum(r1[0].n, r1[0].n, vc, 1.0);
254			fvsum(v1, r2[0].p, r1[0].p, -1.0);
255			fcross(vc, v2, v1);
256			fvsum(r1[0].n, r1[0].n, vc, 1.0);
257			normalize(r1[0].n);
258			r0++; r1++; r2++;
259			}
260			}
261
262
263			int
264			norminterp(resmat, p0, p1, p2, p3) /* compute normal interpolation */
265			register FVECT resmat[4];
266			POINT p0, p1, p2, p3;
267			{
268			#define u ((ax+1)%3)
269			#define v ((ax+2)%3)
270
271			register int ax;
272	greg	1.4	double eqnmat[4][4];
273	greg	1.3	FVECT v1;
274			register int i, j;
275
276			if (!smooth) /* no interpolation if no smoothing */
277			return(-1);
278			/* find dominant axis */
279			VCOPY(v1, p0->n);
280			fvsum(v1, v1, p1->n, 1.0);
281			fvsum(v1, v1, p2->n, 1.0);
282			fvsum(v1, v1, p3->n, 1.0);
283	greg	1.4	ax = ABS(v1[0]) > ABS(v1[1]) ? 0 : 1;
284			ax = ABS(v1[ax]) > ABS(v1[2]) ? ax : 2;
285	greg	1.3	/* assign equation matrix */
286			eqnmat[0][0] = p0->p[u]*p0->p[v];
287			eqnmat[0][1] = p0->p[u];
288			eqnmat[0][2] = p0->p[v];
289			eqnmat[0][3] = 1.0;
290			eqnmat[1][0] = p1->p[u]*p1->p[v];
291			eqnmat[1][1] = p1->p[u];
292			eqnmat[1][2] = p1->p[v];
293			eqnmat[1][3] = 1.0;
294			eqnmat[2][0] = p2->p[u]*p2->p[v];
295			eqnmat[2][1] = p2->p[u];
296			eqnmat[2][2] = p2->p[v];
297			eqnmat[2][3] = 1.0;
298			eqnmat[3][0] = p3->p[u]*p3->p[v];
299			eqnmat[3][1] = p3->p[u];
300			eqnmat[3][2] = p3->p[v];
301			eqnmat[3][3] = 1.0;
302			/* invert matrix (solve system) */
303	greg	1.4	if (!invmat(eqnmat, eqnmat))
304	greg	1.3	return(-1); /* no solution */
305			/* compute result matrix */
306			for (j = 0; j < 4; j++)
307			for (i = 0; i < 3; i++)
308	greg	1.4	resmat[j][i] = eqnmat[j][0]*p0->n[i] +
309			eqnmat[j][1]*p1->n[i] +
310			eqnmat[j][2]*p2->n[i] +
311			eqnmat[j][3]*p3->n[i];
312	greg	1.3	return(ax);
313
314			#undef u
315			#undef v
316			}
317
318
319			/*
320			* invmat - computes the inverse of mat into inverse. Returns 1
321			* if there exists an inverse, 0 otherwise. It uses Gaussian Elimination
322			* method.
323			*/
324
325			invmat(inverse,mat)
326			double mat[4][4],inverse[4][4];
327			{
328			#define SWAP(a,b,t) (t=a,a=b,b=t)
329
330	greg	1.4	double m4tmp[4][4];
331	greg	1.3	register int i,j,k;
332			register double temp;
333
334	greg	1.5	bcopy((char )mat, (char )m4tmp, sizeof(m4tmp));
335	greg	1.4	/* set inverse to identity */
336			for (i = 0; i < 4; i++)
337			for (j = 0; j < 4; j++)
338			inverse[i][j] = i==j ? 1.0 : 0.0;
339	greg	1.3
340			for(i = 0; i < 4; i++) {
341	greg	1.4	/* Look for raw with largest pivot and swap raws */
342			temp = FTINY; j = -1;
343			for(k = i; k < 4; k++)
344			if(ABS(m4tmp[k][i]) > temp) {
345			temp = ABS(m4tmp[k][i]);
346			j = k;
347			}
348			if(j == -1) /* No replacing raw -> no inverse */
349			return(0);
350			if (j != i)
351			for(k = 0; k < 4; k++) {
352			SWAP(m4tmp[i][k],m4tmp[j][k],temp);
353			SWAP(inverse[i][k],inverse[j][k],temp);
354			}
355	greg	1.3
356			temp = m4tmp[i][i];
357			for(k = 0; k < 4; k++) {
358			m4tmp[i][k] /= temp;
359			inverse[i][k] /= temp;
360			}
361			for(j = 0; j < 4; j++) {
362			if(j != i) {
363			temp = m4tmp[j][i];
364			for(k = 0; k < 4; k++) {
365			m4tmp[j][k] -= m4tmp[i][k]*temp;
366			inverse[j][k] -= inverse[i][k]*temp;
367			}
368			}
369			}
370			}
371			return(1);
372	greg	1.4
373	greg	1.3	#undef SWAP
374	greg	1.1	}
375
376
377			eputs(msg)
378			char *msg;
379			{
380			fputs(msg, stderr);
381			}
382
383
384			wputs(msg)
385			char *msg;
386			{
387			eputs(msg);
388			}
389
390
391			quit(code)
392			{
393			exit(code);
394			}
395
396
397			printhead(ac, av) /* print command header */
398			register int ac;
399			register char **av;
400			{
401			putchar('#');
402			while (ac--) {
403			putchar(' ');
404			fputs(*av++, stdout);
405			}
406			putchar('\n');
407			}
408
409
410			double
411			l_hermite()
412			{
413			double t;
414
415			t = argument(5);
416			return( argument(1)((2.0t-3.0)tt+1.0) +
417			argument(2)(-2.0t+3.0)tt +
418			argument(3)((t-2.0)t+1.0)*t +
419			argument(4)(t-1.0)t*t );
420	greg	1.6	}
421
422
423			double
424			l_bezier()
425			{
426			double t;
427
428			t = argument(5);
429			return( argument(1) * (1.+t(-3.+t(3.-t))) +
430			argument(2) * 3.t(1.+t*(-2.+t)) +
431			argument(3) * 3.tt*(1.-t) +
432			argument(4) * ttt );
433	greg	1.7	}
434
435
436			double
437			l_bspline()
438			{
439			double t;
440
441			t = argument(5);
442			return( argument(1) * (1./6.+t(-1./2.+t(1./2.-1./6.*t))) +
443			argument(2) * (2./3.+tt(-1.+1./2.*t)) +
444			argument(3) * (1./6.+t(1./2.+t(1./2.-1./2.*t))) +
445			argument(4) * (1./6.tt*t) );
446	greg	1.1	}