ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/ambient.c
Revision: 2.130
Committed: Thu Jan 23 18:44:20 2025 UTC (3 months, 1 week ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.129: +5 -4 lines
Log Message:
perf: Minor optimization in latest change; still some acceptable contention

File Contents

# Content
1 static const char RCSid[] = "$Id: ambient.c,v 2.129 2025/01/23 02:15:33 greg Exp $";
2 /*
3 * ambient.c - routines dealing with ambient (inter-reflected) component.
4 *
5 * Declarations of external symbols in ambient.h
6 */
7
8 #include "copyright.h"
9
10 #include <string.h>
11
12 #include "platform.h"
13 #include "ray.h"
14 #include "otypes.h"
15 #include "otspecial.h"
16 #include "resolu.h"
17 #include "ambient.h"
18 #include "random.h"
19 #include "pmapamb.h"
20
21 #ifndef OCTSCALE
22 #define OCTSCALE 1.0 /* ceil((valid rad.)/(cube size)) */
23 #endif
24
25 #ifndef MAXASET
26 #define MAXASET 4095 /* maximum number of elements in ambient set */
27 #endif
28 OBJECT ambset[MAXASET+1]={0}; /* ambient include/exclude set */
29
30 double maxarad; /* maximum ambient radius */
31 double minarad; /* minimum ambient radius */
32
33 static AMBTREE atrunk; /* our ambient trunk node */
34
35 static FILE *ambfp = NULL; /* ambient file pointer */
36 static int nunflshed; /* number of unflushed ambient values */
37 static FILE *ambinp = NULL; /* input pointer for ambient i/o */
38
39 static double avsum = 0.; /* computed ambient value sum (log) */
40 static unsigned int navsum = 0; /* number of values in avsum */
41 static unsigned int nambvals = 0; /* total number of indirect values */
42 static off_t lastpos = -1; /* last flush position */
43
44 #define AMBFLUSH (BUFSIZ/AMBVALSIZ)
45
46 #define AVSIZE (sizeof(AMBVAL)-sizeof(SCOLOR)+sizeof(COLORV)*NCSAMP)
47 #define newambval() (AMBVAL *)malloc(AVSIZE)
48
49 #define tfunc(x0, x, x1) (((x)-(x0))/((x1)-(x0)))
50
51 static void initambfile(int cre8);
52 static void avsave(AMBVAL *av);
53 static AMBVAL *avstore(AMBVAL *aval);
54 static AMBTREE *newambtree(void);
55 static void freeambtree(AMBTREE *atp);
56
57 typedef void unloadtf_t(AMBVAL *);
58 static unloadtf_t avinsert;
59 static unloadtf_t avfree;
60 static void unloadatree(AMBTREE *at, unloadtf_t *f);
61
62 static void sortambvals(void);
63
64 static int plugaleak(RAY *r, AMBVAL *ap, FVECT anorm, double ang);
65 static double sumambient(SCOLOR acol, RAY *r, FVECT rn, int al,
66 AMBTREE *at, FVECT c0, double s);
67 static int makeambient(SCOLOR acol, RAY *r, FVECT rn, int al);
68 static int extambient(SCOLOR cr, AMBVAL *ap, FVECT pv, FVECT nv,
69 FVECT uvw[3]);
70
71
72 void
73 setambres( /* set ambient resolution */
74 int ar
75 )
76 {
77 ambres = ar < 0 ? 0 : ar; /* may be done already */
78 /* set min & max radii */
79 if (ar <= 0) {
80 minarad = 0;
81 maxarad = thescene.cusize*0.2;
82 } else {
83 minarad = thescene.cusize / ar;
84 maxarad = 64.0 * minarad; /* heuristic */
85 if (maxarad > thescene.cusize*0.2)
86 maxarad = thescene.cusize*0.2;
87 }
88 if (minarad <= FTINY)
89 minarad = 10.0*FTINY;
90 if (maxarad <= minarad)
91 maxarad = 64.0 * minarad;
92 }
93
94
95 void
96 setambacc( /* set ambient accuracy */
97 double newa
98 )
99 {
100 static double olda; /* remember previous setting here */
101
102 newa *= (newa > 0);
103 if (fabs(newa - olda) >= .05*(newa + olda)) {
104 ambacc = newa;
105 if (ambacc > FTINY && nambvals > 0)
106 sortambvals(); /* rebuild tree */
107 }
108 }
109
110
111 void
112 setambient(void) /* initialize calculation */
113 {
114 int exists;
115 off_t flen;
116 AMBVAL amb;
117 /* make sure we're fresh */
118 ambdone();
119 /* init ambient limits */
120 setambres(ambres);
121 setambacc(ambacc);
122 if (ambfile == NULL || !ambfile[0])
123 return;
124 if (ambacc <= FTINY) {
125 sprintf(errmsg, "zero ambient accuracy so \"%s\" not opened",
126 ambfile);
127 error(WARNING, errmsg);
128 return;
129 }
130 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 lastpos = ftell(ambfp);
150 while (readambval(&amb, ambfp))
151 avstore(&amb); /* load what we can */
152 if (exists < 0) { /* read-only? */
153 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 return;
160 }
161 /* align file pointer */
162 lastpos += (off_t)nambvals*AMBVALSIZ;
163 flen = lseek(fileno(ambfp), 0, SEEK_END);
164 if (flen != lastpos) {
165 sprintf(errmsg,
166 "ignoring last %lu values in ambient file (corrupted)",
167 (unsigned long)((flen - lastpos)/AMBVALSIZ));
168 error(WARNING, errmsg);
169 /* fseek() not needed? */
170 fseeko(ambfp, lastpos, SEEK_SET);
171 ftruncate(fileno(ambfp), lastpos);
172 }
173 } else {
174 initambfile(1); /* else start new file */
175 fflush(ambfp);
176 lastpos = ftell(ambfp);
177 }
178 }
179
180
181 void
182 ambdone(void) /* close ambient file and free memory */
183 {
184 if (ambfp != NULL) { /* close ambient file */
185 fclose(ambfp); /* don't call ambsync() */
186 ambfp = NULL;
187 lastpos = -1;
188 if (ambinp != NULL) {
189 fclose(ambinp);
190 ambinp = NULL;
191 }
192 }
193 /* free ambient tree */
194 unloadatree(&atrunk, avfree);
195 /* reset state variables */
196 avsum = 0.;
197 navsum = 0;
198 nambvals = 0;
199 }
200
201
202 void
203 ambnotify( /* record new modifier */
204 OBJECT obj
205 )
206 {
207 static int hitlimit = 0;
208 OBJREC *o;
209 char **amblp;
210
211 if (obj == OVOID) { /* starting over */
212 ambset[0] = 0;
213 hitlimit = 0;
214 return;
215 }
216 o = objptr(obj);
217 if (hitlimit || !ismodifier(o->otype))
218 return;
219 for (amblp = amblist; *amblp != NULL; amblp++)
220 if (!strcmp(o->oname, *amblp)) {
221 if (ambset[0] >= MAXASET) {
222 error(WARNING, "too many modifiers in ambient list");
223 hitlimit++;
224 return; /* should this be fatal? */
225 }
226 insertelem(ambset, obj);
227 return;
228 }
229 }
230
231
232 void
233 multambient( /* compute ambient component & multiply by coef. */
234 SCOLOR aval,
235 RAY *r,
236 FVECT nrm
237 )
238 {
239 static double logAvgAbsorp = 1;
240 static int rdepth = 0; /* ambient recursion */
241 SCOLOR acol, caustic;
242 int i, ok;
243 double d, l;
244
245 /* 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 if (logAvgAbsorp > 0) /* exclude in -aw to avoid growth */
252 logAvgAbsorp = log(1.-AVGREFL);
253
254 /* 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 {/* 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 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 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 FVECT uvd[2];
278
279 if (sgn*rdot < 0.9999)
280 dgp = dgrad; /* compute rotational grad. */
281 copyscolor(acol, aval);
282 rdepth++;
283 ok = doambient(acol, r, r->rweight*sgn,
284 uvd, NULL, NULL, dgp, NULL);
285 rdepth--;
286 if (!ok)
287 goto dumbamb;
288 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 d += sgn*v1[i] * (dgp[0]*uvd[0][i] + dgp[1]*uvd[1][i]);
294 if (d >= 0.05)
295 scalescolor(acol, d);
296 }
297 copyscolor(aval, acol);
298
299 /* PMAP: add in caustic */
300 saddscolor(aval, caustic);
301 return;
302 }
303 /* interpolate ambient value */
304 scolorblack(acol);
305 d = sumambient(acol, r, nrm, rdepth,
306 &atrunk, thescene.cuorg, thescene.cusize);
307
308 if (d > FTINY) {
309 scalescolor(acol, 1.0/d);
310 smultscolor(aval, acol);
311
312 /* PMAP: add in caustic */
313 saddscolor(aval, caustic);
314 return;
315 }
316
317 rdepth++; /* need to cache new value */
318 ok = makeambient(acol, r, nrm, rdepth-1);
319 rdepth--;
320
321 if (ok) {
322 smultscolor(aval, acol); /* computed new value */
323
324 /* PMAP: add in caustic */
325 saddscolor(aval, caustic);
326 return;
327 }
328
329 dumbamb: /* return global value */
330 if ((ambvwt <= 0) | (navsum == 0)) {
331 smultcolor(aval, ambval);
332
333 /* PMAP: add in caustic */
334 saddscolor(aval, caustic);
335 return;
336 }
337
338 l = bright(ambval); /* average in computations */
339 if (l > FTINY) {
340 d = (log(l)*(double)ambvwt + avsum + logAvgAbsorp*navsum) /
341 (double)(ambvwt + navsum);
342 d = exp(d) / l;
343 scalescolor(aval, d);
344 smultcolor(aval, ambval); /* apply color of ambval */
345 } else {
346 d = exp( avsum/(double)navsum + logAvgAbsorp );
347 scalescolor(aval, d); /* neutral color */
348 }
349 }
350
351
352 /* 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 OBJREC *m = findmaterial(rtst.ro);
386 if (m != NULL && !istransp(m) && !isBSDFproxy(m) &&
387 (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 sumambient( /* get interpolated ambient value */
400 SCOLOR acol,
401 RAY *r,
402 FVECT rn,
403 int al,
404 AMBTREE *at,
405 FVECT c0,
406 double s
407 )
408 { /* initial limit is 10 degrees plus ambacc radians */
409 const double minangle = 10.0 * PI/180.;
410 const int sgn = 1 - 2*(DOT(r->ron,rn) < 0);
411 double maxangle = minangle + ambacc;
412 double wsum = 0.0;
413 FVECT ck0;
414 int i, j;
415 AMBVAL *av;
416
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 if ((wsum >= 0.05) & (s*ambacc > minarad))
435 return(wsum);
436 }
437 /* 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 /* sum this node */
441 for (av = at->alist; av != NULL; av = av->next) {
442 double u, v, d, delta_r2, delta_t2;
443 SCOLOR sct;
444 FVECT uvw[3];
445 /*
446 * Ambient level test
447 */
448 if (av->lvl > al || /* list sorted, so this works */
449 (av->lvl == al) & (av->weight < 0.9*r->rweight))
450 break;
451 /*
452 * Direction test using unperturbed normal
453 */
454 decodedir(uvw[2], av->ndir);
455 d = sgn * DOT(uvw[2], r->ron);
456 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 * Modified ray behind test
463 */
464 VSUB(ck0, r->rop, av->pos);
465 d = DOT(ck0, uvw[2]);
466 if (d < -minarad*ambacc)
467 continue;
468 d /= av->rad[0];
469 delta_t2 = d*d;
470 if (delta_t2 >= ambacc*ambacc)
471 continue;
472 /*
473 * Elliptical radii test based on Hessian
474 */
475 decodedir(uvw[0], av->udir);
476 VCROSS(uvw[1], uvw[2], uvw[0]);
477 d = (u = DOT(ck0, uvw[0])) / av->rad[0];
478 delta_t2 += d*d;
479 d = (v = DOT(ck0, uvw[1])) / av->rad[1];
480 delta_t2 += d*d;
481 if (delta_t2 >= ambacc*ambacc)
482 continue;
483 /*
484 * Test for potential light leak
485 */
486 if (av->corral && plugaleak(r, av, uvw[2], atan2a(v,u)))
487 continue;
488 /*
489 * Extrapolate value and compute final weight (hat function)
490 */
491 if (!extambient(sct, av, r->rop, rn, uvw))
492 continue;
493 d = tfunc(maxangle, sqrt(delta_r2), 0.0) *
494 tfunc(ambacc, sqrt(delta_t2), 0.0);
495 scalescolor(sct, d);
496 saddscolor(acol, sct);
497 wsum += d;
498 }
499 return(wsum);
500 }
501
502
503 static int
504 makeambient( /* make a new ambient value for storage */
505 SCOLOR acol,
506 RAY *r,
507 FVECT rn,
508 int al
509 )
510 {
511 int sgn = 1 - 2*(DOT(r->ron,rn) < 0);
512 AMBVAL amb;
513 FVECT uvw[3];
514 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 setscolor(acol, AVGREFL, AVGREFL, AVGREFL);
522 /* compute ambient */
523 i = doambient(acol, r, amb.weight*sgn,
524 uvw, amb.rad, amb.gpos, amb.gdir, &amb.corral);
525 scalescolor(acol, 1./AVGREFL); /* undo assumed reflectance */
526 if (i <= 0 || amb.rad[0] <= FTINY) /* no Hessian or zero radius */
527 return(i);
528 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 /* store value */
532 VCOPY(amb.pos, r->rop);
533 amb.ndir = encodedir(uvw[2]);
534 amb.udir = encodedir(uvw[0]);
535 amb.lvl = al;
536 copyscolor(amb.val, acol);
537 avsave(&amb); /* insert and save to file */
538 if (DOT(uvw[2],rn) < 0.9999) /* texture? */
539 extambient(acol, &amb, r->rop, rn, uvw);
540 return(1);
541 }
542
543
544 static int
545 extambient( /* extrapolate value at pv, nv */
546 SCOLOR scr,
547 AMBVAL *ap,
548 FVECT pv,
549 FVECT nv,
550 FVECT uvw[3]
551 )
552 {
553 const double min_d = 0.05;
554 const double max_d = 20.;
555 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 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 if (d < min_d) /* clamp min/max scaling */
575 d = min_d;
576 else if (d > max_d)
577 d = max_d;
578 copyscolor(scr, ap->val);
579 scalescolor(scr, d);
580 return(d > min_d);
581 }
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 while (s*(OCTSCALE/2) > av->rad[1]*ambacc) {
603 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 if ( ap->next->lvl > av->lvl ||
618 (ap->next->lvl == av->lvl) &
619 (ap->next->weight <= av->weight) )
620 break;
621 av->next = ap->next;
622 ap->next = (AMBVAL*)av;
623 at->alist = avh.next;
624 }
625
626
627 static void
628 initambfile( /* initialize ambient file */
629 int cre8
630 )
631 {
632 extern char *progname, *octname;
633 static char *mybuf = NULL;
634 int ntries = 3;
635
636 if (!AMBFLUSH)
637 error(INTERNAL, "BUFSIZ too small in initambfile");
638 SET_FILE_BINARY(ambfp);
639 if (mybuf == NULL)
640 mybuf = (char *)bmalloc(BUFSIZ);
641 setbuf(ambfp, mybuf);
642 nunflshed = 0;
643 retry:
644 if (cre8) { /* new file */
645 newheader("RADIANCE", ambfp);
646 fprintf(ambfp, "%s -av %g %g %g -aw %d -ab %d -aa %g ",
647 progname, colval(ambval,RED),
648 colval(ambval,GRN), colval(ambval,BLU),
649 ambvwt, ambounce, ambacc);
650 fprintf(ambfp, "-ad %d -as %d -ar %d ",
651 ambdiv, ambssamp, ambres);
652 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 fprintf(ambfp, "-cw %g %g -cs %d ", WLPART[3], WLPART[0], NCSAMP);
657 if (octname != NULL)
658 fputs(octname, ambfp);
659 fputc('\n', ambfp); /* end of command line, not header! */
660 fprintf(ambfp, "SOFTWARE= %s\n", VersionID);
661 fputnow(ambfp);
662 AMB_CNDX = CNDX; /* use current spectral sampling */
663 AMB_WLPART = WLPART;
664 fputwlsplit(WLPART, ambfp);
665 fputncomp(NCSAMP, ambfp);
666 fputformat(AMBFMT, ambfp);
667 fputc('\n', ambfp);
668 putambmagic(ambfp);
669 } else if (getheader(ambfp, amb_headline, NULL) < 0 || !hasambmagic(ambfp)) {
670 if (--ntries > 0 && ftell(ambfp) == 0) {
671 clearerr(ambfp);
672 sleep(2);
673 goto retry;
674 }
675 error(USER, "bad/incompatible ambient file");
676 }
677 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 }
688
689
690 static void
691 avsave( /* insert and save an ambient value */
692 AMBVAL *av
693 )
694 {
695 avstore(av);
696 if (ambfp == NULL)
697 return;
698 if (writambval(av, ambfp) < 0)
699 goto writerr;
700 if (++nunflshed >= AMBFLUSH)
701 if (ambsync() == EOF)
702 goto writerr;
703 return;
704 writerr:
705 error(SYSTEM, "error writing to ambient file");
706 }
707
708
709 static AMBVAL *
710 avstore( /* allocate memory and save aval */
711 AMBVAL *aval
712 )
713 {
714 AMBVAL *av;
715 double d;
716
717 if ((av = newambval()) == NULL)
718 error(SYSTEM, "out of memory in avstore");
719 memcpy(av, aval, AVSIZE); /* AVSIZE <= sizeof(AMBVAL) */
720 av->next = NULL;
721 nambvals++;
722 d = pbright(av->val);
723 if (d > FTINY) { /* add to log sum for averaging */
724 avsum += log(d);
725 navsum++;
726 }
727 avinsert(av); /* insert in our cache tree */
728 return(av);
729 }
730
731
732 #define ATALLOCSZ 512 /* #/8 trees to allocate at once */
733
734 static AMBTREE *atfreelist = NULL; /* free ambient tree structures */
735
736
737 static AMBTREE *
738 newambtree(void) /* allocate 8 ambient tree structs */
739 {
740 AMBTREE *atp, *upperlim;
741
742 if (atfreelist == NULL) { /* get more nodes */
743 atfreelist = (AMBTREE *)malloc(ATALLOCSZ*8*sizeof(AMBTREE));
744 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 memset(atp, 0, 8*sizeof(AMBTREE));
755 return(atp);
756 }
757
758
759 static void
760 freeambtree( /* free 8 ambient tree structs */
761 AMBTREE *atp
762 )
763 {
764 atp->kid = atfreelist;
765 atfreelist = atp;
766 }
767
768
769 static void
770 unloadatree( /* unload an ambient value tree */
771 AMBTREE *at,
772 unloadtf_t *f
773 )
774 {
775 AMBVAL *av;
776 int i;
777 /* transfer values at this node */
778 for (av = at->alist; av != NULL; av = at->alist) {
779 at->alist = av->next;
780 av->next = NULL;
781 (*f)(av);
782 }
783 if (at->kid == NULL)
784 return;
785 for (i = 0; i < 8; i++) /* transfer and free children */
786 unloadatree(at->kid+i, f);
787 freeambtree(at->kid);
788 at->kid = NULL;
789 }
790
791
792 static void
793 avfree(AMBVAL *av)
794 {
795 free(av);
796 }
797
798
799 static void
800 sortambvals(void) /* resort ambient values */
801 {
802 AMBTREE oldatrunk = atrunk;
803
804 atrunk.alist = NULL;
805 atrunk.kid = NULL;
806 unloadatree(&oldatrunk, avinsert);
807 }
808
809
810 int
811 ambsync(void) /* synchronize ambient file */
812 {
813 off_t newpos;
814 int n;
815 AMBVAL avs;
816
817 if (ambfp == NULL) /* no ambient file? */
818 return(0);
819
820 if (nunflshed > 0) { /* append new values? */
821 if (fflush(ambfp) < 0)
822 return(EOF);
823 newpos = lseek(fileno(ambfp), 0, SEEK_CUR);
824 } else
825 newpos = lseek(fileno(ambfp), 0, SEEK_END);
826
827 if (newpos < 0)
828 goto seekerr;
829 /* how many others added? */
830 n = (newpos - lastpos)/AMBVALSIZ - nunflshed;
831 nunflshed = 0;
832 if (n <= 0) { /* no one helping this time? */
833 lastpos = newpos;
834 return(0);
835 }
836 if (ambinp == NULL) { /* else need to open for input? */
837 ambinp = fopen(ambfile, "r");
838 if (ambinp == NULL) {
839 sprintf(errmsg, "cannot reopen ambient file \"%s\"",
840 ambfile);
841 error(SYSTEM, errmsg);
842 }
843 SET_FILE_BINARY(ambinp);
844 }
845 /* read from last endpoint */
846 if (fseeko(ambinp, lastpos, SEEK_SET) < 0)
847 goto seekerr;
848 while (n-- > 0) { /* load new contributed values */
849 if (!readambval(&avs, ambinp)) {
850 sprintf(errmsg, "ambient file \"%s\" corrupted",
851 ambfile);
852 error(WARNING, errmsg);
853 break;
854 }
855 avstore(&avs);
856 }
857 lastpos = newpos; /* update endpoint */
858 return(0);
859 seekerr:
860 error(SYSTEM, "seek failed in ambsync");
861 return(EOF); /* pro forma return */
862 }