ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/ambient.c
Revision: 2.102
Committed: Sun Apr 24 16:21:32 2016 UTC (8 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.101: +10 -10 lines
Log Message:
Minor changes

File Contents

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