ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/ambient.c
Revision: 2.129
Committed: Thu Jan 23 02:15:33 2025 UTC (3 months, 1 week ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.128: +82 -117 lines
Log Message:
perf: Simplified operation of shared ambient files, more testing needed

File Contents

# User Rev Content
1 greg 2.129 static const char RCSid[] = "$Id: ambient.c,v 2.128 2024/12/19 19:34:31 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.57 #ifndef MAXASET
26 greg 2.67 #define MAXASET 4095 /* maximum number of elements in ambient set */
27 greg 2.57 #endif
28 greg 2.12 OBJECT ambset[MAXASET+1]={0}; /* ambient include/exclude set */
29 greg 1.1
30 greg 2.12 double maxarad; /* maximum ambient radius */
31     double minarad; /* minimum ambient radius */
32 greg 1.1
33 greg 2.12 static AMBTREE atrunk; /* our ambient trunk node */
34 greg 1.1
35     static FILE *ambfp = NULL; /* ambient file pointer */
36 greg 2.129 static int nunflshed; /* number of unflushed ambient values */
37     static FILE *ambinp = NULL; /* input pointer for ambient i/o */
38 greg 1.1
39 gregl 2.43 static double avsum = 0.; /* computed ambient value sum (log) */
40     static unsigned int navsum = 0; /* number of values in avsum */
41 greg 2.35 static unsigned int nambvals = 0; /* total number of indirect values */
42 greg 2.129 static off_t lastpos = -1; /* last flush position */
43 greg 2.26
44 greg 2.12 #define AMBFLUSH (BUFSIZ/AMBVALSIZ)
45 greg 2.6
46 greg 2.118 #define AVSIZE (sizeof(AMBVAL)-sizeof(SCOLOR)+sizeof(COLORV)*NCSAMP)
47     #define newambval() (AMBVAL *)malloc(AVSIZE)
48 greg 1.1
49 greg 2.117 #define tfunc(x0, x, x1) (((x)-(x0))/((x1)-(x0)))
50 greg 2.108
51 greg 2.118 static void initambfile(int cre8);
52 schorsch 2.56 static void avsave(AMBVAL *av);
53     static AMBVAL *avstore(AMBVAL *aval);
54     static AMBTREE *newambtree(void);
55     static void freeambtree(AMBTREE *atp);
56    
57 greg 2.71 typedef void unloadtf_t(AMBVAL *);
58 schorsch 2.56 static unloadtf_t avinsert;
59 greg 2.71 static unloadtf_t avfree;
60 schorsch 2.56 static void unloadatree(AMBTREE *at, unloadtf_t *f);
61    
62 greg 2.112 static void sortambvals(void);
63 schorsch 2.56
64 greg 2.108 static int plugaleak(RAY *r, AMBVAL *ap, FVECT anorm, double ang);
65 greg 2.118 static double sumambient(SCOLOR acol, RAY *r, FVECT rn, int al,
66 greg 2.108 AMBTREE *at, FVECT c0, double s);
67 greg 2.118 static int makeambient(SCOLOR acol, RAY *r, FVECT rn, int al);
68     static int extambient(SCOLOR cr, AMBVAL *ap, FVECT pv, FVECT nv,
69 greg 2.108 FVECT uvw[3]);
70    
71 greg 2.7
72 greg 2.71 void
73 schorsch 2.56 setambres( /* set ambient resolution */
74     int ar
75     )
76 greg 2.3 {
77 greg 2.21 ambres = ar < 0 ? 0 : ar; /* may be done already */
78 greg 2.3 /* set min & max radii */
79     if (ar <= 0) {
80 greg 2.25 minarad = 0;
81 greg 2.89 maxarad = thescene.cusize*0.2;
82 greg 2.3 } else {
83     minarad = thescene.cusize / ar;
84 greg 2.83 maxarad = 64.0 * minarad; /* heuristic */
85 greg 2.89 if (maxarad > thescene.cusize*0.2)
86     maxarad = thescene.cusize*0.2;
87 greg 2.3 }
88 greg 2.25 if (minarad <= FTINY)
89 greg 2.83 minarad = 10.0*FTINY;
90 greg 2.25 if (maxarad <= minarad)
91 greg 2.83 maxarad = 64.0 * minarad;
92 greg 2.3 }
93    
94    
95 greg 2.71 void
96 schorsch 2.56 setambacc( /* set ambient accuracy */
97     double newa
98     )
99 greg 2.20 {
100 greg 2.84 static double olda; /* remember previous setting here */
101 greg 2.82
102     newa *= (newa > 0);
103     if (fabs(newa - olda) >= .05*(newa + olda)) {
104 greg 2.84 ambacc = newa;
105 greg 2.107 if (ambacc > FTINY && nambvals > 0)
106 greg 2.112 sortambvals(); /* rebuild tree */
107 greg 2.74 }
108 greg 2.20 }
109    
110    
111 greg 2.71 void
112 schorsch 2.56 setambient(void) /* initialize calculation */
113 greg 1.1 {
114 greg 2.129 int exists;
115     off_t flen;
116 greg 2.12 AMBVAL amb;
117 greg 2.47 /* make sure we're fresh */
118     ambdone();
119 greg 2.3 /* init ambient limits */
120     setambres(ambres);
121 greg 2.82 setambacc(ambacc);
122 greg 2.47 if (ambfile == NULL || !ambfile[0])
123 greg 2.19 return;
124 greg 2.20 if (ambacc <= FTINY) {
125 greg 2.21 sprintf(errmsg, "zero ambient accuracy so \"%s\" not opened",
126 greg 2.47 ambfile);
127 greg 2.20 error(WARNING, errmsg);
128     return;
129     }
130 greg 2.129 exists = access(ambfile, F_OK) == 0; /* check existence, first */
131     ambfp = fopen(ambfile, "a+"); /* try read/append */
132     if (!exists & (ambfp == NULL)) {
133     sprintf(errmsg, "cannot create ambient file \"%s\"", ambfile);
134     error(SYSTEM, errmsg);
135     }
136     if (ambfp == NULL) { /* try opening read-only? */
137     if ((ambfp = fopen(ambfile, "r")) == NULL) {
138     sprintf(errmsg,
139     "cannot open ambient file \"%s\" for reading",
140     ambfile);
141     error(SYSTEM, errmsg);
142     }
143     exists = -1; /* flag read-only */
144     } else if (exists)
145     rewind(ambfp); /* XXX not necessary? */
146    
147     if (exists) {
148     initambfile(0); /* file already exists */
149 greg 2.66 lastpos = ftell(ambfp);
150 greg 2.19 while (readambval(&amb, ambfp))
151 greg 2.129 avstore(&amb); /* load what we can */
152     if (exists < 0) { /* read-only? */
153 gwlarson 2.46 sprintf(errmsg,
154     "loaded %u values from read-only ambient file",
155     nambvals);
156     error(WARNING, errmsg);
157     fclose(ambfp); /* close file so no writes */
158     ambfp = NULL;
159 greg 2.129 return;
160 gwlarson 2.46 }
161     /* align file pointer */
162 greg 2.129 lastpos += (off_t)nambvals*AMBVALSIZ;
163     flen = lseek(fileno(ambfp), 0, SEEK_END);
164 greg 2.66 if (flen != lastpos) {
165 greg 2.42 sprintf(errmsg,
166 greg 2.129 "ignoring last %lu values in ambient file (corrupted)",
167     (unsigned long)((flen - lastpos)/AMBVALSIZ));
168 greg 2.42 error(WARNING, errmsg);
169 greg 2.129 /* fseek() not needed? */
170     fseeko(ambfp, lastpos, SEEK_SET);
171     ftruncate(fileno(ambfp), lastpos);
172 greg 2.37 }
173 greg 2.129 } else {
174     initambfile(1); /* else start new file */
175 greg 2.68 fflush(ambfp);
176 greg 2.66 lastpos = ftell(ambfp);
177 greg 2.16 }
178 greg 1.8 }
179    
180    
181 greg 2.71 void
182 schorsch 2.56 ambdone(void) /* close ambient file and free memory */
183 greg 2.47 {
184     if (ambfp != NULL) { /* close ambient file */
185 greg 2.129 fclose(ambfp); /* don't call ambsync() */
186 greg 2.47 ambfp = NULL;
187 greg 2.129 lastpos = -1;
188     if (ambinp != NULL) {
189 greg 2.47 fclose(ambinp);
190     ambinp = NULL;
191     }
192     }
193     /* free ambient tree */
194 greg 2.99 unloadatree(&atrunk, avfree);
195 greg 2.47 /* reset state variables */
196     avsum = 0.;
197     navsum = 0;
198     nambvals = 0;
199     }
200    
201    
202 greg 2.71 void
203 schorsch 2.56 ambnotify( /* record new modifier */
204     OBJECT obj
205     )
206 greg 1.8 {
207 greg 1.11 static int hitlimit = 0;
208 greg 2.70 OBJREC *o;
209     char **amblp;
210 greg 1.8
211 greg 2.47 if (obj == OVOID) { /* starting over */
212     ambset[0] = 0;
213     hitlimit = 0;
214     return;
215     }
216     o = objptr(obj);
217 greg 1.11 if (hitlimit || !ismodifier(o->otype))
218 greg 1.8 return;
219     for (amblp = amblist; *amblp != NULL; amblp++)
220     if (!strcmp(o->oname, *amblp)) {
221 greg 1.11 if (ambset[0] >= MAXASET) {
222     error(WARNING, "too many modifiers in ambient list");
223     hitlimit++;
224     return; /* should this be fatal? */
225     }
226 greg 1.8 insertelem(ambset, obj);
227     return;
228 greg 1.1 }
229     }
230    
231 greg 2.71
232     void
233     multambient( /* compute ambient component & multiply by coef. */
234 greg 2.118 SCOLOR aval,
235 greg 2.71 RAY *r,
236     FVECT nrm
237     )
238     {
239 greg 2.111 static double logAvgAbsorp = 1;
240 greg 2.71 static int rdepth = 0; /* ambient recursion */
241 greg 2.118 SCOLOR acol, caustic;
242 greg 2.103 int i, ok;
243 greg 2.71 double d, l;
244    
245 greg 2.94 /* PMAP: Factor in ambient from photon map, if enabled and ray is
246     * ambient. Return as all ambient components accounted for, else
247     * continue. */
248     if (ambPmap(aval, r, rdepth))
249     return;
250    
251 greg 2.111 if (logAvgAbsorp > 0) /* exclude in -aw to avoid growth */
252     logAvgAbsorp = log(1.-AVGREFL);
253    
254 greg 2.94 /* PMAP: Factor in specular-diffuse ambient (caustics) from photon
255     * map, if enabled and ray is primary, else caustic is zero. Continue
256     * with RADIANCE ambient calculation */
257 greg 2.118 {/* XXX TEMPORARY */
258     COLOR pmc;
259     scolor_color(pmc, aval);
260     ambPmapCaustic(pmc, r, rdepth);
261     setscolor(caustic, colval(pmc,RED), colval(pmc,GRN), colval(pmc,BLU));
262     }
263 greg 2.71 if (ambdiv <= 0) /* no ambient calculation */
264     goto dumbamb;
265     /* check number of bounces */
266     if (rdepth >= ambounce)
267     goto dumbamb;
268     /* check ambient list */
269     if (ambincl != -1 && r->ro != NULL &&
270     ambincl != inset(ambset, r->ro->omod))
271     goto dumbamb;
272    
273 greg 2.123 if (ambacc <= FTINY) { /* no ambient storage? */
274     double rdot = DOT(nrm,r->ron);
275     int sgn = 1 - 2*(rdot < 0);
276     float dgrad[2], *dgp = NULL;
277 greg 2.103 FVECT uvd[2];
278    
279 greg 2.123 if (sgn*rdot < 0.9999)
280 greg 2.103 dgp = dgrad; /* compute rotational grad. */
281 greg 2.118 copyscolor(acol, aval);
282 greg 2.71 rdepth++;
283 greg 2.123 ok = doambient(acol, r, r->rweight*sgn,
284 greg 2.103 uvd, NULL, NULL, dgp, NULL);
285 greg 2.71 rdepth--;
286     if (!ok)
287     goto dumbamb;
288 greg 2.103 if ((ok > 0) & (dgp != NULL)) { /* apply texture */
289     FVECT v1;
290     VCROSS(v1, r->ron, nrm);
291     d = 1.0;
292     for (i = 3; i--; )
293 greg 2.123 d += sgn*v1[i] * (dgp[0]*uvd[0][i] + dgp[1]*uvd[1][i]);
294 greg 2.103 if (d >= 0.05)
295 greg 2.118 scalescolor(acol, d);
296 greg 2.103 }
297 greg 2.118 copyscolor(aval, acol);
298 greg 2.94
299     /* PMAP: add in caustic */
300 greg 2.118 saddscolor(aval, caustic);
301 greg 2.71 return;
302     }
303     /* interpolate ambient value */
304 greg 2.118 scolorblack(acol);
305 greg 2.71 d = sumambient(acol, r, nrm, rdepth,
306     &atrunk, thescene.cuorg, thescene.cusize);
307 greg 2.94
308 greg 2.71 if (d > FTINY) {
309 greg 2.118 scalescolor(acol, 1.0/d);
310     smultscolor(aval, acol);
311 greg 2.94
312     /* PMAP: add in caustic */
313 greg 2.118 saddscolor(aval, caustic);
314 greg 2.71 return;
315     }
316 greg 2.94
317 greg 2.71 rdepth++; /* need to cache new value */
318     ok = makeambient(acol, r, nrm, rdepth-1);
319     rdepth--;
320 greg 2.94
321 greg 2.71 if (ok) {
322 greg 2.118 smultscolor(aval, acol); /* computed new value */
323 greg 2.94
324     /* PMAP: add in caustic */
325 greg 2.118 saddscolor(aval, caustic);
326 greg 2.71 return;
327     }
328 greg 2.94
329 greg 2.71 dumbamb: /* return global value */
330     if ((ambvwt <= 0) | (navsum == 0)) {
331 greg 2.118 smultcolor(aval, ambval);
332 greg 2.94
333     /* PMAP: add in caustic */
334 greg 2.118 saddscolor(aval, caustic);
335 greg 2.71 return;
336     }
337 greg 2.94
338     l = bright(ambval); /* average in computations */
339 greg 2.71 if (l > FTINY) {
340 greg 2.111 d = (log(l)*(double)ambvwt + avsum + logAvgAbsorp*navsum) /
341 greg 2.71 (double)(ambvwt + navsum);
342     d = exp(d) / l;
343 greg 2.118 scalescolor(aval, d);
344     smultcolor(aval, ambval); /* apply color of ambval */
345 greg 2.71 } else {
346 greg 2.111 d = exp( avsum/(double)navsum + logAvgAbsorp );
347 greg 2.118 scalescolor(aval, d); /* neutral color */
348 greg 2.71 }
349     }
350    
351    
352 greg 2.86 /* Plug a potential leak where ambient cache value is occluded */
353     static int
354     plugaleak(RAY *r, AMBVAL *ap, FVECT anorm, double ang)
355     {
356     const double cost70sq = 0.1169778; /* cos(70deg)^2 */
357     RAY rtst;
358     FVECT vdif;
359     double normdot, ndotd, nadotd;
360     double a, b, c, t[2];
361    
362     ang += 2.*PI*(ang < 0); /* check direction flags */
363     if ( !(ap->corral>>(int)(ang*(16./PI)) & 1) )
364     return(0);
365     /*
366     * Generate test ray, targeting 20 degrees above sample point plane
367     * along surface normal from cache position. This should be high
368     * enough to miss local geometry we don't really care about.
369     */
370     VSUB(vdif, ap->pos, r->rop);
371     normdot = DOT(anorm, r->ron);
372     ndotd = DOT(vdif, r->ron);
373     nadotd = DOT(vdif, anorm);
374     a = normdot*normdot - cost70sq;
375     b = 2.0*(normdot*ndotd - nadotd*cost70sq);
376     c = ndotd*ndotd - DOT(vdif,vdif)*cost70sq;
377     if (quadratic(t, a, b, c) != 2)
378     return(1); /* should rarely happen */
379     if (t[1] <= FTINY)
380     return(0); /* should fail behind test */
381     rayorigin(&rtst, SHADOW, r, NULL);
382     VSUM(rtst.rdir, vdif, anorm, t[1]); /* further dist. > plane */
383     rtst.rmax = normalize(rtst.rdir); /* short ray test */
384     while (localhit(&rtst, &thescene)) { /* check for occluder */
385 greg 2.105 OBJREC *m = findmaterial(rtst.ro);
386 greg 2.126 if (m != NULL && !istransp(m) && !isBSDFproxy(m) &&
387 greg 2.86 (rtst.clipset == NULL ||
388     !inset(rtst.clipset, rtst.ro->omod)))
389     return(1); /* plug light leak */
390     VCOPY(rtst.rorg, rtst.rop); /* skip invisible surface */
391     rtst.rmax -= rtst.rot;
392     rayclear(&rtst);
393     }
394     return(0); /* seems we're OK */
395     }
396    
397    
398     static double
399 greg 2.72 sumambient( /* get interpolated ambient value */
400 greg 2.118 SCOLOR acol,
401 greg 2.71 RAY *r,
402     FVECT rn,
403     int al,
404     AMBTREE *at,
405     FVECT c0,
406     double s
407     )
408 greg 2.80 { /* initial limit is 10 degrees plus ambacc radians */
409     const double minangle = 10.0 * PI/180.;
410 greg 2.123 const int sgn = 1 - 2*(DOT(r->ron,rn) < 0);
411 greg 2.84 double maxangle = minangle + ambacc;
412 greg 2.71 double wsum = 0.0;
413     FVECT ck0;
414     int i, j;
415     AMBVAL *av;
416 greg 2.83
417     if (at->kid != NULL) { /* sum children first */
418     s *= 0.5;
419     for (i = 0; i < 8; i++) {
420     for (j = 0; j < 3; j++) {
421     ck0[j] = c0[j];
422     if (1<<j & i)
423     ck0[j] += s;
424     if (r->rop[j] < ck0[j] - OCTSCALE*s)
425     break;
426     if (r->rop[j] > ck0[j] + (1.0+OCTSCALE)*s)
427     break;
428     }
429     if (j == 3)
430     wsum += sumambient(acol, r, rn, al,
431     at->kid+i, ck0, s);
432     }
433     /* good enough? */
434 greg 2.117 if ((wsum >= 0.05) & (s*ambacc > minarad))
435 greg 2.83 return(wsum);
436     }
437 greg 2.84 /* adjust maximum angle */
438     if (at->alist != NULL && (at->alist->lvl <= al) & (r->rweight < 0.6))
439     maxangle = (maxangle - PI/2.)*pow(r->rweight,0.13) + PI/2.;
440 greg 2.71 /* sum this node */
441     for (av = at->alist; av != NULL; av = av->next) {
442 greg 2.86 double u, v, d, delta_r2, delta_t2;
443 greg 2.118 SCOLOR sct;
444 greg 2.71 FVECT uvw[3];
445     /*
446 greg 2.72 * Ambient level test
447 greg 2.71 */
448 greg 2.87 if (av->lvl > al || /* list sorted, so this works */
449     (av->lvl == al) & (av->weight < 0.9*r->rweight))
450 greg 2.71 break;
451     /*
452 greg 2.72 * Direction test using unperturbed normal
453 greg 2.71 */
454     decodedir(uvw[2], av->ndir);
455 greg 2.123 d = sgn * DOT(uvw[2], r->ron);
456 greg 2.71 if (d <= 0.0) /* >= 90 degrees */
457     continue;
458     delta_r2 = 2.0 - 2.0*d; /* approx. radians^2 */
459     if (delta_r2 >= maxangle*maxangle)
460     continue;
461     /*
462 greg 2.125 * Modified ray behind test
463 greg 2.81 */
464 greg 2.125 VSUB(ck0, r->rop, av->pos);
465     d = DOT(ck0, uvw[2]);
466 greg 2.116 if (d < -minarad*ambacc)
467 greg 2.81 continue;
468 greg 2.125 d /= av->rad[0];
469     delta_t2 = d*d;
470     if (delta_t2 >= ambacc*ambacc)
471 greg 2.81 continue;
472     /*
473 greg 2.72 * Elliptical radii test based on Hessian
474 greg 2.71 */
475     decodedir(uvw[0], av->udir);
476     VCROSS(uvw[1], uvw[2], uvw[0]);
477 greg 2.86 d = (u = DOT(ck0, uvw[0])) / av->rad[0];
478 greg 2.81 delta_t2 += d*d;
479 greg 2.86 d = (v = DOT(ck0, uvw[1])) / av->rad[1];
480 greg 2.71 delta_t2 += d*d;
481 greg 2.84 if (delta_t2 >= ambacc*ambacc)
482 greg 2.71 continue;
483     /*
484 greg 2.86 * Test for potential light leak
485     */
486     if (av->corral && plugaleak(r, av, uvw[2], atan2a(v,u)))
487     continue;
488     /*
489 greg 2.72 * Extrapolate value and compute final weight (hat function)
490 greg 2.71 */
491 greg 2.118 if (!extambient(sct, av, r->rop, rn, uvw))
492 greg 2.92 continue;
493 greg 2.71 d = tfunc(maxangle, sqrt(delta_r2), 0.0) *
494 greg 2.84 tfunc(ambacc, sqrt(delta_t2), 0.0);
495 greg 2.118 scalescolor(sct, d);
496     saddscolor(acol, sct);
497 greg 2.72 wsum += d;
498 greg 2.71 }
499     return(wsum);
500     }
501    
502    
503 greg 2.86 static int
504 greg 2.71 makeambient( /* make a new ambient value for storage */
505 greg 2.118 SCOLOR acol,
506 greg 2.71 RAY *r,
507     FVECT rn,
508     int al
509     )
510     {
511 greg 2.123 int sgn = 1 - 2*(DOT(r->ron,rn) < 0);
512 greg 2.71 AMBVAL amb;
513 greg 2.72 FVECT uvw[3];
514 greg 2.71 int i;
515    
516     amb.weight = 1.0; /* compute weight */
517     for (i = al; i-- > 0; )
518     amb.weight *= AVGREFL;
519     if (r->rweight < 0.1*amb.weight) /* heuristic override */
520     amb.weight = 1.25*r->rweight;
521 greg 2.118 setscolor(acol, AVGREFL, AVGREFL, AVGREFL);
522 greg 2.71 /* compute ambient */
523 greg 2.123 i = doambient(acol, r, amb.weight*sgn,
524 greg 2.86 uvw, amb.rad, amb.gpos, amb.gdir, &amb.corral);
525 greg 2.118 scalescolor(acol, 1./AVGREFL); /* undo assumed reflectance */
526 greg 2.76 if (i <= 0 || amb.rad[0] <= FTINY) /* no Hessian or zero radius */
527 greg 2.73 return(i);
528 greg 2.123 uvw[2][0] = sgn*r->ron[0]; /* orient unperturbed normal */
529     uvw[2][1] = sgn*r->ron[1];
530     uvw[2][2] = sgn*r->ron[2];
531 greg 2.71 /* store value */
532     VCOPY(amb.pos, r->rop);
533 greg 2.123 amb.ndir = encodedir(uvw[2]);
534 greg 2.72 amb.udir = encodedir(uvw[0]);
535 greg 2.71 amb.lvl = al;
536 greg 2.118 copyscolor(amb.val, acol);
537 greg 2.123 avsave(&amb); /* insert and save to file */
538     if (DOT(uvw[2],rn) < 0.9999) /* texture? */
539 greg 2.72 extambient(acol, &amb, r->rop, rn, uvw);
540 greg 2.71 return(1);
541     }
542    
543    
544 greg 2.92 static int
545 greg 2.71 extambient( /* extrapolate value at pv, nv */
546 greg 2.118 SCOLOR scr,
547 greg 2.71 AMBVAL *ap,
548     FVECT pv,
549 greg 2.72 FVECT nv,
550     FVECT uvw[3]
551 greg 2.71 )
552     {
553 greg 2.93 const double min_d = 0.05;
554 greg 2.110 const double max_d = 20.;
555 greg 2.72 static FVECT my_uvw[3];
556     FVECT v1;
557     int i;
558     double d = 1.0; /* zeroeth order */
559    
560     if (uvw == NULL) { /* need local coordinates? */
561     decodedir(my_uvw[2], ap->ndir);
562     decodedir(my_uvw[0], ap->udir);
563     VCROSS(my_uvw[1], my_uvw[2], my_uvw[0]);
564     uvw = my_uvw;
565     }
566 greg 2.71 for (i = 3; i--; ) /* gradient due to translation */
567     d += (pv[i] - ap->pos[i]) *
568     (ap->gpos[0]*uvw[0][i] + ap->gpos[1]*uvw[1][i]);
569    
570     VCROSS(v1, uvw[2], nv); /* gradient due to rotation */
571     for (i = 3; i--; )
572     d += v1[i] * (ap->gdir[0]*uvw[0][i] + ap->gdir[1]*uvw[1][i]);
573    
574 greg 2.110 if (d < min_d) /* clamp min/max scaling */
575 greg 2.93 d = min_d;
576 greg 2.110 else if (d > max_d)
577     d = max_d;
578 greg 2.118 copyscolor(scr, ap->val);
579     scalescolor(scr, d);
580 greg 2.93 return(d > min_d);
581 greg 2.71 }
582    
583    
584     static void
585     avinsert( /* insert ambient value in our tree */
586     AMBVAL *av
587     )
588     {
589     AMBTREE *at;
590     AMBVAL *ap;
591     AMBVAL avh;
592     FVECT ck0;
593     double s;
594     int branch;
595     int i;
596    
597     if (av->rad[0] <= FTINY)
598     error(CONSISTENCY, "zero ambient radius in avinsert");
599     at = &atrunk;
600     VCOPY(ck0, thescene.cuorg);
601     s = thescene.cusize;
602 greg 2.84 while (s*(OCTSCALE/2) > av->rad[1]*ambacc) {
603 greg 2.71 if (at->kid == NULL)
604     if ((at->kid = newambtree()) == NULL)
605     error(SYSTEM, "out of memory in avinsert");
606     s *= 0.5;
607     branch = 0;
608     for (i = 0; i < 3; i++)
609     if (av->pos[i] > ck0[i] + s) {
610     ck0[i] += s;
611     branch |= 1 << i;
612     }
613     at = at->kid + branch;
614     }
615     avh.next = at->alist; /* order by increasing level */
616     for (ap = &avh; ap->next != NULL; ap = ap->next)
617 greg 2.87 if ( ap->next->lvl > av->lvl ||
618     (ap->next->lvl == av->lvl) &
619     (ap->next->weight <= av->weight) )
620 greg 2.71 break;
621     av->next = ap->next;
622     ap->next = (AMBVAL*)av;
623     at->alist = avh.next;
624     }
625    
626    
627     static void
628 schorsch 2.56 initambfile( /* initialize ambient file */
629 greg 2.68 int cre8
630 schorsch 2.56 )
631 greg 2.9 {
632 greg 2.47 extern char *progname, *octname;
633     static char *mybuf = NULL;
634 greg 2.128 int ntries = 3;
635 greg 2.9
636 greg 2.129 if (!AMBFLUSH)
637     error(INTERNAL, "BUFSIZ too small in initambfile");
638 schorsch 2.50 SET_FILE_BINARY(ambfp);
639 greg 2.47 if (mybuf == NULL)
640 greg 2.129 mybuf = (char *)bmalloc(BUFSIZ);
641 greg 2.47 setbuf(ambfp, mybuf);
642 greg 2.129 nunflshed = 0;
643 greg 2.121 retry:
644 greg 2.68 if (cre8) { /* new file */
645 greg 2.24 newheader("RADIANCE", ambfp);
646 greg 2.41 fprintf(ambfp, "%s -av %g %g %g -aw %d -ab %d -aa %g ",
647 greg 2.9 progname, colval(ambval,RED),
648     colval(ambval,GRN), colval(ambval,BLU),
649 greg 2.41 ambvwt, ambounce, ambacc);
650 greg 2.47 fprintf(ambfp, "-ad %d -as %d -ar %d ",
651     ambdiv, ambssamp, ambres);
652 greg 2.115 fprintf(ambfp, "-dr %d -ds %g -dt %g -dc %g ", directrelay,
653     srcsizerat, shadthresh, shadcert);
654     fprintf(ambfp, "-ss %g -st %g -lr %d -lw %g ", specjitter,
655     specthresh, maxdepth, minweight);
656 greg 2.119 fprintf(ambfp, "-cw %g %g -cs %d ", WLPART[3], WLPART[0], NCSAMP);
657 greg 2.47 if (octname != NULL)
658 greg 2.68 fputs(octname, ambfp);
659 greg 2.120 fputc('\n', ambfp); /* end of command line, not header! */
660 greg 2.9 fprintf(ambfp, "SOFTWARE= %s\n", VersionID);
661 greg 2.47 fputnow(ambfp);
662 greg 2.120 AMB_CNDX = CNDX; /* use current spectral sampling */
663     AMB_WLPART = WLPART;
664 greg 2.118 fputwlsplit(WLPART, ambfp);
665     fputncomp(NCSAMP, ambfp);
666 greg 2.9 fputformat(AMBFMT, ambfp);
667 greg 2.68 fputc('\n', ambfp);
668 greg 2.9 putambmagic(ambfp);
669 greg 2.121 } else if (getheader(ambfp, amb_headline, NULL) < 0 || !hasambmagic(ambfp)) {
670 greg 2.128 if (--ntries > 0 && ftell(ambfp) == 0) {
671 greg 2.122 clearerr(ambfp);
672 greg 2.121 sleep(2);
673     goto retry;
674     }
675 greg 2.118 error(USER, "bad/incompatible ambient file");
676 greg 2.121 }
677 greg 2.118 if ((AMB_CNDX != CNDX) | (AMB_WLPART != WLPART)) {
678     if (setspectrsamp(AMB_CNDX, AMB_WLPART) < 0)
679     error(USER, "bad wavelength sampling in ambient file");
680     if (AMB_CNDX[3] == CNDX[3] && FABSEQ(AMB_WLPART[0],WLPART[0]) &&
681     FABSEQ(AMB_WLPART[3],WLPART[3])) {
682     AMB_CNDX = CNDX;
683     AMB_WLPART = WLPART; /* just the same */
684     } else
685     error(WARNING, "different ambient file wavelength sampling");
686     }
687 greg 2.9 }
688    
689    
690 greg 2.47 static void
691 schorsch 2.56 avsave( /* insert and save an ambient value */
692     AMBVAL *av
693     )
694 greg 1.1 {
695 greg 2.88 avstore(av);
696 greg 1.1 if (ambfp == NULL)
697     return;
698 greg 2.5 if (writambval(av, ambfp) < 0)
699 greg 1.1 goto writerr;
700 greg 2.16 if (++nunflshed >= AMBFLUSH)
701 greg 2.7 if (ambsync() == EOF)
702 greg 1.1 goto writerr;
703     return;
704     writerr:
705 greg 2.47 error(SYSTEM, "error writing to ambient file");
706 greg 1.1 }
707    
708    
709 greg 2.20 static AMBVAL *
710 greg 2.88 avstore( /* allocate memory and save aval */
711 greg 2.70 AMBVAL *aval
712 schorsch 2.56 )
713 greg 2.20 {
714 greg 2.70 AMBVAL *av;
715 gregl 2.43 double d;
716 greg 2.20
717     if ((av = newambval()) == NULL)
718     error(SYSTEM, "out of memory in avstore");
719 greg 2.118 memcpy(av, aval, AVSIZE); /* AVSIZE <= sizeof(AMBVAL) */
720 greg 2.26 av->next = NULL;
721     nambvals++;
722 greg 2.118 d = pbright(av->val);
723 gregl 2.43 if (d > FTINY) { /* add to log sum for averaging */
724     avsum += log(d);
725     navsum++;
726     }
727 greg 2.88 avinsert(av); /* insert in our cache tree */
728 greg 2.20 return(av);
729     }
730    
731    
732 greg 2.26 #define ATALLOCSZ 512 /* #/8 trees to allocate at once */
733    
734     static AMBTREE *atfreelist = NULL; /* free ambient tree structures */
735    
736    
737 greg 2.47 static AMBTREE *
738 schorsch 2.56 newambtree(void) /* allocate 8 ambient tree structs */
739 greg 2.26 {
740 greg 2.70 AMBTREE *atp, *upperlim;
741 greg 2.26
742     if (atfreelist == NULL) { /* get more nodes */
743 greg 2.47 atfreelist = (AMBTREE *)malloc(ATALLOCSZ*8*sizeof(AMBTREE));
744 greg 2.26 if (atfreelist == NULL)
745     return(NULL);
746     /* link new free list */
747     upperlim = atfreelist + 8*(ATALLOCSZ-1);
748     for (atp = atfreelist; atp < upperlim; atp += 8)
749     atp->kid = atp + 8;
750     atp->kid = NULL;
751     }
752     atp = atfreelist;
753     atfreelist = atp->kid;
754 greg 2.102 memset(atp, 0, 8*sizeof(AMBTREE));
755 greg 2.26 return(atp);
756     }
757    
758    
759 greg 2.47 static void
760 schorsch 2.56 freeambtree( /* free 8 ambient tree structs */
761     AMBTREE *atp
762     )
763 greg 2.26 {
764 greg 2.100 atp->kid = atfreelist;
765 greg 2.26 atfreelist = atp;
766     }
767    
768    
769 greg 2.47 static void
770 schorsch 2.56 unloadatree( /* unload an ambient value tree */
771 greg 2.70 AMBTREE *at,
772 schorsch 2.56 unloadtf_t *f
773     )
774 greg 2.20 {
775 greg 2.70 AMBVAL *av;
776     int i;
777 greg 2.20 /* transfer values at this node */
778     for (av = at->alist; av != NULL; av = at->alist) {
779     at->alist = av->next;
780 greg 2.102 av->next = NULL;
781 greg 2.26 (*f)(av);
782 greg 2.20 }
783 greg 2.21 if (at->kid == NULL)
784     return;
785 greg 2.20 for (i = 0; i < 8; i++) /* transfer and free children */
786 greg 2.26 unloadatree(at->kid+i, f);
787 greg 2.20 freeambtree(at->kid);
788 greg 2.26 at->kid = NULL;
789     }
790    
791    
792 schorsch 2.56 static void
793 greg 2.71 avfree(AMBVAL *av)
794     {
795     free(av);
796     }
797    
798 greg 2.112
799 greg 2.71 static void
800 greg 2.112 sortambvals(void) /* resort ambient values */
801 greg 2.27 {
802 greg 2.112 AMBTREE oldatrunk = atrunk;
803 greg 2.27
804 greg 2.112 atrunk.alist = NULL;
805     atrunk.kid = NULL;
806     unloadatree(&oldatrunk, avinsert);
807 greg 2.20 }
808    
809    
810 greg 2.71 int
811 schorsch 2.56 ambsync(void) /* synchronize ambient file */
812 greg 2.7 {
813 greg 2.129 off_t newpos;
814     int n;
815 greg 2.12 AMBVAL avs;
816 greg 2.16
817 greg 2.65 if (ambfp == NULL) /* no ambient file? */
818     return(0);
819 greg 2.129
820     if (nunflshed > 0) { /* append new values? */
821     if (fflush(ambfp) < 0)
822     return(EOF);
823     } else if (fseeko(ambfp, 0, SEEK_END) < 0)
824     goto seekerr;
825    
826     if ((newpos = ftello(ambfp)) < 0)
827 greg 2.21 goto seekerr;
828 greg 2.129 /* how many others added? */
829     n = (newpos - lastpos)/AMBVALSIZ - nunflshed;
830     nunflshed = 0;
831     if (n <= 0) { /* no one helping this time? */
832     lastpos = newpos;
833     return(0);
834     }
835     if (ambinp == NULL) { /* else need to open for input? */
836     ambinp = fopen(ambfile, "r");
837     if (ambinp == NULL) {
838     sprintf(errmsg, "cannot reopen ambient file \"%s\"",
839     ambfile);
840     error(SYSTEM, errmsg);
841 greg 2.14 }
842 greg 2.129 SET_FILE_BINARY(ambinp);
843     }
844     /* read from last endpoint */
845     if (fseeko(ambinp, lastpos, SEEK_SET) < 0)
846     goto seekerr;
847     while (n-- > 0) { /* load new contributed values */
848     if (!readambval(&avs, ambinp)) {
849     sprintf(errmsg, "ambient file \"%s\" corrupted",
850     ambfile);
851     error(WARNING, errmsg);
852     break;
853 greg 2.15 }
854 greg 2.129 avstore(&avs);
855     }
856     lastpos = newpos; /* update endpoint */
857     return(0);
858 greg 2.21 seekerr:
859     error(SYSTEM, "seek failed in ambsync");
860 greg 2.104 return(EOF); /* pro forma return */
861 greg 2.18 }