--- ray/src/gen/gensurf.c 1989/10/18 15:01:23 1.3 +++ ray/src/gen/gensurf.c 1989/10/18 18:49:09 1.4 @@ -27,6 +27,8 @@ static char SCCSid[] = "$SunId$ LBL"; #define FTINY 1e-7 +#define ABS(x) ((x)>=0 ? (x) : -(x)) + #define pvect(p) printf(vformat, (p)[0], (p)[1], (p)[2]) char vformat[] = "%15.9g %15.9g %15.9g\n"; @@ -37,7 +39,7 @@ int smooth = 0; /* apply smoothing? */ char *modname, *surfname; -double funvalue(), l_hermite(), argument(), fabs(); +double funvalue(), l_hermite(), argument(); typedef struct { FVECT p; /* vertex position */ @@ -82,19 +84,21 @@ 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 */ + /* 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 +106,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 +221,15 @@ double s; register POINT *row; int siz; { - double st[2], step; - + double st[2]; + register int i; + /* compute one past each end */ 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; + for (i = -1; i <= siz+1; i++) { + 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); } } @@ -240,20 +239,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 +260,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 +272,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 +283,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 +303,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 +319,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 +330,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(mat, 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 +372,7 @@ double mat[4][4],inverse[4][4]; } } return(1); + #undef SWAP }