--- ray/src/cv/bsdfrep.c	2012/11/22 06:07:17	2.11
+++ ray/src/cv/bsdfrep.c	2013/10/22 05:48:05	2.16
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdfrep.c,v 2.11 2012/11/22 06:07:17 greg Exp $";
+static const char RCSid[] = "$Id: bsdfrep.c,v 2.16 2013/10/22 05:48:05 greg Exp $";
 #endif
 /*
  * Support BSDF representation as radial basis functions.
@@ -26,6 +26,12 @@ int			single_plane_incident = -1;
 int			input_orient = 0;
 int			output_orient = 0;
 
+				/* BSDF histogram */
+unsigned long		bsdf_hist[HISTLEN];
+
+				/* BSDF value for boundary regions */
+double			bsdf_min = 0;
+
 				/* processed incident DSF measurements */
 RBFNODE			*dsf_list = NULL;
 
@@ -199,15 +205,6 @@ rotate_rbf(RBFNODE *rbf, const FVECT invec)
 	VCOPY(rbf->invec, invec);
 }
 
-/* Compute volume associated with Gaussian lobe */
-double
-rbf_volume(const RBFVAL *rbfp)
-{
-	double	rad = R2ANG(rbfp->crad);
-
-	return((2.*M_PI) * rbfp->peak * rad*rad);
-}
-
 /* Compute outgoing vector from grid position */
 void
 ovec_from_pos(FVECT vec, int xpos, int ypos)
@@ -237,26 +234,59 @@ pos_from_vec(int pos[2], const FVECT vec)
 	pos[1] = (int)(sq[1]*grid_res);
 }
 
+/* Compute volume associated with Gaussian lobe */
+double
+rbf_volume(const RBFVAL *rbfp)
+{
+	double	rad = R2ANG(rbfp->crad);
+	FVECT	odir;
+	double	elev, integ;
+				/* infinite integral approximation */
+	integ = (2.*M_PI) * rbfp->peak * rad*rad;
+				/* check if we're near horizon */
+	ovec_from_pos(odir, rbfp->gx, rbfp->gy);
+	elev = output_orient*odir[2];
+				/* apply cut-off correction if > 1% */
+	if (elev < 2.8*rad) {
+		/* elev = asin(elev);	/* this is so crude, anyway... */
+		integ *= 1. - .5*exp(-.5*elev*elev/(rad*rad));
+	}
+	return(integ);
+}
+
 /* Evaluate RBF for DSF at the given normalized outgoing direction */
 double
 eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
 {
+	double		minval = bsdf_min*output_orient*outvec[2];
+	int		pos[2];
 	double		res = 0;
 	const RBFVAL	*rbfp;
 	FVECT		odir;
-	double		sig2;
+	double		rad2;
 	int		n;
-
-	if (rp == NULL)
+				/* check for wrong side */
+	if (outvec[2] > 0 ^ output_orient > 0)
 		return(.0);
+				/* use minimum if no information avail. */
+	if (rp == NULL)
+		return(minval);
+				/* optimization for fast lobe culling */
+	pos_from_vec(pos, outvec);
+				/* sum radial basis function */
 	rbfp = rp->rbfa;
 	for (n = rp->nrbf; n--; rbfp++) {
+		int	d2 = (pos[0]-rbfp->gx)*(pos[0]-rbfp->gx) +
+				(pos[1]-rbfp->gy)*(pos[1]-rbfp->gy);
+		rad2 = R2ANG(rbfp->crad);
+		rad2 *= rad2;
+		if (d2 > (38.*GRIDRES*GRIDRES/M_PI/M_PI)*rad2)
+			continue;
 		ovec_from_pos(odir, rbfp->gx, rbfp->gy);
-		sig2 = R2ANG(rbfp->crad);
-		sig2 = (DOT(odir,outvec) - 1.) / (sig2*sig2);
-		if (sig2 > -19.)
-			res += rbfp->peak * exp(sig2);
+		res += rbfp->peak * exp((DOT(odir,outvec) - 1.) / rad2);
 	}
+	if (res < minval)	/* never return less than minval */
+		return(minval);
 	return(res);
 }
 
@@ -390,6 +420,7 @@ save_bsdf_rep(FILE *ofp)
 	fprintf(ofp, "SYMMETRY=%d\n", !single_plane_incident * inp_coverage);
 	fprintf(ofp, "IO_SIDES= %d %d\n", input_orient, output_orient);
 	fprintf(ofp, "GRIDRES=%d\n", grid_res);
+	fprintf(ofp, "BSDFMIN=%g\n", bsdf_min);
 	fputformat(BSDFREP_FMT, ofp);
 	fputc('\n', ofp);
 					/* write each DSF */
@@ -453,6 +484,10 @@ headline(char *s, void *p)
 	}
 	if (!strncmp(s, "GRIDRES=", 8)) {
 		sscanf(s+8, "%d", &grid_res);
+		return(0);
+	}
+	if (!strncmp(s, "BSDFMIN=", 8)) {
+		sscanf(s+8, "%lf", &bsdf_min);
 		return(0);
 	}
 	if (formatval(fmt, s) && strcmp(fmt, BSDFREP_FMT))