--- ray/src/gen/gensurf.c	1989/02/02 11:16:30	1.1
+++ ray/src/gen/gensurf.c	1991/11/12 17:04:49	2.1
@@ -1,8 +1,10 @@
-/*
-
 #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)
@@ -14,140 +16,106 @@ static char SCCSid[] = "$SunId$ LBL";
  *	4/3/87
  */
 
-#include  <stdio.h>
+#include  "standard.h"
 
-#define  XNAME		"X_"			/* x function name */
-#define  YNAME		"Y_"			/* y function name */
-#define  ZNAME		"Z_"			/* z function name */
+#define  XNAME		"X`SYS`"		/* x function name */
+#define  YNAME		"Y`SYS`"		/* y function name */
+#define  ZNAME		"Z`SYS`"		/* z function name */
 
-#define  PI		3.14159265358979323846
+#define  ABS(x)		((x)>=0 ? (x) : -(x))
 
-#define  FTINY		1e-7
+#define  pvect(p)	printf(vformat, (p)[0], (p)[1], (p)[2])
 
-#define  vertex(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";
 
-double  funvalue(), dist2(), fdot(), l_hermite(), argument();
+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[];
 {
-	static double  *xyz[4];
-	double  *row0, *row1, *dp;
-	double  v1[3], v2[3], vc1[3], vc2[3];
-	double  a1, a2;
+	extern long	eclock;
+	POINT  *row0, *row1, *row2, *rp;
 	int  i, j, m, n;
 	char  stmp[256];
-	double  d;
-	register int  k;
 
-	varset("PI", PI);
-	funset("hermite", 5, l_hermite);
+	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]);
+			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(NULL, stmp);
+	scompile(stmp, NULL, 0);
 	sprintf(stmp, "%s(s,t)=%s;", YNAME, argv[4]);
-	scompile(NULL, stmp);
+	scompile(stmp, NULL, 0);
 	sprintf(stmp, "%s(s,t)=%s;", ZNAME, argv[5]);
-	scompile(NULL, stmp);
+	scompile(stmp, NULL, 0);
 	m = atoi(argv[6]);
 	n = atoi(argv[7]);
 	if (m <= 0 || n <= 0)
 		goto userror;
 
-	row0 = (double *)malloc((n+1)*3*sizeof(double));
-	row1 = (double *)malloc((n+1)*3*sizeof(double));
-	if (row0 == NULL || row1 == NULL) {
+	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);
-
-	comprow(0.0, row1, n);			/* compute zeroeth row */
-
+	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 */
-		dp = row0;
+		rp = row0;
 		row0 = row1;
-		row1 = dp;
-		comprow((double)(i+1)/m, row1, n);
+		row1 = row2;
+		row2 = rp;
+		comprow((double)(i+2)/m, row2, n);
+		compnorms(row0, row1, row2, n);
 
 		for (j = 0; j < n; j++) {
-							/* get vertices */
-			xyz[0] = row0 + 3*j;
-			xyz[1] = row1 + 3*j;
-			xyz[2] = xyz[0] + 3;
-			xyz[3] = xyz[1] + 3;
-							/* rotate vertices */
-			if (dist2(xyz[0],xyz[3]) < dist2(xyz[1],xyz[2])-FTINY) {
-				dp = xyz[0];
-				xyz[0] = xyz[1];
-				xyz[1] = xyz[3];
-				xyz[3] = xyz[2];
-				xyz[2] = dp;
-			}
-							/* get normals */
-			for (k = 0; k < 3; k++) {
-				v1[k] = xyz[1][k] - xyz[0][k];
-				v2[k] = xyz[2][k] - xyz[0][k];
-			}
-			fcross(vc1, v1, v2);
-			a1 = fdot(vc1, vc1);
-			for (k = 0; k < 3; k++) {
-				v1[k] = xyz[2][k] - xyz[3][k];
-				v2[k] = xyz[1][k] - xyz[3][k];
-			}
-			fcross(vc2, v1, v2);
-			a2 = fdot(vc2, vc2);
-							/* check coplanar */
-			if (a1 > FTINY*FTINY && a2 > FTINY*FTINY) {
-				d = fdot(vc1, vc2);
-				if (d*d/a1/a2 >= 1.0-FTINY*FTINY) {
-					if (d > 0.0) {	/* coplanar */
-						printf(
-						"\n%s polygon %s.%d.%d\n",
-						argv[1], argv[2], i+1, j+1);
-						printf("0\n0\n12\n");
-						vertex(xyz[0]);
-						vertex(xyz[1]);
-						vertex(xyz[3]);
-						vertex(xyz[2]);
-					}		/* else overlapped */
-					continue;
-				}			/* else bent */
-			}
-							/* check triangles */
-			if (a1 > FTINY*FTINY) {
-				printf("\n%s polygon %s.%da%d\n",
-					argv[1], argv[2], i+1, j+1);
-				printf("0\n0\n9\n");
-				vertex(xyz[0]);
-				vertex(xyz[1]);
-				vertex(xyz[2]);
-			}
-			if (a2 > FTINY*FTINY) {
-				printf("\n%s polygon %s.%db%d\n",
-					argv[1], argv[2], i+1, j+1);
-				printf("0\n0\n9\n");
-				vertex(xyz[2]);
-				vertex(xyz[1]);
-				vertex(xyz[3]);
-			}
+							/* 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]);
 		}
 	}
 
@@ -155,30 +123,261 @@ char  *argv[];
 
 userror:
 	fprintf(stderr, "Usage: %s material name ", argv[0]);
-	fprintf(stderr, "x(s,t) y(s,t) z(s,t) m n [-e expr] [-f file]\n");
+	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 double  *row;
+register POINT  *row;
 int  siz;
 {
-	double  st[2], step;
-
+	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;
-	st[1] = 0.0;
-	step = 1.0 / siz;
+	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) {
-		*row++ = funvalue(XNAME, 2, st);
-		*row++ = funvalue(YNAME, 2, st);
-		*row++ = funvalue(ZNAME, 2, st);
-		st[1] += step;
+		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;
 {
@@ -222,4 +421,30 @@ l_hermite()			
 		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) );
 }