--- ray/src/px/pf3.c	1993/06/25 17:08:40	2.8
+++ ray/src/px/pf3.c	1996/04/04 13:01:35	2.12
@@ -1,4 +1,4 @@
-/* Copyright (c) 1992 Regents of the University of California */
+/* Copyright (c) 1993 Regents of the University of California */
 
 #ifndef lint
 static char SCCSid[] = "$SunId$ LBL";
@@ -14,7 +14,8 @@ static char SCCSid[] = "$SunId$ LBL";
 
 #include  "color.h"
 
-#define  TEPS		0.25	/* threshold proximity goal */
+#define  TEPS		0.2	/* threshold proximity goal */
+#define  REPS		0.1	/* radius proximity goal */
 
 extern double  CHECKRAD;	/* radius over which gaussian is summed */
 
@@ -42,6 +43,8 @@ extern float  **greybar;	/* grey-averaged input values
 extern int  obarsize;		/* size of output scan bar */
 extern int  orad;		/* output window radius */
 
+extern int  wrapfilt;		/* wrap filter horizontally? */
+
 extern char  *progname;
 
 float  *gausstable;		/* gauss lookup table */
@@ -51,6 +54,8 @@ short  *ringwt;			/* weight (count) of ring values */
 short  *ringndx;		/* ring index table */
 float  *warr;			/* array of pixel weights */
 
+extern double  (*ourbright)();	/* brightness computation function */
+
 double  pickfilt();
 
 #define	 lookgauss(x)		gausstable[(int)(10.*(x)+.5)]
@@ -111,7 +116,7 @@ int  c, r;
 	int  wsum;
 	double  d;
 	int  y;
-	register int  x;
+	register int  x, offs;
 	register COLOR	*scan;
 	
 	wsum = 0;
@@ -124,13 +129,14 @@ int  c, r;
 		if (d >= 0.5) break;
 		scan = scanin[y%barsize];
 		for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
-			if (x < 0) continue;
-			if (x >= xres) break;
+			offs = x < 0 ? xres : x >= xres ? -xres : 0;
+			if (offs && !wrapfilt)
+				continue;
 			d = x_c < 1.0 ? x_c*x - c : (double)(x - xcent);
 			if (d < -0.5) continue;
 			if (d >= 0.5) break;
 			wsum++;
-			addcolor(csum, scan[x]);
+			addcolor(csum, scan[x+offs]);
 		}
 	}
 	if (wsum > 1)
@@ -146,7 +152,7 @@ int  c, r;
 	double  dy, dx, weight, wsum;
 	COLOR  ctmp;
 	int  y;
-	register int  x;
+	register int  x, offs;
 	register COLOR	*scan;
 
 	wsum = FTINY;
@@ -157,12 +163,13 @@ int  c, r;
 		dy = (y_r*(y+.5) - (r+.5))/rad;
 		scan = scanin[y%barsize];
 		for (x = xcent-xrad; x <= xcent+xrad; x++) {
-			if (x < 0) continue;
-			if (x >= xres) break;
+			offs = x < 0 ? xres : x >= xres ? -xres : 0;
+			if (offs && !wrapfilt)
+				continue;
 			dx = (x_c*(x+.5) - (c+.5))/rad;
 			weight = lookgauss(dx*dx + dy*dy);
 			wsum += weight;
-			copycolor(ctmp, scan[x]);
+			copycolor(ctmp, scan[x+offs]);
 			scalecolor(ctmp, weight);
 			addcolor(csum, ctmp);
 		}
@@ -176,10 +183,9 @@ int  xcent, ycent;
 int  ccent, rcent;
 {
 	double  d;
-	int  r, y;
+	int  r, y, offs;
 	register int  c, x;
 	register float  *gscan;
-#define pval gscan
 					/* compute ring sums */
 	bzero((char *)ringsum, (orad+1)*sizeof(float));
 	bzero((char *)ringwt, (orad+1)*sizeof(short));
@@ -188,11 +194,13 @@ int  ccent, rcent;
 		if (rcent+r >= nrows) break;
 		gscan = greybar[(rcent+r)%obarsize];
 		for (c = -orad; c <= orad; c++) {
-			if (ccent+c < 0) continue;
-			if (ccent+c >= ncols) break;
+			offs = ccent+c < 0 ? ncols :
+					ccent+c >= ncols ? -ncols : 0;
+			if (offs && !wrapfilt)
+				continue;
 			x = ringndx[c*c + r*r];
 			if (x < 0) continue;
-			ringsum[x] += gscan[ccent+c];
+			ringsum[x] += gscan[ccent+c+offs];
 			ringwt[x]++;
 		}
 	}
