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

Comparing ray/src/rt/ambient.c (file contents):
Revision 1.13 by greg, Thu Jun 6 10:48:50 1991 UTC vs.
Revision 1.14 by greg, Fri Jun 7 10:04:29 1991 UTC

# Line 1 | Line 1
1 < /* Copyright (c) 1986 Regents of the University of California */
1 > /* Copyright (c) 1991 Regents of the University of California */
2  
3   #ifndef lint
4   static char SCCSid[] = "$SunId$ LBL";
# Line 19 | Line 19 | static char SCCSid[] = "$SunId$ LBL";
19  
20   #include  "otypes.h"
21  
22 + #include  "ambient.h"
23 +
24   #include  "random.h"
25  
26   #define  OCTSCALE       0.5     /* ceil((valid rad.)/(cube size)) */
27  
28 + typedef struct ambtree {
29 +        AMBVAL  *alist;         /* ambient value list */
30 +        struct ambtree  *kid;   /* 8 child nodes */
31 + }  AMBTREE;                     /* ambient octree */
32 +
33   extern CUBE  thescene;          /* contains space boundaries */
34  
28 extern COLOR  ambval;           /* global ambient component */
29 extern double  ambacc;          /* ambient accuracy */
30 extern int  ambres;             /* ambient resolution */
31 extern int  ambdiv;             /* number of divisions for calculation */
32 extern int  ambssamp;           /* number of super-samples */
33 extern int  ambounce;           /* number of ambient bounces */
34 extern char  *amblist[];        /* ambient include/exclude list */
35 extern int  ambincl;            /* include == 1, exclude == 0 */
36
35   #define  MAXASET        511     /* maximum number of elements in ambient set */
36   OBJECT  ambset[MAXASET+1]={0};  /* ambient include/exclude set */
37  
38   double  maxarad;                /* maximum ambient radius */
39   double  minarad;                /* minimum ambient radius */
40  
43 /*
44 * Since we've defined our vectors as float below to save space,
45 * watch out for changes in the definitions of VCOPY() and DOT().
46 */
47 typedef struct ambval {
48        float  pos[3];          /* position in space */
49        float  dir[3];          /* normal direction */
50        int  lvl;               /* recursion level of parent ray */
51        float  weight;          /* weight of parent ray */
52        COLOR  val;             /* computed ambient value */
53        float  rad;             /* validity radius */
54        struct ambval  *next;   /* next in list */
55 }  AMBVAL;                      /* ambient value */
56
57 typedef struct ambtree {
58        AMBVAL  *alist;         /* ambient value list */
59        struct ambtree  *kid;   /* 8 child nodes */
60 }  AMBTREE;                     /* ambient octree */
61
62 typedef struct {
63        float  k;               /* error contribution per sample */
64        COLOR  v;               /* ray sum */
65        int  n;                 /* number of samples */
66        short  t, p;            /* theta, phi indices */
67 }  AMBSAMP;                     /* ambient sample */
68
41   static AMBTREE  atrunk;         /* our ambient trunk node */
42  
43   static FILE  *ambfp = NULL;     /* ambient file pointer */
# Line 74 | Line 46 | static FILE  *ambfp = NULL;    /* ambient file pointer */
46  
47   #define  newambtree()   (AMBTREE *)calloc(8, sizeof(AMBTREE))
48  
77 double  sumambient(), doambient(), makeambient();
49  
79
50   setambient(afile)                       /* initialize calculation */
51   char  *afile;
52   {
# Line 145 | Line 115 | register RAY  *r;
115                  goto dumbamb;
116  
117          if (ambacc <= FTINY) {                  /* no ambient storage */
118 <                if (doambient(acol, r) == 0.0)
118 >                if (doambient(acol, r, NULL, NULL) == 0.0)
119                          goto dumbamb;
120                  goto done;
121          }
# Line 258 | Line 228 | COLOR  acol;
228   register RAY  *r;
229   {
230          AMBVAL  amb;
231 +        FVECT   gp, gd;
232  
233 <        amb.rad = doambient(acol, r);           /* compute ambient */
233 >        amb.rad = doambient(acol, r, gp, gd);   /* compute ambient */
234          if (amb.rad == 0.0)
235                  return(0.0);
236                                                  /* store it */
# Line 268 | Line 239 | register RAY  *r;
239          amb.lvl = r->rlvl;
240          amb.weight = r->rweight;
241          copycolor(amb.val, acol);
242 +        VCOPY(amb.gpos, gp);
243 +        VCOPY(amb.gdir, gd);
244                                                  /* insert into tree */
245          avinsert(&amb, &atrunk, thescene.cuorg, thescene.cusize);
246          avsave(&amb);                           /* write to file */
247          return(amb.rad);
275 }
276
277
278 double
279 doambient(acol, r)                      /* compute ambient component */
280 COLOR  acol;
281 register RAY  *r;
282 {
283        extern int  ambcmp();
284        extern double  sin(), cos(), sqrt();
285        int  hlist[4];
286        double  phi, xd, yd, zd;
287        double  b, b2;
288        register AMBSAMP  *div;
289        AMBSAMP  dnew;
290        RAY  ar;
291        FVECT  ux, uy;
292        double  arad;
293        int  ndivs, nt, np, ns, ne, i, j;
294        register int  k;
295
296        setcolor(acol, 0.0, 0.0, 0.0);
297                                        /* set number of divisions */
298        nt = sqrt(ambdiv * r->rweight * 0.5) + 0.5;
299        np = 2 * nt;
300        ndivs = nt * np;
301                                        /* check first */
302        if (ndivs == 0 || rayorigin(&ar, r, AMBIENT, 0.5) < 0)
303                return(0.0);
304                                        /* set number of super-samples */
305        ns = ambssamp * r->rweight + 0.5;
306        if (ns > 0) {
307                div = (AMBSAMP *)malloc(ndivs*sizeof(AMBSAMP));
308                if (div == NULL)
309                        error(SYSTEM, "out of memory in doambient");
310        } else
311                div = NULL;
312                                        /* make axes */
313        uy[0] = uy[1] = uy[2] = 0.0;
314        for (k = 0; k < 3; k++)
315                if (r->ron[k] < 0.6 && r->ron[k] > -0.6)
316                        break;
317        uy[k] = 1.0;
318        fcross(ux, r->ron, uy);
319        normalize(ux);
320        fcross(uy, ux, r->ron);
321                                        /* set up urand */
322        hlist[0] = r->rno;
323                                                /* sample divisions */
324        arad = 0.0;
325        ne = 0;
326        for (i = 0; i < nt; i++) {
327                hlist[1] = i;
328                for (j = 0; j < np; j++) {
329                        rayorigin(&ar, r, AMBIENT, 0.5);        /* pretested */
330                        hlist[2] = j;
331                        hlist[3] = 0;
332                        zd = sqrt((i+urand(ilhash(hlist,4)))/nt);
333                        hlist[3] = 1;
334                        phi = 2.0*PI * (j+urand(ilhash(hlist,4)))/np;
335                        xd = cos(phi) * zd;
336                        yd = sin(phi) * zd;
337                        zd = sqrt(1.0 - zd*zd);
338                        for (k = 0; k < 3; k++)
339                                ar.rdir[k] = xd*ux[k]+yd*uy[k]+zd*r->ron[k];
340                        dimlist[ndims++] = i*np + j + 38813;
341                        rayvalue(&ar);
342                        ndims--;
343                        if (ar.rot < FHUGE)
344                                arad += 1.0 / ar.rot;
345                        if (div != NULL) {              /* save division */
346                                div[ne].k = 0.0;
347                                copycolor(div[ne].v, ar.rcol);
348                                div[ne].n = 0;
349                                div[ne].t = i; div[ne].p = j;
350                                                        /* sum errors */
351                                b = bright(ar.rcol);
352                                if (i > 0) {            /* from above */
353                                        b2 = bright(div[ne-np].v) - b;
354                                        b2 *= b2 * 0.25;
355                                        div[ne].k += b2;
356                                        div[ne].n++;
357                                        div[ne-np].k += b2;
358                                        div[ne-np].n++;
359                                }
360                                if (j > 0) {            /* from behind */
361                                        b2 = bright(div[ne-1].v) - b;
362                                        b2 *= b2 * 0.25;
363                                        div[ne].k += b2;
364                                        div[ne].n++;
365                                        div[ne-1].k += b2;
366                                        div[ne-1].n++;
367                                }
368                                if (j == np-1) {        /* around */
369                                        b2 = bright(div[ne-(np-1)].v) - b;
370                                        b2 *= b2 * 0.25;
371                                        div[ne].k += b2;
372                                        div[ne].n++;
373                                        div[ne-(np-1)].k += b2;
374                                        div[ne-(np-1)].n++;
375                                }
376                                ne++;
377                        } else
378                                addcolor(acol, ar.rcol);
379                }
380        }
381        for (k = 0; k < ne; k++) {              /* compute errors */
382                if (div[k].n > 1)
383                        div[k].k /= div[k].n;
384                div[k].n = 1;
385        }
386                                                /* sort the divisions */
387        qsort(div, ne, sizeof(AMBSAMP), ambcmp);
388                                                /* skim excess */
389        while (ne > ns) {
390                ne--;
391                addcolor(acol, div[ne].v);
392        }
393                                                /* super-sample */
394        for (i = ns; i > 0; i--) {
395                rayorigin(&ar, r, AMBIENT, 0.5);        /* pretested */
396                hlist[1] = div[0].t;
397                hlist[2] = div[0].p;
398                hlist[3] = 0;
399                zd = sqrt((div[0].t+urand(ilhash(hlist,4)+div[0].n))/nt);
400                hlist[3] = 1;
401                phi = 2.0*PI * (div[0].p+urand(ilhash(hlist,4)+div[0].n))/np;
402                xd = cos(phi) * zd;
403                yd = sin(phi) * zd;
404                zd = sqrt(1.0 - zd*zd);
405                for (k = 0; k < 3; k++)
406                        ar.rdir[k] = xd*ux[k]+yd*uy[k]+zd*r->ron[k];
407                dimlist[ndims++] = div[0].t*np + div[0].p + 38813;
408                rayvalue(&ar);
409                ndims--;
410                rayvalue(&ar);
411                if (ar.rot < FHUGE)
412                        arad += 1.0 / ar.rot;
413                                                /* recompute error */
414                copycolor(dnew.v, div[0].v);
415                addcolor(dnew.v, ar.rcol);
416                dnew.n = div[0].n + 1;
417                dnew.t = div[0].t; dnew.p = div[0].p;
418                b2 = bright(dnew.v)/dnew.n - bright(ar.rcol);
419                b2 = b2*b2 + div[0].k*(div[0].n*div[0].n);
420                dnew.k = b2/(dnew.n*dnew.n);
421                                                /* reinsert */
422                for (k = 0; k < ne-1 && dnew.k < div[k+1].k; k++)
423                        copystruct(&div[k], &div[k+1]);
424                copystruct(&div[k], &dnew);
425
426                if (ne >= i) {          /* extract darkest division */
427                        ne--;
428                        if (div[ne].n > 1) {
429                                b = 1.0/div[ne].n;
430                                scalecolor(div[ne].v, b);
431                                div[ne].n = 1;
432                        }
433                        addcolor(acol, div[ne].v);
434                }
435        }
436        scalecolor(acol, 1.0/ndivs);
437        if (arad <= FTINY)
438                arad = FHUGE;
439        else
440                arad = (ndivs+ns) / arad / sqrt(r->rweight);
441        if (arad > maxarad)
442                arad = maxarad;
443        else if (arad < minarad)
444                arad = minarad;
445        if (div != NULL)
446                free((char *)div);
447        return(arad);
448 }
449
450
451 static int
452 ambcmp(d1, d2)                          /* decreasing order */
453 AMBSAMP  *d1, *d2;
454 {
455        if (d1->k < d2->k)
456                return(1);
457        if (d1->k > d2->k)
458                return(-1);
459        return(0);
248   }
249  
250  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines