--- ray/src/px/pcond4.c	1997/01/28 16:31:17	3.14
+++ ray/src/px/pcond4.c	2004/11/14 16:57:18	3.19
@@ -1,9 +1,6 @@
-/* Copyright (c) 1997 Regents of the University of California */
-
 #ifndef lint
-static char SCCSid[] = "$SunId$ LBL";
+static const char	RCSid[] = "$Id: pcond4.c,v 3.19 2004/11/14 16:57:18 greg Exp $";
 #endif
-
 /*
  * Routines for veiling glare and loss of acuity.
  */
@@ -22,8 +19,14 @@ static float	(*raydir)[3] = NULL;	/* ray direction for
 
 #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;
@@ -36,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:
@@ -47,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:
@@ -72,7 +75,8 @@ compraydir()				/* compute ray directions */
 }
 
 
-compveil()				/* compute veiling image */
+extern void
+compveil(void)				/* compute veiling image */
 {
 	double	t2, t2sum;
 	COLOR	ctmp, vsum;
@@ -108,7 +112,8 @@ compveil()				/* compute veiling image */
 					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,10 +126,140 @@ compveil()				/* compute veiling image */
 }
 
 
-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;
@@ -163,24 +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;
+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;
@@ -213,14 +354,19 @@ 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;
 	while (iy >= fvyr-1) iy--;
@@ -239,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;
@@ -261,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;
@@ -309,11 +461,12 @@ SCANBAR	*sb;
 }
 
 
-SCANBAR	*
-sballoc(se, ns, sl)		/* allocate scanbar */
-int	se;		/* sampling rate exponent */
-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;
@@ -322,8 +475,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)
@@ -336,12 +491,13 @@ int	sl;		/* original scanline length */
 }
 
 
-initacuity()			/* initialize variable acuity sampling */
+extern void
+initacuity(void)			/* initialize variable acuity sampling */
 {
 	FVECT	diffx, diffy, cp;
 	double	omega, maxsr;
 	register int	x, y, i;
-
+	
 	compraydir();			/* compute ray directions */
 
 	inpacuD = (float *)malloc(fvxr*fvyr*sizeof(float));