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 2.80 by greg, Fri Apr 25 18:38:47 2014 UTC vs.
Revision 2.93 by greg, Fri Nov 21 00:53:52 2014 UTC

# Line 52 | Line 52 | static int  nunflshed = 0;     /* number of unflushed ambi
52   #endif
53  
54  
55 – static double  qambacc = 0.;            /* ambient accuracy to the 1/4 power */
55   static double  avsum = 0.;              /* computed ambient value sum (log) */
56   static unsigned int  navsum = 0;        /* number of values in avsum */
57   static unsigned int  nambvals = 0;      /* total number of indirect values */
# Line 110 | Line 109 | setambres(                             /* set ambient resolution */
109                                                  /* set min & max radii */
110          if (ar <= 0) {
111                  minarad = 0;
112 <                maxarad = thescene.cusize / 2.0;
112 >                maxarad = thescene.cusize*0.2;
113          } else {
114                  minarad = thescene.cusize / ar;
115 <                maxarad = 64 * minarad;                 /* heuristic */
116 <                if (maxarad > thescene.cusize / 2.0)
117 <                        maxarad = thescene.cusize / 2.0;
115 >                maxarad = 64.0 * minarad;               /* heuristic */
116 >                if (maxarad > thescene.cusize*0.2)
117 >                        maxarad = thescene.cusize*0.2;
118          }
119          if (minarad <= FTINY)
120 <                minarad = 10*FTINY;
120 >                minarad = 10.0*FTINY;
121          if (maxarad <= minarad)
122 <                maxarad = 64 * minarad;
122 >                maxarad = 64.0 * minarad;
123   }
124  
125  
# Line 129 | Line 128 | setambacc(                             /* set ambient accuracy */
128          double  newa
129   )
130   {
131 <        double  ambdiff;
132 <
133 <        if (newa < 0.0)
134 <                newa = 0.0;
135 <        ambdiff = fabs(newa - ambacc);
137 <        if (ambdiff >= .01 && (ambacc = newa) > FTINY) {
138 <                qambacc = sqrt(sqrt(ambacc));
131 >        static double   olda;           /* remember previous setting here */
132 >        
133 >        newa *= (newa > 0);
134 >        if (fabs(newa - olda) >= .05*(newa + olda)) {
135 >                ambacc = newa;
136                  if (nambvals > 0)
137                          sortambvals(1);         /* rebuild tree */
138          }
# Line 152 | Line 149 | setambient(void)                               /* initialize calculation */
149          ambdone();
150                                                  /* init ambient limits */
151          setambres(ambres);
152 <        qambacc = sqrt(sqrt(ambacc *= (ambacc > FTINY)));
152 >        setambacc(ambacc);
153          if (ambfile == NULL || !ambfile[0])
154                  return;
155          if (ambacc <= FTINY) {
# Line 168 | Line 165 | setambient(void)                               /* initialize calculation */
165                  initambfile(0);                 /* file exists */
166                  lastpos = ftell(ambfp);
167                  while (readambval(&amb, ambfp))
168 <                        avinsert(avstore(&amb));
168 >                        avstore(&amb);
169                  nambshare = nambvals;           /* share loaded values */
170                  if (readonly) {
171                          sprintf(errmsg,
# Line 266 | Line 263 | ambnotify(                     /* record new modifier */
263  
264   /************ THE FOLLOWING ROUTINES DIFFER BETWEEN NEW & OLD ***************/
265  
266 < #ifdef NEWAMB
266 > #ifndef OLDAMB
267  
268   #define tfunc(lwr, x, upr)      (((x)-(lwr))/((upr)-(lwr)))
269  
270 + static int      plugaleak(RAY *r, AMBVAL *ap, FVECT anorm, double ang);
271   static double   sumambient(COLOR acol, RAY *r, FVECT rn, int al,
272                                  AMBTREE *at, FVECT c0, double s);
273   static int      makeambient(COLOR acol, RAY *r, FVECT rn, int al);
274 < static void     extambient(COLOR cr, AMBVAL *ap, FVECT pv, FVECT nv,
274 > static int      extambient(COLOR cr, AMBVAL *ap, FVECT pv, FVECT nv,
275                                  FVECT uvw[3]);
276  
277   void
# Line 301 | Line 299 | multambient(           /* compute ambient component & multiply
299          if (ambacc <= FTINY) {                  /* no ambient storage */
300                  copycolor(acol, aval);
301                  rdepth++;
302 <                ok = doambient(acol, r, r->rweight, NULL, NULL, NULL, NULL);
302 >                ok = doambient(acol, r, r->rweight,
303 >                                NULL, NULL, NULL, NULL, NULL);
304                  rdepth--;
305                  if (!ok)
306                          goto dumbamb;
# Line 347 | Line 346 | dumbamb:                                       /* return global value */
346   }
347  
348  
349 < double
349 > /* Plug a potential leak where ambient cache value is occluded */
350 > static int
351 > plugaleak(RAY *r, AMBVAL *ap, FVECT anorm, double ang)
352 > {
353 >        const double    cost70sq = 0.1169778;   /* cos(70deg)^2 */
354 >        RAY             rtst;
355 >        FVECT           vdif;
356 >        double          normdot, ndotd, nadotd;
357 >        double          a, b, c, t[2];
358 >
359 >        ang += 2.*PI*(ang < 0);                 /* check direction flags */
360 >        if ( !(ap->corral>>(int)(ang*(16./PI)) & 1) )
361 >                return(0);
362 >        /*
363 >         * Generate test ray, targeting 20 degrees above sample point plane
364 >         * along surface normal from cache position.  This should be high
365 >         * enough to miss local geometry we don't really care about.
366 >         */
367 >        VSUB(vdif, ap->pos, r->rop);
368 >        normdot = DOT(anorm, r->ron);
369 >        ndotd = DOT(vdif, r->ron);
370 >        nadotd = DOT(vdif, anorm);
371 >        a = normdot*normdot - cost70sq;
372 >        b = 2.0*(normdot*ndotd - nadotd*cost70sq);
373 >        c = ndotd*ndotd - DOT(vdif,vdif)*cost70sq;
374 >        if (quadratic(t, a, b, c) != 2)
375 >                return(1);                      /* should rarely happen */
376 >        if (t[1] <= FTINY)
377 >                return(0);                      /* should fail behind test */
378 >        rayorigin(&rtst, SHADOW, r, NULL);
379 >        VSUM(rtst.rdir, vdif, anorm, t[1]);     /* further dist. > plane */
380 >        rtst.rmax = normalize(rtst.rdir);       /* short ray test */
381 >        while (localhit(&rtst, &thescene)) {    /* check for occluder */
382 >                if (rtst.ro->omod != OVOID &&
383 >                                (rtst.clipset == NULL ||
384 >                                        !inset(rtst.clipset, rtst.ro->omod)))
385 >                        return(1);              /* plug light leak */
386 >                VCOPY(rtst.rorg, rtst.rop);     /* skip invisible surface */
387 >                rtst.rmax -= rtst.rot;
388 >                rayclear(&rtst);
389 >        }
390 >        return(0);                              /* seems we're OK */
391 > }
392 >
393 >
394 > static double
395   sumambient(             /* get interpolated ambient value */
396          COLOR  acol,
397          RAY  *r,
# Line 359 | Line 403 | sumambient(            /* get interpolated ambient value */
403   )
404   {                       /* initial limit is 10 degrees plus ambacc radians */
405          const double    minangle = 10.0 * PI/180.;
406 <        const double    maxangle = (minangle+ambacc-PI/2.)*pow(r->rweight,0.13)
363 <                                        + PI/2.;
406 >        double          maxangle = minangle + ambacc;
407          double          wsum = 0.0;
408          FVECT           ck0;
409          int             i, j;
410          AMBVAL          *av;
411 +
412 +        if (at->kid != NULL) {          /* sum children first */                                
413 +                s *= 0.5;
414 +                for (i = 0; i < 8; i++) {
415 +                        for (j = 0; j < 3; j++) {
416 +                                ck0[j] = c0[j];
417 +                                if (1<<j & i)
418 +                                        ck0[j] += s;
419 +                                if (r->rop[j] < ck0[j] - OCTSCALE*s)
420 +                                        break;
421 +                                if (r->rop[j] > ck0[j] + (1.0+OCTSCALE)*s)
422 +                                        break;
423 +                        }
424 +                        if (j == 3)
425 +                                wsum += sumambient(acol, r, rn, al,
426 +                                                        at->kid+i, ck0, s);
427 +                }
428 +                                        /* good enough? */
429 +                if (wsum >= 0.05 && s > minarad*10.0)
430 +                        return(wsum);
431 +        }
432 +                                        /* adjust maximum angle */
433 +        if (at->alist != NULL && (at->alist->lvl <= al) & (r->rweight < 0.6))
434 +                maxangle = (maxangle - PI/2.)*pow(r->rweight,0.13) + PI/2.;
435                                          /* sum this node */
436          for (av = at->alist; av != NULL; av = av->next) {
437 <                double  d, delta_r2, delta_t2;
437 >                double  u, v, d, delta_r2, delta_t2;
438                  COLOR   ct;
439                  FVECT   uvw[3];
440                                          /* record access */
# Line 376 | Line 443 | sumambient(            /* get interpolated ambient value */
443                  /*
444                   *  Ambient level test
445                   */
446 <                if (av->lvl > al)       /* list sorted, so this works */
446 >                if (av->lvl > al ||     /* list sorted, so this works */
447 >                                (av->lvl == al) & (av->weight < 0.9*r->rweight))
448                          break;
381 –                if (av->weight < 0.9*r->rweight)
382 –                        continue;
449                  /*
450                   *  Direction test using unperturbed normal
451                   */
# Line 391 | Line 457 | sumambient(            /* get interpolated ambient value */
457                  if (delta_r2 >= maxangle*maxangle)
458                          continue;
459                  /*
460 +                 *  Modified ray behind test
461 +                 */
462 +                VSUB(ck0, r->rop, av->pos);
463 +                d = DOT(ck0, uvw[2]);
464 +                if (d < -minarad*ambacc-.001)
465 +                        continue;
466 +                d /= av->rad[0];
467 +                delta_t2 = d*d;
468 +                if (delta_t2 >= ambacc*ambacc)
469 +                        continue;
470 +                /*
471                   *  Elliptical radii test based on Hessian
472                   */
473                  decodedir(uvw[0], av->udir);
474                  VCROSS(uvw[1], uvw[2], uvw[0]);
475 <                VSUB(ck0, av->pos, r->rop);
399 <                d = DOT(ck0, uvw[0]) / av->rad[0];
400 <                delta_t2 = d*d;
401 <                d = DOT(ck0, uvw[1]) / av->rad[1];
475 >                d = (u = DOT(ck0, uvw[0])) / av->rad[0];
476                  delta_t2 += d*d;
477 <                d = DOT(ck0, uvw[2]) / av->rad[0];
477 >                d = (v = DOT(ck0, uvw[1])) / av->rad[1];
478                  delta_t2 += d*d;
479 <                if (delta_t2 >= qambacc*qambacc)
479 >                if (delta_t2 >= ambacc*ambacc)
480                          continue;
481                  /*
482 +                 *  Test for potential light leak
483 +                 */
484 +                if (av->corral && plugaleak(r, av, uvw[2], atan2a(v,u)))
485 +                        continue;
486 +                /*
487                   *  Extrapolate value and compute final weight (hat function)
488                   */
489 <                extambient(ct, av, r->rop, rn, uvw);
489 >                if (!extambient(ct, av, r->rop, rn, uvw))
490 >                        continue;
491                  d = tfunc(maxangle, sqrt(delta_r2), 0.0) *
492 <                        tfunc(qambacc, sqrt(delta_t2), 0.0);
492 >                        tfunc(ambacc, sqrt(delta_t2), 0.0);
493                  scalecolor(ct, d);
494                  addcolor(acol, ct);
495                  wsum += d;
496          }
417 –        if (at->kid == NULL)
418 –                return(wsum);
419 –                                        /* sum children */
420 –        s *= 0.5;
421 –        for (i = 0; i < 8; i++) {
422 –                for (j = 0; j < 3; j++) {
423 –                        ck0[j] = c0[j];
424 –                        if (1<<j & i)
425 –                                ck0[j] += s;
426 –                        if (r->rop[j] < ck0[j] - OCTSCALE*s)
427 –                                break;
428 –                        if (r->rop[j] > ck0[j] + (1.0+OCTSCALE)*s)
429 –                                break;
430 –                }
431 –                if (j == 3)
432 –                        wsum += sumambient(acol, r, rn, al,
433 –                                                at->kid+i, ck0, s);
434 –        }
497          return(wsum);
498   }
499  
500  
501 < int
501 > static int
502   makeambient(            /* make a new ambient value for storage */
503          COLOR  acol,
504          RAY  *r,
# Line 455 | Line 517 | makeambient(           /* make a new ambient value for storage
517                  amb.weight = 1.25*r->rweight;
518          setcolor(acol, AVGREFL, AVGREFL, AVGREFL);
519                                                  /* compute ambient */
520 <        i = doambient(acol, r, amb.weight, uvw, amb.rad, amb.gpos, amb.gdir);
520 >        i = doambient(acol, r, amb.weight,
521 >                        uvw, amb.rad, amb.gpos, amb.gdir, &amb.corral);
522          scalecolor(acol, 1./AVGREFL);           /* undo assumed reflectance */
523          if (i <= 0 || amb.rad[0] <= FTINY)      /* no Hessian or zero radius */
524                  return(i);
# Line 475 | Line 538 | makeambient(           /* make a new ambient value for storage
538   }
539  
540  
541 < void
541 > static int
542   extambient(             /* extrapolate value at pv, nv */
543          COLOR  cr,
544          AMBVAL   *ap,
# Line 484 | Line 547 | extambient(            /* extrapolate value at pv, nv */
547          FVECT  uvw[3]
548   )
549   {
550 +        const double    min_d = 0.05;
551          static FVECT    my_uvw[3];
552          FVECT           v1;
553          int             i;
# Line 503 | Line 567 | extambient(            /* extrapolate value at pv, nv */
567          for (i = 3; i--; )
568                  d += v1[i] * (ap->gdir[0]*uvw[0][i] + ap->gdir[1]*uvw[1][i]);
569          
570 <        if (d <= 0.0) {
571 <                setcolor(cr, 0.0, 0.0, 0.0);
508 <                return;
509 <        }
570 >        if (d < min_d)                  /* should not use if we can avoid it */
571 >                d = min_d;
572          copycolor(cr, ap->val);
573          scalecolor(cr, d);
574 +        return(d > min_d);
575   }
576  
577  
# Line 530 | Line 593 | avinsert(                              /* insert ambient value in our tree */
593          at = &atrunk;
594          VCOPY(ck0, thescene.cuorg);
595          s = thescene.cusize;
596 <        while (s*(OCTSCALE/2) > av->rad[1]*qambacc) {
596 >        while (s*(OCTSCALE/2) > av->rad[1]*ambacc) {
597                  if (at->kid == NULL)
598                          if ((at->kid = newambtree()) == NULL)
599                                  error(SYSTEM, "out of memory in avinsert");
# Line 545 | Line 608 | avinsert(                              /* insert ambient value in our tree */
608          }
609          avh.next = at->alist;           /* order by increasing level */
610          for (ap = &avh; ap->next != NULL; ap = ap->next)
611 <                if (ap->next->lvl >= av->lvl)
611 >                if ( ap->next->lvl > av->lvl ||
612 >                                (ap->next->lvl == av->lvl) &
613 >                                (ap->next->weight <= av->weight) )
614                          break;
615          av->next = ap->next;
616          ap->next = (AMBVAL*)av;
# Line 658 | Line 723 | sumambient(    /* get interpolated ambient value */
723                  /*
724                   *  Ambient level test.
725                   */
726 <                if (av->lvl > al)       /* list sorted, so this works */
726 >                if (av->lvl > al ||     /* list sorted, so this works */
727 >                                (av->lvl == al) & (av->weight < 0.9*r->rweight))
728                          break;
663 –                if (av->weight < 0.9*r->rweight)
664 –                        continue;
729                  /*
730                   *  Ambient radius test.
731                   */
# Line 843 | Line 907 | avinsert(                              /* insert ambient value in our tree */
907          }
908          avh.next = at->alist;           /* order by increasing level */
909          for (ap = &avh; ap->next != NULL; ap = ap->next)
910 <                if (ap->next->lvl >= av->lvl)
910 >                if ( ap->next->lvl > av->lvl ||
911 >                                (ap->next->lvl == av->lvl) &
912 >                                (ap->next->weight <= av->weight) )
913                          break;
914          av->next = ap->next;
915          ap->next = (AMBVAL*)av;
# Line 895 | Line 961 | avsave(                                /* insert and save an ambient value */
961          AMBVAL  *av
962   )
963   {
964 <        avinsert(avstore(av));
964 >        avstore(av);
965          if (ambfp == NULL)
966                  return;
967          if (writambval(av, ambfp) < 0)
# Line 910 | Line 976 | writerr:
976  
977  
978   static AMBVAL *
979 < avstore(                                /* allocate memory and store aval */
979 > avstore(                                /* allocate memory and save aval */
980          AMBVAL  *aval
981   )
982   {
# Line 928 | Line 994 | avstore(                               /* allocate memory and store aval */
994                  avsum += log(d);
995                  navsum++;
996          }
997 +        avinsert(av);                   /* insert in our cache tree */
998          return(av);
999   }
1000  
# Line 1128 | Line 1195 | sortambvals(                   /* resort ambient values */
1195                  if (i_avlist < nambvals)
1196                          error(CONSISTENCY, "missing ambient values in sortambvals");
1197   #endif
1198 <                qsort((char *)avlist1, nambvals, sizeof(struct avl), &alatcmp);
1199 <                qsort((char *)avlist2, nambvals, sizeof(AMBVAL *), &aposcmp);
1198 >                qsort((char *)avlist1, nambvals, sizeof(struct avl), alatcmp);
1199 >                qsort((char *)avlist2, nambvals, sizeof(AMBVAL *), aposcmp);
1200                  for (i = 0; i < nambvals; i++) {
1201                          if (avlist1[i].p == NULL)
1202                                  continue;
# Line 1210 | Line 1277 | ambsync(void)                  /* synchronize ambient file */
1277                                  error(WARNING, errmsg);
1278                                  break;
1279                          }
1280 <                        avinsert(avstore(&avs));
1280 >                        avstore(&avs);
1281                          n -= AMBVALSIZ;
1282                  }
1283                  lastpos = flen - n;

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines