ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/ambient.c
Revision: 2.108
Committed: Tue May 14 17:39:10 2019 UTC (5 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.107: +10 -343 lines
Log Message:
Stripped out code related to old (pre-Hessian) ambient calculation

File Contents

# User Rev Content
1 greg 2.108 static const char RCSid[] = "$Id: ambient.c,v 2.107 2019/02/16 00:09:53 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     static int rdepth = 0; /* ambient recursion */
276 greg 2.94 COLOR acol, caustic;
277 greg 2.103 int i, ok;
278 greg 2.71 double d, l;
279    
280 greg 2.94 /* PMAP: Factor in ambient from photon map, if enabled and ray is
281     * ambient. Return as all ambient components accounted for, else
282     * continue. */
283     if (ambPmap(aval, r, rdepth))
284     return;
285    
286     /* PMAP: Factor in specular-diffuse ambient (caustics) from photon
287     * map, if enabled and ray is primary, else caustic is zero. Continue
288     * with RADIANCE ambient calculation */
289     copycolor(caustic, aval);
290     ambPmapCaustic(caustic, r, rdepth);
291    
292 greg 2.71 if (ambdiv <= 0) /* no ambient calculation */
293     goto dumbamb;
294     /* check number of bounces */
295     if (rdepth >= ambounce)
296     goto dumbamb;
297     /* check ambient list */
298     if (ambincl != -1 && r->ro != NULL &&
299     ambincl != inset(ambset, r->ro->omod))
300     goto dumbamb;
301    
302     if (ambacc <= FTINY) { /* no ambient storage */
303 greg 2.103 FVECT uvd[2];
304     float dgrad[2], *dgp = NULL;
305    
306     if (nrm != r->ron && DOT(nrm,r->ron) < 0.9999)
307     dgp = dgrad; /* compute rotational grad. */
308 greg 2.71 copycolor(acol, aval);
309     rdepth++;
310 greg 2.86 ok = doambient(acol, r, r->rweight,
311 greg 2.103 uvd, NULL, NULL, dgp, NULL);
312 greg 2.71 rdepth--;
313     if (!ok)
314     goto dumbamb;
315 greg 2.103 if ((ok > 0) & (dgp != NULL)) { /* apply texture */
316     FVECT v1;
317     VCROSS(v1, r->ron, nrm);
318     d = 1.0;
319     for (i = 3; i--; )
320     d += v1[i] * (dgp[0]*uvd[0][i] + dgp[1]*uvd[1][i]);
321     if (d >= 0.05)
322     scalecolor(acol, d);
323     }
324 greg 2.71 copycolor(aval, acol);
325 greg 2.94
326     /* PMAP: add in caustic */
327     addcolor(aval, caustic);
328 greg 2.71 return;
329     }
330    
331     if (tracktime) /* sort to minimize thrashing */
332     sortambvals(0);
333     /* interpolate ambient value */
334     setcolor(acol, 0.0, 0.0, 0.0);
335     d = sumambient(acol, r, nrm, rdepth,
336     &atrunk, thescene.cuorg, thescene.cusize);
337 greg 2.94
338 greg 2.71 if (d > FTINY) {
339     d = 1.0/d;
340     scalecolor(acol, d);
341     multcolor(aval, acol);
342 greg 2.94
343     /* PMAP: add in caustic */
344     addcolor(aval, caustic);
345 greg 2.71 return;
346     }
347 greg 2.94
348 greg 2.71 rdepth++; /* need to cache new value */
349     ok = makeambient(acol, r, nrm, rdepth-1);
350     rdepth--;
351 greg 2.94
352 greg 2.71 if (ok) {
353 greg 2.73 multcolor(aval, acol); /* computed new value */
354 greg 2.94
355     /* PMAP: add in caustic */
356     addcolor(aval, caustic);
357 greg 2.71 return;
358     }
359 greg 2.94
360 greg 2.71 dumbamb: /* return global value */
361     if ((ambvwt <= 0) | (navsum == 0)) {
362     multcolor(aval, ambval);
363 greg 2.94
364     /* PMAP: add in caustic */
365     addcolor(aval, caustic);
366 greg 2.71 return;
367     }
368 greg 2.94
369     l = bright(ambval); /* average in computations */
370 greg 2.71 if (l > FTINY) {
371     d = (log(l)*(double)ambvwt + avsum) /
372     (double)(ambvwt + navsum);
373     d = exp(d) / l;
374     scalecolor(aval, d);
375     multcolor(aval, ambval); /* apply color of ambval */
376     } else {
377     d = exp( avsum / (double)navsum );
378     scalecolor(aval, d); /* neutral color */
379     }
380     }
381    
382    
383 greg 2.86 /* Plug a potential leak where ambient cache value is occluded */
384     static int
385     plugaleak(RAY *r, AMBVAL *ap, FVECT anorm, double ang)
386     {
387     const double cost70sq = 0.1169778; /* cos(70deg)^2 */
388     RAY rtst;
389     FVECT vdif;
390     double normdot, ndotd, nadotd;
391     double a, b, c, t[2];
392    
393     ang += 2.*PI*(ang < 0); /* check direction flags */
394     if ( !(ap->corral>>(int)(ang*(16./PI)) & 1) )
395     return(0);
396     /*
397     * Generate test ray, targeting 20 degrees above sample point plane
398     * along surface normal from cache position. This should be high
399     * enough to miss local geometry we don't really care about.
400     */
401     VSUB(vdif, ap->pos, r->rop);
402     normdot = DOT(anorm, r->ron);
403     ndotd = DOT(vdif, r->ron);
404     nadotd = DOT(vdif, anorm);
405     a = normdot*normdot - cost70sq;
406     b = 2.0*(normdot*ndotd - nadotd*cost70sq);
407     c = ndotd*ndotd - DOT(vdif,vdif)*cost70sq;
408     if (quadratic(t, a, b, c) != 2)
409     return(1); /* should rarely happen */
410     if (t[1] <= FTINY)
411     return(0); /* should fail behind test */
412     rayorigin(&rtst, SHADOW, r, NULL);
413     VSUM(rtst.rdir, vdif, anorm, t[1]); /* further dist. > plane */
414     rtst.rmax = normalize(rtst.rdir); /* short ray test */
415     while (localhit(&rtst, &thescene)) { /* check for occluder */
416 greg 2.105 OBJREC *m = findmaterial(rtst.ro);
417     if (m != NULL && !istransp(m->otype) && !isBSDFproxy(m) &&
418 greg 2.86 (rtst.clipset == NULL ||
419     !inset(rtst.clipset, rtst.ro->omod)))
420     return(1); /* plug light leak */
421     VCOPY(rtst.rorg, rtst.rop); /* skip invisible surface */
422     rtst.rmax -= rtst.rot;
423     rayclear(&rtst);
424     }
425     return(0); /* seems we're OK */
426     }
427    
428    
429     static double
430 greg 2.72 sumambient( /* get interpolated ambient value */
431 greg 2.71 COLOR acol,
432     RAY *r,
433     FVECT rn,
434     int al,
435     AMBTREE *at,
436     FVECT c0,
437     double s
438     )
439 greg 2.80 { /* initial limit is 10 degrees plus ambacc radians */
440     const double minangle = 10.0 * PI/180.;
441 greg 2.84 double maxangle = minangle + ambacc;
442 greg 2.71 double wsum = 0.0;
443     FVECT ck0;
444     int i, j;
445     AMBVAL *av;
446 greg 2.83
447     if (at->kid != NULL) { /* sum children first */
448     s *= 0.5;
449     for (i = 0; i < 8; i++) {
450     for (j = 0; j < 3; j++) {
451     ck0[j] = c0[j];
452     if (1<<j & i)
453     ck0[j] += s;
454     if (r->rop[j] < ck0[j] - OCTSCALE*s)
455     break;
456     if (r->rop[j] > ck0[j] + (1.0+OCTSCALE)*s)
457     break;
458     }
459     if (j == 3)
460     wsum += sumambient(acol, r, rn, al,
461     at->kid+i, ck0, s);
462     }
463     /* good enough? */
464 greg 2.84 if (wsum >= 0.05 && s > minarad*10.0)
465 greg 2.83 return(wsum);
466     }
467 greg 2.84 /* adjust maximum angle */
468     if (at->alist != NULL && (at->alist->lvl <= al) & (r->rweight < 0.6))
469     maxangle = (maxangle - PI/2.)*pow(r->rweight,0.13) + PI/2.;
470 greg 2.71 /* sum this node */
471     for (av = at->alist; av != NULL; av = av->next) {
472 greg 2.86 double u, v, d, delta_r2, delta_t2;
473 greg 2.71 COLOR ct;
474     FVECT uvw[3];
475     /* record access */
476     if (tracktime)
477     av->latick = ambclock;
478     /*
479 greg 2.72 * Ambient level test
480 greg 2.71 */
481 greg 2.87 if (av->lvl > al || /* list sorted, so this works */
482     (av->lvl == al) & (av->weight < 0.9*r->rweight))
483 greg 2.71 break;
484     /*
485 greg 2.72 * Direction test using unperturbed normal
486 greg 2.71 */
487     decodedir(uvw[2], av->ndir);
488     d = DOT(uvw[2], r->ron);
489     if (d <= 0.0) /* >= 90 degrees */
490     continue;
491     delta_r2 = 2.0 - 2.0*d; /* approx. radians^2 */
492     if (delta_r2 >= maxangle*maxangle)
493     continue;
494     /*
495 greg 2.81 * Modified ray behind test
496     */
497 greg 2.85 VSUB(ck0, r->rop, av->pos);
498 greg 2.81 d = DOT(ck0, uvw[2]);
499 greg 2.84 if (d < -minarad*ambacc-.001)
500 greg 2.81 continue;
501     d /= av->rad[0];
502     delta_t2 = d*d;
503 greg 2.84 if (delta_t2 >= ambacc*ambacc)
504 greg 2.81 continue;
505     /*
506 greg 2.72 * Elliptical radii test based on Hessian
507 greg 2.71 */
508     decodedir(uvw[0], av->udir);
509     VCROSS(uvw[1], uvw[2], uvw[0]);
510 greg 2.86 d = (u = DOT(ck0, uvw[0])) / av->rad[0];
511 greg 2.81 delta_t2 += d*d;
512 greg 2.86 d = (v = DOT(ck0, uvw[1])) / av->rad[1];
513 greg 2.71 delta_t2 += d*d;
514 greg 2.84 if (delta_t2 >= ambacc*ambacc)
515 greg 2.71 continue;
516     /*
517 greg 2.86 * Test for potential light leak
518     */
519     if (av->corral && plugaleak(r, av, uvw[2], atan2a(v,u)))
520     continue;
521     /*
522 greg 2.72 * Extrapolate value and compute final weight (hat function)
523 greg 2.71 */
524 greg 2.92 if (!extambient(ct, av, r->rop, rn, uvw))
525     continue;
526 greg 2.71 d = tfunc(maxangle, sqrt(delta_r2), 0.0) *
527 greg 2.84 tfunc(ambacc, sqrt(delta_t2), 0.0);
528 greg 2.71 scalecolor(ct, d);
529     addcolor(acol, ct);
530 greg 2.72 wsum += d;
531 greg 2.71 }
532     return(wsum);
533     }
534    
535    
536 greg 2.86 static int
537 greg 2.71 makeambient( /* make a new ambient value for storage */
538     COLOR acol,
539     RAY *r,
540     FVECT rn,
541     int al
542     )
543     {
544     AMBVAL amb;
545 greg 2.72 FVECT uvw[3];
546 greg 2.71 int i;
547    
548     amb.weight = 1.0; /* compute weight */
549     for (i = al; i-- > 0; )
550     amb.weight *= AVGREFL;
551     if (r->rweight < 0.1*amb.weight) /* heuristic override */
552     amb.weight = 1.25*r->rweight;
553     setcolor(acol, AVGREFL, AVGREFL, AVGREFL);
554     /* compute ambient */
555 greg 2.86 i = doambient(acol, r, amb.weight,
556     uvw, amb.rad, amb.gpos, amb.gdir, &amb.corral);
557 greg 2.71 scalecolor(acol, 1./AVGREFL); /* undo assumed reflectance */
558 greg 2.76 if (i <= 0 || amb.rad[0] <= FTINY) /* no Hessian or zero radius */
559 greg 2.73 return(i);
560 greg 2.71 /* store value */
561     VCOPY(amb.pos, r->rop);
562     amb.ndir = encodedir(r->ron);
563 greg 2.72 amb.udir = encodedir(uvw[0]);
564 greg 2.71 amb.lvl = al;
565     copycolor(amb.val, acol);
566     /* insert into tree */
567     avsave(&amb); /* and save to file */
568 greg 2.72 if (rn != r->ron) { /* texture */
569     VCOPY(uvw[2], r->ron);
570     extambient(acol, &amb, r->rop, rn, uvw);
571     }
572 greg 2.71 return(1);
573     }
574    
575    
576 greg 2.92 static int
577 greg 2.71 extambient( /* extrapolate value at pv, nv */
578     COLOR cr,
579     AMBVAL *ap,
580     FVECT pv,
581 greg 2.72 FVECT nv,
582     FVECT uvw[3]
583 greg 2.71 )
584     {
585 greg 2.93 const double min_d = 0.05;
586 greg 2.72 static FVECT my_uvw[3];
587     FVECT v1;
588     int i;
589     double d = 1.0; /* zeroeth order */
590    
591     if (uvw == NULL) { /* need local coordinates? */
592     decodedir(my_uvw[2], ap->ndir);
593     decodedir(my_uvw[0], ap->udir);
594     VCROSS(my_uvw[1], my_uvw[2], my_uvw[0]);
595     uvw = my_uvw;
596     }
597 greg 2.71 for (i = 3; i--; ) /* gradient due to translation */
598     d += (pv[i] - ap->pos[i]) *
599     (ap->gpos[0]*uvw[0][i] + ap->gpos[1]*uvw[1][i]);
600    
601     VCROSS(v1, uvw[2], nv); /* gradient due to rotation */
602     for (i = 3; i--; )
603     d += v1[i] * (ap->gdir[0]*uvw[0][i] + ap->gdir[1]*uvw[1][i]);
604    
605 greg 2.93 if (d < min_d) /* should not use if we can avoid it */
606     d = min_d;
607 greg 2.71 copycolor(cr, ap->val);
608     scalecolor(cr, d);
609 greg 2.93 return(d > min_d);
610 greg 2.71 }
611    
612    
613     static void
614     avinsert( /* insert ambient value in our tree */
615     AMBVAL *av
616     )
617     {
618     AMBTREE *at;
619     AMBVAL *ap;
620     AMBVAL avh;
621     FVECT ck0;
622     double s;
623     int branch;
624     int i;
625    
626     if (av->rad[0] <= FTINY)
627     error(CONSISTENCY, "zero ambient radius in avinsert");
628     at = &atrunk;
629     VCOPY(ck0, thescene.cuorg);
630     s = thescene.cusize;
631 greg 2.84 while (s*(OCTSCALE/2) > av->rad[1]*ambacc) {
632 greg 2.71 if (at->kid == NULL)
633     if ((at->kid = newambtree()) == NULL)
634     error(SYSTEM, "out of memory in avinsert");
635     s *= 0.5;
636     branch = 0;
637     for (i = 0; i < 3; i++)
638     if (av->pos[i] > ck0[i] + s) {
639     ck0[i] += s;
640     branch |= 1 << i;
641     }
642     at = at->kid + branch;
643     }
644     avh.next = at->alist; /* order by increasing level */
645     for (ap = &avh; ap->next != NULL; ap = ap->next)
646 greg 2.87 if ( ap->next->lvl > av->lvl ||
647     (ap->next->lvl == av->lvl) &
648     (ap->next->weight <= av->weight) )
649 greg 2.71 break;
650     av->next = ap->next;
651     ap->next = (AMBVAL*)av;
652     at->alist = avh.next;
653     }
654    
655    
656     static void
657 schorsch 2.56 initambfile( /* initialize ambient file */
658 greg 2.68 int cre8
659 schorsch 2.56 )
660 greg 2.9 {
661 greg 2.47 extern char *progname, *octname;
662     static char *mybuf = NULL;
663 greg 2.9
664 greg 2.19 #ifdef F_SETLKW
665 greg 2.69 aflock(cre8 ? F_WRLCK : F_RDLCK);
666 greg 2.19 #endif
667 schorsch 2.50 SET_FILE_BINARY(ambfp);
668 greg 2.47 if (mybuf == NULL)
669     mybuf = (char *)bmalloc(BUFSIZ+8);
670     setbuf(ambfp, mybuf);
671 greg 2.68 if (cre8) { /* new file */
672 greg 2.24 newheader("RADIANCE", ambfp);
673 greg 2.41 fprintf(ambfp, "%s -av %g %g %g -aw %d -ab %d -aa %g ",
674 greg 2.9 progname, colval(ambval,RED),
675     colval(ambval,GRN), colval(ambval,BLU),
676 greg 2.41 ambvwt, ambounce, ambacc);
677 greg 2.47 fprintf(ambfp, "-ad %d -as %d -ar %d ",
678     ambdiv, ambssamp, ambres);
679     if (octname != NULL)
680 greg 2.68 fputs(octname, ambfp);
681     fputc('\n', ambfp);
682 greg 2.9 fprintf(ambfp, "SOFTWARE= %s\n", VersionID);
683 greg 2.47 fputnow(ambfp);
684 greg 2.9 fputformat(AMBFMT, ambfp);
685 greg 2.68 fputc('\n', ambfp);
686 greg 2.9 putambmagic(ambfp);
687 greg 2.17 } else if (checkheader(ambfp, AMBFMT, NULL) < 0 || !hasambmagic(ambfp))
688     error(USER, "bad ambient file");
689 greg 2.9 }
690    
691    
692 greg 2.47 static void
693 schorsch 2.56 avsave( /* insert and save an ambient value */
694     AMBVAL *av
695     )
696 greg 1.1 {
697 greg 2.88 avstore(av);
698 greg 1.1 if (ambfp == NULL)
699     return;
700 greg 2.5 if (writambval(av, ambfp) < 0)
701 greg 1.1 goto writerr;
702 greg 2.16 if (++nunflshed >= AMBFLUSH)
703 greg 2.7 if (ambsync() == EOF)
704 greg 1.1 goto writerr;
705     return;
706     writerr:
707 greg 2.47 error(SYSTEM, "error writing to ambient file");
708 greg 1.1 }
709    
710    
711 greg 2.20 static AMBVAL *
712 greg 2.88 avstore( /* allocate memory and save aval */
713 greg 2.70 AMBVAL *aval
714 schorsch 2.56 )
715 greg 2.20 {
716 greg 2.70 AMBVAL *av;
717 gregl 2.43 double d;
718 greg 2.20
719     if ((av = newambval()) == NULL)
720     error(SYSTEM, "out of memory in avstore");
721 schorsch 2.53 *av = *aval;
722 greg 2.26 av->latick = ambclock;
723     av->next = NULL;
724     nambvals++;
725 gregl 2.43 d = bright(av->val);
726     if (d > FTINY) { /* add to log sum for averaging */
727     avsum += log(d);
728     navsum++;
729     }
730 greg 2.88 avinsert(av); /* insert in our cache tree */
731 greg 2.20 return(av);
732     }
733    
734    
735 greg 2.26 #define ATALLOCSZ 512 /* #/8 trees to allocate at once */
736    
737     static AMBTREE *atfreelist = NULL; /* free ambient tree structures */
738    
739    
740 greg 2.47 static AMBTREE *
741 schorsch 2.56 newambtree(void) /* allocate 8 ambient tree structs */
742 greg 2.26 {
743 greg 2.70 AMBTREE *atp, *upperlim;
744 greg 2.26
745     if (atfreelist == NULL) { /* get more nodes */
746 greg 2.47 atfreelist = (AMBTREE *)malloc(ATALLOCSZ*8*sizeof(AMBTREE));
747 greg 2.26 if (atfreelist == NULL)
748     return(NULL);
749     /* link new free list */
750     upperlim = atfreelist + 8*(ATALLOCSZ-1);
751     for (atp = atfreelist; atp < upperlim; atp += 8)
752     atp->kid = atp + 8;
753     atp->kid = NULL;
754     }
755     atp = atfreelist;
756     atfreelist = atp->kid;
757 greg 2.102 memset(atp, 0, 8*sizeof(AMBTREE));
758 greg 2.26 return(atp);
759     }
760    
761    
762 greg 2.47 static void
763 schorsch 2.56 freeambtree( /* free 8 ambient tree structs */
764     AMBTREE *atp
765     )
766 greg 2.26 {
767 greg 2.100 atp->kid = atfreelist;
768 greg 2.26 atfreelist = atp;
769     }
770    
771    
772 greg 2.47 static void
773 schorsch 2.56 unloadatree( /* unload an ambient value tree */
774 greg 2.70 AMBTREE *at,
775 schorsch 2.56 unloadtf_t *f
776     )
777 greg 2.20 {
778 greg 2.70 AMBVAL *av;
779     int i;
780 greg 2.20 /* transfer values at this node */
781     for (av = at->alist; av != NULL; av = at->alist) {
782     at->alist = av->next;
783 greg 2.102 av->next = NULL;
784 greg 2.26 (*f)(av);
785 greg 2.20 }
786 greg 2.21 if (at->kid == NULL)
787     return;
788 greg 2.20 for (i = 0; i < 8; i++) /* transfer and free children */
789 greg 2.26 unloadatree(at->kid+i, f);
790 greg 2.20 freeambtree(at->kid);
791 greg 2.26 at->kid = NULL;
792     }
793    
794    
795 greg 2.39 static struct avl {
796     AMBVAL *p;
797     unsigned long t;
798     } *avlist1; /* ambient value list with ticks */
799     static AMBVAL **avlist2; /* memory positions for sorting */
800 greg 2.26 static int i_avlist; /* index for lists */
801    
802 schorsch 2.56 static int alatcmp(const void *av1, const void *av2);
803 greg 2.26
804 schorsch 2.56 static void
805 greg 2.71 avfree(AMBVAL *av)
806     {
807     free(av);
808     }
809    
810     static void
811 schorsch 2.56 av2list(
812 greg 2.71 AMBVAL *av
813 schorsch 2.56 )
814 greg 2.26 {
815 greg 2.27 #ifdef DEBUG
816 greg 2.26 if (i_avlist >= nambvals)
817     error(CONSISTENCY, "too many ambient values in av2list1");
818 greg 2.27 #endif
819 schorsch 2.56 avlist1[i_avlist].p = avlist2[i_avlist] = (AMBVAL*)av;
820 greg 2.71 avlist1[i_avlist++].t = av->latick;
821 greg 2.26 }
822    
823    
824     static int
825 schorsch 2.56 alatcmp( /* compare ambient values for MRA */
826     const void *av1,
827     const void *av2
828     )
829 greg 2.26 {
830 greg 2.70 long lc = ((struct avl *)av2)->t - ((struct avl *)av1)->t;
831 greg 2.38 return(lc<0 ? -1 : lc>0 ? 1 : 0);
832 greg 2.26 }
833    
834    
835 greg 2.47 /* GW NOTE 2002/10/3:
836     * I used to compare AMBVAL pointers, but found that this was the
837     * cause of a serious consistency error with gcc, since the optimizer
838     * uses some dangerous trick in pointer subtraction that
839     * assumes pointers differ by exact struct size increments.
840     */
841 greg 2.26 static int
842 schorsch 2.56 aposcmp( /* compare ambient value positions */
843     const void *avp1,
844     const void *avp2
845     )
846 greg 2.26 {
847 greg 2.70 long diff = *(char * const *)avp1 - *(char * const *)avp2;
848 greg 2.47 if (diff < 0)
849     return(-1);
850     return(diff > 0);
851 greg 2.26 }
852    
853 greg 2.71
854 greg 2.27 static int
855 schorsch 2.56 avlmemi( /* find list position from address */
856     AMBVAL *avaddr
857     )
858 greg 2.27 {
859 greg 2.70 AMBVAL **avlpp;
860 greg 2.27
861 greg 2.102 avlpp = (AMBVAL **)bsearch(&avaddr, avlist2,
862     nambvals, sizeof(AMBVAL *), aposcmp);
863 greg 2.27 if (avlpp == NULL)
864     error(CONSISTENCY, "address not found in avlmemi");
865     return(avlpp - avlist2);
866     }
867    
868    
869 greg 2.47 static void
870 schorsch 2.56 sortambvals( /* resort ambient values */
871     int always
872     )
873 greg 2.26 {
874     AMBTREE oldatrunk;
875     AMBVAL tav, *tap, *pnext;
876 greg 2.70 int i, j;
877 greg 2.28 /* see if it's time yet */
878 greg 2.40 if (!always && (ambclock++ < lastsort+sortintvl ||
879 greg 2.28 nambvals < SORT_THRESH))
880     return;
881 greg 2.26 /*
882     * The idea here is to minimize memory thrashing
883     * in VM systems by improving reference locality.
884     * We do this by periodically sorting our stored ambient
885     * values in memory in order of most recently to least
886     * recently accessed. This ordering was chosen so that new
887     * ambient values (which tend to be less important) go into
888     * higher memory with the infrequently accessed values.
889     * Since we expect our values to need sorting less
890     * frequently as the process continues, we double our
891     * waiting interval after each call.
892     * This routine is also called by setambacc() with
893     * the "always" parameter set to 1 so that the ambient
894     * tree will be rebuilt with the new accuracy parameter.
895     */
896 greg 2.27 if (tracktime) { /* allocate pointer arrays to sort */
897 greg 2.26 avlist2 = (AMBVAL **)malloc(nambvals*sizeof(AMBVAL *));
898 greg 2.39 avlist1 = (struct avl *)malloc(nambvals*sizeof(struct avl));
899     } else {
900     avlist2 = NULL;
901     avlist1 = NULL;
902     }
903     if (avlist1 == NULL) { /* no time tracking -- rebuild tree? */
904     if (avlist2 != NULL)
905 greg 2.102 free(avlist2);
906 greg 2.27 if (always) { /* rebuild without sorting */
907 schorsch 2.53 oldatrunk = atrunk;
908 greg 2.27 atrunk.alist = NULL;
909     atrunk.kid = NULL;
910 greg 2.99 unloadatree(&oldatrunk, avinsert);
911 greg 2.27 }
912 greg 2.26 } else { /* sort memory by last access time */
913     /*
914     * Sorting memory is tricky because it isn't contiguous.
915     * We have to sort an array of pointers by MRA and also
916     * by memory position. We then copy values in "loops"
917     * to minimize memory hits. Nevertheless, we will visit
918 greg 2.27 * everyone at least twice, and this is an expensive process
919 greg 2.26 * when we're thrashing, which is when we need to do it.
920     */
921 greg 2.27 #ifdef DEBUG
922 greg 2.29 sprintf(errmsg, "sorting %u ambient values at ambclock=%lu...",
923     nambvals, ambclock);
924 greg 2.27 eputs(errmsg);
925     #endif
926     i_avlist = 0;
927 greg 2.99 unloadatree(&atrunk, av2list); /* empty current tree */
928 greg 2.27 #ifdef DEBUG
929     if (i_avlist < nambvals)
930     error(CONSISTENCY, "missing ambient values in sortambvals");
931     #endif
932 greg 2.102 qsort(avlist1, nambvals, sizeof(struct avl), alatcmp);
933     qsort(avlist2, nambvals, sizeof(AMBVAL *), aposcmp);
934 greg 2.26 for (i = 0; i < nambvals; i++) {
935 greg 2.39 if (avlist1[i].p == NULL)
936 greg 2.26 continue;
937     tap = avlist2[i];
938 schorsch 2.53 tav = *tap;
939 greg 2.39 for (j = i; (pnext = avlist1[j].p) != tap;
940 greg 2.27 j = avlmemi(pnext)) {
941 schorsch 2.53 *(avlist2[j]) = *pnext;
942 greg 2.26 avinsert(avlist2[j]);
943 greg 2.39 avlist1[j].p = NULL;
944 greg 2.26 }
945 schorsch 2.53 *(avlist2[j]) = tav;
946 greg 2.26 avinsert(avlist2[j]);
947 greg 2.39 avlist1[j].p = NULL;
948 greg 2.26 }
949 greg 2.102 free(avlist1);
950     free(avlist2);
951 greg 2.28 /* compute new sort interval */
952     sortintvl = ambclock - lastsort;
953 greg 2.32 if (sortintvl >= MAX_SORT_INTVL/2)
954 greg 2.28 sortintvl = MAX_SORT_INTVL;
955     else
956 greg 2.26 sortintvl <<= 1; /* wait twice as long next */
957 greg 2.27 #ifdef DEBUG
958     eputs("done\n");
959     #endif
960 greg 2.26 }
961     if (ambclock >= MAXACLOCK)
962     ambclock = MAXACLOCK/2;
963     lastsort = ambclock;
964 greg 2.20 }
965    
966    
967 greg 2.18 #ifdef F_SETLKW
968 greg 2.10
969 greg 2.47 static void
970 schorsch 2.56 aflock( /* lock/unlock ambient file */
971     int typ
972     )
973 greg 2.19 {
974     static struct flock fls; /* static so initialized to zeroes */
975    
976 greg 2.66 if (typ == fls.l_type) /* already called? */
977     return;
978 greg 2.19 fls.l_type = typ;
979     if (fcntl(fileno(ambfp), F_SETLKW, &fls) < 0)
980     error(SYSTEM, "cannot (un)lock ambient file");
981     }
982    
983    
984 greg 2.71 int
985 schorsch 2.56 ambsync(void) /* synchronize ambient file */
986 greg 2.7 {
987 greg 2.15 long flen;
988 greg 2.12 AMBVAL avs;
989 greg 2.70 int n;
990 greg 2.16
991 greg 2.65 if (ambfp == NULL) /* no ambient file? */
992     return(0);
993 greg 2.63 /* gain appropriate access */
994     aflock(nunflshed ? F_WRLCK : F_RDLCK);
995 greg 2.7 /* see if file has grown */
996 greg 2.55 if ((flen = lseek(fileno(ambfp), (off_t)0, SEEK_END)) < 0)
997 greg 2.21 goto seekerr;
998 greg 2.68 if ((n = flen - lastpos) > 0) { /* file has grown */
999 greg 2.104 if (ambinp == NULL) { /* get new file pointer */
1000     ambinp = fopen(ambfile, "rb");
1001 greg 2.14 if (ambinp == NULL)
1002 greg 2.104 error(SYSTEM, "fopen failed in ambsync");
1003 greg 2.14 }
1004 greg 2.66 if (fseek(ambinp, lastpos, SEEK_SET) < 0)
1005 greg 2.21 goto seekerr;
1006 greg 2.15 while (n >= AMBVALSIZ) { /* load contributed values */
1007 greg 2.37 if (!readambval(&avs, ambinp)) {
1008     sprintf(errmsg,
1009 greg 2.47 "ambient file \"%s\" corrupted near character %ld",
1010     ambfile, flen - n);
1011 greg 2.37 error(WARNING, errmsg);
1012     break;
1013     }
1014 greg 2.88 avstore(&avs);
1015 greg 2.15 n -= AMBVALSIZ;
1016     }
1017 greg 2.104 lastpos = flen - n; /* check alignment */
1018     if (n && lseek(fileno(ambfp), (off_t)lastpos, SEEK_SET) < 0)
1019     goto seekerr;
1020 greg 2.7 }
1021     n = fflush(ambfp); /* calls write() at last */
1022 greg 2.104 lastpos += (long)nunflshed*AMBVALSIZ;
1023 greg 2.19 aflock(F_UNLCK); /* release file */
1024 greg 2.16 nunflshed = 0;
1025 greg 2.7 return(n);
1026 greg 2.21 seekerr:
1027     error(SYSTEM, "seek failed in ambsync");
1028 greg 2.104 return(EOF); /* pro forma return */
1029 greg 2.18 }
1030    
1031 greg 2.71 #else /* ! F_SETLKW */
1032 greg 2.18
1033 greg 2.71 int
1034 schorsch 2.56 ambsync(void) /* flush ambient file */
1035 greg 2.18 {
1036 greg 2.65 if (ambfp == NULL)
1037     return(0);
1038 greg 2.18 nunflshed = 0;
1039     return(fflush(ambfp));
1040 greg 1.1 }
1041 greg 2.10
1042 greg 2.71 #endif /* ! F_SETLKW */