--- ray/src/cv/bsdfmesh.c	2012/11/08 22:05:04	2.3
+++ ray/src/cv/bsdfmesh.c	2013/06/28 23:18:51	2.9
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdfmesh.c,v 2.3 2012/11/08 22:05:04 greg Exp $";
+static const char RCSid[] = "$Id: bsdfmesh.c,v 2.9 2013/06/28 23:18:51 greg Exp $";
 #endif
 /*
  * Create BSDF advection mesh from radial basis functions.
@@ -119,6 +119,7 @@ run_subprocess(void)
 		if (pid < 0) {
 			fprintf(stderr, "%s: cannot fork subprocess\n",
 					progname);
+			await_children(nchild);
 			exit(1);
 		}
 		++nchild;			/* subprocess started */
@@ -172,7 +173,8 @@ price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, co
 	    FVECT		vfrom;
 	    ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
 	    for (j = to_rbf->nrbf; j--; ) {
-		pricerow(pm,i)[j] = acos(DOT(vfrom, vto[j])) +
+		double		dprod = DOT(vfrom, vto[j]);
+		pricerow(pm,i)[j] = ((dprod >= 1.) ? .0 : acos(dprod)) +
 				fabs(R2ANG(to_rbf->rbfa[j].crad) - from_ang);
 		psortrow(pm,i)[j] = j;
 	    }
