--- ray/src/px/pcond3.c	1996/10/03 16:52:50	3.1
+++ ray/src/px/pcond3.c	2003/02/22 02:07:27	3.11
@@ -1,9 +1,6 @@
-/* Copyright (c) 1996 Regents of the University of California */
-
 #ifndef lint
-static char SCCSid[] = "$SunId$ LBL";
+static const char	RCSid[] = "$Id: pcond3.c,v 3.11 2003/02/22 02:07:27 greg Exp $";
 #endif
-
 /*
  * Routines for computing and applying brightness mapping.
  */
@@ -11,26 +8,238 @@ static char SCCSid[] = "$SunId$ LBL";
 #include "pcond.h"
 
 
-#define CVRATIO		0.025		/* fraction of pixels allowed > env. */
+#define CVRATIO		0.025		/* fraction of samples allowed > env. */
 
-#define BotMesopic	5.62e-3		/* top of scotopic range */
-#define TopMesopic	5.62		/* bottom of photopic range */
+#define LN_10		2.30258509299404568402
+#define exp10(x)	exp(LN_10*(x))
 
-#define exp10(x)	exp(2.302585093*(x))
+double	modhist[HISTRES];		/* modified histogram */
+double	mhistot;			/* modified histogram total */
+double	cumf[HISTRES+1];		/* cumulative distribution function */
 
-int	modhist[HISTRES];		/* modified histogram */
-float	cumf[HISTRES+1];		/* cumulative distribution function */
 
+getfixations(fp)		/* 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.rt=inpres.rt)&YMAJOR) {
+		fvres.xr = fvxr;
+		fvres.yr = fvyr;
+	} else {
+		fvres.xr = fvyr;
+		fvres.yr = fvxr;
+	}
+				/* read each picture position */
+	while (fscanf(fp, "%d %d", &pos[0], &pos[1]) == 2) {
+				/* convert to closest index in foveal image */
+		loc2pix(pos, &fvres,
+				(pos[0]+.5)/inpres.xr, (pos[1]+.5)/inpres.yr);
+				/* include nine neighborhood samples */
+		for (px = pos[0]-1; px <= pos[0]+1; px++) {
+			if (px < 0 || px >= fvxr)
+				continue;
+			for (py = pos[1]-1; py <= pos[1]+1; py++) {
+				if (py < 0 || py >= fvyr)
+					continue;
+				for (i = nfixations; i-- > 0; )
+					if (fixlst[i][0] == px &&
+							fixlst[i][1] == py)
+						break;
+				if (i >= 0)
+					continue;	/* already there */
+				if (nfixations % FIXHUNK == 0) {
+					if (nfixations)
+						fixlst = (short (*)[2])
+							realloc((char *)fixlst,
+							(nfixations+FIXHUNK)*
+							2*sizeof(short));
+					else
+						fixlst = (short (*)[2])malloc(
+							FIXHUNK*2*sizeof(short)
+							);
+					if (fixlst == NULL)
+						syserror("malloc");
+				}
+				fixlst[nfixations][0] = px;
+				fixlst[nfixations][1] = py;
+				nfixations++;
+			}
+		}
+	}
+	if (!feof(fp)) {
+		fprintf(stderr, "%s: format error reading fixation data\n",
+				progname);
+		exit(1);
+	}
+#undef	FIXHUNK
+}
 
+
+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;
+{
+	double	xr, yr, p;
+				/* paraboloid, 0 at 90 degrees from center */
+	xr = (x - .5*(fvxr-1))/90.;	/* 180 degree fisheye has fv?r == 90 */
+	yr = (y - .5*(fvyr-1))/90.;
+	p = 1. - xr*xr - yr*yr;
+	return(p < 0. ? 0. : p);
+}
+
+
+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++)
+		for (x = 0; x < fvxr; x++) {
+			l = plum(fovscan(y)[x]);
+			if (l < lwmin) lwmin = l;
+			if (l > lwmax) lwmax = l;
+		}
+	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.;
+	for (x = 0; x < HISTRES; x++)
+		bwhist[x] = 0.;
+					/* global average */
+	if (!(what2do&DO_FIXHIST) || fixfrac < 1.-FTINY)
+		for (y = 0; y < fvyr; y++)
+			for (x = 0; x < fvxr; x++) {
+				l = plum(fovscan(y)[x]);
+				if (l < lwmin+FTINY) continue;
+				if (l >= lwmax-FTINY) continue;
+				b = Bl(l);
+				w = what2do&DO_CWEIGHT ? centprob(x,y) : 1.;
+				bwavg += w*b;
+				bwhist[bwhi(b)] += w;
+				histot += w;
+			}
+					/* average fixation points */
+	if (what2do&DO_FIXHIST && nfixations > 0) {
+		if (histot > FTINY)
+			w = fixfrac/(1.-fixfrac)*histot/nfixations;
+		else
+			w = 1.;
+		for (x = 0; x < nfixations; x++) {
+			l = plum(fovscan(fixlst[x][1])[fixlst[x][0]]);
+			if (l < lwmin+FTINY) continue;
+			if (l >= lwmax-FTINY) continue;
+			b = Bl(l);
+			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 */
 {
 	register int	i;
-	register long	sum;
+	register double	sum;
 
-	cumf[0] = 0.;
-	sum = modhist[0];
-	for (i = 1; i < HISTRES; i++) {
-		cumf[i] = (double)sum/histot;
+	mhistot = 0.;		/* compute modified total */
+	for (i = 0; i < HISTRES; i++)
+		mhistot += modhist[i];
+
+	sum = 0.;		/* compute cumulative function */
+	for (i = 0; i < HISTRES; i++) {
+		cumf[i] = sum/mhistot;
 		sum += modhist[i];
 	}
 	cumf[HISTRES] = 1.;
@@ -60,7 +269,7 @@ double	Lw;
 		return(Bldmin);
 	if (b >= bwmax-FTINY)
 		return(Bldmax);
-	return(Bldmin + cf(Bl(Lw))*(Bldmax-Bldmin));
+	return(Bldmin + cf(b)*(Bldmax-Bldmin));
 }
 
 
@@ -87,7 +296,7 @@ double	La;
 
 
 double
-clampf(Lw)		/* derivative clamping function */
+clampf(Lw)			/* histogram clamping function */
 double	Lw;
 {
 	double	bLw, ratio;
@@ -97,68 +306,73 @@ 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;
 
-int
-shiftdir(bw)		/* compute shift direction for histogram */
-double	bw;
+	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 */
 {
-	if (what2do&DO_HSENS && cf(bw) - (bw - bwmin)/(Bldmax - bwmin))
-		return(1);
-	return(-1);
+	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 */
 {
-	int	hdiffs[HISTRES], above, below;
-	double	T, b, s;
-	int	maxd, maxi, sd;
+	double	Tdb, b, s;
+	double	ceiling, trimmings;
 	register int	i;
 					/* copy initial histogram */
-	for (i = 0; i < HISTRES; i++)
-		modhist[i] = bwhist[i];
-	T = histot * (bwmax - bwmin) / HISTRES;
-	s = (bwmax - bwmin)/HISTRES;
+	bcopy((char *)bwhist, (char *)modhist, sizeof(modhist));
+	s = (bwmax - bwmin)/HISTRES;	/* s is delta b */
 					/* loop until satisfactory */
-	for ( ; ; ) {
-		mkcumf();		/* sync brightness mapping */
-		above = below = 0;	/* compute visibility overflow */
+	do {
+		mkcumf();			/* sync brightness mapping */
+		if (mhistot <= histot*CVRATIO)
+			return(-1);		/* no compression needed! */
+		Tdb = mhistot * s;
+		trimmings = 0.;			/* clip to envelope */
 		for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
-			hdiffs[i] = modhist[i] - (int)(T*clampf(Lb(b)) + .5);
-			if (hdiffs[i] > 0) above += hdiffs[i];
-			else below -= hdiffs[i];
+			ceiling = Tdb*clampf(Lb(b));
+			if (modhist[i] > ceiling) {
+				trimmings += modhist[i] - ceiling;
+				modhist[i] = ceiling;
+			}
 		}
-		if (above <= histot*CVRATIO)
-			break;		/* close enough */
-		if (above-below >= 0)
-			return(-1);	/* Houston, we have a problem.... */
-		/* original looped here as well (BEGIN_L2) */
-		maxd = 0;		/* find largest overvis */
-		for (i = 0; i < HISTRES; i++)
-			if (hdiffs[i] > maxd)
-				maxd = hdiffs[maxi=i];
-		/* broke loop here when (maxd == 0) (BREAK_L2) */
-		for (sd = shiftdir((maxi+.5)/HISTRES*(bwmax-bwmin)+bwmin);
-				hdiffs[maxi] == maxd; sd = -sd)
-			for (i = maxi+sd; i >= 0 & i < HISTRES; i += sd)
-				if (hdiffs[i] < 0) {
-					if (hdiffs[i] <= -maxd) {
-						modhist[i] += maxd;
-						modhist[maxi] -= maxd;
-						hdiffs[i] += maxd;
-						hdiffs[maxi] = 0;
-					} else {
-						modhist[maxi] += hdiffs[i];
-						modhist[i] -= hdiffs[i];
-						hdiffs[maxi] += hdiffs[i];
-						hdiffs[i] = 0;
-					}
-					break;
-				}
-		/* (END_L2) */
-	}
-	return(0);
+	} while (trimmings > histot*CVRATIO);
+
+#if ADJ_VEIL
+	mkcrfimage();			/* contrast reduction image */
+#endif
+
+	return(0);			/* we got it */
 }
 
 
@@ -180,13 +394,13 @@ int	xres;
 			incolor = (Lw - BotMesopic) /
 					(TopMesopic - BotMesopic);
 		if (incolor < 1.-FTINY) {
+			b = (1.-incolor)*slum(scan[i])*inpexp/SWNORM;
+			if (lumf == rgblum) b /= WHTEFFICACY;
+			setcolor(ctmp, b, b, b);
 			if (incolor <= FTINY)
 				setcolor(scan[i], 0., 0., 0.);
 			else
 				scalecolor(scan[i], incolor);
-			b = (1.-incolor)*slum(scan[i])*inpexp/SWNORM;
-			if (lumf == rgblum) b /= WHTEFFICACY;
-			setcolor(ctmp, b, b, b);
 			addcolor(scan[i], ctmp);
 		}
 	}
@@ -198,53 +412,40 @@ 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);
 	}
 }
 
 
