--- ray/src/gen/gensurf.c 1989/10/18 15:01:23 1.3 +++ ray/src/gen/gensurf.c 1990/07/09 09:55:47 1.9 @@ -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 -#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,12 +46,15 @@ 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); funset("hermite", 5, l_hermite); + funset("bezier", 5, l_bezier); + funset("bspline", 5, l_bspline); if (argc < 8) goto userror; @@ -82,19 +82,22 @@ char *argv[]; 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++; } } @@ -240,20 +249,12 @@ register POINT *r0, *r1, *r2; int siz; { FVECT v1, v2, vc; + 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) { + 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); @@ -269,14 +270,6 @@ int siz; 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 +282,7 @@ POINT *p0, *p1, *p2, *p3; #define v ((ax+2)%3) register int ax; - double eqnmat[4][4], solmat[4][4]; + double eqnmat[4][4]; FVECT v1; register int i, j; @@ -300,8 +293,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 +313,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 +329,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 @@ -364,26 +340,31 @@ double mat[4][4],inverse[4][4]; { #define SWAP(a,b,t) (t=a,a=b,b=t) + double m4tmp[4][4]; register int i,j,k; register double temp; - setident4(inverse); - copymat4(m4tmp, mat); + 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++) { - 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 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++) { @@ -401,6 +382,7 @@ double mat[4][4],inverse[4][4]; } } return(1); + #undef SWAP } @@ -448,4 +430,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) ); }