--- ray/src/cv/bsdfmesh.c	2013/11/08 03:42:13	2.14
+++ 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.14 2013/11/08 03:42:13 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.
@@ -484,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
@@ -492,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)