@@ -204,16 +212,16 @@ int  ccent, rcent;
 		if (d < -0.5) continue;
 		if (d >= 0.5) break;
 		for (x = xcent+1-xbrad; x <= xcent+xbrad; x++) {
-			if (x < 0) continue;
-			if (x >= xres) break;
+			offs = x < 0 ? xres : x >= xres ? -xres : 0;
+			if (offs && !wrapfilt)
+				continue;
 			d = x_c < 1.0 ? x_c*x - ccent : (double)(x - xcent);
 			if (d < -0.5) continue;
 			if (d >= 0.5) break;
-			pval = scanin[y%barsize][x];
-			sumans(x, y, rcent, ccent, pickfilt(bright(pval)));
+			sumans(x, y, rcent, ccent,
+			pickfilt((*ourbright)(scanin[y%barsize][x+offs])));
 		}
 	}
-#undef pval
 }
 
 
@@ -222,10 +230,12 @@ pickfilt(p0)			/* find filter multiplier for p0 */
 double  p0;
 {
 	double  m = 1.0;
-	double  t, avg, wsum;
-	int  ilimit = 5/TEPS;
+	double  t, num, denom, avg, wsum;
+	double  mlimit[2];
+	int  ilimit = 4/TEPS;
 	register int  i;
 				/* iterative search for m */
+	mlimit[0] = 1.0; mlimit[1] = orad/rad/CHECKRAD;
 	do {
 					/* compute grey weighted average */
 		i = 1.25*CHECKRAD*rad*m + .5;
@@ -237,21 +247,39 @@ double  p0;
 			avg += t*ringsum[i];
 			wsum += t*ringwt[i];
 		}
+		if (avg < 1e-20)		/* zero inclusive average */
+			return(1.0);
 		avg /= wsum;
 					/* check threshold */
-		t = p0 - avg;
-		if (t < 0.0) t = -t;
-		t *= gausstable[0]/(m*m*avg);
-		if (t <= thresh && (m <= 1.0+FTINY ||
-				(thresh-t)/thresh <= TEPS))
-			break;
-		if (t > thresh && (m*CHECKRAD*rad >= orad-FTINY ||
-				(t-thresh)/thresh <= TEPS))
-			break;
+		denom = m*m/gausstable[0] - p0/avg;
+		if (denom <= FTINY) {		/* zero exclusive average */
+			if (m >= mlimit[1]-REPS)
+				break;
+			m = mlimit[1];
+			continue;
+		}
+		num = p0/avg - 1.0;
+		if (num < 0.0) num = -num;
+		t = num/denom;
+		if (t <= thresh) {
+			if (m <= mlimit[0]+REPS || (thresh-t)/thresh <= TEPS)
+				break;
+		} else {
+			if (m >= mlimit[1]-REPS || (t-thresh)/thresh <= TEPS)
+				break;
+		}
+		t = m;			/* remember current m */
 					/* next guesstimate */
-		m *= sqrt(t/thresh);
-		if (m < 1.0) m = 1.0;
-		else if (m*CHECKRAD*rad > orad) m = orad/rad/CHECKRAD;
+		m = sqrt(gausstable[0]*(num/thresh + p0/avg));
+		if (m < t) {		/* bound it */
+			if (m <= mlimit[0]+FTINY)
+				m = 0.5*(mlimit[0] + t);
+			mlimit[1] = t;
+		} else {
+			if (m >= mlimit[1]-FTINY)
+				m = 0.5*(mlimit[1] + t);
+			mlimit[0] = t;
+		}
 	} while (--ilimit > 0);
 	return(m);
 }
@@ -264,12 +292,16 @@ double  m;
 {
 	double  dy, dx;
 	COLOR  pval, ctmp;
-	int  ksiz, r;
+	int  ksiz, r, offs;
 	double  pc, pr, norm;
 	register int  i, c;
 	register COLOR	*scan;
-
-	copycolor(pval, scanin[py%barsize][px]);
+	/*
+	 * This normalization method fails at the picture borders because
+	 * a different number of input pixels contribute there.
+	 */
+	scan = scanin[py%barsize] + (px < 0 ? xres : px >= xres ? -xres : 0);
+	copycolor(pval, scan[px]);
 	pc = x_c*(px+.5);
 	pr = y_r*(py+.5);
 	ksiz = CHECKRAD*m*rad + 1;
@@ -282,8 +314,10 @@ double  m;
 		if (r >= nrows) break;
 		dy = (pr - (r+.5))/(m*rad);
 		for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
-			if (c < 0) continue;
-			if (c >= ncols) break;
+			if (!wrapfilt) {
+				if (c < 0) continue;
+				if (c >= ncols) break;
+			}
 			dx = (pc - (c+.5))/(m*rad);
 			norm += warr[i++] = lookgauss(dx*dx + dy*dy);
 		}
@@ -298,12 +332,13 @@ double  m;
 		if (r >= nrows) break;
 		scan = scoutbar[r%obarsize];
 		for (c = ccent-ksiz; c <= ccent+ksiz; c++) {
-			if (c < 0) continue;
-			if (c >= ncols) break;
+			offs = c < 0 ? ncols : c >= ncols ? -ncols : 0;
+			if (offs && !wrapfilt)
+				continue;
 			copycolor(ctmp, pval);
 			dx = norm*warr[i++];
 			scalecolor(ctmp, dx);
-			addcolor(scan[c], ctmp);
+			addcolor(scan[c+offs], ctmp);
 		}
 	}
 }