--- ray/src/px/pcond3.c	1998/03/12 15:47:34	3.10
+++ ray/src/px/pcond3.c	2021/04/13 02:56:42	3.19
@@ -1,14 +1,14 @@
-/* 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.19 2021/04/13 02:56:42 greg Exp $";
 #endif
-
 /*
  * Routines for computing and applying brightness mapping.
  */
 
+#include <string.h>
+
 #include "pcond.h"
+#include "random.h"
 
 
 #define CVRATIO		0.025		/* fraction of samples allowed > env. */
@@ -16,20 +16,32 @@ static char SCCSid[] = "$SunId$ LBL";
 #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);
+static float *getlumsamp(int n);
+#if ADJ_VEIL
+static void mkcrfimage(void);
+#endif
 
-getfixations(fp)		/* load fixation history list */
-FILE	*fp;
+
+
+void
+getfixations(		/* load fixation history list */
+FILE	*fp
+)
 {
 #define	FIXHUNK		128
 	RESOLU	fvres;
 	int	pos[2];
-	register int	px, py, i;
+	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 {
@@ -57,7 +69,7 @@ FILE	*fp;
 				if (nfixations % FIXHUNK == 0) {
 					if (nfixations)
 						fixlst = (short (*)[2])
-							realloc((char *)fixlst,
+							realloc((void *)fixlst,
 							(nfixations+FIXHUNK)*
 							2*sizeof(short));
 					else
@@ -82,25 +94,32 @@ FILE	*fp;
 }
 
 
-gethisto(fp)			/* load precomputed luminance histogram */
-FILE	*fp;
+void
+gethisto(			/* load precomputed luminance histogram */
+	FILE	*fp
+)
 {
-	float	histo[MAXPREHIST];
+	double	histo[MAXPREHIST];
 	double	histart, histep;
-	double	l, b, lastb, w;
+	double	b, lastb, w;
 	int	n;
-	register int	i;
+	int	i;
 					/* load data */
 	for (i = 0; i < MAXPREHIST &&
-			fscanf(fp, "%lf %f", &b, &histo[i]) == 2; i++) {
-		if (i > 1 && fabs(b - lastb - histep) > 1e-3) {
+			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)
-			histep = b - (histart = lastb);
+			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)) {
@@ -115,10 +134,10 @@ FILE	*fp;
 					/* find extrema */
 	for (i = 0; i < n && histo[i] <= FTINY; i++)
 		;
-	bwmin = histart + i*histep;
+	bwmin = histart + (i-.001)*histep;
 	for (i = n; i-- && histo[i] <= FTINY; )
 		;
-	bwmax = histart + i*histep;
+	bwmax = histart + (i+1.001)*histep;
 	if (bwmax > Bl(LMAX))
 		bwmax = Bl(LMAX);
 	if (bwmin < Bl(LMIN))
@@ -130,8 +149,8 @@ FILE	*fp;
 	for (i = 0; i < HISTRES; i++)
 		bwhist[i] = 0.;
 	for (i = 0, b = histart; i < n; i++, b += histep) {
-		if (b < bwmin) continue;
-		if (b > bwmax) break;
+		if (b < bwmin+FTINY) continue;
+		if (b >= bwmax-FTINY) break;
 		w = histo[i];
 		bwavg += w*b;
 		bwhist[bwhi(b)] += w;
@@ -145,9 +164,11 @@ FILE	*fp;
 }
 
 
-double
-centprob(x, y)			/* center-weighting probability function */
-int	x, y;
+static double
+centprob(			/* center-weighting probability function */
+	int	x,
+	int	y
+)
 {
 	double	xr, yr, p;
 				/* paraboloid, 0 at 90 degrees from center */
@@ -158,10 +179,55 @@ int	x, y;
 }
 
 
-comphist()			/* create foveal sampling histogram */
+static float *
+getlumsamp(int n)		/* get set of random point sample luminances */
 {
+	float	*ls = (float *)malloc(n*sizeof(float));
+	COLR	*cscan = (COLR *)malloc(scanlen(&inpres)*sizeof(COLR));
+	long	startpos = ftell(infp);
+	long	npleft = (long)inpres.xr*inpres.yr;
+	int	x;
+	
+	if ((ls == NULL) | (cscan == NULL))
+		syserror("malloc");
+	x = 0;				/* read/convert samples */
+	while (n > 0) {
+		COLOR	col;
+		int	sv = 0;
+		double	rval, cumprob = 0;
+
+		if (x <= 0 && freadcolrs(cscan, x=scanlen(&inpres), infp) < 0) {
+			fprintf(stderr, "%s: %s: scanline read error\n",
+					progname, infn);
+			exit(1);
+		}
+		rval = frandom();	/* random distance to next sample */
+		while ((cumprob += (1.-cumprob)*n/(npleft-sv)) < rval)
+			sv++;
+		if (x < ++sv) {		/* out of pixels in this scanline */
+			npleft -= x;
+			x = 0;
+			continue;
+		}
+		x -= sv;
+		colr_color(col, cscan[x]);
+		ls[--n] = plum(col);
+		npleft -= sv;
+	}
+	free(cscan);			/* clean up and reset file pointer */
+	if (fseek(infp, startpos, SEEK_SET) < 0)
+		syserror("fseek");
+	return(ls);
+}
+
+
+void
+comphist(void)			/* create foveal sampling histogram */
+{
 	double	l, b, w, lwmin, lwmax;
-	register int	x, y;
+	float	*lumsamp;
+	int	nlumsamp;
+	int	x, y;
 					/* check for precalculated histogram */
 	if (what2do&DO_PREHIST)
 		return;
@@ -173,6 +239,16 @@ comphist()			/* create foveal sampling histogram */
 			if (l < lwmin) lwmin = l;
 			if (l > lwmax) lwmax = l;
 		}
+					/* sample luminance pixels */
+	nlumsamp = fvxr*fvyr*16;
+	if (nlumsamp > inpres.xr*inpres.yr)
+		nlumsamp = inpres.xr*inpres.yr;
+	lumsamp = getlumsamp(nlumsamp);
+	for (x = nlumsamp; x--; ) {
+		l = lumsamp[x];
+		if (l < lwmin) lwmin = l;
+		if (l > lwmax) lwmax = l;
+	}
 	lwmax *= 1.01;
 	if (lwmax > LMAX)
 		lwmax = LMAX;
@@ -187,21 +263,32 @@ comphist()			/* create foveal sampling histogram */
 					/* (re)compute histogram */
 	bwavg = 0.;
 	histot = 0.;
-	for (x = 0; x < HISTRES; x++)
-		bwhist[x] = 0.;
+	memset(bwhist, 0, sizeof(bwhist));
 					/* global average */
-	if (!(what2do&DO_FIXHIST) || fixfrac < 1.-FTINY)
+	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) continue;
-				if (l > lwmax) continue;
+				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;
 			}
+					/* weight for point luminances */
+		w = 1. * histot / nlumsamp;
+		for (x = nlumsamp; x--; ) {
+			l = lumsamp[x];
+			if (l < lwmin+FTINY) continue;
+			if (l >= lwmax-FTINY) continue;
+			b = Bl(l);
+			bwavg += w*b;
+			bwhist[bwhi(b)] += w;
+			histot += w;
+		}
+	}
 					/* average fixation points */
 	if (what2do&DO_FIXHIST && nfixations > 0) {
 		if (histot > FTINY)
@@ -210,8 +297,8 @@ 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 += w*b;
 			bwhist[bwhi(b)] += w;
@@ -223,13 +310,15 @@ comphist()			/* create foveal sampling histogram */
 		bwhist[1] *= 0.5;
 		bwhist[0] += bwhist[1];
 	}
+	free(lumsamp);
 }
 
 
