--- ray/src/gen/gensurf.c	1990/01/19 00:03:08	1.5
+++ ray/src/gen/gensurf.c	2004/08/21 11:54:06	2.16
@@ -1,9 +1,6 @@
-/* Copyright (c) 1989 Regents of the University of California */
-
 #ifndef lint
-static char SCCSid[] = "$SunId$ LBL";
+static const char RCSid[] = "$Id: gensurf.c,v 2.16 2004/08/21 11:54:06 greg Exp $";
 #endif
-
 /*
  *  gensurf.c - program to generate functional surfaces
  *
@@ -14,16 +11,27 @@ static char SCCSid[] = "$SunId$ LBL";
  *  rule applied to (s,t).
  *
  *	4/3/87
+ *
+ *	4/16/02	Added conditional vertex output
  */
 
 #include  "standard.h"
 
-#define  XNAME		"X_"			/* x function name */
-#define  YNAME		"Y_"			/* y function name */
-#define  ZNAME		"Z_"			/* z function name */
+#include  "paths.h"
+#include  "resolu.h"
+#include  "rterror.h"
+#include  "calcomp.h"
 
+char  XNAME[] =		"X`SYS";		/* x function name */
+char  YNAME[] =		"Y`SYS";		/* y function name */
+char  ZNAME[] =		"Z`SYS";		/* z function name */
+
+char  VNAME[] = 	"valid";		/* valid vertex name */
+
 #define  ABS(x)		((x)>=0 ? (x) : -(x))
 
+#define  ZEROVECT(v)	(DOT(v,v) <= FTINY*FTINY)
+
 #define  pvect(p)	printf(vformat, (p)[0], (p)[1], (p)[2])
 
 char  vformat[] = "%15.9g %15.9g %15.9g\n";
@@ -31,17 +39,41 @@ char  tsargs[] = "4 surf_dx surf_dy surf_dz surf.cal\n
 char  texname[] = "Phong";
 
 int  smooth = 0;		/* apply smoothing? */
+int  objout = 0;		/* output .OBJ format? */
 
 char  *modname, *surfname;
 
-double  funvalue(), l_hermite(), argument();
+				/* recorded data flags */
+#define  HASBORDER	01
+#define  TRIPLETS	02
+				/* a data structure */
+struct {
+	int	flags;			/* data type */
+	short	m, n;			/* number of s and t values */
+	RREAL	*data;			/* the data itself, s major sort */
+} datarec;			/* our recorded data */
 
+/* XXX this is redundant with rt/noise3.c, should go to a library */
+double  l_hermite(), l_bezier(), l_bspline(), l_dataval();
+
 typedef struct {
+	int  valid;	/* point is valid (vertex number) */
 	FVECT  p;	/* vertex position */
 	FVECT  n;	/* average normal */
+	RREAL  uv[2];	/* (u,v) position */
 } POINT;
 
 
+void loaddata(char *file, int m, int n, int pointsize);
+double l_dataval(char *nam);
+void putobjrow(POINT *rp, int n);
+void putsquare(POINT *p0, POINT *p1, POINT *p2, POINT *p3);
+void comprow(double s, POINT *row, int siz);
+void compnorms(POINT *r0, POINT *r1, POINT *r2, int siz);
+int norminterp(FVECT resmat[4], POINT *p0, POINT *p1, POINT *p2, POINT *p3);
+
+
+int
 main(argc, argv)
 int  argc;
 char  *argv[];
@@ -50,35 +82,53 @@ char  *argv[];
 	int  i, j, m, n;
 	char  stmp[256];
 
-	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 if (!strcmp(argv[i], "-o"))
+			objout++;
 		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;
-
+	if (!strcmp(argv[5], "-") || access(argv[5], 4) == 0) {	/* file? */
+		funset(ZNAME, 2, ':', l_dataval);
+		if (!strcmp(argv[5],argv[3]) && !strcmp(argv[5],argv[4])) {
+			loaddata(argv[5], m, n, 3);
+			funset(XNAME, 2, ':', l_dataval);
+			funset(YNAME, 2, ':', l_dataval);
+		} else {
+			loaddata(argv[5], m, n, 1);
+			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);
+		}
+	} else {
+		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);
+	}
 	row0 = (POINT *)malloc((n+3)*sizeof(POINT));
 	row1 = (POINT *)malloc((n+3)*sizeof(POINT));
 	row2 = (POINT *)malloc((n+3)*sizeof(POINT));
