--- ray/src/px/pcond3.c	1997/01/29 13:22:13	3.8
+++ ray/src/px/pcond3.c	2004/03/28 20:33:14	3.15
@@ -1,34 +1,45 @@
-/* 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.15 2004/03/28 20:33:14 schorsch Exp $";
 #endif
-
 /*
  * Routines for computing and applying brightness mapping.
  */
 
+#include <string.h>
+
 #include "pcond.h"
 
 
 #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 */
 
+static double centprob(int	x, int	y);
+static void mkcumf(void);
+static double cf(double	b);
+static double BLw(double	Lw);
+#if ADJ_VEIL
+static void mkcrfimage(void);
+#endif
 
-getfixations(fp)		/* load fixation history list */
-FILE	*fp;
+
+
+extern void
+getfixations(		/* load fixation history list */
+FILE	*fp
+)
 {
 #define	FIXHUNK		128
 	RESOLU	fvres;
 	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 +67,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,10 +92,82 @@ FILE	*fp;
 }
 
 
-double
-centprob(x, y)			/* center-weighting probability function */
-int	x, y;
+extern void
+gethisto(			/* load precomputed luminance histogram */
+	FILE	*fp
+)
 {
+	double	histo[MAXPREHIST];
+	double	histart, histep;
+	double	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];
+	}
+}
+
+
+static double
+centprob(			/* center-weighting probability function */
+	int	x,
+	int	y
+)
+{
 	double	xr, yr, p;
 				/* paraboloid, 0 at 90 degrees from center */
 	xr = (x - .5*(fvxr-1))/90.;	/* 180 degree fisheye has fv?r == 90 */
@@ -94,11 +177,14 @@ int	x, y;
 }
 
 
-comphist()			/* create foveal sampling histogram */
+extern void
+comphist(void)			/* 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 +193,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 +214,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,19 +230,24 @@ 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];
+	}
 }
 
 
-mkcumf()			/* make cumulative distribution function */
+static void
+mkcumf(void)			/* make cumulative distribution function */
 {
 	register int	i;
 	register double	sum;
@@ -169,9 +265,10 @@ mkcumf()			/* make cumulative distribution function */
 }
 
 
-double
-cf(b)				/* return cumulative function at b */
-double	b;
+static double
+cf(				/* return cumulative function at b */
+	double	b
+)
 {
 	double	x;
 	register int	i;
@@ -182,9 +279,10 @@ double	b;
 }
 
 
-double
-BLw(Lw)				/* map world luminance to display brightness */
-double	Lw;
+static double
+BLw(				/* map world luminance to display brightness */
+	double	Lw
+)
 {
 	double	b;
 
@@ -196,9 +294,10 @@ double	Lw;
 }
 
 
-double
-htcontrs(La)		/* human threshold contrast sensitivity, dL(La) */
-double	La;
+extern double
+htcontrs(		/* human threshold contrast sensitivity, dL(La) */
+	double	La
+)
 {
 	double	l10La, l10dL;
 				/* formula taken from Ferwerda et al. [SG96] */
@@ -218,9 +317,10 @@ double	La;
 }
 
 
-double
-clampf(Lw)		/* derivative clamping function */
-double	Lw;
+extern double
+clampf(			/* histogram clamping function */
+	double	Lw
+)
 {
 	double	bLw, ratio;
 
@@ -229,25 +329,63 @@ double	Lw;
 	return(ratio/(Lb1(bLw)*(Bldmax-Bldmin)*Bl1(Lw)));
 }
 
+extern double
+crfactor(			/* contrast reduction factor */
+	double	Lw
+)
+{
+	int	i = HISTRES*(Bl(Lw) - bwmin)/(bwmax - bwmin);
+	double	bLw, ratio, Tdb;
 
-int
-mkbrmap()			/* make dynamic range map */
+	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
+static void
+mkcrfimage(void)			/* compute contrast reduction factor image */
 {
-	double	T, b, s;
+	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
+
+
+extern int
+mkbrmap(void)			/* make dynamic range map */
+{
+	double	Tdb, b, s;
 	double	ceiling, trimmings;
 	register int	i;
 					/* copy initial histogram */
-	bcopy((char *)bwhist, (char *)modhist, sizeof(modhist));
-	s = (bwmax - bwmin)/HISTRES;
+	memcpy((void *)modhist, (void *)bwhist, 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,13 +393,19 @@ mkbrmap()			/* make dynamic range map */
 		}
 	} while (trimmings > histot*CVRATIO);
 
+#if ADJ_VEIL
+	mkcrfimage();			/* contrast reduction image */
+#endif
+
 	return(0);			/* we got it */
 }
 
 
-scotscan(scan, xres)		/* apply scotopic color sensitivity loss */
-COLOR	*scan;
-int	xres;
+extern void
+scotscan(		/* apply scotopic color sensitivity loss */
+	COLOR	*scan,
+	int	xres
+)
 {
 	COLOR	ctmp;
 	double	incolor, b, Lw;
@@ -290,29 +434,33 @@ int	xres;
 }
 
 
-mapscan(scan, xres)		/* apply tone mapping operator to scanline */
-COLOR	*scan;
-int	xres;
+extern void
+mapscan(		/* apply tone mapping operator to scanline */
+	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);
 	}
 }
 
 
-putmapping(fp)			/* put out mapping function */
-FILE	*fp;
+extern void
+putmapping(			/* put out mapping function */
+	FILE	*fp
+)
 {
 	double	b, s;
 	register int	i;