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.83 by greg, Sat Apr 26 05:15:17 2014 UTC vs.
Revision 2.98 by greg, Sun Aug 23 00:17:12 2015 UTC

# Line 1 | Line 1
1 – #ifndef lint
1   static const char       RCSid[] = "$Id$";
3 – #endif
2   /*
3   *  ambient.c - routines dealing with ambient (inter-reflected) component.
4   *
# Line 17 | Line 15 | static const char      RCSid[] = "$Id$";
15   #include  "resolu.h"
16   #include  "ambient.h"
17   #include  "random.h"
18 + #include  "pmapamb.h"
19  
20   #ifndef  OCTSCALE
21   #define  OCTSCALE       1.0     /* ceil((valid rad.)/(cube size)) */
# Line 52 | Line 51 | static int  nunflshed = 0;     /* number of unflushed ambi
51   #endif
52  
53  
55 – static double  qambacc = 0.;            /* ambient accuracy to the 1/4 power */
54   static double  avsum = 0.;              /* computed ambient value sum (log) */
55   static unsigned int  navsum = 0;        /* number of values in avsum */
56   static unsigned int  nambvals = 0;      /* total number of indirect values */
# Line 110 | Line 108 | setambres(                             /* set ambient resolution */
108                                                  /* set min & max radii */
109          if (ar <= 0) {
110                  minarad = 0;
111 <                maxarad = thescene.cusize*0.5;
111 >                maxarad = thescene.cusize*0.2;
112          } else {
113                  minarad = thescene.cusize / ar;
114                  maxarad = 64.0 * minarad;               /* heuristic */
115 <                if (maxarad > thescene.cusize*0.5)
116 <                        maxarad = thescene.cusize*0.5;
115 >                if (maxarad > thescene.cusize*0.2)
116 >                        maxarad = thescene.cusize*0.2;
117          }
118          if (minarad <= FTINY)
119                  minarad = 10.0*FTINY;
# Line 129 | Line 127 | setambacc(                             /* set ambient accuracy */
127          double  newa
128   )
129   {
130 <        double  olda = qambacc*qambacc*qambacc*qambacc;
130 >        static double   olda;           /* remember previous setting here */
131          
132          newa *= (newa > 0);
133          if (fabs(newa - olda) >= .05*(newa + olda)) {
134 <                qambacc = sqrt(sqrt(ambacc = newa));
134 >                ambacc = newa;
135                  if (nambvals > 0)
136                          sortambvals(1);         /* rebuild tree */
137          }
# Line 166 | Line 164 | setambient(void)                               /* initialize calculation */
164                  initambfile(0);                 /* file exists */
165                  lastpos = ftell(ambfp);
166                  while (readambval(&amb, ambfp))
167 <                        avinsert(avstore(&amb));
167 >                        avstore(&amb);
168                  nambshare = nambvals;           /* share loaded values */
169                  if (readonly) {
170                          sprintf(errmsg,
# Line 198 | Line 196 | setambient(void)                               /* initialize calculation */
196                  sprintf(errmsg, "cannot open ambient file \"%s\"", ambfile);
197                  error(SYSTEM, errmsg);
198          }
201 – #ifdef getc_unlocked
202 –        flockfile(ambfp);                       /* application-level lock */
203 – #endif
199   #ifdef  F_SETLKW
200          aflock(F_UNLCK);                        /* release file */
201   #endif
# Line 264 | Line 259 | ambnotify(                     /* record new modifier */
259  
260   /************ THE FOLLOWING ROUTINES DIFFER BETWEEN NEW & OLD ***************/
261  
262 < #ifdef NEWAMB
262 > #ifndef OLDAMB
263  
264   #define tfunc(lwr, x, upr)      (((x)-(lwr))/((upr)-(lwr)))
265  
266 + static int      plugaleak(RAY *r, AMBVAL *ap, FVECT anorm, double ang);
267   static double   sumambient(COLOR acol, RAY *r, FVECT rn, int al,
268                                  AMBTREE *at, FVECT c0, double s);
269   static int      makeambient(COLOR acol, RAY *r, FVECT rn, int al);
270 < static void     extambient(COLOR cr, AMBVAL *ap, FVECT pv, FVECT nv,
270 > static int      extambient(COLOR cr, AMBVAL *ap, FVECT pv, FVECT nv,
271                                  FVECT uvw[3]);
272  
273   void
# Line 282 | Line 278 | multambient(           /* compute ambient component & multiply
278   )
279   {
280          static int  rdepth = 0;                 /* ambient recursion */
281 <        COLOR   acol;
281 >        COLOR   acol, caustic;
282          int     ok;
283          double  d, l;
284  
285 +        /* PMAP: Factor in ambient from photon map, if enabled and ray is
286 +         * ambient. Return as all ambient components accounted for, else
287 +         * continue. */
288 +        if (ambPmap(aval, r, rdepth))
289 +                return;
290 +
291 +        /* PMAP: Factor in specular-diffuse ambient (caustics) from photon
292 +         * map, if enabled and ray is primary, else caustic is zero.  Continue
293 +         * with RADIANCE ambient calculation */
294 +        copycolor(caustic, aval);
295 +        ambPmapCaustic(caustic, r, rdepth);
296 +        
297          if (ambdiv <= 0)                        /* no ambient calculation */
298                  goto dumbamb;
299                                                  /* check number of bounces */
# Line 299 | Line 307 | multambient(           /* compute ambient component & multiply
307          if (ambacc <= FTINY) {                  /* no ambient storage */
308                  copycolor(acol, aval);
309                  rdepth++;
310 <                ok = doambient(acol, r, r->rweight, NULL, NULL, NULL, NULL);
310 >                ok = doambient(acol, r, r->rweight,
311 >                                NULL, NULL, NULL, NULL, NULL);
312                  rdepth--;
313                  if (!ok)
314                          goto dumbamb;
315                  copycolor(aval, acol);
316 +
317 +                /* PMAP: add in caustic */
318 +                addcolor(aval, caustic);
319                  return;
320          }
321  
# Line 313 | Line 325 | multambient(           /* compute ambient component & multiply
325          setcolor(acol, 0.0, 0.0, 0.0);
326          d = sumambient(acol, r, nrm, rdepth,
327                          &atrunk, thescene.cuorg, thescene.cusize);
328 +                        
329          if (d > FTINY) {
330                  d = 1.0/d;
331                  scalecolor(acol, d);
332                  multcolor(aval, acol);
333 +
334 +                /* PMAP: add in caustic */
335 +                addcolor(aval, caustic);
336                  return;
337          }
338 +        
339          rdepth++;                               /* need to cache new value */
340          ok = makeambient(acol, r, nrm, rdepth-1);
341          rdepth--;
342 +        
343          if (ok) {
344                  multcolor(aval, acol);          /* computed new value */
345 +
346 +                /* PMAP: add in caustic */
347 +                addcolor(aval, caustic);
348                  return;
349          }
350 +        
351   dumbamb:                                        /* return global value */
352          if ((ambvwt <= 0) | (navsum == 0)) {
353                  multcolor(aval, ambval);
354 +                
355 +                /* PMAP: add in caustic */
356 +                addcolor(aval, caustic);
357                  return;
358          }
359 <        l = bright(ambval);                     /* average in computations */
359 >        
360 >        l = bright(ambval);                     /* average in computations */  
361          if (l > FTINY) {
362                  d = (log(l)*(double)ambvwt + avsum) /
363                                  (double)(ambvwt + navsum);
# Line 345 | Line 371 | dumbamb:                                       /* return global value */
371   }
372  
373  
374 < double
374 > /* Plug a potential leak where ambient cache value is occluded */
375 > static int
376 > plugaleak(RAY *r, AMBVAL *ap, FVECT anorm, double ang)
377 > {
378 >        const double    cost70sq = 0.1169778;   /* cos(70deg)^2 */
379 >        RAY             rtst;
380 >        FVECT           vdif;
381 >        double          normdot, ndotd, nadotd;
382 >        double          a, b, c, t[2];
383 >
384 >        ang += 2.*PI*(ang < 0);                 /* check direction flags */
385 >        if ( !(ap->corral>>(int)(ang*(16./PI)) & 1) )
386 >                return(0);
387 >        /*
388 >         * Generate test ray, targeting 20 degrees above sample point plane
389 >         * along surface normal from cache position.  This should be high
390 >         * enough to miss local geometry we don't really care about.
391 >         */
392 >        VSUB(vdif, ap->pos, r->rop);
393 >        normdot = DOT(anorm, r->ron);
394 >        ndotd = DOT(vdif, r->ron);
395 >        nadotd = DOT(vdif, anorm);
396 >        a = normdot*normdot - cost70sq;
397 >        b = 2.0*(normdot*ndotd - nadotd*cost70sq);
398 >        c = ndotd*ndotd - DOT(vdif,vdif)*cost70sq;
399 >        if (quadratic(t, a, b, c) != 2)
400 >                return(1);                      /* should rarely happen */
401 >        if (t[1] <= FTINY)
402 >                return(0);                      /* should fail behind test */
403 >        rayorigin(&rtst, SHADOW, r, NULL);
404 >        VSUM(rtst.rdir, vdif, anorm, t[1]);     /* further dist. > plane */
405 >        rtst.rmax = normalize(rtst.rdir);       /* short ray test */
406 >        while (localhit(&rtst, &thescene)) {    /* check for occluder */
407 >                if (rtst.ro->omod != OVOID &&
408 >                                (rtst.clipset == NULL ||
409 >                                        !inset(rtst.clipset, rtst.ro->omod)))
410 >                        return(1);              /* plug light leak */
411 >                VCOPY(rtst.rorg, rtst.rop);     /* skip invisible surface */
412 >                rtst.rmax -= rtst.rot;
413 >                rayclear(&rtst);
414 >        }
415 >        return(0);                              /* seems we're OK */
416 > }
417 >
418 >
419 > static double
420   sumambient(             /* get interpolated ambient value */
421          COLOR  acol,
422          RAY  *r,
# Line 357 | Line 428 | sumambient(            /* get interpolated ambient value */
428   )
429   {                       /* initial limit is 10 degrees plus ambacc radians */
430          const double    minangle = 10.0 * PI/180.;
431 <        const double    maxangle = (minangle+ambacc-PI/2.)*pow(r->rweight,0.13)
361 <                                        + PI/2.;
431 >        double          maxangle = minangle + ambacc;
432          double          wsum = 0.0;
433          FVECT           ck0;
434          int             i, j;
# Line 381 | Line 451 | sumambient(            /* get interpolated ambient value */
451                                                          at->kid+i, ck0, s);
452                  }
453                                          /* good enough? */
454 <                if (wsum > 0.04 && s > (minarad*0.8+maxarad*0.2))
454 >                if (wsum >= 0.05 && s > minarad*10.0)
455                          return(wsum);
456          }
457 +                                        /* adjust maximum angle */
458 +        if (at->alist != NULL && (at->alist->lvl <= al) & (r->rweight < 0.6))
459 +                maxangle = (maxangle - PI/2.)*pow(r->rweight,0.13) + PI/2.;
460                                          /* sum this node */
461          for (av = at->alist; av != NULL; av = av->next) {
462 <                double  d, delta_r2, delta_t2;
462 >                double  u, v, d, delta_r2, delta_t2;
463                  COLOR   ct;
464                  FVECT   uvw[3];
465                                          /* record access */
# Line 395 | Line 468 | sumambient(            /* get interpolated ambient value */
468                  /*
469                   *  Ambient level test
470                   */
471 <                if (av->lvl > al)       /* list sorted, so this works */
471 >                if (av->lvl > al ||     /* list sorted, so this works */
472 >                                (av->lvl == al) & (av->weight < 0.9*r->rweight))
473                          break;
400 –                if (av->weight < 0.9*r->rweight)
401 –                        continue;
474                  /*
475                   *  Direction test using unperturbed normal
476                   */
# Line 412 | Line 484 | sumambient(            /* get interpolated ambient value */
484                  /*
485                   *  Modified ray behind test
486                   */
487 <                VSUB(ck0, av->pos, r->rop);
487 >                VSUB(ck0, r->rop, av->pos);
488                  d = DOT(ck0, uvw[2]);
489 <                if (d < -minarad*qambacc-.001)
489 >                if (d < -minarad*ambacc-.001)
490                          continue;
491                  d /= av->rad[0];
492                  delta_t2 = d*d;
493 <                if (delta_t2 >= qambacc*qambacc)
493 >                if (delta_t2 >= ambacc*ambacc)
494                          continue;
495                  /*
496                   *  Elliptical radii test based on Hessian
497                   */
498                  decodedir(uvw[0], av->udir);
499                  VCROSS(uvw[1], uvw[2], uvw[0]);
500 <                d = DOT(ck0, uvw[0]) / av->rad[0];
500 >                d = (u = DOT(ck0, uvw[0])) / av->rad[0];
501                  delta_t2 += d*d;
502 <                d = DOT(ck0, uvw[1]) / av->rad[1];
502 >                d = (v = DOT(ck0, uvw[1])) / av->rad[1];
503                  delta_t2 += d*d;
504 <                if (delta_t2 >= qambacc*qambacc)
504 >                if (delta_t2 >= ambacc*ambacc)
505                          continue;
506                  /*
507 +                 *  Test for potential light leak
508 +                 */
509 +                if (av->corral && plugaleak(r, av, uvw[2], atan2a(v,u)))
510 +                        continue;
511 +                /*
512                   *  Extrapolate value and compute final weight (hat function)
513                   */
514 <                extambient(ct, av, r->rop, rn, uvw);
514 >                if (!extambient(ct, av, r->rop, rn, uvw))
515 >                        continue;
516                  d = tfunc(maxangle, sqrt(delta_r2), 0.0) *
517 <                        tfunc(qambacc, sqrt(delta_t2), 0.0);
517 >                        tfunc(ambacc, sqrt(delta_t2), 0.0);
518                  scalecolor(ct, d);
519                  addcolor(acol, ct);
520                  wsum += d;
# Line 445 | Line 523 | sumambient(            /* get interpolated ambient value */
523   }
524  
525  
526 < int
526 > static int
527   makeambient(            /* make a new ambient value for storage */
528          COLOR  acol,
529          RAY  *r,
# Line 464 | Line 542 | makeambient(           /* make a new ambient value for storage
542                  amb.weight = 1.25*r->rweight;
543          setcolor(acol, AVGREFL, AVGREFL, AVGREFL);
544                                                  /* compute ambient */
545 <        i = doambient(acol, r, amb.weight, uvw, amb.rad, amb.gpos, amb.gdir);
545 >        i = doambient(acol, r, amb.weight,
546 >                        uvw, amb.rad, amb.gpos, amb.gdir, &amb.corral);
547          scalecolor(acol, 1./AVGREFL);           /* undo assumed reflectance */
548          if (i <= 0 || amb.rad[0] <= FTINY)      /* no Hessian or zero radius */
549                  return(i);
# Line 484 | Line 563 | makeambient(           /* make a new ambient value for storage
563   }
564  
565  
566 < void
566 > static int
567   extambient(             /* extrapolate value at pv, nv */
568          COLOR  cr,
569          AMBVAL   *ap,
# Line 493 | Line 572 | extambient(            /* extrapolate value at pv, nv */
572          FVECT  uvw[3]
573   )
574   {
575 +        const double    min_d = 0.05;
576          static FVECT    my_uvw[3];
577          FVECT           v1;
578          int             i;
# Line 512 | Line 592 | extambient(            /* extrapolate value at pv, nv */
592          for (i = 3; i--; )
593                  d += v1[i] * (ap->gdir[0]*uvw[0][i] + ap->gdir[1]*uvw[1][i]);
594          
595 <        if (d <= 0.0) {
596 <                setcolor(cr, 0.0, 0.0, 0.0);
517 <                return;
518 <        }
595 >        if (d < min_d)                  /* should not use if we can avoid it */
596 >                d = min_d;
597          copycolor(cr, ap->val);
598          scalecolor(cr, d);
599 +        return(d > min_d);
600   }
601  
602  
# Line 539 | Line 618 | avinsert(                              /* insert ambient value in our tree */
618          at = &atrunk;
619          VCOPY(ck0, thescene.cuorg);
620          s = thescene.cusize;
621 <        while (s*(OCTSCALE/2) > av->rad[1]*qambacc) {
621 >        while (s*(OCTSCALE/2) > av->rad[1]*ambacc) {
622                  if (at->kid == NULL)
623                          if ((at->kid = newambtree()) == NULL)
624                                  error(SYSTEM, "out of memory in avinsert");
# Line 554 | Line 633 | avinsert(                              /* insert ambient value in our tree */
633          }
634          avh.next = at->alist;           /* order by increasing level */
635          for (ap = &avh; ap->next != NULL; ap = ap->next)
636 <                if (ap->next->lvl >= av->lvl)
636 >                if ( ap->next->lvl > av->lvl ||
637 >                                (ap->next->lvl == av->lvl) &
638 >                                (ap->next->weight <= av->weight) )
639                          break;
640          av->next = ap->next;
641          ap->next = (AMBVAL*)av;
# Line 578 | Line 659 | multambient(           /* compute ambient component & multiply
659   )
660   {
661          static int  rdepth = 0;                 /* ambient recursion */
662 <        COLOR   acol;
662 >        COLOR   acol, caustic;
663          double  d, l;
664  
665 +        /* PMAP: Factor in ambient from global photon map (if enabled) and return
666 +         * as all ambient components accounted for */
667 +        if (ambPmap(aval, r, rdepth))
668 +                return;
669 +
670 +        /* PMAP: Otherwise factor in ambient from caustic photon map
671 +         * (ambPmapCaustic() returns zero if caustic photons disabled) and
672 +         * continue with RADIANCE ambient calculation */
673 +        copycolor(caustic, aval);
674 +        ambPmapCaustic(caustic, r, rdepth);
675 +        
676          if (ambdiv <= 0)                        /* no ambient calculation */
677                  goto dumbamb;
678                                                  /* check number of bounces */
# Line 598 | Line 690 | multambient(           /* compute ambient component & multiply
690                  rdepth--;
691                  if (d <= FTINY)
692                          goto dumbamb;
693 <                copycolor(aval, acol);
693 >                copycolor(aval, acol);          
694 >        
695 >           /* PMAP: add in caustic */
696 >                addcolor(aval, caustic);        
697                  return;
698          }
699  
# Line 608 | Line 703 | multambient(           /* compute ambient component & multiply
703          setcolor(acol, 0.0, 0.0, 0.0);
704          d = sumambient(acol, r, nrm, rdepth,
705                          &atrunk, thescene.cuorg, thescene.cusize);
706 +                        
707          if (d > FTINY) {
708                  d = 1.0/d;
709                  scalecolor(acol, d);
710                  multcolor(aval, acol);
711 +                
712 +                /* PMAP: add in caustic */
713 +                addcolor(aval, caustic);        
714                  return;
715          }
716 +        
717          rdepth++;                               /* need to cache new value */
718          d = makeambient(acol, r, nrm, rdepth-1);
719          rdepth--;
720 +        
721          if (d > FTINY) {
722                  multcolor(aval, acol);          /* got new value */
723 +
724 +                /* PMAP: add in caustic */
725 +                addcolor(aval, caustic);                        
726                  return;
727          }
728 +        
729   dumbamb:                                        /* return global value */
730          if ((ambvwt <= 0) | (navsum == 0)) {
731                  multcolor(aval, ambval);
732 +
733 +                /* PMAP: add in caustic */
734 +                addcolor(aval, caustic);        
735                  return;
736          }
737 +        
738          l = bright(ambval);                     /* average in computations */
739          if (l > FTINY) {
740                  d = (log(l)*(double)ambvwt + avsum) /
# Line 667 | Line 776 | sumambient(    /* get interpolated ambient value */
776                  /*
777                   *  Ambient level test.
778                   */
779 <                if (av->lvl > al)       /* list sorted, so this works */
779 >                if (av->lvl > al ||     /* list sorted, so this works */
780 >                                (av->lvl == al) & (av->weight < 0.9*r->rweight))
781                          break;
672 –                if (av->weight < 0.9*r->rweight)
673 –                        continue;
782                  /*
783                   *  Ambient radius test.
784                   */
# Line 852 | Line 960 | avinsert(                              /* insert ambient value in our tree */
960          }
961          avh.next = at->alist;           /* order by increasing level */
962          for (ap = &avh; ap->next != NULL; ap = ap->next)
963 <                if (ap->next->lvl >= av->lvl)
963 >                if ( ap->next->lvl > av->lvl ||
964 >                                (ap->next->lvl == av->lvl) &
965 >                                (ap->next->weight <= av->weight) )
966                          break;
967          av->next = ap->next;
968          ap->next = (AMBVAL*)av;
# Line 904 | Line 1014 | avsave(                                /* insert and save an ambient value */
1014          AMBVAL  *av
1015   )
1016   {
1017 <        avinsert(avstore(av));
1017 >        avstore(av);
1018          if (ambfp == NULL)
1019                  return;
1020          if (writambval(av, ambfp) < 0)
# Line 919 | Line 1029 | writerr:
1029  
1030  
1031   static AMBVAL *
1032 < avstore(                                /* allocate memory and store aval */
1032 > avstore(                                /* allocate memory and save aval */
1033          AMBVAL  *aval
1034   )
1035   {
# Line 937 | Line 1047 | avstore(                               /* allocate memory and store aval */
1047                  avsum += log(d);
1048                  navsum++;
1049          }
1050 +        avinsert(av);                   /* insert in our cache tree */
1051          return(av);
1052   }
1053  
# Line 1137 | Line 1248 | sortambvals(                   /* resort ambient values */
1248                  if (i_avlist < nambvals)
1249                          error(CONSISTENCY, "missing ambient values in sortambvals");
1250   #endif
1251 <                qsort((char *)avlist1, nambvals, sizeof(struct avl), &alatcmp);
1252 <                qsort((char *)avlist2, nambvals, sizeof(AMBVAL *), &aposcmp);
1251 >                qsort((char *)avlist1, nambvals, sizeof(struct avl), alatcmp);
1252 >                qsort((char *)avlist2, nambvals, sizeof(AMBVAL *), aposcmp);
1253                  for (i = 0; i < nambvals; i++) {
1254                          if (avlist1[i].p == NULL)
1255                                  continue;
# Line 1219 | Line 1330 | ambsync(void)                  /* synchronize ambient file */
1330                                  error(WARNING, errmsg);
1331                                  break;
1332                          }
1333 <                        avinsert(avstore(&avs));
1333 >                        avstore(&avs);
1334                          n -= AMBVALSIZ;
1335                  }
1336                  lastpos = flen - n;

Diff Legend

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