ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/rt/ambient.c
(Generate patch)

Comparing ray/src/rt/ambient.c (file contents):
Revision 2.34 by greg, Mon Nov 6 12:03:13 1995 UTC vs.
Revision 2.103 by greg, Tue Nov 1 20:39:39 2016 UTC

# Line 1 | Line 1
1 < /* Copyright (c) 1995 Regents of the University of California */
2 <
3 < #ifndef lint
4 < static char SCCSid[] = "$SunId$ LBL";
5 < #endif
6 <
1 > static const char       RCSid[] = "$Id$";
2   /*
3   *  ambient.c - routines dealing with ambient (inter-reflected) component.
4 + *
5 + *  Declarations of external symbols in ambient.h
6   */
7  
8 < #include  "ray.h"
8 > #include "copyright.h"
9  
10 < #include  "octree.h"
10 > #include <string.h>
11  
12 + #include  "platform.h"
13 + #include  "ray.h"
14   #include  "otypes.h"
15 <
15 > #include  "resolu.h"
16   #include  "ambient.h"
18
17   #include  "random.h"
18 + #include  "pmapamb.h"
19  
20   #ifndef  OCTSCALE
21   #define  OCTSCALE       1.0     /* ceil((valid rad.)/(cube size)) */
22   #endif
24 #ifndef  AMBVWT
25 #define  AMBVWT         250     /* relative ambient value weight (# calcs) */
26 #endif
23  
28 typedef struct ambtree {
29        AMBVAL  *alist;         /* ambient value list */
30        struct ambtree  *kid;   /* 8 child nodes */
31 }  AMBTREE;                     /* ambient octree */
32
33 extern CUBE  thescene;          /* contains space boundaries */
34
24   extern char  *shm_boundary;     /* memory sharing boundary */
25  
26 < #define  MAXASET        511     /* maximum number of elements in ambient set */
26 > #ifndef  MAXASET
27 > #define  MAXASET        4095    /* maximum number of elements in ambient set */
28 > #endif
29   OBJECT  ambset[MAXASET+1]={0};  /* ambient include/exclude set */
30  
31   double  maxarad;                /* maximum ambient radius */
# Line 46 | Line 37 | static FILE  *ambfp = NULL;    /* ambient file pointer */
37   static int  nunflshed = 0;      /* number of unflushed ambient values */
38  
39   #ifndef SORT_THRESH
40 < #ifdef BIGMEM
41 < #define SORT_THRESH     ((9L<<20)/sizeof(AMBVAL))
40 > #ifdef SMLMEM
41 > #define SORT_THRESH     ((16L<<20)/sizeof(AMBVAL))
42   #else
43 < #define SORT_THRESH     ((3L<<20)/sizeof(AMBVAL))
43 > #define SORT_THRESH     ((64L<<20)/sizeof(AMBVAL))
44   #endif
45   #endif
46   #ifndef SORT_INTVL
47 < #define SORT_INTVL      (SORT_THRESH*256)
47 > #define SORT_INTVL      (SORT_THRESH<<1)
48   #endif
49   #ifndef MAX_SORT_INTVL
50 < #define MAX_SORT_INTVL  (SORT_INTVL<<4)
50 > #define MAX_SORT_INTVL  (SORT_INTVL<<6)
51   #endif
52  
53 < static COLOR  avsum = BLKCOLOR;         /* computed ambient value sum */
54 < static unsigned int  nambvals = 0;      /* number of computed ambient values */
53 >
54 > static double  avsum = 0.;              /* computed ambient value sum (log) */
55 > static unsigned int  navsum = 0;        /* number of values in avsum */
56 > static unsigned int  nambvals = 0;      /* total number of indirect values */
57 > static unsigned int  nambshare = 0;     /* number of values from file */
58   static unsigned long  ambclock = 0;     /* ambient access clock */
59   static unsigned long  lastsort = 0;     /* time of last value sort */
60   static long  sortintvl = SORT_INTVL;    /* time until next sort */
61 + static FILE  *ambinp = NULL;            /* auxiliary file for input */
62 + static long  lastpos = -1;              /* last flush position */
63  
64   #define MAXACLOCK       (1L<<30)        /* clock turnover value */
65          /*
66           * Track access times unless we are sharing ambient values
67           * through memory on a multiprocessor, when we want to avoid
68 <         * claiming our own memory (copy on write).
68 >         * claiming our own memory (copy on write).  Go ahead anyway
69 >         * if more than two thirds of our values are unshared.
70 >         * Compile with -Dtracktime=0 to turn this code off.
71           */
72 < #define tracktime       (shm_boundary == NULL || ambfp == NULL)
72 > #ifndef tracktime
73 > #define tracktime       (shm_boundary == NULL || nambvals > 3*nambshare)
74 > #endif
75  
76   #define  AMBFLUSH       (BUFSIZ/AMBVALSIZ)
77  
78 < #define  newambval()    (AMBVAL *)bmalloc(sizeof(AMBVAL))
78 > #define  newambval()    (AMBVAL *)malloc(sizeof(AMBVAL))
79  
80 < extern long  ftell(), lseek();
81 < static int  initambfile(), avsave(), avinsert(), sortambvals();
82 < static AMBVAL  *avstore();
80 > static void initambfile(int creat);
81 > static void avsave(AMBVAL *av);
82 > static AMBVAL *avstore(AMBVAL  *aval);
83 > static AMBTREE *newambtree(void);
84 > static void freeambtree(AMBTREE  *atp);
85 >
86 > typedef void unloadtf_t(AMBVAL *);
87 > static unloadtf_t avinsert;
88 > static unloadtf_t av2list;
89 > static unloadtf_t avfree;
90 > static void unloadatree(AMBTREE  *at, unloadtf_t *f);
91 >
92 > static int aposcmp(const void *avp1, const void *avp2);
93 > static int avlmemi(AMBVAL *avaddr);
94 > static void sortambvals(int always);
95 >
96   #ifdef  F_SETLKW
97 < static  aflock();
97 > static void aflock(int  typ);
98   #endif
99  
100  
101 < setambres(ar)                           /* set ambient resolution */
102 < int  ar;
101 > void
102 > setambres(                              /* set ambient resolution */
103 >        int  ar
104 > )
105   {
106          ambres = ar < 0 ? 0 : ar;               /* may be done already */
107                                                  /* set min & max radii */
108          if (ar <= 0) {
109                  minarad = 0;
110 <                maxarad = thescene.cusize / 2.0;
110 >                maxarad = thescene.cusize*0.2;
111          } else {
112                  minarad = thescene.cusize / ar;
113 <                maxarad = 64 * minarad;                 /* heuristic */
114 <                if (maxarad > thescene.cusize / 2.0)
115 <                        maxarad = thescene.cusize / 2.0;
113 >                maxarad = 64.0 * minarad;               /* heuristic */
114 >                if (maxarad > thescene.cusize*0.2)
115 >                        maxarad = thescene.cusize*0.2;
116          }
117          if (minarad <= FTINY)
118 <                minarad = 10*FTINY;
118 >                minarad = 10.0*FTINY;
119          if (maxarad <= minarad)
120 <                maxarad = 64 * minarad;
120 >                maxarad = 64.0 * minarad;
121   }
122  
123  
124 < setambacc(newa)                         /* set ambient accuracy */
125 < double  newa;
124 > void
125 > setambacc(                              /* set ambient accuracy */
126 >        double  newa
127 > )
128   {
129 <        double  ambdiff;
130 <
131 <        if (newa < 0.0)
132 <                newa = 0.0;
133 <        ambdiff = fabs(newa - ambacc);
134 <        if (ambdiff >= .01 && (ambacc = newa) > FTINY && nambvals > 0)
135 <                sortambvals(1);                 /* rebuild tree */
129 >        static double   olda;           /* remember previous setting here */
130 >        
131 >        newa *= (newa > 0);
132 >        if (fabs(newa - olda) >= .05*(newa + olda)) {
133 >                ambacc = newa;
134 >                if (nambvals > 0)
135 >                        sortambvals(1);         /* rebuild tree */
136 >        }
137   }
138  
139  
140 < setambient(afile)                       /* initialize calculation */
141 < char  *afile;
140 > void
141 > setambient(void)                                /* initialize calculation */
142   {
143 <        long  headlen;
143 >        int     readonly = 0;
144 >        long    flen;
145          AMBVAL  amb;
146 +                                                /* make sure we're fresh */
147 +        ambdone();
148                                                  /* init ambient limits */
149          setambres(ambres);
150          setambacc(ambacc);
151 <        if (afile == NULL)
151 >        if (ambfile == NULL || !ambfile[0])
152                  return;
153          if (ambacc <= FTINY) {
154                  sprintf(errmsg, "zero ambient accuracy so \"%s\" not opened",
155 <                                afile);
155 >                                ambfile);
156                  error(WARNING, errmsg);
157                  return;
158          }
159                                                  /* open ambient file */
160 <        if ((ambfp = fopen(afile, "r+")) != NULL) {
161 <                initambfile(0);
162 <                headlen = ftell(ambfp);
160 >        if ((ambfp = fopen(ambfile, "r+")) == NULL)
161 >                readonly = (ambfp = fopen(ambfile, "r")) != NULL;
162 >        if (ambfp != NULL) {
163 >                initambfile(0);                 /* file exists */
164 >                lastpos = ftell(ambfp);
165                  while (readambval(&amb, ambfp))
166 <                        avinsert(avstore(&amb));
167 <                                                /* align */
168 <                fseek(ambfp, -((ftell(ambfp)-headlen)%AMBVALSIZ), 1);
169 <        } else if ((ambfp = fopen(afile, "w+")) != NULL)
170 <                initambfile(1);
171 <        else {
172 <                sprintf(errmsg, "cannot open ambient file \"%s\"", afile);
166 >                        avstore(&amb);
167 >                nambshare = nambvals;           /* share loaded values */
168 >                if (readonly) {
169 >                        sprintf(errmsg,
170 >                                "loaded %u values from read-only ambient file",
171 >                                        nambvals);
172 >                        error(WARNING, errmsg);
173 >                        fclose(ambfp);          /* close file so no writes */
174 >                        ambfp = NULL;
175 >                        return;                 /* avoid ambsync() */
176 >                }
177 >                                                /* align file pointer */
178 >                lastpos += (long)nambvals*AMBVALSIZ;
179 >                flen = lseek(fileno(ambfp), (off_t)0, SEEK_END);
180 >                if (flen != lastpos) {
181 >                        sprintf(errmsg,
182 >                        "ignoring last %ld values in ambient file (corrupted)",
183 >                                        (flen - lastpos)/AMBVALSIZ);
184 >                        error(WARNING, errmsg);
185 >                        fseek(ambfp, lastpos, SEEK_SET);
186 >                        ftruncate(fileno(ambfp), (off_t)lastpos);
187 >                }
188 >        } else if ((ambfp = fopen(ambfile, "w+")) != NULL) {
189 >                initambfile(1);                 /* else create new file */
190 >                fflush(ambfp);
191 >                lastpos = ftell(ambfp);
192 >        } else {
193 >                sprintf(errmsg, "cannot open ambient file \"%s\"", ambfile);
194                  error(SYSTEM, errmsg);
195          }
196 <        nunflshed++;    /* lie */
197 <        ambsync();
196 > #ifdef  F_SETLKW
197 >        aflock(F_UNLCK);                        /* release file */
198 > #endif
199   }
200  
201  
202 < ambnotify(obj)                  /* record new modifier */
203 < OBJECT  obj;
202 > void
203 > ambdone(void)                   /* close ambient file and free memory */
204   {
205 +        if (ambfp != NULL) {            /* close ambient file */
206 +                ambsync();
207 +                fclose(ambfp);
208 +                ambfp = NULL;
209 +                if (ambinp != NULL) {  
210 +                        fclose(ambinp);
211 +                        ambinp = NULL;
212 +                }
213 +                lastpos = -1;
214 +        }
215 +                                        /* free ambient tree */
216 +        unloadatree(&atrunk, avfree);
217 +                                        /* reset state variables */
218 +        avsum = 0.;
219 +        navsum = 0;
220 +        nambvals = 0;
221 +        nambshare = 0;
222 +        ambclock = 0;
223 +        lastsort = 0;
224 +        sortintvl = SORT_INTVL;
225 + }
226 +
227 +
228 + void
229 + ambnotify(                      /* record new modifier */
230 +        OBJECT  obj
231 + )
232 + {
233          static int  hitlimit = 0;
234 <        register OBJREC  *o = objptr(obj);
235 <        register char  **amblp;
234 >        OBJREC   *o;
235 >        char  **amblp;
236  
237 +        if (obj == OVOID) {             /* starting over */
238 +                ambset[0] = 0;
239 +                hitlimit = 0;
240 +                return;
241 +        }
242 +        o = objptr(obj);
243          if (hitlimit || !ismodifier(o->otype))
244                  return;
245          for (amblp = amblist; *amblp != NULL; amblp++)
# Line 175 | Line 254 | OBJECT obj;
254                  }
255   }
256  
257 + /************ THE FOLLOWING ROUTINES DIFFER BETWEEN NEW & OLD ***************/
258  
259 < ambient(acol, r, nrm)           /* compute ambient component for ray */
260 < COLOR  acol;
261 < register RAY  *r;
262 < FVECT  nrm;
259 > #ifndef OLDAMB
260 >
261 > #define tfunc(lwr, x, upr)      (((x)-(lwr))/((upr)-(lwr)))
262 >
263 > static int      plugaleak(RAY *r, AMBVAL *ap, FVECT anorm, double ang);
264 > static double   sumambient(COLOR acol, RAY *r, FVECT rn, int al,
265 >                                AMBTREE *at, FVECT c0, double s);
266 > static int      makeambient(COLOR acol, RAY *r, FVECT rn, int al);
267 > static int      extambient(COLOR cr, AMBVAL *ap, FVECT pv, FVECT nv,
268 >                                FVECT uvw[3]);
269 >
270 > void
271 > multambient(            /* compute ambient component & multiply by coef. */
272 >        COLOR  aval,
273 >        RAY  *r,
274 >        FVECT  nrm
275 > )
276   {
277          static int  rdepth = 0;                 /* ambient recursion */
278 <        double  d;
278 >        COLOR   acol, caustic;
279 >        int     i, ok;
280 >        double  d, l;
281  
282 +        /* PMAP: Factor in ambient from photon map, if enabled and ray is
283 +         * ambient. Return as all ambient components accounted for, else
284 +         * continue. */
285 +        if (ambPmap(aval, r, rdepth))
286 +                return;
287 +
288 +        /* PMAP: Factor in specular-diffuse ambient (caustics) from photon
289 +         * map, if enabled and ray is primary, else caustic is zero.  Continue
290 +         * with RADIANCE ambient calculation */
291 +        copycolor(caustic, aval);
292 +        ambPmapCaustic(caustic, r, rdepth);
293 +        
294          if (ambdiv <= 0)                        /* no ambient calculation */
295                  goto dumbamb;
296                                                  /* check number of bounces */
# Line 195 | Line 302 | FVECT  nrm;
302                  goto dumbamb;
303  
304          if (ambacc <= FTINY) {                  /* no ambient storage */
305 +                FVECT   uvd[2];
306 +                float   dgrad[2], *dgp = NULL;
307 +
308 +                if (nrm != r->ron && DOT(nrm,r->ron) < 0.9999)
309 +                        dgp = dgrad;            /* compute rotational grad. */
310 +                copycolor(acol, aval);
311                  rdepth++;
312 +                ok = doambient(acol, r, r->rweight,
313 +                                uvd, NULL, NULL, dgp, NULL);
314 +                rdepth--;
315 +                if (!ok)
316 +                        goto dumbamb;
317 +                if ((ok > 0) & (dgp != NULL)) { /* apply texture */
318 +                        FVECT   v1;
319 +                        VCROSS(v1, r->ron, nrm);
320 +                        d = 1.0;
321 +                        for (i = 3; i--; )
322 +                                d += v1[i] * (dgp[0]*uvd[0][i] + dgp[1]*uvd[1][i]);
323 +                        if (d >= 0.05)
324 +                                scalecolor(acol, d);
325 +                }
326 +                copycolor(aval, acol);
327 +
328 +                /* PMAP: add in caustic */
329 +                addcolor(aval, caustic);
330 +                return;
331 +        }
332 +
333 +        if (tracktime)                          /* sort to minimize thrashing */
334 +                sortambvals(0);
335 +                                                /* interpolate ambient value */
336 +        setcolor(acol, 0.0, 0.0, 0.0);
337 +        d = sumambient(acol, r, nrm, rdepth,
338 +                        &atrunk, thescene.cuorg, thescene.cusize);
339 +                        
340 +        if (d > FTINY) {
341 +                d = 1.0/d;
342 +                scalecolor(acol, d);
343 +                multcolor(aval, acol);
344 +
345 +                /* PMAP: add in caustic */
346 +                addcolor(aval, caustic);
347 +                return;
348 +        }
349 +        
350 +        rdepth++;                               /* need to cache new value */
351 +        ok = makeambient(acol, r, nrm, rdepth-1);
352 +        rdepth--;
353 +        
354 +        if (ok) {
355 +                multcolor(aval, acol);          /* computed new value */
356 +
357 +                /* PMAP: add in caustic */
358 +                addcolor(aval, caustic);
359 +                return;
360 +        }
361 +        
362 + dumbamb:                                        /* return global value */
363 +        if ((ambvwt <= 0) | (navsum == 0)) {
364 +                multcolor(aval, ambval);
365 +                
366 +                /* PMAP: add in caustic */
367 +                addcolor(aval, caustic);
368 +                return;
369 +        }
370 +        
371 +        l = bright(ambval);                     /* average in computations */  
372 +        if (l > FTINY) {
373 +                d = (log(l)*(double)ambvwt + avsum) /
374 +                                (double)(ambvwt + navsum);
375 +                d = exp(d) / l;
376 +                scalecolor(aval, d);
377 +                multcolor(aval, ambval);        /* apply color of ambval */
378 +        } else {
379 +                d = exp( avsum / (double)navsum );
380 +                scalecolor(aval, d);            /* neutral color */
381 +        }
382 + }
383 +
384 +
385 + /* Plug a potential leak where ambient cache value is occluded */
386 + static int
387 + plugaleak(RAY *r, AMBVAL *ap, FVECT anorm, double ang)
388 + {
389 +        const double    cost70sq = 0.1169778;   /* cos(70deg)^2 */
390 +        RAY             rtst;
391 +        FVECT           vdif;
392 +        double          normdot, ndotd, nadotd;
393 +        double          a, b, c, t[2];
394 +
395 +        ang += 2.*PI*(ang < 0);                 /* check direction flags */
396 +        if ( !(ap->corral>>(int)(ang*(16./PI)) & 1) )
397 +                return(0);
398 +        /*
399 +         * Generate test ray, targeting 20 degrees above sample point plane
400 +         * along surface normal from cache position.  This should be high
401 +         * enough to miss local geometry we don't really care about.
402 +         */
403 +        VSUB(vdif, ap->pos, r->rop);
404 +        normdot = DOT(anorm, r->ron);
405 +        ndotd = DOT(vdif, r->ron);
406 +        nadotd = DOT(vdif, anorm);
407 +        a = normdot*normdot - cost70sq;
408 +        b = 2.0*(normdot*ndotd - nadotd*cost70sq);
409 +        c = ndotd*ndotd - DOT(vdif,vdif)*cost70sq;
410 +        if (quadratic(t, a, b, c) != 2)
411 +                return(1);                      /* should rarely happen */
412 +        if (t[1] <= FTINY)
413 +                return(0);                      /* should fail behind test */
414 +        rayorigin(&rtst, SHADOW, r, NULL);
415 +        VSUM(rtst.rdir, vdif, anorm, t[1]);     /* further dist. > plane */
416 +        rtst.rmax = normalize(rtst.rdir);       /* short ray test */
417 +        while (localhit(&rtst, &thescene)) {    /* check for occluder */
418 +                if (rtst.ro->omod != OVOID &&
419 +                                (rtst.clipset == NULL ||
420 +                                        !inset(rtst.clipset, rtst.ro->omod)))
421 +                        return(1);              /* plug light leak */
422 +                VCOPY(rtst.rorg, rtst.rop);     /* skip invisible surface */
423 +                rtst.rmax -= rtst.rot;
424 +                rayclear(&rtst);
425 +        }
426 +        return(0);                              /* seems we're OK */
427 + }
428 +
429 +
430 + static double
431 + sumambient(             /* get interpolated ambient value */
432 +        COLOR  acol,
433 +        RAY  *r,
434 +        FVECT  rn,
435 +        int  al,
436 +        AMBTREE  *at,
437 +        FVECT  c0,
438 +        double  s
439 + )
440 + {                       /* initial limit is 10 degrees plus ambacc radians */
441 +        const double    minangle = 10.0 * PI/180.;
442 +        double          maxangle = minangle + ambacc;
443 +        double          wsum = 0.0;
444 +        FVECT           ck0;
445 +        int             i, j;
446 +        AMBVAL          *av;
447 +
448 +        if (at->kid != NULL) {          /* sum children first */                                
449 +                s *= 0.5;
450 +                for (i = 0; i < 8; i++) {
451 +                        for (j = 0; j < 3; j++) {
452 +                                ck0[j] = c0[j];
453 +                                if (1<<j & i)
454 +                                        ck0[j] += s;
455 +                                if (r->rop[j] < ck0[j] - OCTSCALE*s)
456 +                                        break;
457 +                                if (r->rop[j] > ck0[j] + (1.0+OCTSCALE)*s)
458 +                                        break;
459 +                        }
460 +                        if (j == 3)
461 +                                wsum += sumambient(acol, r, rn, al,
462 +                                                        at->kid+i, ck0, s);
463 +                }
464 +                                        /* good enough? */
465 +                if (wsum >= 0.05 && s > minarad*10.0)
466 +                        return(wsum);
467 +        }
468 +                                        /* adjust maximum angle */
469 +        if (at->alist != NULL && (at->alist->lvl <= al) & (r->rweight < 0.6))
470 +                maxangle = (maxangle - PI/2.)*pow(r->rweight,0.13) + PI/2.;
471 +                                        /* sum this node */
472 +        for (av = at->alist; av != NULL; av = av->next) {
473 +                double  u, v, d, delta_r2, delta_t2;
474 +                COLOR   ct;
475 +                FVECT   uvw[3];
476 +                                        /* record access */
477 +                if (tracktime)
478 +                        av->latick = ambclock;
479 +                /*
480 +                 *  Ambient level test
481 +                 */
482 +                if (av->lvl > al ||     /* list sorted, so this works */
483 +                                (av->lvl == al) & (av->weight < 0.9*r->rweight))
484 +                        break;
485 +                /*
486 +                 *  Direction test using unperturbed normal
487 +                 */
488 +                decodedir(uvw[2], av->ndir);
489 +                d = DOT(uvw[2], r->ron);
490 +                if (d <= 0.0)           /* >= 90 degrees */
491 +                        continue;
492 +                delta_r2 = 2.0 - 2.0*d; /* approx. radians^2 */
493 +                if (delta_r2 >= maxangle*maxangle)
494 +                        continue;
495 +                /*
496 +                 *  Modified ray behind test
497 +                 */
498 +                VSUB(ck0, r->rop, av->pos);
499 +                d = DOT(ck0, uvw[2]);
500 +                if (d < -minarad*ambacc-.001)
501 +                        continue;
502 +                d /= av->rad[0];
503 +                delta_t2 = d*d;
504 +                if (delta_t2 >= ambacc*ambacc)
505 +                        continue;
506 +                /*
507 +                 *  Elliptical radii test based on Hessian
508 +                 */
509 +                decodedir(uvw[0], av->udir);
510 +                VCROSS(uvw[1], uvw[2], uvw[0]);
511 +                d = (u = DOT(ck0, uvw[0])) / av->rad[0];
512 +                delta_t2 += d*d;
513 +                d = (v = DOT(ck0, uvw[1])) / av->rad[1];
514 +                delta_t2 += d*d;
515 +                if (delta_t2 >= ambacc*ambacc)
516 +                        continue;
517 +                /*
518 +                 *  Test for potential light leak
519 +                 */
520 +                if (av->corral && plugaleak(r, av, uvw[2], atan2a(v,u)))
521 +                        continue;
522 +                /*
523 +                 *  Extrapolate value and compute final weight (hat function)
524 +                 */
525 +                if (!extambient(ct, av, r->rop, rn, uvw))
526 +                        continue;
527 +                d = tfunc(maxangle, sqrt(delta_r2), 0.0) *
528 +                        tfunc(ambacc, sqrt(delta_t2), 0.0);
529 +                scalecolor(ct, d);
530 +                addcolor(acol, ct);
531 +                wsum += d;
532 +        }
533 +        return(wsum);
534 + }
535 +
536 +
537 + static int
538 + makeambient(            /* make a new ambient value for storage */
539 +        COLOR  acol,
540 +        RAY  *r,
541 +        FVECT  rn,
542 +        int  al
543 + )
544 + {
545 +        AMBVAL  amb;
546 +        FVECT   uvw[3];
547 +        int     i;
548 +
549 +        amb.weight = 1.0;                       /* compute weight */
550 +        for (i = al; i-- > 0; )
551 +                amb.weight *= AVGREFL;
552 +        if (r->rweight < 0.1*amb.weight)        /* heuristic override */
553 +                amb.weight = 1.25*r->rweight;
554 +        setcolor(acol, AVGREFL, AVGREFL, AVGREFL);
555 +                                                /* compute ambient */
556 +        i = doambient(acol, r, amb.weight,
557 +                        uvw, amb.rad, amb.gpos, amb.gdir, &amb.corral);
558 +        scalecolor(acol, 1./AVGREFL);           /* undo assumed reflectance */
559 +        if (i <= 0 || amb.rad[0] <= FTINY)      /* no Hessian or zero radius */
560 +                return(i);
561 +                                                /* store value */
562 +        VCOPY(amb.pos, r->rop);
563 +        amb.ndir = encodedir(r->ron);
564 +        amb.udir = encodedir(uvw[0]);
565 +        amb.lvl = al;
566 +        copycolor(amb.val, acol);
567 +                                                /* insert into tree */
568 +        avsave(&amb);                           /* and save to file */
569 +        if (rn != r->ron) {                     /* texture */
570 +                VCOPY(uvw[2], r->ron);
571 +                extambient(acol, &amb, r->rop, rn, uvw);
572 +        }
573 +        return(1);
574 + }
575 +
576 +
577 + static int
578 + extambient(             /* extrapolate value at pv, nv */
579 +        COLOR  cr,
580 +        AMBVAL   *ap,
581 +        FVECT  pv,
582 +        FVECT  nv,
583 +        FVECT  uvw[3]
584 + )
585 + {
586 +        const double    min_d = 0.05;
587 +        static FVECT    my_uvw[3];
588 +        FVECT           v1;
589 +        int             i;
590 +        double          d = 1.0;        /* zeroeth order */
591 +
592 +        if (uvw == NULL) {              /* need local coordinates? */
593 +                decodedir(my_uvw[2], ap->ndir);
594 +                decodedir(my_uvw[0], ap->udir);
595 +                VCROSS(my_uvw[1], my_uvw[2], my_uvw[0]);
596 +                uvw = my_uvw;
597 +        }
598 +        for (i = 3; i--; )              /* gradient due to translation */
599 +                d += (pv[i] - ap->pos[i]) *
600 +                        (ap->gpos[0]*uvw[0][i] + ap->gpos[1]*uvw[1][i]);
601 +
602 +        VCROSS(v1, uvw[2], nv);         /* gradient due to rotation */
603 +        for (i = 3; i--; )
604 +                d += v1[i] * (ap->gdir[0]*uvw[0][i] + ap->gdir[1]*uvw[1][i]);
605 +        
606 +        if (d < min_d)                  /* should not use if we can avoid it */
607 +                d = min_d;
608 +        copycolor(cr, ap->val);
609 +        scalecolor(cr, d);
610 +        return(d > min_d);
611 + }
612 +
613 +
614 + static void
615 + avinsert(                               /* insert ambient value in our tree */
616 +        AMBVAL *av
617 + )
618 + {
619 +        AMBTREE  *at;
620 +        AMBVAL  *ap;
621 +        AMBVAL  avh;
622 +        FVECT  ck0;
623 +        double  s;
624 +        int  branch;
625 +        int  i;
626 +
627 +        if (av->rad[0] <= FTINY)
628 +                error(CONSISTENCY, "zero ambient radius in avinsert");
629 +        at = &atrunk;
630 +        VCOPY(ck0, thescene.cuorg);
631 +        s = thescene.cusize;
632 +        while (s*(OCTSCALE/2) > av->rad[1]*ambacc) {
633 +                if (at->kid == NULL)
634 +                        if ((at->kid = newambtree()) == NULL)
635 +                                error(SYSTEM, "out of memory in avinsert");
636 +                s *= 0.5;
637 +                branch = 0;
638 +                for (i = 0; i < 3; i++)
639 +                        if (av->pos[i] > ck0[i] + s) {
640 +                                ck0[i] += s;
641 +                                branch |= 1 << i;
642 +                        }
643 +                at = at->kid + branch;
644 +        }
645 +        avh.next = at->alist;           /* order by increasing level */
646 +        for (ap = &avh; ap->next != NULL; ap = ap->next)
647 +                if ( ap->next->lvl > av->lvl ||
648 +                                (ap->next->lvl == av->lvl) &
649 +                                (ap->next->weight <= av->weight) )
650 +                        break;
651 +        av->next = ap->next;
652 +        ap->next = (AMBVAL*)av;
653 +        at->alist = avh.next;
654 + }
655 +
656 +
657 + #else /* ! NEWAMB */
658 +
659 + static double   sumambient(COLOR acol, RAY *r, FVECT rn, int al,
660 +                                AMBTREE *at, FVECT c0, double s);
661 + static double   makeambient(COLOR acol, RAY *r, FVECT rn, int al);
662 + static void     extambient(COLOR cr, AMBVAL *ap, FVECT pv, FVECT nv);
663 +
664 +
665 + void
666 + multambient(            /* compute ambient component & multiply by coef. */
667 +        COLOR  aval,
668 +        RAY  *r,
669 +        FVECT  nrm
670 + )
671 + {
672 +        static int  rdepth = 0;                 /* ambient recursion */
673 +        COLOR   acol, caustic;
674 +        double  d, l;
675 +
676 +        /* PMAP: Factor in ambient from global photon map (if enabled) and return
677 +         * as all ambient components accounted for */
678 +        if (ambPmap(aval, r, rdepth))
679 +                return;
680 +
681 +        /* PMAP: Otherwise factor in ambient from caustic photon map
682 +         * (ambPmapCaustic() returns zero if caustic photons disabled) and
683 +         * continue with RADIANCE ambient calculation */
684 +        copycolor(caustic, aval);
685 +        ambPmapCaustic(caustic, r, rdepth);
686 +        
687 +        if (ambdiv <= 0)                        /* no ambient calculation */
688 +                goto dumbamb;
689 +                                                /* check number of bounces */
690 +        if (rdepth >= ambounce)
691 +                goto dumbamb;
692 +                                                /* check ambient list */
693 +        if (ambincl != -1 && r->ro != NULL &&
694 +                        ambincl != inset(ambset, r->ro->omod))
695 +                goto dumbamb;
696 +
697 +        if (ambacc <= FTINY) {                  /* no ambient storage */
698 +                copycolor(acol, aval);
699 +                rdepth++;
700                  d = doambient(acol, r, r->rweight, NULL, NULL);
701                  rdepth--;
702                  if (d <= FTINY)
703                          goto dumbamb;
704 +                copycolor(aval, acol);          
705 +        
706 +           /* PMAP: add in caustic */
707 +                addcolor(aval, caustic);        
708                  return;
709          }
710 <                                                /* resort memory? */
711 <        sortambvals(0);
712 <                                                /* get ambient value */
710 >
711 >        if (tracktime)                          /* sort to minimize thrashing */
712 >                sortambvals(0);
713 >                                                /* interpolate ambient value */
714          setcolor(acol, 0.0, 0.0, 0.0);
715          d = sumambient(acol, r, nrm, rdepth,
716                          &atrunk, thescene.cuorg, thescene.cusize);
717 +                        
718          if (d > FTINY) {
719 <                scalecolor(acol, 1.0/d);
719 >                d = 1.0/d;
720 >                scalecolor(acol, d);
721 >                multcolor(aval, acol);
722 >                
723 >                /* PMAP: add in caustic */
724 >                addcolor(aval, caustic);        
725                  return;
726          }
727 +        
728          rdepth++;                               /* need to cache new value */
729          d = makeambient(acol, r, nrm, rdepth-1);
730          rdepth--;
731 <        if (d > FTINY)
731 >        
732 >        if (d > FTINY) {
733 >                multcolor(aval, acol);          /* got new value */
734 >
735 >                /* PMAP: add in caustic */
736 >                addcolor(aval, caustic);                        
737                  return;
738 +        }
739 +        
740   dumbamb:                                        /* return global value */
741 <        copycolor(acol, ambval);
742 < #if  AMBVWT
743 <        if (nambvals == 0)
741 >        if ((ambvwt <= 0) | (navsum == 0)) {
742 >                multcolor(aval, ambval);
743 >
744 >                /* PMAP: add in caustic */
745 >                addcolor(aval, caustic);        
746                  return;
747 <        scalecolor(acol, (double)AMBVWT);
748 <        addcolor(acol, avsum);                  /* average in computations */
749 <        d = 1.0/(AMBVWT+nambvals);
750 <        scalecolor(acol, d);
751 < #endif
747 >        }
748 >        
749 >        l = bright(ambval);                     /* average in computations */
750 >        if (l > FTINY) {
751 >                d = (log(l)*(double)ambvwt + avsum) /
752 >                                (double)(ambvwt + navsum);
753 >                d = exp(d) / l;
754 >                scalecolor(aval, d);
755 >                multcolor(aval, ambval);        /* apply color of ambval */
756 >        } else {
757 >                d = exp( avsum / (double)navsum );
758 >                scalecolor(aval, d);            /* neutral color */
759 >        }
760   }
761  
762  
763 < double
764 < sumambient(acol, r, rn, al, at, c0, s)  /* get interpolated ambient value */
765 < COLOR  acol;
766 < register RAY  *r;
767 < FVECT  rn;
768 < int  al;
769 < AMBTREE  *at;
770 < FVECT  c0;
771 < double  s;
763 > static double
764 > sumambient(     /* get interpolated ambient value */
765 >        COLOR  acol,
766 >        RAY  *r,
767 >        FVECT  rn,
768 >        int  al,
769 >        AMBTREE  *at,
770 >        FVECT  c0,
771 >        double  s
772 > )
773   {
774          double  d, e1, e2, wt, wsum;
775          COLOR  ct;
776          FVECT  ck0;
777          int  i;
778 <        register int  j;
779 <        register AMBVAL  *av;
778 >        int  j;
779 >        AMBVAL   *av;
780  
781          wsum = 0.0;
782                                          /* do this node */
783          for (av = at->alist; av != NULL; av = av->next) {
784 +                double  rn_dot = -2.0;
785                  if (tracktime)
786 <                        av->latick = ambclock++;
786 >                        av->latick = ambclock;
787                  /*
788                   *  Ambient level test.
789                   */
790 <                if (av->lvl > al)       /* list sorted, so this works */
790 >                if (av->lvl > al ||     /* list sorted, so this works */
791 >                                (av->lvl == al) & (av->weight < 0.9*r->rweight))
792                          break;
260                if (av->weight < r->rweight-FTINY)
261                        continue;
793                  /*
794                   *  Ambient radius test.
795                   */
796 <                d = av->pos[0] - r->rop[0];
797 <                e1 = d * d;
267 <                d = av->pos[1] - r->rop[1];
268 <                e1 += d * d;
269 <                d = av->pos[2] - r->rop[2];
270 <                e1 += d * d;
271 <                e1 /= av->rad * av->rad;
796 >                VSUB(ck0, av->pos, r->rop);
797 >                e1 = DOT(ck0, ck0) / (av->rad * av->rad);
798                  if (e1 > ambacc*ambacc*1.21)
799                          continue;
800                  /*
801 <                 *  Normal direction test.
801 >                 *  Direction test using closest normal.
802                   */
803 <                e2 = (1.0 - DOT(av->dir, r->ron)) * r->rweight;
804 <                if (e2 < 0.0) e2 = 0.0;
805 <                if (e1 + e2 > ambacc*ambacc*1.21)
803 >                d = DOT(av->dir, r->ron);
804 >                if (rn != r->ron) {
805 >                        rn_dot = DOT(av->dir, rn);
806 >                        if (rn_dot > 1.0-FTINY)
807 >                                rn_dot = 1.0-FTINY;
808 >                        if (rn_dot >= d-FTINY) {
809 >                                d = rn_dot;
810 >                                rn_dot = -2.0;
811 >                        }
812 >                }
813 >                e2 = (1.0 - d) * r->rweight;
814 >                if (e2 < 0.0)
815 >                        e2 = 0.0;
816 >                else if (e1 + e2 > ambacc*ambacc*1.21)
817                          continue;
818                  /*
819                   *  Ray behind test.
# Line 290 | Line 827 | double s;
827                  /*
828                   *  Jittering final test reduces image artifacts.
829                   */
830 <                wt = sqrt(e1) + sqrt(e2);
830 >                e1 = sqrt(e1);
831 >                e2 = sqrt(e2);
832 >                wt = e1 + e2;
833                  if (wt > ambacc*(.9+.2*urand(9015+samplendx)))
834                          continue;
835 +                /*
836 +                 *  Recompute directional error using perturbed normal
837 +                 */
838 +                if (rn_dot > 0.0) {
839 +                        e2 = sqrt((1.0 - rn_dot)*r->rweight);
840 +                        wt = e1 + e2;
841 +                }
842                  if (wt <= 1e-3)
843                          wt = 1e3;
844                  else
# Line 317 | Line 863 | double s;
863                                  break;
864                  }
865                  if (j == 3)
866 <                        wsum += sumambient(acol, r, rn, al, at->kid+i, ck0, s);
866 >                        wsum += sumambient(acol, r, rn, al,
867 >                                                at->kid+i, ck0, s);
868          }
869          return(wsum);
870   }
871  
872  
873 < double
874 < makeambient(acol, r, rn, al)    /* make a new ambient value */
875 < COLOR  acol;
876 < register RAY  *r;
877 < FVECT  rn;
878 < int  al;
873 > static double
874 > makeambient(            /* make a new ambient value for storage */
875 >        COLOR  acol,
876 >        RAY  *r,
877 >        FVECT  rn,
878 >        int  al
879 > )
880   {
881          AMBVAL  amb;
882          FVECT   gp, gd;
883 <                                                /* compute weight */
884 <        amb.weight = pow(AVGREFL, (double)al);
885 <        if (r->rweight < 0.1*amb.weight)        /* heuristic */
886 <                amb.weight = r->rweight;
883 >        int     i;
884 >
885 >        amb.weight = 1.0;                       /* compute weight */
886 >        for (i = al; i-- > 0; )
887 >                amb.weight *= AVGREFL;
888 >        if (r->rweight < 0.1*amb.weight)        /* heuristic override */
889 >                amb.weight = 1.25*r->rweight;
890 >        setcolor(acol, AVGREFL, AVGREFL, AVGREFL);
891                                                  /* compute ambient */
892          amb.rad = doambient(acol, r, amb.weight, gp, gd);
893 <        if (amb.rad <= FTINY)
893 >        if (amb.rad <= FTINY) {
894 >                setcolor(acol, 0.0, 0.0, 0.0);
895                  return(0.0);
896 <                                                /* store it */
896 >        }
897 >        scalecolor(acol, 1./AVGREFL);           /* undo assumed reflectance */
898 >                                                /* store value */
899          VCOPY(amb.pos, r->rop);
900          VCOPY(amb.dir, r->ron);
901          amb.lvl = al;
# Line 355 | Line 910 | int  al;
910   }
911  
912  
913 < extambient(cr, ap, pv, nv)              /* extrapolate value at pv, nv */
914 < COLOR  cr;
915 < register AMBVAL  *ap;
916 < FVECT  pv, nv;
913 > static void
914 > extambient(             /* extrapolate value at pv, nv */
915 >        COLOR  cr,
916 >        AMBVAL   *ap,
917 >        FVECT  pv,
918 >        FVECT  nv
919 > )
920   {
921 <        FVECT  v1, v2;
922 <        register int  i;
921 >        FVECT  v1;
922 >        int  i;
923          double  d;
924  
925          d = 1.0;                        /* zeroeth order */
# Line 369 | Line 927 | FVECT  pv, nv;
927          for (i = 0; i < 3; i++)
928                  d += ap->gpos[i]*(pv[i]-ap->pos[i]);
929                                          /* gradient due to rotation */
930 <        VCOPY(v1, ap->dir);
931 <        fcross(v2, v1, nv);
374 <        d += DOT(ap->gdir, v2);
930 >        VCROSS(v1, ap->dir, nv);
931 >        d += DOT(ap->gdir, v1);
932          if (d <= 0.0) {
933                  setcolor(cr, 0.0, 0.0, 0.0);
934                  return;
# Line 381 | Line 938 | FVECT  pv, nv;
938   }
939  
940  
941 < static
942 < initambfile(creat)              /* initialize ambient file */
943 < int  creat;
941 > static void
942 > avinsert(                               /* insert ambient value in our tree */
943 >        AMBVAL *av
944 > )
945   {
946 <        extern char  *progname, *octname, VersionID[];
946 >        AMBTREE  *at;
947 >        AMBVAL  *ap;
948 >        AMBVAL  avh;
949 >        FVECT  ck0;
950 >        double  s;
951 >        int  branch;
952 >        int  i;
953  
954 +        if (av->rad <= FTINY)
955 +                error(CONSISTENCY, "zero ambient radius in avinsert");
956 +        at = &atrunk;
957 +        VCOPY(ck0, thescene.cuorg);
958 +        s = thescene.cusize;
959 +        while (s*(OCTSCALE/2) > av->rad*ambacc) {
960 +                if (at->kid == NULL)
961 +                        if ((at->kid = newambtree()) == NULL)
962 +                                error(SYSTEM, "out of memory in avinsert");
963 +                s *= 0.5;
964 +                branch = 0;
965 +                for (i = 0; i < 3; i++)
966 +                        if (av->pos[i] > ck0[i] + s) {
967 +                                ck0[i] += s;
968 +                                branch |= 1 << i;
969 +                        }
970 +                at = at->kid + branch;
971 +        }
972 +        avh.next = at->alist;           /* order by increasing level */
973 +        for (ap = &avh; ap->next != NULL; ap = ap->next)
974 +                if ( ap->next->lvl > av->lvl ||
975 +                                (ap->next->lvl == av->lvl) &
976 +                                (ap->next->weight <= av->weight) )
977 +                        break;
978 +        av->next = ap->next;
979 +        ap->next = (AMBVAL*)av;
980 +        at->alist = avh.next;
981 + }
982 +
983 + #endif  /* ! NEWAMB */
984 +
985 + /************* FOLLOWING ROUTINES SAME FOR NEW & OLD METHODS ***************/
986 +
987 + static void
988 + initambfile(            /* initialize ambient file */
989 +        int  cre8
990 + )
991 + {
992 +        extern char  *progname, *octname;
993 +        static char  *mybuf = NULL;
994 +
995   #ifdef  F_SETLKW
996 <        aflock(creat ? F_WRLCK : F_RDLCK);
996 >        aflock(cre8 ? F_WRLCK : F_RDLCK);
997   #endif
998 < #ifdef MSDOS
999 <        setmode(fileno(ambfp), O_BINARY);
1000 < #endif
1001 <        setbuf(ambfp, bmalloc(BUFSIZ+8));
1002 <        if (creat) {                    /* new file */
998 >        SET_FILE_BINARY(ambfp);
999 >        if (mybuf == NULL)
1000 >                mybuf = (char *)bmalloc(BUFSIZ+8);
1001 >        setbuf(ambfp, mybuf);
1002 >        if (cre8) {                     /* new file */
1003                  newheader("RADIANCE", ambfp);
1004 <                fprintf(ambfp, "%s -av %g %g %g -ab %d -aa %g ",
1004 >                fprintf(ambfp, "%s -av %g %g %g -aw %d -ab %d -aa %g ",
1005                                  progname, colval(ambval,RED),
1006                                  colval(ambval,GRN), colval(ambval,BLU),
1007 <                                ambounce, ambacc);
1008 <                fprintf(ambfp, "-ad %d -as %d -ar %d %s\n",
1009 <                                ambdiv, ambssamp, ambres,
1010 <                                octname==NULL ? "" : octname);
1007 >                                ambvwt, ambounce, ambacc);
1008 >                fprintf(ambfp, "-ad %d -as %d -ar %d ",
1009 >                                ambdiv, ambssamp, ambres);
1010 >                if (octname != NULL)
1011 >                        fputs(octname, ambfp);
1012 >                fputc('\n', ambfp);
1013                  fprintf(ambfp, "SOFTWARE= %s\n", VersionID);
1014 +                fputnow(ambfp);
1015                  fputformat(AMBFMT, ambfp);
1016 <                putc('\n', ambfp);
1016 >                fputc('\n', ambfp);
1017                  putambmagic(ambfp);
1018          } else if (checkheader(ambfp, AMBFMT, NULL) < 0 || !hasambmagic(ambfp))
1019                  error(USER, "bad ambient file");
1020   }
1021  
1022  
1023 < static
1024 < avsave(av)                              /* insert and save an ambient value */
1025 < AMBVAL  *av;
1023 > static void
1024 > avsave(                         /* insert and save an ambient value */
1025 >        AMBVAL  *av
1026 > )
1027   {
1028 <        avinsert(avstore(av));
1028 >        avstore(av);
1029          if (ambfp == NULL)
1030                  return;
1031          if (writambval(av, ambfp) < 0)
# Line 426 | Line 1035 | AMBVAL *av;
1035                          goto writerr;
1036          return;
1037   writerr:
1038 <        error(SYSTEM, "error writing ambient file");
1038 >        error(SYSTEM, "error writing to ambient file");
1039   }
1040  
1041  
1042   static AMBVAL *
1043 < avstore(aval)                           /* allocate memory and store aval */
1044 < register AMBVAL  *aval;
1043 > avstore(                                /* allocate memory and save aval */
1044 >        AMBVAL  *aval
1045 > )
1046   {
1047 <        register AMBVAL  *av;
1047 >        AMBVAL  *av;
1048 >        double  d;
1049  
1050          if ((av = newambval()) == NULL)
1051                  error(SYSTEM, "out of memory in avstore");
1052 <        copystruct(av, aval);
1052 >        *av = *aval;
1053          av->latick = ambclock;
1054          av->next = NULL;
444        addcolor(avsum, av->val);       /* add to sum for averaging */
1055          nambvals++;
1056 +        d = bright(av->val);
1057 +        if (d > FTINY) {                /* add to log sum for averaging */
1058 +                avsum += log(d);
1059 +                navsum++;
1060 +        }
1061 +        avinsert(av);                   /* insert in our cache tree */
1062          return(av);
1063   }
1064  
# Line 452 | Line 1068 | register AMBVAL  *aval;
1068   static AMBTREE  *atfreelist = NULL;     /* free ambient tree structures */
1069  
1070  
1071 < static
1072 < AMBTREE *
457 < newambtree()                            /* allocate 8 ambient tree structs */
1071 > static AMBTREE *
1072 > newambtree(void)                                /* allocate 8 ambient tree structs */
1073   {
1074 <        register AMBTREE  *atp, *upperlim;
1074 >        AMBTREE  *atp, *upperlim;
1075  
1076          if (atfreelist == NULL) {       /* get more nodes */
1077 <                atfreelist = (AMBTREE *)bmalloc(ATALLOCSZ*8*sizeof(AMBTREE));
1077 >                atfreelist = (AMBTREE *)malloc(ATALLOCSZ*8*sizeof(AMBTREE));
1078                  if (atfreelist == NULL)
1079                          return(NULL);
1080                                          /* link new free list */
# Line 470 | Line 1085 | newambtree()                           /* allocate 8 ambient tree structs */
1085          }
1086          atp = atfreelist;
1087          atfreelist = atp->kid;
1088 <        bzero((char *)atp, 8*sizeof(AMBTREE));
1088 >        memset(atp, 0, 8*sizeof(AMBTREE));
1089          return(atp);
1090   }
1091  
1092  
1093 < static
1094 < freeambtree(atp)                        /* free 8 ambient tree structs */
1095 < AMBTREE  *atp;
1093 > static void
1094 > freeambtree(                    /* free 8 ambient tree structs */
1095 >        AMBTREE  *atp
1096 > )
1097   {
1098          atp->kid = atfreelist;
1099          atfreelist = atp;
1100   }
1101  
1102  
1103 < static
1104 < avinsert(av)                            /* insert ambient value in our tree */
1105 < register AMBVAL  *av;
1103 > static void
1104 > unloadatree(                    /* unload an ambient value tree */
1105 >        AMBTREE  *at,
1106 >        unloadtf_t *f
1107 > )
1108   {
1109 <        register AMBTREE  *at;
1110 <        register AMBVAL  *ap;
493 <        AMBVAL  avh;
494 <        FVECT  ck0;
495 <        double  s;
496 <        int  branch;
497 <        register int  i;
498 <
499 <        if (av->rad <= FTINY)
500 <                error(CONSISTENCY, "zero ambient radius in avinsert");
501 <        at = &atrunk;
502 <        VCOPY(ck0, thescene.cuorg);
503 <        s = thescene.cusize;
504 <        while (s*(OCTSCALE/2) > av->rad*ambacc) {
505 <                if (at->kid == NULL)
506 <                        if ((at->kid = newambtree()) == NULL)
507 <                                error(SYSTEM, "out of memory in avinsert");
508 <                s *= 0.5;
509 <                branch = 0;
510 <                for (i = 0; i < 3; i++)
511 <                        if (av->pos[i] > ck0[i] + s) {
512 <                                ck0[i] += s;
513 <                                branch |= 1 << i;
514 <                        }
515 <                at = at->kid + branch;
516 <        }
517 <        avh.next = at->alist;           /* order by increasing level */
518 <        for (ap = &avh; ap->next != NULL; ap = ap->next)
519 <                if (ap->next->lvl >= av->lvl)
520 <                        break;
521 <        av->next = ap->next;
522 <        ap->next = av;
523 <        at->alist = avh.next;
524 < }
525 <
526 <
527 < static
528 < unloadatree(at, f)                      /* unload an ambient value tree */
529 < register AMBTREE  *at;
530 < int     (*f)();
531 < {
532 <        register AMBVAL  *av;
533 <        register int  i;
1109 >        AMBVAL  *av;
1110 >        int  i;
1111                                          /* transfer values at this node */
1112          for (av = at->alist; av != NULL; av = at->alist) {
1113                  at->alist = av->next;
1114 +                av->next = NULL;
1115                  (*f)(av);
1116          }
1117          if (at->kid == NULL)
# Line 545 | Line 1123 | int    (*f)();
1123   }
1124  
1125  
1126 < static AMBVAL   **avlist1, **avlist2;   /* ambient value lists for sorting */
1126 > static struct avl {
1127 >        AMBVAL  *p;
1128 >        unsigned long   t;
1129 > }       *avlist1;                       /* ambient value list with ticks */
1130 > static AMBVAL   **avlist2;              /* memory positions for sorting */
1131   static int      i_avlist;               /* index for lists */
1132  
1133 + static int alatcmp(const void *av1, const void *av2);
1134  
1135 < static
1136 < av2list(av)
554 < AMBVAL  *av;
1135 > static void
1136 > avfree(AMBVAL *av)
1137   {
1138 +        free(av);
1139 + }
1140 +
1141 + static void
1142 + av2list(
1143 +        AMBVAL *av
1144 + )
1145 + {
1146   #ifdef DEBUG
1147          if (i_avlist >= nambvals)
1148                  error(CONSISTENCY, "too many ambient values in av2list1");
1149   #endif
1150 <        avlist1[i_avlist] = avlist2[i_avlist] = av;
1151 <        i_avlist++;
1150 >        avlist1[i_avlist].p = avlist2[i_avlist] = (AMBVAL*)av;
1151 >        avlist1[i_avlist++].t = av->latick;
1152   }
1153  
1154  
1155   static int
1156 < alatcmp(avp1, avp2)                     /* compare ambient values for MRA */
1157 < AMBVAL  **avp1, **avp2;
1156 > alatcmp(                        /* compare ambient values for MRA */
1157 >        const void *av1,
1158 >        const void *av2
1159 > )
1160   {
1161 <        return((**avp2).latick - (**avp1).latick);
1161 >        long  lc = ((struct avl *)av2)->t - ((struct avl *)av1)->t;
1162 >        return(lc<0 ? -1 : lc>0 ? 1 : 0);
1163   }
1164  
1165  
1166 + /* GW NOTE 2002/10/3:
1167 + * I used to compare AMBVAL pointers, but found that this was the
1168 + * cause of a serious consistency error with gcc, since the optimizer
1169 + * uses some dangerous trick in pointer subtraction that
1170 + * assumes pointers differ by exact struct size increments.
1171 + */
1172   static int
1173 < aposcmp(avp1, avp2)                     /* compare ambient value positions */
1174 < AMBVAL  **avp1, **avp2;
1173 > aposcmp(                        /* compare ambient value positions */
1174 >        const void      *avp1,
1175 >        const void      *avp2
1176 > )
1177   {
1178 <        return(*avp1 - *avp2);
1178 >        long    diff = *(char * const *)avp1 - *(char * const *)avp2;
1179 >        if (diff < 0)
1180 >                return(-1);
1181 >        return(diff > 0);
1182   }
1183  
1184  
581 #ifdef DEBUG
1185   static int
1186 < avlmemi(avaddr)                         /* find list position from address */
1187 < AMBVAL  *avaddr;
1186 > avlmemi(                                /* find list position from address */
1187 >        AMBVAL  *avaddr
1188 > )
1189   {
1190 <        register AMBVAL  **avlpp;
1190 >        AMBVAL  **avlpp;
1191  
1192 <        avlpp = (AMBVAL **)bsearch((char *)&avaddr, (char *)avlist2,
1192 >        avlpp = (AMBVAL **)bsearch(&avaddr, avlist2,
1193                          nambvals, sizeof(AMBVAL *), aposcmp);
1194          if (avlpp == NULL)
1195                  error(CONSISTENCY, "address not found in avlmemi");
1196          return(avlpp - avlist2);
1197   }
594 #else
595 #define avlmemi(avaddr) ((AMBVAL **)bsearch((char *)&avaddr,(char *)avlist2, \
596                                nambvals,sizeof(AMBVAL *),aposcmp) - avlist2)
597 #endif
1198  
1199  
1200 < static
1201 < sortambvals(always)                     /* resort ambient values */
1202 < int     always;
1200 > static void
1201 > sortambvals(                    /* resort ambient values */
1202 >        int     always
1203 > )
1204   {
1205          AMBTREE  oldatrunk;
1206          AMBVAL  tav, *tap, *pnext;
1207 <        register int    i, j;
1207 >        int     i, j;
1208                                          /* see if it's time yet */
1209 <        if (!always && (ambclock < lastsort+sortintvl ||
1209 >        if (!always && (ambclock++ < lastsort+sortintvl ||
1210                          nambvals < SORT_THRESH))
1211                  return;
1212          /*
# Line 624 | Line 1225 | int    always;
1225           * tree will be rebuilt with the new accuracy parameter.
1226           */
1227          if (tracktime) {                /* allocate pointer arrays to sort */
627                avlist1 = (AMBVAL **)malloc(nambvals*sizeof(AMBVAL *));
1228                  avlist2 = (AMBVAL **)malloc(nambvals*sizeof(AMBVAL *));
1229 <        } else
1230 <                avlist1 = avlist2 = NULL;
1231 <        if (avlist2 == NULL) {          /* no time tracking -- rebuild tree? */
1232 <                if (avlist1 != NULL)
1233 <                        free((char *)avlist1);
1229 >                avlist1 = (struct avl *)malloc(nambvals*sizeof(struct avl));
1230 >        } else {
1231 >                avlist2 = NULL;
1232 >                avlist1 = NULL;
1233 >        }
1234 >        if (avlist1 == NULL) {          /* no time tracking -- rebuild tree? */
1235 >                if (avlist2 != NULL)
1236 >                        free(avlist2);
1237                  if (always) {           /* rebuild without sorting */
1238 <                        copystruct(&oldatrunk, &atrunk);
1238 >                        oldatrunk = atrunk;
1239                          atrunk.alist = NULL;
1240                          atrunk.kid = NULL;
1241                          unloadatree(&oldatrunk, avinsert);
# Line 657 | Line 1260 | int    always;
1260                  if (i_avlist < nambvals)
1261                          error(CONSISTENCY, "missing ambient values in sortambvals");
1262   #endif
1263 <                qsort((char *)avlist1, nambvals, sizeof(AMBVAL *), alatcmp);
1264 <                qsort((char *)avlist2, nambvals, sizeof(AMBVAL *), aposcmp);
1263 >                qsort(avlist1, nambvals, sizeof(struct avl), alatcmp);
1264 >                qsort(avlist2, nambvals, sizeof(AMBVAL *), aposcmp);
1265                  for (i = 0; i < nambvals; i++) {
1266 <                        if (avlist1[i] == NULL)
1266 >                        if (avlist1[i].p == NULL)
1267                                  continue;
1268                          tap = avlist2[i];
1269 <                        copystruct(&tav, tap);
1270 <                        for (j = i; (pnext = avlist1[j]) != tap;
1269 >                        tav = *tap;
1270 >                        for (j = i; (pnext = avlist1[j].p) != tap;
1271                                          j = avlmemi(pnext)) {
1272 <                                copystruct(avlist2[j], pnext);
1272 >                                *(avlist2[j]) = *pnext;
1273                                  avinsert(avlist2[j]);
1274 <                                avlist1[j] = NULL;
1274 >                                avlist1[j].p = NULL;
1275                          }
1276 <                        copystruct(avlist2[j], &tav);
1276 >                        *(avlist2[j]) = tav;
1277                          avinsert(avlist2[j]);
1278 <                        avlist1[j] = NULL;
1278 >                        avlist1[j].p = NULL;
1279                  }
1280 <                free((char *)avlist1);
1281 <                free((char *)avlist2);
1280 >                free(avlist1);
1281 >                free(avlist2);
1282                                                  /* compute new sort interval */
1283                  sortintvl = ambclock - lastsort;
1284                  if (sortintvl >= MAX_SORT_INTVL/2)
# Line 694 | Line 1297 | int    always;
1297  
1298   #ifdef  F_SETLKW
1299  
1300 < static
1301 < aflock(typ)                     /* lock/unlock ambient file */
1302 < int  typ;
1300 > static void
1301 > aflock(                 /* lock/unlock ambient file */
1302 >        int  typ
1303 > )
1304   {
1305          static struct flock  fls;       /* static so initialized to zeroes */
1306  
1307 +        if (typ == fls.l_type)          /* already called? */
1308 +                return;
1309          fls.l_type = typ;
1310          if (fcntl(fileno(ambfp), F_SETLKW, &fls) < 0)
1311                  error(SYSTEM, "cannot (un)lock ambient file");
# Line 707 | Line 1313 | int  typ;
1313  
1314  
1315   int
1316 < ambsync()                       /* synchronize ambient file */
1316 > ambsync(void)                   /* synchronize ambient file */
1317   {
712        static FILE  *ambinp = NULL;
713        static long  lastpos = -1;
1318          long  flen;
1319          AMBVAL  avs;
1320 <        register int  n;
1320 >        int  n;
1321  
1322 <        if (nunflshed == 0)
1322 >        if (ambfp == NULL)      /* no ambient file? */
1323                  return(0);
1324 <        if (lastpos < 0)        /* initializing (locked in initambfile) */
1325 <                goto syncend;
722 <                                /* gain exclusive access */
723 <        aflock(F_WRLCK);
1324 >                                /* gain appropriate access */
1325 >        aflock(nunflshed ? F_WRLCK : F_RDLCK);
1326                                  /* see if file has grown */
1327 <        if ((flen = lseek(fileno(ambfp), 0L, 2)) < 0)
1327 >        if ((flen = lseek(fileno(ambfp), (off_t)0, SEEK_END)) < 0)
1328                  goto seekerr;
1329 <        if (n = flen - lastpos) {               /* file has grown */
1329 >        if ((n = flen - lastpos) > 0) {         /* file has grown */
1330                  if (ambinp == NULL) {           /* use duplicate filedes */
1331                          ambinp = fdopen(dup(fileno(ambfp)), "r");
1332                          if (ambinp == NULL)
1333                                  error(SYSTEM, "fdopen failed in ambsync");
1334                  }
1335 <                if (fseek(ambinp, lastpos, 0) < 0)
1335 >                if (fseek(ambinp, lastpos, SEEK_SET) < 0)
1336                          goto seekerr;
1337                  while (n >= AMBVALSIZ) {        /* load contributed values */
1338 <                        readambval(&avs, ambinp);
1339 <                        avinsert(avstore(&avs));
1338 >                        if (!readambval(&avs, ambinp)) {
1339 >                                sprintf(errmsg,
1340 >                        "ambient file \"%s\" corrupted near character %ld",
1341 >                                                ambfile, flen - n);
1342 >                                error(WARNING, errmsg);
1343 >                                break;
1344 >                        }
1345 >                        avstore(&avs);
1346                          n -= AMBVALSIZ;
1347                  }
1348 +                lastpos = flen - n;
1349                  /*** seek always as safety measure
1350                  if (n) ***/                     /* alignment */
1351 <                        if (lseek(fileno(ambfp), flen-n, 0) < 0)
1351 >                        if (lseek(fileno(ambfp), (off_t)lastpos, SEEK_SET) < 0)
1352                                  goto seekerr;
1353          }
745 #ifdef  DEBUG
746        if (ambfp->_ptr - ambfp->_base != nunflshed*AMBVALSIZ) {
747                sprintf(errmsg, "ambient file buffer at %d rather than %d",
748                                ambfp->_ptr - ambfp->_base,
749                                nunflshed*AMBVALSIZ);
750                error(CONSISTENCY, errmsg);
751        }
752 #endif
753 syncend:
1354          n = fflush(ambfp);                      /* calls write() at last */
1355 <        if ((lastpos = lseek(fileno(ambfp), 0L, 1)) < 0)
1355 >        if (n != EOF)
1356 >                lastpos += (long)nunflshed*AMBVALSIZ;
1357 >        else if ((lastpos = lseek(fileno(ambfp), (off_t)0, SEEK_CUR)) < 0)
1358                  goto seekerr;
1359 +                
1360          aflock(F_UNLCK);                        /* release file */
1361          nunflshed = 0;
1362          return(n);
1363   seekerr:
1364          error(SYSTEM, "seek failed in ambsync");
1365 +        return -1; /* pro forma return */
1366   }
1367  
1368 < #else
1368 > #else   /* ! F_SETLKW */
1369  
1370   int
1371 < ambsync()                       /* flush ambient file */
1371 > ambsync(void)                   /* flush ambient file */
1372   {
1373 <        if (nunflshed == 0)
1373 >        if (ambfp == NULL)
1374                  return(0);
1375          nunflshed = 0;
1376          return(fflush(ambfp));
1377   }
1378  
1379 < #endif
1379 > #endif  /* ! F_SETLKW */

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines