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.72 by greg, Fri Apr 11 22:54:34 2014 UTC vs.
Revision 2.97 by greg, Fri Aug 21 18:21:05 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 51 | Line 50 | static int  nunflshed = 0;     /* number of unflushed ambi
50   #define MAX_SORT_INTVL  (SORT_INTVL<<6)
51   #endif
52  
53 +
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 108 | Line 108 | setambres(                             /* set ambient resolution */
108                                                  /* set min & max radii */
109          if (ar <= 0) {
110                  minarad = 0;
111 <                maxarad = thescene.cusize / 2.0;
111 >                maxarad = thescene.cusize*0.2;
112          } else {
113                  minarad = thescene.cusize / ar;
114 <                maxarad = 64 * minarad;                 /* heuristic */
115 <                if (maxarad > thescene.cusize / 2.0)
116 <                        maxarad = thescene.cusize / 2.0;
114 >                maxarad = 64.0 * minarad;               /* heuristic */
115 >                if (maxarad > thescene.cusize*0.2)
116 >                        maxarad = thescene.cusize*0.2;
117          }
118          if (minarad <= FTINY)
119 <                minarad = 10*FTINY;
119 >                minarad = 10.0*FTINY;
120          if (maxarad <= minarad)
121 <                maxarad = 64 * minarad;
121 >                maxarad = 64.0 * minarad;
122   }
123  
124  
# Line 127 | Line 127 | setambacc(                             /* set ambient accuracy */
127          double  newa
128   )
129   {
130 <        double  ambdiff;
131 <
132 <        if (newa < 0.0)
133 <                newa = 0.0;
134 <        ambdiff = fabs(newa - ambacc);
135 <        if (ambdiff >= .01 && (ambacc = newa) > FTINY && nambvals > 0)
136 <                sortambvals(1);                 /* rebuild tree */
130 >        static double   olda;           /* remember previous setting here */
131 >        
132 >        newa *= (newa > 0);
133 >        if (fabs(newa - olda) >= .05*(newa + olda)) {
134 >                ambacc = newa;
135 >                if (nambvals > 0)
136 >                        sortambvals(1);         /* rebuild tree */
137 >        }
138   }
139  
140  
# Line 163 | 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 195 | Line 196 | setambient(void)                               /* initialize calculation */
196                  sprintf(errmsg, "cannot open ambient file \"%s\"", ambfile);
197                  error(SYSTEM, errmsg);
198          }
199 < #ifdef getc_unlocked
199 > #if 0
200          flockfile(ambfp);                       /* application-level lock */
201   #endif
202   #ifdef  F_SETLKW
# Line 261 | Line 262 | ambnotify(                     /* record new modifier */
262  
263   /************ THE FOLLOWING ROUTINES DIFFER BETWEEN NEW & OLD ***************/
264  
265 < #ifdef NEWAMB
265 > #ifndef OLDAMB
266  
267   #define tfunc(lwr, x, upr)      (((x)-(lwr))/((upr)-(lwr)))
268  
269 + static int      plugaleak(RAY *r, AMBVAL *ap, FVECT anorm, double ang);
270   static double   sumambient(COLOR acol, RAY *r, FVECT rn, int al,
271                                  AMBTREE *at, FVECT c0, double s);
272   static int      makeambient(COLOR acol, RAY *r, FVECT rn, int al);
273 < static void     extambient(COLOR cr, AMBVAL *ap, FVECT pv, FVECT nv,
273 > static int      extambient(COLOR cr, AMBVAL *ap, FVECT pv, FVECT nv,
274                                  FVECT uvw[3]);
275  
276   void
# Line 279 | Line 281 | multambient(           /* compute ambient component & multiply
281   )
282   {
283          static int  rdepth = 0;                 /* ambient recursion */
284 <        COLOR   acol;
284 >        COLOR   acol, caustic;
285          int     ok;
286          double  d, l;
287  
288 +        /* PMAP: Factor in ambient from photon map, if enabled and ray is
289 +         * ambient. Return as all ambient components accounted for, else
290 +         * continue. */
291 +        if (ambPmap(aval, r, rdepth))
292 +                return;
293 +
294 +        /* PMAP: Factor in specular-diffuse ambient (caustics) from photon
295 +         * map, if enabled and ray is primary, else caustic is zero.  Continue
296 +         * with RADIANCE ambient calculation */
297 +        copycolor(caustic, aval);
298 +        ambPmapCaustic(caustic, r, rdepth);
299 +        
300          if (ambdiv <= 0)                        /* no ambient calculation */
301                  goto dumbamb;
302                                                  /* check number of bounces */
# Line 296 | Line 310 | multambient(           /* compute ambient component & multiply
310          if (ambacc <= FTINY) {                  /* no ambient storage */
311                  copycolor(acol, aval);
312                  rdepth++;
313 <                ok = doambient(acol, r, r->rweight, NULL, NULL, NULL, NULL);
313 >                ok = doambient(acol, r, r->rweight,
314 >                                NULL, NULL, NULL, NULL, NULL);
315                  rdepth--;
316                  if (!ok)
317                          goto dumbamb;
318                  copycolor(aval, acol);
319 +
320 +                /* PMAP: add in caustic */
321 +                addcolor(aval, caustic);
322                  return;
323          }
324  
# Line 310 | Line 328 | multambient(           /* compute ambient component & multiply
328          setcolor(acol, 0.0, 0.0, 0.0);
329          d = sumambient(acol, r, nrm, rdepth,
330                          &atrunk, thescene.cuorg, thescene.cusize);
331 +                        
332          if (d > FTINY) {
333                  d = 1.0/d;
334                  scalecolor(acol, d);
335                  multcolor(aval, acol);
336 +
337 +                /* PMAP: add in caustic */
338 +                addcolor(aval, caustic);
339                  return;
340          }
341 +        
342          rdepth++;                               /* need to cache new value */
343          ok = makeambient(acol, r, nrm, rdepth-1);
344          rdepth--;
345 +        
346          if (ok) {
347 <                multcolor(aval, acol);          /* got new value */
347 >                multcolor(aval, acol);          /* computed new value */
348 >
349 >                /* PMAP: add in caustic */
350 >                addcolor(aval, caustic);
351                  return;
352          }
353 +        
354   dumbamb:                                        /* return global value */
355          if ((ambvwt <= 0) | (navsum == 0)) {
356                  multcolor(aval, ambval);
357 +                
358 +                /* PMAP: add in caustic */
359 +                addcolor(aval, caustic);
360                  return;
361          }
362 <        l = bright(ambval);                     /* average in computations */
362 >        
363 >        l = bright(ambval);                     /* average in computations */  
364          if (l > FTINY) {
365                  d = (log(l)*(double)ambvwt + avsum) /
366                                  (double)(ambvwt + navsum);
# Line 342 | Line 374 | dumbamb:                                       /* return global value */
374   }
375  
376  
377 < double
377 > /* Plug a potential leak where ambient cache value is occluded */
378 > static int
379 > plugaleak(RAY *r, AMBVAL *ap, FVECT anorm, double ang)
380 > {
381 >        const double    cost70sq = 0.1169778;   /* cos(70deg)^2 */
382 >        RAY             rtst;
383 >        FVECT           vdif;
384 >        double          normdot, ndotd, nadotd;
385 >        double          a, b, c, t[2];
386 >
387 >        ang += 2.*PI*(ang < 0);                 /* check direction flags */
388 >        if ( !(ap->corral>>(int)(ang*(16./PI)) & 1) )
389 >                return(0);
390 >        /*
391 >         * Generate test ray, targeting 20 degrees above sample point plane
392 >         * along surface normal from cache position.  This should be high
393 >         * enough to miss local geometry we don't really care about.
394 >         */
395 >        VSUB(vdif, ap->pos, r->rop);
396 >        normdot = DOT(anorm, r->ron);
397 >        ndotd = DOT(vdif, r->ron);
398 >        nadotd = DOT(vdif, anorm);
399 >        a = normdot*normdot - cost70sq;
400 >        b = 2.0*(normdot*ndotd - nadotd*cost70sq);
401 >        c = ndotd*ndotd - DOT(vdif,vdif)*cost70sq;
402 >        if (quadratic(t, a, b, c) != 2)
403 >                return(1);                      /* should rarely happen */
404 >        if (t[1] <= FTINY)
405 >                return(0);                      /* should fail behind test */
406 >        rayorigin(&rtst, SHADOW, r, NULL);
407 >        VSUM(rtst.rdir, vdif, anorm, t[1]);     /* further dist. > plane */
408 >        rtst.rmax = normalize(rtst.rdir);       /* short ray test */
409 >        while (localhit(&rtst, &thescene)) {    /* check for occluder */
410 >                if (rtst.ro->omod != OVOID &&
411 >                                (rtst.clipset == NULL ||
412 >                                        !inset(rtst.clipset, rtst.ro->omod)))
413 >                        return(1);              /* plug light leak */
414 >                VCOPY(rtst.rorg, rtst.rop);     /* skip invisible surface */
415 >                rtst.rmax -= rtst.rot;
416 >                rayclear(&rtst);
417 >        }
418 >        return(0);                              /* seems we're OK */
419 > }
420 >
421 >
422 > static double
423   sumambient(             /* get interpolated ambient value */
424          COLOR  acol,
425          RAY  *r,
# Line 352 | Line 429 | sumambient(            /* get interpolated ambient value */
429          FVECT  c0,
430          double  s
431   )
432 < {                                       /* initial limit is ambacc radians */
433 <        const double    maxangle = (ambacc-PI/2.)*pow(r->rweight,0.13) + PI/2.;
432 > {                       /* initial limit is 10 degrees plus ambacc radians */
433 >        const double    minangle = 10.0 * PI/180.;
434 >        double          maxangle = minangle + ambacc;
435          double          wsum = 0.0;
436          FVECT           ck0;
437          int             i, j;
438          AMBVAL          *av;
439 +
440 +        if (at->kid != NULL) {          /* sum children first */                                
441 +                s *= 0.5;
442 +                for (i = 0; i < 8; i++) {
443 +                        for (j = 0; j < 3; j++) {
444 +                                ck0[j] = c0[j];
445 +                                if (1<<j & i)
446 +                                        ck0[j] += s;
447 +                                if (r->rop[j] < ck0[j] - OCTSCALE*s)
448 +                                        break;
449 +                                if (r->rop[j] > ck0[j] + (1.0+OCTSCALE)*s)
450 +                                        break;
451 +                        }
452 +                        if (j == 3)
453 +                                wsum += sumambient(acol, r, rn, al,
454 +                                                        at->kid+i, ck0, s);
455 +                }
456 +                                        /* good enough? */
457 +                if (wsum >= 0.05 && s > minarad*10.0)
458 +                        return(wsum);
459 +        }
460 +                                        /* adjust maximum angle */
461 +        if (at->alist != NULL && (at->alist->lvl <= al) & (r->rweight < 0.6))
462 +                maxangle = (maxangle - PI/2.)*pow(r->rweight,0.13) + PI/2.;
463                                          /* sum this node */
464          for (av = at->alist; av != NULL; av = av->next) {
465 <                double  d, delta_r2, delta_t2;
465 >                double  u, v, d, delta_r2, delta_t2;
466                  COLOR   ct;
467                  FVECT   uvw[3];
468                                          /* record access */
# Line 369 | Line 471 | sumambient(            /* get interpolated ambient value */
471                  /*
472                   *  Ambient level test
473                   */
474 <                if (av->lvl > al)       /* list sorted, so this works */
474 >                if (av->lvl > al ||     /* list sorted, so this works */
475 >                                (av->lvl == al) & (av->weight < 0.9*r->rweight))
476                          break;
374 –                if (av->weight < 0.9*r->rweight)
375 –                        continue;
477                  /*
478                   *  Direction test using unperturbed normal
479                   */
# Line 384 | Line 485 | sumambient(            /* get interpolated ambient value */
485                  if (delta_r2 >= maxangle*maxangle)
486                          continue;
487                  /*
488 +                 *  Modified ray behind test
489 +                 */
490 +                VSUB(ck0, r->rop, av->pos);
491 +                d = DOT(ck0, uvw[2]);
492 +                if (d < -minarad*ambacc-.001)
493 +                        continue;
494 +                d /= av->rad[0];
495 +                delta_t2 = d*d;
496 +                if (delta_t2 >= ambacc*ambacc)
497 +                        continue;
498 +                /*
499                   *  Elliptical radii test based on Hessian
500                   */
501                  decodedir(uvw[0], av->udir);
502                  VCROSS(uvw[1], uvw[2], uvw[0]);
503 <                VSUB(ck0, av->pos, r->rop);
392 <                d = DOT(ck0, uvw[0]) / av->rad[0];
393 <                delta_t2 = d*d;
394 <                d = DOT(ck0, uvw[1]) / av->rad[1];
503 >                d = (u = DOT(ck0, uvw[0])) / av->rad[0];
504                  delta_t2 += d*d;
505 +                d = (v = DOT(ck0, uvw[1])) / av->rad[1];
506 +                delta_t2 += d*d;
507                  if (delta_t2 >= ambacc*ambacc)
508                          continue;
509                  /*
510 <                 *  Intersection behind test
510 >                 *  Test for potential light leak
511                   */
512 <                d = 0.0;
402 <                for (j = 0; j < 3; j++)
403 <                        d += (r->rop[j] - av->pos[j])*(uvw[2][j] + r->ron[j]);
404 <                if (d*0.5 < -minarad*ambacc-.001)
512 >                if (av->corral && plugaleak(r, av, uvw[2], atan2a(v,u)))
513                          continue;
514                  /*
515                   *  Extrapolate value and compute final weight (hat function)
516                   */
517 <                extambient(ct, av, r->rop, rn, uvw);
517 >                if (!extambient(ct, av, r->rop, rn, uvw))
518 >                        continue;
519                  d = tfunc(maxangle, sqrt(delta_r2), 0.0) *
520                          tfunc(ambacc, sqrt(delta_t2), 0.0);
521                  scalecolor(ct, d);
522                  addcolor(acol, ct);
523                  wsum += d;
524          }
416 –        if (at->kid == NULL)
417 –                return(wsum);
418 –                                        /* sum children */
419 –        s *= 0.5;
420 –        for (i = 0; i < 8; i++) {
421 –                for (j = 0; j < 3; j++) {
422 –                        ck0[j] = c0[j];
423 –                        if (1<<j & i)
424 –                                ck0[j] += s;
425 –                        if (r->rop[j] < ck0[j] - OCTSCALE*s)
426 –                                break;
427 –                        if (r->rop[j] > ck0[j] + (1.0+OCTSCALE)*s)
428 –                                break;
429 –                }
430 –                if (j == 3)
431 –                        wsum += sumambient(acol, r, rn, al,
432 –                                                at->kid+i, ck0, s);
433 –        }
525          return(wsum);
526   }
527  
528  
529 < int
529 > static int
530   makeambient(            /* make a new ambient value for storage */
531          COLOR  acol,
532          RAY  *r,
# Line 454 | Line 545 | makeambient(           /* make a new ambient value for storage
545                  amb.weight = 1.25*r->rweight;
546          setcolor(acol, AVGREFL, AVGREFL, AVGREFL);
547                                                  /* compute ambient */
548 <        if (!doambient(acol, r, amb.weight, uvw, amb.rad, amb.gpos, amb.gdir)) {
549 <                setcolor(acol, 0.0, 0.0, 0.0);
459 <                return(0);
460 <        }
548 >        i = doambient(acol, r, amb.weight,
549 >                        uvw, amb.rad, amb.gpos, amb.gdir, &amb.corral);
550          scalecolor(acol, 1./AVGREFL);           /* undo assumed reflectance */
551 +        if (i <= 0 || amb.rad[0] <= FTINY)      /* no Hessian or zero radius */
552 +                return(i);
553                                                  /* store value */
554          VCOPY(amb.pos, r->rop);
555          amb.ndir = encodedir(r->ron);
# Line 475 | Line 566 | makeambient(           /* make a new ambient value for storage
566   }
567  
568  
569 < void
569 > static int
570   extambient(             /* extrapolate value at pv, nv */
571          COLOR  cr,
572          AMBVAL   *ap,
# Line 484 | Line 575 | extambient(            /* extrapolate value at pv, nv */
575          FVECT  uvw[3]
576   )
577   {
578 +        const double    min_d = 0.05;
579          static FVECT    my_uvw[3];
580          FVECT           v1;
581          int             i;
# Line 503 | Line 595 | extambient(            /* extrapolate value at pv, nv */
595          for (i = 3; i--; )
596                  d += v1[i] * (ap->gdir[0]*uvw[0][i] + ap->gdir[1]*uvw[1][i]);
597          
598 <        if (d <= 0.0) {
599 <                setcolor(cr, 0.0, 0.0, 0.0);
508 <                return;
509 <        }
598 >        if (d < min_d)                  /* should not use if we can avoid it */
599 >                d = min_d;
600          copycolor(cr, ap->val);
601          scalecolor(cr, d);
602 +        return(d > min_d);
603   }
604  
605  
# Line 545 | Line 636 | avinsert(                              /* insert ambient value in our tree */
636          }
637          avh.next = at->alist;           /* order by increasing level */
638          for (ap = &avh; ap->next != NULL; ap = ap->next)
639 <                if (ap->next->lvl >= av->lvl)
639 >                if ( ap->next->lvl > av->lvl ||
640 >                                (ap->next->lvl == av->lvl) &
641 >                                (ap->next->weight <= av->weight) )
642                          break;
643          av->next = ap->next;
644          ap->next = (AMBVAL*)av;
# Line 569 | Line 662 | multambient(           /* compute ambient component & multiply
662   )
663   {
664          static int  rdepth = 0;                 /* ambient recursion */
665 <        COLOR   acol;
665 >        COLOR   acol, caustic;
666          double  d, l;
667  
668 +        /* PMAP: Factor in ambient from global photon map (if enabled) and return
669 +         * as all ambient components accounted for */
670 +        if (ambPmap(aval, r, rdepth))
671 +                return;
672 +
673 +        /* PMAP: Otherwise factor in ambient from caustic photon map
674 +         * (ambPmapCaustic() returns zero if caustic photons disabled) and
675 +         * continue with RADIANCE ambient calculation */
676 +        copycolor(caustic, aval);
677 +        ambPmapCaustic(caustic, r, rdepth);
678 +        
679          if (ambdiv <= 0)                        /* no ambient calculation */
680                  goto dumbamb;
681                                                  /* check number of bounces */
# Line 589 | Line 693 | multambient(           /* compute ambient component & multiply
693                  rdepth--;
694                  if (d <= FTINY)
695                          goto dumbamb;
696 <                copycolor(aval, acol);
696 >                copycolor(aval, acol);          
697 >        
698 >           /* PMAP: add in caustic */
699 >                addcolor(aval, caustic);        
700                  return;
701          }
702  
# Line 599 | Line 706 | multambient(           /* compute ambient component & multiply
706          setcolor(acol, 0.0, 0.0, 0.0);
707          d = sumambient(acol, r, nrm, rdepth,
708                          &atrunk, thescene.cuorg, thescene.cusize);
709 +                        
710          if (d > FTINY) {
711                  d = 1.0/d;
712                  scalecolor(acol, d);
713                  multcolor(aval, acol);
714 +                
715 +                /* PMAP: add in caustic */
716 +                addcolor(aval, caustic);        
717                  return;
718          }
719 +        
720          rdepth++;                               /* need to cache new value */
721          d = makeambient(acol, r, nrm, rdepth-1);
722          rdepth--;
723 +        
724          if (d > FTINY) {
725                  multcolor(aval, acol);          /* got new value */
726 +
727 +                /* PMAP: add in caustic */
728 +                addcolor(aval, caustic);                        
729                  return;
730          }
731 +        
732   dumbamb:                                        /* return global value */
733          if ((ambvwt <= 0) | (navsum == 0)) {
734                  multcolor(aval, ambval);
735 +
736 +                /* PMAP: add in caustic */
737 +                addcolor(aval, caustic);        
738                  return;
739          }
740 +        
741          l = bright(ambval);                     /* average in computations */
742          if (l > FTINY) {
743                  d = (log(l)*(double)ambvwt + avsum) /
# Line 658 | Line 779 | sumambient(    /* get interpolated ambient value */
779                  /*
780                   *  Ambient level test.
781                   */
782 <                if (av->lvl > al)       /* list sorted, so this works */
782 >                if (av->lvl > al ||     /* list sorted, so this works */
783 >                                (av->lvl == al) & (av->weight < 0.9*r->rweight))
784                          break;
663 –                if (av->weight < 0.9*r->rweight)
664 –                        continue;
785                  /*
786                   *  Ambient radius test.
787                   */
# Line 843 | Line 963 | avinsert(                              /* insert ambient value in our tree */
963          }
964          avh.next = at->alist;           /* order by increasing level */
965          for (ap = &avh; ap->next != NULL; ap = ap->next)
966 <                if (ap->next->lvl >= av->lvl)
966 >                if ( ap->next->lvl > av->lvl ||
967 >                                (ap->next->lvl == av->lvl) &
968 >                                (ap->next->weight <= av->weight) )
969                          break;
970          av->next = ap->next;
971          ap->next = (AMBVAL*)av;
# Line 895 | Line 1017 | avsave(                                /* insert and save an ambient value */
1017          AMBVAL  *av
1018   )
1019   {
1020 <        avinsert(avstore(av));
1020 >        avstore(av);
1021          if (ambfp == NULL)
1022                  return;
1023          if (writambval(av, ambfp) < 0)
# Line 910 | Line 1032 | writerr:
1032  
1033  
1034   static AMBVAL *
1035 < avstore(                                /* allocate memory and store aval */
1035 > avstore(                                /* allocate memory and save aval */
1036          AMBVAL  *aval
1037   )
1038   {
# Line 928 | Line 1050 | avstore(                               /* allocate memory and store aval */
1050                  avsum += log(d);
1051                  navsum++;
1052          }
1053 +        avinsert(av);                   /* insert in our cache tree */
1054          return(av);
1055   }
1056  
# Line 1128 | Line 1251 | sortambvals(                   /* resort ambient values */
1251                  if (i_avlist < nambvals)
1252                          error(CONSISTENCY, "missing ambient values in sortambvals");
1253   #endif
1254 <                qsort((char *)avlist1, nambvals, sizeof(struct avl), &alatcmp);
1255 <                qsort((char *)avlist2, nambvals, sizeof(AMBVAL *), &aposcmp);
1254 >                qsort((char *)avlist1, nambvals, sizeof(struct avl), alatcmp);
1255 >                qsort((char *)avlist2, nambvals, sizeof(AMBVAL *), aposcmp);
1256                  for (i = 0; i < nambvals; i++) {
1257                          if (avlist1[i].p == NULL)
1258                                  continue;
# Line 1210 | Line 1333 | ambsync(void)                  /* synchronize ambient file */
1333                                  error(WARNING, errmsg);
1334                                  break;
1335                          }
1336 <                        avinsert(avstore(&avs));
1336 >                        avstore(&avs);
1337                          n -= AMBVALSIZ;
1338                  }
1339                  lastpos = flen - n;

Diff Legend

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