--- ray/src/px/pcond.c	1996/10/10 17:09:24	3.4
+++ ray/src/px/pcond.c	2003/02/22 02:07:27	3.14
@@ -1,40 +1,41 @@
-/* Copyright (c) 1996 Regents of the University of California */
-
 #ifndef lint
-static char SCCSid[] = "$SunId$ LBL";
+static const char	RCSid[] = "$Id: pcond.c,v 3.14 2003/02/22 02:07:27 greg Exp $";
 #endif
-
 /*
  * Condition Radiance picture for display/output
+ *  Added white-balance adjustment 10/01 (GW).
  */
 
 #include "pcond.h"
 
-#include "random.h"
 
-
 #define	LDMAX		100		/* default max. display luminance */
-#define LDMINF		0.01		/* default min. display lum. factor */
+#define LDDYN		32		/* default dynamic range */
 
 int	what2do = 0;			/* desired adjustments */
 
 double	ldmax = LDMAX;			/* maximum output luminance */
-double	ldmin = 0.;			/* minimum output luminance */
-double	Bldmin, Bldmax;			/* Bl(ldmin) and Bl(ldmax) */
+double	lddyn = LDDYN;			/* display dynamic range */
+double	Bldmin, Bldmax;			/* Bl(ldmax/lddyn) and Bl(ldmax) */
 
 char	*progname;			/* global argv[0] */
 
 char	*infn;				/* input file name */
 FILE	*infp;				/* input stream */
+FILE	*mapfp = NULL;			/* tone-mapping function stream */
 VIEW	ourview = STDVIEW;		/* picture view */
 int	gotview = 0;			/* picture has view */
 double	pixaspect = 1.0;		/* pixel aspect ratio */
+double	fixfrac = 0.;			/* histogram share due to fixations */
 RESOLU	inpres;				/* input picture resolution */
 
 COLOR	*fovimg;			/* foveal (1 degree) averaged image */
-short	fvxr, fvyr;			/* foveal image resolution */
-int	bwhist[HISTRES];		/* luminance histogram */
-long	histot;				/* total count of histogram */
+int	fvxr, fvyr;			/* foveal image resolution */
+float	*crfimg;			/* contrast reduction factors */
+short	(*fixlst)[2];			/* fixation history list */
+int	nfixations;			/* number of fixation points */
+double	bwhist[HISTRES];		/* luminance histogram */
+double	histot;				/* total count of histogram */
 double	bwmin, bwmax;			/* histogram limits */
 double	bwavg;				/* mean brightness */
 
@@ -77,6 +78,15 @@ char	*argv[];
 		case 'w':
 			bool(DO_CWEIGHT);
 			break;
+		case 'i':
+			if (i+1 >= argc) goto userr;
+			fixfrac = atof(argv[++i]);
+			if (fixfrac > FTINY) what2do |= DO_FIXHIST;
+			else what2do &= ~DO_FIXHIST;
+			break;
+		case 'I':
+			bool(DO_PREHIST);
+			break;
 		case 'l':
 			bool(DO_LINEAR);
 			break;
@@ -103,39 +113,53 @@ char	*argv[];
 			if (i+1 >= argc) goto userr;
 			mbcalfile = argv[++i];
 			break;
-		case 't':
+		case 'm':
 			if (i+1 >= argc) goto userr;
+			cwarpfile = argv[++i];
+			break;
+		case 'u':
+			if (i+1 >= argc) goto userr;
 			ldmax = atof(argv[++i]);
 			if (ldmax <= FTINY)
 				goto userr;
 			break;
-		case 'b':
+		case 'd':
 			if (i+1 >= argc) goto userr;
-			ldmin = atof(argv[++i]);
+			lddyn = atof(argv[++i]);
 			break;
+		case 'x':
+			if (i+1 >= argc) goto userr;
+			if ((mapfp = fopen(argv[++i], "w")) == NULL) {
+				fprintf(stderr,
+					"%s: cannot open for writing\n",
+						argv[i]);
+				exit(1);
+			}
+			break;
 		default:
 			goto userr;
 		}
