--- ray/src/cv/bsdfrep.c	2012/11/18 03:56:39	2.10
+++ ray/src/cv/bsdfrep.c	2013/09/20 00:21:39	2.14
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdfrep.c,v 2.10 2012/11/18 03:56:39 greg Exp $";
+static const char RCSid[] = "$Id: bsdfrep.c,v 2.14 2013/09/20 00:21:39 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 */
+int			bsdf_hist[HISTLEN];
+
+				/* BSDF value for boundary regions */
+double			bsdf_min = 0;
+
 				/* processed incident DSF measurements */
 RBFNODE			*dsf_list = NULL;
 
@@ -82,7 +88,7 @@ new_input_direction(double new_theta, double new_phi)
 int
 use_symmetry(FVECT vec)
 {
-	double	phi = get_phi360(vec);
+	const double	phi = get_phi360(vec);
 
 	switch (inp_coverage) {
 	case INP_QUAD1|INP_QUAD2|INP_QUAD3|INP_QUAD4:
@@ -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,18 +234,43 @@ 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];
 	double		res = 0;
 	const RBFVAL	*rbfp;
 	FVECT		odir;
 	double		sig2;
 	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);
+				/* sum radial basis function */
 	rbfp = rp->rbfa;
 	for (n = rp->nrbf; n--; rbfp++) {
 		ovec_from_pos(odir, rbfp->gx, rbfp->gy);
@@ -257,6 +279,8 @@ eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
 		if (sig2 > -19.)
 			res += rbfp->peak * exp(sig2);
 	}
+	if (res < minval)	/* never return less than minval */
+		return(minval);
 	return(res);
 }
 
@@ -390,6 +414,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 +478,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))