--- ray/src/rt/ambient.c	1997/01/02 09:37:13	2.42
+++ ray/src/rt/ambient.c	2023/11/17 20:02:07	2.119
@@ -1,37 +1,30 @@
-/* Copyright (c) 1996 Regents of the University of California */
-
-#ifndef lint
-static char SCCSid[] = "$SunId$ LBL";
-#endif
-
+static const char	RCSid[] = "$Id: ambient.c,v 2.119 2023/11/17 20:02:07 greg Exp $";
 /*
  *  ambient.c - routines dealing with ambient (inter-reflected) component.
+ *
+ *  Declarations of external symbols in ambient.h
  */
 
-#include  "ray.h"
+#include "copyright.h"
 
-#include  "octree.h"
+#include <string.h>
 
+#include  "platform.h"
+#include  "ray.h"
 #include  "otypes.h"
-
+#include  "otspecial.h"
+#include  "resolu.h"
 #include  "ambient.h"
-
 #include  "random.h"
+#include  "pmapamb.h"
 
 #ifndef  OCTSCALE
 #define	 OCTSCALE	1.0	/* ceil((valid rad.)/(cube size)) */
 #endif
 
-typedef struct ambtree {
-	AMBVAL	*alist;		/* ambient value list */
-	struct ambtree	*kid;	/* 8 child nodes */
-}  AMBTREE;			/* ambient octree */
-
-extern CUBE  thescene;		/* contains space boundaries */
-
-extern char  *shm_boundary;	/* memory sharing boundary */
-
-#define	 MAXASET	511	/* maximum number of elements in ambient set */
+#ifndef  MAXASET
+#define	 MAXASET	4095	/* maximum number of elements in ambient set */
+#endif
 OBJECT	ambset[MAXASET+1]={0};	/* ambient include/exclude set */
 
 double	maxarad;		/* maximum ambient radius */
@@ -42,134 +35,184 @@ static AMBTREE	atrunk;		/* our ambient trunk node */
 static FILE  *ambfp = NULL;	/* ambient file pointer */
 static int  nunflshed = 0;	/* number of unflushed ambient values */
 
-#ifndef SORT_THRESH
-#ifdef BIGMEM
-#define SORT_THRESH	((9L<<20)/sizeof(AMBVAL))
-#else
-#define SORT_THRESH	((3L<<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
-
-static COLOR  avsum = BLKCOLOR;		/* computed ambient value sum */
+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 */
 
-#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.
-	 */
-#define tracktime	(shm_boundary == NULL || nambvals > 3*nambshare)
-
 #define	 AMBFLUSH	(BUFSIZ/AMBVALSIZ)
 
-#define	 newambval()	(AMBVAL *)bmalloc(sizeof(AMBVAL))
+#define  AVSIZE		(sizeof(AMBVAL)-sizeof(SCOLOR)+sizeof(COLORV)*NCSAMP)
+#define	 newambval()	(AMBVAL *)malloc(AVSIZE)
 
-extern long  ftell(), lseek();
-static int  initambfile(), avsave(), avinsert(), sortambvals(), avlmemi();
-static AMBVAL  *avstore();
+#define  tfunc(x0, x, x1)	(((x)-(x0))/((x1)-(x0)))
+
+static void initambfile(int cre8);
+static void avsave(AMBVAL *av);
+static AMBVAL *avstore(AMBVAL  *aval);
+static AMBTREE *newambtree(void);
+static void freeambtree(AMBTREE  *atp);
+
+typedef void unloadtf_t(AMBVAL *);
+static unloadtf_t avinsert;
+static unloadtf_t avfree;
+static void unloadatree(AMBTREE  *at, unloadtf_t *f);
+
+static void sortambvals(void);
+
+static int	plugaleak(RAY *r, AMBVAL *ap, FVECT anorm, double ang);
+static double	sumambient(SCOLOR acol, RAY *r, FVECT rn, int al,
+				AMBTREE *at, FVECT c0, double s);
+static int	makeambient(SCOLOR acol, RAY *r, FVECT rn, int al);
+static int	extambient(SCOLOR cr, AMBVAL *ap, FVECT pv, FVECT nv,
+				FVECT uvw[3]);
+
 #ifdef  F_SETLKW
-static  aflock();
+static void aflock(int  typ);
 #endif
 
 
-setambres(ar)				/* set ambient resolution */
-int  ar;
+void
+setambres(				/* 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;
+		maxarad = thescene.cusize*0.2;
 	} else {
 		minarad = thescene.cusize / ar;
-		maxarad = 64 * minarad;			/* heuristic */
-		if (maxarad > thescene.cusize / 2.0)
-			maxarad = thescene.cusize / 2.0;
+		maxarad = 64.0 * minarad;		/* heuristic */
+		if (maxarad > thescene.cusize*0.2)
+			maxarad = thescene.cusize*0.2;
 	}
 	if (minarad <= FTINY)
-		minarad = 10*FTINY;
+		minarad = 10.0*FTINY;
 	if (maxarad <= minarad)
-		maxarad = 64 * minarad;
+		maxarad = 64.0 * minarad;
 }
 
 
-setambacc(newa)				/* set ambient accuracy */
-double  newa;
+void
+setambacc(				/* set ambient accuracy */
+	double  newa
+)
 {
-	double  ambdiff;
-
-	if (newa < 0.0)
-		newa = 0.0;
-	ambdiff = fabs(newa - ambacc);
-	if (ambdiff >= .01 && (ambacc = newa) > FTINY && nambvals > 0)
-		sortambvals(1);			/* rebuild tree */
+	static double	olda;		/* remember previous setting here */
+	
+	newa *= (newa > 0);
+	if (fabs(newa - olda) >= .05*(newa + olda)) {
+		ambacc = newa;
+		if (ambacc > FTINY && nambvals > 0)
+			sortambvals();		/* rebuild tree */
+	}
 }
 
 
-setambient(afile)			/* initialize calculation */
-char  *afile;
+void
+setambient(void)				/* initialize calculation */
 {
-	long  pos, flen;
+	int	readonly = 0;
+	long	flen;
 	AMBVAL	amb;
+						/* make sure we're fresh */
+	ambdone();
 						/* init ambient limits */
 	setambres(ambres);
 	setambacc(ambacc);
-	if (afile == NULL)
+	if (ambfile == NULL || !ambfile[0])
 		return;
 	if (ambacc <= FTINY) {
 		sprintf(errmsg, "zero ambient accuracy so \"%s\" not opened",
-				afile);
+				ambfile);
 		error(WARNING, errmsg);
 		return;
 	}
 						/* open ambient file */
-	if ((ambfp = fopen(afile, "r+")) != NULL) {
-		initambfile(0);
-		pos = ftell(ambfp);
+	if ((ambfp = fopen(ambfile, "r+")) == NULL)
+		readonly = (ambfp = fopen(ambfile, "r")) != NULL;
+	if (ambfp != NULL) {
+		initambfile(0);			/* file exists */
+		lastpos = ftell(ambfp);
 		while (readambval(&amb, ambfp))
-			avinsert(avstore(&amb));
-						/* align */
-		pos += (long)nambvals*AMBVALSIZ;
-		flen = lseek(fileno(ambfp), 0L, 2);
-		if (flen != pos) {
+			avstore(&amb);
+		nambshare = nambvals;		/* share loaded values */
+		if (readonly) {
 			sprintf(errmsg,
+				"loaded %u values from read-only ambient file",
+					nambvals);
+			error(WARNING, errmsg);
+			fclose(ambfp);		/* close file so no writes */
+			ambfp = NULL;
+			return;			/* avoid ambsync() */
+		}
+						/* align file pointer */
+		lastpos += (long)nambvals*AMBVALSIZ;
+		flen = lseek(fileno(ambfp), (off_t)0, SEEK_END);
+		if (flen != lastpos) {
+			sprintf(errmsg,
 			"ignoring last %ld values in ambient file (corrupted)",
-					(flen - pos)/AMBVALSIZ);
+					(flen - lastpos)/AMBVALSIZ);
 			error(WARNING, errmsg);
-			fseek(ambfp, pos, 0);
-			ftruncate(fileno(ambfp), pos);
+			fseek(ambfp, lastpos, SEEK_SET);
+			ftruncate(fileno(ambfp), (off_t)lastpos);
 		}
-		nambshare = nambvals;
-	} else if ((ambfp = fopen(afile, "w+")) != NULL)
-		initambfile(1);
-	else {
-		sprintf(errmsg, "cannot open ambient file \"%s\"", afile);
+	} else if ((ambfp = fopen(ambfile, "w+")) != NULL) {
+		initambfile(1);			/* else create new file */
+		fflush(ambfp);
+		lastpos = ftell(ambfp);
+	} else {
+		sprintf(errmsg, "cannot open ambient file \"%s\"", ambfile);
 		error(SYSTEM, errmsg);
 	}
-	nunflshed++;	/* lie */
-	ambsync();
+#ifdef  F_SETLKW
+	aflock(F_UNLCK);			/* release file */
+#endif
 }
 
 
-ambnotify(obj)			/* record new modifier */
-OBJECT	obj;
+void
+ambdone(void)			/* 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, avfree);
+					/* reset state variables */
+	avsum = 0.;
+	navsum = 0;
+	nambvals = 0;
+	nambshare = 0;
+}
+
+
+void
+ambnotify(			/* record new modifier */
+	OBJECT	obj
+)
+{
 	static int  hitlimit = 0;
-	register OBJREC	 *o = objptr(obj);
-	register char  **amblp;
+	OBJREC	 *o;
+	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++)
@@ -185,14 +228,37 @@ OBJECT	obj;
 }
 
 
-ambient(acol, r, nrm)		/* compute ambient component for ray */
-COLOR  acol;
-register RAY  *r;
-FVECT  nrm;
+void
+multambient(		/* compute ambient component & multiply by coef. */
+	SCOLOR  aval,
+	RAY  *r,
+	FVECT  nrm
+)
 {
+	static double  logAvgAbsorp = 1;
 	static int  rdepth = 0;			/* ambient recursion */
-	double	d;
+	SCOLOR	acol, caustic;
+	int	i, ok;
+	double	d, l;
 
+	/* PMAP: Factor in ambient from photon map, if enabled and ray is
+	 * ambient. Return as all ambient components accounted for, else
+	 * continue. */
+	if (ambPmap(aval, r, rdepth))
+		return;
+
+	if (logAvgAbsorp > 0)			/* exclude in -aw to avoid growth */
+		logAvgAbsorp = log(1.-AVGREFL);
+
+	/* PMAP: Factor in specular-diffuse ambient (caustics) from photon
+	 * map, if enabled and ray is primary, else caustic is zero.  Continue
+	 * with RADIANCE ambient calculation */
+{/* XXX TEMPORARY */
+	COLOR	pmc;
+	scolor_color(pmc, aval);
+	ambPmapCaustic(pmc, r, rdepth);
+	setscolor(caustic, colval(pmc,RED), colval(pmc,GRN), colval(pmc,BLU));
+}
 	if (ambdiv <= 0)			/* no ambient calculation */
 		goto dumbamb;
 						/* check number of bounces */
@@ -204,227 +270,415 @@ FVECT  nrm;
 		goto dumbamb;
 
 	if (ambacc <= FTINY) {			/* no ambient storage */
+		FVECT	uvd[2];
+		float	dgrad[2], *dgp = NULL;
+
+		if (nrm != r->ron && DOT(nrm,r->ron) < 0.9999)
+			dgp = dgrad;		/* compute rotational grad. */
+		copyscolor(acol, aval);
 		rdepth++;
-		d = doambient(acol, r, r->rweight, NULL, NULL);
+		ok = doambient(acol, r, r->rweight,
+				uvd, NULL, NULL, dgp, NULL);
 		rdepth--;
-		if (d <= FTINY)
+		if (!ok)
 			goto dumbamb;
+		if ((ok > 0) & (dgp != NULL)) {	/* apply texture */
+			FVECT	v1;
+			VCROSS(v1, r->ron, nrm);
+			d = 1.0;
+			for (i = 3; i--; )
+				d += v1[i] * (dgp[0]*uvd[0][i] + dgp[1]*uvd[1][i]);
+			if (d >= 0.05)
+				scalescolor(acol, d);
+		}
+		copyscolor(aval, acol);
+
+		/* PMAP: add in caustic */
+		saddscolor(aval, caustic);
 		return;
 	}
-
-	if (tracktime)				/* sort to minimize thrashing */
-		sortambvals(0);
-						/* get ambient value */
-	setcolor(acol, 0.0, 0.0, 0.0);
+						/* interpolate ambient value */
+	scolorblack(acol);
 	d = sumambient(acol, r, nrm, rdepth,
 			&atrunk, thescene.cuorg, thescene.cusize);
+			
 	if (d > FTINY) {
-		scalecolor(acol, 1.0/d);
+		scalescolor(acol, 1.0/d);
+		smultscolor(aval, acol);
+
+		/* PMAP: add in caustic */
+		saddscolor(aval, caustic);
 		return;
 	}
+	
 	rdepth++;				/* need to cache new value */
-	d = makeambient(acol, r, nrm, rdepth-1);
+	ok = makeambient(acol, r, nrm, rdepth-1);
 	rdepth--;
-	if (d > FTINY)
+	
+	if (ok) {
+		smultscolor(aval, acol);	/* computed new value */
+
+		/* PMAP: add in caustic */
+		saddscolor(aval, caustic);
 		return;
+	}
+	
 dumbamb:					/* return global value */
-	copycolor(acol, ambval);
-	if (ambvwt <= 0 | nambvals == 0)
+	if ((ambvwt <= 0) | (navsum == 0)) {
+		smultcolor(aval, ambval);
+		
+		/* PMAP: add in caustic */
+		saddscolor(aval, caustic);
 		return;
-	scalecolor(acol, (double)ambvwt);
-	addcolor(acol, avsum);			/* average in computations */
-	d = 1.0/(ambvwt+nambvals);
-	scalecolor(acol, d);
+	}
+	
+	l = bright(ambval);			/* average in computations */	
+	if (l > FTINY) {
+		d = (log(l)*(double)ambvwt + avsum + logAvgAbsorp*navsum) /
+				(double)(ambvwt + navsum);
+		d = exp(d) / l;
+		scalescolor(aval, d);
+		smultcolor(aval, ambval);	/* apply color of ambval */
+	} else {
+		d = exp( avsum/(double)navsum + logAvgAbsorp );
+		scalescolor(aval, d);		/* neutral color */
+	}
 }
 
 
-double
-sumambient(acol, r, rn, al, at, c0, s)	/* get interpolated ambient value */
-COLOR  acol;
-register RAY  *r;
-FVECT  rn;
-int  al;
-AMBTREE	 *at;
-FVECT  c0;
-double	s;
+/* Plug a potential leak where ambient cache value is occluded */
+static int
+plugaleak(RAY *r, AMBVAL *ap, FVECT anorm, double ang)
 {
-	double	d, e1, e2, wt, wsum;
-	COLOR  ct;
-	FVECT  ck0;
-	int  i;
-	register int  j;
-	register AMBVAL	 *av;
+	const double	cost70sq = 0.1169778;	/* cos(70deg)^2 */
+	RAY		rtst;
+	FVECT		vdif;
+	double		normdot, ndotd, nadotd;
+	double		a, b, c, t[2];
 
-	wsum = 0.0;
-					/* do this node */
+	ang += 2.*PI*(ang < 0);			/* check direction flags */
+	if ( !(ap->corral>>(int)(ang*(16./PI)) & 1) )
+		return(0);
+	/*
+	 * Generate test ray, targeting 20 degrees above sample point plane
+	 * along surface normal from cache position.  This should be high
+	 * enough to miss local geometry we don't really care about.
+	 */
+	VSUB(vdif, ap->pos, r->rop);
+	normdot = DOT(anorm, r->ron);
+	ndotd = DOT(vdif, r->ron);
+	nadotd = DOT(vdif, anorm);
+	a = normdot*normdot - cost70sq;
+	b = 2.0*(normdot*ndotd - nadotd*cost70sq);
+	c = ndotd*ndotd - DOT(vdif,vdif)*cost70sq;
+	if (quadratic(t, a, b, c) != 2)
+		return(1);			/* should rarely happen */
+	if (t[1] <= FTINY)
+		return(0);			/* should fail behind test */
+	rayorigin(&rtst, SHADOW, r, NULL);
+	VSUM(rtst.rdir, vdif, anorm, t[1]);	/* further dist. > plane */
+	rtst.rmax = normalize(rtst.rdir);	/* short ray test */
+	while (localhit(&rtst, &thescene)) {	/* check for occluder */
+		OBJREC	*m = findmaterial(rtst.ro);
+		if (m != NULL && !istransp(m->otype) && !isBSDFproxy(m) &&
+				(rtst.clipset == NULL ||
+					!inset(rtst.clipset, rtst.ro->omod)))
+			return(1);		/* plug light leak */
+		VCOPY(rtst.rorg, rtst.rop);	/* skip invisible surface */
+		rtst.rmax -= rtst.rot;
+		rayclear(&rtst);
+	}
+	return(0);				/* seems we're OK */
+}
+
+
+static double
+sumambient(		/* get interpolated ambient value */
+	SCOLOR  acol,
+	RAY  *r,
+	FVECT  rn,
+	int  al,
+	AMBTREE	 *at,
+	FVECT  c0,
+	double	s
+)
+{			/* initial limit is 10 degrees plus ambacc radians */
+	const double	minangle = 10.0 * PI/180.;
+	double		maxangle = minangle + ambacc;
+	double		wsum = 0.0;
+	FVECT		ck0;
+	int		i, j;
+	AMBVAL		*av;
+
+	if (at->kid != NULL) {		/* sum children first */				
+		s *= 0.5;
+		for (i = 0; i < 8; i++) {
+			for (j = 0; j < 3; j++) {
+				ck0[j] = c0[j];
+				if (1<<j & i)
+					ck0[j] += s;
+				if (r->rop[j] < ck0[j] - OCTSCALE*s)
+					break;
+				if (r->rop[j] > ck0[j] + (1.0+OCTSCALE)*s)
+					break;
+			}
+			if (j == 3)
+				wsum += sumambient(acol, r, rn, al,
+							at->kid+i, ck0, s);
+		}
+					/* good enough? */
+		if ((wsum >= 0.05) & (s*ambacc > minarad))
+			return(wsum);
+	}
+					/* adjust maximum angle */
+	if (at->alist != NULL && (at->alist->lvl <= al) & (r->rweight < 0.6))
+		maxangle = (maxangle - PI/2.)*pow(r->rweight,0.13) + PI/2.;
+					/* sum this node */
 	for (av = at->alist; av != NULL; av = av->next) {
-		if (tracktime)
-			av->latick = ambclock;
+		double	u, v, d, delta_r2, delta_t2;
+		SCOLOR	sct;
+		FVECT	uvw[3];
 		/*
-		 *  Ambient level test.
+		 *  Ambient level test
 		 */
-		if (av->lvl > al)	/* list sorted, so this works */
+		if (av->lvl > al ||	/* list sorted, so this works */
+				(av->lvl == al) & (av->weight < 0.9*r->rweight))
 			break;
-		if (av->weight < r->rweight-FTINY)
+		/*
+		 *  Direction test using unperturbed normal
+		 */
+		decodedir(uvw[2], av->ndir);
+		d = DOT(uvw[2], r->ron);
+		if (d <= 0.0)		/* >= 90 degrees */
 			continue;
+		delta_r2 = 2.0 - 2.0*d;	/* approx. radians^2 */
+		if (delta_r2 >= maxangle*maxangle)
+			continue;
 		/*
-		 *  Ambient radius test.
+		 *  Modified ray behind test
 		 */
-		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)
+		VSUB(ck0, r->rop, av->pos);
+		d = DOT(ck0, uvw[2]);
+		if (d < -minarad*ambacc)
 			continue;
+		d /= av->rad[0];
+		delta_t2 = d*d;
+		if (delta_t2 >= ambacc*ambacc)
+			continue;
 		/*
-		 *  Normal direction test.
+		 *  Elliptical radii test based on Hessian
 		 */
-		e2 = (1.0 - DOT(av->dir, r->ron)) * r->rweight;
-		if (e2 < 0.0) e2 = 0.0;
-		if (e1 + e2 > ambacc*ambacc*1.21)
+		decodedir(uvw[0], av->udir);
+		VCROSS(uvw[1], uvw[2], uvw[0]);
+		d = (u = DOT(ck0, uvw[0])) / av->rad[0];
+		delta_t2 += d*d;
+		d = (v = DOT(ck0, uvw[1])) / av->rad[1];
+		delta_t2 += d*d;
+		if (delta_t2 >= ambacc*ambacc)
 			continue;
 		/*
-		 *  Ray behind test.
+		 *  Test for potential light leak
 		 */
-		d = 0.0;
-		for (j = 0; j < 3; j++)
-			d += (r->rop[j] - av->pos[j]) *
-					(av->dir[j] + r->ron[j]);
-		if (d*0.5 < -minarad*ambacc-.001)
+		if (av->corral && plugaleak(r, av, uvw[2], atan2a(v,u)))
 			continue;
 		/*
-		 *  Jittering final test reduces image artifacts.
+		 *  Extrapolate value and compute final weight (hat function)
 		 */
-		wt = sqrt(e1) + sqrt(e2);
-		if (wt > ambacc*(.9+.2*urand(9015+samplendx)))
+		if (!extambient(sct, av, r->rop, rn, uvw))
 			continue;
-		if (wt <= 1e-3)
-			wt = 1e3;
-		else
-			wt = 1.0 / wt;
-		wsum += wt;
-		extambient(ct, av, r->rop, rn);
-		scalecolor(ct, wt);
-		addcolor(acol, ct);
+		d = tfunc(maxangle, sqrt(delta_r2), 0.0) *
+			tfunc(ambacc, sqrt(delta_t2), 0.0);
+		scalescolor(sct, d);
+		saddscolor(acol, sct);
+		wsum += d;
 	}
-	if (at->kid == NULL)
-		return(wsum);
-					/* do children */
-	s *= 0.5;
-	for (i = 0; i < 8; i++) {
-		for (j = 0; j < 3; j++) {
-			ck0[j] = c0[j];
-			if (1<<j & i)
-				ck0[j] += s;
-			if (r->rop[j] < ck0[j] - OCTSCALE*s)
-				break;
-			if (r->rop[j] > ck0[j] + (1.0+OCTSCALE)*s)
-				break;
-		}
-		if (j == 3)
-			wsum += sumambient(acol, r, rn, al, at->kid+i, ck0, s);
-	}
 	return(wsum);
 }
 
 
-double
-makeambient(acol, r, rn, al)	/* make a new ambient value */
-COLOR  acol;
-register RAY  *r;
-FVECT  rn;
-int  al;
+static int
+makeambient(		/* make a new ambient value for storage */
+	SCOLOR  acol,
+	RAY  *r,
+	FVECT  rn,
+	int  al
+)
 {
 	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;
+	FVECT	uvw[3];
+	int	i;
+
+	amb.weight = 1.0;			/* compute weight */
+	for (i = al; i-- > 0; )
+		amb.weight *= AVGREFL;
+	if (r->rweight < 0.1*amb.weight)	/* heuristic override */
+		amb.weight = 1.25*r->rweight;
+	setscolor(acol, AVGREFL, AVGREFL, AVGREFL);
 						/* compute ambient */
-	amb.rad = doambient(acol, r, amb.weight, gp, gd);
-	if (amb.rad <= FTINY)
-		return(0.0);
-						/* store it */
+	i = doambient(acol, r, amb.weight,
+			uvw, amb.rad, amb.gpos, amb.gdir, &amb.corral);
+	scalescolor(acol, 1./AVGREFL);		/* undo assumed reflectance */
+	if (i <= 0 || amb.rad[0] <= FTINY)	/* no Hessian or zero radius */
+		return(i);
+						/* store value */
 	VCOPY(amb.pos, r->rop);
-	VCOPY(amb.dir, r->ron);
+	amb.ndir = encodedir(r->ron);
+	amb.udir = encodedir(uvw[0]);
 	amb.lvl = al;
-	copycolor(amb.val, acol);
-	VCOPY(amb.gpos, gp);
-	VCOPY(amb.gdir, gd);
+	copyscolor(amb.val, acol);
 						/* insert into tree */
 	avsave(&amb);				/* and save to file */
-	if (rn != r->ron)
-		extambient(acol, &amb, r->rop, rn);	/* texture */
-	return(amb.rad);
+	if (rn != r->ron) {			/* texture */
+		VCOPY(uvw[2], r->ron);
+		extambient(acol, &amb, r->rop, rn, uvw);
+	}
+	return(1);
 }
 
 
-extambient(cr, ap, pv, nv)		/* extrapolate value at pv, nv */
-COLOR  cr;
-register AMBVAL	 *ap;
-FVECT  pv, nv;
+static int
+extambient(		/* extrapolate value at pv, nv */
+	SCOLOR  scr,
+	AMBVAL	 *ap,
+	FVECT  pv,
+	FVECT  nv,
+	FVECT  uvw[3]
+)
 {
-	FVECT  v1, v2;
-	register int  i;
-	double	d;
+	const double	min_d = 0.05;
+	const double	max_d = 20.;
+	static FVECT	my_uvw[3];
+	FVECT		v1;
+	int		i;
+	double		d = 1.0;	/* zeroeth order */
 
-	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 */
-	VCOPY(v1, ap->dir);
-	fcross(v2, v1, nv);
-	d += DOT(ap->gdir, v2);
-	if (d <= 0.0) {
-		setcolor(cr, 0.0, 0.0, 0.0);
-		return;
+	if (uvw == NULL) {		/* need local coordinates? */
+		decodedir(my_uvw[2], ap->ndir);
+		decodedir(my_uvw[0], ap->udir);
+		VCROSS(my_uvw[1], my_uvw[2], my_uvw[0]);
+		uvw = my_uvw;
 	}
-	copycolor(cr, ap->val);
-	scalecolor(cr, d);
+	for (i = 3; i--; )		/* gradient due to translation */
+		d += (pv[i] - ap->pos[i]) *
+			(ap->gpos[0]*uvw[0][i] + ap->gpos[1]*uvw[1][i]);
+
+	VCROSS(v1, uvw[2], nv);		/* gradient due to rotation */
+	for (i = 3; i--; )
+		d += v1[i] * (ap->gdir[0]*uvw[0][i] + ap->gdir[1]*uvw[1][i]);
+	
+	if (d < min_d)			/* clamp min/max scaling */
+		d = min_d;
+	else if (d > max_d)
+		d = max_d;
+	copyscolor(scr, ap->val);
+	scalescolor(scr, d);
+	return(d > min_d);
 }
 
 
-static
-initambfile(creat)		/* initialize ambient file */
-int  creat;
+static void
+avinsert(				/* insert ambient value in our tree */
+	AMBVAL *av
+)
 {
-	extern char  *progname, *octname, VersionID[];
+	AMBTREE  *at;
+	AMBVAL  *ap;
+	AMBVAL  avh;
+	FVECT  ck0;
+	double	s;
+	int  branch;
+	int  i;
 
+	if (av->rad[0] <= FTINY)
+		error(CONSISTENCY, "zero ambient radius in avinsert");
+	at = &atrunk;
+	VCOPY(ck0, thescene.cuorg);
+	s = thescene.cusize;
+	while (s*(OCTSCALE/2) > av->rad[1]*ambacc) {
+		if (at->kid == NULL)
+			if ((at->kid = newambtree()) == NULL)
+				error(SYSTEM, "out of memory in avinsert");
+		s *= 0.5;
+		branch = 0;
+		for (i = 0; i < 3; i++)
+			if (av->pos[i] > ck0[i] + s) {
+				ck0[i] += s;
+				branch |= 1 << i;
+			}
+		at = at->kid + branch;
+	}
+	avh.next = at->alist;		/* order by increasing level */
+	for (ap = &avh; ap->next != NULL; ap = ap->next)
+		if ( ap->next->lvl > av->lvl ||
+				(ap->next->lvl == av->lvl) &
+				(ap->next->weight <= av->weight) )
+			break;
+	av->next = ap->next;
+	ap->next = (AMBVAL*)av;
+	at->alist = avh.next;
+}
+
+
+static void
+initambfile(		/* initialize ambient file */
+	int  cre8
+)
+{
+	extern char  *progname, *octname;
+	static char  *mybuf = NULL;
+
 #ifdef	F_SETLKW
-	aflock(creat ? F_WRLCK : F_RDLCK);
+	aflock(cre8 ? F_WRLCK : F_RDLCK);
 #endif
-#ifdef MSDOS
-	setmode(fileno(ambfp), O_BINARY);
-#endif
-	setbuf(ambfp, bmalloc(BUFSIZ+8));
-	if (creat) {			/* new file */
+	SET_FILE_BINARY(ambfp);
+	if (mybuf == NULL)
+		mybuf = (char *)bmalloc(BUFSIZ+8);
+	setbuf(ambfp, mybuf);
+	if (cre8) {			/* 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 %s\n",
-				ambdiv, ambssamp, ambres,
-				octname==NULL ? "" : octname);
+		fprintf(ambfp, "-ad %d -as %d -ar %d ",
+				ambdiv, ambssamp, ambres);
+		fprintf(ambfp, "-dr %d -ds %g -dt %g -dc %g ", directrelay,
+				srcsizerat, shadthresh, shadcert);
+		fprintf(ambfp, "-ss %g -st %g -lr %d -lw %g ", specjitter,
+				specthresh, maxdepth, minweight);
+		fprintf(ambfp, "-cw %g %g -cs %d ", WLPART[3], WLPART[0], NCSAMP);
+		if (octname != NULL)
+			fputs(octname, ambfp);
+		fputc('\n', ambfp);
 		fprintf(ambfp, "SOFTWARE= %s\n", VersionID);
+		fputnow(ambfp);
+		fputwlsplit(WLPART, ambfp);
+		fputncomp(NCSAMP, ambfp);
 		fputformat(AMBFMT, ambfp);
-		putc('\n', ambfp);
+		fputc('\n', ambfp);
 		putambmagic(ambfp);
-	} else if (checkheader(ambfp, AMBFMT, NULL) < 0 || !hasambmagic(ambfp))
-		error(USER, "bad ambient file");
+	} else if (getheader(ambfp, amb_headline, NULL) < 0 || !hasambmagic(ambfp))
+		error(USER, "bad/incompatible ambient file");
+
+	if ((AMB_CNDX != CNDX) | (AMB_WLPART != WLPART)) {
+		if (setspectrsamp(AMB_CNDX, AMB_WLPART) < 0)
+			error(USER, "bad wavelength sampling in ambient file");
+		if (AMB_CNDX[3] == CNDX[3] && FABSEQ(AMB_WLPART[0],WLPART[0]) &&
+					FABSEQ(AMB_WLPART[3],WLPART[3])) {
+			AMB_CNDX = CNDX;
+			AMB_WLPART = WLPART;		/* just the same */
+		} else
+			error(WARNING, "different ambient file wavelength sampling");
+	}
 }
 
 
-static
-avsave(av)				/* insert and save an ambient value */
-AMBVAL	*av;
+static void
+avsave(				/* insert and save an ambient value */
+	AMBVAL	*av
+)
 {
-	avinsert(avstore(av));
+	avstore(av);
 	if (ambfp == NULL)
 		return;
 	if (writambval(av, ambfp) < 0)
@@ -434,23 +688,29 @@ AMBVAL	*av;
 			goto writerr;
 	return;
 writerr:
-	error(SYSTEM, "error writing ambient file");
+	error(SYSTEM, "error writing to ambient file");
 }
 
 
 static AMBVAL *
-avstore(aval)				/* allocate memory and store aval */
-register AMBVAL  *aval;
+avstore(				/* allocate memory and save aval */
+	AMBVAL  *aval
+)
 {
-	register AMBVAL  *av;
+	AMBVAL  *av;
+	double	d;
 
 	if ((av = newambval()) == NULL)
 		error(SYSTEM, "out of memory in avstore");
-	copystruct(av, aval);
-	av->latick = ambclock;
+	memcpy(av, aval, AVSIZE);	/* AVSIZE <= sizeof(AMBVAL) */
 	av->next = NULL;
-	addcolor(avsum, av->val);	/* add to sum for averaging */
 	nambvals++;
+	d = pbright(av->val);
+	if (d > FTINY) {		/* add to log sum for averaging */
+		avsum += log(d);
+		navsum++;
+	}
+	avinsert(av);			/* insert in our cache tree */
 	return(av);
 }
 
@@ -460,14 +720,13 @@ register AMBVAL  *aval;
 static AMBTREE  *atfreelist = NULL;	/* free ambient tree structures */
 
 
-static
-AMBTREE *
-newambtree()				/* allocate 8 ambient tree structs */
+static AMBTREE *
+newambtree(void)				/* allocate 8 ambient tree structs */
 {
-	register AMBTREE  *atp, *upperlim;
+	AMBTREE  *atp, *upperlim;
 
 	if (atfreelist == NULL) {	/* get more nodes */
-		atfreelist = (AMBTREE *)bmalloc(ATALLOCSZ*8*sizeof(AMBTREE));
+		atfreelist = (AMBTREE *)malloc(ATALLOCSZ*8*sizeof(AMBTREE));
 		if (atfreelist == NULL)
 			return(NULL);
 					/* link new free list */
@@ -478,70 +737,33 @@ newambtree()				/* allocate 8 ambient tree structs */
 	}
 	atp = atfreelist;
 	atfreelist = atp->kid;
-	bzero((char *)atp, 8*sizeof(AMBTREE));
+	memset(atp, 0, 8*sizeof(AMBTREE));
 	return(atp);
 }
 
 
-static
-freeambtree(atp)			/* free 8 ambient tree structs */
-AMBTREE  *atp;
+static void
+freeambtree(			/* free 8 ambient tree structs */
+	AMBTREE  *atp
+)
 {
 	atp->kid = atfreelist;
 	atfreelist = atp;
 }
 
 
-static
-avinsert(av)				/* insert ambient value in our tree */
-register AMBVAL	 *av;
+static void
+unloadatree(			/* unload an ambient value tree */
+	AMBTREE  *at,
+	unloadtf_t *f
+)
 {
-	register AMBTREE  *at;
-	register AMBVAL  *ap;
-	AMBVAL  avh;
-	FVECT  ck0;
-	double	s;
-	int  branch;
-	register int  i;
-
-	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)
-				error(SYSTEM, "out of memory in avinsert");
-		s *= 0.5;
-		branch = 0;
-		for (i = 0; i < 3; i++)
-			if (av->pos[i] > ck0[i] + s) {
-				ck0[i] += s;
-				branch |= 1 << i;
-			}
-		at = at->kid + branch;
-	}
-	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
-unloadatree(at, f)			/* unload an ambient value tree */
-register AMBTREE  *at;
-int	(*f)();
-{
-	register AMBVAL  *av;
-	register int  i;
+	AMBVAL  *av;
+	int  i;
 					/* transfer values at this node */
 	for (av = at->alist; av != NULL; av = at->alist) {
 		at->alist = av->next;
+		av->next = NULL;
 		(*f)(av);
 	}
 	if (at->kid == NULL)
@@ -553,244 +775,102 @@ int	(*f)();
 }
 
 