-	if (mbcalfile != NULL & outprims != stdprims) {
-		fprintf(stderr, "%s: only one of -p or -f option supported\n",
+	if ((what2do & (DO_FIXHIST|DO_PREHIST)) == (DO_FIXHIST|DO_PREHIST)) {
+		fprintf(stderr, "%s: only one of -i or -I option\n", progname);
+		exit(1);
+	}
+	if ((mbcalfile != NULL) + (cwarpfile != NULL) +
+			(outprims != stdprims) > 1) {
+		fprintf(stderr,
+			"%s: only one of -p, -m or -f option supported\n",
 				progname);
 		exit(1);
 	}
 	if (outprims == stdprims & inprims != stdprims)
 		outprims = inprims;
-	if (ldmin <= FTINY)
-		ldmin = ldmax*LDMINF;
-	else if (ldmin >= ldmax) {
-		fprintf(stderr, "%s: Ldmin (%f) >= Ldmax (%f)!\n", progname,
-				ldmin, ldmax);
-		exit(1);
-	}
-	Bldmin = Bl(ldmin);
+	Bldmin = Bl(ldmax/lddyn);
 	Bldmax = Bl(ldmax);
 	if (i >= argc || i+2 < argc)
 		goto userr;
+					/* open input file */
 	if ((infp = fopen(infn=argv[i], "r")) == NULL)
 		syserror(infn);
+					/* open output file */
 	if (i+2 == argc && freopen(argv[i+1], "w", stdout) == NULL)
 		syserror(argv[i+1]);
 #ifdef MSDOS
@@ -144,12 +168,20 @@ char	*argv[];
 #endif
 	getahead();			/* load input header */
 	printargs(argc, argv, stdout);	/* add to output header */
-	if (outprims != inprims)
+	if (mbcalfile == NULL & outprims != stdprims)
 		fputprims(outprims, stdout);
+	if ((what2do & (DO_PREHIST|DO_VEIL|DO_ACUITY)) != DO_PREHIST)
+		getfovimg();		/* get foveal sample image? */
+	if (what2do&DO_PREHIST)		/* get histogram? */
+		gethisto(stdin);
+	else if (what2do&DO_FIXHIST)	/* get fixation history? */
+		getfixations(stdin);
 	mapimage();			/* map the picture */
+	if (mapfp != NULL)		/* write out basic mapping */
+		putmapping(mapfp);
 	exit(0);
 userr:
-	fprintf(stderr, "Usage: %s [-{h|a|v|s|c|l|w}[+-]][-e ev][-p xr yr xg yg xb yb xw yw|-f mbf.cal][-t Ldmax][-b Ldmin] inpic [outpic]\n",
+	fprintf(stderr, "Usage: %s [-{h|a|v|s|c|l|w}[+-]][-I|-i ffrac][-e ev][-p xr yr xg yg xb yb xw yw|-f mbf.cal|-m rgb.cwp][-u Ldmax][-d Lddyn][-x mapfile] inpic [outpic]\n",
 			progname);
 	exit(1);
 #undef bool
@@ -175,22 +207,22 @@ char	*s;
 		if (!strcmp(fmt,COLRFMT)) lumf = rgblum;
 		else if (!strcmp(fmt,CIEFMT)) lumf = cielum;
 		else lumf = NULL;
-		return;			/* don't echo */
+		return(0);		/* don't echo */
 	}
 	if (isprims(s)) {		/* get input primaries */
 		primsval(inprimS, s);
 		inprims= inprimS;
-		return;			/* don't echo */
+		return(0);		/* don't echo */
 	}
 	if (isexpos(s)) {		/* picture exposure */
 		inpexp *= exposval(s);
-		return;			/* don't echo */
+		return(0);		/* don't echo */
 	}
 	if (isaspect(s))		/* pixel aspect ratio */
 		pixaspect *= aspectval(s);
 	if (isview(s))			/* image view */
 		gotview += sscanview(&ourview, s);
-	fputs(s, stdout);
+	return(fputs(s, stdout));
 }
 
 
@@ -205,7 +237,7 @@ getahead()			/* load picture header */
 		exit(1);
 	}
 	if (lumf == rgblum)
