--- ray/src/cv/bsdfrbf.c	2012/11/13 04:23:38	2.2
+++ ray/src/cv/bsdfrbf.c	2013/09/25 17:42:45	2.7
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdfrbf.c,v 2.2 2012/11/13 04:23:38 greg Exp $";
+static const char RCSid[] = "$Id: bsdfrbf.c,v 2.7 2013/09/25 17:42:45 greg Exp $";
 #endif
 /*
  * Radial basis function representation for BSDF data.
@@ -17,6 +17,12 @@ static const char RCSid[] = "$Id: bsdfrbf.c,v 2.2 2012
 #ifndef RSCA
 #define RSCA		2.7		/* radius scaling factor (empirical) */
 #endif
+#ifndef MAXFRAC
+#define MAXFRAC		0.5		/* maximum contribution to neighbor */
+#endif
+#ifndef NNEIGH
+#define NNEIGH		10		/* number of neighbors to consider */
+#endif
 				/* our loaded grid for this incident angle */
 GRIDVAL			dsf_grid[GRIDRES][GRIDRES];
 
@@ -51,6 +57,10 @@ add_bsdf_data(double theta_out, double phi_out, double
 	if (!isDSF)
 		val *= ovec[2];		/* convert from BSDF to DSF */
 
+					/* update BSDF histogram */
+	if (val < BSDF2BIG*ovec[2] && val > BSDF2SML*ovec[2])
+		++bsdf_hist[histndx(val/ovec[2])];
+
 	pos_from_vec(pos, ovec);
 
 	dsf_grid[pos[0]][pos[1]].vsum += val;
@@ -185,6 +195,89 @@ cull_values(void)
 		}
 }
 
+/* Compute minimum BSDF from histogram and clear it */
+static void
+comp_bsdf_min()
+{
+	int	cnt;
+	int	i, target;
+
+	cnt = 0;
+	for (i = HISTLEN; i--; )
+		cnt += bsdf_hist[i];
+	if (!cnt) {				/* shouldn't happen */
+		bsdf_min = 0;
+		return;
+	}
+	target = cnt/100;			/* ignore bottom 1% */
+	cnt = 0;
+	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 */
+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	rad0 = R2ANG(dsf_grid[i][j].crad);
+	double		currad = RSCA * rad0;
+	int		neigh[NNEIGH][2];
+	int		n;
+	FVECT		our_dir;
+
+	ovec_from_pos(our_dir, i, j);
+	n = get_neighbors(neigh, NNEIGH, 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 = MAXFRAC * 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)
@@ -203,6 +296,8 @@ make_rbfrep(void)
 	for (i = 0; i < GRIDRES; i++)
 	    for (j = 0; j < GRIDRES; j++)
 		nn += dsf_grid[i][j].nval;
+				/* compute minimum BSDF */
+	comp_bsdf_min();
 				/* allocate RBF array */
 	newnode = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(nn-1));
 	if (newnode == NULL)
@@ -221,7 +316,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;
@@ -260,7 +355,11 @@ make_rbfrep(void)
 	nn = 0;			/* compute sum for normalization */
 	while (nn < newnode->nrbf)
 		newnode->vtotal += rbf_volume(&newnode->rbfa[nn++]);
-
+#ifdef DEBUG
+	fprintf(stderr, "Integrated DSF at (%.1f,%.1f) deg. is %.2f\n",
+			get_theta180(newnode->invec), get_phi360(newnode->invec),
+			newnode->vtotal);
+#endif
 	insert_dsf(newnode);
 
 	return(newnode);