--- ray/src/cv/bsdfrbf.c	2012/10/19 04:14:29	2.1
+++ ray/src/cv/bsdfrbf.c	2013/10/17 19:09:11	2.9
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdfrbf.c,v 2.1 2012/10/19 04:14:29 greg Exp $";
+static const char RCSid[] = "$Id: bsdfrbf.c,v 2.9 2013/10/17 19:09:11 greg Exp $";
 #endif
 /*
  * Radial basis function representation for BSDF data.
@@ -14,9 +14,18 @@ static const char RCSid[] = "$Id: bsdfrbf.c,v 2.1 2012
 #include <math.h>
 #include "bsdfrep.h"
 
-#ifndef RSCA
-#define RSCA		2.7		/* radius scaling factor (empirical) */
+#ifndef MINRSCA
+#define MINRSCA		0.15		/* minimum radius scaling factor */
 #endif
+#ifndef MAXRSCA
+#define MAXRSCA		2.7		/* maximum radius scaling factor */
+#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,9 +57,15 @@ 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 */
+	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;
@@ -58,16 +73,22 @@ 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;
 	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++) {
@@ -111,9 +132,9 @@ compute_radii(void)
 	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;
@@ -158,12 +179,12 @@ cull_values(void)
 		    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;
@@ -185,15 +206,120 @@ 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);
+	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 */
@@ -202,12 +328,12 @@ 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) {
-		fprintf(stderr, "%s: Out of memory in make_rbfrep()\n", progname);
-		exit(1);
-	}
+	if (newnode == NULL)
+		goto memerr;
 	newnode->ord = -1;
 	newnode->next = NULL;
 	newnode->ejl = NULL;
@@ -222,12 +348,24 @@ 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)
+		goto memerr;
+	memcpy(itera, newnode->rbfa, sizeof(RBFVAL)*newnode->nrbf);
 	do {
 		double	dsum = 0, dsum2 = 0;
 		nn = 0;
@@ -237,12 +375,13 @@ make_rbfrep(void)
 				FVECT	odir;
 				double	corr;
 				ovec_from_pos(odir, i, j);
-				newnode->rbfa[nn++].peak *= corr =
+				itera[nn++].peak *= corr =
 					dsf_grid[i][j].vsum /
 						eval_rbfrep(newnode, odir);
-				dsum += corr - 1.;
-				dsum2 += (corr-1.)*(corr-1.);
+				dsum += 1. - corr;
+				dsum2 += (1.-corr)*(1.-corr);
 			}
+		memcpy(newnode->rbfa, itera, sizeof(RBFVAL)*newnode->nrbf);
 		lastVar = thisVar;
 		thisVar = dsum2/(double)nn;
 #ifdef DEBUG
@@ -252,11 +391,19 @@ make_rbfrep(void)
 #endif
 	} while (--niter > 0 && lastVar-thisVar > 0.02*lastVar);
 
+	free(itera);
 	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);
+memerr:
+	fprintf(stderr, "%s: Out of memory in make_rbfrep()\n", progname);
+	exit(1);
 }