--- ray/src/cv/bsdfmesh.c	2014/03/05 22:47:16	2.17
+++ ray/src/cv/bsdfmesh.c	2014/03/26 00:11:30	2.26
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdfmesh.c,v 2.17 2014/03/05 22:47:16 greg Exp $";
+static const char RCSid[] = "$Id: bsdfmesh.c,v 2.26 2014/03/26 00:11:30 greg Exp $";
 #endif
 /*
  * Create BSDF advection mesh from radial basis functions.
@@ -18,6 +18,10 @@ static const char RCSid[] = "$Id: bsdfmesh.c,v 2.17 20
 #include <string.h>
 #include <math.h>
 #include "bsdfrep.h"
+
+#ifndef NEIGH_FACT2
+#define NEIGH_FACT2	0.1	/* empirical neighborhood distance weight */
+#endif
 				/* number of processes to run */
 int			nprocs = 1;
 				/* number of children (-1 in child) */
@@ -135,6 +139,47 @@ run_subprocess(void)
 
 #endif	/* ! _WIN32 */
 
+/* Compute normalized distribution scattering functions for comparison */
+static void
+compute_nDSFs(const RBFNODE *rbf0, const RBFNODE *rbf1)
+{
+	const double	nf0 = (GRIDRES*GRIDRES) / rbf0->vtotal;
+	const double	nf1 = (GRIDRES*GRIDRES) / rbf1->vtotal;
+	int		x, y;
+	FVECT		dv;
+
+	for (x = GRIDRES; x--; )
+	    for (y = GRIDRES; y--; ) {
+		ovec_from_pos(dv, x, y);	/* cube root (brightness) */
+		dsf_grid[x][y].val[0] = pow(nf0*eval_rbfrep(rbf0, dv), .3333);
+		dsf_grid[x][y].val[1] = pow(nf1*eval_rbfrep(rbf1, dv), .3333);
+	    }
+}	
+
+/* Compute neighborhood distance-squared (dissimilarity) */
+static double
+neighborhood_dist2(int x0, int y0, int x1, int y1)
+{
+	int	rad = GRIDRES>>5;
+	double	sum2 = 0.;
+	double	d;
+	int	p[4];
+	int	i, j;
+						/* check radius */
+	p[0] = x0; p[1] = y0; p[2] = x1; p[3] = y1;
+	for (i = 4; i--; ) {
+		if (p[i] < rad) rad = p[i];
+		if (GRIDRES-1-p[i] < rad) rad = GRIDRES-1-p[i];
+	}
+	for (i = -rad; i <= rad; i++)
+	    for (j = -rad; j <= rad; j++) {
+	    	d = dsf_grid[x0+i][y0+j].val[0] -
+			dsf_grid[x1+i][y1+j].val[1];
+		sum2 += d*d;
+	    }
+	return(sum2 / (4*rad*(rad+1) + 1));
+}
+
 /* Comparison routine needed for sorting price row */
 static int
 msrt_cmp(void *b, const void *p1, const void *p2)
@@ -155,6 +200,7 @@ price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, co
 	FVECT	*vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
 	int	i, j;
 
+	compute_nDSFs(from_rbf, to_rbf);
 	pm->nrows = from_rbf->nrbf;
 	pm->ncols = to_rbf->nrbf;
 	pm->price = (float *)malloc(sizeof(float) * pm->nrows*pm->ncols);
@@ -176,12 +222,15 @@ 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		d;		/* quadratic cost function */
-		d = DOT(vfrom, vto[j]);
-		d = (d >= 1.) ? .0 : acos(d);
+		double	d;			/* quadratic cost function */
+		d = Acos(DOT(vfrom, vto[j]));
 		pm->prow[j] = d*d;
 		d = R2ANG(to_rbf->rbfa[j].crad) - from_ang;
