--- ray/src/rt/ambient.c 2014/04/25 23:04:16 2.81 +++ ray/src/rt/ambient.c 2016/04/24 16:21:32 2.102 @@ -1,6 +1,4 @@ -#ifndef lint -static const char RCSid[] = "$Id: ambient.c,v 2.81 2014/04/25 23:04:16 greg Exp $"; -#endif +static const char RCSid[] = "$Id: ambient.c,v 2.102 2016/04/24 16:21:32 greg Exp $"; /* * ambient.c - routines dealing with ambient (inter-reflected) component. * @@ -17,6 +15,7 @@ static const char RCSid[] = "$Id: ambient.c,v 2.81 201 #include "resolu.h" #include "ambient.h" #include "random.h" +#include "pmapamb.h" #ifndef OCTSCALE #define OCTSCALE 1.0 /* ceil((valid rad.)/(cube size)) */ @@ -52,7 +51,6 @@ static int nunflshed = 0; /* number of unflushed ambi #endif -static double qambacc = 0.; /* ambient accuracy to the 1/4 power */ static double avsum = 0.; /* computed ambient value sum (log) */ static unsigned int navsum = 0; /* number of values in avsum */ static unsigned int nambvals = 0; /* total number of indirect values */ @@ -78,7 +76,6 @@ static long lastpos = -1; /* last flush position */ #define AMBFLUSH (BUFSIZ/AMBVALSIZ) #define newambval() (AMBVAL *)malloc(sizeof(AMBVAL)) -#define freeav(av) free((void *)av); static void initambfile(int creat); static void avsave(AMBVAL *av); @@ -110,17 +107,17 @@ setambres( /* set ambient resolution */ /* set min & max radii */ if (ar <= 0) { minarad = 0; - maxarad = thescene.cusize / 2.0; + maxarad = thescene.cusize*0.2; } else { minarad = thescene.cusize / ar; - maxarad = 64 * minarad; /* heuristic */ - if (maxarad > thescene.cusize / 2.0) - maxarad = thescene.cusize / 2.0; + maxarad = 64.0 * minarad; /* heuristic */ + if (maxarad > thescene.cusize*0.2) + maxarad = thescene.cusize*0.2; } if (minarad <= FTINY) - minarad = 10*FTINY; + minarad = 10.0*FTINY; if (maxarad <= minarad) - maxarad = 64 * minarad; + maxarad = 64.0 * minarad; } @@ -129,13 +126,11 @@ setambacc( /* set ambient accuracy */ double newa ) { - double ambdiff; - - if (newa < 0.0) - newa = 0.0; - ambdiff = fabs(newa - ambacc); - if (ambdiff >= .01 && (ambacc = newa) > FTINY) { - qambacc = sqrt(sqrt(ambacc)); + static double olda; /* remember previous setting here */ + + newa *= (newa > 0); + if (fabs(newa - olda) >= .05*(newa + olda)) { + ambacc = newa; if (nambvals > 0) sortambvals(1); /* rebuild tree */ } @@ -152,7 +147,7 @@ setambient(void) /* initialize calculation */ ambdone(); /* init ambient limits */ setambres(ambres); - qambacc = sqrt(sqrt(ambacc *= (ambacc > FTINY))); + setambacc(ambacc); if (ambfile == NULL || !ambfile[0]) return; if (ambacc <= FTINY) { @@ -168,7 +163,7 @@ setambient(void) /* initialize calculation */ initambfile(0); /* file exists */ lastpos = ftell(ambfp); while (readambval(&amb, ambfp)) - avinsert(avstore(&amb)); + avstore(&amb); nambshare = nambvals; /* share loaded values */ if (readonly) { sprintf(errmsg, @@ -188,9 +183,7 @@ setambient(void) /* initialize calculation */ (flen - lastpos)/AMBVALSIZ); error(WARNING, errmsg); fseek(ambfp, lastpos, SEEK_SET); -#ifndef _WIN32 /* XXX we need a replacement for that one */ ftruncate(fileno(ambfp), (off_t)lastpos); -#endif } } else if ((ambfp = fopen(ambfile, "w+")) != NULL) { initambfile(1); /* else create new file */ @@ -200,9 +193,6 @@ setambient(void) /* initialize calculation */ sprintf(errmsg, "cannot open ambient file \"%s\"", ambfile); error(SYSTEM, errmsg); } -#ifdef getc_unlocked - flockfile(ambfp); /* application-level lock */ -#endif #ifdef F_SETLKW aflock(F_UNLCK); /* release file */ #endif @@ -223,7 +213,7 @@ ambdone(void) /* close ambient file and free memory lastpos = -1; } /* free ambient tree */ - unloadatree(&atrunk, &avfree); + unloadatree(&atrunk, avfree); /* reset state variables */ avsum = 0.; navsum = 0; @@ -266,14 +256,15 @@ ambnotify( /* record new modifier */ /************ THE FOLLOWING ROUTINES DIFFER BETWEEN NEW & OLD ***************/ -#ifdef NEWAMB +#ifndef OLDAMB #define tfunc(lwr, x, upr) (((x)-(lwr))/((upr)-(lwr))) +static int plugaleak(RAY *r, AMBVAL *ap, FVECT anorm, double ang); static double sumambient(COLOR acol, RAY *r, FVECT rn, int al, AMBTREE *at, FVECT c0, double s); static int makeambient(COLOR acol, RAY *r, FVECT rn, int al); -static void extambient(COLOR cr, AMBVAL *ap, FVECT pv, FVECT nv, +static int extambient(COLOR cr, AMBVAL *ap, FVECT pv, FVECT nv, FVECT uvw[3]); void @@ -284,10 +275,22 @@ multambient( /* compute ambient component & multiply ) { static int rdepth = 0; /* ambient recursion */ - COLOR acol; + COLOR acol, caustic; int ok; double d, l; + /* PMAP: Factor in ambient from photon map, if enabled and ray is + * ambient. Return as all ambient components accounted for, else + * continue. */ + if (ambPmap(aval, r, rdepth)) + return; + + /* PMAP: Factor in specular-diffuse ambient (caustics) from photon + * map, if enabled and ray is primary, else caustic is zero. Continue + * with RADIANCE ambient calculation */ + copycolor(caustic, aval); + ambPmapCaustic(caustic, r, rdepth); + if (ambdiv <= 0) /* no ambient calculation */ goto dumbamb; /* check number of bounces */ @@ -301,11 +304,15 @@ multambient( /* compute ambient component & multiply if (ambacc <= FTINY) { /* no ambient storage */ copycolor(acol, aval); rdepth++; - ok = doambient(acol, r, r->rweight, NULL, NULL, NULL, NULL); + ok = doambient(acol, r, r->rweight, + NULL, NULL, NULL, NULL, NULL); rdepth--; if (!ok) goto dumbamb; copycolor(aval, acol); + + /* PMAP: add in caustic */ + addcolor(aval, caustic); return; } @@ -315,25 +322,39 @@ multambient( /* compute ambient component & multiply setcolor(acol, 0.0, 0.0, 0.0); d = sumambient(acol, r, nrm, rdepth, &atrunk, thescene.cuorg, thescene.cusize); + if (d > FTINY) { d = 1.0/d; scalecolor(acol, d); multcolor(aval, acol); + + /* PMAP: add in caustic */ + addcolor(aval, caustic); return; } + rdepth++; /* need to cache new value */ ok = makeambient(acol, r, nrm, rdepth-1); rdepth--; + if (ok) { multcolor(aval, acol); /* computed new value */ + + /* PMAP: add in caustic */ + addcolor(aval, caustic); return; } + dumbamb: /* return global value */ if ((ambvwt <= 0) | (navsum == 0)) { multcolor(aval, ambval); + + /* PMAP: add in caustic */ + addcolor(aval, caustic); return; } - l = bright(ambval); /* average in computations */ + + l = bright(ambval); /* average in computations */ if (l > FTINY) { d = (log(l)*(double)ambvwt + avsum) / (double)(ambvwt + navsum); @@ -347,7 +368,52 @@ dumbamb: /* return global value */ } -double +/* Plug a potential leak where ambient cache value is occluded */ +static int +plugaleak(RAY *r, AMBVAL *ap, FVECT anorm, double ang) +{ + const double cost70sq = 0.1169778; /* cos(70deg)^2 */ + RAY rtst; + FVECT vdif; + double normdot, ndotd, nadotd; + double a, b, c, t[2]; + + ang += 2.*PI*(ang < 0); /* check direction flags */ + if ( !(ap->corral>>(int)(ang*(16./PI)) & 1) ) + return(0); + /* + * Generate test ray, targeting 20 degrees above sample point plane + * along surface normal from cache position. This should be high + * enough to miss local geometry we don't really care about. + */ + VSUB(vdif, ap->pos, r->rop); + normdot = DOT(anorm, r->ron); + ndotd = DOT(vdif, r->ron); + nadotd = DOT(vdif, anorm); + a = normdot*normdot - cost70sq; + b = 2.0*(normdot*ndotd - nadotd*cost70sq); + c = ndotd*ndotd - DOT(vdif,vdif)*cost70sq; + if (quadratic(t, a, b, c) != 2) + return(1); /* should rarely happen */ + if (t[1] <= FTINY) + return(0); /* should fail behind test */ + rayorigin(&rtst, SHADOW, r, NULL); + VSUM(rtst.rdir, vdif, anorm, t[1]); /* further dist. > plane */ + rtst.rmax = normalize(rtst.rdir); /* short ray test */ + while (localhit(&rtst, &thescene)) { /* check for occluder */ + if (rtst.ro->omod != OVOID && + (rtst.clipset == NULL || + !inset(rtst.clipset, rtst.ro->omod))) + return(1); /* plug light leak */ + VCOPY(rtst.rorg, rtst.rop); /* skip invisible surface */ + rtst.rmax -= rtst.rot; + rayclear(&rtst); + } + return(0); /* seems we're OK */ +} + + +static double sumambient( /* get interpolated ambient value */ COLOR acol, RAY *r, @@ -359,15 +425,38 @@ sumambient( /* get interpolated ambient value */ ) { /* initial limit is 10 degrees plus ambacc radians */ const double minangle = 10.0 * PI/180.; - const double maxangle = (minangle+ambacc-PI/2.)*pow(r->rweight,0.13) - + PI/2.; + double maxangle = minangle + ambacc; double wsum = 0.0; FVECT ck0; int i, j; AMBVAL *av; + + if (at->kid != NULL) { /* sum children first */ + s *= 0.5; + for (i = 0; i < 8; i++) { + for (j = 0; j < 3; j++) { + ck0[j] = c0[j]; + if (1<rop[j] < ck0[j] - OCTSCALE*s) + break; + if (r->rop[j] > ck0[j] + (1.0+OCTSCALE)*s) + break; + } + if (j == 3) + wsum += sumambient(acol, r, rn, al, + at->kid+i, ck0, s); + } + /* good enough? */ + if (wsum >= 0.05 && s > minarad*10.0) + return(wsum); + } + /* adjust maximum angle */ + if (at->alist != NULL && (at->alist->lvl <= al) & (r->rweight < 0.6)) + maxangle = (maxangle - PI/2.)*pow(r->rweight,0.13) + PI/2.; /* sum this node */ for (av = at->alist; av != NULL; av = av->next) { - double d, delta_r2, delta_t2; + double u, v, d, delta_r2, delta_t2; COLOR ct; FVECT uvw[3]; /* record access */ @@ -376,10 +465,9 @@ sumambient( /* get interpolated ambient value */ /* * Ambient level test */ - if (av->lvl > al) /* list sorted, so this works */ + if (av->lvl > al || /* list sorted, so this works */ + (av->lvl == al) & (av->weight < 0.9*r->rweight)) break; - if (av->weight < 0.9*r->rweight) - continue; /* * Direction test using unperturbed normal */ @@ -393,58 +481,46 @@ sumambient( /* get interpolated ambient value */ /* * Modified ray behind test */ - VSUB(ck0, av->pos, r->rop); + VSUB(ck0, r->rop, av->pos); d = DOT(ck0, uvw[2]); - if (d < -minarad*qambacc-.001) + if (d < -minarad*ambacc-.001) continue; d /= av->rad[0]; delta_t2 = d*d; - if (delta_t2 >= qambacc*qambacc) + if (delta_t2 >= ambacc*ambacc) continue; /* * Elliptical radii test based on Hessian */ decodedir(uvw[0], av->udir); VCROSS(uvw[1], uvw[2], uvw[0]); - d = DOT(ck0, uvw[0]) / av->rad[0]; + d = (u = DOT(ck0, uvw[0])) / av->rad[0]; delta_t2 += d*d; - d = DOT(ck0, uvw[1]) / av->rad[1]; + d = (v = DOT(ck0, uvw[1])) / av->rad[1]; delta_t2 += d*d; - if (delta_t2 >= qambacc*qambacc) + if (delta_t2 >= ambacc*ambacc) continue; /* + * Test for potential light leak + */ + if (av->corral && plugaleak(r, av, uvw[2], atan2a(v,u))) + continue; + /* * Extrapolate value and compute final weight (hat function) */ - extambient(ct, av, r->rop, rn, uvw); + if (!extambient(ct, av, r->rop, rn, uvw)) + continue; d = tfunc(maxangle, sqrt(delta_r2), 0.0) * - tfunc(qambacc, sqrt(delta_t2), 0.0); + tfunc(ambacc, sqrt(delta_t2), 0.0); scalecolor(ct, d); addcolor(acol, ct); wsum += d; } - if (at->kid == NULL) - return(wsum); - /* sum children */ - s *= 0.5; - for (i = 0; i < 8; i++) { - for (j = 0; j < 3; j++) { - ck0[j] = c0[j]; - if (1<rop[j] < ck0[j] - OCTSCALE*s) - break; - if (r->rop[j] > ck0[j] + (1.0+OCTSCALE)*s) - break; - } - if (j == 3) - wsum += sumambient(acol, r, rn, al, - at->kid+i, ck0, s); - } return(wsum); } -int +static int makeambient( /* make a new ambient value for storage */ COLOR acol, RAY *r, @@ -463,7 +539,8 @@ makeambient( /* make a new ambient value for storage amb.weight = 1.25*r->rweight; setcolor(acol, AVGREFL, AVGREFL, AVGREFL); /* compute ambient */ - i = doambient(acol, r, amb.weight, uvw, amb.rad, amb.gpos, amb.gdir); + i = doambient(acol, r, amb.weight, + uvw, amb.rad, amb.gpos, amb.gdir, &amb.corral); scalecolor(acol, 1./AVGREFL); /* undo assumed reflectance */ if (i <= 0 || amb.rad[0] <= FTINY) /* no Hessian or zero radius */ return(i); @@ -483,7 +560,7 @@ makeambient( /* make a new ambient value for storage } -void +static int extambient( /* extrapolate value at pv, nv */ COLOR cr, AMBVAL *ap, @@ -492,6 +569,7 @@ extambient( /* extrapolate value at pv, nv */ FVECT uvw[3] ) { + const double min_d = 0.05; static FVECT my_uvw[3]; FVECT v1; int i; @@ -511,12 +589,11 @@ extambient( /* extrapolate value at pv, nv */ for (i = 3; i--; ) d += v1[i] * (ap->gdir[0]*uvw[0][i] + ap->gdir[1]*uvw[1][i]); - if (d <= 0.0) { - setcolor(cr, 0.0, 0.0, 0.0); - return; - } + if (d < min_d) /* should not use if we can avoid it */ + d = min_d; copycolor(cr, ap->val); scalecolor(cr, d); + return(d > min_d); } @@ -538,7 +615,7 @@ avinsert( /* insert ambient value in our tree */ at = &atrunk; VCOPY(ck0, thescene.cuorg); s = thescene.cusize; - while (s*(OCTSCALE/2) > av->rad[1]*qambacc) { + while (s*(OCTSCALE/2) > av->rad[1]*ambacc) { if (at->kid == NULL) if ((at->kid = newambtree()) == NULL) error(SYSTEM, "out of memory in avinsert"); @@ -553,7 +630,9 @@ avinsert( /* insert ambient value in our tree */ } avh.next = at->alist; /* order by increasing level */ for (ap = &avh; ap->next != NULL; ap = ap->next) - if (ap->next->lvl >= av->lvl) + if ( ap->next->lvl > av->lvl || + (ap->next->lvl == av->lvl) & + (ap->next->weight <= av->weight) ) break; av->next = ap->next; ap->next = (AMBVAL*)av; @@ -577,9 +656,20 @@ multambient( /* compute ambient component & multiply ) { static int rdepth = 0; /* ambient recursion */ - COLOR acol; + COLOR acol, caustic; double d, l; + /* PMAP: Factor in ambient from global photon map (if enabled) and return + * as all ambient components accounted for */ + if (ambPmap(aval, r, rdepth)) + return; + + /* PMAP: Otherwise factor in ambient from caustic photon map + * (ambPmapCaustic() returns zero if caustic photons disabled) and + * continue with RADIANCE ambient calculation */ + copycolor(caustic, aval); + ambPmapCaustic(caustic, r, rdepth); + if (ambdiv <= 0) /* no ambient calculation */ goto dumbamb; /* check number of bounces */ @@ -597,7 +687,10 @@ multambient( /* compute ambient component & multiply rdepth--; if (d <= FTINY) goto dumbamb; - copycolor(aval, acol); + copycolor(aval, acol); + + /* PMAP: add in caustic */ + addcolor(aval, caustic); return; } @@ -607,24 +700,38 @@ multambient( /* compute ambient component & multiply setcolor(acol, 0.0, 0.0, 0.0); d = sumambient(acol, r, nrm, rdepth, &atrunk, thescene.cuorg, thescene.cusize); + if (d > FTINY) { d = 1.0/d; scalecolor(acol, d); multcolor(aval, acol); + + /* PMAP: add in caustic */ + addcolor(aval, caustic); return; } + rdepth++; /* need to cache new value */ d = makeambient(acol, r, nrm, rdepth-1); rdepth--; + if (d > FTINY) { multcolor(aval, acol); /* got new value */ + + /* PMAP: add in caustic */ + addcolor(aval, caustic); return; } + dumbamb: /* return global value */ if ((ambvwt <= 0) | (navsum == 0)) { multcolor(aval, ambval); + + /* PMAP: add in caustic */ + addcolor(aval, caustic); return; } + l = bright(ambval); /* average in computations */ if (l > FTINY) { d = (log(l)*(double)ambvwt + avsum) / @@ -666,10 +773,9 @@ sumambient( /* get interpolated ambient value */ /* * Ambient level test. */ - if (av->lvl > al) /* list sorted, so this works */ + if (av->lvl > al || /* list sorted, so this works */ + (av->lvl == al) & (av->weight < 0.9*r->rweight)) break; - if (av->weight < 0.9*r->rweight) - continue; /* * Ambient radius test. */ @@ -851,7 +957,9 @@ avinsert( /* insert ambient value in our tree */ } avh.next = at->alist; /* order by increasing level */ for (ap = &avh; ap->next != NULL; ap = ap->next) - if (ap->next->lvl >= av->lvl) + if ( ap->next->lvl > av->lvl || + (ap->next->lvl == av->lvl) & + (ap->next->weight <= av->weight) ) break; av->next = ap->next; ap->next = (AMBVAL*)av; @@ -903,7 +1011,7 @@ avsave( /* insert and save an ambient value */ AMBVAL *av ) { - avinsert(avstore(av)); + avstore(av); if (ambfp == NULL) return; if (writambval(av, ambfp) < 0) @@ -918,7 +1026,7 @@ writerr: static AMBVAL * -avstore( /* allocate memory and store aval */ +avstore( /* allocate memory and save aval */ AMBVAL *aval ) { @@ -936,6 +1044,7 @@ avstore( /* allocate memory and store aval */ avsum += log(d); navsum++; } + avinsert(av); /* insert in our cache tree */ return(av); } @@ -962,7 +1071,7 @@ newambtree(void) /* allocate 8 ambient tree structs } atp = atfreelist; atfreelist = atp->kid; - memset((char *)atp, '\0', 8*sizeof(AMBTREE)); + memset(atp, 0, 8*sizeof(AMBTREE)); return(atp); } @@ -988,6 +1097,7 @@ unloadatree( /* unload an ambient value tree */ /* transfer values at this node */ for (av = at->alist; av != NULL; av = at->alist) { at->alist = av->next; + av->next = NULL; (*f)(av); } if (at->kid == NULL) @@ -1065,8 +1175,8 @@ avlmemi( /* find list position from address */ { AMBVAL **avlpp; - avlpp = (AMBVAL **)bsearch((char *)&avaddr, (char *)avlist2, - nambvals, sizeof(AMBVAL *), &aposcmp); + avlpp = (AMBVAL **)bsearch(&avaddr, avlist2, + nambvals, sizeof(AMBVAL *), aposcmp); if (avlpp == NULL) error(CONSISTENCY, "address not found in avlmemi"); return(avlpp - avlist2); @@ -1109,12 +1219,12 @@ sortambvals( /* resort ambient values */ } if (avlist1 == NULL) { /* no time tracking -- rebuild tree? */ if (avlist2 != NULL) - free((void *)avlist2); + free(avlist2); if (always) { /* rebuild without sorting */ oldatrunk = atrunk; atrunk.alist = NULL; atrunk.kid = NULL; - unloadatree(&oldatrunk, &avinsert); + unloadatree(&oldatrunk, avinsert); } } else { /* sort memory by last access time */ /* @@ -1131,13 +1241,13 @@ sortambvals( /* resort ambient values */ eputs(errmsg); #endif i_avlist = 0; - unloadatree(&atrunk, &av2list); /* empty current tree */ + unloadatree(&atrunk, av2list); /* empty current tree */ #ifdef DEBUG if (i_avlist < nambvals) error(CONSISTENCY, "missing ambient values in sortambvals"); #endif - qsort((char *)avlist1, nambvals, sizeof(struct avl), &alatcmp); - qsort((char *)avlist2, nambvals, sizeof(AMBVAL *), &aposcmp); + qsort(avlist1, nambvals, sizeof(struct avl), alatcmp); + qsort(avlist2, nambvals, sizeof(AMBVAL *), aposcmp); for (i = 0; i < nambvals; i++) { if (avlist1[i].p == NULL) continue; @@ -1153,8 +1263,8 @@ sortambvals( /* resort ambient values */ avinsert(avlist2[j]); avlist1[j].p = NULL; } - free((void *)avlist1); - free((void *)avlist2); + free(avlist1); + free(avlist2); /* compute new sort interval */ sortintvl = ambclock - lastsort; if (sortintvl >= MAX_SORT_INTVL/2) @@ -1218,7 +1328,7 @@ ambsync(void) /* synchronize ambient file */ error(WARNING, errmsg); break; } - avinsert(avstore(&avs)); + avstore(&avs); n -= AMBVALSIZ; } lastpos = flen - n;