--- ray/src/gen/gensurf.c 1992/01/30 14:11:39 2.2 +++ ray/src/gen/gensurf.c 2003/02/22 02:07:23 2.6 @@ -1,9 +1,6 @@ #ifndef lint -static char SCCSid[] = "$SunId$ LBL"; +static const char RCSid[] = "$Id: gensurf.c,v 2.6 2003/02/22 02:07:23 greg Exp $"; #endif - -/* Copyright (c) 1989 Regents of the University of California */ - /* * gensurf.c - program to generate functional surfaces * @@ -14,14 +11,18 @@ static char SCCSid[] = "$SunId$ LBL"; * rule applied to (s,t). * * 4/3/87 + * + * 4/16/02 Added conditional vertex output */ #include "standard.h" -char XNAME[] = "X`SYS`"; /* x function name */ -char YNAME[] = "Y`SYS`"; /* y function name */ -char ZNAME[] = "Z`SYS`"; /* z function name */ +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 pvect(p) printf(vformat, (p)[0], (p)[1], (p)[2]) @@ -48,6 +49,7 @@ double l_hermite(), l_bezier(), l_bspline(), l_datava extern double funvalue(), argument(); typedef struct { + int valid; /* point is valid */ FVECT p; /* vertex position */ FVECT n; /* average normal */ } POINT; @@ -134,8 +136,13 @@ char *argv[]; compnorms(row0, row1, row2, 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 @@ -158,14 +165,15 @@ char *file; int m, n; int pointsize; { - extern char *fgetword(); FILE *fp; char word[64]; register int size; register FLOAT *dp; datarec.flags = HASBORDER; /* assume border values */ - size = (m+1)*(n+1)*pointsize; + datarec.m = m+1; + datarec.n = n+1; + size = datarec.m*datarec.n*pointsize; if (pointsize == 3) datarec.flags |= TRIPLETS; dp = (FLOAT *)malloc(size*sizeof(FLOAT)); @@ -196,9 +204,12 @@ int pointsize; if (dp != NULL) datarec.data = dp; datarec.flags &= ~HASBORDER; + datarec.m = m; + datarec.n = n; size = 0; } - if (size || fgetword(word, sizeof(word), fp) != NULL) { + 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); @@ -218,11 +229,11 @@ char *nam; /* compute coordinates */ u = argument(1); v = argument(2); if (datarec.flags & HASBORDER) { - i = u *= datarec.m; - j = v *= datarec.n; + i = u *= datarec.m-1; + j = v *= datarec.n-1; } else { - i = u = u*(datarec.m+1) - .5; - j = v = v*(datarec.n+1) - .5; + 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; @@ -230,18 +241,18 @@ char *nam; else if (j > datarec.n-2) j = datarec.n-2; /* compute value */ if (datarec.flags & TRIPLETS) { - dp = datarec.data + 3*(j*datarec.n + i); - if (nam == YNAME) - dp++; - else if (nam == ZNAME) + 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.n; + dp += 3*datarec.m; d10 = dp[0]; d11 = dp[3]; } else { - dp = datarec.data + j*datarec.n + i; + dp = datarec.data + j*datarec.m + i; d00 = dp[0]; d01 = dp[1]; - dp += datarec.n; + dp += datarec.m; d10 = dp[0]; d11 = dp[1]; } /* bilinear interpolation */ @@ -258,14 +269,20 @@ POINT *p0, *p1, *p2, *p3; 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) { + 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); + } + ok2 = (p1->valid & p2->valid & p3->valid); + if (ok2) { + 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 */ @@ -343,6 +360,7 @@ int siz; { double st[2]; int end; + int checkvalid; register int i; if (smooth) { @@ -355,11 +373,18 @@ int siz; end = siz; } st[0] = s; + 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; + } 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); + } i++; } } @@ -370,14 +395,33 @@ register POINT *r0, *r1, *r2; int siz; { FVECT v1, v2; - register int i; if (!smooth) /* not needed if no smoothing */ return; - /* compute middle points */ + /* compute row 1 normals */ while (siz-- >= 0) { - fvsum(v1, r2[0].p, r0[0].p, -1.0); - fvsum(v2, r1[1].p, r1[-1].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); normalize(r1[0].n); r0++; r1++; r2++; @@ -425,7 +469,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++) @@ -441,64 +485,7 @@ 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) -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 -} - - +void eputs(msg) char *msg; { @@ -506,6 +493,7 @@ char *msg; } +void wputs(msg) char *msg; { @@ -513,7 +501,9 @@ char *msg; } +void quit(code) +int code; { exit(code); }