--- ray/src/cv/bsdfrbf.c	2013/06/28 23:18:51	2.5
+++ 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.5 2013/06/28 23:18:51 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.
@@ -14,9 +14,24 @@ static const char RCSid[] = "$Id: bsdfrbf.c,v 2.5 2013
 #include <math.h>
 #include "bsdfrep.h"
 
-#ifndef RSCA
-#define RSCA		2.7		/* radius scaling factor (empirical) */
+#ifndef MINRSCA
+#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
+#ifndef NNEIGH
+#define NNEIGH		10		/* number of neighbors to consider */
+#endif
 				/* our loaded grid for this incident angle */
 GRIDVAL			dsf_grid[GRIDRES][GRIDRES];
 
@@ -48,7 +63,9 @@ add_bsdf_data(double theta_out, double phi_out, double
 	ovec[1] = sin((M_PI/180.)*phi_out) * ovec[2];
 	ovec[2] = sqrt(1. - ovec[2]*ovec[2]);
 
-	if (!isDSF)
+	if (val <= 0)			/* truncate to zero */
+		val = 0;
+	else if (!isDSF)
 		val *= ovec[2];		/* convert from BSDF to DSF */
 
 					/* update BSDF histogram */
@@ -62,10 +79,11 @@ add_bsdf_data(double theta_out, double phi_out, double
 }
 
 /* Compute radii for non-empty bins */
-/* (distance to furthest empty bin for which non-empty bin is the closest) */
+/* (distance to furthest empty bin for which non-empty test bin is closest) */
 static void
 compute_radii(void)
 {
+	const int	cradmin = ANG2R(.5*M_PI/GRIDRES);
 	unsigned int	fill_grid[GRIDRES][GRIDRES];
 	unsigned short	fill_cnt[GRIDRES][GRIDRES];
 	FVECT		ovec0, ovec1;
@@ -110,14 +128,66 @@ 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));
 	for (i = 0; i < GRIDRES; i++)
 	    for (j = 0; j < GRIDRES; j++) {
-		if (!dsf_grid[i][j].crad)
-			continue;		/* missing distance */
-		r = R2ANG(dsf_grid[i][j].crad)*(2.*RSCA*GRIDRES/M_PI);
+		if (!dsf_grid[i][j].nval)
+			continue;		/* not part of this */
+		r = R2ANG(dsf_grid[i][j].crad)*(2.*MAXRSCA*GRIDRES/M_PI);
 		for (ii = i-r; ii <= i+r; ii++) {
 		    if (ii < 0) continue;
 		    if (ii >= GRIDRES) break;
@@ -138,36 +208,51 @@ 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++) {
 		    if (ii < 0) continue;
 		    if (ii >= GRIDRES) break;
 		    for (jj = j-r; jj <= j+r; jj++) {
+			if ((ii == i) & (jj == j))
+				continue;	/* don't get self-absorbed */
 			if (jj < 0) continue;
 			if (jj >= GRIDRES) break;
 			if (!dsf_grid[ii][jj].nval)
 				continue;
-			if ((ii == i) & (jj == j))
-				continue;	/* don't get self-absorbed */
 			ovec_from_pos(ovec1, ii, jj);
 			if (2. - 2.*DOT(ovec0,ovec1) >= maxang2)
 				continue;
@@ -179,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++)
@@ -189,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()
 {
@@ -208,19 +293,100 @@ 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 */
+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);
+	const double	minrad = MINRSCA * rad0;
+	double		currad = MAXRSCA * 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 <= minrad)		/* limit how small we'll go */
+			return(ANG2R(minrad));
+	}
+	return(ANG2R(currad));			/* encode selected radius */
+}
+
 /* Count up filled nodes and build RBF representation from current grid */ 
 RBFNODE *
 make_rbfrep(void)
 {
+	long	cradsum = 0, ocradsum = 0;
 	int	niter = 16;
 	double	lastVar, thisVar = 100.;
 	int	nn;
 	RBFNODE	*newnode;
 	RBFVAL	*itera;
 	int	i, j;
+
+#ifdef DEBUG
+{
+	int	maxcnt = 0, nempty = 0;
+	long	cntsum = 0;
+	for (i = 0; i < GRIDRES; i++)
+	    for (j = 0; j < GRIDRES; j++)
+		if (!dsf_grid[i][j].nval) {
+			++nempty;
+		} else {
+			if (dsf_grid[i][j].nval > maxcnt)
+				maxcnt = dsf_grid[i][j].nval;
+			cntsum += dsf_grid[i][j].nval;
+		}
+	fprintf(stderr, "Average, maximum bin count: %d, %d (%.1f%% empty)\n",
+			(int)(cntsum/((GRIDRES*GRIDRES)-nempty)), maxcnt,
+			100./(GRIDRES*GRIDRES)*nempty);
+}
+#endif
 				/* compute RBF radii */
 	compute_radii();
 				/* coagulate lobes */
@@ -249,11 +415,19 @@ 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;
+			ocradsum += dsf_grid[i][j].crad;
+			cradsum +=
+			newnode->rbfa[nn].crad = adj_coded_radius(i, j);
 			newnode->rbfa[nn].gx = i;
 			newnode->rbfa[nn].gy = j;
 			++nn;
 		}
+#ifdef DEBUG
+	fprintf(stderr,
+	"Average radius reduced from %.2f to %.2f degrees for %d lobes\n",
+			180./M_PI*MAXRSCA*R2ANG(ocradsum/newnode->nrbf),
+			180./M_PI*R2ANG(cradsum/newnode->nrbf), newnode->nrbf);
+#endif
 				/* iterate to improve interpolation accuracy */
 	itera = (RBFVAL *)malloc(sizeof(RBFVAL)*newnode->nrbf);
 	if (itera == NULL)