src/gen/gensurf.c

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

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

/*
 *  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  <stdio.h>
#include  "fvect.h"

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

#define  PI             3.14159265358979323846

#define  FTINY          1e-7

#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(), argument(), fabs();

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);

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

        st[0] = s;
        st[1] = 0.0;
        step = 1.0 / siz;
        while (siz-- >= 0) {
                row->p[0] = funvalue(XNAME, 2, st);
                row->p[1] = funvalue(YNAME, 2, st);
                row->p[2] = funvalue(ZNAME, 2, st);
                row++;
                st[1] += step;
        }
}


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

        if (!smooth)                    /* not needed if no smoothing */
                return;
                                        /* compute first point */
        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);
        normalize(r1[0].n);
        r0++; r1++; r2++;
                                        /* 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++;
        }
                                        /* compute end point */
        fvsum(v1, r0[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, 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);
}


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], solmat[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 = fabs(v1[0]) > fabs(v1[1]) ? 0 : 1;
        ax = fabs(v1[ax]) > fabs(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(solmat, eqnmat))
                return(-1);                     /* no solution */
                                        /* compute result matrix */
        for (j = 0; j < 4; j++)
                for (i = 0; i < 3; i++)
                        resmat[j][i] =  solmat[j][0]*p0->n[i] +
                                        solmat[j][1]*p1->n[i] +
                                        solmat[j][2]*p2->n[i] +
                                        solmat[j][3]*p3->n[i];
        return(ax);

#undef u
#undef v
}


static double  m4tmp[4][4];             /* for efficiency */

#define  copymat4(m4a,m4b)      bcopy((char *)m4b,(char *)m4a,sizeof(m4tmp))


setident4(m4)
double  m4[4][4];
{
        static double  ident[4][4] = {
                1.,0.,0.,0.,
                0.,1.,0.,0.,
                0.,0.,1.,0.,
                0.,0.,0.,1.,
        };
        copymat4(m4, ident);
}