-mkcumf()			/* make cumulative distribution function */
+static void
+mkcumf(void)			/* make cumulative distribution function */
 {
-	register int	i;
-	register double	sum;
+	int	i;
+	double	sum;
 
 	mhistot = 0.;		/* compute modified total */
 	for (i = 0; i < HISTRES; i++)
@@ -244,12 +333,13 @@ 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;
+	int	i;
 
 	i = x = HISTRES*(b - bwmin)/(bwmax - bwmin);
 	x -= (double)i;
@@ -257,9 +347,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;
 
@@ -272,8 +363,9 @@ double	Lw;
 
 
 double
-htcontrs(La)		/* human threshold contrast sensitivity, dL(La) */
-double	La;
+htcontrs(		/* human threshold contrast sensitivity, dL(La) */
+	double	La
+)
 {
 	double	l10La, l10dL;
 				/* formula taken from Ferwerda et al. [SG96] */
@@ -294,8 +386,9 @@ double	La;
 
 
 double
-clampf(Lw)		/* derivative clamping function */
-double	Lw;
+clampf(			/* histogram clamping function */
+	double	Lw
+)
 {
 	double	bLw, ratio;
 
@@ -304,43 +397,87 @@ double	Lw;
 	return(ratio/(Lb1(bLw)*(Bldmax-Bldmin)*Bl1(Lw)));
 }
 
+double
+crfactor(			/* 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
+static void
+mkcrfimage(void)			/* 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 */
+mkbrmap(void)			/* make dynamic range map */
 {
-	double	T, b, s;
+	double	Tdb, b, s;
 	double	ceiling, trimmings;
-	register int	i;
+	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;
 			}
 		}
-	} while (trimmings > histot*CVRATIO);
+	} while (trimmings > mhistot*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;
+void
+scotscan(		/* apply scotopic color sensitivity loss */
+	COLOR	*scan,
+	int	xres
+)
 {
 	COLOR	ctmp;
 	double	incolor, b, Lw;
-	register int	i;
+	int	i;
 
 	for (i = 0; i < xres; i++) {
 		Lw = plum(scan[i]);
@@ -365,32 +502,36 @@ int	xres;
 }
 
 
-mapscan(scan, xres)		/* apply tone mapping operator to scanline */
-COLOR	*scan;
-int	xres;
+void
+mapscan(		/* apply tone mapping operator to scanline */
+	COLOR	*scan,
+	int	xres
+)
 {
 	double	mult, Lw, b;
-	register int	i;
+	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;
+void
+putmapping(			/* put out mapping function */
+	FILE	*fp
+)
 {
 	double	b, s;
-	register int	i;
+	int	i;
 	double	wlum, sf, dlum;
 
 	sf = scalef*inpexp;