--- ray/src/cv/bsdfmesh.c	2013/10/24 16:11:37	2.11
+++ ray/src/cv/bsdfmesh.c	2014/02/19 05:16:06	2.16
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdfmesh.c,v 2.11 2013/10/24 16:11:37 greg Exp $";
+static const char RCSid[] = "$Id: bsdfmesh.c,v 2.16 2014/02/19 05:16:06 greg Exp $";
 #endif
 /*
  * Create BSDF advection mesh from radial basis functions.
@@ -176,9 +176,12 @@ price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, co
 	    pm->prow = pricerow(pm,i);
 	    srow = psortrow(pm,i);
 	    for (j = to_rbf->nrbf; j--; ) {
-		double		dprod = DOT(vfrom, vto[j]);
-		pm->prow[j] = ((dprod >= 1.) ? .0 : acos(dprod)) +
-				fabs(R2ANG(to_rbf->rbfa[j].crad) - from_ang);
+		double		d;		/* quadratic cost function */
+		d = DOT(vfrom, vto[j]);
+		d = (d >= 1.) ? .0 : acos(d);
+		pm->prow[j] = d*d;
+		d = R2ANG(to_rbf->rbfa[j].crad) - from_ang;
+		pm->prow[j] += d*d;	
 		srow[j] = j;
 	    }
 	    qsort_r(srow, pm->ncols, sizeof(short), pm, &msrt_cmp);
@@ -236,7 +239,7 @@ migration_step(MIGRATION *mig, double *src_rem, double
 {
 	const int	max2check = 100;
 	const double	maxamt = 1./(double)pm->ncols;
-	const double	minamt = maxamt*5e-6;
+	const double	minamt = maxamt*1e-4;
 	double		*src_cost;
 	short		(*rord)[2];
 	struct {
@@ -333,10 +336,10 @@ create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
 			return(NULL);
 						/* else allocate */
 #ifdef DEBUG
-	fprintf(stderr, "Building path from (theta,phi) (%.0f,%.0f) ",
+	fprintf(stderr, "Building path from (theta,phi) (%.1f,%.1f) ",
 			get_theta180(from_rbf->invec),
 			get_phi360(from_rbf->invec));
-	fprintf(stderr, "to (%.0f,%.0f) with %d x %d matrix\n",
+	fprintf(stderr, "to (%.1f,%.1f) with %d x %d matrix\n",
 			get_theta180(to_rbf->invec),
 			get_phi360(to_rbf->invec), 
 			from_rbf->nrbf, to_rbf->nrbf);
@@ -404,7 +407,7 @@ overlaps_tri(const RBFNODE *bv0, const RBFNODE *bv1, c
 	return(vother[im_rev] != NULL);
 }
 
-/* Find context hull vertex to complete triangle (oriented call) */
+/* Find convex hull vertex to complete triangle (oriented call) */
 static RBFNODE *
 find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rbf1)
 {
@@ -425,7 +428,7 @@ find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rb
 		if (DOT(vp, vmid) <= FTINY)
 			continue;		/* wrong orientation */
 		area2 = .25*DOT(vp,vp);
-		VSUB(vp, rbf->invec, rbf0->invec);
+		VSUB(vp, rbf->invec, vmid);
 		dprod = -DOT(vp, vejn);
 		VSUM(vp, vp, vejn, dprod);	/* above guarantees non-zero */
 		dprod = DOT(vp, vmid) / VLEN(vp);
@@ -481,6 +484,78 @@ mesh_from_edge(MIGRATION *edge)
 		}
 	}
 }
+
+/* Add normal direction if missing */
+static void
+check_normal_incidence(void)
+{
+	static const FVECT	norm_vec = {.0, .0, 1.};
+	const int		saved_nprocs = nprocs;
+	RBFNODE			*near_rbf, *mir_rbf, *rbf;
+	double			bestd;
+	int			n;
+
+	if (dsf_list == NULL)
+		return;				/* XXX should be error? */
+	near_rbf = dsf_list;
+	bestd = input_orient*near_rbf->invec[2];
+	if (single_plane_incident) {		/* ordered plane incidence? */
+		if (bestd >= 1.-2.*FTINY)
+			return;			/* already have normal */
+	} else {
+		switch (inp_coverage) {
+		case INP_QUAD1:
+		case INP_QUAD2:
+		case INP_QUAD3:
+		case INP_QUAD4:
+			break;			/* quadrilateral symmetry? */
+		default:
+			return;			/* else we can interpolate */
+		}
+		for (rbf = near_rbf->next; rbf != NULL; rbf = rbf->next) {
+			const double	d = input_orient*rbf->invec[2];
+			if (d >= 1.-2.*FTINY)
+				return;		/* seems we have normal */
+			if (d > bestd) {
+				near_rbf = rbf;
+				bestd = d;
+			}
+		}
+	}
+	if (mig_list != NULL) {			/* need to be called first */
+		fprintf(stderr, "%s: Late call to check_normal_incidence()\n",
+				progname);
+		exit(1);
+	}
+#ifdef DEBUG
+	fprintf(stderr, "Interpolating normal incidence by mirroring (%.1f,%.1f)\n",
+			get_theta180(near_rbf->invec), get_phi360(near_rbf->invec));
+#endif
+						/* mirror nearest incidence */
+	n = sizeof(RBFNODE) + sizeof(RBFVAL)*(near_rbf->nrbf-1);
+	mir_rbf = (RBFNODE *)malloc(n);
+	if (mir_rbf == NULL)
+		goto memerr;
+	memcpy(mir_rbf, near_rbf, n);
+	mir_rbf->ord = near_rbf->ord - 1;	/* not used, I think */
+	mir_rbf->next = NULL;
+	rev_rbf_symmetry(mir_rbf, MIRROR_X|MIRROR_Y);
+	nprocs = 1;				/* compute migration matrix */
+	if (mig_list != create_migration(mir_rbf, near_rbf))
+		exit(1);			/* XXX should never happen! */
+						/* interpolate normal dist. */
+	rbf = e_advect_rbf(mig_list, norm_vec, 2*near_rbf->nrbf);
+	nprocs = saved_nprocs;			/* final clean-up */
+	free(mir_rbf);
+	free(mig_list);
+	mig_list = near_rbf->ejl = NULL;
+	insert_dsf(rbf);			/* insert interpolated normal */
+	return;
+memerr:
+	fprintf(stderr, "%s: Out of memory in check_normal_incidence()\n",
+				progname);
+	exit(1);
+}
 	
 /* Build our triangle mesh from recorded RBFs */
 void
@@ -489,6 +564,8 @@ build_mesh(void)
 	double		best2 = M_PI*M_PI;
 	RBFNODE		*shrt_edj[2];
 	RBFNODE		*rbf0, *rbf1;
+						/* add normal if needed */
+	check_normal_incidence();
 						/* check if isotropic */
 	if (single_plane_incident) {
 		for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)