--- ray/src/rt/ambient.c	1989/05/26 17:58:01	1.3
+++ ray/src/rt/ambient.c	2003/06/26 00:58:10	2.51
@@ -1,127 +1,251 @@
-/* Copyright (c) 1986 Regents of the University of California */
-
 #ifndef lint
-static char SCCSid[] = "$SunId$ LBL";
+static const char	RCSid[] = "$Id: ambient.c,v 2.51 2003/06/26 00:58:10 schorsch Exp $";
 #endif
-
 /*
  *  ambient.c - routines dealing with ambient (inter-reflected) component.
  *
- *  The macro AMBFLUSH (if defined) is the number of ambient values
- *	to wait before flushing to the ambient file.
- *
- *     5/9/86
+ *  Declarations of external symbols in ambient.h
  */
 
+#include "copyright.h"
+
+#include  "platform.h"
 #include  "ray.h"
+#include  "otypes.h"
+#include  "ambient.h"
+#include  "random.h"
 
-#include  "octree.h"
+#ifndef  OCTSCALE
+#define	 OCTSCALE	1.0	/* ceil((valid rad.)/(cube size)) */
+#endif
 
-#include  "random.h"
+extern char  *shm_boundary;	/* memory sharing boundary */
 
-#define  OCTSCALE	0.5	/* ceil((valid rad.)/(cube size)) */
+#define	 MAXASET	511	/* maximum number of elements in ambient set */
+OBJECT	ambset[MAXASET+1]={0};	/* ambient include/exclude set */
 
-extern CUBE  thescene;		/* contains space boundaries */
+double	maxarad;		/* maximum ambient radius */
+double	minarad;		/* minimum ambient radius */
 
-extern COLOR  ambval;		/* global ambient component */
-extern double  ambacc;		/* ambient accuracy */
-extern int  ambres;		/* ambient resolution */
-extern int  ambdiv;		/* number of divisions for calculation */
-extern int  ambssamp;		/* number of super-samples */
-extern int  ambounce;		/* number of ambient bounces */
-extern char  *amblist[];	/* ambient include/exclude list */
-extern int  ambincl;		/* include == 1, exclude == 0 */
+static AMBTREE	atrunk;		/* our ambient trunk node */
 
-OBJECT  ambset[128];		/* ambient include/exclude set */
+static FILE  *ambfp = NULL;	/* ambient file pointer */
+static int  nunflshed = 0;	/* number of unflushed ambient values */
 
-double  maxarad;		/* maximum ambient radius */
-double  minarad;		/* minimum ambient radius */
+#ifndef SORT_THRESH
+#ifdef SMLMEM
+#define SORT_THRESH	((3L<<20)/sizeof(AMBVAL))
+#else
+#define SORT_THRESH	((9L<<20)/sizeof(AMBVAL))
+#endif
+#endif
+#ifndef SORT_INTVL
+#define SORT_INTVL	(SORT_THRESH<<1)
+#endif
+#ifndef MAX_SORT_INTVL
+#define MAX_SORT_INTVL	(SORT_INTVL<<6)
+#endif
 
-typedef struct ambval {
-	FVECT  pos;		/* position in space */
-	FVECT  dir;		/* normal direction */
-	int  lvl;		/* recursion level of parent ray */
-	float  weight;		/* weight of parent ray */
-	COLOR  val;		/* computed ambient value */
-	float  rad;		/* validity radius */
-	struct ambval  *next;	/* next in list */
-}  AMBVAL;			/* ambient value */
+static double  avsum = 0.;		/* computed ambient value sum (log) */
+static unsigned int  navsum = 0;	/* number of values in avsum */
+static unsigned int  nambvals = 0;	/* total number of indirect values */
+static unsigned int  nambshare = 0;	/* number of values from file */
+static unsigned long  ambclock = 0;	/* ambient access clock */
+static unsigned long  lastsort = 0;	/* time of last value sort */
+static long  sortintvl = SORT_INTVL;	/* time until next sort */
+static FILE  *ambinp = NULL;		/* auxiliary file for input */
+static long  lastpos = -1;		/* last flush position */
 
-typedef struct ambtree {
-	AMBVAL  *alist;		/* ambient value list */
-	struct ambtree  *kid;	/* 8 child nodes */
-}  AMBTREE;			/* ambient octree */
+#define MAXACLOCK	(1L<<30)	/* clock turnover value */
+	/*
+	 * Track access times unless we are sharing ambient values
+	 * through memory on a multiprocessor, when we want to avoid
+	 * claiming our own memory (copy on write).  Go ahead anyway
+	 * if more than two thirds of our values are unshared.
+	 * Compile with -Dtracktime=0 to turn this code off.
+	 */
+#ifndef tracktime
+#define tracktime	(shm_boundary == NULL || nambvals > 3*nambshare)
+#endif
 
-typedef struct {
-	float  k;		/* error contribution per sample */
-	COLOR  v;		/* ray sum */
-	int  n;			/* number of samples */
-	short  t, p;		/* theta, phi indices */
-}  AMBSAMP;			/* ambient sample */
+#define	 AMBFLUSH	(BUFSIZ/AMBVALSIZ)
 
-static AMBTREE  atrunk;		/* our ambient trunk node */
+#define	 newambval()	(AMBVAL *)malloc(sizeof(AMBVAL))
+#define  freeav(av)	free((void *)av);
 
