--- ray/src/gen/gensurf.c	1989/10/18 15:01:23	1.3
+++ ray/src/gen/gensurf.c	1991/04/23 13:04:58	1.13
@@ -1,9 +1,9 @@
-/* Copyright (c) 1989 Regents of the University of California */
-
 #ifndef lint
 static char SCCSid[] = "$SunId$ LBL";
 #endif
 
+/* Copyright (c) 1989 Regents of the University of California */
+
 /*
  *  gensurf.c - program to generate functional surfaces
  *
@@ -16,17 +16,14 @@ static char SCCSid[] = "$SunId$ LBL";
  *	4/3/87
  */
 
-#include  <stdio.h>
-#include  "fvect.h"
+#include  "standard.h"
 
 #define  XNAME		"X_"			/* x function name */
 #define  YNAME		"Y_"			/* y function name */
 #define  ZNAME		"Z_"			/* 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])
 
 char  vformat[] = "%15.9g %15.9g %15.9g\n";
@@ -37,7 +34,7 @@ int  smooth = 0;		/* apply smoothing? */
 
 char  *modname, *surfname;
 
-double  funvalue(), l_hermite(), argument(), fabs();
+double  funvalue(), l_hermite(), l_bezier(), l_bspline(), argument();
 
 typedef struct {
 	FVECT  p;	/* vertex position */
@@ -49,19 +46,22 @@ 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);
+	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"))
@@ -72,29 +72,32 @@ char  *argv[];
 	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 = (POINT *)malloc((n+1)*sizeof(POINT));
-	row1 = (POINT *)malloc((n+1)*sizeof(POINT));
-	row2 = (POINT *)malloc((n+1)*sizeof(POINT));
+	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);
-						/* compute first two rows */
+	eclock = 0;
+						/* initialize */
+	comprow(-1.0/m, row0, n);
 	comprow(0.0, row1, n);
 	comprow(1.0/m, row2, n);
-	compnorms(row1, row1, row2, n);
+	compnorms(row0, row1, row2, n);
 						/* for each row */
 	for (i = 0; i < m; i++) {
 						/* compute next row */
@@ -102,11 +105,8 @@ char  *argv[];
 		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);
+		comprow((double)(i+2)/m, row2, n);
+		compnorms(row0, row1, row2, n);
 
 		for (j = 0; j < n; j++) {
 							/* put polygons */
@@ -220,17 +220,26 @@ double  s;
 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 (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;
+	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++;
 	}
 }
 
@@ -239,44 +248,19 @@ compnorms(r0, r1, r2, siz)		/* compute row of averaged
 register POINT  *r0, *r1, *r2;
 int  siz;
 {
-	FVECT  v1, v2, vc;
+	FVECT  v1, v2;
+	register int  i;
 
 	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);
+	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);
-		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);
 }
 
 
@@ -289,7 +273,7 @@ POINT  *p0, *p1, *p2, *p3;
 #define v  ((ax+2)%3)
 
 	register int  ax;
-	double  eqnmat[4][4], solmat[4][4];
+	MAT4  eqnmat;
 	FVECT  v1;
 	register int  i, j;
 
@@ -300,8 +284,8 @@ POINT  *p0, *p1, *p2, *p3;
 	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;
+	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];
@@ -320,15 +304,15 @@ POINT  *p0, *p1, *p2, *p3;
 	eqnmat[3][2] = p3->p[v];
 	eqnmat[3][3] = 1.0;
 					/* invert matrix (solve system) */
-	if (!invmat(solmat, eqnmat))
+	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] =	solmat[j][0]*p0->n[i] +
-					solmat[j][1]*p1->n[i] +
-					solmat[j][2]*p2->n[i] +
-					solmat[j][3]*p3->n[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
@@ -336,23 +320,6 @@ POINT  *p0, *p1, *p2, *p3;
 }
 
 
-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
@@ -360,30 +327,35 @@ double  m4[4][4];
  */
 
 invmat(inverse,mat)
-double mat[4][4],inverse[4][4];
+MAT4  inverse, mat;
 {
 #define SWAP(a,b,t) (t=a,a=b,b=t)
 
+	MAT4  m4tmp;
 	register int i,j,k;
 	register double temp;
 
-	setident4(inverse);
 	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++) {
-		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);
-			}
+		/* 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++) {
@@ -401,6 +373,7 @@ double mat[4][4],inverse[4][4];
 			}
 		}
 	return(1);
+
 #undef SWAP
 }
 
@@ -448,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) );
 }