ViewVC Help

View File | Revision Log | Show Annotations | Download File | Root Listing

root/Development/ray/src/cv/bsdfmesh.c

Comparing ray/src/cv/bsdfmesh.c (file contents):
Revision 2.1 by greg, Fri Oct 19 04:14:29 2012 UTC vs.
Revision 2.26 by greg, Wed Mar 26 00:11:30 2014 UTC

#	Line 18 | Line 18 | static const char RCSid[] = "$Id$";
18		#include <string.h>
19		#include <math.h>
20		#include "bsdfrep.h"
21	+
22	+	#ifndef NEIGH_FACT2
23	+	#define NEIGH_FACT2     0.1     /* empirical neighborhood distance weight */
24	+	#endif
25		/* number of processes to run */
26		int                     nprocs = 1;
27		/* number of children (-1 in child) */
28		static int              nchild = 0;
29
30	<	/* Compute (and allocate) migration price matrix for optimization */
31	<	static float *
32	<	price_routes(const RBFNODE *from_rbf, const RBFNODE *to_rbf)
33	<	{
34	<	float   *pmtx = (float *)malloc(sizeof(float) *
35	<	from_rbf->nrbf * to_rbf->nrbf);
32	<	FVECT   *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
33	<	int     i, j;
30	>	typedef struct {
31	>	int             nrows, ncols;   /* array size (matches migration) */
32	>	float           *price;         /* migration prices */
33	>	short           *sord;          /* sort for each row, low to high */
34	>	float           *prow;          /* current price row */
35	>	} PRICEMAT;                     /* sorted pricing matrix */
36
37	<	if ((pmtx == NULL) | (vto == NULL)) {
38	<	fprintf(stderr, "%s: Out of memory in migration_costs()\n",
37	<	progname);
38	<	exit(1);
39	<	}
40	<	for (j = to_rbf->nrbf; j--; )           /* save repetitive ops. */
41	<	ovec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);
37	>	#define pricerow(p,i)   ((p)->price + (i)*(p)->ncols)
38	>	#define psortrow(p,i)   ((p)->sord + (i)*(p)->ncols)
39
43	–	for (i = from_rbf->nrbf; i--; ) {
44	–	const double        from_ang = R2ANG(from_rbf->rbfa[i].crad);
45	–	FVECT               vfrom;
46	–	ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
47	–	for (j = to_rbf->nrbf; j--; )
48	–	pmtx[i*to_rbf->nrbf + j] = acos(DOT(vfrom, vto[j])) +
49	–	fabs(R2ANG(to_rbf->rbfa[j].crad) - from_ang);
50	–	}
51	–	free(vto);
52	–	return(pmtx);
53	–	}
54	–
55	–	/* Comparison routine needed for sorting price row */
56	–	static const float      *price_arr;
57	–	static int
58	–	msrt_cmp(const void *p1, const void *p2)
59	–	{
60	–	float   c1 = price_arr[*(const int *)p1];
61	–	float   c2 = price_arr[*(const int *)p2];
62	–
63	–	if (c1 > c2) return(1);
64	–	if (c1 < c2) return(-1);
65	–	return(0);
66	–	}
67	–
68	–	/* Compute minimum (optimistic) cost for moving the given source material */
69	–	static double
70	–	min_cost(double amt2move, const double *avail, const float *price, int n)
71	–	{
72	–	static int      *price_sort = NULL;
73	–	static int      n_alloc = 0;
74	–	double          total_cost = 0;
75	–	int             i;
76	–
77	–	if (amt2move <= FTINY)                  /* pre-emptive check */
78	–	return(0.);
79	–	if (n > n_alloc) {                      /* (re)allocate sort array */
80	–	if (n_alloc) free(price_sort);
81	–	price_sort = (int *)malloc(sizeof(int)*n);
82	–	if (price_sort == NULL) {
83	–	fprintf(stderr, "%s: Out of memory in min_cost()\n",
84	–	progname);
85	–	exit(1);
86	–	}
87	–	n_alloc = n;
88	–	}
89	–	for (i = n; i--; )
90	–	price_sort[i] = i;
91	–	price_arr = price;
92	–	qsort(price_sort, n, sizeof(int), &msrt_cmp);
93	–	/* move cheapest first */
94	–	for (i = 0; i < n && amt2move > FTINY; i++) {
95	–	int     d = price_sort[i];
96	–	double  amt = (amt2move < avail[d]) ? amt2move : avail[d];
97	–
98	–	total_cost += amt * price[d];
99	–	amt2move -= amt;
100	–	}
101	–	return(total_cost);
102	–	}
103	–
104	–	/* Take a step in migration by choosing optimal bucket to transfer */
105	–	static double
106	–	migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, const float *pmtx)
107	–	{
108	–	const double    maxamt = .1;
109	–	const double    minamt = maxamt*.0001;
110	–	static double   *src_cost = NULL;
111	–	static int      n_alloc = 0;
112	–	struct {
113	–	int     s, d;   /* source and destination */
114	–	double  price;  /* price estimate per amount moved */
115	–	double  amt;    /* amount we can move */
116	–	} cur, best;
117	–	int             i;
118	–
119	–	if (mtx_nrows(mig) > n_alloc) {         /* allocate cost array */
120	–	if (n_alloc)
121	–	free(src_cost);
122	–	src_cost = (double *)malloc(sizeof(double)*mtx_nrows(mig));
123	–	if (src_cost == NULL) {
124	–	fprintf(stderr, "%s: Out of memory in migration_step()\n",
125	–	progname);
126	–	exit(1);
127	–	}
128	–	n_alloc = mtx_nrows(mig);
129	–	}
130	–	for (i = mtx_nrows(mig); i--; )         /* starting costs for diff. */
131	–	src_cost[i] = min_cost(src_rem[i], dst_rem,
132	–	pmtx+i*mtx_ncols(mig), mtx_ncols(mig));
133	–
134	–	/* find best source & dest. */
135	–	best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
136	–	for (cur.s = mtx_nrows(mig); cur.s--; ) {
137	–	const float *price = pmtx + cur.s*mtx_ncols(mig);
138	–	double      cost_others = 0;
139	–	if (src_rem[cur.s] < minamt)
140	–	continue;
141	–	cur.d = -1;                         /* examine cheapest dest. */
142	–	for (i = mtx_ncols(mig); i--; )
143	–	if (dst_rem[i] > minamt &&
144	–	(cur.d < 0 || price[i] < price[cur.d]))
145	–	cur.d = i;
146	–	if (cur.d < 0)
147	–	return(.0);
148	–	if ((cur.price = price[cur.d]) >= best.price)
149	–	continue;                   /* no point checking further */
150	–	cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
151	–	src_rem[cur.s] : dst_rem[cur.d];
152	–	if (cur.amt > maxamt) cur.amt = maxamt;
153	–	dst_rem[cur.d] -= cur.amt;          /* add up differential costs */
154	–	for (i = mtx_nrows(mig); i--; )
155	–	if (i != cur.s)
156	–	cost_others += min_cost(src_rem[i], dst_rem,
157	–	price, mtx_ncols(mig))
158	–	- src_cost[i];
159	–	dst_rem[cur.d] += cur.amt;          /* undo trial move */
160	–	cur.price += cost_others/cur.amt;   /* adjust effective price */
161	–	if (cur.price < best.price)         /* are we better than best? */
162	–	best = cur;
163	–	}
164	–	if ((best.s < 0) | (best.d < 0))
165	–	return(.0);
166	–	/* make the actual move */
167	–	mig->mtx[mtx_ndx(mig,best.s,best.d)] += best.amt;
168	–	src_rem[best.s] -= best.amt;
169	–	dst_rem[best.d] -= best.amt;
170	–	return(best.amt);
171	–	}
172	–
173	–	#ifdef DEBUG
174	–	static char *
175	–	thetaphi(const FVECT v)
176	–	{
177	–	static char     buf[128];
178	–	double          theta, phi;
179	–
180	–	theta = 180./M_PI*acos(v[2]);
181	–	phi = 180./M_PI*atan2(v[1],v[0]);
182	–	sprintf(buf, "(%.0f,%.0f)", theta, phi);
183	–
184	–	return(buf);
185	–	}
186	–	#endif
187	–
40		/* Create a new migration holder (sharing memory for multiprocessing) */
41		static MIGRATION *
42		new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
#	Line 272 | Line 124 | run_subprocess(void)
124		if (pid < 0) {
125		fprintf(stderr, "%s: cannot fork subprocess\n",
126		progname);
127	+	await_children(nchild);
128		exit(1);
129		}
130		++nchild;                       /* subprocess started */
#	Line 286 | Line 139 | run_subprocess(void)
139
140		#endif  /* ! _WIN32 */
141
142	+	/* Compute normalized distribution scattering functions for comparison */
143	+	static void
144	+	compute_nDSFs(const RBFNODE *rbf0, const RBFNODE *rbf1)
145	+	{
146	+	const double    nf0 = (GRIDRES*GRIDRES) / rbf0->vtotal;
147	+	const double    nf1 = (GRIDRES*GRIDRES) / rbf1->vtotal;
148	+	int             x, y;
149	+	FVECT           dv;
150	+
151	+	for (x = GRIDRES; x--; )
152	+	for (y = GRIDRES; y--; ) {
153	+	ovec_from_pos(dv, x, y);        /* cube root (brightness) */
154	+	dsf_grid[x][y].val[0] = pow(nf0*eval_rbfrep(rbf0, dv), .3333);
155	+	dsf_grid[x][y].val[1] = pow(nf1*eval_rbfrep(rbf1, dv), .3333);
156	+	}
157	+	}
158	+
159	+	/* Compute neighborhood distance-squared (dissimilarity) */
160	+	static double
161	+	neighborhood_dist2(int x0, int y0, int x1, int y1)
162	+	{
163	+	int     rad = GRIDRES>>5;
164	+	double  sum2 = 0.;
165	+	double  d;
166	+	int     p[4];
167	+	int     i, j;
168	+	/* check radius */
169	+	p[0] = x0; p[1] = y0; p[2] = x1; p[3] = y1;
170	+	for (i = 4; i--; ) {
171	+	if (p[i] < rad) rad = p[i];
172	+	if (GRIDRES-1-p[i] < rad) rad = GRIDRES-1-p[i];
173	+	}
174	+	for (i = -rad; i <= rad; i++)
175	+	for (j = -rad; j <= rad; j++) {
176	+	d = dsf_grid[x0+i][y0+j].val[0] -
177	+	dsf_grid[x1+i][y1+j].val[1];
178	+	sum2 += d*d;
179	+	}
180	+	return(sum2 / (4*rad*(rad+1) + 1));
181	+	}
182	+
183	+	/* Comparison routine needed for sorting price row */
184	+	static int
185	+	msrt_cmp(void *b, const void *p1, const void *p2)
186	+	{
187	+	PRICEMAT        *pm = (PRICEMAT *)b;
188	+	float           c1 = pm->prow[*(const short *)p1];
189	+	float           c2 = pm->prow[*(const short *)p2];
190	+
191	+	if (c1 > c2) return(1);
192	+	if (c1 < c2) return(-1);
193	+	return(0);
194	+	}
195	+
196	+	/* Compute (and allocate) migration price matrix for optimization */
197	+	static void
198	+	price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, const RBFNODE *to_rbf)
199	+	{
200	+	FVECT   *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
201	+	int     i, j;
202	+
203	+	compute_nDSFs(from_rbf, to_rbf);
204	+	pm->nrows = from_rbf->nrbf;
205	+	pm->ncols = to_rbf->nrbf;
206	+	pm->price = (float *)malloc(sizeof(float) * pm->nrows*pm->ncols);
207	+	pm->sord = (short *)malloc(sizeof(short) * pm->nrows*pm->ncols);
208	+
209	+	if ((pm->price == NULL) | (pm->sord == NULL) | (vto == NULL)) {
210	+	fprintf(stderr, "%s: Out of memory in migration_costs()\n",
211	+	progname);
212	+	exit(1);
213	+	}
214	+	for (j = to_rbf->nrbf; j--; )           /* save repetitive ops. */
215	+	ovec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);
216	+
217	+	for (i = from_rbf->nrbf; i--; ) {
218	+	const double        from_ang = R2ANG(from_rbf->rbfa[i].crad);
219	+	FVECT               vfrom;
220	+	short               *srow;
221	+	ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
222	+	pm->prow = pricerow(pm,i);
223	+	srow = psortrow(pm,i);
224	+	for (j = to_rbf->nrbf; j--; ) {
225	+	double  d;                      /* quadratic cost function */
226	+	d = Acos(DOT(vfrom, vto[j]));
227	+	pm->prow[j] = d*d;
228	+	d = R2ANG(to_rbf->rbfa[j].crad) - from_ang;
229	+	pm->prow[j] += d*d;
230	+	/* neighborhood difference */
231	+	pm->prow[j] += NEIGH_FACT2 * neighborhood_dist2(
232	+	from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy,
233	+	to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy );
234	+	srow[j] = j;
235	+	}
236	+	qsort_r(srow, pm->ncols, sizeof(short), pm, &msrt_cmp);
237	+	}
238	+	free(vto);
239	+	}
240	+
241	+	/* Free price matrix */
242	+	static void
243	+	free_routes(PRICEMAT *pm)
244	+	{
245	+	free(pm->price); pm->price = NULL;
246	+	free(pm->sord); pm->sord = NULL;
247	+	}
248	+
249	+	/* Compute minimum (optimistic) cost for moving the given source material */
250	+	static double
251	+	min_cost(double amt2move, const double *avail, const PRICEMAT *pm, int s)
252	+	{
253	+	const short     *srow = psortrow(pm,s);
254	+	const float     *prow = pricerow(pm,s);
255	+	double          total_cost = 0;
256	+	int             j;
257	+	/* move cheapest first */
258	+	for (j = 0; (j < pm->ncols) & (amt2move > FTINY); j++) {
259	+	int     d = srow[j];
260	+	double  amt = (amt2move < avail[d]) ? amt2move : avail[d];
261	+
262	+	total_cost += amt * prow[d];
263	+	amt2move -= amt;
264	+	}
265	+	return(total_cost);
266	+	}
267	+
268	+	typedef struct {
269	+	short   s, d;           /* source and destination */
270	+	float   dc;             /* discount to push inventory */
271	+	} ROWSENT;              /* row sort entry */
272	+
273	+	/* Compare entries by discounted moving price */
274	+	static int
275	+	rmovcmp(void *b, const void *p1, const void *p2)
276	+	{
277	+	PRICEMAT        *pm = (PRICEMAT *)b;
278	+	const ROWSENT   *re1 = (const ROWSENT *)p1;
279	+	const ROWSENT   *re2 = (const ROWSENT *)p2;
280	+	double          price_diff;
281	+
282	+	if (re1->d < 0) return(re2->d >= 0);
283	+	if (re2->d < 0) return(-1);
284	+	price_diff = re1->dc*pricerow(pm,re1->s)[re1->d] -
285	+	re2->dc*pricerow(pm,re2->s)[re2->d];
286	+	if (price_diff > 0) return(1);
287	+	if (price_diff < 0) return(-1);
288	+	return(0);
289	+	}
290	+
291	+	/* Take a step in migration by choosing reasonable bucket to transfer */
292	+	static double
293	+	migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, PRICEMAT *pm)
294	+	{
295	+	const int       max2check = 100;
296	+	const double    maxamt = 1./(double)pm->ncols;
297	+	const double    minamt = maxamt*1e-4;
298	+	double          *src_cost;
299	+	ROWSENT         *rord;
300	+	struct {
301	+	int     s, d;   /* source and destination */
302	+	double  price;  /* cost per amount moved */
303	+	double  amt;    /* amount we can move */
304	+	} cur, best;
305	+	int             r2check, i, ri;
306	+	/*
307	+	* Check cheapest available routes only -- a higher adjusted
308	+	* destination price implies that another source is closer, so
309	+	* we can hold off considering more expensive options until
310	+	* some other (hopefully better) moves have been made.
311	+	* A discount based on source remaining is supposed to prioritize
312	+	* movement from large lobes, but it doesn't seem to do much,
313	+	* so we have it set to 1.0 at the moment.
314	+	*/
315	+	#define discount(qr)    1.0
316	+	/* most promising row order */
317	+	rord = (ROWSENT *)malloc(sizeof(ROWSENT)*pm->nrows);
318	+	if (rord == NULL)
319	+	goto memerr;
320	+	for (ri = pm->nrows; ri--; ) {
321	+	rord[ri].s = ri;
322	+	rord[ri].d = -1;
323	+	rord[ri].dc = 1.f;
324	+	if (src_rem[ri] <= minamt)          /* enough source material? */
325	+	continue;
326	+	for (i = 0; i < pm->ncols; i++)
327	+	if (dst_rem[ rord[ri].d = psortrow(pm,ri)[i] ] > minamt)
328	+	break;
329	+	if (i >= pm->ncols) {               /* moved all we can? */
330	+	free(rord);
331	+	return(.0);
332	+	}
333	+	rord[ri].dc = discount(src_rem[ri]);
334	+	}
335	+	if (pm->nrows > max2check)              /* sort if too many sources */
336	+	qsort_r(rord, pm->nrows, sizeof(ROWSENT), pm, &rmovcmp);
337	+	/* allocate cost array */
338	+	src_cost = (double *)malloc(sizeof(double)*pm->nrows);
339	+	if (src_cost == NULL)
340	+	goto memerr;
341	+	for (i = pm->nrows; i--; )              /* starting costs for diff. */
342	+	src_cost[i] = min_cost(src_rem[i], dst_rem, pm, i);
343	+	/* find best source & dest. */
344	+	best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
345	+	if ((r2check = pm->nrows) > max2check)
346	+	r2check = max2check;            /* put a limit on search */
347	+	for (ri = 0; ri < r2check; ri++) {      /* check each source row */
348	+	double      cost_others = 0;
349	+	cur.s = rord[ri].s;
350	+	if ((cur.d = rord[ri].d) < 0 ||
351	+	rord[ri].dc*pricerow(pm,cur.s)[cur.d] >= best.price) {
352	+	if (pm->nrows > max2check) break;       /* sorted end */
353	+	continue;                       /* else skip this one */
354	+	}
355	+	cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
356	+	src_rem[cur.s] : dst_rem[cur.d];
357	+	/* don't just leave smidgen */
358	+	if (cur.amt > maxamt*1.02) cur.amt = maxamt;
359	+	dst_rem[cur.d] -= cur.amt;          /* add up opportunity costs */
360	+	for (i = pm->nrows; i--; )
361	+	if (i != cur.s)
362	+	cost_others += min_cost(src_rem[i], dst_rem, pm, i)
363	+	- src_cost[i];
364	+	dst_rem[cur.d] += cur.amt;          /* undo trial move */
365	+	/* discount effective price */
366	+	cur.price = ( pricerow(pm,cur.s)[cur.d] + cost_others/cur.amt ) *
367	+	rord[ri].dc;
368	+	if (cur.price < best.price)         /* are we better than best? */
369	+	best = cur;
370	+	}
371	+	free(src_cost);                         /* clean up */
372	+	free(rord);
373	+	if ((best.s < 0) | (best.d < 0))        /* nothing left to move? */
374	+	return(.0);
375	+	/* else make the actual move */
376	+	mtx_coef(mig,best.s,best.d) += best.amt;
377	+	src_rem[best.s] -= best.amt;
378	+	dst_rem[best.d] -= best.amt;
379	+	return(best.amt);
380	+	memerr:
381	+	fprintf(stderr, "%s: Out of memory in migration_step()\n", progname);
382	+	exit(1);
383	+	#undef discount
384	+	}
385	+
386	+	#ifdef DUMP_MATRIX
387	+	/* Dump transport plan and corresponding price matrix to a text file */
388	+	static void
389	+	dump_matrix(const MIGRATION *me, const PRICEMAT *pm)
390	+	{
391	+	char    fname[256];
392	+	FILE    *fp;
393	+	int     i, j;
394	+
395	+	sprintf(fname, "edge_%d-%d.txt", me->rbfv[0]->ord, me->rbfv[1]->ord);
396	+	if ((fp = fopen(fname, "w")) == NULL)
397	+	return;
398	+	for (j = 0; j < 2; j++) {
399	+	fprintf(fp, "Available from %d source RBF lobes in node %d:\n",
400	+	me->rbfv[j]->nrbf, me->rbfv[j]->ord);
401	+	for (i = 0; i < me->rbfv[j]->nrbf; i++)
402	+	fprintf(fp, " %.4e", rbf_volume(&me->rbfv[j]->rbfa[i]) /
403	+	me->rbfv[j]->vtotal);
404	+	fputc('\n', fp);
405	+	}
406	+	fprintf(fp, "Price (quadratic distance metric) matrix:\n");
407	+	for (i = 0; i < pm->nrows; i++) {
408	+	for (j = 0; j < pm->ncols; j++)
409	+	fprintf(fp, " %.4e", pricerow(pm,i)[j]);
410	+	fputc('\n', fp);
411	+	}
412	+	fprintf(fp, "Solution matrix (transport plan):\n");
413	+	for (i = 0; i < mtx_nrows(me); i++) {
414	+	for (j = 0; j < mtx_ncols(me); j++)
415	+	fprintf(fp, " %.4e", mtx_coef(me,i,j));
416	+	fputc('\n', fp);
417	+	}
418	+	fclose(fp);
419	+	}
420	+	#endif
421	+
422		/* Compute and insert migration along directed edge (may fork child) */
423		static MIGRATION *
424		create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
425		{
426	<	const double    end_thresh = 0.1/(from_rbf->nrbf*to_rbf->nrbf);
427	<	const double    check_thresh = 0.01;
295	<	const double    rel_thresh = 5e-6;
296	<	float           *pmtx;
426	>	const double    end_thresh = 5e-6;
427	>	PRICEMAT        pmtx;
428		MIGRATION       *newmig;
429		double          *src_rem, *dst_rem;
430		double          total_rem = 1., move_amt;
431	<	int             i;
431	>	int             i, j;
432		/* check if exists already */
433		for (newmig = from_rbf->ejl; newmig != NULL;
434		newmig = nextedge(from_rbf,newmig))
435		if (newmig->rbfv[1] == to_rbf)
436		return(NULL);
437		/* else allocate */
438	+	#ifdef DEBUG
439	+	fprintf(stderr, "Building path from (theta,phi) (%.1f,%.1f) ",
440	+	get_theta180(from_rbf->invec),
441	+	get_phi360(from_rbf->invec));
442	+	fprintf(stderr, "to (%.1f,%.1f) with %d x %d matrix\n",
443	+	get_theta180(to_rbf->invec),
444	+	get_phi360(to_rbf->invec),
445	+	from_rbf->nrbf, to_rbf->nrbf);
446	+	#endif
447		newmig = new_migration(from_rbf, to_rbf);
448		if (run_subprocess())
449		return(newmig);                 /* child continues */
450	<	pmtx = price_routes(from_rbf, to_rbf);
450	>	price_routes(&pmtx, from_rbf, to_rbf);
451		src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
452		dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
453		if ((src_rem == NULL) | (dst_rem == NULL)) {
#	Line 315 | Line 455 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
455		progname);
456		exit(1);
457		}
318	–	#ifdef DEBUG
319	–	fprintf(stderr, "Building path from (theta,phi) %s ",
320	–	thetaphi(from_rbf->invec));
321	–	fprintf(stderr, "to %s", thetaphi(to_rbf->invec));
322	–	/* if (nchild) */ fputc('\n', stderr);
323	–	#endif
458		/* starting quantities */
459		memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
460		for (i = from_rbf->nrbf; i--; )
461		src_rem[i] = rbf_volume(&from_rbf->rbfa[i]) / from_rbf->vtotal;
462	<	for (i = to_rbf->nrbf; i--; )
463	<	dst_rem[i] = rbf_volume(&to_rbf->rbfa[i]) / to_rbf->vtotal;
462	>	for (j = to_rbf->nrbf; j--; )
463	>	dst_rem[j] = rbf_volume(&to_rbf->rbfa[j]) / to_rbf->vtotal;
464	>
465		do {                                    /* move a bit at a time */
466	<	move_amt = migration_step(newmig, src_rem, dst_rem, pmtx);
466	>	move_amt = migration_step(newmig, src_rem, dst_rem, &pmtx);
467		total_rem -= move_amt;
468	<	#ifdef DEBUG
469	<	if (!nchild)
335	<	/* fputc('.', stderr); */
336	<	fprintf(stderr, "%.9f remaining...\r", total_rem);
337	<	#endif
338	<	} while (total_rem > end_thresh && (total_rem > check_thresh) |
339	<	(move_amt > rel_thresh*total_rem));
340	<	#ifdef DEBUG
341	<	if (!nchild) fputs("\ndone.\n", stderr);
342	<	else fprintf(stderr, "finished with %.9f remaining\n", total_rem);
343	<	#endif
468	>	} while ((total_rem > end_thresh) & (move_amt > 0));
469	>
470		for (i = from_rbf->nrbf; i--; ) {       /* normalize final matrix */
471	<	float       nf = rbf_volume(&from_rbf->rbfa[i]);
346	<	int         j;
471	>	double      nf = rbf_volume(&from_rbf->rbfa[i]);
472		if (nf <= FTINY) continue;
473		nf = from_rbf->vtotal / nf;
474		for (j = to_rbf->nrbf; j--; )
475	<	newmig->mtx[mtx_ndx(newmig,i,j)] *= nf;
475	>	mtx_coef(newmig,i,j) *= nf;     /* row now sums to 1.0 */
476		}
477	+	#ifdef DUMP_MATRIX
478	+	dump_matrix(newmig, &pmtx);
479	+	#endif
480		end_subprocess();                       /* exit here if subprocess */
481	<	free(pmtx);                             /* free working arrays */
481	>	free_routes(&pmtx);                     /* free working arrays */
482		free(src_rem);
483		free(dst_rem);
484		return(newmig);
#	Line 382 | Line 510 | overlaps_tri(const RBFNODE *bv0, const RBFNODE *bv1, c
510		return(vother[im_rev] != NULL);
511		}
512
513	<	/* Find context hull vertex to complete triangle (oriented call) */
513	>	/* Find convex hull vertex to complete triangle (oriented call) */
514		static RBFNODE *
515		find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rbf1)
516		{
#	Line 403 | Line 531 | find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rb
531		if (DOT(vp, vmid) <= FTINY)
532		continue;               /* wrong orientation */
533		area2 = .25*DOT(vp,vp);
534	<	VSUB(vp, rbf->invec, rbf0->invec);
534	>	VSUB(vp, rbf->invec, vmid);
535		dprod = -DOT(vp, vejn);
536		VSUM(vp, vp, vejn, dprod);      /* above guarantees non-zero */
537		dprod = DOT(vp, vmid) / VLEN(vp);
#	Line 459 | Line 587 | mesh_from_edge(MIGRATION *edge)
587		}
588		}
589		}
590	+
591	+	/* Add normal direction if missing */
592	+	static void
593	+	check_normal_incidence(void)
594	+	{
595	+	static FVECT            norm_vec = {.0, .0, 1.};
596	+	const int               saved_nprocs = nprocs;
597	+	RBFNODE                 *near_rbf, *mir_rbf, *rbf;
598	+	double                  bestd;
599	+	int                     n;
600	+
601	+	if (dsf_list == NULL)
602	+	return;                         /* XXX should be error? */
603	+	near_rbf = dsf_list;
604	+	bestd = input_orient*near_rbf->invec[2];
605	+	if (single_plane_incident) {            /* ordered plane incidence? */
606	+	if (bestd >= 1.-2.*FTINY)
607	+	return;                 /* already have normal */
608	+	} else {
609	+	switch (inp_coverage) {
610	+	case INP_QUAD1:
611	+	case INP_QUAD2:
612	+	case INP_QUAD3:
613	+	case INP_QUAD4:
614	+	break;                  /* quadrilateral symmetry? */
615	+	default:
616	+	return;                 /* else we can interpolate */
617	+	}
618	+	for (rbf = near_rbf->next; rbf != NULL; rbf = rbf->next) {
619	+	const double    d = input_orient*rbf->invec[2];
620	+	if (d >= 1.-2.*FTINY)
621	+	return;         /* seems we have normal */
622	+	if (d > bestd) {
623	+	near_rbf = rbf;
624	+	bestd = d;
625	+	}
626	+	}
627	+	}
628	+	if (mig_list != NULL) {                 /* need to be called first */
629	+	fprintf(stderr, "%s: Late call to check_normal_incidence()\n",
630	+	progname);
631	+	exit(1);
632	+	}
633	+	#ifdef DEBUG
634	+	fprintf(stderr, "Interpolating normal incidence by mirroring (%.1f,%.1f)\n",
635	+	get_theta180(near_rbf->invec), get_phi360(near_rbf->invec));
636	+	#endif
637	+	/* mirror nearest incidence */
638	+	n = sizeof(RBFNODE) + sizeof(RBFVAL)*(near_rbf->nrbf-1);
639	+	mir_rbf = (RBFNODE *)malloc(n);
640	+	if (mir_rbf == NULL)
641	+	goto memerr;
642	+	memcpy(mir_rbf, near_rbf, n);
643	+	mir_rbf->ord = near_rbf->ord - 1;       /* not used, I think */
644	+	mir_rbf->next = NULL;
645	+	mir_rbf->ejl = NULL;
646	+	rev_rbf_symmetry(mir_rbf, MIRROR_X|MIRROR_Y);
647	+	nprocs = 1;                             /* compute migration matrix */
648	+	if (create_migration(mir_rbf, near_rbf) == NULL)
649	+	exit(1);                        /* XXX should never happen! */
650	+	norm_vec[2] = input_orient;             /* interpolate normal dist. */
651	+	rbf = e_advect_rbf(mig_list, norm_vec, 2*near_rbf->nrbf);
652	+	nprocs = saved_nprocs;                  /* final clean-up */
653	+	free(mir_rbf);
654	+	free(mig_list);
655	+	mig_list = near_rbf->ejl = NULL;
656	+	insert_dsf(rbf);                        /* insert interpolated normal */
657	+	return;
658	+	memerr:
659	+	fprintf(stderr, "%s: Out of memory in check_normal_incidence()\n",
660	+	progname);
661	+	exit(1);
662	+	}
663
664		/* Build our triangle mesh from recorded RBFs */
665		void
#	Line 467 | Line 668 | build_mesh(void)
668		double          best2 = M_PI*M_PI;
669		RBFNODE         *shrt_edj[2];
670		RBFNODE         *rbf0, *rbf1;
671	+	/* add normal if needed */
672	+	check_normal_incidence();
673		/* check if isotropic */
674		if (single_plane_incident) {
675		for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines (old)
>	Changed lines (new)