-		pm->prow[j] += d*d;	
+		pm->prow[j] += d*d;
+						/* neighborhood difference */
+		pm->prow[j] += NEIGH_FACT2 * neighborhood_dist2(
+				from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy,
+				to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy );
 		srow[j] = j;
 	    }
 	    qsort_r(srow, pm->ncols, sizeof(short), pm, &msrt_cmp);
@@ -334,6 +383,42 @@ memerr:
 #undef discount
 }
 
+#ifdef DUMP_MATRIX
+/* Dump transport plan and corresponding price matrix to a text file */
+static void
+dump_matrix(const MIGRATION *me, const PRICEMAT *pm)
+{
+	char	fname[256];
+	FILE	*fp;
+	int	i, j;
+
+	sprintf(fname, "edge_%d-%d.txt", me->rbfv[0]->ord, me->rbfv[1]->ord);
+	if ((fp = fopen(fname, "w")) == NULL)
+		return;
+	for (j = 0; j < 2; j++) {
+		fprintf(fp, "Available from %d source RBF lobes in node %d:\n",
+					me->rbfv[j]->nrbf, me->rbfv[j]->ord);
+		for (i = 0; i < me->rbfv[j]->nrbf; i++)
+			fprintf(fp, " %.4e", rbf_volume(&me->rbfv[j]->rbfa[i]) /
+							me->rbfv[j]->vtotal);
+		fputc('\n', fp);
+	}
+	fprintf(fp, "Price (quadratic distance metric) matrix:\n");
+	for (i = 0; i < pm->nrows; i++) {
+	    for (j = 0; j < pm->ncols; j++)
+		fprintf(fp, " %.4e", pricerow(pm,i)[j]);
+	    fputc('\n', fp);
+	}
+	fprintf(fp, "Solution matrix (transport plan):\n");
+	for (i = 0; i < mtx_nrows(me); i++) {
+	    for (j = 0; j < mtx_ncols(me); j++)
+		fprintf(fp, " %.4e", mtx_coef(me,i,j));
+	    fputc('\n', fp);
+	}
+	fclose(fp);
+}
+#endif
+
 /* Compute and insert migration along directed edge (may fork child) */
 static MIGRATION *
 create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
@@ -389,6 +474,9 @@ create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
 	    for (j = to_rbf->nrbf; j--; )
 		mtx_coef(newmig,i,j) *= nf;	/* row now sums to 1.0 */
 	}
+#ifdef DUMP_MATRIX
+	dump_matrix(newmig, &pmtx);
+#endif
 	end_subprocess();			/* exit here if subprocess */
 	free_routes(&pmtx);			/* free working arrays */
 	free(src_rem);
@@ -504,7 +592,7 @@ mesh_from_edge(MIGRATION *edge)
 static void
 check_normal_incidence(void)
 {
-	static const FVECT	norm_vec = {.0, .0, 1.};
+	static FVECT		norm_vec = {.0, .0, 1.};
 	const int		saved_nprocs = nprocs;
 	RBFNODE			*near_rbf, *mir_rbf, *rbf;
 	double			bestd;
@@ -554,11 +642,12 @@ check_normal_incidence(void)
 	memcpy(mir_rbf, near_rbf, n);
 	mir_rbf->ord = near_rbf->ord - 1;	/* not used, I think */
 	mir_rbf->next = NULL;
+	mir_rbf->ejl = NULL;
 	rev_rbf_symmetry(mir_rbf, MIRROR_X|MIRROR_Y);
 	nprocs = 1;				/* compute migration matrix */
-	if (mig_list != create_migration(mir_rbf, near_rbf))
+	if (create_migration(mir_rbf, near_rbf) == NULL)
 		exit(1);			/* XXX should never happen! */
-						/* interpolate normal dist. */
+	norm_vec[2] = input_orient;		/* interpolate normal dist. */
 	rbf = e_advect_rbf(mig_list, norm_vec, 2*near_rbf->nrbf);
 	nprocs = saved_nprocs;			/* final clean-up */
 	free(mir_rbf);