ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/gensurf.c
(Generate patch)

Comparing ray/src/gen/gensurf.c (file contents):
Revision 1.3 by greg, Wed Oct 18 15:01:23 1989 UTC vs.
Revision 1.4 by greg, Wed Oct 18 18:49:09 1989 UTC

# Line 27 | Line 27 | static char SCCSid[] = "$SunId$ LBL";
27  
28   #define  FTINY          1e-7
29  
30 + #define  ABS(x)         ((x)>=0 ? (x) : -(x))
31 +
32   #define  pvect(p)       printf(vformat, (p)[0], (p)[1], (p)[2])
33  
34   char  vformat[] = "%15.9g %15.9g %15.9g\n";
# Line 37 | Line 39 | int  smooth = 0;               /* apply smoothing? */
39  
40   char  *modname, *surfname;
41  
42 < double  funvalue(), l_hermite(), argument(), fabs();
42 > double  funvalue(), l_hermite(), argument();
43  
44   typedef struct {
45          FVECT  p;       /* vertex position */
# Line 82 | Line 84 | char  *argv[];
84          if (m <= 0 || n <= 0)
85                  goto userror;
86  
87 <        row0 = (POINT *)malloc((n+1)*sizeof(POINT));
88 <        row1 = (POINT *)malloc((n+1)*sizeof(POINT));
89 <        row2 = (POINT *)malloc((n+1)*sizeof(POINT));
87 >        row0 = (POINT *)malloc((n+3)*sizeof(POINT));
88 >        row1 = (POINT *)malloc((n+3)*sizeof(POINT));
89 >        row2 = (POINT *)malloc((n+3)*sizeof(POINT));
90          if (row0 == NULL || row1 == NULL || row2 == NULL) {
91                  fprintf(stderr, "%s: out of memory\n", argv[0]);
92                  quit(1);
93          }
94 +        row0++; row1++; row2++;
95                                                  /* print header */
96          printhead(argc, argv);
97 <                                                /* compute first two rows */
97 >                                                /* initialize */
98 >        comprow(-1.0/m, row0, n);
99          comprow(0.0, row1, n);
100          comprow(1.0/m, row2, n);
101 <        compnorms(row1, row1, row2, n);
101 >        compnorms(row0, row1, row2, n);
102                                                  /* for each row */
103          for (i = 0; i < m; i++) {
104                                                  /* compute next row */
# Line 102 | Line 106 | char  *argv[];
106                  row0 = row1;
107                  row1 = row2;
108                  row2 = rp;
109 <                if (i+2 <= m) {
110 <                        comprow((double)(i+2)/m, row2, n);
107 <                        compnorms(row0, row1, row2, n);
108 <                } else
109 <                        compnorms(row0, row1, row1, n);
109 >                comprow((double)(i+2)/m, row2, n);
110 >                compnorms(row0, row1, row2, n);
111  
112                  for (j = 0; j < n; j++) {
113                                                          /* put polygons */
# Line 220 | Line 221 | double  s;
221   register POINT  *row;
222   int  siz;
223   {
224 <        double  st[2], step;
225 <
224 >        double  st[2];
225 >        register int  i;
226 >                                        /* compute one past each end */
227          st[0] = s;
228 <        st[1] = 0.0;
229 <        step = 1.0 / siz;
230 <        while (siz-- >= 0) {
231 <                row->p[0] = funvalue(XNAME, 2, st);
232 <                row->p[1] = funvalue(YNAME, 2, st);
231 <                row->p[2] = funvalue(ZNAME, 2, st);
232 <                row++;
233 <                st[1] += step;
228 >        for (i = -1; i <= siz+1; i++) {
229 >                st[1] = (double)i/siz;
230 >                row[i].p[0] = funvalue(XNAME, 2, st);
231 >                row[i].p[1] = funvalue(YNAME, 2, st);
232 >                row[i].p[2] = funvalue(ZNAME, 2, st);
233          }
234   }
235  
# Line 240 | Line 239 | register POINT  *r0, *r1, *r2;
239   int  siz;
240   {
241          FVECT  v1, v2, vc;
242 +        register int  i;
243  
244          if (!smooth)                    /* not needed if no smoothing */
245                  return;
246                                        /* compute first point */
247        fvsum(v1, r2[0].p, r1[0].p, -1.0);
248        fvsum(v2, r1[1].p, r1[0].p, -1.0);
249        fcross(r1[0].n, v1, v2);
250        fvsum(v1, r0[0].p, r1[0].p, -1.0);
251        fcross(vc, v2, v1);
252        fvsum(r1[0].n, r1[0].n, vc, 1.0);
253        normalize(r1[0].n);
254        r0++; r1++; r2++;
246                                          /* compute middle points */
247 <        while (--siz > 0) {
247 >        while (siz-- >= 0) {
248                  fvsum(v1, r2[0].p, r1[0].p, -1.0);
249                  fvsum(v2, r1[1].p, r1[0].p, -1.0);
250                  fcross(r1[0].n, v1, v2);
# Line 269 | Line 260 | int  siz;
260                  normalize(r1[0].n);
261                  r0++; r1++; r2++;
262          }
272                                        /* compute end point */
273        fvsum(v1, r0[0].p, r1[0].p, -1.0);
274        fvsum(v2, r1[-1].p, r1[0].p, -1.0);
275        fcross(r1[0].n, v1, v2);
276        fvsum(v1, r2[0].p, r1[0].p, -1.0);
277        fcross(vc, v2, v1);
278        fvsum(r1[0].n, r1[0].n, vc, 1.0);
279        normalize(r1[0].n);
263   }
264  
265  
# Line 289 | Line 272 | POINT  *p0, *p1, *p2, *p3;
272   #define v  ((ax+2)%3)
273  
274          register int  ax;
275 <        double  eqnmat[4][4], solmat[4][4];
275 >        double  eqnmat[4][4];
276          FVECT  v1;
277          register int  i, j;
278  
# Line 300 | Line 283 | POINT  *p0, *p1, *p2, *p3;
283          fvsum(v1, v1, p1->n, 1.0);
284          fvsum(v1, v1, p2->n, 1.0);
285          fvsum(v1, v1, p3->n, 1.0);
286 <        ax = fabs(v1[0]) > fabs(v1[1]) ? 0 : 1;
287 <        ax = fabs(v1[ax]) > fabs(v1[2]) ? ax : 2;
286 >        ax = ABS(v1[0]) > ABS(v1[1]) ? 0 : 1;
287 >        ax = ABS(v1[ax]) > ABS(v1[2]) ? ax : 2;
288                                          /* assign equation matrix */
289          eqnmat[0][0] = p0->p[u]*p0->p[v];
290          eqnmat[0][1] = p0->p[u];
# Line 320 | Line 303 | POINT  *p0, *p1, *p2, *p3;
303          eqnmat[3][2] = p3->p[v];
304          eqnmat[3][3] = 1.0;
305                                          /* invert matrix (solve system) */
306 <        if (!invmat(solmat, eqnmat))
306 >        if (!invmat(eqnmat, eqnmat))
307                  return(-1);                     /* no solution */
308                                          /* compute result matrix */
309          for (j = 0; j < 4; j++)
310                  for (i = 0; i < 3; i++)
311 <                        resmat[j][i] =  solmat[j][0]*p0->n[i] +
312 <                                        solmat[j][1]*p1->n[i] +
313 <                                        solmat[j][2]*p2->n[i] +
314 <                                        solmat[j][3]*p3->n[i];
311 >                        resmat[j][i] =  eqnmat[j][0]*p0->n[i] +
312 >                                        eqnmat[j][1]*p1->n[i] +
313 >                                        eqnmat[j][2]*p2->n[i] +
314 >                                        eqnmat[j][3]*p3->n[i];
315          return(ax);
316  
317   #undef u
# Line 336 | Line 319 | POINT  *p0, *p1, *p2, *p3;
319   }
320  
321  
339 static double  m4tmp[4][4];             /* for efficiency */
340
341 #define  copymat4(m4a,m4b)      bcopy((char *)m4b,(char *)m4a,sizeof(m4tmp))
342
343
344 setident4(m4)
345 double  m4[4][4];
346 {
347        static double  ident[4][4] = {
348                1.,0.,0.,0.,
349                0.,1.,0.,0.,
350                0.,0.,1.,0.,
351                0.,0.,0.,1.,
352        };
353        copymat4(m4, ident);
354 }
355
322   /*
323   * invmat - computes the inverse of mat into inverse.  Returns 1
324   * if there exists an inverse, 0 otherwise.  It uses Gaussian Elimination
# Line 364 | Line 330 | double mat[4][4],inverse[4][4];
330   {
331   #define SWAP(a,b,t) (t=a,a=b,b=t)
332  
333 +        double  m4tmp[4][4];
334          register int i,j,k;
335          register double temp;
336  
337 <        setident4(inverse);
338 <        copymat4(m4tmp, mat);
337 >        bcopy(mat, m4tmp, sizeof(m4tmp));
338 >                                        /* set inverse to identity */
339 >        for (i = 0; i < 4; i++)
340 >                for (j = 0; j < 4; j++)
341 >                        inverse[i][j] = i==j ? 1.0 : 0.0;
342  
343          for(i = 0; i < 4; i++) {
344 <                if(m4tmp[i][i] == 0) {    /* Pivot is zero */
345 <                        /* Look for a raw with pivot != 0 and swap raws */
346 <                        for(j = i + 1; j < 4; j++)
347 <                                if(m4tmp[j][i] != 0) {
348 <                                        for( k = 0; k < 4; k++) {
349 <                                                SWAP(m4tmp[i][k],m4tmp[j][k],temp);
350 <                                                SWAP(inverse[i][k],inverse[j][k],temp);
351 <                                                }
352 <                                        break;
353 <                                        }
354 <                        if(j == 4)      /* No replacing raw -> no inverse */
355 <                                return(0);
356 <                        }
344 >                /* Look for raw with largest pivot and swap raws */
345 >                temp = FTINY; j = -1;
346 >                for(k = i; k < 4; k++)
347 >                        if(ABS(m4tmp[k][i]) > temp) {
348 >                                temp = ABS(m4tmp[k][i]);
349 >                                j = k;
350 >                                }
351 >                if(j == -1)     /* No replacing raw -> no inverse */
352 >                        return(0);
353 >                if (j != i)
354 >                        for(k = 0; k < 4; k++) {
355 >                                SWAP(m4tmp[i][k],m4tmp[j][k],temp);
356 >                                SWAP(inverse[i][k],inverse[j][k],temp);
357 >                                }
358  
359                  temp = m4tmp[i][i];
360                  for(k = 0; k < 4; k++) {
# Line 401 | Line 372 | double mat[4][4],inverse[4][4];
372                          }
373                  }
374          return(1);
375 +
376   #undef SWAP
377   }
378  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines