--- ray/src/px/pcond3.c	1997/01/29 13:22:13	3.8
+++ ray/src/px/pcond3.c	2003/05/13 17:58:33	3.13
@@ -1,9 +1,6 @@
-/* Copyright (c) 1997 Regents of the University of California */
-
 #ifndef lint
-static char SCCSid[] = "$SunId$ LBL";
+static const char	RCSid[] = "$Id: pcond3.c,v 3.13 2003/05/13 17:58:33 greg Exp $";
 #endif
-
 /*
  * Routines for computing and applying brightness mapping.
  */
@@ -13,11 +10,12 @@ static char SCCSid[] = "$SunId$ LBL";
 
 #define CVRATIO		0.025		/* fraction of samples allowed > env. */
 
-#define exp10(x)	exp(2.302585093*(x))
+#define LN_10		2.30258509299404568402
+#define exp10(x)	exp(LN_10*(x))
 
-float	modhist[HISTRES];		/* modified histogram */
+double	modhist[HISTRES];		/* modified histogram */
 double	mhistot;			/* modified histogram total */
-float	cumf[HISTRES+1];		/* cumulative distribution function */
+double	cumf[HISTRES+1];		/* cumulative distribution function */
 
 
 getfixations(fp)		/* load fixation history list */
@@ -28,7 +26,7 @@ FILE	*fp;
 	int	pos[2];
 	register int	px, py, i;
 				/* initialize our resolution struct */
-	if ((fvres.or=inpres.or)&YMAJOR) {
+	if ((fvres.rt=inpres.rt)&YMAJOR) {
 		fvres.xr = fvxr;
 		fvres.yr = fvyr;
 	} else {
@@ -56,7 +54,7 @@ FILE	*fp;
 				if (nfixations % FIXHUNK == 0) {
 					if (nfixations)
 						fixlst = (short (*)[2])
-							realloc((char *)fixlst,
+							realloc((void *)fixlst,
 							(nfixations+FIXHUNK)*
 							2*sizeof(short));
 					else
@@ -81,6 +79,74 @@ FILE	*fp;
 }
 
 
+gethisto(fp)			/* load precomputed luminance histogram */
+FILE	*fp;
+{
+	double	histo[MAXPREHIST];
+	double	histart, histep;
+	double	l, b, lastb, w;
+	int	n;
+	register int	i;
+					/* load data */
+	for (i = 0; i < MAXPREHIST &&
+			fscanf(fp, "%lf %lf", &b, &histo[i]) == 2; i++) {
+		if (i > 1 && fabs(b - lastb - histep) > .001) {
+			fprintf(stderr,
+				"%s: uneven step size in histogram data\n",
+					progname);
+			exit(1);
+		}
+		if (i == 1)
+			if ((histep = b - (histart = lastb)) <= FTINY) {
+				fprintf(stderr,
+					"%s: illegal step in histogram data\n",
+						progname);
+				exit(1);
+			}
+		lastb = b;
+	}
+	if (i < 2 || !feof(fp)) {
+		fprintf(stderr,
+		"%s: format/length error loading histogram (log10L %f at %d)\n",
+				progname, b, i);
+		exit(1);
+	}
+	n = i;
+	histart *= LN_10;
+	histep *= LN_10;
+					/* find extrema */
+	for (i = 0; i < n && histo[i] <= FTINY; i++)
+		;
+	bwmin = histart + (i-.001)*histep;
+	for (i = n; i-- && histo[i] <= FTINY; )
+		;
+	bwmax = histart + (i+1.001)*histep;
+	if (bwmax > Bl(LMAX))
+		bwmax = Bl(LMAX);
+	if (bwmin < Bl(LMIN))
+		bwmin = Bl(LMIN);
+	else				/* duplicate bottom bin */
+		bwmin = bwmax - (bwmax-bwmin)*HISTRES/(HISTRES-1);
+					/* convert histogram */
+	bwavg = 0.; histot = 0.;
+	for (i = 0; i < HISTRES; i++)
+		bwhist[i] = 0.;
+	for (i = 0, b = histart; i < n; i++, b += histep) {
+		if (b < bwmin+FTINY) continue;
+		if (b >= bwmax-FTINY) break;
+		w = histo[i];
+		bwavg += w*b;
+		bwhist[bwhi(b)] += w;
+		histot += w;
+	}
+	bwavg /= histot;
+	if (bwmin > Bl(LMIN)+FTINY) {	/* add false samples at bottom */
+		bwhist[1] *= 0.5;
+		bwhist[0] += bwhist[1];
+	}
+}
+
+
 double
 centprob(x, y)			/* center-weighting probability function */
 int	x, y;
@@ -98,7 +164,9 @@ comphist()			/* create foveal sampling histogram */
 {
 	double	l, b, w, lwmin, lwmax;
 	register int	x, y;
-
+					/* check for precalculated histogram */
+	if (what2do&DO_PREHIST)
+		return;
 	lwmin = 1e10;			/* find extrema */
 	lwmax = 0.;
 	for (y = 0; y < fvyr; y++)
@@ -107,12 +175,17 @@ comphist()			/* create foveal sampling histogram */
 			if (l < lwmin) lwmin = l;
 			if (l > lwmax) lwmax = l;
 		}
-	lwmin -= FTINY;
-	lwmax += FTINY;
-	if (lwmin < LMIN) lwmin = LMIN;
-	if (lwmax > LMAX) lwmax = LMAX;
-	bwmin = Bl(lwmin);
+	lwmax *= 1.01;
+	if (lwmax > LMAX)
+		lwmax = LMAX;
 	bwmax = Bl(lwmax);
+	if (lwmin < LMIN) {
+		lwmin = LMIN;
+		bwmin = Bl(LMIN);
+	} else {			/* duplicate bottom bin */
+		bwmin = bwmax - (bwmax-Bl(lwmin))*HISTRES/(HISTRES-1);
+		lwmin = Lb(bwmin);
+	}
 					/* (re)compute histogram */
 	bwavg = 0.;
 	histot = 0.;
@@ -123,11 +196,11 @@ comphist()			/* create foveal sampling histogram */
 		for (y = 0; y < fvyr; y++)
 			for (x = 0; x < fvxr; x++) {
 				l = plum(fovscan(y)[x]);
-				if (l < lwmin) continue;
-				if (l > lwmax) continue;
+				if (l < lwmin+FTINY) continue;
+				if (l >= lwmax-FTINY) continue;
 				b = Bl(l);
-				bwavg += b;
 				w = what2do&DO_CWEIGHT ? centprob(x,y) : 1.;
+				bwavg += w*b;
 				bwhist[bwhi(b)] += w;
 				histot += w;
 			}
@@ -139,15 +212,19 @@ comphist()			/* create foveal sampling histogram */
 			w = 1.;
 		for (x = 0; x < nfixations; x++) {
 			l = plum(fovscan(fixlst[x][1])[fixlst[x][0]]);
-			if (l < lwmin) continue;
-			if (l > lwmax) continue;
+			if (l < lwmin+FTINY) continue;
+			if (l >= lwmax-FTINY) continue;
 			b = Bl(l);
-			bwavg += b;
+			bwavg += w*b;
 			bwhist[bwhi(b)] += w;
 			histot += w;
 		}
 	}
 	bwavg /= histot;
+	if (lwmin > LMIN+FTINY) {	/* add false samples at bottom */
+		bwhist[1] *= 0.5;
+		bwhist[0] += bwhist[1];
+	}
 }
 
 
@@ -219,7 +296,7 @@ double	La;
 
 
 double
-clampf(Lw)		/* derivative clamping function */
+clampf(Lw)			/* histogram clamping function */
 double	Lw;
 {
 	double	bLw, ratio;
@@ -229,25 +306,61 @@ double	Lw;
 	return(ratio/(Lb1(bLw)*(Bldmax-Bldmin)*Bl1(Lw)));
 }
 
+double
+crfactor(Lw)			/* contrast reduction factor */
+double	Lw;
+{
+	int	i = HISTRES*(Bl(Lw) - bwmin)/(bwmax - bwmin);
+	double	bLw, ratio, Tdb;
 
+	if (i <= 0)
+		return(1.0);
+	if (i >= HISTRES)
+		return(1.0);
+	bLw = BLw(Lw);
+	ratio = what2do&DO_HSENS ? htcontrs(Lb(bLw))/htcontrs(Lw) : Lb(bLw)/Lw;
+	Tdb = mhistot * (bwmax - bwmin) / HISTRES;
+	return(modhist[i]*Lb1(bLw)*(Bldmax-Bldmin)*Bl1(Lw)/(Tdb*ratio));
+}
+
+
+#if ADJ_VEIL
+mkcrfimage()			/* compute contrast reduction factor image */
+{
+	int	i;
+	float	*crfptr;
+	COLOR	*fovptr;
+
+	if (crfimg == NULL)
+		crfimg = (float *)malloc(fvxr*fvyr*sizeof(float));
+	if (crfimg == NULL)
+		syserror("malloc");
+	crfptr = crfimg;
+	fovptr = fovimg;
+	for (i = fvxr*fvyr; i--; crfptr++, fovptr++)
+		crfptr[0] = crfactor(plum(fovptr[0]));
+}
+#endif
+
+
 int
 mkbrmap()			/* make dynamic range map */
 {
-	double	T, b, s;
+	double	Tdb, b, s;
 	double	ceiling, trimmings;
 	register int	i;
 					/* copy initial histogram */
-	bcopy((char *)bwhist, (char *)modhist, sizeof(modhist));
-	s = (bwmax - bwmin)/HISTRES;
+	bcopy((void *)bwhist, (void *)modhist, sizeof(modhist));
+	s = (bwmax - bwmin)/HISTRES;	/* s is delta b */
 					/* loop until satisfactory */
 	do {
 		mkcumf();			/* sync brightness mapping */
 		if (mhistot <= histot*CVRATIO)
 			return(-1);		/* no compression needed! */
-		T = mhistot * (bwmax - bwmin) / HISTRES;
+		Tdb = mhistot * s;
 		trimmings = 0.;			/* clip to envelope */
 		for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
-			ceiling = T*clampf(Lb(b));
+			ceiling = Tdb*clampf(Lb(b));
 			if (modhist[i] > ceiling) {
 				trimmings += modhist[i] - ceiling;
 				modhist[i] = ceiling;
@@ -255,6 +368,10 @@ mkbrmap()			/* make dynamic range map */
 		}
 	} while (trimmings > histot*CVRATIO);
 
+#if ADJ_VEIL
+	mkcrfimage();			/* contrast reduction image */
+#endif
+
 	return(0);			/* we got it */
 }
 
@@ -295,18 +412,18 @@ COLOR	*scan;
 int	xres;
 {
 	double	mult, Lw, b;
-	register int	i;
+	register int	x;
 
-	for (i = 0; i < xres; i++) {
-		Lw = plum(scan[i]);
+	for (x = 0; x < xres; x++) {
+		Lw = plum(scan[x]);
 		if (Lw < LMIN) {
-			setcolor(scan[i], 0., 0., 0.);
+			setcolor(scan[x], 0., 0., 0.);
 			continue;
 		}
-		b = BLw(Lw);
-		mult = (Lb(b) - ldmin)/(ldmax - ldmin) / (Lw*inpexp);
+		b = BLw(Lw);		/* apply brightness mapping */
+		mult = (Lb(b) - ldmin)/(ldmax - ldmin)/(Lw*inpexp);
 		if (lumf == rgblum) mult *= WHTEFFICACY;
-		scalecolor(scan[i], mult);
+		scalecolor(scan[x], mult);
 	}
 }