ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/ambient.c
Revision: 2.111
Committed: Thu Nov 18 19:38:21 2021 UTC (3 years, 5 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.110: +7 -3 lines
Log Message:
perf(rpict, rtrace, rvu): improved function of -aw rendering parameter

File Contents

# User Rev Content
1 greg 2.111 static const char RCSid[] = "$Id: ambient.c,v 2.110 2021/02/19 22:05:46 greg Exp $";
2 greg 1.1 /*
3     * ambient.c - routines dealing with ambient (inter-reflected) component.
4 greg 2.47 *
5     * Declarations of external symbols in ambient.h
6     */
7    
8 greg 2.48 #include "copyright.h"
9 greg 1.1
10 schorsch 2.52 #include <string.h>
11    
12 schorsch 2.50 #include "platform.h"
13 greg 1.1 #include "ray.h"
14 greg 1.8 #include "otypes.h"
15 greg 2.105 #include "otspecial.h"
16 schorsch 2.56 #include "resolu.h"
17 greg 1.14 #include "ambient.h"
18 greg 1.1 #include "random.h"
19 greg 2.94 #include "pmapamb.h"
20 greg 1.1
21 greg 2.32 #ifndef OCTSCALE
22 greg 2.29 #define OCTSCALE 1.0 /* ceil((valid rad.)/(cube size)) */
23 greg 2.32 #endif
24 greg 1.1
25 greg 2.27 extern char *shm_boundary; /* memory sharing boundary */
26 greg 2.26
27 greg 2.57 #ifndef MAXASET
28 greg 2.67 #define MAXASET 4095 /* maximum number of elements in ambient set */
29 greg 2.57 #endif
30 greg 2.12 OBJECT ambset[MAXASET+1]={0}; /* ambient include/exclude set */
31 greg 1.1
32 greg 2.12 double maxarad; /* maximum ambient radius */
33     double minarad; /* minimum ambient radius */
34 greg 1.1
35 greg 2.12 static AMBTREE atrunk; /* our ambient trunk node */
36 greg 1.1
37     static FILE *ambfp = NULL; /* ambient file pointer */
38 greg 2.16 static int nunflshed = 0; /* number of unflushed ambient values */
39 greg 1.1
40 greg 2.26 #ifndef SORT_THRESH
41 greg 2.49 #ifdef SMLMEM
42 greg 2.71 #define SORT_THRESH ((16L<<20)/sizeof(AMBVAL))
43 greg 2.49 #else
44 greg 2.71 #define SORT_THRESH ((64L<<20)/sizeof(AMBVAL))
45 greg 2.26 #endif
46     #endif
47     #ifndef SORT_INTVL
48 greg 2.41 #define SORT_INTVL (SORT_THRESH<<1)
49 greg 2.26 #endif
50 greg 2.28 #ifndef MAX_SORT_INTVL
51 greg 2.41 #define MAX_SORT_INTVL (SORT_INTVL<<6)
52 greg 2.28 #endif
53 greg 2.26
54 greg 2.74
55 gregl 2.43 static double avsum = 0.; /* computed ambient value sum (log) */
56     static unsigned int navsum = 0; /* number of values in avsum */
57 greg 2.35 static unsigned int nambvals = 0; /* total number of indirect values */
58     static unsigned int nambshare = 0; /* number of values from file */
59 greg 2.27 static unsigned long ambclock = 0; /* ambient access clock */
60     static unsigned long lastsort = 0; /* time of last value sort */
61 greg 2.26 static long sortintvl = SORT_INTVL; /* time until next sort */
62 greg 2.47 static FILE *ambinp = NULL; /* auxiliary file for input */
63     static long lastpos = -1; /* last flush position */
64 greg 2.26
65     #define MAXACLOCK (1L<<30) /* clock turnover value */
66 greg 2.27 /*
67     * Track access times unless we are sharing ambient values
68     * through memory on a multiprocessor, when we want to avoid
69 greg 2.35 * claiming our own memory (copy on write). Go ahead anyway
70     * if more than two thirds of our values are unshared.
71 gregl 2.43 * Compile with -Dtracktime=0 to turn this code off.
72 greg 2.27 */
73 gregl 2.43 #ifndef tracktime
74 greg 2.35 #define tracktime (shm_boundary == NULL || nambvals > 3*nambshare)
75 gregl 2.43 #endif
76 greg 2.26
77 greg 2.12 #define AMBFLUSH (BUFSIZ/AMBVALSIZ)
78 greg 2.6
79 greg 2.47 #define newambval() (AMBVAL *)malloc(sizeof(AMBVAL))
80 greg 1.1
81 greg 2.108 #define tfunc(lwr, x, upr) (((x)-(lwr))/((upr)-(lwr)))
82    
83 schorsch 2.56 static void initambfile(int creat);
84     static void avsave(AMBVAL *av);
85     static AMBVAL *avstore(AMBVAL *aval);
86     static AMBTREE *newambtree(void);
87     static void freeambtree(AMBTREE *atp);
88    
89 greg 2.71 typedef void unloadtf_t(AMBVAL *);
90 schorsch 2.56 static unloadtf_t avinsert;
91     static unloadtf_t av2list;
92 greg 2.71 static unloadtf_t avfree;
93 schorsch 2.56 static void unloadatree(AMBTREE *at, unloadtf_t *f);
94    
95     static int aposcmp(const void *avp1, const void *avp2);
96     static int avlmemi(AMBVAL *avaddr);
97     static void sortambvals(int always);
98    
99 greg 2.108 static int plugaleak(RAY *r, AMBVAL *ap, FVECT anorm, double ang);
100     static double sumambient(COLOR acol, RAY *r, FVECT rn, int al,
101     AMBTREE *at, FVECT c0, double s);
102     static int makeambient(COLOR acol, RAY *r, FVECT rn, int al);
103     static int extambient(COLOR cr, AMBVAL *ap, FVECT pv, FVECT nv,
104     FVECT uvw[3]);
105    
106 greg 2.19 #ifdef F_SETLKW
107 schorsch 2.56 static void aflock(int typ);
108 greg 2.19 #endif
109 greg 1.1
110 greg 2.7
111 greg 2.71 void
112 schorsch 2.56 setambres( /* set ambient resolution */
113     int ar
114     )
115 greg 2.3 {
116 greg 2.21 ambres = ar < 0 ? 0 : ar; /* may be done already */
117 greg 2.3 /* set min & max radii */
118     if (ar <= 0) {
119 greg 2.25 minarad = 0;
120 greg 2.89 maxarad = thescene.cusize*0.2;
121 greg 2.3 } else {
122     minarad = thescene.cusize / ar;
123 greg 2.83 maxarad = 64.0 * minarad; /* heuristic */
124 greg 2.89 if (maxarad > thescene.cusize*0.2)
125     maxarad = thescene.cusize*0.2;
126 greg 2.3 }
127 greg 2.25 if (minarad <= FTINY)
128 greg 2.83 minarad = 10.0*FTINY;
129 greg 2.25 if (maxarad <= minarad)
130 greg 2.83 maxarad = 64.0 * minarad;
131 greg 2.3 }
132    
133    
134 greg 2.71 void
135 schorsch 2.56 setambacc( /* set ambient accuracy */
136     double newa
137     )
138 greg 2.20 {
139 greg 2.84 static double olda; /* remember previous setting here */
140 greg 2.82
141     newa *= (newa > 0);
142     if (fabs(newa - olda) >= .05*(newa + olda)) {
143 greg 2.84 ambacc = newa;
144 greg 2.107 if (ambacc > FTINY && nambvals > 0)
145 greg 2.74 sortambvals(1); /* rebuild tree */
146     }
147 greg 2.20 }
148    
149    
150 greg 2.71 void
151 schorsch 2.56 setambient(void) /* initialize calculation */
152 greg 1.1 {
153 gwlarson 2.46 int readonly = 0;
154 greg 2.66 long flen;
155 greg 2.12 AMBVAL amb;
156 greg 2.47 /* make sure we're fresh */
157     ambdone();
158 greg 2.3 /* init ambient limits */
159     setambres(ambres);
160 greg 2.82 setambacc(ambacc);
161 greg 2.47 if (ambfile == NULL || !ambfile[0])
162 greg 2.19 return;
163 greg 2.20 if (ambacc <= FTINY) {
164 greg 2.21 sprintf(errmsg, "zero ambient accuracy so \"%s\" not opened",
165 greg 2.47 ambfile);
166 greg 2.20 error(WARNING, errmsg);
167     return;
168     }
169 greg 2.3 /* open ambient file */
170 greg 2.47 if ((ambfp = fopen(ambfile, "r+")) == NULL)
171     readonly = (ambfp = fopen(ambfile, "r")) != NULL;
172 gwlarson 2.46 if (ambfp != NULL) {
173     initambfile(0); /* file exists */
174 greg 2.66 lastpos = ftell(ambfp);
175 greg 2.19 while (readambval(&amb, ambfp))
176 greg 2.88 avstore(&amb);
177 gwlarson 2.46 nambshare = nambvals; /* share loaded values */
178     if (readonly) {
179     sprintf(errmsg,
180     "loaded %u values from read-only ambient file",
181     nambvals);
182     error(WARNING, errmsg);
183     fclose(ambfp); /* close file so no writes */
184     ambfp = NULL;
185     return; /* avoid ambsync() */
186     }
187     /* align file pointer */
188 greg 2.66 lastpos += (long)nambvals*AMBVALSIZ;
189 greg 2.55 flen = lseek(fileno(ambfp), (off_t)0, SEEK_END);
190 greg 2.66 if (flen != lastpos) {
191 greg 2.42 sprintf(errmsg,
192 greg 2.37 "ignoring last %ld values in ambient file (corrupted)",
193 greg 2.66 (flen - lastpos)/AMBVALSIZ);
194 greg 2.42 error(WARNING, errmsg);
195 greg 2.66 fseek(ambfp, lastpos, SEEK_SET);
196     ftruncate(fileno(ambfp), (off_t)lastpos);
197 greg 2.37 }
198 greg 2.47 } else if ((ambfp = fopen(ambfile, "w+")) != NULL) {
199 gwlarson 2.46 initambfile(1); /* else create new file */
200 greg 2.68 fflush(ambfp);
201 greg 2.66 lastpos = ftell(ambfp);
202 gwlarson 2.46 } else {
203 greg 2.47 sprintf(errmsg, "cannot open ambient file \"%s\"", ambfile);
204 greg 2.19 error(SYSTEM, errmsg);
205 greg 2.16 }
206 greg 2.68 #ifdef F_SETLKW
207     aflock(F_UNLCK); /* release file */
208     #endif
209 greg 1.8 }
210    
211    
212 greg 2.71 void
213 schorsch 2.56 ambdone(void) /* close ambient file and free memory */
214 greg 2.47 {
215     if (ambfp != NULL) { /* close ambient file */
216     ambsync();
217     fclose(ambfp);
218     ambfp = NULL;
219     if (ambinp != NULL) {
220     fclose(ambinp);
221     ambinp = NULL;
222     }
223     lastpos = -1;
224     }
225     /* free ambient tree */
226 greg 2.99 unloadatree(&atrunk, avfree);
227 greg 2.47 /* reset state variables */
228     avsum = 0.;
229     navsum = 0;
230     nambvals = 0;
231     nambshare = 0;
232     ambclock = 0;
233     lastsort = 0;
234     sortintvl = SORT_INTVL;
235     }
236    
237    
238 greg 2.71 void
239 schorsch 2.56 ambnotify( /* record new modifier */
240     OBJECT obj
241     )
242 greg 1.8 {
243 greg 1.11 static int hitlimit = 0;
244 greg 2.70 OBJREC *o;
245     char **amblp;
246 greg 1.8
247 greg 2.47 if (obj == OVOID) { /* starting over */
248     ambset[0] = 0;
249     hitlimit = 0;
250     return;
251     }
252     o = objptr(obj);
253 greg 1.11 if (hitlimit || !ismodifier(o->otype))
254 greg 1.8 return;
255     for (amblp = amblist; *amblp != NULL; amblp++)
256     if (!strcmp(o->oname, *amblp)) {
257 greg 1.11 if (ambset[0] >= MAXASET) {
258     error(WARNING, "too many modifiers in ambient list");
259     hitlimit++;
260     return; /* should this be fatal? */
261     }
262 greg 1.8 insertelem(ambset, obj);
263     return;
264 greg 1.1 }
265     }
266    
267 greg 2.71
268     void
269     multambient( /* compute ambient component & multiply by coef. */
270     COLOR aval,
271     RAY *r,
272     FVECT nrm
273     )
274     {
275 greg 2.111 static double logAvgAbsorp = 1;
276 greg 2.71 static int rdepth = 0; /* ambient recursion */
277 greg 2.94 COLOR acol, caustic;
278 greg 2.103 int i, ok;
279 greg 2.71 double d, l;
280    
281 greg 2.94 /* PMAP: Factor in ambient from photon map, if enabled and ray is
282     * ambient. Return as all ambient components accounted for, else
283     * continue. */
284     if (ambPmap(aval, r, rdepth))
285     return;
286    
287 greg 2.111 if (logAvgAbsorp > 0) /* exclude in -aw to avoid growth */
288     logAvgAbsorp = log(1.-AVGREFL);
289    
290 greg 2.94 /* PMAP: Factor in specular-diffuse ambient (caustics) from photon
291     * map, if enabled and ray is primary, else caustic is zero. Continue
292     * with RADIANCE ambient calculation */
293     copycolor(caustic, aval);
294     ambPmapCaustic(caustic, r, rdepth);
295    
296 greg 2.71 if (ambdiv <= 0) /* no ambient calculation */
297     goto dumbamb;
298     /* check number of bounces */
299     if (rdepth >= ambounce)
300     goto dumbamb;
301     /* check ambient list */
302     if (ambincl != -1 && r->ro != NULL &&
303     ambincl != inset(ambset, r->ro->omod))
304     goto dumbamb;
305    
306     if (ambacc <= FTINY) { /* no ambient storage */
307 greg 2.103 FVECT uvd[2];
308     float dgrad[2], *dgp = NULL;
309    
310     if (nrm != r->ron && DOT(nrm,r->ron) < 0.9999)
311     dgp = dgrad; /* compute rotational grad. */
312 greg 2.71 copycolor(acol, aval);
313     rdepth++;
314 greg 2.86 ok = doambient(acol, r, r->rweight,
315 greg 2.103 uvd, NULL, NULL, dgp, NULL);
316 greg 2.71 rdepth--;
317     if (!ok)
318     goto dumbamb;
319 greg 2.103 if ((ok > 0) & (dgp != NULL)) { /* apply texture */
320     FVECT v1;
321     VCROSS(v1, r->ron, nrm);
322     d = 1.0;
323     for (i = 3; i--; )
324     d += v1[i] * (dgp[0]*uvd[0][i] + dgp[1]*uvd[1][i]);
325     if (d >= 0.05)
326     scalecolor(acol, d);
327     }
328 greg 2.71 copycolor(aval, acol);
329 greg 2.94
330     /* PMAP: add in caustic */
331     addcolor(aval, caustic);
332 greg 2.71 return;
333     }
334    
335     if (tracktime) /* sort to minimize thrashing */
336     sortambvals(0);
337     /* interpolate ambient value */
338     setcolor(acol, 0.0, 0.0, 0.0);
339     d = sumambient(acol, r, nrm, rdepth,
340     &atrunk, thescene.cuorg, thescene.cusize);
341 greg 2.94
342 greg 2.71 if (d > FTINY) {
343     d = 1.0/d;
344     scalecolor(acol, d);
345     multcolor(aval, acol);
346 greg 2.94
347     /* PMAP: add in caustic */
348     addcolor(aval, caustic);
349 greg 2.71 return;
350     }
351 greg 2.94
352 greg 2.71 rdepth++; /* need to cache new value */
353     ok = makeambient(acol, r, nrm, rdepth-1);
354     rdepth--;
355 greg 2.94
356 greg 2.71 if (ok) {
357 greg 2.73 multcolor(aval, acol); /* computed new value */
358 greg 2.94
359     /* PMAP: add in caustic */
360     addcolor(aval, caustic);
361 greg 2.71 return;
362     }
363 greg 2.94
364 greg 2.71 dumbamb: /* return global value */
365     if ((ambvwt <= 0) | (navsum == 0)) {
366     multcolor(aval, ambval);
367 greg 2.94
368     /* PMAP: add in caustic */
369     addcolor(aval, caustic);
370 greg 2.71 return;
371     }
372 greg 2.94
373     l = bright(ambval); /* average in computations */
374 greg 2.71 if (l > FTINY) {
375 greg 2.111 d = (log(l)*(double)ambvwt + avsum + logAvgAbsorp*navsum) /
376 greg 2.71 (double)(ambvwt + navsum);
377     d = exp(d) / l;
378     scalecolor(aval, d);
379     multcolor(aval, ambval); /* apply color of ambval */
380     } else {
381 greg 2.111 d = exp( avsum/(double)navsum + logAvgAbsorp );
382 greg 2.71 scalecolor(aval, d); /* neutral color */
383     }
384     }
385    
386    
387 greg 2.86 /* Plug a potential leak where ambient cache value is occluded */
388     static int
389     plugaleak(RAY *r, AMBVAL *ap, FVECT anorm, double ang)
390     {
391     const double cost70sq = 0.1169778; /* cos(70deg)^2 */
392     RAY rtst;
393     FVECT vdif;
394     double normdot, ndotd, nadotd;
395     double a, b, c, t[2];
396    
397     ang += 2.*PI*(ang < 0); /* check direction flags */
398     if ( !(ap->corral>>(int)(ang*(16./PI)) & 1) )
399     return(0);
400     /*
401     * Generate test ray, targeting 20 degrees above sample point plane
402     * along surface normal from cache position. This should be high
403     * enough to miss local geometry we don't really care about.
404     */
405     VSUB(vdif, ap->pos, r->rop);
406     normdot = DOT(anorm, r->ron);
407     ndotd = DOT(vdif, r->ron);
408     nadotd = DOT(vdif, anorm);
409     a = normdot*normdot - cost70sq;
410     b = 2.0*(normdot*ndotd - nadotd*cost70sq);
411     c = ndotd*ndotd - DOT(vdif,vdif)*cost70sq;
412     if (quadratic(t, a, b, c) != 2)
413     return(1); /* should rarely happen */
414     if (t[1] <= FTINY)
415     return(0); /* should fail behind test */
416     rayorigin(&rtst, SHADOW, r, NULL);
417     VSUM(rtst.rdir, vdif, anorm, t[1]); /* further dist. > plane */
418     rtst.rmax = normalize(rtst.rdir); /* short ray test */
419     while (localhit(&rtst, &thescene)) { /* check for occluder */
420 greg 2.105 OBJREC *m = findmaterial(rtst.ro);
421     if (m != NULL && !istransp(m->otype) && !isBSDFproxy(m) &&
422 greg 2.86 (rtst.clipset == NULL ||
423     !inset(rtst.clipset, rtst.ro->omod)))
424     return(1); /* plug light leak */
425     VCOPY(rtst.rorg, rtst.rop); /* skip invisible surface */
426     rtst.rmax -= rtst.rot;
427     rayclear(&rtst);
428     }
429     return(0); /* seems we're OK */
430     }
431    
432    
433     static double
434 greg 2.72 sumambient( /* get interpolated ambient value */
435 greg 2.71 COLOR acol,
436     RAY *r,
437     FVECT rn,
438     int al,
439     AMBTREE *at,
440     FVECT c0,
441     double s
442     )
443 greg 2.80 { /* initial limit is 10 degrees plus ambacc radians */
444     const double minangle = 10.0 * PI/180.;
445 greg 2.84 double maxangle = minangle + ambacc;
446 greg 2.71 double wsum = 0.0;
447     FVECT ck0;
448     int i, j;
449     AMBVAL *av;
450 greg 2.83
451     if (at->kid != NULL) { /* sum children first */
452     s *= 0.5;
453     for (i = 0; i < 8; i++) {
454     for (j = 0; j < 3; j++) {
455     ck0[j] = c0[j];
456     if (1<<j & i)
457     ck0[j] += s;
458     if (r->rop[j] < ck0[j] - OCTSCALE*s)
459     break;
460     if (r->rop[j] > ck0[j] + (1.0+OCTSCALE)*s)
461     break;
462     }
463     if (j == 3)
464     wsum += sumambient(acol, r, rn, al,
465     at->kid+i, ck0, s);
466     }
467     /* good enough? */
468 greg 2.84 if (wsum >= 0.05 && s > minarad*10.0)
469 greg 2.83 return(wsum);
470     }
471 greg 2.84 /* adjust maximum angle */
472     if (at->alist != NULL && (at->alist->lvl <= al) & (r->rweight < 0.6))
473     maxangle = (maxangle - PI/2.)*pow(r->rweight,0.13) + PI/2.;
474 greg 2.71 /* sum this node */
475     for (av = at->alist; av != NULL; av = av->next) {
476 greg 2.86 double u, v, d, delta_r2, delta_t2;
477 greg 2.71 COLOR ct;
478     FVECT uvw[3];
479     /* record access */
480     if (tracktime)
481     av->latick = ambclock;
482     /*
483 greg 2.72 * Ambient level test
484 greg 2.71 */
485 greg 2.87 if (av->lvl > al || /* list sorted, so this works */
486     (av->lvl == al) & (av->weight < 0.9*r->rweight))
487 greg 2.71 break;
488     /*
489 greg 2.72 * Direction test using unperturbed normal
490 greg 2.71 */
491     decodedir(uvw[2], av->ndir);
492     d = DOT(uvw[2], r->ron);
493     if (d <= 0.0) /* >= 90 degrees */
494     continue;
495     delta_r2 = 2.0 - 2.0*d; /* approx. radians^2 */
496     if (delta_r2 >= maxangle*maxangle)
497     continue;
498     /*
499 greg 2.81 * Modified ray behind test
500     */
501 greg 2.85 VSUB(ck0, r->rop, av->pos);
502 greg 2.81 d = DOT(ck0, uvw[2]);
503 greg 2.84 if (d < -minarad*ambacc-.001)
504 greg 2.81 continue;
505     d /= av->rad[0];
506     delta_t2 = d*d;
507 greg 2.84 if (delta_t2 >= ambacc*ambacc)
508 greg 2.81 continue;
509     /*
510 greg 2.72 * Elliptical radii test based on Hessian
511 greg 2.71 */
512     decodedir(uvw[0], av->udir);
513     VCROSS(uvw[1], uvw[2], uvw[0]);
514 greg 2.86 d = (u = DOT(ck0, uvw[0])) / av->rad[0];
515 greg 2.81 delta_t2 += d*d;
516 greg 2.86 d = (v = DOT(ck0, uvw[1])) / av->rad[1];
517 greg 2.71 delta_t2 += d*d;
518 greg 2.84 if (delta_t2 >= ambacc*ambacc)
519 greg 2.71 continue;
520     /*
521 greg 2.86 * Test for potential light leak
522     */
523     if (av->corral && plugaleak(r, av, uvw[2], atan2a(v,u)))
524     continue;
525     /*
526 greg 2.72 * Extrapolate value and compute final weight (hat function)
527 greg 2.71 */
528 greg 2.92 if (!extambient(ct, av, r->rop, rn, uvw))
529     continue;
530 greg 2.71 d = tfunc(maxangle, sqrt(delta_r2), 0.0) *
531 greg 2.84 tfunc(ambacc, sqrt(delta_t2), 0.0);
532 greg 2.71 scalecolor(ct, d);
533     addcolor(acol, ct);
534 greg 2.72 wsum += d;
535 greg 2.71 }
536     return(wsum);
537     }
538    
539    
540 greg 2.86 static int
541 greg 2.71 makeambient( /* make a new ambient value for storage */
542     COLOR acol,
543     RAY *r,
544     FVECT rn,
545     int al
546     )
547     {
548     AMBVAL amb;
549 greg 2.72 FVECT uvw[3];
550 greg 2.71 int i;
551    
552     amb.weight = 1.0; /* compute weight */
553     for (i = al; i-- > 0; )
554     amb.weight *= AVGREFL;
555     if (r->rweight < 0.1*amb.weight) /* heuristic override */
556     amb.weight = 1.25*r->rweight;
557     setcolor(acol, AVGREFL, AVGREFL, AVGREFL);
558     /* compute ambient */
559 greg 2.86 i = doambient(acol, r, amb.weight,
560     uvw, amb.rad, amb.gpos, amb.gdir, &amb.corral);
561 greg 2.71 scalecolor(acol, 1./AVGREFL); /* undo assumed reflectance */
562 greg 2.76 if (i <= 0 || amb.rad[0] <= FTINY) /* no Hessian or zero radius */
563 greg 2.73 return(i);
564 greg 2.71 /* store value */
565     VCOPY(amb.pos, r->rop);
566     amb.ndir = encodedir(r->ron);
567 greg 2.72 amb.udir = encodedir(uvw[0]);
568 greg 2.71 amb.lvl = al;
569     copycolor(amb.val, acol);
570     /* insert into tree */
571     avsave(&amb); /* and save to file */
572 greg 2.72 if (rn != r->ron) { /* texture */
573     VCOPY(uvw[2], r->ron);
574     extambient(acol, &amb, r->rop, rn, uvw);
575     }
576 greg 2.71 return(1);
577     }
578    
579    
580 greg 2.92 static int
581 greg 2.71 extambient( /* extrapolate value at pv, nv */
582     COLOR cr,
583     AMBVAL *ap,
584     FVECT pv,
585 greg 2.72 FVECT nv,
586     FVECT uvw[3]
587 greg 2.71 )
588     {
589 greg 2.93 const double min_d = 0.05;
590 greg 2.110 const double max_d = 20.;
591 greg 2.72 static FVECT my_uvw[3];
592     FVECT v1;
593     int i;
594     double d = 1.0; /* zeroeth order */
595    
596     if (uvw == NULL) { /* need local coordinates? */
597     decodedir(my_uvw[2], ap->ndir);
598     decodedir(my_uvw[0], ap->udir);
599     VCROSS(my_uvw[1], my_uvw[2], my_uvw[0]);
600     uvw = my_uvw;
601     }
602 greg 2.71 for (i = 3; i--; ) /* gradient due to translation */
603     d += (pv[i] - ap->pos[i]) *
604     (ap->gpos[0]*uvw[0][i] + ap->gpos[1]*uvw[1][i]);
605    
606     VCROSS(v1, uvw[2], nv); /* gradient due to rotation */
607     for (i = 3; i--; )
608     d += v1[i] * (ap->gdir[0]*uvw[0][i] + ap->gdir[1]*uvw[1][i]);
609    
610 greg 2.110 if (d < min_d) /* clamp min/max scaling */
611 greg 2.93 d = min_d;
612 greg 2.110 else if (d > max_d)
613     d = max_d;
614 greg 2.71 copycolor(cr, ap->val);
615     scalecolor(cr, d);
616 greg 2.93 return(d > min_d);
617 greg 2.71 }
618    
619    
620     static void
621     avinsert( /* insert ambient value in our tree */
622     AMBVAL *av
623     )
624     {
625     AMBTREE *at;
626     AMBVAL *ap;
627     AMBVAL avh;
628     FVECT ck0;
629     double s;
630     int branch;
631     int i;
632    
633     if (av->rad[0] <= FTINY)
634     error(CONSISTENCY, "zero ambient radius in avinsert");
635     at = &atrunk;
636     VCOPY(ck0, thescene.cuorg);
637     s = thescene.cusize;
638 greg 2.84 while (s*(OCTSCALE/2) > av->rad[1]*ambacc) {
639 greg 2.71 if (at->kid == NULL)
640     if ((at->kid = newambtree()) == NULL)
641     error(SYSTEM, "out of memory in avinsert");
642     s *= 0.5;
643     branch = 0;
644     for (i = 0; i < 3; i++)
645     if (av->pos[i] > ck0[i] + s) {
646     ck0[i] += s;
647     branch |= 1 << i;
648     }
649     at = at->kid + branch;
650     }
651     avh.next = at->alist; /* order by increasing level */
652     for (ap = &avh; ap->next != NULL; ap = ap->next)
653 greg 2.87 if ( ap->next->lvl > av->lvl ||
654     (ap->next->lvl == av->lvl) &
655     (ap->next->weight <= av->weight) )
656 greg 2.71 break;
657     av->next = ap->next;
658     ap->next = (AMBVAL*)av;
659     at->alist = avh.next;
660     }
661    
662    
663     static void
664 schorsch 2.56 initambfile( /* initialize ambient file */
665 greg 2.68 int cre8
666 schorsch 2.56 )
667 greg 2.9 {
668 greg 2.47 extern char *progname, *octname;
669     static char *mybuf = NULL;
670 greg 2.9
671 greg 2.19 #ifdef F_SETLKW
672 greg 2.69 aflock(cre8 ? F_WRLCK : F_RDLCK);
673 greg 2.19 #endif
674 schorsch 2.50 SET_FILE_BINARY(ambfp);
675 greg 2.47 if (mybuf == NULL)
676     mybuf = (char *)bmalloc(BUFSIZ+8);
677     setbuf(ambfp, mybuf);
678 greg 2.68 if (cre8) { /* new file */
679 greg 2.24 newheader("RADIANCE", ambfp);
680 greg 2.41 fprintf(ambfp, "%s -av %g %g %g -aw %d -ab %d -aa %g ",
681 greg 2.9 progname, colval(ambval,RED),
682     colval(ambval,GRN), colval(ambval,BLU),
683 greg 2.41 ambvwt, ambounce, ambacc);
684 greg 2.47 fprintf(ambfp, "-ad %d -as %d -ar %d ",
685     ambdiv, ambssamp, ambres);
686     if (octname != NULL)
687 greg 2.68 fputs(octname, ambfp);
688     fputc('\n', ambfp);
689 greg 2.9 fprintf(ambfp, "SOFTWARE= %s\n", VersionID);
690 greg 2.47 fputnow(ambfp);
691 greg 2.9 fputformat(AMBFMT, ambfp);
692 greg 2.68 fputc('\n', ambfp);
693 greg 2.9 putambmagic(ambfp);
694 greg 2.17 } else if (checkheader(ambfp, AMBFMT, NULL) < 0 || !hasambmagic(ambfp))
695     error(USER, "bad ambient file");
696 greg 2.9 }
697    
698    
699 greg 2.47 static void
700 schorsch 2.56 avsave( /* insert and save an ambient value */
701     AMBVAL *av
702     )
703 greg 1.1 {
704 greg 2.88 avstore(av);
705 greg 1.1 if (ambfp == NULL)
706     return;
707 greg 2.5 if (writambval(av, ambfp) < 0)
708 greg 1.1 goto writerr;
709 greg 2.16 if (++nunflshed >= AMBFLUSH)
710 greg 2.7 if (ambsync() == EOF)
711 greg 1.1 goto writerr;
712     return;
713     writerr:
714 greg 2.47 error(SYSTEM, "error writing to ambient file");
715 greg 1.1 }
716    
717    
718 greg 2.20 static AMBVAL *
719 greg 2.88 avstore( /* allocate memory and save aval */
720 greg 2.70 AMBVAL *aval
721 schorsch 2.56 )
722 greg 2.20 {
723 greg 2.70 AMBVAL *av;
724 gregl 2.43 double d;
725 greg 2.20
726     if ((av = newambval()) == NULL)
727     error(SYSTEM, "out of memory in avstore");
728 schorsch 2.53 *av = *aval;
729 greg 2.26 av->latick = ambclock;
730     av->next = NULL;
731     nambvals++;
732 gregl 2.43 d = bright(av->val);
733     if (d > FTINY) { /* add to log sum for averaging */
734     avsum += log(d);
735     navsum++;
736     }
737 greg 2.88 avinsert(av); /* insert in our cache tree */
738 greg 2.20 return(av);
739     }
740    
741    
742 greg 2.26 #define ATALLOCSZ 512 /* #/8 trees to allocate at once */
743    
744     static AMBTREE *atfreelist = NULL; /* free ambient tree structures */
745    
746    
747 greg 2.47 static AMBTREE *
748 schorsch 2.56 newambtree(void) /* allocate 8 ambient tree structs */
749 greg 2.26 {
750 greg 2.70 AMBTREE *atp, *upperlim;
751 greg 2.26
752     if (atfreelist == NULL) { /* get more nodes */
753 greg 2.47 atfreelist = (AMBTREE *)malloc(ATALLOCSZ*8*sizeof(AMBTREE));
754 greg 2.26 if (atfreelist == NULL)
755     return(NULL);
756     /* link new free list */
757     upperlim = atfreelist + 8*(ATALLOCSZ-1);
758     for (atp = atfreelist; atp < upperlim; atp += 8)
759     atp->kid = atp + 8;
760     atp->kid = NULL;
761     }
762     atp = atfreelist;
763     atfreelist = atp->kid;
764 greg 2.102 memset(atp, 0, 8*sizeof(AMBTREE));
765 greg 2.26 return(atp);
766     }
767    
768    
769 greg 2.47 static void
770 schorsch 2.56 freeambtree( /* free 8 ambient tree structs */
771     AMBTREE *atp
772     )
773 greg 2.26 {
774 greg 2.100 atp->kid = atfreelist;
775 greg 2.26 atfreelist = atp;
776     }
777    
778    
779 greg 2.47 static void
780 schorsch 2.56 unloadatree( /* unload an ambient value tree */
781 greg 2.70 AMBTREE *at,
782 schorsch 2.56 unloadtf_t *f
783     )
784 greg 2.20 {
785 greg 2.70 AMBVAL *av;
786     int i;
787 greg 2.20 /* transfer values at this node */
788     for (av = at->alist; av != NULL; av = at->alist) {
789     at->alist = av->next;
790 greg 2.102 av->next = NULL;
791 greg 2.26 (*f)(av);
792 greg 2.20 }
793 greg 2.21 if (at->kid == NULL)
794     return;
795 greg 2.20 for (i = 0; i < 8; i++) /* transfer and free children */
796 greg 2.26 unloadatree(at->kid+i, f);
797 greg 2.20 freeambtree(at->kid);
798 greg 2.26 at->kid = NULL;
799     }
800    
801    
802 greg 2.39 static struct avl {
803     AMBVAL *p;
804     unsigned long t;
805     } *avlist1; /* ambient value list with ticks */
806     static AMBVAL **avlist2; /* memory positions for sorting */
807 greg 2.26 static int i_avlist; /* index for lists */
808    
809 schorsch 2.56 static int alatcmp(const void *av1, const void *av2);
810 greg 2.26
811 schorsch 2.56 static void
812 greg 2.71 avfree(AMBVAL *av)
813     {
814     free(av);
815     }
816    
817     static void
818 schorsch 2.56 av2list(
819 greg 2.71 AMBVAL *av
820 schorsch 2.56 )
821 greg 2.26 {
822 greg 2.27 #ifdef DEBUG
823 greg 2.26 if (i_avlist >= nambvals)
824     error(CONSISTENCY, "too many ambient values in av2list1");
825 greg 2.27 #endif
826 schorsch 2.56 avlist1[i_avlist].p = avlist2[i_avlist] = (AMBVAL*)av;
827 greg 2.71 avlist1[i_avlist++].t = av->latick;
828 greg 2.26 }
829    
830    
831     static int
832 schorsch 2.56 alatcmp( /* compare ambient values for MRA */
833     const void *av1,
834     const void *av2
835     )
836 greg 2.26 {
837 greg 2.70 long lc = ((struct avl *)av2)->t - ((struct avl *)av1)->t;
838 greg 2.38 return(lc<0 ? -1 : lc>0 ? 1 : 0);
839 greg 2.26 }
840    
841    
842 greg 2.47 /* GW NOTE 2002/10/3:
843     * I used to compare AMBVAL pointers, but found that this was the
844     * cause of a serious consistency error with gcc, since the optimizer
845     * uses some dangerous trick in pointer subtraction that
846     * assumes pointers differ by exact struct size increments.
847     */
848 greg 2.26 static int
849 schorsch 2.56 aposcmp( /* compare ambient value positions */
850     const void *avp1,
851     const void *avp2
852     )
853 greg 2.26 {
854 greg 2.70 long diff = *(char * const *)avp1 - *(char * const *)avp2;
855 greg 2.47 if (diff < 0)
856     return(-1);
857     return(diff > 0);
858 greg 2.26 }
859    
860 greg 2.71
861 greg 2.27 static int
862 schorsch 2.56 avlmemi( /* find list position from address */
863     AMBVAL *avaddr
864     )
865 greg 2.27 {
866 greg 2.70 AMBVAL **avlpp;
867 greg 2.27
868 greg 2.102 avlpp = (AMBVAL **)bsearch(&avaddr, avlist2,
869     nambvals, sizeof(AMBVAL *), aposcmp);
870 greg 2.27 if (avlpp == NULL)
871     error(CONSISTENCY, "address not found in avlmemi");
872     return(avlpp - avlist2);
873     }
874    
875    
876 greg 2.47 static void
877 schorsch 2.56 sortambvals( /* resort ambient values */
878     int always
879     )
880 greg 2.26 {
881     AMBTREE oldatrunk;
882     AMBVAL tav, *tap, *pnext;
883 greg 2.70 int i, j;
884 greg 2.28 /* see if it's time yet */
885 greg 2.40 if (!always && (ambclock++ < lastsort+sortintvl ||
886 greg 2.28 nambvals < SORT_THRESH))
887     return;
888 greg 2.26 /*
889     * The idea here is to minimize memory thrashing
890     * in VM systems by improving reference locality.
891     * We do this by periodically sorting our stored ambient
892     * values in memory in order of most recently to least
893     * recently accessed. This ordering was chosen so that new
894     * ambient values (which tend to be less important) go into
895     * higher memory with the infrequently accessed values.
896     * Since we expect our values to need sorting less
897     * frequently as the process continues, we double our
898     * waiting interval after each call.
899     * This routine is also called by setambacc() with
900     * the "always" parameter set to 1 so that the ambient
901     * tree will be rebuilt with the new accuracy parameter.
902     */
903 greg 2.27 if (tracktime) { /* allocate pointer arrays to sort */
904 greg 2.26 avlist2 = (AMBVAL **)malloc(nambvals*sizeof(AMBVAL *));
905 greg 2.39 avlist1 = (struct avl *)malloc(nambvals*sizeof(struct avl));
906     } else {
907     avlist2 = NULL;
908     avlist1 = NULL;
909     }
910     if (avlist1 == NULL) { /* no time tracking -- rebuild tree? */
911     if (avlist2 != NULL)
912 greg 2.102 free(avlist2);
913 greg 2.27 if (always) { /* rebuild without sorting */
914 schorsch 2.53 oldatrunk = atrunk;
915 greg 2.27 atrunk.alist = NULL;
916     atrunk.kid = NULL;
917 greg 2.99 unloadatree(&oldatrunk, avinsert);
918 greg 2.27 }
919 greg 2.26 } else { /* sort memory by last access time */
920     /*
921     * Sorting memory is tricky because it isn't contiguous.
922     * We have to sort an array of pointers by MRA and also
923     * by memory position. We then copy values in "loops"
924     * to minimize memory hits. Nevertheless, we will visit
925 greg 2.27 * everyone at least twice, and this is an expensive process
926 greg 2.26 * when we're thrashing, which is when we need to do it.
927     */
928 greg 2.27 #ifdef DEBUG
929 greg 2.29 sprintf(errmsg, "sorting %u ambient values at ambclock=%lu...",
930     nambvals, ambclock);
931 greg 2.27 eputs(errmsg);
932     #endif
933     i_avlist = 0;
934 greg 2.99 unloadatree(&atrunk, av2list); /* empty current tree */
935 greg 2.27 #ifdef DEBUG
936     if (i_avlist < nambvals)
937     error(CONSISTENCY, "missing ambient values in sortambvals");
938     #endif
939 greg 2.102 qsort(avlist1, nambvals, sizeof(struct avl), alatcmp);
940     qsort(avlist2, nambvals, sizeof(AMBVAL *), aposcmp);
941 greg 2.26 for (i = 0; i < nambvals; i++) {
942 greg 2.39 if (avlist1[i].p == NULL)
943 greg 2.26 continue;
944     tap = avlist2[i];
945 schorsch 2.53 tav = *tap;
946 greg 2.39 for (j = i; (pnext = avlist1[j].p) != tap;
947 greg 2.27 j = avlmemi(pnext)) {
948 schorsch 2.53 *(avlist2[j]) = *pnext;
949 greg 2.26 avinsert(avlist2[j]);
950 greg 2.39 avlist1[j].p = NULL;
951 greg 2.26 }
952 schorsch 2.53 *(avlist2[j]) = tav;
953 greg 2.26 avinsert(avlist2[j]);
954 greg 2.39 avlist1[j].p = NULL;
955 greg 2.26 }
956 greg 2.102 free(avlist1);
957     free(avlist2);
958 greg 2.28 /* compute new sort interval */
959     sortintvl = ambclock - lastsort;
960 greg 2.32 if (sortintvl >= MAX_SORT_INTVL/2)
961 greg 2.28 sortintvl = MAX_SORT_INTVL;
962     else
963 greg 2.26 sortintvl <<= 1; /* wait twice as long next */
964 greg 2.27 #ifdef DEBUG
965     eputs("done\n");
966     #endif
967 greg 2.26 }
968     if (ambclock >= MAXACLOCK)
969     ambclock = MAXACLOCK/2;
970     lastsort = ambclock;
971 greg 2.20 }
972    
973    
974 greg 2.18 #ifdef F_SETLKW
975 greg 2.10
976 greg 2.47 static void
977 schorsch 2.56 aflock( /* lock/unlock ambient file */
978     int typ
979     )
980 greg 2.19 {
981     static struct flock fls; /* static so initialized to zeroes */
982    
983 greg 2.66 if (typ == fls.l_type) /* already called? */
984     return;
985 greg 2.109
986 greg 2.19 fls.l_type = typ;
987 greg 2.109 do
988     if (fcntl(fileno(ambfp), F_SETLKW, &fls) != -1)
989     return;
990     while (errno == EINTR);
991    
992     error(SYSTEM, "cannot (un)lock ambient file");
993 greg 2.19 }
994    
995    
996 greg 2.71 int
997 schorsch 2.56 ambsync(void) /* synchronize ambient file */
998 greg 2.7 {
999 greg 2.15 long flen;
1000 greg 2.12 AMBVAL avs;
1001 greg 2.70 int n;
1002 greg 2.16
1003 greg 2.65 if (ambfp == NULL) /* no ambient file? */
1004     return(0);
1005 greg 2.63 /* gain appropriate access */
1006     aflock(nunflshed ? F_WRLCK : F_RDLCK);
1007 greg 2.7 /* see if file has grown */
1008 greg 2.55 if ((flen = lseek(fileno(ambfp), (off_t)0, SEEK_END)) < 0)
1009 greg 2.21 goto seekerr;
1010 greg 2.68 if ((n = flen - lastpos) > 0) { /* file has grown */
1011 greg 2.104 if (ambinp == NULL) { /* get new file pointer */
1012     ambinp = fopen(ambfile, "rb");
1013 greg 2.14 if (ambinp == NULL)
1014 greg 2.104 error(SYSTEM, "fopen failed in ambsync");
1015 greg 2.14 }
1016 greg 2.66 if (fseek(ambinp, lastpos, SEEK_SET) < 0)
1017 greg 2.21 goto seekerr;
1018 greg 2.15 while (n >= AMBVALSIZ) { /* load contributed values */
1019 greg 2.37 if (!readambval(&avs, ambinp)) {
1020     sprintf(errmsg,
1021 greg 2.47 "ambient file \"%s\" corrupted near character %ld",
1022     ambfile, flen - n);
1023 greg 2.37 error(WARNING, errmsg);
1024     break;
1025     }
1026 greg 2.88 avstore(&avs);
1027 greg 2.15 n -= AMBVALSIZ;
1028     }
1029 greg 2.104 lastpos = flen - n; /* check alignment */
1030     if (n && lseek(fileno(ambfp), (off_t)lastpos, SEEK_SET) < 0)
1031     goto seekerr;
1032 greg 2.7 }
1033     n = fflush(ambfp); /* calls write() at last */
1034 greg 2.104 lastpos += (long)nunflshed*AMBVALSIZ;
1035 greg 2.19 aflock(F_UNLCK); /* release file */
1036 greg 2.16 nunflshed = 0;
1037 greg 2.7 return(n);
1038 greg 2.21 seekerr:
1039     error(SYSTEM, "seek failed in ambsync");
1040 greg 2.104 return(EOF); /* pro forma return */
1041 greg 2.18 }
1042    
1043 greg 2.71 #else /* ! F_SETLKW */
1044 greg 2.18
1045 greg 2.71 int
1046 schorsch 2.56 ambsync(void) /* flush ambient file */
1047 greg 2.18 {
1048 greg 2.65 if (ambfp == NULL)
1049     return(0);
1050 greg 2.18 nunflshed = 0;
1051     return(fflush(ambfp));
1052 greg 1.1 }
1053 greg 2.10
1054 greg 2.71 #endif /* ! F_SETLKW */