--- ray/src/cv/bsdfmesh.c	2014/08/21 13:44:05	2.31
+++ ray/src/cv/bsdfmesh.c	2019/04/23 14:30:36	2.40
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdfmesh.c,v 2.31 2014/08/21 13:44:05 greg Exp $";
+static const char RCSid[] = "$Id: bsdfmesh.c,v 2.40 2019/04/23 14:30:36 greg Exp $";
 #endif
 /*
  * Create BSDF advection mesh from radial basis functions.
@@ -7,7 +7,7 @@ static const char RCSid[] = "$Id: bsdfmesh.c,v 2.31 20
  *	G. Ward
  */
 
-#ifndef _WIN32
+#if !defined(_WIN32) && !defined(_WIN64)
 #include <unistd.h>
 #include <sys/wait.h>
 #include <sys/mman.h>
@@ -51,7 +51,7 @@ eval_DSFsurround(const RBFNODE *rbf, const FVECT outve
 		if (i) spinvector(tvec, tvec, outvec, phinc);
 		if (tvec[2] > 0 ^ output_orient > 0)
 			continue;
-		sum += eval_rbfrep(rbf, tvec) * output_orient*tvec[2];
+		sum += eval_rbfrep(rbf, tvec) * COSF(tvec[2]);
 		++n;
 	}
 	if (n < 2)				/* should never happen! */
@@ -63,56 +63,93 @@ eval_DSFsurround(const RBFNODE *rbf, const FVECT outve
 static double
 est_DSFrad(const RBFNODE *rbf, const FVECT outvec)
 {
-	const double	rad_epsilon = 0.03;
-	const double	DSFtarget = 0.60653066 * eval_rbfrep(rbf,outvec)
-						* output_orient*outvec[2];
+	const double	rad_epsilon = 0.01;
+	const double	DSFtarget = 0.60653066 * eval_rbfrep(rbf,outvec) *
+							COSF(outvec[2]);
 	double		inside_rad = rad_epsilon;
 	double		outside_rad = 0.5;
 	double		DSFinside = eval_DSFsurround(rbf, outvec, inside_rad);
 	double		DSFoutside = eval_DSFsurround(rbf, outvec, outside_rad);
 #define	interp_rad	inside_rad + (outside_rad-inside_rad) * \
 				(DSFtarget-DSFinside) / (DSFoutside-DSFinside)
-						/* interpolation search */
-	while (outside_rad-inside_rad > rad_epsilon) {
+						/* Newton's method (sort of) */
+	do {
 		double	test_rad = interp_rad;
-		double	DSFtest = eval_DSFsurround(rbf, outvec, test_rad);
-		if (DSFtarget < DSFtest) {
+		double	DSFtest;
+		if ((test_rad >= outside_rad) | (test_rad <= inside_rad))
+			test_rad = .5*(inside_rad + outside_rad);
+		DSFtest = eval_DSFsurround(rbf, outvec, test_rad);
+		if (DSFtest > DSFtarget) {
 			inside_rad = test_rad;
 			DSFinside = DSFtest;
 		} else {
 			outside_rad = test_rad;
 			DSFoutside = DSFtest;
 		}
-	}
-	return(interp_rad);
+	} while (outside_rad-inside_rad > rad_epsilon);
+
+	return(.5*(inside_rad + outside_rad));
 #undef interp_rad
 }
 
-/* Compute average BSDF peak from current DSF's */
+static int
+dbl_cmp(const void *p1, const void *p2)
+{
+	double	d1 = *(const double *)p1;
+	double	d2 = *(const double *)p2;
+
+	if (d1 > d2) return(1);
+	if (d1 < d2) return(-1);
+	return(0);
+}
+
+/* Conservative estimate of average BSDF value from current DSF's */
 static void
 comp_bsdf_spec(void)
 {
-	double		peak_sum = 0;
+	double		vmod_sum = 0;
 	double		rad_sum = 0;
 	int		n = 0;
+	double		*cost_list = NULL;
+	double		max_cost = 1.;
 	RBFNODE		*rbf;
 	FVECT		sdv;
-
-	if (dsf_list == NULL) {
-		bsdf_spec_peak = 0;
+						/* sort by incident altitude */
+	for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
+		n++;
+	if (n >= 10)
+		cost_list = (double *)malloc(sizeof(double)*n);
+	if (cost_list == NULL) {
+		bsdf_spec_val = 0;
 		bsdf_spec_rad = 0;
 		return;
 	}
+	n = 0;
+	for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
+		cost_list[n++] = rbf->invec[2]*input_orient;
+	qsort(cost_list, n, sizeof(double), dbl_cmp);
+	max_cost = cost_list[(n+3)/4];		/* accept 25% nearest grazing */
+	free(cost_list);
+	n = 0;
 	for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
+		double	this_rad, cosfact, vest;
+		if (rbf->invec[2]*input_orient > max_cost)
+			continue;
 		sdv[0] = -rbf->invec[0];
 		sdv[1] = -rbf->invec[1];
 		sdv[2] = rbf->invec[2]*(2*(input_orient==output_orient) - 1);
-		peak_sum += eval_rbfrep(rbf, sdv);
-		rad_sum += est_DSFrad(rbf, sdv);
+		cosfact = COSF(sdv[2]);
+		this_rad = est_DSFrad(rbf, sdv);
+		vest = eval_rbfrep(rbf, sdv) * cosfact *
+				(2.*M_PI) * this_rad*this_rad;
+		if (vest > rbf->vtotal)		/* don't over-estimate energy */
+			vest = rbf->vtotal;
+		vmod_sum += vest / cosfact;	/* remove cosine factor */
+		rad_sum += this_rad;
 		++n;
 	}
-	bsdf_spec_peak = peak_sum/(double)n;
 	bsdf_spec_rad = rad_sum/(double)n;
+	bsdf_spec_val = vmod_sum/(2.*M_PI*n*bsdf_spec_rad*bsdf_spec_rad);
 }
 
 /* Create a new migration holder (sharing memory for multiprocessing) */
@@ -122,7 +159,7 @@ new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
 	size_t		memlen = sizeof(MIGRATION) +
 				sizeof(float)*(from_rbf->nrbf*to_rbf->nrbf - 1);
 	MIGRATION	*newmig;
-#ifdef _WIN32
+#if defined(_WIN32) || defined(_WIN64)
 	if (nprocs > 1)
 		fprintf(stderr, "%s: warning - multiprocessing not supported\n",
 				progname);
@@ -153,7 +190,7 @@ new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
 	return(mig_list = newmig);
 }
 
-#ifdef _WIN32
+#if defined(_WIN32) || defined(_WIN64)
 #define await_children(n)	(void)(n)
 #define run_subprocess()	0
 #define end_subprocess()	(void)0
@@ -451,7 +488,7 @@ check_normal_incidence(void)
 		default:
 			return;			/* else we can interpolate */
 		}
