--- ray/src/px/pcond4.c	1996/10/08 12:48:00	3.8
+++ ray/src/px/pcond4.c	2004/11/08 15:50:59	3.18
@@ -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: pcond4.c,v 3.18 2004/11/08 15:50:59 greg Exp $";
 #endif
-
 /*
  * Routines for veiling glare and loss of acuity.
  */
@@ -14,21 +11,22 @@ static char SCCSid[] = "$SunId$ LBL";
 
 #define	VADAPT		0.08		/* fraction of adaptation from veil */
 
-extern COLOR	*fovimg;		/* foveal (1 degree) averaged image */
-extern short	fvxr, fvyr;		/* foveal image resolution */
+static COLOR	*veilimg = NULL;	/* veiling image */
 
-#define fovscan(y)	(fovimg+(y)*fvxr)
-
-static COLOR	*veilimg;		/* veiling image */
-
 #define veilscan(y)	(veilimg+(y)*fvxr)
 
 static float	(*raydir)[3] = NULL;	/* ray direction for each pixel */
 
 #define rdirscan(y)	(raydir+(y)*fvxr)
 
+static void compraydir(void);
+#if ADJ_VEIL
+static void smoothveil(void);
+#endif
 
-compraydir()				/* compute ray directions */
+
+static void
+compraydir(void)				/* compute ray directions */
 {
 	FVECT	rorg, rdir;
 	double	h, v;
@@ -41,7 +39,7 @@ compraydir()				/* compute ray directions */
 		syserror("malloc");
 
 	for (y = 0; y < fvyr; y++) {
-		switch (inpres.or) {
+		switch (inpres.rt) {
 		case YMAJOR: case YMAJOR|XDECR:
 			v = (y+.5)/fvyr; break;
 		case YMAJOR|YDECR: case YMAJOR|YDECR|XDECR:
@@ -52,7 +50,7 @@ compraydir()				/* compute ray directions */
 			h = 1. - (y+.5)/fvyr; break;
 		}
 		for (x = 0; x < fvxr; x++) {
-			switch (inpres.or) {
+			switch (inpres.rt) {
 			case YMAJOR: case YMAJOR|YDECR:
 				h = (x+.5)/fvxr; break;
 			case YMAJOR|XDECR: case YMAJOR|XDECR|YDECR:
@@ -77,12 +75,16 @@ compraydir()				/* compute ray directions */
 }
 
 
-compveil()				/* compute veiling image */
+extern void
+compveil(void)				/* compute veiling image */
 {
 	double	t2, t2sum;
 	COLOR	ctmp, vsum;
 	int	px, py;
 	register int	x, y;
+
+	if (veilimg != NULL)		/* already done? */
+		return;
 					/* compute ray directions */
 	compraydir();
 					/* compute veil image */
@@ -100,26 +102,164 @@ compveil()				/* compute veiling image */
 							rdirscan(y)[x]);
 					if (t2 <= FTINY) continue;
 					/*	use approximation instead
-					t2 = acos(t2);
-					t2 = 1./(t2*t2);
+					t3 = acos(t2);
+					t2 = t2/(t3*t3);
 					*/
-					t2 = .5 / (1. - t2);
+					t2 *= .5 / (1. - t2);
 					copycolor(ctmp, fovscan(y)[x]);
 					scalecolor(ctmp, t2);
 					addcolor(vsum, ctmp);
 					t2sum += t2;
 				}
 			/* VADAPT of original is subtracted in addveil() */
-			scalecolor(vsum, VADAPT/t2sum);
+			if (t2sum > FTINY)
+				scalecolor(vsum, VADAPT/t2sum);
 			copycolor(veilscan(py)[px], vsum);
 		}
+					/* modify FOV sample image */
+	for (y = 0; y < fvyr; y++)
+		for (x = 0; x < fvxr; x++) {
+			scalecolor(fovscan(y)[x], 1.-VADAPT);
+			addcolor(fovscan(y)[x], veilscan(y)[x]);
+		}
+	comphist();			/* recompute histogram */
 }
 
 
-addveil(sl, y)				/* add veil to scanline */
-COLOR	*sl;
-int	y;
+#if ADJ_VEIL
+/*
+ * The following veil adjustment was added to compensate for
+ * the fact that contrast reduction gets confused with veil
+ * in the human visual system.  Therefore, we reduce the
+ * veil in portions of the image where our mapping has
+ * already reduced contrast below the target value.
+ * This gets called after the intial veil has been computed
+ * and added to the foveal image, and the mapping has been
+ * determined.
+ */
+extern void
+adjveil(void)				/* adjust veil image */
 {
+	float	*crfptr = crfimg;
+	COLOR	*fovptr = fovimg;
+	COLOR	*veilptr = veilimg;
+	double	s2nits = 1./inpexp;
+	double	vl, vl2, fovl, vlsum;
+	double	deltavc[3];
+	int	i, j;
+
+	if (lumf == rgblum)
+		s2nits *= WHTEFFICACY;
+
+	for (i = fvxr*fvyr; i--; crfptr++, fovptr++, veilptr++) {
+		if (crfptr[0] >= 0.95)
+			continue;
+		vl = plum(veilptr[0]);
+		fovl = (plum(fovptr[0]) - vl) * (1./(1.-VADAPT));
+		if (vl <= 0.05*fovl)
+			continue;
+		vlsum = vl;
+		for (j = 2; j < 11; j++) {
+			vlsum += crfptr[0]*vl - (1.0 - crfptr[0])*fovl;
+			vl2 = vlsum / (double)j;
+			if (vl2 < 0.0)
+				vl2 = 0.0;
+			crfptr[0] = crfactor(fovl + vl2);
+		}
+		/* desaturation code causes color fringes at this level */
+		for (j = 3; j--; ) {
+			double	vc = colval(veilptr[0],j);
+			double	fovc = (colval(fovptr[0],j) - vc) *
+						(1./(1.-VADAPT));
+			deltavc[j] = (1.-crfptr[0])*(fovl/s2nits - fovc);
+			if (vc + deltavc[j] < 0.0)
+				break;
+		}
+		if (j < 0)
+			addcolor(veilptr[0], deltavc);
+		else
+			scalecolor(veilptr[0], vl2/vl);
+	}
+	smoothveil();			/* smooth our result */
+}
+
+
+static void
+smoothveil(void)				/* smooth veil image */
+{
+	COLOR 	*nveilimg;
+	COLOR	*ovptr, *nvptr;
+	int	x, y, i;
+
+	nveilimg = (COLOR *)malloc(fvxr*fvyr*sizeof(COLOR));
+	if (nveilimg == NULL)
+		return;
+	for (y = 1; y < fvyr-1; y++) {
+		ovptr = veilimg + y*fvxr + 1;
+		nvptr = nveilimg + y*fvxr + 1;
+		for (x = 1; x < fvxr-1; x++, ovptr++, nvptr++)
+			for (i = 3; i--; )
+				nvptr[0][i] = 0.5 * ovptr[0][i]
+					+ (1./12.) *
+				(ovptr[-1][i] + ovptr[-fvxr][i] +
+					ovptr[1][i] + ovptr[fvxr][i])
+					+ (1./24.) *
+				(ovptr[-fvxr-1][i] + ovptr[-fvxr+1][i] +
+					ovptr[fvxr-1][i] + ovptr[fvxr+1][i]);
+	}
+	ovptr = veilimg + 1;
+	nvptr = nveilimg + 1;
+	for (x = 1; x < fvxr-1; x++, ovptr++, nvptr++)
+		for (i = 3; i--; )
+			nvptr[0][i] = 0.5 * ovptr[0][i]
+				+ (1./9.) *
+			(ovptr[-1][i] + ovptr[1][i] + ovptr[fvxr][i])
+				+ (1./12.) *
+			(ovptr[fvxr-1][i] + ovptr[fvxr+1][i]);
+	ovptr = veilimg + (fvyr-1)*fvxr + 1;
+	nvptr = nveilimg + (fvyr-1)*fvxr + 1;
+	for (x = 1; x < fvxr-1; x++, ovptr++, nvptr++)
+		for (i = 3; i--; )
+			nvptr[0][i] = 0.5 * ovptr[0][i]
+				+ (1./9.) *
+			(ovptr[-1][i] + ovptr[1][i] + ovptr[-fvxr][i])
+				+ (1./12.) *
+			(ovptr[-fvxr-1][i] + ovptr[-fvxr+1][i]);
+	ovptr = veilimg + fvxr;
+	nvptr = nveilimg + fvxr;
+	for (y = 1; y < fvyr-1; y++, ovptr += fvxr, nvptr += fvxr)
+		for (i = 3; i--; )
+			nvptr[0][i] = 0.5 * ovptr[0][i]
+				+ (1./9.) *
+			(ovptr[-fvxr][i] + ovptr[1][i] + ovptr[fvxr][i])
+				+ (1./12.) *
+			(ovptr[-fvxr+1][i] + ovptr[fvxr+1][i]);
+	ovptr = veilimg + fvxr - 1;
+	nvptr = nveilimg + fvxr - 1;
+	for (y = 1; y < fvyr-1; y++, ovptr += fvxr, nvptr += fvxr)
+		for (i = 3; i--; )
+			nvptr[0][i] = 0.5 * ovptr[0][i]
+				+ (1./9.) *
+			(ovptr[-fvxr][i] + ovptr[-1][i] + ovptr[fvxr][i])
+				+ (1./12.) *
+			(ovptr[-fvxr-1][i] + ovptr[fvxr-1][i]);
+	for (i = 3; i--; ) {
+		nveilimg[0][i] = veilimg[0][i];
+		nveilimg[fvxr-1][i] = veilimg[fvxr-1][i];
+		nveilimg[(fvyr-1)*fvxr][i] = veilimg[(fvyr-1)*fvxr][i];
+		nveilimg[fvyr*fvxr-1][i] = veilimg[fvyr*fvxr-1][i];
+	}
+	free((void *)veilimg);
+	veilimg = nveilimg;
+}
+#endif
+
+extern void
+addveil(				/* add veil to scanline */
+	COLOR	*sl,
+	int	y
+)
+{
 	int	vx, vy;
 	double	dx, dy;
 	double	lv, uv;
@@ -158,47 +298,30 @@ typedef struct {
 
 SCANBAR	*rootbar;		/* root scan bar (lowest resolution) */
 
-float	*inpacuD;		/* input acuity data (cycles/degree) */
+float	*inpacuD = NULL;	/* input acuity data (cycles/degree) */
 
 #define tsampr(x,y)	inpacuD[(y)*fvxr+(x)]
 
+static COLOR * getascan(SCANBAR	*sb, int	y);
+static void acusample(COLOR	col, int	x, int	y, double	sr);
+static void ascanval(COLOR	col, int	x, int	y, SCANBAR	*sb);
+static SCANBAR	*sballoc(int	se, int	ns, int	sl);
 
-double
-hacuity(La)			/* return visual acuity in cycles/degree */
-double	La;
-{				/* data due to S. Shaler (we should fit it!) */
-#define NPOINTS	20
-	static float	l10lum[NPOINTS] = {
-		-3.10503,-2.66403,-2.37703,-2.09303,-1.64403,-1.35803,
-		-1.07403,-0.67203,-0.38503,-0.10103,0.29397,0.58097,0.86497,
-		1.25697,1.54397,1.82797,2.27597,2.56297,2.84697,3.24897
-	};
-	static float	resfreq[NPOINTS] = {
-		2.09,3.28,3.79,4.39,6.11,8.83,10.94,18.66,23.88,31.05,37.42,
-		37.68,41.60,43.16,45.30,47.00,48.43,48.32,51.06,51.09
-	};
-	double	l10La;
-	register int	i;
-					/* check limits */
-	if (La <= 7.85e-4)
-		return(resfreq[0]);
-	if (La >= 1.78e3)
-		return(resfreq[NPOINTS-1]);
-					/* interpolate data */
-	l10La = log10(La);
-	for (i = 0; i < NPOINTS-2 && l10lum[i+1] <= l10La; i++)
-		;
-	return( ( (l10lum[i+1] - l10La)*resfreq[i] +
-			(l10La - l10lum[i])*resfreq[i+1] ) /
-			(l10lum[i+1] - l10lum[i]) );
-#undef NPOINTS
+extern double
+hacuity(			/* return visual acuity in cycles/degree */
+	double	La
+)
+{
+					/* functional fit */
+	return(17.25*atan(1.4*log10(La) + 0.35) + 25.72);
 }
 
 
-COLOR *
-getascan(sb, y)			/* find/read scanline y for scanbar sb */
-register SCANBAR	*sb;
-int	y;
+static COLOR *
+getascan(			/* find/read scanline y for scanbar sb */
+	register SCANBAR	*sb,
+	int	y
+)
 {
 	register COLOR	*sl0, *sl1, *mysl;
 	register int	i;
@@ -231,22 +354,27 @@ int	y;
 }
 
 
-acuscan(scln, y)		/* get acuity-sampled scanline */
-COLOR	*scln;
-int	y;
+extern void
+acuscan(		/* get acuity-sampled scanline */
+	COLOR	*scln,
+	int	y
+)
 {
 	double	sr;
 	double	dx, dy;
 	int	ix, iy;
 	register int	x;
+	
+	if (inpacuD == NULL)
+		return;
 					/* compute foveal y position */
 	iy = dy = (y+.5)/numscans(&inpres)*fvyr - .5;
-	if (iy >= fvyr-1) iy--;
+	while (iy >= fvyr-1) iy--;
 	dy -= (double)iy;
 	for (x = 0; x < scanlen(&inpres); x++) {
 					/* compute foveal x position */
 		ix = dx = (x+.5)/scanlen(&inpres)*fvxr - .5;
-		if (ix >= fvxr-1) ix--;
+		while (ix >= fvxr-1) ix--;
 		dx -= (double)ix;
 					/* interpolate sample rate */
 		sr = (1.-dy)*((1.-dx)*tsampr(ix,iy) + dx*tsampr(ix+1,iy)) +
@@ -257,10 +385,13 @@ int	y;
 }
 
 
-acusample(col, x, y, sr)	/* interpolate sample at (x,y) using rate sr */
-COLOR	col;
-int	x, y;
-double	sr;
+static void
+acusample(	/* interpolate sample at (x,y) using rate sr */
+	COLOR	col,
+	int	x,
+	int	y,
+	double	sr
+)
 {
 	COLOR	c1;
 	double	d;
@@ -279,10 +410,13 @@ double	sr;
 }
 
 
-ascanval(col, x, y, sb)		/* interpolate scanbar at orig. coords (x,y) */
-COLOR	col;
-int	x, y;
-SCANBAR	*sb;
+static void
+ascanval(		/* interpolate scanbar at orig. coords (x,y) */
+	COLOR	col,
+	int	x,
+	int	y,
+	SCANBAR	*sb
+)
 {
 	COLOR	*sl0, *sl1, c1, c1y;
 	double	dx, dy;
@@ -302,11 +436,10 @@ SCANBAR	*sb;
 	dy -= (double)iy;
 					/* get scanlines */
 	sl0 = getascan(sb, iy);
-	if (sl0 == NULL) {
-		fprintf(stderr, "%s: internal - cannot backspace in ascanval\n",
-				progname);
-		exit(1);
-	}
+#ifdef DEBUG
+	if (sl0 == NULL)
+		error(INTERNAL, "cannot backspace in ascanval");
+#endif
 	sl1 = getascan(sb, iy+1);
 					/* 2D linear interpolation */
 	copycolor(col, sl0[ix]);
@@ -322,41 +455,51 @@ SCANBAR	*sb;
 	scalecolor(col, 1.-dy);
 	scalecolor(c1y, dy);
 	addcolor(col, c1y);
+	for (ix = 0; ix < 3; ix++)	/* make sure no negative */
+		if (colval(col,ix) < 0.)
+			colval(col,ix) = 0.;
 }
 
 
-SCANBAR	*
-sballoc(sr, ns, sl)		/* allocate scanbar */
-int	sr;		/* sampling rate */
-int	ns;		/* number of scanlines */
-int	sl;		/* original scanline length */
+static SCANBAR	*
+sballoc(		/* allocate scanbar */
+	int	se,		/* sampling rate exponent */
+	int	ns,		/* number of scanlines */
+	int	sl		/* original scanline length */
+)
 {
 	SCANBAR	*sbarr;
 	register SCANBAR	*sb;
 
-	sbarr = sb = (SCANBAR *)malloc((sr+1)*sizeof(SCANBAR));
+	sbarr = sb = (SCANBAR *)malloc((se+1)*sizeof(SCANBAR));
 	if (sb == NULL)
 		syserror("malloc");
 	do {
-		sb->sdata = (COLOR *)malloc((sl>>sr)*ns*sizeof(COLOR));
+		sb->len = sl>>se;
+		if (sb->len <= 0)
+			continue;
+		sb->sampe = se;
+		sb->nscans = ns;
+		sb->sdata = (COLOR *)malloc(sb->len*ns*sizeof(COLOR));
 		if (sb->sdata == NULL)
 			syserror("malloc");
-		sb->sampe = sr;
-		sb->nscans = ns;
-		sb->len = sl>>sr;
 		sb->nread = 0;
 		ns <<= 1;
 		sb++;
-	} while (--sr >= 0);
+	} while (--se >= 0);
 	return(sbarr);
 }
 
 
-initacuity()			/* initialize variable acuity sampling */
+extern int
+initacuity(void)			/* initialize variable acuity sampling */
 {
 	FVECT	diffx, diffy, cp;
 	double	omega, maxsr;
 	register int	x, y, i;
+	
+	if (fvxr < 3 || fvyr < 3)
+		return(0);		/* too small to work with */
 
 	compraydir();			/* compute ray directions */
 
@@ -376,7 +519,7 @@ initacuity()			/* initialize variable acuity sampling 
 			}
 			fcross(cp, diffx, diffy);
 			omega = 0.5 * sqrt(DOT(cp,cp));
-			if (omega <= FTINY)
+			if (omega <= FTINY*FTINY)
 				tsampr(x,y) = 1.;
 			else if ((tsampr(x,y) = PI/180. / sqrt(omega) /
 					hacuity(plum(fovscan(y)[x]))) > maxsr)
@@ -388,9 +531,10 @@ initacuity()			/* initialize variable acuity sampling 
 		tsampr(x,fvyr-1) = tsampr(x,fvyr-2);
 	}
 	for (y = 0; y < fvyr; y++) {
-		tsampr(y,0) = tsampr(y,1);
-		tsampr(y,fvxr-1) = tsampr(y,fvxr-2);
+		tsampr(0,y) = tsampr(1,y);
+		tsampr(fvxr-1,y) = tsampr(fvxr-2,y);
 	}
 					/* initialize with next power of two */
 	rootbar = sballoc((int)(log(maxsr)/log(2.))+1, 2, scanlen(&inpres));
+	return(1);
 }