@@ -197,7 +199,7 @@ min_cost(double amt2move, const double *avail, const P
 	int		j;
 
 	if (amt2move <= FTINY)			/* pre-emptive check */
-		return(0.);
+		return(.0);
 						/* move cheapest first */
 	for (j = 0; j < pm->ncols && amt2move > FTINY; j++) {
 		int	d = psortrow(pm,s)[j];
@@ -213,27 +215,21 @@ min_cost(double amt2move, const double *avail, const P
 static double
 migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, const PRICEMAT *pm)
 {
-	const double	maxamt = .1;
+	const double	maxamt = 1./(double)pm->ncols;
 	const double	minamt = maxamt*5e-6;
-	static double	*src_cost = NULL;
-	static int	n_alloc = 0;
+	double		*src_cost;
 	struct {
 		int	s, d;	/* source and destination */
 		double	price;	/* price estimate per amount moved */
 		double	amt;	/* amount we can move */
 	} cur, best;
 	int		i;
-
-	if (pm->nrows > n_alloc) {		/* allocate cost array */
-		if (n_alloc)
-			free(src_cost);
-		src_cost = (double *)malloc(sizeof(double)*pm->nrows);
-		if (src_cost == NULL) {
-			fprintf(stderr, "%s: Out of memory in migration_step()\n",
-					progname);
-			exit(1);
-		}
-		n_alloc = pm->nrows;
+						/* allocate cost array */
+	src_cost = (double *)malloc(sizeof(double)*pm->nrows);
+	if (src_cost == NULL) {
+		fprintf(stderr, "%s: Out of memory in migration_step()\n",
+				progname);
+		exit(1);
 	}
 	for (i = pm->nrows; i--; )		/* starting costs for diff. */
 		src_cost[i] = min_cost(src_rem[i], dst_rem, pm, i);
@@ -241,19 +237,18 @@ migration_step(MIGRATION *mig, double *src_rem, double
 						/* find best source & dest. */
 	best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
 	for (cur.s = pm->nrows; cur.s--; ) {
-	    const float	*price = pricerow(pm,cur.s);
 	    double	cost_others = 0;
+
 	    if (src_rem[cur.s] <= minamt)
 		    continue;
-	    cur.d = -1;				/* examine cheapest dest. */
-	    for (i = pm->ncols; i--; )
-		if (dst_rem[i] > minamt &&
-				(cur.d < 0 || price[i] < price[cur.d]))
-			cur.d = i;
-	    if (cur.d < 0)
-		    return(.0);
-	    if ((cur.price = price[cur.d]) >= best.price)
-		    continue;			/* no point checking further */
+						/* examine cheapest dest. */
+	    for (i = 0; i < pm->ncols; i++)
+		if (dst_rem[ cur.d = psortrow(pm,cur.s)[i] ] > minamt)
+			break;
+	    if (i >= pm->ncols)
+		break;
+	    if ((cur.price = pricerow(pm,cur.s)[cur.d]) >= best.price)
+		continue;			/* no point checking further */
 	    cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
 				src_rem[cur.s] : dst_rem[cur.d];
 	    if (cur.amt > maxamt) cur.amt = maxamt;
@@ -267,30 +262,17 @@ migration_step(MIGRATION *mig, double *src_rem, double
 	    if (cur.price < best.price)		/* are we better than best? */
 		    best = cur;
 	}
-	if ((best.s < 0) | (best.d < 0))
+	free(src_cost);				/* finish up */
+
+	if ((best.s < 0) | (best.d < 0))	/* nothing left to move? */
 		return(.0);
-						/* make the actual move */
+						/* else make the actual move */
 	mtx_coef(mig,best.s,best.d) += best.amt;
 	src_rem[best.s] -= best.amt;
 	dst_rem[best.d] -= best.amt;
 	return(best.amt);
 }
 
-#ifdef DEBUG
-static char *
-thetaphi(const FVECT v)
-{
-	static char	buf[128];
-	double		theta, phi;
-
-	theta = 180./M_PI*acos(v[2]);
-	phi = 180./M_PI*atan2(v[1],v[0]);
-	sprintf(buf, "(%.0f,%.0f)", theta, phi);
-
-	return(buf);
-}
-#endif
-
 /* Compute and insert migration along directed edge (may fork child) */
 static MIGRATION *
 create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
@@ -300,13 +282,22 @@ create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
 	MIGRATION	*newmig;
 	double		*src_rem, *dst_rem;
 	double		total_rem = 1., move_amt;
-	int		i;
+	int		i, j;
 						/* check if exists already */
 	for (newmig = from_rbf->ejl; newmig != NULL;
 			newmig = nextedge(from_rbf,newmig))
 		if (newmig->rbfv[1] == to_rbf)
 			return(NULL);
 						/* else allocate */
+#ifdef DEBUG
+	fprintf(stderr, "Building path from (theta,phi) (%.0f,%.0f) ",
+			get_theta180(from_rbf->invec),
+			get_phi360(from_rbf->invec));
+	fprintf(stderr, "to (%.0f,%.0f) with %d x %d matrix\n",
+			get_theta180(to_rbf->invec),
+			get_phi360(to_rbf->invec), 
+			from_rbf->nrbf, to_rbf->nrbf);
+#endif
 	newmig = new_migration(from_rbf, to_rbf);
 	if (run_subprocess())
 		return(newmig);			/* child continues */
@@ -318,38 +309,24 @@ create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
 				progname);
 		exit(1);
 	}
-#ifdef DEBUG
-	fprintf(stderr, "Building path from (theta,phi) %s ",
-			thetaphi(from_rbf->invec));
-	fprintf(stderr, "to %s with %d x %d matrix\n",
-			thetaphi(to_rbf->invec), 
-			from_rbf->nrbf, to_rbf->nrbf);
-#endif
 						/* starting quantities */
 	memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
 	for (i = from_rbf->nrbf; i--; )
 		src_rem[i] = rbf_volume(&from_rbf->rbfa[i]) / from_rbf->vtotal;
-	for (i = to_rbf->nrbf; i--; )
-		dst_rem[i] = rbf_volume(&to_rbf->rbfa[i]) / to_rbf->vtotal;
+	for (j = to_rbf->nrbf; j--; )
+		dst_rem[j] = rbf_volume(&to_rbf->rbfa[j]) / to_rbf->vtotal;
+
 	do {					/* move a bit at a time */
 		move_amt = migration_step(newmig, src_rem, dst_rem, &pmtx);
 		total_rem -= move_amt;
-#ifdef DEBUG
-		if (!nchild)
-			fprintf(stderr, "\r%.9f remaining...", total_rem);
-#endif
 	} while ((total_rem > end_thresh) & (move_amt > 0));
-#ifdef DEBUG
-	if (!nchild) fputs("done.\n", stderr);
-	else fprintf(stderr, "finished with %.9f remaining\n", total_rem);
-#endif
+
 	for (i = from_rbf->nrbf; i--; ) {	/* normalize final matrix */
-	    float	nf = rbf_volume(&from_rbf->rbfa[i]);
-	    int		j;
+	    double	nf = rbf_volume(&from_rbf->rbfa[i]);
 	    if (nf <= FTINY) continue;
 	    nf = from_rbf->vtotal / nf;
 	    for (j = to_rbf->nrbf; j--; )
-		mtx_coef(newmig,i,j) *= nf;
+		mtx_coef(newmig,i,j) *= nf;	/* row now sums to 1.0 */
 	}
 	end_subprocess();			/* exit here if subprocess */
 	free_routes(&pmtx);			/* free working arrays */