ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/ambient.c
Revision: 2.105
Committed: Tue Jan 9 05:01:15 2018 UTC (6 years, 4 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R2
Changes since 2.104: +5 -2 lines
Log Message:
Made proximity test (plugaleak) more robust with transparent surfaces, etc.

File Contents

# Content
1 static const char RCSid[] = "$Id: ambient.c,v 2.104 2017/01/27 22:00:49 greg Exp $";
2 /*
3 * ambient.c - routines dealing with ambient (inter-reflected) component.
4 *
5 * Declarations of external symbols in ambient.h
6 */
7
8 #include "copyright.h"
9
10 #include <string.h>
11
12 #include "platform.h"
13 #include "ray.h"
14 #include "otypes.h"
15 #include "otspecial.h"
16 #include "resolu.h"
17 #include "ambient.h"
18 #include "source.h"
19 #include "random.h"
20 #include "pmapamb.h"
21
22 #ifndef OCTSCALE
23 #define OCTSCALE 1.0 /* ceil((valid rad.)/(cube size)) */
24 #endif
25
26 extern char *shm_boundary; /* memory sharing boundary */
27
28 #ifndef MAXASET
29 #define MAXASET 4095 /* maximum number of elements in ambient set */
30 #endif
31 OBJECT ambset[MAXASET+1]={0}; /* ambient include/exclude set */
32
33 double maxarad; /* maximum ambient radius */
34 double minarad; /* minimum ambient radius */
35
36 static AMBTREE atrunk; /* our ambient trunk node */
37
38 static FILE *ambfp = NULL; /* ambient file pointer */
39 static int nunflshed = 0; /* number of unflushed ambient values */
40
41 #ifndef SORT_THRESH
42 #ifdef SMLMEM
43 #define SORT_THRESH ((16L<<20)/sizeof(AMBVAL))
44 #else
45 #define SORT_THRESH ((64L<<20)/sizeof(AMBVAL))
46 #endif
47 #endif
48 #ifndef SORT_INTVL
49 #define SORT_INTVL (SORT_THRESH<<1)
50 #endif
51 #ifndef MAX_SORT_INTVL
52 #define MAX_SORT_INTVL (SORT_INTVL<<6)
53 #endif
54
55
56 static double avsum = 0.; /* computed ambient value sum (log) */
57 static unsigned int navsum = 0; /* number of values in avsum */
58 static unsigned int nambvals = 0; /* total number of indirect values */
59 static unsigned int nambshare = 0; /* number of values from file */
60 static unsigned long ambclock = 0; /* ambient access clock */
61 static unsigned long lastsort = 0; /* time of last value sort */
62 static long sortintvl = SORT_INTVL; /* time until next sort */
63 static FILE *ambinp = NULL; /* auxiliary file for input */
64 static long lastpos = -1; /* last flush position */
65
66 #define MAXACLOCK (1L<<30) /* clock turnover value */
67 /*
68 * Track access times unless we are sharing ambient values
69 * through memory on a multiprocessor, when we want to avoid
70 * claiming our own memory (copy on write). Go ahead anyway
71 * if more than two thirds of our values are unshared.
72 * Compile with -Dtracktime=0 to turn this code off.
73 */
74 #ifndef tracktime
75 #define tracktime (shm_boundary == NULL || nambvals > 3*nambshare)
76 #endif
77
78 #define AMBFLUSH (BUFSIZ/AMBVALSIZ)
79
80 #define newambval() (AMBVAL *)malloc(sizeof(AMBVAL))
81
82 static void initambfile(int creat);
83 static void avsave(AMBVAL *av);
84 static AMBVAL *avstore(AMBVAL *aval);
85 static AMBTREE *newambtree(void);
86 static void freeambtree(AMBTREE *atp);
87
88 typedef void unloadtf_t(AMBVAL *);
89 static unloadtf_t avinsert;
90 static unloadtf_t av2list;
91 static unloadtf_t avfree;
92 static void unloadatree(AMBTREE *at, unloadtf_t *f);
93
94 static int aposcmp(const void *avp1, const void *avp2);
95 static int avlmemi(AMBVAL *avaddr);
96 static void sortambvals(int always);
97
98 #ifdef F_SETLKW
99 static void aflock(int typ);
100 #endif
101
102
103 void
104 setambres( /* set ambient resolution */
105 int ar
106 )
107 {
108 ambres = ar < 0 ? 0 : ar; /* may be done already */
109 /* set min & max radii */
110 if (ar <= 0) {
111 minarad = 0;
112 maxarad = thescene.cusize*0.2;
113 } else {
114 minarad = thescene.cusize / ar;
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.0*FTINY;
121 if (maxarad <= minarad)
122 maxarad = 64.0 * minarad;
123 }
124
125
126 void
127 setambacc( /* set ambient accuracy */
128 double newa
129 )
130 {
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
142 void
143 setambient(void) /* initialize calculation */
144 {
145 int readonly = 0;
146 long flen;
147 AMBVAL amb;
148 /* make sure we're fresh */
149 ambdone();
150 /* init ambient limits */
151 setambres(ambres);
152 setambacc(ambacc);
153 if (ambfile == NULL || !ambfile[0])
154 return;
155 if (ambacc <= FTINY) {
156 sprintf(errmsg, "zero ambient accuracy so \"%s\" not opened",
157 ambfile);
158 error(WARNING, errmsg);
159 return;
160 }
161 /* open ambient file */
162 if ((ambfp = fopen(ambfile, "r+")) == NULL)
163 readonly = (ambfp = fopen(ambfile, "r")) != NULL;
164 if (ambfp != NULL) {
165 initambfile(0); /* file exists */
166 lastpos = ftell(ambfp);
167 while (readambval(&amb, ambfp))
168 avstore(&amb);
169 nambshare = nambvals; /* share loaded values */
170 if (readonly) {
171 sprintf(errmsg,
172 "loaded %u values from read-only ambient file",
173 nambvals);
174 error(WARNING, errmsg);
175 fclose(ambfp); /* close file so no writes */
176 ambfp = NULL;
177 return; /* avoid ambsync() */
178 }
179 /* align file pointer */
180 lastpos += (long)nambvals*AMBVALSIZ;
181 flen = lseek(fileno(ambfp), (off_t)0, SEEK_END);
182 if (flen != lastpos) {
183 sprintf(errmsg,
184 "ignoring last %ld values in ambient file (corrupted)",
185 (flen - lastpos)/AMBVALSIZ);
186 error(WARNING, errmsg);
187 fseek(ambfp, lastpos, SEEK_SET);
188 ftruncate(fileno(ambfp), (off_t)lastpos);
189 }
190 } else if ((ambfp = fopen(ambfile, "w+")) != NULL) {
191 initambfile(1); /* else create new file */
192 fflush(ambfp);
193 lastpos = ftell(ambfp);
194 } else {
195 sprintf(errmsg, "cannot open ambient file \"%s\"", ambfile);
196 error(SYSTEM, errmsg);
197 }
198 #ifdef F_SETLKW
199 aflock(F_UNLCK); /* release file */
200 #endif
201 }
202
203
204 void
205 ambdone(void) /* close ambient file and free memory */
206 {
207 if (ambfp != NULL) { /* close ambient file */
208 ambsync();
209 fclose(ambfp);
210 ambfp = NULL;
211 if (ambinp != NULL) {
212 fclose(ambinp);
213 ambinp = NULL;
214 }
215 lastpos = -1;
216 }
217 /* free ambient tree */
218 unloadatree(&atrunk, avfree);
219 /* reset state variables */
220 avsum = 0.;
221 navsum = 0;
222 nambvals = 0;
223 nambshare = 0;
224 ambclock = 0;
225 lastsort = 0;
226 sortintvl = SORT_INTVL;
227 }
228
229
230 void
231 ambnotify( /* record new modifier */
232 OBJECT obj
233 )
234 {
235 static int hitlimit = 0;
236 OBJREC *o;
237 char **amblp;
238
239 if (obj == OVOID) { /* starting over */
240 ambset[0] = 0;
241 hitlimit = 0;
242 return;
243 }
244 o = objptr(obj);
245 if (hitlimit || !ismodifier(o->otype))
246 return;
247 for (amblp = amblist; *amblp != NULL; amblp++)
248 if (!strcmp(o->oname, *amblp)) {
249 if (ambset[0] >= MAXASET) {
250 error(WARNING, "too many modifiers in ambient list");
251 hitlimit++;
252 return; /* should this be fatal? */
253 }
254 insertelem(ambset, obj);
255 return;
256 }
257 }
258
259 /************ THE FOLLOWING ROUTINES DIFFER BETWEEN NEW & OLD ***************/
260
261 #ifndef OLDAMB
262
263 #define tfunc(lwr, x, upr) (((x)-(lwr))/((upr)-(lwr)))
264
265 static int plugaleak(RAY *r, AMBVAL *ap, FVECT anorm, double ang);
266 static double sumambient(COLOR acol, RAY *r, FVECT rn, int al,
267 AMBTREE *at, FVECT c0, double s);
268 static int makeambient(COLOR acol, RAY *r, FVECT rn, int al);
269 static int extambient(COLOR cr, AMBVAL *ap, FVECT pv, FVECT nv,
270 FVECT uvw[3]);
271
272 void
273 multambient( /* compute ambient component & multiply by coef. */
274 COLOR aval,
275 RAY *r,
276 FVECT nrm
277 )
278 {
279 static int rdepth = 0; /* ambient recursion */
280 COLOR acol, caustic;
281 int i, ok;
282 double d, l;
283
284 /* PMAP: Factor in ambient from photon map, if enabled and ray is
285 * ambient. Return as all ambient components accounted for, else
286 * continue. */
287 if (ambPmap(aval, r, rdepth))
288 return;
289
290 /* PMAP: Factor in specular-diffuse ambient (caustics) from photon
291 * map, if enabled and ray is primary, else caustic is zero. Continue
292 * with RADIANCE ambient calculation */
293 copycolor(caustic, aval);
294 ambPmapCaustic(caustic, r, rdepth);
295
296 if (ambdiv <= 0) /* no ambient calculation */
297 goto dumbamb;
298 /* check number of bounces */
299 if (rdepth >= ambounce)
300 goto dumbamb;
301 /* check ambient list */
302 if (ambincl != -1 && r->ro != NULL &&
303 ambincl != inset(ambset, r->ro->omod))
304 goto dumbamb;
305
306 if (ambacc <= FTINY) { /* no ambient storage */
307 FVECT uvd[2];
308 float dgrad[2], *dgp = NULL;
309
310 if (nrm != r->ron && DOT(nrm,r->ron) < 0.9999)
311 dgp = dgrad; /* compute rotational grad. */
312 copycolor(acol, aval);
313 rdepth++;
314 ok = doambient(acol, r, r->rweight,
315 uvd, NULL, NULL, dgp, NULL);
316 rdepth--;
317 if (!ok)
318 goto dumbamb;
319 if ((ok > 0) & (dgp != NULL)) { /* apply texture */
320 FVECT v1;
321 VCROSS(v1, r->ron, nrm);
322 d = 1.0;
323 for (i = 3; i--; )
324 d += v1[i] * (dgp[0]*uvd[0][i] + dgp[1]*uvd[1][i]);
325 if (d >= 0.05)
326 scalecolor(acol, d);
327 }
328 copycolor(aval, acol);
329
330 /* PMAP: add in caustic */
331 addcolor(aval, caustic);
332 return;
333 }
334
335 if (tracktime) /* sort to minimize thrashing */
336 sortambvals(0);
337 /* interpolate ambient value */
338 setcolor(acol, 0.0, 0.0, 0.0);
339 d = sumambient(acol, r, nrm, rdepth,
340 &atrunk, thescene.cuorg, thescene.cusize);
341
342 if (d > FTINY) {
343 d = 1.0/d;
344 scalecolor(acol, d);
345 multcolor(aval, acol);
346
347 /* PMAP: add in caustic */
348 addcolor(aval, caustic);
349 return;
350 }
351
352 rdepth++; /* need to cache new value */
353 ok = makeambient(acol, r, nrm, rdepth-1);
354 rdepth--;
355
356 if (ok) {
357 multcolor(aval, acol); /* computed new value */
358
359 /* PMAP: add in caustic */
360 addcolor(aval, caustic);
361 return;
362 }
363
364 dumbamb: /* return global value */
365 if ((ambvwt <= 0) | (navsum == 0)) {
366 multcolor(aval, ambval);
367
368 /* PMAP: add in caustic */
369 addcolor(aval, caustic);
370 return;
371 }
372
373 l = bright(ambval); /* average in computations */
374 if (l > FTINY) {
375 d = (log(l)*(double)ambvwt + avsum) /
376 (double)(ambvwt + navsum);
377 d = exp(d) / l;
378 scalecolor(aval, d);
379 multcolor(aval, ambval); /* apply color of ambval */
380 } else {
381 d = exp( avsum / (double)navsum );
382 scalecolor(aval, d); /* neutral color */
383 }
384 }
385
386
387 /* Plug a potential leak where ambient cache value is occluded */
388 static int
389 plugaleak(RAY *r, AMBVAL *ap, FVECT anorm, double ang)
390 {
391 const double cost70sq = 0.1169778; /* cos(70deg)^2 */
392 RAY rtst;
393 FVECT vdif;
394 double normdot, ndotd, nadotd;
395 double a, b, c, t[2];
396
397 ang += 2.*PI*(ang < 0); /* check direction flags */
398 if ( !(ap->corral>>(int)(ang*(16./PI)) & 1) )
399 return(0);
400 /*
401 * Generate test ray, targeting 20 degrees above sample point plane
402 * along surface normal from cache position. This should be high
403 * enough to miss local geometry we don't really care about.
404 */
405 VSUB(vdif, ap->pos, r->rop);
406 normdot = DOT(anorm, r->ron);
407 ndotd = DOT(vdif, r->ron);
408 nadotd = DOT(vdif, anorm);
409 a = normdot*normdot - cost70sq;
410 b = 2.0*(normdot*ndotd - nadotd*cost70sq);
411 c = ndotd*ndotd - DOT(vdif,vdif)*cost70sq;
412 if (quadratic(t, a, b, c) != 2)
413 return(1); /* should rarely happen */
414 if (t[1] <= FTINY)
415 return(0); /* should fail behind test */
416 rayorigin(&rtst, SHADOW, r, NULL);
417 VSUM(rtst.rdir, vdif, anorm, t[1]); /* further dist. > plane */
418 rtst.rmax = normalize(rtst.rdir); /* short ray test */
419 while (localhit(&rtst, &thescene)) { /* check for occluder */
420 OBJREC *m = findmaterial(rtst.ro);
421 if (m != NULL && !istransp(m->otype) && !isBSDFproxy(m) &&
422 (rtst.clipset == NULL ||
423 !inset(rtst.clipset, rtst.ro->omod)))
424 return(1); /* plug light leak */
425 VCOPY(rtst.rorg, rtst.rop); /* skip invisible surface */
426 rtst.rmax -= rtst.rot;
427 rayclear(&rtst);
428 }
429 return(0); /* seems we're OK */
430 }
431
432
433 static double
434 sumambient( /* get interpolated ambient value */
435 COLOR acol,
436 RAY *r,
437 FVECT rn,
438 int al,
439 AMBTREE *at,
440 FVECT c0,
441 double s
442 )
443 { /* initial limit is 10 degrees plus ambacc radians */
444 const double minangle = 10.0 * PI/180.;
445 double maxangle = minangle + ambacc;
446 double wsum = 0.0;
447 FVECT ck0;
448 int i, j;
449 AMBVAL *av;
450
451 if (at->kid != NULL) { /* sum children first */
452 s *= 0.5;
453 for (i = 0; i < 8; i++) {
454 for (j = 0; j < 3; j++) {
455 ck0[j] = c0[j];
456 if (1<<j & i)
457 ck0[j] += s;
458 if (r->rop[j] < ck0[j] - OCTSCALE*s)
459 break;
460 if (r->rop[j] > ck0[j] + (1.0+OCTSCALE)*s)
461 break;
462 }
463 if (j == 3)
464 wsum += sumambient(acol, r, rn, al,
465 at->kid+i, ck0, s);
466 }
467 /* good enough? */
468 if (wsum >= 0.05 && s > minarad*10.0)
469 return(wsum);
470 }
471 /* adjust maximum angle */
472 if (at->alist != NULL && (at->alist->lvl <= al) & (r->rweight < 0.6))
473 maxangle = (maxangle - PI/2.)*pow(r->rweight,0.13) + PI/2.;
474 /* sum this node */
475 for (av = at->alist; av != NULL; av = av->next) {
476 double u, v, d, delta_r2, delta_t2;
477 COLOR ct;
478 FVECT uvw[3];
479 /* record access */
480 if (tracktime)
481 av->latick = ambclock;
482 /*
483 * Ambient level test
484 */
485 if (av->lvl > al || /* list sorted, so this works */
486 (av->lvl == al) & (av->weight < 0.9*r->rweight))
487 break;
488 /*
489 * Direction test using unperturbed normal
490 */
491 decodedir(uvw[2], av->ndir);
492 d = DOT(uvw[2], r->ron);
493 if (d <= 0.0) /* >= 90 degrees */
494 continue;
495 delta_r2 = 2.0 - 2.0*d; /* approx. radians^2 */
496 if (delta_r2 >= maxangle*maxangle)
497 continue;
498 /*
499 * Modified ray behind test
500 */
501 VSUB(ck0, r->rop, av->pos);
502 d = DOT(ck0, uvw[2]);
503 if (d < -minarad*ambacc-.001)
504 continue;
505 d /= av->rad[0];
506 delta_t2 = d*d;
507 if (delta_t2 >= ambacc*ambacc)
508 continue;
509 /*
510 * Elliptical radii test based on Hessian
511 */
512 decodedir(uvw[0], av->udir);
513 VCROSS(uvw[1], uvw[2], uvw[0]);
514 d = (u = DOT(ck0, uvw[0])) / av->rad[0];
515 delta_t2 += d*d;
516 d = (v = DOT(ck0, uvw[1])) / av->rad[1];
517 delta_t2 += d*d;
518 if (delta_t2 >= ambacc*ambacc)
519 continue;
520 /*
521 * Test for potential light leak
522 */
523 if (av->corral && plugaleak(r, av, uvw[2], atan2a(v,u)))
524 continue;
525 /*
526 * Extrapolate value and compute final weight (hat function)
527 */
528 if (!extambient(ct, av, r->rop, rn, uvw))
529 continue;
530 d = tfunc(maxangle, sqrt(delta_r2), 0.0) *
531 tfunc(ambacc, sqrt(delta_t2), 0.0);
532 scalecolor(ct, d);
533 addcolor(acol, ct);
534 wsum += d;
535 }
536 return(wsum);
537 }
538
539
540 static int
541 makeambient( /* make a new ambient value for storage */
542 COLOR acol,
543 RAY *r,
544 FVECT rn,
545 int al
546 )
547 {
548 AMBVAL amb;
549 FVECT uvw[3];
550 int i;
551
552 amb.weight = 1.0; /* compute weight */
553 for (i = al; i-- > 0; )
554 amb.weight *= AVGREFL;
555 if (r->rweight < 0.1*amb.weight) /* heuristic override */
556 amb.weight = 1.25*r->rweight;
557 setcolor(acol, AVGREFL, AVGREFL, AVGREFL);
558 /* compute ambient */
559 i = doambient(acol, r, amb.weight,
560 uvw, amb.rad, amb.gpos, amb.gdir, &amb.corral);
561 scalecolor(acol, 1./AVGREFL); /* undo assumed reflectance */
562 if (i <= 0 || amb.rad[0] <= FTINY) /* no Hessian or zero radius */
563 return(i);
564 /* store value */
565 VCOPY(amb.pos, r->rop);
566 amb.ndir = encodedir(r->ron);
567 amb.udir = encodedir(uvw[0]);
568 amb.lvl = al;
569 copycolor(amb.val, acol);
570 /* insert into tree */
571 avsave(&amb); /* and save to file */
572 if (rn != r->ron) { /* texture */
573 VCOPY(uvw[2], r->ron);
574 extambient(acol, &amb, r->rop, rn, uvw);
575 }
576 return(1);
577 }
578
579
580 static int
581 extambient( /* extrapolate value at pv, nv */
582 COLOR cr,
583 AMBVAL *ap,
584 FVECT pv,
585 FVECT nv,
586 FVECT uvw[3]
587 )
588 {
589 const double min_d = 0.05;
590 static FVECT my_uvw[3];
591 FVECT v1;
592 int i;
593 double d = 1.0; /* zeroeth order */
594
595 if (uvw == NULL) { /* need local coordinates? */
596 decodedir(my_uvw[2], ap->ndir);
597 decodedir(my_uvw[0], ap->udir);
598 VCROSS(my_uvw[1], my_uvw[2], my_uvw[0]);
599 uvw = my_uvw;
600 }
601 for (i = 3; i--; ) /* gradient due to translation */
602 d += (pv[i] - ap->pos[i]) *
603 (ap->gpos[0]*uvw[0][i] + ap->gpos[1]*uvw[1][i]);
604
605 VCROSS(v1, uvw[2], nv); /* gradient due to rotation */
606 for (i = 3; i--; )
607 d += v1[i] * (ap->gdir[0]*uvw[0][i] + ap->gdir[1]*uvw[1][i]);
608
609 if (d < min_d) /* should not use if we can avoid it */
610 d = min_d;
611 copycolor(cr, ap->val);
612 scalecolor(cr, d);
613 return(d > min_d);
614 }
615
616
617 static void
618 avinsert( /* insert ambient value in our tree */
619 AMBVAL *av
620 )
621 {
622 AMBTREE *at;
623 AMBVAL *ap;
624 AMBVAL avh;
625 FVECT ck0;
626 double s;
627 int branch;
628 int i;
629
630 if (av->rad[0] <= FTINY)
631 error(CONSISTENCY, "zero ambient radius in avinsert");
632 at = &atrunk;
633 VCOPY(ck0, thescene.cuorg);
634 s = thescene.cusize;
635 while (s*(OCTSCALE/2) > av->rad[1]*ambacc) {
636 if (at->kid == NULL)
637 if ((at->kid = newambtree()) == NULL)
638 error(SYSTEM, "out of memory in avinsert");
639 s *= 0.5;
640 branch = 0;
641 for (i = 0; i < 3; i++)
642 if (av->pos[i] > ck0[i] + s) {
643 ck0[i] += s;
644 branch |= 1 << i;
645 }
646 at = at->kid + branch;
647 }
648 avh.next = at->alist; /* order by increasing level */
649 for (ap = &avh; ap->next != NULL; ap = ap->next)
650 if ( ap->next->lvl > av->lvl ||
651 (ap->next->lvl == av->lvl) &
652 (ap->next->weight <= av->weight) )
653 break;
654 av->next = ap->next;
655 ap->next = (AMBVAL*)av;
656 at->alist = avh.next;
657 }
658
659
660 #else /* ! NEWAMB */
661
662 static double sumambient(COLOR acol, RAY *r, FVECT rn, int al,
663 AMBTREE *at, FVECT c0, double s);
664 static double makeambient(COLOR acol, RAY *r, FVECT rn, int al);
665 static void extambient(COLOR cr, AMBVAL *ap, FVECT pv, FVECT nv);
666
667
668 void
669 multambient( /* compute ambient component & multiply by coef. */
670 COLOR aval,
671 RAY *r,
672 FVECT nrm
673 )
674 {
675 static int rdepth = 0; /* ambient recursion */
676 COLOR acol, caustic;
677 double d, l;
678
679 /* PMAP: Factor in ambient from global photon map (if enabled) and return
680 * as all ambient components accounted for */
681 if (ambPmap(aval, r, rdepth))
682 return;
683
684 /* PMAP: Otherwise factor in ambient from caustic photon map
685 * (ambPmapCaustic() returns zero if caustic photons disabled) and
686 * continue with RADIANCE ambient calculation */
687 copycolor(caustic, aval);
688 ambPmapCaustic(caustic, r, rdepth);
689
690 if (ambdiv <= 0) /* no ambient calculation */
691 goto dumbamb;
692 /* check number of bounces */
693 if (rdepth >= ambounce)
694 goto dumbamb;
695 /* check ambient list */
696 if (ambincl != -1 && r->ro != NULL &&
697 ambincl != inset(ambset, r->ro->omod))
698 goto dumbamb;
699
700 if (ambacc <= FTINY) { /* no ambient storage */
701 copycolor(acol, aval);
702 rdepth++;
703 d = doambient(acol, r, r->rweight, NULL, NULL);
704 rdepth--;
705 if (d <= FTINY)
706 goto dumbamb;
707 copycolor(aval, acol);
708
709 /* PMAP: add in caustic */
710 addcolor(aval, caustic);
711 return;
712 }
713
714 if (tracktime) /* sort to minimize thrashing */
715 sortambvals(0);
716 /* interpolate ambient value */
717 setcolor(acol, 0.0, 0.0, 0.0);
718 d = sumambient(acol, r, nrm, rdepth,
719 &atrunk, thescene.cuorg, thescene.cusize);
720
721 if (d > FTINY) {
722 d = 1.0/d;
723 scalecolor(acol, d);
724 multcolor(aval, acol);
725
726 /* PMAP: add in caustic */
727 addcolor(aval, caustic);
728 return;
729 }
730
731 rdepth++; /* need to cache new value */
732 d = makeambient(acol, r, nrm, rdepth-1);
733 rdepth--;
734
735 if (d > FTINY) {
736 multcolor(aval, acol); /* got new value */
737
738 /* PMAP: add in caustic */
739 addcolor(aval, caustic);
740 return;
741 }
742
743 dumbamb: /* return global value */
744 if ((ambvwt <= 0) | (navsum == 0)) {
745 multcolor(aval, ambval);
746
747 /* PMAP: add in caustic */
748 addcolor(aval, caustic);
749 return;
750 }
751
752 l = bright(ambval); /* average in computations */
753 if (l > FTINY) {
754 d = (log(l)*(double)ambvwt + avsum) /
755 (double)(ambvwt + navsum);
756 d = exp(d) / l;
757 scalecolor(aval, d);
758 multcolor(aval, ambval); /* apply color of ambval */
759 } else {
760 d = exp( avsum / (double)navsum );
761 scalecolor(aval, d); /* neutral color */
762 }
763 }
764
765
766 static double
767 sumambient( /* get interpolated ambient value */
768 COLOR acol,
769 RAY *r,
770 FVECT rn,
771 int al,
772 AMBTREE *at,
773 FVECT c0,
774 double s
775 )
776 {
777 double d, e1, e2, wt, wsum;
778 COLOR ct;
779 FVECT ck0;
780 int i;
781 int j;
782 AMBVAL *av;
783
784 wsum = 0.0;
785 /* do this node */
786 for (av = at->alist; av != NULL; av = av->next) {
787 double rn_dot = -2.0;
788 if (tracktime)
789 av->latick = ambclock;
790 /*
791 * Ambient level test.
792 */
793 if (av->lvl > al || /* list sorted, so this works */
794 (av->lvl == al) & (av->weight < 0.9*r->rweight))
795 break;
796 /*
797 * Ambient radius test.
798 */
799 VSUB(ck0, av->pos, r->rop);
800 e1 = DOT(ck0, ck0) / (av->rad * av->rad);
801 if (e1 > ambacc*ambacc*1.21)
802 continue;
803 /*
804 * Direction test using closest normal.
805 */
806 d = DOT(av->dir, r->ron);
807 if (rn != r->ron) {
808 rn_dot = DOT(av->dir, rn);
809 if (rn_dot > 1.0-FTINY)
810 rn_dot = 1.0-FTINY;
811 if (rn_dot >= d-FTINY) {
812 d = rn_dot;
813 rn_dot = -2.0;
814 }
815 }
816 e2 = (1.0 - d) * r->rweight;
817 if (e2 < 0.0)
818 e2 = 0.0;
819 else if (e1 + e2 > ambacc*ambacc*1.21)
820 continue;
821 /*
822 * Ray behind test.
823 */
824 d = 0.0;
825 for (j = 0; j < 3; j++)
826 d += (r->rop[j] - av->pos[j]) *
827 (av->dir[j] + r->ron[j]);
828 if (d*0.5 < -minarad*ambacc-.001)
829 continue;
830 /*
831 * Jittering final test reduces image artifacts.
832 */
833 e1 = sqrt(e1);
834 e2 = sqrt(e2);
835 wt = e1 + e2;
836 if (wt > ambacc*(.9+.2*urand(9015+samplendx)))
837 continue;
838 /*
839 * Recompute directional error using perturbed normal
840 */
841 if (rn_dot > 0.0) {
842 e2 = sqrt((1.0 - rn_dot)*r->rweight);
843 wt = e1 + e2;
844 }
845 if (wt <= 1e-3)
846 wt = 1e3;
847 else
848 wt = 1.0 / wt;
849 wsum += wt;
850 extambient(ct, av, r->rop, rn);
851 scalecolor(ct, wt);
852 addcolor(acol, ct);
853 }
854 if (at->kid == NULL)
855 return(wsum);
856 /* do children */
857 s *= 0.5;
858 for (i = 0; i < 8; i++) {
859 for (j = 0; j < 3; j++) {
860 ck0[j] = c0[j];
861 if (1<<j & i)
862 ck0[j] += s;
863 if (r->rop[j] < ck0[j] - OCTSCALE*s)
864 break;
865 if (r->rop[j] > ck0[j] + (1.0+OCTSCALE)*s)
866 break;
867 }
868 if (j == 3)
869 wsum += sumambient(acol, r, rn, al,
870 at->kid+i, ck0, s);
871 }
872 return(wsum);
873 }
874
875
876 static double
877 makeambient( /* make a new ambient value for storage */
878 COLOR acol,
879 RAY *r,
880 FVECT rn,
881 int al
882 )
883 {
884 AMBVAL amb;
885 FVECT gp, gd;
886 int i;
887
888 amb.weight = 1.0; /* compute weight */
889 for (i = al; i-- > 0; )
890 amb.weight *= AVGREFL;
891 if (r->rweight < 0.1*amb.weight) /* heuristic override */
892 amb.weight = 1.25*r->rweight;
893 setcolor(acol, AVGREFL, AVGREFL, AVGREFL);
894 /* compute ambient */
895 amb.rad = doambient(acol, r, amb.weight, gp, gd);
896 if (amb.rad <= FTINY) {
897 setcolor(acol, 0.0, 0.0, 0.0);
898 return(0.0);
899 }
900 scalecolor(acol, 1./AVGREFL); /* undo assumed reflectance */
901 /* store value */
902 VCOPY(amb.pos, r->rop);
903 VCOPY(amb.dir, r->ron);
904 amb.lvl = al;
905 copycolor(amb.val, acol);
906 VCOPY(amb.gpos, gp);
907 VCOPY(amb.gdir, gd);
908 /* insert into tree */
909 avsave(&amb); /* and save to file */
910 if (rn != r->ron)
911 extambient(acol, &amb, r->rop, rn); /* texture */
912 return(amb.rad);
913 }
914
915
916 static void
917 extambient( /* extrapolate value at pv, nv */
918 COLOR cr,
919 AMBVAL *ap,
920 FVECT pv,
921 FVECT nv
922 )
923 {
924 FVECT v1;
925 int i;
926 double d;
927
928 d = 1.0; /* zeroeth order */
929 /* gradient due to translation */
930 for (i = 0; i < 3; i++)
931 d += ap->gpos[i]*(pv[i]-ap->pos[i]);
932 /* gradient due to rotation */
933 VCROSS(v1, ap->dir, nv);
934 d += DOT(ap->gdir, v1);
935 if (d <= 0.0) {
936 setcolor(cr, 0.0, 0.0, 0.0);
937 return;
938 }
939 copycolor(cr, ap->val);
940 scalecolor(cr, d);
941 }
942
943
944 static void
945 avinsert( /* insert ambient value in our tree */
946 AMBVAL *av
947 )
948 {
949 AMBTREE *at;
950 AMBVAL *ap;
951 AMBVAL avh;
952 FVECT ck0;
953 double s;
954 int branch;
955 int i;
956
957 if (av->rad <= FTINY)
958 error(CONSISTENCY, "zero ambient radius in avinsert");
959 at = &atrunk;
960 VCOPY(ck0, thescene.cuorg);
961 s = thescene.cusize;
962 while (s*(OCTSCALE/2) > av->rad*ambacc) {
963 if (at->kid == NULL)
964 if ((at->kid = newambtree()) == NULL)
965 error(SYSTEM, "out of memory in avinsert");
966 s *= 0.5;
967 branch = 0;
968 for (i = 0; i < 3; i++)
969 if (av->pos[i] > ck0[i] + s) {
970 ck0[i] += s;
971 branch |= 1 << i;
972 }
973 at = at->kid + branch;
974 }
975 avh.next = at->alist; /* order by increasing level */
976 for (ap = &avh; ap->next != NULL; ap = ap->next)
977 if ( ap->next->lvl > av->lvl ||
978 (ap->next->lvl == av->lvl) &
979 (ap->next->weight <= av->weight) )
980 break;
981 av->next = ap->next;
982 ap->next = (AMBVAL*)av;
983 at->alist = avh.next;
984 }
985
986 #endif /* ! NEWAMB */
987
988 /************* FOLLOWING ROUTINES SAME FOR NEW & OLD METHODS ***************/
989
990 static void
991 initambfile( /* initialize ambient file */
992 int cre8
993 )
994 {
995 extern char *progname, *octname;
996 static char *mybuf = NULL;
997
998 #ifdef F_SETLKW
999 aflock(cre8 ? F_WRLCK : F_RDLCK);
1000 #endif
1001 SET_FILE_BINARY(ambfp);
1002 if (mybuf == NULL)
1003 mybuf = (char *)bmalloc(BUFSIZ+8);
1004 setbuf(ambfp, mybuf);
1005 if (cre8) { /* new file */
1006 newheader("RADIANCE", ambfp);
1007 fprintf(ambfp, "%s -av %g %g %g -aw %d -ab %d -aa %g ",
1008 progname, colval(ambval,RED),
1009 colval(ambval,GRN), colval(ambval,BLU),
1010 ambvwt, ambounce, ambacc);
1011 fprintf(ambfp, "-ad %d -as %d -ar %d ",
1012 ambdiv, ambssamp, ambres);
1013 if (octname != NULL)
1014 fputs(octname, ambfp);
1015 fputc('\n', ambfp);
1016 fprintf(ambfp, "SOFTWARE= %s\n", VersionID);
1017 fputnow(ambfp);
1018 fputformat(AMBFMT, ambfp);
1019 fputc('\n', ambfp);
1020 putambmagic(ambfp);
1021 } else if (checkheader(ambfp, AMBFMT, NULL) < 0 || !hasambmagic(ambfp))
1022 error(USER, "bad ambient file");
1023 }
1024
1025
1026 static void
1027 avsave( /* insert and save an ambient value */
1028 AMBVAL *av
1029 )
1030 {
1031 avstore(av);
1032 if (ambfp == NULL)
1033 return;
1034 if (writambval(av, ambfp) < 0)
1035 goto writerr;
1036 if (++nunflshed >= AMBFLUSH)
1037 if (ambsync() == EOF)
1038 goto writerr;
1039 return;
1040 writerr:
1041 error(SYSTEM, "error writing to ambient file");
1042 }
1043
1044
1045 static AMBVAL *
1046 avstore( /* allocate memory and save aval */
1047 AMBVAL *aval
1048 )
1049 {
1050 AMBVAL *av;
1051 double d;
1052
1053 if ((av = newambval()) == NULL)
1054 error(SYSTEM, "out of memory in avstore");
1055 *av = *aval;
1056 av->latick = ambclock;
1057 av->next = NULL;
1058 nambvals++;
1059 d = bright(av->val);
1060 if (d > FTINY) { /* add to log sum for averaging */
1061 avsum += log(d);
1062 navsum++;
1063 }
1064 avinsert(av); /* insert in our cache tree */
1065 return(av);
1066 }
1067
1068
1069 #define ATALLOCSZ 512 /* #/8 trees to allocate at once */
1070
1071 static AMBTREE *atfreelist = NULL; /* free ambient tree structures */
1072
1073
1074 static AMBTREE *
1075 newambtree(void) /* allocate 8 ambient tree structs */
1076 {
1077 AMBTREE *atp, *upperlim;
1078
1079 if (atfreelist == NULL) { /* get more nodes */
1080 atfreelist = (AMBTREE *)malloc(ATALLOCSZ*8*sizeof(AMBTREE));
1081 if (atfreelist == NULL)
1082 return(NULL);
1083 /* link new free list */
1084 upperlim = atfreelist + 8*(ATALLOCSZ-1);
1085 for (atp = atfreelist; atp < upperlim; atp += 8)
1086 atp->kid = atp + 8;
1087 atp->kid = NULL;
1088 }
1089 atp = atfreelist;
1090 atfreelist = atp->kid;
1091 memset(atp, 0, 8*sizeof(AMBTREE));
1092 return(atp);
1093 }
1094
1095
1096 static void
1097 freeambtree( /* free 8 ambient tree structs */
1098 AMBTREE *atp
1099 )
1100 {
1101 atp->kid = atfreelist;
1102 atfreelist = atp;
1103 }
1104
1105
1106 static void
1107 unloadatree( /* unload an ambient value tree */
1108 AMBTREE *at,
1109 unloadtf_t *f
1110 )
1111 {
1112 AMBVAL *av;
1113 int i;
1114 /* transfer values at this node */
1115 for (av = at->alist; av != NULL; av = at->alist) {
1116 at->alist = av->next;
1117 av->next = NULL;
1118 (*f)(av);
1119 }
1120 if (at->kid == NULL)
1121 return;
1122 for (i = 0; i < 8; i++) /* transfer and free children */
1123 unloadatree(at->kid+i, f);
1124 freeambtree(at->kid);
1125 at->kid = NULL;
1126 }
1127
1128
1129 static struct avl {
1130 AMBVAL *p;
1131 unsigned long t;
1132 } *avlist1; /* ambient value list with ticks */
1133 static AMBVAL **avlist2; /* memory positions for sorting */
1134 static int i_avlist; /* index for lists */
1135
1136 static int alatcmp(const void *av1, const void *av2);
1137
1138 static void
1139 avfree(AMBVAL *av)
1140 {
1141 free(av);
1142 }
1143
1144 static void
1145 av2list(
1146 AMBVAL *av
1147 )
1148 {
1149 #ifdef DEBUG
1150 if (i_avlist >= nambvals)
1151 error(CONSISTENCY, "too many ambient values in av2list1");
1152 #endif
1153 avlist1[i_avlist].p = avlist2[i_avlist] = (AMBVAL*)av;
1154 avlist1[i_avlist++].t = av->latick;
1155 }
1156
1157
1158 static int
1159 alatcmp( /* compare ambient values for MRA */
1160 const void *av1,
1161 const void *av2
1162 )
1163 {
1164 long lc = ((struct avl *)av2)->t - ((struct avl *)av1)->t;
1165 return(lc<0 ? -1 : lc>0 ? 1 : 0);
1166 }
1167
1168
1169 /* GW NOTE 2002/10/3:
1170 * I used to compare AMBVAL pointers, but found that this was the
1171 * cause of a serious consistency error with gcc, since the optimizer
1172 * uses some dangerous trick in pointer subtraction that
1173 * assumes pointers differ by exact struct size increments.
1174 */
1175 static int
1176 aposcmp( /* compare ambient value positions */
1177 const void *avp1,
1178 const void *avp2
1179 )
1180 {
1181 long diff = *(char * const *)avp1 - *(char * const *)avp2;
1182 if (diff < 0)
1183 return(-1);
1184 return(diff > 0);
1185 }
1186
1187
1188 static int
1189 avlmemi( /* find list position from address */
1190 AMBVAL *avaddr
1191 )
1192 {
1193 AMBVAL **avlpp;
1194
1195 avlpp = (AMBVAL **)bsearch(&avaddr, avlist2,
1196 nambvals, sizeof(AMBVAL *), aposcmp);
1197 if (avlpp == NULL)
1198 error(CONSISTENCY, "address not found in avlmemi");
1199 return(avlpp - avlist2);
1200 }
1201
1202
1203 static void
1204 sortambvals( /* resort ambient values */
1205 int always
1206 )
1207 {
1208 AMBTREE oldatrunk;
1209 AMBVAL tav, *tap, *pnext;
1210 int i, j;
1211 /* see if it's time yet */
1212 if (!always && (ambclock++ < lastsort+sortintvl ||
1213 nambvals < SORT_THRESH))
1214 return;
1215 /*
1216 * The idea here is to minimize memory thrashing
1217 * in VM systems by improving reference locality.
1218 * We do this by periodically sorting our stored ambient
1219 * values in memory in order of most recently to least
1220 * recently accessed. This ordering was chosen so that new
1221 * ambient values (which tend to be less important) go into
1222 * higher memory with the infrequently accessed values.
1223 * Since we expect our values to need sorting less
1224 * frequently as the process continues, we double our
1225 * waiting interval after each call.
1226 * This routine is also called by setambacc() with
1227 * the "always" parameter set to 1 so that the ambient
1228 * tree will be rebuilt with the new accuracy parameter.
1229 */
1230 if (tracktime) { /* allocate pointer arrays to sort */
1231 avlist2 = (AMBVAL **)malloc(nambvals*sizeof(AMBVAL *));
1232 avlist1 = (struct avl *)malloc(nambvals*sizeof(struct avl));
1233 } else {
1234 avlist2 = NULL;
1235 avlist1 = NULL;
1236 }
1237 if (avlist1 == NULL) { /* no time tracking -- rebuild tree? */
1238 if (avlist2 != NULL)
1239 free(avlist2);
1240 if (always) { /* rebuild without sorting */
1241 oldatrunk = atrunk;
1242 atrunk.alist = NULL;
1243 atrunk.kid = NULL;
1244 unloadatree(&oldatrunk, avinsert);
1245 }
1246 } else { /* sort memory by last access time */
1247 /*
1248 * Sorting memory is tricky because it isn't contiguous.
1249 * We have to sort an array of pointers by MRA and also
1250 * by memory position. We then copy values in "loops"
1251 * to minimize memory hits. Nevertheless, we will visit
1252 * everyone at least twice, and this is an expensive process
1253 * when we're thrashing, which is when we need to do it.
1254 */
1255 #ifdef DEBUG
1256 sprintf(errmsg, "sorting %u ambient values at ambclock=%lu...",
1257 nambvals, ambclock);
1258 eputs(errmsg);
1259 #endif
1260 i_avlist = 0;
1261 unloadatree(&atrunk, av2list); /* empty current tree */
1262 #ifdef DEBUG
1263 if (i_avlist < nambvals)
1264 error(CONSISTENCY, "missing ambient values in sortambvals");
1265 #endif
1266 qsort(avlist1, nambvals, sizeof(struct avl), alatcmp);
1267 qsort(avlist2, nambvals, sizeof(AMBVAL *), aposcmp);
1268 for (i = 0; i < nambvals; i++) {
1269 if (avlist1[i].p == NULL)
1270 continue;
1271 tap = avlist2[i];
1272 tav = *tap;
1273 for (j = i; (pnext = avlist1[j].p) != tap;
1274 j = avlmemi(pnext)) {
1275 *(avlist2[j]) = *pnext;
1276 avinsert(avlist2[j]);
1277 avlist1[j].p = NULL;
1278 }
1279 *(avlist2[j]) = tav;
1280 avinsert(avlist2[j]);
1281 avlist1[j].p = NULL;
1282 }
1283 free(avlist1);
1284 free(avlist2);
1285 /* compute new sort interval */
1286 sortintvl = ambclock - lastsort;
1287 if (sortintvl >= MAX_SORT_INTVL/2)
1288 sortintvl = MAX_SORT_INTVL;
1289 else
1290 sortintvl <<= 1; /* wait twice as long next */
1291 #ifdef DEBUG
1292 eputs("done\n");
1293 #endif
1294 }
1295 if (ambclock >= MAXACLOCK)
1296 ambclock = MAXACLOCK/2;
1297 lastsort = ambclock;
1298 }
1299
1300
1301 #ifdef F_SETLKW
1302
1303 static void
1304 aflock( /* lock/unlock ambient file */
1305 int typ
1306 )
1307 {
1308 static struct flock fls; /* static so initialized to zeroes */
1309
1310 if (typ == fls.l_type) /* already called? */
1311 return;
1312 fls.l_type = typ;
1313 if (fcntl(fileno(ambfp), F_SETLKW, &fls) < 0)
1314 error(SYSTEM, "cannot (un)lock ambient file");
1315 }
1316
1317
1318 int
1319 ambsync(void) /* synchronize ambient file */
1320 {
1321 long flen;
1322 AMBVAL avs;
1323 int n;
1324
1325 if (ambfp == NULL) /* no ambient file? */
1326 return(0);
1327 /* gain appropriate access */
1328 aflock(nunflshed ? F_WRLCK : F_RDLCK);
1329 /* see if file has grown */
1330 if ((flen = lseek(fileno(ambfp), (off_t)0, SEEK_END)) < 0)
1331 goto seekerr;
1332 if ((n = flen - lastpos) > 0) { /* file has grown */
1333 if (ambinp == NULL) { /* get new file pointer */
1334 ambinp = fopen(ambfile, "rb");
1335 if (ambinp == NULL)
1336 error(SYSTEM, "fopen failed in ambsync");
1337 }
1338 if (fseek(ambinp, lastpos, SEEK_SET) < 0)
1339 goto seekerr;
1340 while (n >= AMBVALSIZ) { /* load contributed values */
1341 if (!readambval(&avs, ambinp)) {
1342 sprintf(errmsg,
1343 "ambient file \"%s\" corrupted near character %ld",
1344 ambfile, flen - n);
1345 error(WARNING, errmsg);
1346 break;
1347 }
1348 avstore(&avs);
1349 n -= AMBVALSIZ;
1350 }
1351 lastpos = flen - n; /* check alignment */
1352 if (n && lseek(fileno(ambfp), (off_t)lastpos, SEEK_SET) < 0)
1353 goto seekerr;
1354 }
1355 n = fflush(ambfp); /* calls write() at last */
1356 lastpos += (long)nunflshed*AMBVALSIZ;
1357 aflock(F_UNLCK); /* release file */
1358 nunflshed = 0;
1359 return(n);
1360 seekerr:
1361 error(SYSTEM, "seek failed in ambsync");
1362 return(EOF); /* pro forma return */
1363 }
1364
1365 #else /* ! F_SETLKW */
1366
1367 int
1368 ambsync(void) /* flush ambient file */
1369 {
1370 if (ambfp == NULL)
1371 return(0);
1372 nunflshed = 0;
1373 return(fflush(ambfp));
1374 }
1375
1376 #endif /* ! F_SETLKW */