-		comprgb2xyzmat(inrgb2xyz, inprims);
+		comprgb2xyzWBmat(inrgb2xyz, inprims);
 	else if (mbcalfile != NULL) {
 		fprintf(stderr, "%s: macbethcal only works with RGB pictures\n",
 				progname);
@@ -236,35 +268,17 @@ mapimage()				/* map picture and send to stdout */
 {
 	COLOR	*scan;
 
-#ifdef DEBUG
-	fprintf(stderr, "%s: generating histogram...", progname);
-#endif
-	fovhist();			/* generate adaptation histogram */
-#ifdef DEBUG
-	fputs("done\n", stderr);
-#endif
+	comphist();			/* generate adaptation histogram */
 	check2do();			/* modify what2do flags */
-	if (what2do&DO_VEIL) {
-#ifdef DEBUG
-		fprintf(stderr, "%s: computing veiling...", progname);
+	if (what2do&DO_VEIL)
+		compveil();		/* compute veil image */
+	if (!(what2do&DO_LINEAR))
+		if (mkbrmap() < 0)	/* make tone map */
+			what2do |= DO_LINEAR;	/* failed! -- use linear */
+#if ADJ_VEIL
+		else if (what2do&DO_VEIL)
+			adjveil();	/* else adjust veil image */
 #endif
-		compveil();
-#ifdef DEBUG
-		fputs("done\n", stderr);
-#endif
-	}
-#ifdef DEBUG
-	fprintf(stderr, "%s: computing brightness mapping...", progname);
-#endif
-	if (!(what2do&DO_LINEAR) && mkbrmap() < 0) {	/* make tone map */
-		what2do |= DO_LINEAR;		/* use linear scaling */
-#ifdef DEBUG
-		fputs("failed!\n", stderr);
-	} else
-		fputs("done\n", stderr);
-#else
-	}
-#endif
 	if (what2do&DO_LINEAR) {
 		if (scalef <= FTINY) {
 			if (what2do&DO_HSENS)
@@ -274,16 +288,13 @@ mapimage()				/* map picture and send to stdout */
 				scalef = Lb(0.5*(Bldmax+Bldmin)) / Lb(bwavg);
 			scalef *= WHTEFFICACY/(inpexp*ldmax);
 		}
-#ifdef DEBUG
-		fprintf(stderr, "%s: linear scaling factor = %f\n",
-				progname, scalef);
-#endif
 		fputexpos(inpexp*scalef, stdout);	/* record exposure */
 		if (lumf == cielum) scalef /= WHTEFFICACY;
 	}
-	putchar('\n');			/* complete header */
+	fputformat(COLRFMT, stdout);	/* complete header */
+	putchar('\n');
 	fputsresolu(&inpres, stdout);	/* resolution doesn't change */
-
+					/* condition our image */
 	for (scan = firstscan(); scan != NULL; scan = nextscan())
 		if (fwritescan(scan, scanlen(&inpres), stdout) < 0) {
 			fprintf(stderr, "%s: scanline write error\n",
@@ -293,38 +304,25 @@ mapimage()				/* map picture and send to stdout */
 }
 
 
-double
-centprob(x, y)			/* center-weighting probability function */
-int	x, y;
+getfovimg()			/* load foveal sampled image */
 {
-	double	xr, yr;
-
-	xr = (x+.5)/fvxr - .5;
-	yr = (y+.5)/fvyr - .5;
-	return(1. - xr*xr - yr*yr);	/* radial, == 0.5 at corners */
-}
-
-
-fovhist()			/* create foveal sampled image and histogram */
-{
 	extern FILE	*popen();
 	char	combuf[128];
-	double	l, b, lwmin, lwmax;
 	FILE	*fp;
 	int	x, y;
-
+						/* compute image size */
 	fvxr = sqrt(ourview.hn2)/FOVDIA + 0.5;
 	if (fvxr < 2) fvxr = 2;
 	fvyr = sqrt(ourview.vn2)/FOVDIA + 0.5;
 	if (fvyr < 2) fvyr = 2;
-	if (!(inpres.or & YMAJOR)) {		/* picture is rotated? */
+	if (!(inpres.rt & YMAJOR)) {		/* picture is rotated? */
 		y = fvyr;
 		fvyr = fvxr;
 		fvxr = y;
 	}
 	if ((fovimg = (COLOR *)malloc(fvxr*fvyr*sizeof(COLOR))) == NULL)
 		syserror("malloc");
-	sprintf(combuf, "pfilt -1 -b -x %d -y %d %s", fvxr, fvyr, infn);
+	sprintf(combuf, "pfilt -1 -b -pa 0 -x %d -y %d %s", fvxr, fvyr, infn);
 	if ((fp = popen(combuf, "r")) == NULL)
 		syserror("popen");
 	getheader(fp, NULL, NULL);	/* skip header */
@@ -334,34 +332,6 @@ fovhist()			/* create foveal sampled image and histogr
 		if (freadscan(fovscan(y), fvxr, fp) < 0)
 			goto readerr;
 	pclose(fp);
-	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;
-		}
-	if (lwmin < LMIN) lwmin = LMIN;
-	if (lwmax > LMAX) lwmax = LMAX;
-					/* compute histogram */
-	bwmin = Bl(lwmin)*(1. - .01/HISTRES);
-	bwmax = Bl(lwmax)*(1. + .01/HISTRES);
-	bwavg = 0.;
-	for (y = 0; y < fvyr; y++)
-		for (x = 0; x < fvxr; x++) {
-			if (what2do & DO_CWEIGHT &&
-					frandom() > centprob(x,y))
-				continue;
-			l = plum(fovscan(y)[x]);
-			b = Bl(l);
-			if (b < bwmin) continue;
-			if (b > bwmax) continue;
-			bwavg += b;
-			bwhc(b)++;
-			histot++;
-		}
-	bwavg /= (double)histot;
 	return;
 readerr:
 	fprintf(stderr, "%s: error reading from pfilt process in fovimage\n",
@@ -372,22 +342,21 @@ readerr:
 
 check2do()		/* check histogram to see what isn't worth doing */
 {
-	long	sum;
+	double	sum;
 	double	b, l;
 	register int	i;
 
 					/* check for within display range */
-	l = Lb(bwmax)/Lb(bwmin);
-	if (l <= ldmax/ldmin)
+	if (bwmax - bwmin <= Bldmax - Bldmin)
 		what2do |= DO_LINEAR;
 					/* determine if veiling significant */
-	if (l < 100.)			/* heuristic */
+	if (bwmax - bwmin < 4.5)		/* heuristic */
 		what2do &= ~DO_VEIL;
 
 	if (!(what2do & (DO_ACUITY|DO_COLOR)))
 		return;
 					/* find 5th percentile */
-	sum = histot*0.05 + .5;
+	sum = histot*0.05;
 	for (i = 0; i < HISTRES; i++)
 		if ((sum -= bwhist[i]) <= 0)
 			break;
@@ -396,9 +365,9 @@ check2do()		/* check histogram to see what isn't worth
 					/* determine if acuity adj. useful */
 	if (what2do&DO_ACUITY &&
 			hacuity(l) >= (inpres.xr/sqrt(ourview.hn2) +
-			inpres.yr/sqrt(ourview.vn2))/(2.*180./PI*2.))
+			inpres.yr/sqrt(ourview.vn2))/(2.*180./PI))
 		what2do &= ~DO_ACUITY;
 					/* color sensitivity loss? */
-	if (l >= 6.0)
+	if (l >= TopMesopic)
 		what2do &= ~DO_COLOR;
 }