--- ray/src/px/pcond3.c	1996/10/03 16:52:50	3.1
+++ ray/src/px/pcond3.c	1997/01/10 16:48:21	3.6
@@ -11,26 +11,27 @@ 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 exp10(x)	exp(2.302585093*(x))
 
-int	modhist[HISTRES];		/* modified histogram */
+float	modhist[HISTRES];		/* modified histogram */
+double	mhistot;			/* modified histogram total */
 float	cumf[HISTRES+1];		/* cumulative distribution function */
 
 
 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 +61,7 @@ double	Lw;
 		return(Bldmin);
 	if (b >= bwmax-FTINY)
 		return(Bldmax);
-	return(Bldmin + cf(Bl(Lw))*(Bldmax-Bldmin));
+	return(Bldmin + cf(b)*(Bldmax-Bldmin));
 }
 
 
@@ -99,66 +100,31 @@ double	Lw;
 
 
 int
-shiftdir(bw)		/* compute shift direction for histogram */
-double	bw;
-{
-	if (what2do&DO_HSENS && cf(bw) - (bw - bwmin)/(Bldmax - bwmin))
-		return(1);
-	return(-1);
-}
-
-
-int
 mkbrmap()			/* make dynamic range map */
 {
-	int	hdiffs[HISTRES], above, below;
 	double	T, b, s;
-	int	maxd, maxi, sd;
+	double	ceiling, trimmings;
 	register int	i;
 					/* copy initial histogram */
-	for (i = 0; i < HISTRES; i++)
-		modhist[i] = bwhist[i];
-	T = histot * (bwmax - bwmin) / HISTRES;
+	bcopy((char *)bwhist, (char *)modhist, sizeof(modhist));
 	s = (bwmax - bwmin)/HISTRES;
 					/* 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! */
+		T = mhistot * (bwmax - bwmin) / HISTRES;
+		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 = T*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);
+
+	return(0);			/* we got it */
 }
 
 
@@ -180,13 +146,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);
 		}
 	}
@@ -214,37 +180,21 @@ int	xres;
 }
 
 
-#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;
 
-	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)
+			fprintf(fp, "%e %e\n", wlum, sf*wlum);
+		else
+			fprintf(fp, "%e %e\n", wlum, Lb(BLw(wlum)));
+	}
 }
-#endif