ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/ambient.c
Revision: 2.104
Committed: Fri Jan 27 22:00:49 2017 UTC (7 years, 3 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R1
Changes since 2.103: +9 -15 lines
Log Message:
Created independent ambient file pointer to resolve some race issues

File Contents

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