-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
-av2list(av)
-register AMBVAL	*av;
+static void
+avfree(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;
+	free(av);
 }
 
 
-static int
-alatcmp(av1, av2)			/* compare ambient values for MRA */
-struct avl	*av1, *av2;
+static void
+sortambvals(void)			/* resort ambient values */
 {
-	register long  lc = av2->t - av1->t;
-	return(lc<0 ? -1 : lc>0 ? 1 : 0);
-}
+	AMBTREE  oldatrunk = atrunk;
 
-
-static int
-aposcmp(avp1, avp2)			/* compare ambient value positions */
-AMBVAL	**avp1, **avp2;
-{
-	return(*avp1 - *avp2);
+	atrunk.alist = NULL;
+	atrunk.kid = NULL;
+	unloadatree(&oldatrunk, avinsert);
 }
 
 
-#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
-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((char *)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((char *)avlist1);
-		free((char *)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
-aflock(typ)			/* lock/unlock ambient file */
-int  typ;
+static void
+aflock(			/* lock/unlock ambient file */
+	int  typ
+)
 {
 	static struct flock  fls;	/* static so initialized to zeroes */
 
+	if (typ == fls.l_type)		/* already called? */
+		return;
+
 	fls.l_type = typ;
-	if (fcntl(fileno(ambfp), F_SETLKW, &fls) < 0)
-		error(SYSTEM, "cannot (un)lock ambient file");
+	do
+		if (fcntl(fileno(ambfp), F_SETLKW, &fls) != -1)
+			return;
+	while (errno == EINTR);
+	
+	error(SYSTEM, "cannot (un)lock ambient file");
 }
 
 
 int
-ambsync()			/* synchronize ambient file */
+ambsync(void)			/* synchronize ambient file */
 {
-	static FILE  *ambinp = NULL;
-	static long  lastpos = -1;
 	long  flen;
 	AMBVAL	avs;
-	register int  n;
+	int  n;
 
-	if (nunflshed == 0)
+	if (ambfp == NULL)	/* no ambient file? */
 		return(0);
-	if (lastpos < 0)	/* initializing (locked in initambfile) */
-		goto syncend;
-				/* gain exclusive access */
-	aflock(F_WRLCK);
+				/* gain appropriate access */
+	aflock(nunflshed ? F_WRLCK : F_RDLCK);
 				/* see if file has grown */
-	if ((flen = lseek(fileno(ambfp), 0L, 2)) < 0)
+	if ((flen = lseek(fileno(ambfp), (off_t)0, SEEK_END)) < 0)
 		goto seekerr;
-	if (n = flen - lastpos) {		/* file has grown */
-		if (ambinp == NULL) {		/* use duplicate filedes */
-			ambinp = fdopen(dup(fileno(ambfp)), "r");
+	if ((n = flen - lastpos) > 0) {		/* file has grown */
+		if (ambinp == NULL) {		/* get new file pointer */
+			ambinp = fopen(ambfile, "rb");
 			if (ambinp == NULL)
-				error(SYSTEM, "fdopen failed in ambsync");
+				error(SYSTEM, "fopen failed in ambsync");
 		}
-		if (fseek(ambinp, lastpos, 0) < 0)
+		if (fseek(ambinp, lastpos, SEEK_SET) < 0)
 			goto seekerr;
 		while (n >= AMBVALSIZ) {	/* load contributed values */
 			if (!readambval(&avs, ambinp)) {
 				sprintf(errmsg,
-				"ambient file corrupted near character %ld",
-						flen - n);
+			"ambient file \"%s\" corrupted near character %ld",
+						ambfile, flen - n);
 				error(WARNING, errmsg);
 				break;
 			}
-			avinsert(avstore(&avs));
+			avstore(&avs);
 			n -= AMBVALSIZ;
 		}
-		/*** seek always as safety measure
-		if (n) ***/			/* alignment */
-			if (lseek(fileno(ambfp), flen-n, 0) < 0)
-				goto seekerr;
+		lastpos = flen - n;		/* check alignment */
+		if (n && lseek(fileno(ambfp), (off_t)lastpos, SEEK_SET) < 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), 0L, 1)) < 0)
-		goto seekerr;
+	lastpos += (long)nunflshed*AMBVALSIZ;
 	aflock(F_UNLCK);			/* release file */
 	nunflshed = 0;
 	return(n);
 seekerr:
 	error(SYSTEM, "seek failed in ambsync");
+	return(EOF);	/* pro forma return */
 }
 
-#else
+#else	/* ! F_SETLKW */
 
 int
-ambsync()			/* flush ambient file */
+ambsync(void)			/* flush ambient file */
 {
-	if (nunflshed == 0)
+	if (ambfp == NULL)
 		return(0);
 	nunflshed = 0;
 	return(fflush(ambfp));
 }
 
-#endif
+#endif	/* ! F_SETLKW */