-		for (rbf = near_rbf->next; rbf != NULL; rbf = rbf->next) {
+		for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
 			const double	d = input_orient*rbf->invec[2];
 			if (d >= 1.-2.*FTINY)
 				return;		/* seems we have normal */
@@ -501,9 +538,11 @@ memerr:
 void
 build_mesh(void)
 {
+	int		nrbfs = 0, nmigs = 0;
 	double		best2 = M_PI*M_PI;
 	RBFNODE		*shrt_edj[2];
 	RBFNODE		*rbf0, *rbf1;
+	const MIGRATION	*ej;
 						/* average specular peak */
 	comp_bsdf_spec();
 						/* add normal if needed */
@@ -517,7 +556,7 @@ build_mesh(void)
 		return;
 	}
 	shrt_edj[0] = shrt_edj[1] = NULL;	/* start w/ shortest edge */
-	for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
+	for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next) {
 	    for (rbf1 = rbf0->next; rbf1 != NULL; rbf1 = rbf1->next) {
 		double	dist2 = 2. - 2.*DOT(rbf0->invec,rbf1->invec);
 		if (dist2 < best2) {
@@ -525,6 +564,8 @@ build_mesh(void)
 			shrt_edj[1] = rbf1;
 			best2 = dist2;
 		}
+	    }
+	    ++nrbfs;
 	}
 	if (shrt_edj[0] == NULL) {
 		fprintf(stderr, "%s: Cannot find shortest edge\n", progname);
@@ -535,6 +576,13 @@ build_mesh(void)
 		mesh_from_edge(create_migration(shrt_edj[0], shrt_edj[1]));
 	else
 		mesh_from_edge(create_migration(shrt_edj[1], shrt_edj[0]));
+						/* count up edges */
+	for (ej = mig_list; ej != NULL; ej = ej->next)
+		++nmigs;
+	if (nmigs < nrbfs-1)			/* did meshing fail? */
+		fprintf(stderr,
+	    "%s: warning - %d incident directions but only %d interpolant(s)\n",
+				progname, nrbfs, nmigs);
 						/* complete migrations */
 	await_children(nchild);
 }