--- ray/src/cv/bsdfinterp.c	2012/10/23 21:09:29	2.4
+++ ray/src/cv/bsdfinterp.c	2012/12/12 04:49:59	2.8
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdfinterp.c,v 2.4 2012/10/23 21:09:29 greg Exp $";
+static const char RCSid[] = "$Id: bsdfinterp.c,v 2.8 2012/12/12 04:49:59 greg Exp $";
 #endif
 /*
  * Interpolate BSDF data from radial basis functions in advection mesh.
@@ -126,7 +126,7 @@ in_tri(const RBFNODE *v1, const RBFNODE *v2, const RBF
 	return(sgn2 == sgn3);
 }
 
-/* Test and set for edge */
+/* Test (and set) bitmap for edge */
 static int
 check_edge(unsigned char *emap, int nedges, const MIGRATION *mig, int mark)
 {
@@ -202,8 +202,15 @@ get_interp(MIGRATION *miga[3], FVECT invec)
 							input_orient*invec[2]) {
 				for (miga[0] = rbf->ejl; miga[0] != NULL;
 						miga[0] = nextedge(rbf,miga[0]))
-					if (opp_rbf(rbf,miga[0]) == rbf->next)
+					if (opp_rbf(rbf,miga[0]) == rbf->next) {
+						double	nf = 1.-rbf->invec[2]*rbf->invec[2];
+						if (nf > FTINY) {
+							nf = sqrt((1.-invec[2]*invec[2])/nf);
+							invec[0] = nf*rbf->invec[0];
+							invec[1] = nf*rbf->invec[1];
+						}
 						return(0);
+					}
 				break;
 			}
 		}
@@ -247,21 +254,23 @@ e_advect_rbf(const MIGRATION *mig, const FVECT invec)
 	double		t, full_dist;
 						/* get relative position */
 	t = acos(DOT(invec, mig->rbfv[0]->invec));
-	if (t < M_PI/GRIDRES) {			/* near first DSF */
+	if (t < M_PI/grid_res) {		/* near first DSF */
 		n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[0]->nrbf-1);
 		rbf = (RBFNODE *)malloc(n);
 		if (rbf == NULL)
 			goto memerr;
 		memcpy(rbf, mig->rbfv[0], n);	/* just duplicate */
+		rbf->next = NULL; rbf->ejl = NULL;
 		return(rbf);
 	}
 	full_dist = acos(DOT(mig->rbfv[0]->invec, mig->rbfv[1]->invec));
-	if (t > full_dist-M_PI/GRIDRES) {	/* near second DSF */
+	if (t > full_dist-M_PI/grid_res) {	/* near second DSF */
 		n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[1]->nrbf-1);
 		rbf = (RBFNODE *)malloc(n);
 		if (rbf == NULL)
 			goto memerr;
 		memcpy(rbf, mig->rbfv[1], n);	/* just duplicate */
+		rbf->next = NULL; rbf->ejl = NULL;
 		return(rbf);
 	}
 	t /= full_dist;	
@@ -332,7 +341,10 @@ advect_rbf(const FVECT invec)
 		return(NULL);
 	if (miga[1] == NULL) {			/* advect along edge? */
 		rbf = e_advect_rbf(miga[0], sivec);
-		rev_rbf_symmetry(rbf, sym);
+		if (single_plane_incident)
+			rotate_rbf(rbf, invec);
+		else
+			rev_rbf_symmetry(rbf, sym);
 		return(rbf);
 	}
 #ifdef DEBUG
@@ -400,7 +412,6 @@ advect_rbf(const FVECT invec)
 		    float		mc = mtx_coef(miga[2],i,k);
 		    const RBFVAL	*rbf2k;
 		    double		rad2k;
-		    FVECT		vout;
 		    int			pos[2];
 		    if ((mb <= FTINY) & (mc <= FTINY))
 			continue;
@@ -409,8 +420,8 @@ advect_rbf(const FVECT invec)
 		    rad2k = R2ANG(rbf2k->crad);
 		    rbf->rbfa[n].crad = ANG2R(sqrt(srad2 + t*rad2k*rad2k));
 		    ovec_from_pos(v2, rbf2k->gx, rbf2k->gy);
-		    geodesic(vout, v1, v2, t, GEOD_REL);
-		    pos_from_vec(pos, vout);
+		    geodesic(v2, v1, v2, t, GEOD_REL);
+		    pos_from_vec(pos, v2);
 		    rbf->rbfa[n].gx = pos[0];
 		    rbf->rbfa[n].gy = pos[1];
 		    ++n;