/*
 * 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)

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

        setident4(inverse);
        copymat4(m4tmp, mat);

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

                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 );
}
Revision:	1.3
Committed:	Wed Oct 18 15:01:23 1989 UTC (36 years ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 1.2:	+307 -83 lines
Log Message:	added smoothing option
#	User	Rev	Content
1	greg	1.2	/* Copyright (c) 1989 Regents of the University of California */
2	greg	1.1
3			#ifndef lint
4			static char SCCSid[] = "$SunId$ LBL";
5			#endif
6	greg	1.2
7			/*
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			#include <stdio.h>
20	greg	1.3	#include "fvect.h"
21	greg	1.1
22			#define XNAME "X_" /* x function name */
23			#define YNAME "Y_" /* y function name */
24			#define ZNAME "Z_" /* z function name */
25
26			#define PI 3.14159265358979323846
27
28			#define FTINY 1e-7
29
30	greg	1.3	#define pvect(p) printf(vformat, (p)[0], (p)[1], (p)[2])
31	greg	1.1
32			char vformat[] = "%15.9g %15.9g %15.9g\n";
33	greg	1.3	char tsargs[] = "4 surf_dx surf_dy surf_dz surf.cal\n";
34			char texname[] = "Phong";
35	greg	1.1
36	greg	1.3	int smooth = 0; /* apply smoothing? */
37	greg	1.1
38	greg	1.3	char modname, surfname;
39	greg	1.1
40	greg	1.3	double funvalue(), l_hermite(), argument(), fabs();
41
42			typedef struct {
43			FVECT p; /* vertex position */
44			FVECT n; /* average normal */
45			} POINT;
46
47
48	greg	1.1	main(argc, argv)
49			int argc;
50			char *argv[];
51			{
52	greg	1.3	POINT row0, row1, row2, rp;
53	greg	1.1	int i, j, m, n;
54			char stmp[256];
55
56			varset("PI", PI);
57			funset("hermite", 5, l_hermite);
58
59			if (argc < 8)
60			goto userror;
61
62			for (i = 8; i < argc; i++)
63			if (!strcmp(argv[i], "-e"))
64			scompile(NULL, argv[++i]);
65			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			scompile(NULL, stmp);
76			sprintf(stmp, "%s(s,t)=%s;", YNAME, argv[4]);
77			scompile(NULL, stmp);
78			sprintf(stmp, "%s(s,t)=%s;", ZNAME, argv[5]);
79			scompile(NULL, stmp);
80			m = atoi(argv[6]);
81			n = atoi(argv[7]);
82			if (m <= 0 \|\| n <= 0)
83			goto userror;
84
85	greg	1.3	row0 = (POINT )malloc((n+1)sizeof(POINT));
86			row1 = (POINT )malloc((n+1)sizeof(POINT));
87			row2 = (POINT )malloc((n+1)sizeof(POINT));
88			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.3	/* print header */
93	greg	1.1	printhead(argc, argv);
94	greg	1.3	/* compute first two rows */
95			comprow(0.0, row1, n);
96			comprow(1.0/m, row2, n);
97			compnorms(row1, row1, row2, n);
98			/* for each row */
99	greg	1.1	for (i = 0; i < m; i++) {
100			/* compute next row */
101	greg	1.3	rp = row0;
102	greg	1.1	row0 = row1;
103	greg	1.3	row1 = row2;
104			row2 = rp;
105			if (i+2 <= m) {
106			comprow((double)(i+2)/m, row2, n);
107			compnorms(row0, row1, row2, n);
108			} else
109			compnorms(row0, row1, row1, 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			double st[2], step;
224
225			st[0] = s;
226			st[1] = 0.0;
227			step = 1.0 / siz;
228			while (siz-- >= 0) {
229	greg	1.3	row->p[0] = funvalue(XNAME, 2, st);
230			row->p[1] = funvalue(YNAME, 2, st);
231			row->p[2] = funvalue(ZNAME, 2, st);
232			row++;
233	greg	1.1	st[1] += step;
234			}
235	greg	1.3	}
236
237
238			compnorms(r0, r1, r2, siz) /* compute row of averaged normals */
239			register POINT r0, r1, *r2;
240			int siz;
241			{
242			FVECT v1, v2, vc;
243
244			if (!smooth) /* not needed if no smoothing */
245			return;
246			/* compute first point */
247			fvsum(v1, r2[0].p, r1[0].p, -1.0);
248			fvsum(v2, r1[1].p, r1[0].p, -1.0);
249			fcross(r1[0].n, v1, v2);
250			fvsum(v1, r0[0].p, r1[0].p, -1.0);
251			fcross(vc, v2, v1);
252			fvsum(r1[0].n, r1[0].n, vc, 1.0);
253			normalize(r1[0].n);
254			r0++; r1++; r2++;
255			/* compute middle points */
256			while (--siz > 0) {
257			fvsum(v1, r2[0].p, r1[0].p, -1.0);
258			fvsum(v2, r1[1].p, r1[0].p, -1.0);
259			fcross(r1[0].n, v1, v2);
260			fvsum(v1, r0[0].p, r1[0].p, -1.0);
261			fcross(vc, v2, v1);
262			fvsum(r1[0].n, r1[0].n, vc, 1.0);
263			fvsum(v2, r1[-1].p, r1[0].p, -1.0);
264			fcross(vc, v1, v2);
265			fvsum(r1[0].n, r1[0].n, vc, 1.0);
266			fvsum(v1, r2[0].p, r1[0].p, -1.0);
267			fcross(vc, v2, v1);
268			fvsum(r1[0].n, r1[0].n, vc, 1.0);
269			normalize(r1[0].n);
270			r0++; r1++; r2++;
271			}
272			/* compute end point */
273			fvsum(v1, r0[0].p, r1[0].p, -1.0);
274			fvsum(v2, r1[-1].p, r1[0].p, -1.0);
275			fcross(r1[0].n, v1, v2);
276			fvsum(v1, r2[0].p, r1[0].p, -1.0);
277			fcross(vc, v2, v1);
278			fvsum(r1[0].n, r1[0].n, vc, 1.0);
279			normalize(r1[0].n);
280			}
281
282
283			int
284			norminterp(resmat, p0, p1, p2, p3) /* compute normal interpolation */
285			register FVECT resmat[4];
286			POINT p0, p1, p2, p3;
287			{
288			#define u ((ax+1)%3)
289			#define v ((ax+2)%3)
290
291			register int ax;
292			double eqnmat[4][4], solmat[4][4];
293			FVECT v1;
294			register int i, j;
295
296			if (!smooth) /* no interpolation if no smoothing */
297			return(-1);
298			/* find dominant axis */
299			VCOPY(v1, p0->n);
300			fvsum(v1, v1, p1->n, 1.0);
301			fvsum(v1, v1, p2->n, 1.0);
302			fvsum(v1, v1, p3->n, 1.0);
303			ax = fabs(v1[0]) > fabs(v1[1]) ? 0 : 1;
304			ax = fabs(v1[ax]) > fabs(v1[2]) ? ax : 2;
305			/* assign equation matrix */
306			eqnmat[0][0] = p0->p[u]*p0->p[v];
307			eqnmat[0][1] = p0->p[u];
308			eqnmat[0][2] = p0->p[v];
309			eqnmat[0][3] = 1.0;
310			eqnmat[1][0] = p1->p[u]*p1->p[v];
311			eqnmat[1][1] = p1->p[u];
312			eqnmat[1][2] = p1->p[v];
313			eqnmat[1][3] = 1.0;
314			eqnmat[2][0] = p2->p[u]*p2->p[v];
315			eqnmat[2][1] = p2->p[u];
316			eqnmat[2][2] = p2->p[v];
317			eqnmat[2][3] = 1.0;
318			eqnmat[3][0] = p3->p[u]*p3->p[v];
319			eqnmat[3][1] = p3->p[u];
320			eqnmat[3][2] = p3->p[v];
321			eqnmat[3][3] = 1.0;
322			/* invert matrix (solve system) */
323			if (!invmat(solmat, eqnmat))
324			return(-1); /* no solution */
325			/* compute result matrix */
326			for (j = 0; j < 4; j++)
327			for (i = 0; i < 3; i++)
328			resmat[j][i] = solmat[j][0]*p0->n[i] +
329			solmat[j][1]*p1->n[i] +
330			solmat[j][2]*p2->n[i] +
331			solmat[j][3]*p3->n[i];
332			return(ax);
333
334			#undef u
335			#undef v
336			}
337
338
339			static double m4tmp[4][4]; /* for efficiency */
340
341			#define copymat4(m4a,m4b) bcopy((char )m4b,(char )m4a,sizeof(m4tmp))
342
343
344			setident4(m4)
345			double m4[4][4];
346			{
347			static double ident[4][4] = {
348			1.,0.,0.,0.,
349			0.,1.,0.,0.,
350			0.,0.,1.,0.,
351			0.,0.,0.,1.,
352			};
353			copymat4(m4, ident);
354			}
355
356			/*
357			* invmat - computes the inverse of mat into inverse. Returns 1
358			* if there exists an inverse, 0 otherwise. It uses Gaussian Elimination
359			* method.
360			*/
361
362			invmat(inverse,mat)
363			double mat[4][4],inverse[4][4];
364			{
365			#define SWAP(a,b,t) (t=a,a=b,b=t)
366
367			register int i,j,k;
368			register double temp;
369
370			setident4(inverse);
371			copymat4(m4tmp, mat);
372
373			for(i = 0; i < 4; i++) {
374			if(m4tmp[i][i] == 0) { /* Pivot is zero */
375			/* Look for a raw with pivot != 0 and swap raws */
376			for(j = i + 1; j < 4; j++)
377			if(m4tmp[j][i] != 0) {
378			for( k = 0; k < 4; k++) {
379			SWAP(m4tmp[i][k],m4tmp[j][k],temp);
380			SWAP(inverse[i][k],inverse[j][k],temp);
381			}
382			break;
383			}
384			if(j == 4) /* No replacing raw -> no inverse */
385			return(0);
386			}
387
388			temp = m4tmp[i][i];
389			for(k = 0; k < 4; k++) {
390			m4tmp[i][k] /= temp;
391			inverse[i][k] /= temp;
392			}
393			for(j = 0; j < 4; j++) {
394			if(j != i) {
395			temp = m4tmp[j][i];
396			for(k = 0; k < 4; k++) {
397			m4tmp[j][k] -= m4tmp[i][k]*temp;
398			inverse[j][k] -= inverse[i][k]*temp;
399			}
400			}
401			}
402			}
403			return(1);
404			#undef SWAP
405	greg	1.1	}
406
407
408			eputs(msg)
409			char *msg;
410			{
411			fputs(msg, stderr);
412			}
413
414
415			wputs(msg)
416			char *msg;
417			{
418			eputs(msg);
419			}
420
421
422			quit(code)
423			{
424			exit(code);
425			}
426
427
428			printhead(ac, av) /* print command header */
429			register int ac;
430			register char **av;
431			{
432			putchar('#');
433			while (ac--) {
434			putchar(' ');
435			fputs(*av++, stdout);
436			}
437			putchar('\n');
438			}
439
440
441			double
442			l_hermite()
443			{
444			double t;
445
446			t = argument(5);
447			return( argument(1)((2.0t-3.0)tt+1.0) +
448			argument(2)(-2.0t+3.0)tt +
449			argument(3)((t-2.0)t+1.0)*t +
450			argument(4)(t-1.0)t*t );
451			}