ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/ambient.c
Revision: 2.123
Committed: Fri Apr 5 01:10:26 2024 UTC (4 weeks ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 2.122: +18 -14 lines
Log Message:
fix: Improved tracking of reflected vs. transmitted rays for antimatter

File Contents

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