--- ray/src/px/pcond4.c	1997/01/09 13:56:32	3.11
+++ ray/src/px/pcond4.c	2003/02/22 02:07:27	3.16
@@ -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.16 2003/02/22 02:07:27 greg Exp $";
 #endif
-
 /*
  * Routines for veiling glare and loss of acuity.
  */
@@ -36,7 +33,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:
@@ -47,7 +44,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:
@@ -98,17 +95,18 @@ 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 */
@@ -121,6 +119,132 @@ compveil()				/* compute veiling image */
 }
 
 
+#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.
+ */
+adjveil()				/* 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 */
+}
+
+
+smoothveil()				/* 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
+
 addveil(sl, y)				/* add veil to scanline */
 COLOR	*sl;
 int	y;
@@ -171,32 +295,9 @@ float	*inpacuD;		/* input acuity data (cycles/degree) 
 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
+{
+					/* functional fit */
+	return(17.25*atan(1.4*log10(La) + 0.35) + 25.72);
 }
 
 
@@ -308,11 +409,8 @@ SCANBAR	*sb;
 					/* get scanlines */
 	sl0 = getascan(sb, iy);
 #ifdef DEBUG
-	if (sl0 == NULL) {
-		fprintf(stderr, "%s: internal - cannot backspace in ascanval\n",
-				progname);
-		abort();
-	}
+	if (sl0 == NULL)
+		error(INTERNAL, "cannot backspace in ascanval");
 #endif
 	sl1 = getascan(sb, iy+1);
 					/* 2D linear interpolation */
@@ -348,8 +446,10 @@ int	sl;		/* original scanline length */
 	if (sb == NULL)
 		syserror("malloc");
 	do {
-		sb->sampe = se;
 		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)