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.71 by greg, Fri Apr 11 20:31:37 2014 UTC vs.
Revision 2.92 by greg, Fri Nov 21 00:30:11 2014 UTC

# Line 51 | Line 51 | static int  nunflshed = 0;     /* number of unflushed ambi
51   #define MAX_SORT_INTVL  (SORT_INTVL<<6)
52   #endif
53  
54 +
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 108 | 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 127 | 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);
136 <        if (ambdiff >= .01 && (ambacc = newa) > FTINY && nambvals > 0)
137 <                sortambvals(1);                 /* rebuild tree */
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 >        }
139   }
140  
141  
# Line 163 | 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 261 | 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)))
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 int      extambient(COLOR cr, AMBVAL *ap, FVECT pv, FVECT nv,
275 +                                FVECT uvw[3]);
276  
277   void
278   multambient(            /* compute ambient component & multiply by coef. */
# Line 291 | 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 315 | Line 324 | multambient(           /* compute ambient component & multiply
324          ok = makeambient(acol, r, nrm, rdepth-1);
325          rdepth--;
326          if (ok) {
327 <                multcolor(aval, acol);          /* got new value */
327 >                multcolor(aval, acol);          /* computed new value */
328                  return;
329          }
330   dumbamb:                                        /* return global value */
# Line 337 | Line 346 | dumbamb:                                       /* return global value */
346   }
347  
348  
349 < double
350 < sumambient(     /* get interpolated ambient value */
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,
398          FVECT  rn,
# Line 347 | Line 401 | sumambient(    /* get interpolated ambient value */
401          FVECT  c0,
402          double  s
403   )
404 < {                                       /* initial limit is ambacc radians */
405 <        const double    maxangle = (ambacc-PI/2.)*pow(r->rweight,0.13) + PI/2.;
404 > {                       /* initial limit is 10 degrees plus ambacc radians */
405 >        const double    minangle = 10.0 * PI/180.;
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 */
441                  if (tracktime)
442                          av->latick = ambclock;
443                  /*
444 <                 *  Ambient level test.
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;
369 –                if (av->weight < 0.9*r->rweight)
370 –                        continue;
449                  /*
450 <                 *  Direction test using unperturbed normal.
450 >                 *  Direction test using unperturbed normal
451                   */
452                  decodedir(uvw[2], av->ndir);
453                  d = DOT(uvw[2], r->ron);
# Line 379 | Line 457 | sumambient(    /* get interpolated ambient value */
457                  if (delta_r2 >= maxangle*maxangle)
458                          continue;
459                  /*
460 <                 *  Ambient radius test.
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);
387 <                d = DOT(ck0, uvw[0]) / av->rad[0];
388 <                delta_t2 = d*d;
389 <                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 = (v = DOT(ck0, uvw[1])) / av->rad[1];
478 +                delta_t2 += d*d;
479                  if (delta_t2 >= ambacc*ambacc)
480                          continue;
481                  /*
482 <                 *  Ray behind test.
482 >                 *  Test for potential light leak
483                   */
484 <                d = 0.0;
397 <                for (j = 0; j < 3; j++)
398 <                        d += (r->rop[j] - av->pos[j])*(uvw[2][j] + r->ron[j]);
399 <                if (d*0.5 < -minarad*ambacc-.001)
484 >                if (av->corral && plugaleak(r, av, uvw[2], atan2a(v,u)))
485                          continue;
486                  /*
487 <                 *  Convert to final weight (hat function)
487 >                 *  Extrapolate value and compute final weight (hat function)
488                   */
489 +                if (!extambient(ct, av, r->rop, rn, uvw))
490 +                        continue;
491                  d = tfunc(maxangle, sqrt(delta_r2), 0.0) *
492                          tfunc(ambacc, sqrt(delta_t2), 0.0);
406 –                wsum += d;
407 –                extambient(ct, av, uvw, r->rop, rn);
493                  scalecolor(ct, d);
494                  addcolor(acol, ct);
495 +                wsum += d;
496          }
411 –        if (at->kid == NULL)
412 –                return(wsum);
413 –                                        /* sum children */
414 –        s *= 0.5;
415 –        for (i = 0; i < 8; i++) {
416 –                for (j = 0; j < 3; j++) {
417 –                        ck0[j] = c0[j];
418 –                        if (1<<j & i)
419 –                                ck0[j] += s;
420 –                        if (r->rop[j] < ck0[j] - OCTSCALE*s)
421 –                                break;
422 –                        if (r->rop[j] > ck0[j] + (1.0+OCTSCALE)*s)
423 –                                break;
424 –                }
425 –                if (j == 3)
426 –                        wsum += sumambient(acol, r, rn, al,
427 –                                                at->kid+i, ck0, s);
428 –        }
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 439 | Line 507 | makeambient(           /* make a new ambient value for storage
507   )
508   {
509          AMBVAL  amb;
510 <        FVECT   uv[2];
510 >        FVECT   uvw[3];
511          int     i;
512  
513          amb.weight = 1.0;                       /* compute weight */
# Line 449 | 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 <        if (!doambient(acol, r, amb.weight, uv, amb.rad, amb.gpos, amb.gdir)) {
521 <                setcolor(acol, 0.0, 0.0, 0.0);
454 <                return(0);
455 <        }
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);
525                                                  /* store value */
526          VCOPY(amb.pos, r->rop);
527          amb.ndir = encodedir(r->ron);
528 <        amb.udir = encodedir(uv[0]);
528 >        amb.udir = encodedir(uvw[0]);
529          amb.lvl = al;
530          copycolor(amb.val, acol);
531                                                  /* insert into tree */
532          avsave(&amb);                           /* and save to file */
533 <        if (rn != r->ron)
534 <                extambient(acol, &amb, r->rop, rn);     /* texture */
533 >        if (rn != r->ron) {                     /* texture */
534 >                VCOPY(uvw[2], r->ron);
535 >                extambient(acol, &amb, r->rop, rn, uvw);
536 >        }
537          return(1);
538   }
539  
540  
541 < void
541 > static int
542   extambient(             /* extrapolate value at pv, nv */
543          COLOR  cr,
544          AMBVAL   *ap,
475 –        FVECT  uvw[3],
545          FVECT  pv,
546 <        FVECT  nv
546 >        FVECT  nv,
547 >        FVECT  uvw[3]
548   )
549   {
550 <        FVECT  v1;
551 <        int  i;
552 <        double  d = 1.0;                /* zeroeth order */
550 >        static FVECT    my_uvw[3];
551 >        FVECT           v1;
552 >        int             i;
553 >        double          d = 1.0;        /* zeroeth order */
554  
555 +        if (uvw == NULL) {              /* need local coordinates? */
556 +                decodedir(my_uvw[2], ap->ndir);
557 +                decodedir(my_uvw[0], ap->udir);
558 +                VCROSS(my_uvw[1], my_uvw[2], my_uvw[0]);
559 +                uvw = my_uvw;
560 +        }
561          for (i = 3; i--; )              /* gradient due to translation */
562                  d += (pv[i] - ap->pos[i]) *
563                          (ap->gpos[0]*uvw[0][i] + ap->gpos[1]*uvw[1][i]);
# Line 491 | Line 568 | extambient(            /* extrapolate value at pv, nv */
568          
569          if (d <= 0.0) {
570                  setcolor(cr, 0.0, 0.0, 0.0);
571 <                return;
571 >                return(0);              /* should not use if we can avoid it */
572          }
573          copycolor(cr, ap->val);
574          scalecolor(cr, d);
575 +        return(1);
576   }
577  
578  
# Line 531 | Line 609 | avinsert(                              /* insert ambient value in our tree */
609          }
610          avh.next = at->alist;           /* order by increasing level */
611          for (ap = &avh; ap->next != NULL; ap = ap->next)
612 <                if (ap->next->lvl >= av->lvl)
612 >                if ( ap->next->lvl > av->lvl ||
613 >                                (ap->next->lvl == av->lvl) &
614 >                                (ap->next->weight <= av->weight) )
615                          break;
616          av->next = ap->next;
617          ap->next = (AMBVAL*)av;
# Line 541 | Line 621 | avinsert(                              /* insert ambient value in our tree */
621  
622   #else /* ! NEWAMB */
623  
624 + static double   sumambient(COLOR acol, RAY *r, FVECT rn, int al,
625 +                                AMBTREE *at, FVECT c0, double s);
626 + static double   makeambient(COLOR acol, RAY *r, FVECT rn, int al);
627 + static void     extambient(COLOR cr, AMBVAL *ap, FVECT pv, FVECT nv);
628  
629 +
630   void
631   multambient(            /* compute ambient component & multiply by coef. */
632          COLOR  aval,
# Line 612 | Line 697 | dumbamb:                                       /* return global value */
697   }
698  
699  
700 < double
700 > static double
701   sumambient(     /* get interpolated ambient value */
702          COLOR  acol,
703          RAY  *r,
# Line 639 | Line 724 | sumambient(    /* get interpolated ambient value */
724                  /*
725                   *  Ambient level test.
726                   */
727 <                if (av->lvl > al)       /* list sorted, so this works */
727 >                if (av->lvl > al ||     /* list sorted, so this works */
728 >                                (av->lvl == al) & (av->weight < 0.9*r->rweight))
729                          break;
644 –                if (av->weight < 0.9*r->rweight)
645 –                        continue;
730                  /*
731                   *  Ambient radius test.
732                   */
# Line 723 | Line 807 | sumambient(    /* get interpolated ambient value */
807   }
808  
809  
810 < double
810 > static double
811   makeambient(            /* make a new ambient value for storage */
812          COLOR  acol,
813          RAY  *r,
# Line 763 | Line 847 | makeambient(           /* make a new ambient value for storage
847   }
848  
849  
850 < void
850 > static void
851   extambient(             /* extrapolate value at pv, nv */
852          COLOR  cr,
853          AMBVAL   *ap,
# Line 824 | Line 908 | avinsert(                              /* insert ambient value in our tree */
908          }
909          avh.next = at->alist;           /* order by increasing level */
910          for (ap = &avh; ap->next != NULL; ap = ap->next)
911 <                if (ap->next->lvl >= av->lvl)
911 >                if ( ap->next->lvl > av->lvl ||
912 >                                (ap->next->lvl == av->lvl) &
913 >                                (ap->next->weight <= av->weight) )
914                          break;
915          av->next = ap->next;
916          ap->next = (AMBVAL*)av;
# Line 876 | Line 962 | avsave(                                /* insert and save an ambient value */
962          AMBVAL  *av
963   )
964   {
965 <        avinsert(avstore(av));
965 >        avstore(av);
966          if (ambfp == NULL)
967                  return;
968          if (writambval(av, ambfp) < 0)
# Line 891 | Line 977 | writerr:
977  
978  
979   static AMBVAL *
980 < avstore(                                /* allocate memory and store aval */
980 > avstore(                                /* allocate memory and save aval */
981          AMBVAL  *aval
982   )
983   {
# Line 909 | Line 995 | avstore(                               /* allocate memory and store aval */
995                  avsum += log(d);
996                  navsum++;
997          }
998 +        avinsert(av);                   /* insert in our cache tree */
999          return(av);
1000   }
1001  
# Line 1109 | Line 1196 | sortambvals(                   /* resort ambient values */
1196                  if (i_avlist < nambvals)
1197                          error(CONSISTENCY, "missing ambient values in sortambvals");
1198   #endif
1199 <                qsort((char *)avlist1, nambvals, sizeof(struct avl), &alatcmp);
1200 <                qsort((char *)avlist2, nambvals, sizeof(AMBVAL *), &aposcmp);
1199 >                qsort((char *)avlist1, nambvals, sizeof(struct avl), alatcmp);
1200 >                qsort((char *)avlist2, nambvals, sizeof(AMBVAL *), aposcmp);
1201                  for (i = 0; i < nambvals; i++) {
1202                          if (avlist1[i].p == NULL)
1203                                  continue;
# Line 1191 | Line 1278 | ambsync(void)                  /* synchronize ambient file */
1278                                  error(WARNING, errmsg);
1279                                  break;
1280                          }
1281 <                        avinsert(avstore(&avs));
1281 >                        avstore(&avs);
1282                          n -= AMBVALSIZ;
1283                  }
1284                  lastpos = flen - n;

Diff Legend

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