-static FILE  *ambfp = NULL;	/* ambient file pointer */
+static void  initambfile(), avsave(), avinsert(), sortambvals(), unloadatree();
+static int  avlmemi();
+static AMBVAL  *avstore();
+#ifdef  F_SETLKW
+static void  aflock();
+#endif
 
-#define  newambval()	(AMBVAL *)bmalloc(sizeof(AMBVAL))
 
-#define  newambtree()	(AMBTREE *)calloc(8, sizeof(AMBTREE))
+void
+setambres(ar)				/* set ambient resolution */
+int  ar;
+{
+	ambres = ar < 0 ? 0 : ar;		/* may be done already */
+						/* set min & max radii */
+	if (ar <= 0) {
+		minarad = 0;
+		maxarad = thescene.cusize / 2.0;
+	} else {
+		minarad = thescene.cusize / ar;
+		maxarad = 64 * minarad;			/* heuristic */
+		if (maxarad > thescene.cusize / 2.0)
+			maxarad = thescene.cusize / 2.0;
+	}
+	if (minarad <= FTINY)
+		minarad = 10*FTINY;
+	if (maxarad <= minarad)
+		maxarad = 64 * minarad;
+}
 
-double  sumambient(), doambient(), makeambient();
 
+void
+setambacc(newa)				/* set ambient accuracy */
+double  newa;
+{
+	double  ambdiff;
 
-setambient(afile)			/* initialize calculation */
-char  *afile;
+	if (newa < 0.0)
+		newa = 0.0;
+	ambdiff = fabs(newa - ambacc);
+	if (ambdiff >= .01 && (ambacc = newa) > FTINY && nambvals > 0)
+		sortambvals(1);			/* rebuild tree */
+}
+
+
+void
+setambient()				/* initialize calculation */
 {
-	long  ftell();
-	char  **amblp;
-	OBJECT  obj;
-	AMBVAL  amb;
-					/* set up ambient set */
-	ambset[0] = 0;
-	for (amblp = amblist; *amblp != NULL; amblp++) {
-		if ((obj = modifier(*amblp)) == OVOID) {
-			sprintf(errmsg, "unknown %s modifier \"%s\"", 
-				ambincl ? "include" : "exclude", *amblp);
+	int	readonly = 0;
+	long  pos, flen;
+	AMBVAL	amb;
+						/* make sure we're fresh */
+	ambdone();
+						/* init ambient limits */
+	setambres(ambres);
+	setambacc(ambacc);
+	if (ambfile == NULL || !ambfile[0])
+		return;
+	if (ambacc <= FTINY) {
+		sprintf(errmsg, "zero ambient accuracy so \"%s\" not opened",
+				ambfile);
+		error(WARNING, errmsg);
+		return;
+	}
+						/* open ambient file */
+	if ((ambfp = fopen(ambfile, "r+")) == NULL)
+		readonly = (ambfp = fopen(ambfile, "r")) != NULL;
+	if (ambfp != NULL) {
+		initambfile(0);			/* file exists */
+		pos = ftell(ambfp);
+		while (readambval(&amb, ambfp))
+			avinsert(avstore(&amb));
+		nambshare = nambvals;		/* share loaded values */
+		if (readonly) {
+			sprintf(errmsg,
+				"loaded %u values from read-only ambient file",
+					nambvals);
 			error(WARNING, errmsg);
-			continue;
+			fclose(ambfp);		/* close file so no writes */
+			ambfp = NULL;
+			return;			/* avoid ambsync() */
 		}
-		if (!inset(ambset, obj))
-			insertelem(ambset, obj);
+						/* align file pointer */
+		pos += (long)nambvals*AMBVALSIZ;
+		flen = lseek(fileno(ambfp), (off_t)0L, 2);
+		if (flen != pos) {
+			sprintf(errmsg,
+			"ignoring last %ld values in ambient file (corrupted)",
+					(flen - pos)/AMBVALSIZ);
+			error(WARNING, errmsg);
+			fseek(ambfp, pos, 0);
+#ifndef _WIN32 /* XXX we need a replacement for that one */
+			ftruncate(fileno(ambfp), (off_t)pos);
+#endif
+		}
+	} else if ((ambfp = fopen(ambfile, "w+")) != NULL) {
+		initambfile(1);			/* else create new file */
+	} else {
+		sprintf(errmsg, "cannot open ambient file \"%s\"", ambfile);
+		error(SYSTEM, errmsg);
 	}
-	maxarad = thescene.cusize / 2.0;		/* maximum radius */
-							/* minimum radius */
-	minarad = ambres > 0 ? thescene.cusize/ambres : 0.0;
+	nunflshed++;	/* lie */
+	ambsync();
+}
 
-					/* open ambient file */
-	if (afile != NULL)
-		if ((ambfp = fopen(afile, "r+")) != NULL) {
-			while (fread(&amb, sizeof(AMBVAL), 1, ambfp) == 1)
-				avinsert(&amb, &atrunk, thescene.cuorg,
-						thescene.cusize);
-							/* align */
-			fseek(ambfp, -(ftell(ambfp)%sizeof(AMBVAL)), 1);
-		} else if ((ambfp = fopen(afile, "w")) == NULL) {
-			sprintf(errmsg, "cannot open ambient file \"%s\"",
-					afile);
-			error(SYSTEM, errmsg);
+
+void
+ambdone()			/* close ambient file and free memory */
+{
+	if (ambfp != NULL) {		/* close ambient file */
+		ambsync();
+		fclose(ambfp);
+		ambfp = NULL;
+		if (ambinp != NULL) {	
+			fclose(ambinp);
+			ambinp = NULL;
 		}
+		lastpos = -1;
+	}
+					/* free ambient tree */
+	unloadatree(&atrunk, free);
+					/* reset state variables */
+	avsum = 0.;
+	navsum = 0;
+	nambvals = 0;
+	nambshare = 0;
+	ambclock = 0;
+	lastsort = 0;
+	sortintvl = SORT_INTVL;
 }
 
 
-ambient(acol, r)		/* compute ambient component for ray */
+void
+ambnotify(obj)			/* record new modifier */
+OBJECT	obj;
+{
+	static int  hitlimit = 0;
+	register OBJREC	 *o;
+	register char  **amblp;
+
+	if (obj == OVOID) {		/* starting over */
+		ambset[0] = 0;
+		hitlimit = 0;
+		return;
+	}
+	o = objptr(obj);
+	if (hitlimit || !ismodifier(o->otype))
+		return;
+	for (amblp = amblist; *amblp != NULL; amblp++)
+		if (!strcmp(o->oname, *amblp)) {
+			if (ambset[0] >= MAXASET) {
+				error(WARNING, "too many modifiers in ambient list");
+				hitlimit++;
+				return;		/* should this be fatal? */
+			}
+			insertelem(ambset, obj);
+			return;
+		}
+}
+
+
+void
+ambient(acol, r, nrm)		/* compute ambient component for ray */
 COLOR  acol;
 register RAY  *r;
+FVECT  nrm;
 {
 	static int  rdepth = 0;			/* ambient recursion */
-	double  wsum;
+	double	d, l;
 
-	rdepth++;				/* increment level */
-
 	if (ambdiv <= 0)			/* no ambient calculation */
 		goto dumbamb;
 						/* check number of bounces */
-	if (rdepth > ambounce)
+	if (rdepth >= ambounce)
 		goto dumbamb;
 						/* check ambient list */
 	if (ambincl != -1 && r->ro != NULL &&
@@ -129,64 +253,102 @@ register RAY  *r;
 		goto dumbamb;
 
 	if (ambacc <= FTINY) {			/* no ambient storage */
-		if (doambient(acol, r) == 0.0)
+		rdepth++;
+		d = doambient(acol, r, r->rweight, NULL, NULL);
+		rdepth--;
+		if (d <= FTINY)
 			goto dumbamb;
-		goto done;
+		return;
 	}
+
+	if (tracktime)				/* sort to minimize thrashing */
+		sortambvals(0);
 						/* get ambient value */
 	setcolor(acol, 0.0, 0.0, 0.0);
-	wsum = sumambient(acol, r, &atrunk, thescene.cuorg, thescene.cusize);
-	if (wsum > FTINY)
-		scalecolor(acol, 1.0/wsum);
-	else if (makeambient(acol, r) == 0.0)
-		goto dumbamb;
-	goto done;
-
+	d = sumambient(acol, r, nrm, rdepth,
+			&atrunk, thescene.cuorg, thescene.cusize);
+	if (d > FTINY) {
+		scalecolor(acol, 1.0/d);
+		return;
+	}
+	rdepth++;				/* need to cache new value */
+	d = makeambient(acol, r, nrm, rdepth-1);
+	rdepth--;
+	if (d > FTINY)
+		return;
 dumbamb:					/* return global value */
 	copycolor(acol, ambval);
-done:						/* must finish here! */
-	rdepth--;
+	if (ambvwt <= 0 | navsum == 0)
+		return;
+	l = bright(ambval);			/* average in computations */
+	if (l > FTINY) {
+		d = (log(l)*(double)ambvwt + avsum) /
+				(double)(ambvwt + navsum);
+		d = exp(d) / l;
+		scalecolor(acol, d);		/* apply color of ambval */
+	} else {
+		d = exp( avsum / (double)navsum );
+		setcolor(acol, d, d, d);	/* neutral color */
+	}
 }
 
 
 double
-sumambient(acol, r, at, c0, s)		/* get interpolated ambient value */
+sumambient(acol, r, rn, al, at, c0, s)	/* get interpolated ambient value */
 COLOR  acol;
 register RAY  *r;
-AMBTREE  *at;
+FVECT  rn;
+int  al;
+AMBTREE	 *at;
 FVECT  c0;
-double  s;
+double	s;
 {
-	extern double  sqrt();
-	double  d, e1, e2, wt, wsum;
+	double	d, e1, e2, wt, wsum;
 	COLOR  ct;
 	FVECT  ck0;
 	int  i;
 	register int  j;
-	register AMBVAL  *av;
-					/* do this node */
+	register AMBVAL	 *av;
+
 	wsum = 0.0;
+					/* do this node */
 	for (av = at->alist; av != NULL; av = av->next) {
+		double	rn_dot = -2.0;
+		if (tracktime)
+			av->latick = ambclock;
 		/*
-		 *  Ray strength test.
+		 *  Ambient level test.
 		 */
-		if (av->lvl > r->rlvl || av->weight < r->rweight-FTINY)
+		if (av->lvl > al)	/* list sorted, so this works */
+			break;
+		if (av->weight < r->rweight-FTINY)
 			continue;
 		/*
 		 *  Ambient radius test.
 		 */
-		e1 = 0.0;
-		for (j = 0; j < 3; j++) {
-			d = av->pos[j] - r->rop[j];
-			e1 += d * d;
-		}
+		d = av->pos[0] - r->rop[0];
+		e1 = d * d;
+		d = av->pos[1] - r->rop[1];
+		e1 += d * d;
+		d = av->pos[2] - r->rop[2];
+		e1 += d * d;
 		e1 /= av->rad * av->rad;
 		if (e1 > ambacc*ambacc*1.21)
 			continue;
 		/*
-		 *  Normal direction test.
+		 *  Direction test using closest normal.
 		 */
-		e2 = (1.0 - DOT(av->dir, r->ron)) * r->rweight;
+		d = DOT(av->dir, r->ron);
+		if (rn != r->ron) {
+			rn_dot = DOT(av->dir, rn);
+			if (rn_dot > 1.0-FTINY)
+				rn_dot = 1.0-FTINY;
+			if (rn_dot >= d-FTINY) {
+				d = rn_dot;
+				rn_dot = -2.0;
+			}
+		}
+		e2 = (1.0 - d) * r->rweight;
 		if (e2 < 0.0) e2 = 0.0;
 		if (e1 + e2 > ambacc*ambacc*1.21)
 			continue;
@@ -197,20 +359,29 @@ double  s;
 		for (j = 0; j < 3; j++)
 			d += (r->rop[j] - av->pos[j]) *
 					(av->dir[j] + r->ron[j]);
-		if (d < -minarad)
+		if (d*0.5 < -minarad*ambacc-.001)
 			continue;
 		/*
 		 *  Jittering final test reduces image artifacts.
 		 */
-		wt = sqrt(e1) + sqrt(e2);
-		if (wt > ambacc*(0.9 + 0.2*frandom()))
+		e1 = sqrt(e1);
+		e2 = sqrt(e2);
+		wt = e1 + e2;
+		if (wt > ambacc*(.9+.2*urand(9015+samplendx)))
 			continue;
+		/*
+		 *  Recompute directional error using perturbed normal
+		 */
+		if (rn_dot > 0.0) {
+			e2 = sqrt((1.0 - rn_dot)*r->rweight);
+			wt = e1 + e2;
+		}
 		if (wt <= 1e-3)
 			wt = 1e3;
 		else
 			wt = 1.0 / wt;
 		wsum += wt;
-		copycolor(ct, av->val);
+		extambient(ct, av, r->rop, rn);
 		scalecolor(ct, wt);
 		addcolor(acol, ct);
 	}
@@ -229,243 +400,203 @@ double  s;
 				break;
 		}
 		if (j == 3)
-			wsum += sumambient(acol, r, at->kid+i, ck0, s);
+			wsum += sumambient(acol, r, rn, al, at->kid+i, ck0, s);
 	}
 	return(wsum);
 }
 
 
 double
-makeambient(acol, r)		/* make a new ambient value */
+makeambient(acol, r, rn, al)	/* make a new ambient value */
 COLOR  acol;
 register RAY  *r;
+FVECT  rn;
+int  al;
 {
-	AMBVAL  amb;
-
-	amb.rad = doambient(acol, r);		/* compute ambient */
-	if (amb.rad == 0.0)
+	AMBVAL	amb;
+	FVECT	gp, gd;
+						/* compute weight */
+	amb.weight = pow(AVGREFL, (double)al);
+	if (r->rweight < 0.1*amb.weight)	/* heuristic */
+		amb.weight = r->rweight;
+						/* compute ambient */
+	amb.rad = doambient(acol, r, amb.weight, gp, gd);
+	if (amb.rad <= FTINY)
 		return(0.0);
 						/* store it */
 	VCOPY(amb.pos, r->rop);
 	VCOPY(amb.dir, r->ron);
-	amb.lvl = r->rlvl;
-	amb.weight = r->rweight;
+	amb.lvl = al;
 	copycolor(amb.val, acol);
+	VCOPY(amb.gpos, gp);
+	VCOPY(amb.gdir, gd);
 						/* insert into tree */
-	avinsert(&amb, &atrunk, thescene.cuorg, thescene.cusize);
-	avsave(&amb);				/* write to file */
+	avsave(&amb);				/* and save to file */
+	if (rn != r->ron)
+		extambient(acol, &amb, r->rop, rn);	/* texture */
 	return(amb.rad);
 }
 
 
-double
-doambient(acol, r)			/* compute ambient component */
-COLOR  acol;
-register RAY  *r;
+void
+extambient(cr, ap, pv, nv)		/* extrapolate value at pv, nv */
+COLOR  cr;
+register AMBVAL	 *ap;
+FVECT  pv, nv;
 {
-	extern int  ambcmp();
-	extern double  sin(), cos(), sqrt();
-	double  phi, xd, yd, zd;
-	double  b, b2;
-	register AMBSAMP  *div;
-	AMBSAMP  dnew;
-	RAY  ar;
-	FVECT  ux, uy;
-	double  arad;
-	int  ndivs, nt, np, ns, ne, i, j;
-	register int  k;
+	FVECT  v1;
+	register int  i;
+	double	d;
 
-	setcolor(acol, 0.0, 0.0, 0.0);
-					/* set number of divisions */
-	nt = sqrt(ambdiv * r->rweight * 0.5) + 0.5;
-	np = 2 * nt;
-	ndivs = nt * np;
-					/* check first */
-	if (ndivs == 0 || rayorigin(&ar, r, AMBIENT, 0.5) < 0)
-		return(0.0);
-					/* set number of super-samples */
-	ns = ambssamp * r->rweight + 0.5;
-	if (ns > 0) {
-		div = (AMBSAMP *)malloc(ndivs*sizeof(AMBSAMP));
-		if (div == NULL)
-			error(SYSTEM, "out of memory in doambient");
+	d = 1.0;			/* zeroeth order */
+					/* gradient due to translation */
+	for (i = 0; i < 3; i++)
+		d += ap->gpos[i]*(pv[i]-ap->pos[i]);
+					/* gradient due to rotation */
+	VCROSS(v1, ap->dir, nv);
+	d += DOT(ap->gdir, v1);
+	if (d <= 0.0) {
+		setcolor(cr, 0.0, 0.0, 0.0);
+		return;
 	}
-					/* make axes */
-	uy[0] = uy[1] = uy[2] = 0.0;
-	for (k = 0; k < 3; k++)
-		if (r->ron[k] < 0.6 && r->ron[k] > -0.6)
-			break;
-	uy[k] = 1.0;
-	fcross(ux, r->ron, uy);
-	normalize(ux);
-	fcross(uy, ux, r->ron);
-						/* sample divisions */
-	arad = 0.0;
-	ne = 0;
-	for (i = 0; i < nt; i++)
-		for (j = 0; j < np; j++) {
-			rayorigin(&ar, r, AMBIENT, 0.5);	/* pretested */
-			zd = sqrt((i+frandom())/nt);
-			phi = 2.0*PI * (j+frandom())/np;
-			xd = cos(phi) * zd;
-			yd = sin(phi) * zd;
-			zd = sqrt(1.0 - zd*zd);
-			for (k = 0; k < 3; k++)
-				ar.rdir[k] = xd*ux[k]+yd*uy[k]+zd*r->ron[k];
-			rayvalue(&ar);
-			if (ar.rot < FHUGE)
-				arad += 1.0 / ar.rot;
-			if (ns > 0) {			/* save division */
-				div[ne].k = 0.0;
-				copycolor(div[ne].v, ar.rcol);
-				div[ne].n = 0;
-				div[ne].t = i; div[ne].p = j;
-							/* sum errors */
-				b = bright(ar.rcol);
-				if (i > 0) {		/* from above */
-					b2 = bright(div[ne-np].v) - b;
-					b2 *= b2 * 0.25;
-					div[ne].k += b2;
-					div[ne].n++;
-					div[ne-np].k += b2;
-					div[ne-np].n++;
-				}
-				if (j > 0) {		/* from behind */
-					b2 = bright(div[ne-1].v) - b;
-					b2 *= b2 * 0.25;
-					div[ne].k += b2;
-					div[ne].n++;
-					div[ne-1].k += b2;
-					div[ne-1].n++;
-				}
-				if (j == np-1) {	/* around */
-					b2 = bright(div[ne-(np-1)].v) - b;
-					b2 *= b2 * 0.25;
-					div[ne].k += b2;
-					div[ne].n++;
-					div[ne-(np-1)].k += b2;
-					div[ne-(np-1)].n++;
-				}
-				ne++;
-			} else
-				addcolor(acol, ar.rcol);
-		}
-	for (k = 0; k < ne; k++) {		/* compute errors */
-		if (div[k].n > 1)
-			div[k].k /= div[k].n;
-		div[k].n = 1;
-	}
-						/* sort the divisions */
-	qsort(div, ne, sizeof(AMBSAMP), ambcmp);
-						/* skim excess */
-	while (ne > ns) {
-		ne--;
-		addcolor(acol, div[ne].v);
-	}
-						/* super-sample */
-	for (i = ns; i > 0; i--) {
-		rayorigin(&ar, r, AMBIENT, 0.5);	/* pretested */
-		zd = sqrt((div[0].t+frandom())/nt);
-		phi = 2.0*PI * (div[0].p+frandom())/np;
-		xd = cos(phi) * zd;
-		yd = sin(phi) * zd;
-		zd = sqrt(1.0 - zd*zd);
-		for (k = 0; k < 3; k++)
-			ar.rdir[k] = xd*ux[k]+yd*uy[k]+zd*r->ron[k];
-		rayvalue(&ar);
-		if (ar.rot < FHUGE)
-			arad += 1.0 / ar.rot;
-						/* recompute error */
-		copycolor(dnew.v, div[0].v);
-		addcolor(dnew.v, ar.rcol);
-		dnew.n = div[0].n + 1;
-		dnew.t = div[0].t; dnew.p = div[0].p;
-		b2 = bright(dnew.v)/dnew.n - bright(ar.rcol);
-		b2 = b2*b2 + div[0].k*(div[0].n*div[0].n);
-		dnew.k = b2/(dnew.n*dnew.n);
-						/* reinsert */
-		for (k = 0; k < ne-1 && dnew.k < div[k+1].k; k++)
-			bcopy(&div[k+1], &div[k], sizeof(AMBSAMP));
-		bcopy(&dnew, &div[k], sizeof(AMBSAMP));
-
-		if (ne >= i) {		/* extract darkest division */
-			ne--;
-			if (div[ne].n > 1)
-				scalecolor(div[ne].v, 1.0/div[ne].n);
-			addcolor(acol, div[ne].v);
-		}
-	}
-	scalecolor(acol, 1.0/ndivs);
-	if (arad <= FTINY)
-		arad = FHUGE;
-	else
-		arad = (ndivs+ns) / arad / sqrt(r->rweight);
-	if (arad > maxarad)
-		arad = maxarad;
-	else if (arad < minarad)
-		arad = minarad;
-	if (ns > 0)
-		free((char *)div);
-	return(arad);
+	copycolor(cr, ap->val);
+	scalecolor(cr, d);
 }
 
 
-static int
-ambcmp(d1, d2)				/* decreasing order */
-AMBSAMP  *d1, *d2;
+static void
+initambfile(creat)		/* initialize ambient file */
+int  creat;
 {
-	if (d1->k < d2->k)
-		return(1);
-	if (d1->k > d2->k)
-		return(-1);
-	return(0);
+	extern char  *progname, *octname;
+	static char  *mybuf = NULL;
+
+#ifdef	F_SETLKW
+	aflock(creat ? F_WRLCK : F_RDLCK);
+#endif
+	SET_FILE_BINARY(ambfp);
+	if (mybuf == NULL)
+		mybuf = (char *)bmalloc(BUFSIZ+8);
+	setbuf(ambfp, mybuf);
+	if (creat) {			/* new file */
+		newheader("RADIANCE", ambfp);
+		fprintf(ambfp, "%s -av %g %g %g -aw %d -ab %d -aa %g ",
+				progname, colval(ambval,RED),
+				colval(ambval,GRN), colval(ambval,BLU),
+				ambvwt, ambounce, ambacc);
+		fprintf(ambfp, "-ad %d -as %d -ar %d ",
+				ambdiv, ambssamp, ambres);
+		if (octname != NULL)
+			printargs(1, &octname, ambfp);
+		else
+			fputc('\n', ambfp);
+		fprintf(ambfp, "SOFTWARE= %s\n", VersionID);
+		fputnow(ambfp);
+		fputformat(AMBFMT, ambfp);
+		putc('\n', ambfp);
+		putambmagic(ambfp);
+	} else if (checkheader(ambfp, AMBFMT, NULL) < 0 || !hasambmagic(ambfp))
+		error(USER, "bad ambient file");
 }
 
 
-static
-avsave(av)				/* save an ambient value */
-AMBVAL  *av;
+static void
+avsave(av)				/* insert and save an ambient value */
+AMBVAL	*av;
 {
-#ifdef  AMBFLUSH
-	static int  nunflshed = 0;
-#endif
+	avinsert(avstore(av));
 	if (ambfp == NULL)
 		return;
-	if (fwrite(av, sizeof(AMBVAL), 1, ambfp) != 1)
+	if (writambval(av, ambfp) < 0)
 		goto writerr;
-#ifdef  AMBFLUSH
-	if (++nunflshed >= AMBFLUSH) {
-		if (fflush(ambfp) == EOF)
+	if (++nunflshed >= AMBFLUSH)
+		if (ambsync() == EOF)
 			goto writerr;
-		nunflshed = 0;
-	}
-#endif
 	return;
 writerr:
-	error(SYSTEM, "error writing ambient file");
+	error(SYSTEM, "error writing to ambient file");
 }
 
 
-static
-avinsert(aval, at, c0, s)		/* insert ambient value in a tree */
-AMBVAL  *aval;
-register AMBTREE  *at;
-FVECT  c0;
-double  s;
+static AMBVAL *
+avstore(aval)				/* allocate memory and store aval */
+register AMBVAL  *aval;
 {
+	register AMBVAL  *av;
+	double	d;
+
+	if ((av = newambval()) == NULL)
+		error(SYSTEM, "out of memory in avstore");
+	copystruct(av, aval);
+	av->latick = ambclock;
+	av->next = NULL;
+	nambvals++;
+	d = bright(av->val);
+	if (d > FTINY) {		/* add to log sum for averaging */
+		avsum += log(d);
+		navsum++;
+	}
+	return(av);
+}
+
+
+#define ATALLOCSZ	512		/* #/8 trees to allocate at once */
+
+static AMBTREE  *atfreelist = NULL;	/* free ambient tree structures */
+
+
+static AMBTREE *
+newambtree()				/* allocate 8 ambient tree structs */
+{
+	register AMBTREE  *atp, *upperlim;
+
+	if (atfreelist == NULL) {	/* get more nodes */
+		atfreelist = (AMBTREE *)malloc(ATALLOCSZ*8*sizeof(AMBTREE));
+		if (atfreelist == NULL)
+			return(NULL);
+					/* link new free list */
+		upperlim = atfreelist + 8*(ATALLOCSZ-1);
+		for (atp = atfreelist; atp < upperlim; atp += 8)
+			atp->kid = atp + 8;
+		atp->kid = NULL;
+	}
+	atp = atfreelist;
+	atfreelist = atp->kid;
+	bzero((char *)atp, 8*sizeof(AMBTREE));
+	return(atp);
+}
+
+
+static void
+freeambtree(atp)			/* free 8 ambient tree structs */
+AMBTREE  *atp;
+{
+	atp->kid = atfreelist;
+	atfreelist = atp;
+}
+
+
+static void
+avinsert(av)				/* insert ambient value in our tree */
+register AMBVAL	 *av;
+{
+	register AMBTREE  *at;
+	register AMBVAL  *ap;
+	AMBVAL  avh;
 	FVECT  ck0;
+	double	s;
 	int  branch;
-	register AMBVAL  *av;
 	register int  i;
 
-	if ((av = newambval()) == NULL)
-		goto memerr;
-	bcopy(aval, av, sizeof(AMBVAL));
-	VCOPY(ck0, c0);
+	if (av->rad <= FTINY)
+		error(CONSISTENCY, "zero ambient radius in avinsert");
+	at = &atrunk;
+	VCOPY(ck0, thescene.cuorg);
+	s = thescene.cusize;
 	while (s*(OCTSCALE/2) > av->rad*ambacc) {
 		if (at->kid == NULL)
 			if ((at->kid = newambtree()) == NULL)
-				goto memerr;
+				error(SYSTEM, "out of memory in avinsert");
 		s *= 0.5;
 		branch = 0;
 		for (i = 0; i < 3; i++)
@@ -475,9 +606,281 @@ double  s;
 			}
 		at = at->kid + branch;
 	}
-	av->next = at->alist;
-	at->alist = av;
-	return;
-memerr:
-	error(SYSTEM, "out of memory in avinsert");
+	avh.next = at->alist;		/* order by increasing level */
+	for (ap = &avh; ap->next != NULL; ap = ap->next)
+		if (ap->next->lvl >= av->lvl)
+			break;
+	av->next = ap->next;
+	ap->next = av;
+	at->alist = avh.next;
 }
+
+
+static void
+unloadatree(at, f)			/* unload an ambient value tree */
+register AMBTREE  *at;
+void	(*f)();
+{
+	register AMBVAL  *av;
+	register int  i;
+					/* transfer values at this node */
+	for (av = at->alist; av != NULL; av = at->alist) {
+		at->alist = av->next;
+		(*f)(av);
+	}
+	if (at->kid == NULL)
+		return;
+	for (i = 0; i < 8; i++)		/* transfer and free children */
+		unloadatree(at->kid+i, f);
+	freeambtree(at->kid);
+	at->kid = NULL;
+}
+
+
+static struct avl {
+	AMBVAL	*p;
+	unsigned long	t;
+}	*avlist1;			/* ambient value list with ticks */
+static AMBVAL	**avlist2;		/* memory positions for sorting */
+static int	i_avlist;		/* index for lists */
+
+
+static int
+av2list(av)
+register AMBVAL	*av;
+{
+#ifdef DEBUG
+	if (i_avlist >= nambvals)
+		error(CONSISTENCY, "too many ambient values in av2list1");
+#endif
+	avlist1[i_avlist].p = avlist2[i_avlist] = av;
+	avlist1[i_avlist++].t = av->latick;
+}
+
+
+static int
+alatcmp(av1, av2)			/* compare ambient values for MRA */
+struct avl	*av1, *av2;
+{
+	register long  lc = av2->t - av1->t;
+	return(lc<0 ? -1 : lc>0 ? 1 : 0);
+}
+
+
+/* GW NOTE 2002/10/3:
+ * I used to compare AMBVAL pointers, but found that this was the
+ * cause of a serious consistency error with gcc, since the optimizer
+ * uses some dangerous trick in pointer subtraction that
+ * assumes pointers differ by exact struct size increments.
+ */
+static int
+aposcmp(avp1, avp2)			/* compare ambient value positions */
+const void	*avp1, *avp2;
+{
+	register long	diff = *(char * const *)avp1 - *(char * const *)avp2;
+	if (diff < 0)
+		return(-1);
+	return(diff > 0);
+}
+
+#if 1
+static int
+avlmemi(avaddr)				/* find list position from address */
+AMBVAL	*avaddr;
+{
+	register AMBVAL  **avlpp;
+
+	avlpp = (AMBVAL **)bsearch((char *)&avaddr, (char *)avlist2,
+			nambvals, sizeof(AMBVAL *), aposcmp);
+	if (avlpp == NULL)
+		error(CONSISTENCY, "address not found in avlmemi");
+	return(avlpp - avlist2);
+}
+#else
+#define avlmemi(avaddr)	((AMBVAL **)bsearch((char *)&avaddr,(char *)avlist2, \
+				nambvals,sizeof(AMBVAL *),aposcmp) - avlist2)
+#endif
+
+
+static void
+sortambvals(always)			/* resort ambient values */
+int	always;
+{
+	AMBTREE  oldatrunk;
+	AMBVAL	tav, *tap, *pnext;
+	register int	i, j;
+					/* see if it's time yet */
+	if (!always && (ambclock++ < lastsort+sortintvl ||
+			nambvals < SORT_THRESH))
+		return;
+	/*
+	 * The idea here is to minimize memory thrashing
+	 * in VM systems by improving reference locality.
+	 * We do this by periodically sorting our stored ambient
+	 * values in memory in order of most recently to least
+	 * recently accessed.  This ordering was chosen so that new
+	 * ambient values (which tend to be less important) go into
+	 * higher memory with the infrequently accessed values.
+	 *	Since we expect our values to need sorting less
+	 * frequently as the process continues, we double our
+	 * waiting interval after each call.
+	 * 	This routine is also called by setambacc() with
+	 * the "always" parameter set to 1 so that the ambient
+	 * tree will be rebuilt with the new accuracy parameter.
+	 */
+	if (tracktime) {		/* allocate pointer arrays to sort */
+		avlist2 = (AMBVAL **)malloc(nambvals*sizeof(AMBVAL *));
+		avlist1 = (struct avl *)malloc(nambvals*sizeof(struct avl));
+	} else {
+		avlist2 = NULL;
+		avlist1 = NULL;
+	}
+	if (avlist1 == NULL) {		/* no time tracking -- rebuild tree? */
+		if (avlist2 != NULL)
+			free((void *)avlist2);
+		if (always) {		/* rebuild without sorting */
+			copystruct(&oldatrunk, &atrunk);
+			atrunk.alist = NULL;
+			atrunk.kid = NULL;
+			unloadatree(&oldatrunk, avinsert);
+		}
+	} else {			/* sort memory by last access time */
+		/*
+		 * Sorting memory is tricky because it isn't contiguous.
+		 * We have to sort an array of pointers by MRA and also
+		 * by memory position.  We then copy values in "loops"
+		 * to minimize memory hits.  Nevertheless, we will visit
+		 * everyone at least twice, and this is an expensive process
+		 * when we're thrashing, which is when we need to do it.
+		 */
+#ifdef DEBUG
+		sprintf(errmsg, "sorting %u ambient values at ambclock=%lu...",
+				nambvals, ambclock);
+		eputs(errmsg);
+#endif
+		i_avlist = 0;
+		unloadatree(&atrunk, av2list);	/* empty current tree */
+#ifdef DEBUG
+		if (i_avlist < nambvals)
+			error(CONSISTENCY, "missing ambient values in sortambvals");
+#endif
+		qsort((char *)avlist1, nambvals, sizeof(struct avl), alatcmp);
+		qsort((char *)avlist2, nambvals, sizeof(AMBVAL *), aposcmp);
+		for (i = 0; i < nambvals; i++) {
+			if (avlist1[i].p == NULL)
+				continue;
+			tap = avlist2[i];
+			copystruct(&tav, tap);
+			for (j = i; (pnext = avlist1[j].p) != tap;
+					j = avlmemi(pnext)) {
+				copystruct(avlist2[j], pnext);
+				avinsert(avlist2[j]);
+				avlist1[j].p = NULL;
+			}
+			copystruct(avlist2[j], &tav);
+			avinsert(avlist2[j]);
+			avlist1[j].p = NULL;
+		}
+		free((void *)avlist1);
+		free((void *)avlist2);
+						/* compute new sort interval */
+		sortintvl = ambclock - lastsort;
+		if (sortintvl >= MAX_SORT_INTVL/2)
+			sortintvl = MAX_SORT_INTVL;
+		else
+			sortintvl <<= 1;	/* wait twice as long next */
+#ifdef DEBUG
+		eputs("done\n");
+#endif
+	}
+	if (ambclock >= MAXACLOCK)
+		ambclock = MAXACLOCK/2;
+	lastsort = ambclock;
+}
+
+
+#ifdef	F_SETLKW
+
+static void
+aflock(typ)			/* lock/unlock ambient file */
+int  typ;
+{
+	static struct flock  fls;	/* static so initialized to zeroes */
+
+	fls.l_type = typ;
+	if (fcntl(fileno(ambfp), F_SETLKW, &fls) < 0)
+		error(SYSTEM, "cannot (un)lock ambient file");
+}
+
+
+int
+ambsync()			/* synchronize ambient file */
+{
+	long  flen;
+	AMBVAL	avs;
+	register int  n;
+
+	if (nunflshed == 0)
+		return(0);
+	if (lastpos < 0)	/* initializing (locked in initambfile) */
+		goto syncend;
+				/* gain exclusive access */
+	aflock(F_WRLCK);
+				/* see if file has grown */
+	if ((flen = lseek(fileno(ambfp), (off_t)0L, 2)) < 0)
+		goto seekerr;
+	if (n = flen - lastpos) {		/* file has grown */
+		if (ambinp == NULL) {		/* use duplicate filedes */
+			ambinp = fdopen(dup(fileno(ambfp)), "r");
+			if (ambinp == NULL)
+				error(SYSTEM, "fdopen failed in ambsync");
+		}
+		if (fseek(ambinp, lastpos, 0) < 0)
+			goto seekerr;
+		while (n >= AMBVALSIZ) {	/* load contributed values */
+			if (!readambval(&avs, ambinp)) {
+				sprintf(errmsg,
+			"ambient file \"%s\" corrupted near character %ld",
+						ambfile, flen - n);
+				error(WARNING, errmsg);
+				break;
+			}
+			avinsert(avstore(&avs));
+			n -= AMBVALSIZ;
+		}
+		/*** seek always as safety measure
+		if (n) ***/			/* alignment */
+			if (lseek(fileno(ambfp), (off_t)(flen-n), 0) < 0)
+				goto seekerr;
+	}
+#ifdef  DEBUG
+	if (ambfp->_ptr - ambfp->_base != nunflshed*AMBVALSIZ) {
+		sprintf(errmsg, "ambient file buffer at %d rather than %d",
+				ambfp->_ptr - ambfp->_base,
+				nunflshed*AMBVALSIZ);
+		error(CONSISTENCY, errmsg);
+	}
+#endif
+syncend:
+	n = fflush(ambfp);			/* calls write() at last */
+	if ((lastpos = lseek(fileno(ambfp), (off_t)0L, 1)) < 0)
+		goto seekerr;
+	aflock(F_UNLCK);			/* release file */
+	nunflshed = 0;
+	return(n);
+seekerr:
+	error(SYSTEM, "seek failed in ambsync");
+}
+
+#else
+
+int
+ambsync()			/* flush ambient file */
+{
+	if (nunflshed == 0)
+		return(0);
+	nunflshed = 0;
+	return(fflush(ambfp));
+}
+
+#endif