@@ -88,12 +138,18 @@ char  *argv[];
 	}
 	row0++; row1++; row2++;
 						/* print header */
-	printhead(argc, argv);
+	fputs("# ", stdout);
+	printargs(argc, argv, stdout);
+	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);
+	if (objout) {
+		printf("\nusemtl %s\n\n", modname);
+		putobjrow(row1, n);
+	}
 						/* for each row */
 	for (i = 0; i < m; i++) {
 						/* compute next row */
@@ -103,10 +159,17 @@ char  *argv[];
 		row2 = rp;
 		comprow((double)(i+2)/m, row2, n);
 		compnorms(row0, row1, row2, n);
+		if (objout)
+			putobjrow(row1, n);
 
 		for (j = 0; j < n; j++) {
+			int  orient = (j & 1);
 							/* put polygons */
-			if ((i+j) & 1)
+			if (!(row0[j].valid && row1[j+1].valid))
+				orient = 1;
+			else if (!(row1[j].valid && row0[j+1].valid))
+				orient = 0;
+			if (orient)
 				putsquare(&row0[j], &row1[j],
 						&row0[j+1], &row1[j+1]);
 			else
@@ -115,34 +178,203 @@ char  *argv[];
 		}
 	}
 
-	quit(0);
+	return 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);
+	return 1;
 }
 
 
