--- ray/src/cv/bsdfrbf.c	2013/06/28 23:18:51	2.5
+++ ray/src/cv/bsdfrbf.c	2013/09/25 05:03:10	2.6
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdfrbf.c,v 2.5 2013/06/28 23:18:51 greg Exp $";
+static const char RCSid[] = "$Id: bsdfrbf.c,v 2.6 2013/09/25 05:03:10 greg Exp $";
 #endif
 /*
  * Radial basis function representation for BSDF data.
@@ -211,6 +211,68 @@ comp_bsdf_min()
 	memset(bsdf_hist, 0, sizeof(bsdf_hist));
 }
 
+/* Find n nearest sub-sampled neighbors to the given grid position */
+static int
+get_neighbors(int neigh[][2], int n, const int i, const int j)
+{
+	int	k = 0;
+	int	r;
+						/* search concentric squares */
+	for (r = 1; r < GRIDRES; r++) {
+		int	ii, jj;
+		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) {
+				if (jj < 0) continue;
+				if (jj >= GRIDRES) break;
+				if (dsf_grid[ii][jj].nval) {
+					neigh[k][0] = ii;
+					neigh[k][1] = jj;
+					if (++k >= n)
+						return(n);
+				}
+			}
+		}
+	}
+	return(k);
+}
+
+/* Adjust coded radius for the given grid position based on neighborhood */
+static int
+adj_coded_radius(const int i, const int j)
+{
+	const double	max_frac = 0.33;
+	const double	rad0 = R2ANG(dsf_grid[i][j].crad);
+	double		currad = RSCA * rad0;
+	int		neigh[5][2];
+	int		n;
+	FVECT		our_dir;
+
+	ovec_from_pos(our_dir, i, j);
+	n = get_neighbors(neigh, 5, i, j);
+	while (n--) {
+		FVECT	their_dir;
+		double	max_ratio, rad_ok2;
+						/* check our value at neighbor */
+		ovec_from_pos(their_dir, neigh[n][0], neigh[n][1]);
+		max_ratio = max_frac * dsf_grid[neigh[n][0]][neigh[n][1]].vsum
+				/ dsf_grid[i][j].vsum;
+		if (max_ratio >= 1)
+			continue;
+		rad_ok2 = (DOT(their_dir,our_dir) - 1.)/log(max_ratio);
+		if (rad_ok2 >= currad*currad)
+			continue;		/* value fraction OK */
+		currad = sqrt(rad_ok2);		/* else reduce lobe radius */
+		if (currad <= rad0)		/* limit how small we'll go */
+			return(dsf_grid[i][j].crad);
+	}
+	return(ANG2R(currad));			/* encode selected radius */
+}
+
 /* Count up filled nodes and build RBF representation from current grid */ 
 RBFNODE *
 make_rbfrep(void)
@@ -249,7 +311,7 @@ make_rbfrep(void)
 	    for (j = 0; j < GRIDRES; j++)
 		if (dsf_grid[i][j].nval) {
 			newnode->rbfa[nn].peak = dsf_grid[i][j].vsum;
-			newnode->rbfa[nn].crad = RSCA*dsf_grid[i][j].crad + .5;
+			newnode->rbfa[nn].crad = adj_coded_radius(i, j);
 			newnode->rbfa[nn].gx = i;
 			newnode->rbfa[nn].gy = j;
 			++nn;