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[];
{
        extern long     eclock;
        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(argv[++i], NULL, 0);
                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(stmp, NULL, 0);
        sprintf(stmp, "%s(s,t)=%s;", YNAME, argv[4]);
        scompile(stmp, NULL, 0);
        sprintf(stmp, "%s(s,t)=%s;", ZNAME, argv[5]);
        scompile(stmp, NULL, 0);
        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);
        eclock = 0;
                                                /* 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];
        int  end;
        register int  i;
        
        if (smooth) {
                i = -1;                 /* compute one past each end */
                end = siz+1;
        } else {
                if (s < -FTINY || s > 1.0+FTINY)
                        return;
                i = 0;
                end = siz;
        }
        st[0] = s;
        while (i <= end) {
                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);
                i++;
        }
}


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

        if (!smooth)                    /* not needed if no smoothing */
                return;
                                        /* compute middle points */
        while (siz-- >= 0) {
                fvsum(v1, r2[0].p, r0[0].p, -1.0);
                fvsum(v2, r1[1].p, r1[-1].p, -1.0);
                fcross(r1[0].n, v1, v2);
                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;
        MAT4  eqnmat;
        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)
MAT4  inverse, mat;
{
#define SWAP(a,b,t) (t=a,a=b,b=t)

        MAT4  m4tmp;
        register int i,j,k;
        register double temp;

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