-#ifdef DEBUG
-doplots()			/* generate debugging plots */
+putmapping(fp)			/* put out mapping function */
+FILE	*fp;
 {
-	double	T, b, s;
-	FILE	*fp;
-	char	fname[128];
+	double	b, s;
 	register int	i;
+	double	wlum, sf, dlum;
 
-	T = histot * (bwmax - bwmin) / HISTRES;
+	sf = scalef*inpexp;
+	if (lumf == cielum) sf *= WHTEFFICACY;
 	s = (bwmax - bwmin)/HISTRES;
-
-	sprintf(fname, "%s_hist.plt", infn);
-	if ((fp = fopen(fname, "w")) == NULL)
-		syserror(fname);
-	fputs("include=curve.plt\n", fp);
-	fputs("title=\"Brightness Frequency Distribution\"\n", fp);
-	fprintf(fp, "subtitle=%s\n", infn);
-	fputs("ymin=0\n", fp);
-	fputs("xlabel=\"Perceptual Brightness B(Lw)\"\n", fp);
-	fputs("ylabel=\"Frequency Count\"\n", fp);
-	fputs("Alabel=\"Histogram\"\n", fp);
-	fputs("Alintype=0\n", fp);
-	fputs("Blabel=\"Envelope\"\n", fp);
-	fputs("Bsymsize=0\n", fp);
-	fputs("Adata=\n", fp);
-	for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s)
-		fprintf(fp, "\t%f %d\n", b, modhist[i]);
-	fputs(";\nBdata=\n", fp);
-	for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s)
-		fprintf(fp, "\t%f %f\n", b, T*clampf(Lb(b)));
-	fputs(";\n", fp);
-	fclose(fp);
+	for (i = 0, b = bwmin + .5*s; i < HISTRES; i++, b += s) {
+		wlum = Lb(b);
+		if (what2do&DO_LINEAR) {
+			dlum = sf*wlum;
+			if (dlum > ldmax) dlum = ldmax;
+			else if (dlum < ldmin) dlum = ldmin;
+			fprintf(fp, "%e %e\n", wlum, dlum);
+		} else
+			fprintf(fp, "%e %e\n", wlum, Lb(BLw(wlum)));
+	}
 }
-#endif