ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/ambient.c
Revision: 2.101
Committed: Sun Mar 6 01:13:17 2016 UTC (8 years, 2 months ago) by schorsch
Content type: text/plain
Branch: MAIN
Changes since 2.100: +1 -3 lines
Log Message:
Prepare for SCons build on Win32 and Win64

File Contents

# Content
1 static const char RCSid[] = "$Id: ambient.c,v 2.100 2016/03/03 22:47:21 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 #define freeav(av) free((void *)av);
80
81 static void initambfile(int creat);
82 static void avsave(AMBVAL *av);
83 static AMBVAL *avstore(AMBVAL *aval);
84 static AMBTREE *newambtree(void);
85 static void freeambtree(AMBTREE *atp);
86
87 typedef void unloadtf_t(AMBVAL *);
88 static unloadtf_t avinsert;
89 static unloadtf_t av2list;
90 static unloadtf_t avfree;
91 static void unloadatree(AMBTREE *at, unloadtf_t *f);
92
93 static int aposcmp(const void *avp1, const void *avp2);
94 static int avlmemi(AMBVAL *avaddr);
95 static void sortambvals(int always);
96
97 #ifdef F_SETLKW
98 static void aflock(int typ);
99 #endif
100
101
102 void
103 setambres( /* set ambient resolution */
104 int ar
105 )
106 {
107 ambres = ar < 0 ? 0 : ar; /* may be done already */
108 /* set min & max radii */
109 if (ar <= 0) {
110 minarad = 0;
111 maxarad = thescene.cusize*0.2;
112 } else {
113 minarad = thescene.cusize / ar;
114 maxarad = 64.0 * minarad; /* heuristic */
115 if (maxarad > thescene.cusize*0.2)
116 maxarad = thescene.cusize*0.2;
117 }
118 if (minarad <= FTINY)
119 minarad = 10.0*FTINY;
120 if (maxarad <= minarad)
121 maxarad = 64.0 * minarad;
122 }
123
124
125 void
126 setambacc( /* set ambient accuracy */
127 double newa
128 )
129 {
130 static double olda; /* remember previous setting here */
131
132 newa *= (newa > 0);
133 if (fabs(newa - olda) >= .05*(newa + olda)) {
134 ambacc = newa;
135 if (nambvals > 0)
136 sortambvals(1); /* rebuild tree */
137 }
138 }
139
140
141 void
142 setambient(void) /* initialize calculation */
143 {
144 int readonly = 0;
145 long flen;
146 AMBVAL amb;
147 /* make sure we're fresh */
148 ambdone();
149 /* init ambient limits */
150 setambres(ambres);
151 setambacc(ambacc);
152 if (ambfile == NULL || !ambfile[0])
153 return;
154 if (ambacc <= FTINY) {
155 sprintf(errmsg, "zero ambient accuracy so \"%s\" not opened",
156 ambfile);
157 error(WARNING, errmsg);
158 return;
159 }
160 /* open ambient file */
161 if ((ambfp = fopen(ambfile, "r+")) == NULL)
162 readonly = (ambfp = fopen(ambfile, "r")) != NULL;
163 if (ambfp != NULL) {
164 initambfile(0); /* file exists */
165 lastpos = ftell(ambfp);
166 while (readambval(&amb, ambfp))
167 avstore(&amb);
168 nambshare = nambvals; /* share loaded values */
169 if (readonly) {
170 sprintf(errmsg,
171 "loaded %u values from read-only ambient file",
172 nambvals);
173 error(WARNING, errmsg);
174 fclose(ambfp); /* close file so no writes */
175 ambfp = NULL;
176 return; /* avoid ambsync() */
177 }
178 /* align file pointer */
179 lastpos += (long)nambvals*AMBVALSIZ;
180 flen = lseek(fileno(ambfp), (off_t)0, SEEK_END);
181 if (flen != lastpos) {
182 sprintf(errmsg,
183 "ignoring last %ld values in ambient file (corrupted)",
184 (flen - lastpos)/AMBVALSIZ);
185 error(WARNING, errmsg);
186 fseek(ambfp, lastpos, SEEK_SET);
187 ftruncate(fileno(ambfp), (off_t)lastpos);
188 }
189 } else if ((ambfp = fopen(ambfile, "w+")) != NULL) {
190 initambfile(1); /* else create new file */
191 fflush(ambfp);
192 lastpos = ftell(ambfp);
193 } else {
194 sprintf(errmsg, "cannot open ambient file \"%s\"", ambfile);
195 error(SYSTEM, errmsg);
196 }
197 #ifdef F_SETLKW
198 aflock(F_UNLCK); /* release file */
199 #endif
200 }
201
202
203 void
204 ambdone(void) /* close ambient file and free memory */
205 {
206 if (ambfp != NULL) { /* close ambient file */
207 ambsync();
208 fclose(ambfp);
209 ambfp = NULL;
210 if (ambinp != NULL) {
211 fclose(ambinp);
212 ambinp = NULL;
213 }
214 lastpos = -1;
215 }
216 /* free ambient tree */
217 unloadatree(&atrunk, avfree);
218 /* reset state variables */
219 avsum = 0.;
220 navsum = 0;
221 nambvals = 0;
222 nambshare = 0;
223 ambclock = 0;
224 lastsort = 0;
225 sortintvl = SORT_INTVL;
226 }
227
228
229 void
230 ambnotify( /* record new modifier */
231 OBJECT obj
232 )
233 {
234 static int hitlimit = 0;
235 OBJREC *o;
236 char **amblp;
237
238 if (obj == OVOID) { /* starting over */
239 ambset[0] = 0;
240 hitlimit = 0;
241 return;
242 }
243 o = objptr(obj);
244 if (hitlimit || !ismodifier(o->otype))
245 return;
246 for (amblp = amblist; *amblp != NULL; amblp++)
247 if (!strcmp(o->oname, *amblp)) {
248 if (ambset[0] >= MAXASET) {
249 error(WARNING, "too many modifiers in ambient list");
250 hitlimit++;
251 return; /* should this be fatal? */
252 }
253 insertelem(ambset, obj);
254 return;
255 }
256 }
257
258 /************ THE FOLLOWING ROUTINES DIFFER BETWEEN NEW & OLD ***************/
259
260 #ifndef OLDAMB
261
262 #define tfunc(lwr, x, upr) (((x)-(lwr))/((upr)-(lwr)))
263
264 static int plugaleak(RAY *r, AMBVAL *ap, FVECT anorm, double ang);
265 static double sumambient(COLOR acol, RAY *r, FVECT rn, int al,
266 AMBTREE *at, FVECT c0, double s);
267 static int makeambient(COLOR acol, RAY *r, FVECT rn, int al);
268 static int extambient(COLOR cr, AMBVAL *ap, FVECT pv, FVECT nv,
269 FVECT uvw[3]);
270
271 void
272 multambient( /* compute ambient component & multiply by coef. */
273 COLOR aval,
274 RAY *r,
275 FVECT nrm
276 )
277 {
278 static int rdepth = 0; /* ambient recursion */
279 COLOR acol, caustic;
280 int ok;
281 double d, l;
282
283 /* PMAP: Factor in ambient from photon map, if enabled and ray is
284 * ambient. Return as all ambient components accounted for, else
285 * continue. */
286 if (ambPmap(aval, r, rdepth))
287 return;
288
289 /* PMAP: Factor in specular-diffuse ambient (caustics) from photon
290 * map, if enabled and ray is primary, else caustic is zero. Continue
291 * with RADIANCE ambient calculation */
292 copycolor(caustic, aval);
293 ambPmapCaustic(caustic, r, rdepth);
294
295 if (ambdiv <= 0) /* no ambient calculation */
296 goto dumbamb;
297 /* check number of bounces */
298 if (rdepth >= ambounce)
299 goto dumbamb;
300 /* check ambient list */
301 if (ambincl != -1 && r->ro != NULL &&
302 ambincl != inset(ambset, r->ro->omod))
303 goto dumbamb;
304
305 if (ambacc <= FTINY) { /* no ambient storage */
306 copycolor(acol, aval);
307 rdepth++;
308 ok = doambient(acol, r, r->rweight,
309 NULL, NULL, NULL, NULL, NULL);
310 rdepth--;
311 if (!ok)
312 goto dumbamb;
313 copycolor(aval, acol);
314
315 /* PMAP: add in caustic */
316 addcolor(aval, caustic);
317 return;
318 }
319
320 if (tracktime) /* sort to minimize thrashing */
321 sortambvals(0);
322 /* interpolate ambient value */
323 setcolor(acol, 0.0, 0.0, 0.0);
324 d = sumambient(acol, r, nrm, rdepth,
325 &atrunk, thescene.cuorg, thescene.cusize);
326
327 if (d > FTINY) {
328 d = 1.0/d;
329 scalecolor(acol, d);
330 multcolor(aval, acol);
331
332 /* PMAP: add in caustic */
333 addcolor(aval, caustic);
334 return;
335 }
336
337 rdepth++; /* need to cache new value */
338 ok = makeambient(acol, r, nrm, rdepth-1);
339 rdepth--;
340
341 if (ok) {
342 multcolor(aval, acol); /* computed new value */
343
344 /* PMAP: add in caustic */
345 addcolor(aval, caustic);
346 return;
347 }
348
349 dumbamb: /* return global value */
350 if ((ambvwt <= 0) | (navsum == 0)) {
351 multcolor(aval, ambval);
352
353 /* PMAP: add in caustic */
354 addcolor(aval, caustic);
355 return;
356 }
357
358 l = bright(ambval); /* average in computations */
359 if (l > FTINY) {
360 d = (log(l)*(double)ambvwt + avsum) /
361 (double)(ambvwt + navsum);
362 d = exp(d) / l;
363 scalecolor(aval, d);
364 multcolor(aval, ambval); /* apply color of ambval */
365 } else {
366 d = exp( avsum / (double)navsum );
367 scalecolor(aval, d); /* neutral color */
368 }
369 }
370
371
372 /* Plug a potential leak where ambient cache value is occluded */
373 static int
374 plugaleak(RAY *r, AMBVAL *ap, FVECT anorm, double ang)
375 {
376 const double cost70sq = 0.1169778; /* cos(70deg)^2 */
377 RAY rtst;
378 FVECT vdif;
379 double normdot, ndotd, nadotd;
380 double a, b, c, t[2];
381
382 ang += 2.*PI*(ang < 0); /* check direction flags */
383 if ( !(ap->corral>>(int)(ang*(16./PI)) & 1) )
384 return(0);
385 /*
386 * Generate test ray, targeting 20 degrees above sample point plane
387 * along surface normal from cache position. This should be high
388 * enough to miss local geometry we don't really care about.
389 */
390 VSUB(vdif, ap->pos, r->rop);
391 normdot = DOT(anorm, r->ron);
392 ndotd = DOT(vdif, r->ron);
393 nadotd = DOT(vdif, anorm);
394 a = normdot*normdot - cost70sq;
395 b = 2.0*(normdot*ndotd - nadotd*cost70sq);
396 c = ndotd*ndotd - DOT(vdif,vdif)*cost70sq;
397 if (quadratic(t, a, b, c) != 2)
398 return(1); /* should rarely happen */
399 if (t[1] <= FTINY)
400 return(0); /* should fail behind test */
401 rayorigin(&rtst, SHADOW, r, NULL);
402 VSUM(rtst.rdir, vdif, anorm, t[1]); /* further dist. > plane */
403 rtst.rmax = normalize(rtst.rdir); /* short ray test */
404 while (localhit(&rtst, &thescene)) { /* check for occluder */
405 if (rtst.ro->omod != OVOID &&
406 (rtst.clipset == NULL ||
407 !inset(rtst.clipset, rtst.ro->omod)))
408 return(1); /* plug light leak */
409 VCOPY(rtst.rorg, rtst.rop); /* skip invisible surface */
410 rtst.rmax -= rtst.rot;
411 rayclear(&rtst);
412 }
413 return(0); /* seems we're OK */
414 }
415
416
417 static double
418 sumambient( /* get interpolated ambient value */
419 COLOR acol,
420 RAY *r,
421 FVECT rn,
422 int al,
423 AMBTREE *at,
424 FVECT c0,
425 double s
426 )
427 { /* initial limit is 10 degrees plus ambacc radians */
428 const double minangle = 10.0 * PI/180.;
429 double maxangle = minangle + ambacc;
430 double wsum = 0.0;
431 FVECT ck0;
432 int i, j;
433 AMBVAL *av;
434
435 if (at->kid != NULL) { /* sum children first */
436 s *= 0.5;
437 for (i = 0; i < 8; i++) {
438 for (j = 0; j < 3; j++) {
439 ck0[j] = c0[j];
440 if (1<<j & i)
441 ck0[j] += s;
442 if (r->rop[j] < ck0[j] - OCTSCALE*s)
443 break;
444 if (r->rop[j] > ck0[j] + (1.0+OCTSCALE)*s)
445 break;
446 }
447 if (j == 3)
448 wsum += sumambient(acol, r, rn, al,
449 at->kid+i, ck0, s);
450 }
451 /* good enough? */
452 if (wsum >= 0.05 && s > minarad*10.0)
453 return(wsum);
454 }
455 /* adjust maximum angle */
456 if (at->alist != NULL && (at->alist->lvl <= al) & (r->rweight < 0.6))
457 maxangle = (maxangle - PI/2.)*pow(r->rweight,0.13) + PI/2.;
458 /* sum this node */
459 for (av = at->alist; av != NULL; av = av->next) {
460 double u, v, d, delta_r2, delta_t2;
461 COLOR ct;
462 FVECT uvw[3];
463 /* record access */
464 if (tracktime)
465 av->latick = ambclock;
466 /*
467 * Ambient level test
468 */
469 if (av->lvl > al || /* list sorted, so this works */
470 (av->lvl == al) & (av->weight < 0.9*r->rweight))
471 break;
472 /*
473 * Direction test using unperturbed normal
474 */
475 decodedir(uvw[2], av->ndir);
476 d = DOT(uvw[2], r->ron);
477 if (d <= 0.0) /* >= 90 degrees */
478 continue;
479 delta_r2 = 2.0 - 2.0*d; /* approx. radians^2 */
480 if (delta_r2 >= maxangle*maxangle)
481 continue;
482 /*
483 * Modified ray behind test
484 */
485 VSUB(ck0, r->rop, av->pos);
486 d = DOT(ck0, uvw[2]);
487 if (d < -minarad*ambacc-.001)
488 continue;
489 d /= av->rad[0];
490 delta_t2 = d*d;
491 if (delta_t2 >= ambacc*ambacc)
492 continue;
493 /*
494 * Elliptical radii test based on Hessian
495 */
496 decodedir(uvw[0], av->udir);
497 VCROSS(uvw[1], uvw[2], uvw[0]);
498 d = (u = DOT(ck0, uvw[0])) / av->rad[0];
499 delta_t2 += d*d;
500 d = (v = DOT(ck0, uvw[1])) / av->rad[1];
501 delta_t2 += d*d;
502 if (delta_t2 >= ambacc*ambacc)
503 continue;
504 /*
505 * Test for potential light leak
506 */
507 if (av->corral && plugaleak(r, av, uvw[2], atan2a(v,u)))
508 continue;
509 /*
510 * Extrapolate value and compute final weight (hat function)
511 */
512 if (!extambient(ct, av, r->rop, rn, uvw))
513 continue;
514 d = tfunc(maxangle, sqrt(delta_r2), 0.0) *
515 tfunc(ambacc, sqrt(delta_t2), 0.0);
516 scalecolor(ct, d);
517 addcolor(acol, ct);
518 wsum += d;
519 }
520 return(wsum);
521 }
522
523
524 static int
525 makeambient( /* make a new ambient value for storage */
526 COLOR acol,
527 RAY *r,
528 FVECT rn,
529 int al
530 )
531 {
532 AMBVAL amb;
533 FVECT uvw[3];
534 int i;
535
536 amb.weight = 1.0; /* compute weight */
537 for (i = al; i-- > 0; )
538 amb.weight *= AVGREFL;
539 if (r->rweight < 0.1*amb.weight) /* heuristic override */
540 amb.weight = 1.25*r->rweight;
541 setcolor(acol, AVGREFL, AVGREFL, AVGREFL);
542 /* compute ambient */
543 i = doambient(acol, r, amb.weight,
544 uvw, amb.rad, amb.gpos, amb.gdir, &amb.corral);
545 scalecolor(acol, 1./AVGREFL); /* undo assumed reflectance */
546 if (i <= 0 || amb.rad[0] <= FTINY) /* no Hessian or zero radius */
547 return(i);
548 /* store value */
549 VCOPY(amb.pos, r->rop);
550 amb.ndir = encodedir(r->ron);
551 amb.udir = encodedir(uvw[0]);
552 amb.lvl = al;
553 copycolor(amb.val, acol);
554 /* insert into tree */
555 avsave(&amb); /* and save to file */
556 if (rn != r->ron) { /* texture */
557 VCOPY(uvw[2], r->ron);
558 extambient(acol, &amb, r->rop, rn, uvw);
559 }
560 return(1);
561 }
562
563
564 static int
565 extambient( /* extrapolate value at pv, nv */
566 COLOR cr,
567 AMBVAL *ap,
568 FVECT pv,
569 FVECT nv,
570 FVECT uvw[3]
571 )
572 {
573 const double min_d = 0.05;
574 static FVECT my_uvw[3];
575 FVECT v1;
576 int i;
577 double d = 1.0; /* zeroeth order */
578
579 if (uvw == NULL) { /* need local coordinates? */
580 decodedir(my_uvw[2], ap->ndir);
581 decodedir(my_uvw[0], ap->udir);
582 VCROSS(my_uvw[1], my_uvw[2], my_uvw[0]);
583 uvw = my_uvw;
584 }
585 for (i = 3; i--; ) /* gradient due to translation */
586 d += (pv[i] - ap->pos[i]) *
587 (ap->gpos[0]*uvw[0][i] + ap->gpos[1]*uvw[1][i]);
588
589 VCROSS(v1, uvw[2], nv); /* gradient due to rotation */
590 for (i = 3; i--; )
591 d += v1[i] * (ap->gdir[0]*uvw[0][i] + ap->gdir[1]*uvw[1][i]);
592
593 if (d < min_d) /* should not use if we can avoid it */
594 d = min_d;
595 copycolor(cr, ap->val);
596 scalecolor(cr, d);
597 return(d > min_d);
598 }
599
600
601 static void
602 avinsert( /* insert ambient value in our tree */
603 AMBVAL *av
604 )
605 {
606 AMBTREE *at;
607 AMBVAL *ap;
608 AMBVAL avh;
609 FVECT ck0;
610 double s;
611 int branch;
612 int i;
613
614 if (av->rad[0] <= FTINY)
615 error(CONSISTENCY, "zero ambient radius in avinsert");
616 at = &atrunk;
617 VCOPY(ck0, thescene.cuorg);
618 s = thescene.cusize;
619 while (s*(OCTSCALE/2) > av->rad[1]*ambacc) {
620 if (at->kid == NULL)
621 if ((at->kid = newambtree()) == NULL)
622 error(SYSTEM, "out of memory in avinsert");
623 s *= 0.5;
624 branch = 0;
625 for (i = 0; i < 3; i++)
626 if (av->pos[i] > ck0[i] + s) {
627 ck0[i] += s;
628 branch |= 1 << i;
629 }
630 at = at->kid + branch;
631 }
632 avh.next = at->alist; /* order by increasing level */
633 for (ap = &avh; ap->next != NULL; ap = ap->next)
634 if ( ap->next->lvl > av->lvl ||
635 (ap->next->lvl == av->lvl) &
636 (ap->next->weight <= av->weight) )
637 break;
638 av->next = ap->next;
639 ap->next = (AMBVAL*)av;
640 at->alist = avh.next;
641 }
642
643
644 #else /* ! NEWAMB */
645
646 static double sumambient(COLOR acol, RAY *r, FVECT rn, int al,
647 AMBTREE *at, FVECT c0, double s);
648 static double makeambient(COLOR acol, RAY *r, FVECT rn, int al);
649 static void extambient(COLOR cr, AMBVAL *ap, FVECT pv, FVECT nv);
650
651
652 void
653 multambient( /* compute ambient component & multiply by coef. */
654 COLOR aval,
655 RAY *r,
656 FVECT nrm
657 )
658 {
659 static int rdepth = 0; /* ambient recursion */
660 COLOR acol, caustic;
661 double d, l;
662
663 /* PMAP: Factor in ambient from global photon map (if enabled) and return
664 * as all ambient components accounted for */
665 if (ambPmap(aval, r, rdepth))
666 return;
667
668 /* PMAP: Otherwise factor in ambient from caustic photon map
669 * (ambPmapCaustic() returns zero if caustic photons disabled) and
670 * continue with RADIANCE ambient calculation */
671 copycolor(caustic, aval);
672 ambPmapCaustic(caustic, r, rdepth);
673
674 if (ambdiv <= 0) /* no ambient calculation */
675 goto dumbamb;
676 /* check number of bounces */
677 if (rdepth >= ambounce)
678 goto dumbamb;
679 /* check ambient list */
680 if (ambincl != -1 && r->ro != NULL &&
681 ambincl != inset(ambset, r->ro->omod))
682 goto dumbamb;
683
684 if (ambacc <= FTINY) { /* no ambient storage */
685 copycolor(acol, aval);
686 rdepth++;
687 d = doambient(acol, r, r->rweight, NULL, NULL);
688 rdepth--;
689 if (d <= FTINY)
690 goto dumbamb;
691 copycolor(aval, acol);
692
693 /* PMAP: add in caustic */
694 addcolor(aval, caustic);
695 return;
696 }
697
698 if (tracktime) /* sort to minimize thrashing */
699 sortambvals(0);
700 /* interpolate ambient value */
701 setcolor(acol, 0.0, 0.0, 0.0);
702 d = sumambient(acol, r, nrm, rdepth,
703 &atrunk, thescene.cuorg, thescene.cusize);
704
705 if (d > FTINY) {
706 d = 1.0/d;
707 scalecolor(acol, d);
708 multcolor(aval, acol);
709
710 /* PMAP: add in caustic */
711 addcolor(aval, caustic);
712 return;
713 }
714
715 rdepth++; /* need to cache new value */
716 d = makeambient(acol, r, nrm, rdepth-1);
717 rdepth--;
718
719 if (d > FTINY) {
720 multcolor(aval, acol); /* got new value */
721
722 /* PMAP: add in caustic */
723 addcolor(aval, caustic);
724 return;
725 }
726
727 dumbamb: /* return global value */
728 if ((ambvwt <= 0) | (navsum == 0)) {
729 multcolor(aval, ambval);
730
731 /* PMAP: add in caustic */
732 addcolor(aval, caustic);
733 return;
734 }
735
736 l = bright(ambval); /* average in computations */
737 if (l > FTINY) {
738 d = (log(l)*(double)ambvwt + avsum) /
739 (double)(ambvwt + navsum);
740 d = exp(d) / l;
741 scalecolor(aval, d);
742 multcolor(aval, ambval); /* apply color of ambval */
743 } else {
744 d = exp( avsum / (double)navsum );
745 scalecolor(aval, d); /* neutral color */
746 }
747 }
748
749
750 static double
751 sumambient( /* get interpolated ambient value */
752 COLOR acol,
753 RAY *r,
754 FVECT rn,
755 int al,
756 AMBTREE *at,
757 FVECT c0,
758 double s
759 )
760 {
761 double d, e1, e2, wt, wsum;
762 COLOR ct;
763 FVECT ck0;
764 int i;
765 int j;
766 AMBVAL *av;
767
768 wsum = 0.0;
769 /* do this node */
770 for (av = at->alist; av != NULL; av = av->next) {
771 double rn_dot = -2.0;
772 if (tracktime)
773 av->latick = ambclock;
774 /*
775 * Ambient level test.
776 */
777 if (av->lvl > al || /* list sorted, so this works */
778 (av->lvl == al) & (av->weight < 0.9*r->rweight))
779 break;
780 /*
781 * Ambient radius test.
782 */
783 VSUB(ck0, av->pos, r->rop);
784 e1 = DOT(ck0, ck0) / (av->rad * av->rad);
785 if (e1 > ambacc*ambacc*1.21)
786 continue;
787 /*
788 * Direction test using closest normal.
789 */
790 d = DOT(av->dir, r->ron);
791 if (rn != r->ron) {
792 rn_dot = DOT(av->dir, rn);
793 if (rn_dot > 1.0-FTINY)
794 rn_dot = 1.0-FTINY;
795 if (rn_dot >= d-FTINY) {
796 d = rn_dot;
797 rn_dot = -2.0;
798 }
799 }
800 e2 = (1.0 - d) * r->rweight;
801 if (e2 < 0.0)
802 e2 = 0.0;
803 else if (e1 + e2 > ambacc*ambacc*1.21)
804 continue;
805 /*
806 * Ray behind test.
807 */
808 d = 0.0;
809 for (j = 0; j < 3; j++)
810 d += (r->rop[j] - av->pos[j]) *
811 (av->dir[j] + r->ron[j]);
812 if (d*0.5 < -minarad*ambacc-.001)
813 continue;
814 /*
815 * Jittering final test reduces image artifacts.
816 */
817 e1 = sqrt(e1);
818 e2 = sqrt(e2);
819 wt = e1 + e2;
820 if (wt > ambacc*(.9+.2*urand(9015+samplendx)))
821 continue;
822 /*
823 * Recompute directional error using perturbed normal
824 */
825 if (rn_dot > 0.0) {
826 e2 = sqrt((1.0 - rn_dot)*r->rweight);
827 wt = e1 + e2;
828 }
829 if (wt <= 1e-3)
830 wt = 1e3;
831 else
832 wt = 1.0 / wt;
833 wsum += wt;
834 extambient(ct, av, r->rop, rn);
835 scalecolor(ct, wt);
836 addcolor(acol, ct);
837 }
838 if (at->kid == NULL)
839 return(wsum);
840 /* do children */
841 s *= 0.5;
842 for (i = 0; i < 8; i++) {
843 for (j = 0; j < 3; j++) {
844 ck0[j] = c0[j];
845 if (1<<j & i)
846 ck0[j] += s;
847 if (r->rop[j] < ck0[j] - OCTSCALE*s)
848 break;
849 if (r->rop[j] > ck0[j] + (1.0+OCTSCALE)*s)
850 break;
851 }
852 if (j == 3)
853 wsum += sumambient(acol, r, rn, al,
854 at->kid+i, ck0, s);
855 }
856 return(wsum);
857 }
858
859
860 static double
861 makeambient( /* make a new ambient value for storage */
862 COLOR acol,
863 RAY *r,
864 FVECT rn,
865 int al
866 )
867 {
868 AMBVAL amb;
869 FVECT gp, gd;
870 int i;
871
872 amb.weight = 1.0; /* compute weight */
873 for (i = al; i-- > 0; )
874 amb.weight *= AVGREFL;
875 if (r->rweight < 0.1*amb.weight) /* heuristic override */
876 amb.weight = 1.25*r->rweight;
877 setcolor(acol, AVGREFL, AVGREFL, AVGREFL);
878 /* compute ambient */
879 amb.rad = doambient(acol, r, amb.weight, gp, gd);
880 if (amb.rad <= FTINY) {
881 setcolor(acol, 0.0, 0.0, 0.0);
882 return(0.0);
883 }
884 scalecolor(acol, 1./AVGREFL); /* undo assumed reflectance */
885 /* store value */
886 VCOPY(amb.pos, r->rop);
887 VCOPY(amb.dir, r->ron);
888 amb.lvl = al;
889 copycolor(amb.val, acol);
890 VCOPY(amb.gpos, gp);
891 VCOPY(amb.gdir, gd);
892 /* insert into tree */
893 avsave(&amb); /* and save to file */
894 if (rn != r->ron)
895 extambient(acol, &amb, r->rop, rn); /* texture */
896 return(amb.rad);
897 }
898
899
900 static void
901 extambient( /* extrapolate value at pv, nv */
902 COLOR cr,
903 AMBVAL *ap,
904 FVECT pv,
905 FVECT nv
906 )
907 {
908 FVECT v1;
909 int i;
910 double d;
911
912 d = 1.0; /* zeroeth order */
913 /* gradient due to translation */
914 for (i = 0; i < 3; i++)
915 d += ap->gpos[i]*(pv[i]-ap->pos[i]);
916 /* gradient due to rotation */
917 VCROSS(v1, ap->dir, nv);
918 d += DOT(ap->gdir, v1);
919 if (d <= 0.0) {
920 setcolor(cr, 0.0, 0.0, 0.0);
921 return;
922 }
923 copycolor(cr, ap->val);
924 scalecolor(cr, d);
925 }
926
927
928 static void
929 avinsert( /* insert ambient value in our tree */
930 AMBVAL *av
931 )
932 {
933 AMBTREE *at;
934 AMBVAL *ap;
935 AMBVAL avh;
936 FVECT ck0;
937 double s;
938 int branch;
939 int i;
940
941 if (av->rad <= FTINY)
942 error(CONSISTENCY, "zero ambient radius in avinsert");
943 at = &atrunk;
944 VCOPY(ck0, thescene.cuorg);
945 s = thescene.cusize;
946 while (s*(OCTSCALE/2) > av->rad*ambacc) {
947 if (at->kid == NULL)
948 if ((at->kid = newambtree()) == NULL)
949 error(SYSTEM, "out of memory in avinsert");
950 s *= 0.5;
951 branch = 0;
952 for (i = 0; i < 3; i++)
953 if (av->pos[i] > ck0[i] + s) {
954 ck0[i] += s;
955 branch |= 1 << i;
956 }
957 at = at->kid + branch;
958 }
959 avh.next = at->alist; /* order by increasing level */
960 for (ap = &avh; ap->next != NULL; ap = ap->next)
961 if ( ap->next->lvl > av->lvl ||
962 (ap->next->lvl == av->lvl) &
963 (ap->next->weight <= av->weight) )
964 break;
965 av->next = ap->next;
966 ap->next = (AMBVAL*)av;
967 at->alist = avh.next;
968 }
969
970 #endif /* ! NEWAMB */
971
972 /************* FOLLOWING ROUTINES SAME FOR NEW & OLD METHODS ***************/
973
974 static void
975 initambfile( /* initialize ambient file */
976 int cre8
977 )
978 {
979 extern char *progname, *octname;
980 static char *mybuf = NULL;
981
982 #ifdef F_SETLKW
983 aflock(cre8 ? F_WRLCK : F_RDLCK);
984 #endif
985 SET_FILE_BINARY(ambfp);
986 if (mybuf == NULL)
987 mybuf = (char *)bmalloc(BUFSIZ+8);
988 setbuf(ambfp, mybuf);
989 if (cre8) { /* new file */
990 newheader("RADIANCE", ambfp);
991 fprintf(ambfp, "%s -av %g %g %g -aw %d -ab %d -aa %g ",
992 progname, colval(ambval,RED),
993 colval(ambval,GRN), colval(ambval,BLU),
994 ambvwt, ambounce, ambacc);
995 fprintf(ambfp, "-ad %d -as %d -ar %d ",
996 ambdiv, ambssamp, ambres);
997 if (octname != NULL)
998 fputs(octname, ambfp);
999 fputc('\n', ambfp);
1000 fprintf(ambfp, "SOFTWARE= %s\n", VersionID);
1001 fputnow(ambfp);
1002 fputformat(AMBFMT, ambfp);
1003 fputc('\n', ambfp);
1004 putambmagic(ambfp);
1005 } else if (checkheader(ambfp, AMBFMT, NULL) < 0 || !hasambmagic(ambfp))
1006 error(USER, "bad ambient file");
1007 }
1008
1009
1010 static void
1011 avsave( /* insert and save an ambient value */
1012 AMBVAL *av
1013 )
1014 {
1015 avstore(av);
1016 if (ambfp == NULL)
1017 return;
1018 if (writambval(av, ambfp) < 0)
1019 goto writerr;
1020 if (++nunflshed >= AMBFLUSH)
1021 if (ambsync() == EOF)
1022 goto writerr;
1023 return;
1024 writerr:
1025 error(SYSTEM, "error writing to ambient file");
1026 }
1027
1028
1029 static AMBVAL *
1030 avstore( /* allocate memory and save aval */
1031 AMBVAL *aval
1032 )
1033 {
1034 AMBVAL *av;
1035 double d;
1036
1037 if ((av = newambval()) == NULL)
1038 error(SYSTEM, "out of memory in avstore");
1039 *av = *aval;
1040 av->latick = ambclock;
1041 av->next = NULL;
1042 nambvals++;
1043 d = bright(av->val);
1044 if (d > FTINY) { /* add to log sum for averaging */
1045 avsum += log(d);
1046 navsum++;
1047 }
1048 avinsert(av); /* insert in our cache tree */
1049 return(av);
1050 }
1051
1052
1053 #define ATALLOCSZ 512 /* #/8 trees to allocate at once */
1054
1055 static AMBTREE *atfreelist = NULL; /* free ambient tree structures */
1056
1057
1058 static AMBTREE *
1059 newambtree(void) /* allocate 8 ambient tree structs */
1060 {
1061 AMBTREE *atp, *upperlim;
1062
1063 if (atfreelist == NULL) { /* get more nodes */
1064 atfreelist = (AMBTREE *)malloc(ATALLOCSZ*8*sizeof(AMBTREE));
1065 if (atfreelist == NULL)
1066 return(NULL);
1067 /* link new free list */
1068 upperlim = atfreelist + 8*(ATALLOCSZ-1);
1069 for (atp = atfreelist; atp < upperlim; atp += 8)
1070 atp->kid = atp + 8;
1071 atp->kid = NULL;
1072 }
1073 atp = atfreelist;
1074 atfreelist = atp->kid;
1075 memset((char *)atp, '\0', 8*sizeof(AMBTREE));
1076 return(atp);
1077 }
1078
1079
1080 static void
1081 freeambtree( /* free 8 ambient tree structs */
1082 AMBTREE *atp
1083 )
1084 {
1085 atp->kid = atfreelist;
1086 atfreelist = atp;
1087 }
1088
1089
1090 static void
1091 unloadatree( /* unload an ambient value tree */
1092 AMBTREE *at,
1093 unloadtf_t *f
1094 )
1095 {
1096 AMBVAL *av;
1097 int i;
1098 /* transfer values at this node */
1099 for (av = at->alist; av != NULL; av = at->alist) {
1100 at->alist = av->next;
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((char *)&avaddr, (char *)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((void *)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((char *)avlist1, nambvals, sizeof(struct avl), alatcmp);
1250 qsort((char *)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((void *)avlist1);
1267 free((void *)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 */