-putsquare(p0, p1, p2, p3)		/* put out a square */
-POINT  *p0, *p1, *p2, *p3;
+void
+loaddata(		/* load point data from file */
+	char  *file,
+	int  m,
+	int  n,
+	int  pointsize
+)
 {
+	FILE  *fp;
+	char  word[64];
+	register int  size;
+	register RREAL  *dp;
+
+	datarec.flags = HASBORDER;		/* assume border values */
+	datarec.m = m+1;
+	datarec.n = n+1;
+	size = datarec.m*datarec.n*pointsize;
+	if (pointsize == 3)
+		datarec.flags |= TRIPLETS;
+	dp = (RREAL *)malloc(size*sizeof(RREAL));
+	if ((datarec.data = dp) == NULL) {
+		fputs("Out of memory\n", stderr);
+		exit(1);
+	}
+	if (!strcmp(file, "-")) {
+		file = "<stdin>";
+		fp = stdin;
+	} else if ((fp = fopen(file, "r")) == NULL) {
+		fputs(file, stderr);
+		fputs(": cannot open\n", stderr);
+		exit(1);
+	}
+	while (size > 0 && fgetword(word, sizeof(word), fp) != NULL) {
+		if (!isflt(word)) {
+			fprintf(stderr, "%s: garbled data value: %s\n",
+					file, word);
+			exit(1);
+		}
+		*dp++ = atof(word);
+		size--;
+	}
+	if (size == (m+n+1)*pointsize) {	/* no border after all */
+		dp = (RREAL *)realloc((void *)datarec.data,
+				m*n*pointsize*sizeof(RREAL));
+		if (dp != NULL)
+			datarec.data = dp;
+		datarec.flags &= ~HASBORDER;
+		datarec.m = m;
+		datarec.n = n;
+		size = 0;
+	}
+	if (datarec.m < 2 || datarec.n < 2 || size != 0 ||
+			fgetword(word, sizeof(word), fp) != NULL) {
+		fputs(file, stderr);
+		fputs(": bad number of data points\n", stderr);
+		exit(1);
+	}
+	fclose(fp);
+}
+
+
+double
+l_dataval(				/* return recorded data value */
+	char  *nam
+)
+{
+	double  u, v;
+	register int  i, j;
+	register RREAL  *dp;
+	double  d00, d01, d10, d11;
+						/* compute coordinates */
+	u = argument(1); v = argument(2);
+	if (datarec.flags & HASBORDER) {
+		i = u *= datarec.m-1;
+		j = v *= datarec.n-1;
+	} else {
+		i = u = u*datarec.m - .5;
+		j = v = v*datarec.n - .5;
+	}
+	if (i < 0) i = 0;
+	else if (i > datarec.m-2) i = datarec.m-2;
+	if (j < 0) j = 0;
+	else if (j > datarec.n-2) j = datarec.n-2;
+						/* compute value */
+	if (datarec.flags & TRIPLETS) {
+		dp = datarec.data + 3*(j*datarec.m + i);
+		if (nam == ZNAME)
+			dp += 2;
+		else if (nam == YNAME)
+			dp++;
+		d00 = dp[0]; d01 = dp[3];
+		dp += 3*datarec.m;
+		d10 = dp[0]; d11 = dp[3];
+	} else {
+		dp = datarec.data + j*datarec.m + i;
+		d00 = dp[0]; d01 = dp[1];
+		dp += datarec.m;
+		d10 = dp[0]; d11 = dp[1];
+	}
+						/* bilinear interpolation */
+	return((j+1-v)*((i+1-u)*d00+(u-i)*d01)+(v-j)*((i+1-u)*d10+(u-i)*d11));
+}
+
+
+void
+putobjrow(			/* output vertex row to .OBJ */
+	register POINT  *rp,
+	int  n
+)
+{
+	static int	nverts = 0;
+
+	for ( ; n-- >= 0; rp++) {
+		if (!rp->valid)
+			continue;
+		fputs("v ", stdout);
+		pvect(rp->p);
+		if (smooth && !ZEROVECT(rp->n))
+			printf("\tvn %.9g %.9g %.9g\n",
+					rp->n[0], rp->n[1], rp->n[2]);
+		printf("\tvt %.9g %.9g\n", rp->uv[0], rp->uv[1]);
+		rp->valid = ++nverts;
+	}
+}
+
+
+void
+putsquare(		/* put out a square */
+	POINT *p0,
+	POINT *p1,
+	POINT *p2,
+	POINT *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;
+	ok1 = (p0->valid && p1->valid && p2->valid);
+	if (ok1) {
+		VSUB(v1, p1->p, p0->p);
+		VSUB(v2, p2->p, p0->p);
+		fcross(vc1, v1, v2);
+		ok1 = (normalize(vc1) != 0.0);
+	}
+	ok2 = (p1->valid && p2->valid && p3->valid);
+	if (ok2) {
+		VSUB(v1, p2->p, p3->p);
+		VSUB(v2, p1->p, p3->p);
+		fcross(vc2, v1, v2);
+		ok2 = (normalize(vc2) != 0.0);
+	}
 	if (!(ok1 | ok2))
 		return;
+	if (objout) {			/* output .OBJ faces */
+		int	p0n=0, p1n=0, p2n=0, p3n=0;
+		if (smooth) {
+			if (!ZEROVECT(p0->n))
+				p0n = p0->valid;
+			if (!ZEROVECT(p1->n))
+				p1n = p1->valid;
+			if (!ZEROVECT(p2->n))
+				p2n = p2->valid;
+			if (!ZEROVECT(p3->n))
+				p3n = p3->valid;
+		}
+		if (ok1 & ok2 && fdot(vc1,vc2) >= 1.0-FTINY*FTINY) {
+			printf("f %d/%d/%d %d/%d/%d %d/%d/%d %d/%d/%d\n",
+					p0->valid, p0->valid, p0n,
+					p1->valid, p1->valid, p1n,
+					p3->valid, p3->valid, p3n,
+					p2->valid, p2->valid, p2n);
+			return;
+		}
+		if (ok1)
+			printf("f %d/%d/%d %d/%d/%d %d/%d/%d\n",
+					p0->valid, p0->valid, p0n,
+					p1->valid, p1->valid, p1n,
+					p2->valid, p2->valid, p2n);
+		if (ok2)
+			printf("f %d/%d/%d %d/%d/%d %d/%d/%d\n",
+					p2->valid, p2->valid, p2n,
+					p1->valid, p1->valid, p1n,
+					p3->valid, p3->valid, p3n);
+		return;
+	}
 					/* compute normal interpolation */
 	axis = norminterp(norm, p0, p1, p2, p3);
 
@@ -211,47 +443,85 @@ POINT  *p0, *p1, *p2, *p3;
 }
 
 
