--- ray/src/cv/bsdfrbf.c	2012/10/19 04:14:29	2.1
+++ ray/src/cv/bsdfrbf.c	2013/06/28 23:18:51	2.5
@@ -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.5 2013/06/28 23:18:51 greg Exp $";
 #endif
 /*
  * Radial basis function representation for BSDF data.
@@ -51,6 +51,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 +189,28 @@ 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));
+}
+
 /* Count up filled nodes and build RBF representation from current grid */ 
 RBFNODE *
 make_rbfrep(void)
@@ -193,6 +219,7 @@ make_rbfrep(void)
 	double	lastVar, thisVar = 100.;
 	int	nn;
 	RBFNODE	*newnode;
+	RBFVAL	*itera;
 	int	i, j;
 				/* compute RBF radii */
 	compute_radii();
@@ -202,12 +229,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;
@@ -228,6 +255,10 @@ make_rbfrep(void)
 			++nn;
 		}
 				/* 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 +268,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 +284,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);
 }