--- ray/src/cv/bsdfrbf.c	2013/10/17 19:09:11	2.9
+++ ray/src/cv/bsdfrbf.c	2013/10/19 00:11:50	2.11
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdfrbf.c,v 2.9 2013/10/17 19:09:11 greg Exp $";
+static const char RCSid[] = "$Id: bsdfrbf.c,v 2.11 2013/10/19 00:11:50 greg Exp $";
 #endif
 /*
  * Radial basis function representation for BSDF data.
@@ -15,11 +15,17 @@ static const char RCSid[] = "$Id: bsdfrbf.c,v 2.9 2013
 #include "bsdfrep.h"
 
 #ifndef MINRSCA
-#define MINRSCA		0.15		/* minimum radius scaling factor */
+#define MINRSCA		1.0		/* minimum radius scaling factor */
 #endif
 #ifndef MAXRSCA
 #define MAXRSCA		2.7		/* maximum radius scaling factor */
 #endif
+#ifndef VARTHRESH
+#define VARTHRESH	0.0015		/* culling variance threshold */
+#endif
+#ifndef DIFFMAX2
+#define DIFFMAX2	(16.*VARTHRESH)	/* maximum ignored sample variance */
+#endif
 #ifndef MAXFRAC
 #define MAXFRAC		0.5		/* maximum contribution to neighbor */
 #endif
@@ -84,11 +90,6 @@ compute_radii(void)
 	double		ang2, lastang2;
 	int		r, i, j, jn, ii, jj, inear, jnear;
 
-	for (i = 0; i < GRIDRES; i++)		/* initialize minimum radii */
-	    for (j = 0; j < GRIDRES; j++)
-		if (dsf_grid[i][j].nval)
-		    dsf_grid[i][j].crad = cradmin;
-
 	r = GRIDRES/2;				/* proceed in zig-zag */
 	for (i = 0; i < GRIDRES; i++)
 	    for (jn = 0; jn < GRIDRES; jn++) {
@@ -127,6 +128,58 @@ compute_radii(void)
 						/* next search radius */
 		r = ang2*(2.*GRIDRES/M_PI) + 3;
 	    }
+	for (i = 0; i < GRIDRES; i++)		/* grow radii where uniform */
+	    for (j = 0; j < GRIDRES; j++) {
+		double	midmean = 0.0;
+		int	nsum = 0;
+		if (!dsf_grid[i][j].nval)
+			continue;
+		r = 1;				/* avg. immediate neighbors */
+		for (ii = i-r; ii <= i+r; ii++) {
+		    if (ii < 0) continue;
+		    if (ii >= GRIDRES) break;
+		    for (jj = j-r; jj <= j+r; jj++) {
+			if (jj < 0) continue;
+			if (jj >= GRIDRES) break;
+			midmean += dsf_grid[ii][jj].vsum;
+			nsum += dsf_grid[ii][jj].nval;
+		    }
+		}
+		midmean /= (double)nsum;
+		while (++r < GRIDRES) {		/* attempt to grow perimeter */
+		    double	diff2sum = 0.0;
+		    nsum = 0;
+		    for (ii = i-r; ii <= i+r; ii++) {
+			int	jstep = 1;
+			if (ii < 0) continue;
+			if (ii >= GRIDRES) break;
+			if ((i-r < ii) & (ii < i+r))
+				jstep = r<<1;
+			for (jj = j-r; jj <= j+r; jj += jstep) {
+				double	d2;
+				if (jj < 0) continue;
+				if (jj >= GRIDRES) break;
+				if (!dsf_grid[ii][jj].nval)
+					continue;
+				d2 = midmean - dsf_grid[ii][jj].vsum /
+						(double)dsf_grid[ii][jj].nval;
+				d2 *= d2;
+				if (d2 > DIFFMAX2*midmean*midmean)
+					goto escape;
+				diff2sum += d2;
+				++nsum;
+			}
+		    }
+		    if (diff2sum > VARTHRESH*midmean*midmean*(double)nsum)
+			break;
+		}
+escape:		--r;
+		r = ANG2R(r*(M_PI/MAXRSCA/GRIDRES));
+		if (r < cradmin)
+			r = cradmin;
+		if (dsf_grid[i][j].crad < r)
+			dsf_grid[i][j].crad = r;
+	    }
 						/* blur radii over hemisphere */
 	memset(fill_grid, 0, sizeof(fill_grid));
 	memset(fill_cnt, 0, sizeof(fill_cnt));
@@ -155,24 +208,39 @@ compute_radii(void)
 			dsf_grid[i][j].crad = fill_grid[i][j]/fill_cnt[i][j];
 }
 
+/* Radius comparison for qsort() */
+static int
+radius_cmp(const void *p1, const void *p2)
+{
+	return( (int)dsf_grid[0][*(const int *)p1].crad -
+		(int)dsf_grid[0][*(const int *)p2].crad );
+}
+
 /* Cull points for more uniform distribution, leave all nval 0 or 1 */
 static void
 cull_values(void)
 {
+	int	indx[GRIDRES*GRIDRES];
 	FVECT	ovec0, ovec1;
 	double	maxang, maxang2;
-	int	i, j, ii, jj, r;
+	int	i, j, k, ii, jj, r;
+						/* sort by radius first */
+	for (k = GRIDRES*GRIDRES; k--; )
+		indx[k] = k;
+	qsort(indx, GRIDRES*GRIDRES, sizeof(int), &radius_cmp);
 						/* simple greedy algorithm */
-	for (i = 0; i < GRIDRES; i++)
-	    for (j = 0; j < GRIDRES; j++) {
+	for (k = GRIDRES*GRIDRES; k--; ) {
+		i = indx[k]/GRIDRES;		/* from biggest radius down */
+		j = indx[k] - i*GRIDRES;
 		if (!dsf_grid[i][j].nval)
 			continue;
 		if (!dsf_grid[i][j].crad)
-			continue;		/* shouldn't happen */
+			break;		/* shouldn't happen */
 		ovec_from_pos(ovec0, i, j);
 		maxang = 2.*R2ANG(dsf_grid[i][j].crad);
-		if (maxang > ovec0[2])		/* clamp near horizon */
-			maxang = ovec0[2];
+						/* clamp near horizon */
+		if (maxang > output_orient*ovec0[2])
+			maxang = output_orient*ovec0[2];
 		r = maxang*(2.*GRIDRES/M_PI) + 1;
 		maxang2 = maxang*maxang;
 		for (ii = i-r; ii <= i+r; ii++) {
@@ -196,7 +264,7 @@ cull_values(void)
 			dsf_grid[ii][jj].nval = 0;
 		    }
 		}
-	    }
+	}
 						/* final averaging pass */
 	for (i = 0; i < GRIDRES; i++)
 	    for (j = 0; j < GRIDRES; j++)
@@ -206,7 +274,7 @@ cull_values(void)
 		}
 }
 
-/* Compute minimum BSDF from histogram and clear it */
+/* Compute minimum BSDF from histogram (does not clear) */
 static void
 comp_bsdf_min()
 {
@@ -225,7 +293,6 @@ comp_bsdf_min()
 	for (i = 0; cnt <= target; i++)
 		cnt += bsdf_hist[i];
 	bsdf_min = histval(i-1);
-	memset(bsdf_hist, 0, sizeof(bsdf_hist));
 }
 
 /* Find n nearest sub-sampled neighbors to the given grid position */