-comprow(s, row, siz)			/* compute row of values */
-double  s;
-register POINT  *row;
-int  siz;
+void
+comprow(			/* compute row of values */
+	double  s,
+	register POINT  *row,
+	int  siz
+)
 {
 	double  st[2];
+	int  end;
+	int  checkvalid;
 	register int  i;
-					/* compute one past each end */
+	
+	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;
-	for (i = -1; i <= siz+1; i++) {
+	checkvalid = (fundefined(VNAME) == 2);
+	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);
+		if (checkvalid && funvalue(VNAME, 2, st) <= 0.0) {
+			row[i].valid = 0;
+			row[i].p[0] = row[i].p[1] = row[i].p[2] = 0.0;
+			row[i].uv[0] = row[i].uv[1] = 0.0;
+		} else {
+			row[i].valid = 1;
+			row[i].p[0] = funvalue(XNAME, 2, st);
+			row[i].p[1] = funvalue(YNAME, 2, st);
+			row[i].p[2] = funvalue(ZNAME, 2, st);
+			row[i].uv[0] = st[0];
+			row[i].uv[1] = st[1];
+		}
+		i++;
 	}
 }
 
 
-compnorms(r0, r1, r2, siz)		/* compute row of averaged normals */
-register POINT  *r0, *r1, *r2;
-int  siz;
+void
+compnorms(		/* compute row of averaged normals */
+	register POINT  *r0,
+	register POINT  *r1,
+	register POINT  *r2,
+	int  siz
+)
 {
-	FVECT  v1, v2, vc;
-	register int  i;
+	FVECT  v1, v2;
 
 	if (!smooth)			/* not needed if no smoothing */
 		return;
-					/* compute middle points */
+					/* compute row 1 normals */
 	while (siz-- >= 0) {
-		fvsum(v1, r2[0].p, r1[0].p, -1.0);
-		fvsum(v2, r1[1].p, r1[0].p, -1.0);
+		if (!r1[0].valid)
+			continue;
+		if (!r0[0].valid) {
+			if (!r2[0].valid) {
+				r1[0].n[0] = r1[0].n[1] = r1[0].n[2] = 0.0;
+				continue;
+			}
+			fvsum(v1, r2[0].p, r1[0].p, -1.0);
+		} else if (!r2[0].valid)
+			fvsum(v1, r1[0].p, r0[0].p, -1.0);
+		else
+			fvsum(v1, r2[0].p, r0[0].p, -1.0);
+		if (!r1[-1].valid) {
+			if (!r1[1].valid) {
+				r1[0].n[0] = r1[0].n[1] = r1[0].n[2] = 0.0;
+				continue;
+			}
+			fvsum(v2, r1[1].p, r1[0].p, -1.0);
+		} else if (!r1[1].valid)
+			fvsum(v2, r1[0].p, r1[-1].p, -1.0);
+		else
+			fvsum(v2, r1[1].p, r1[-1].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++;
 	}
@@ -259,15 +529,19 @@ int  siz;
 
 
 int
-norminterp(resmat, p0, p1, p2, p3)	/* compute normal interpolation */
-register FVECT  resmat[4];
-POINT  *p0, *p1, *p2, *p3;
+norminterp(	/* compute normal interpolation */
+	register FVECT  resmat[4],
+	POINT  *p0,
+	POINT  *p1,
+	POINT  *p2,
+	POINT  *p3
+)
 {
 #define u  ((ax+1)%3)
 #define v  ((ax+2)%3)
 
 	register int  ax;
-	double  eqnmat[4][4];
+	MAT4  eqnmat;
 	FVECT  v1;
 	register int  i, j;
 
@@ -298,7 +572,7 @@ POINT  *p0, *p1, *p2, *p3;
 	eqnmat[3][2] = p3->p[v];
 	eqnmat[3][3] = 1.0;
 					/* invert matrix (solve system) */
-	if (!invmat(eqnmat, eqnmat))
+	if (!invmat4(eqnmat, eqnmat))
 		return(-1);			/* no solution */
 					/* compute result matrix */
 	for (j = 0; j < 4; j++)
@@ -314,105 +588,40 @@ POINT  *p0, *p1, *p2, *p3;
 }
 
 
-/*
- * 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];
+double
+l_hermite(char *nm)
 {
-#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
+	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 );
 }
 
 
-eputs(msg)
-char  *msg;
+double
+l_bezier(char *nm)
 {
-	fputs(msg, stderr);
-}
+	double  t;
 
-
-wputs(msg)
-char  *msg;
-{
-	eputs(msg);
+	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 );
 }
 
 
-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()			
+l_bspline(char *nm)
 {
 	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 );
+	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) );
 }