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.15 by greg, Tue Feb 18 16:06:51 2014 UTC

#	Line 23 | Line 23 | int                    nprocs = 1;
23		/* number of children (-1 in child) */
24		static int              nchild = 0;
25
26	<	/* Compute (and allocate) migration price matrix for optimization */
27	<	static float *
28	<	price_routes(const RBFNODE *from_rbf, const RBFNODE *to_rbf)
29	<	{
30	<	float   *pmtx = (float *)malloc(sizeof(float) *
31	<	from_rbf->nrbf * to_rbf->nrbf);
32	<	FVECT   *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
33	<	int     i, j;
26	>	typedef struct {
27	>	int             nrows, ncols;   /* array size (matches migration) */
28	>	float           *price;         /* migration prices */
29	>	short           *sord;          /* sort for each row, low to high */
30	>	float           *prow;          /* current price row */
31	>	} PRICEMAT;                     /* sorted pricing matrix */
32
33	<	if ((pmtx == NULL) | (vto == NULL)) {
34	<	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);
33	>	#define pricerow(p,i)   ((p)->price + (i)*(p)->ncols)
34	>	#define psortrow(p,i)   ((p)->sord + (i)*(p)->ncols)
35
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	–
36		/* Create a new migration holder (sharing memory for multiprocessing) */
37		static MIGRATION *
38		new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
#	Line 272 | Line 120 | run_subprocess(void)
120		if (pid < 0) {
121		fprintf(stderr, "%s: cannot fork subprocess\n",
122		progname);
123	+	await_children(nchild);
124		exit(1);
125		}
126		++nchild;                       /* subprocess started */
#	Line 286 | Line 135 | run_subprocess(void)
135
136		#endif  /* ! _WIN32 */
137
138	+	/* Comparison routine needed for sorting price row */
139	+	static int
140	+	msrt_cmp(void *b, const void *p1, const void *p2)
141	+	{
142	+	PRICEMAT        *pm = (PRICEMAT *)b;
143	+	float           c1 = pm->prow[*(const short *)p1];
144	+	float           c2 = pm->prow[*(const short *)p2];
145	+
146	+	if (c1 > c2) return(1);
147	+	if (c1 < c2) return(-1);
148	+	return(0);
149	+	}
150	+
151	+	/* Compute (and allocate) migration price matrix for optimization */
152	+	static void
153	+	price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, const RBFNODE *to_rbf)
154	+	{
155	+	FVECT   *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
156	+	int     i, j;
157	+
158	+	pm->nrows = from_rbf->nrbf;
159	+	pm->ncols = to_rbf->nrbf;
160	+	pm->price = (float *)malloc(sizeof(float) * pm->nrows*pm->ncols);
161	+	pm->sord = (short *)malloc(sizeof(short) * pm->nrows*pm->ncols);
162	+
163	+	if ((pm->price == NULL) | (pm->sord == NULL) | (vto == NULL)) {
164	+	fprintf(stderr, "%s: Out of memory in migration_costs()\n",
165	+	progname);
166	+	exit(1);
167	+	}
168	+	for (j = to_rbf->nrbf; j--; )           /* save repetitive ops. */
169	+	ovec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);
170	+
171	+	for (i = from_rbf->nrbf; i--; ) {
172	+	const double        from_ang = R2ANG(from_rbf->rbfa[i].crad);
173	+	FVECT               vfrom;
174	+	short               *srow;
175	+	ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
176	+	pm->prow = pricerow(pm,i);
177	+	srow = psortrow(pm,i);
178	+	for (j = to_rbf->nrbf; j--; ) {
179	+	double          d;              /* quadratic cost function */
180	+	d = DOT(vfrom, vto[j]);
181	+	d = (d >= 1.) ? .0 : acos(d);
182	+	pm->prow[j] = d*d;
183	+	d = R2ANG(to_rbf->rbfa[j].crad) - from_ang;
184	+	pm->prow[j] += d*d;
185	+	srow[j] = j;
186	+	}
187	+	qsort_r(srow, pm->ncols, sizeof(short), pm, &msrt_cmp);
188	+	}
189	+	free(vto);
190	+	}
191	+
192	+	/* Free price matrix */
193	+	static void
194	+	free_routes(PRICEMAT *pm)
195	+	{
196	+	free(pm->price); pm->price = NULL;
197	+	free(pm->sord); pm->sord = NULL;
198	+	}
199	+
200	+	/* Compute minimum (optimistic) cost for moving the given source material */
201	+	static double
202	+	min_cost(double amt2move, const double *avail, const PRICEMAT *pm, int s)
203	+	{
204	+	const short     *srow = psortrow(pm,s);
205	+	const float     *prow = pricerow(pm,s);
206	+	double          total_cost = 0;
207	+	int             j;
208	+	/* move cheapest first */
209	+	for (j = 0; (j < pm->ncols) & (amt2move > FTINY); j++) {
210	+	int     d = srow[j];
211	+	double  amt = (amt2move < avail[d]) ? amt2move : avail[d];
212	+
213	+	total_cost += amt * prow[d];
214	+	amt2move -= amt;
215	+	}
216	+	return(total_cost);
217	+	}
218	+
219	+	/* Compare entries by moving price */
220	+	static int
221	+	rmovcmp(void *b, const void *p1, const void *p2)
222	+	{
223	+	PRICEMAT        *pm = (PRICEMAT *)b;
224	+	const short     *ij1 = (const short *)p1;
225	+	const short     *ij2 = (const short *)p2;
226	+	float           price_diff;
227	+
228	+	if (ij1[1] < 0) return(ij2[1] >= 0);
229	+	if (ij2[1] < 0) return(-1);
230	+	price_diff = pricerow(pm,ij1[0])[ij1[1]] - pricerow(pm,ij2[0])[ij2[1]];
231	+	if (price_diff > 0) return(1);
232	+	if (price_diff < 0) return(-1);
233	+	return(0);
234	+	}
235	+
236	+	/* Take a step in migration by choosing reasonable bucket to transfer */
237	+	static double
238	+	migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, PRICEMAT *pm)
239	+	{
240	+	const int       max2check = 100;
241	+	const double    maxamt = 1./(double)pm->ncols;
242	+	const double    minamt = maxamt*1e-4;
243	+	double          *src_cost;
244	+	short           (*rord)[2];
245	+	struct {
246	+	int     s, d;   /* source and destination */
247	+	double  price;  /* price estimate per amount moved */
248	+	double  amt;    /* amount we can move */
249	+	} cur, best;
250	+	int             r2check, i, ri;
251	+	/*
252	+	* Check cheapest available routes only -- a higher adjusted
253	+	* destination price implies that another source is closer, so
254	+	* we can hold off considering more expensive options until
255	+	* some other (hopefully better) moves have been made.
256	+	*/
257	+	/* most promising row order */
258	+	rord = (short (*)[2])malloc(sizeof(short)*2*pm->nrows);
259	+	if (rord == NULL)
260	+	goto memerr;
261	+	for (ri = pm->nrows; ri--; ) {
262	+	rord[ri][0] = ri;
263	+	rord[ri][1] = -1;
264	+	if (src_rem[ri] <= minamt)          /* enough source material? */
265	+	continue;
266	+	for (i = 0; i < pm->ncols; i++)
267	+	if (dst_rem[ rord[ri][1] = psortrow(pm,ri)[i] ] > minamt)
268	+	break;
269	+	if (i >= pm->ncols) {               /* moved all we can? */
270	+	free(rord);
271	+	return(.0);
272	+	}
273	+	}
274	+	if (pm->nrows > max2check)              /* sort if too many sources */
275	+	qsort_r(rord, pm->nrows, sizeof(short)*2, pm, &rmovcmp);
276	+	/* allocate cost array */
277	+	src_cost = (double *)malloc(sizeof(double)*pm->nrows);
278	+	if (src_cost == NULL)
279	+	goto memerr;
280	+	for (i = pm->nrows; i--; )              /* starting costs for diff. */
281	+	src_cost[i] = min_cost(src_rem[i], dst_rem, pm, i);
282	+	/* find best source & dest. */
283	+	best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
284	+	if ((r2check = pm->nrows) > max2check)
285	+	r2check = max2check;            /* put a limit on search */
286	+	for (ri = 0; ri < r2check; ri++) {      /* check each source row */
287	+	double      cost_others = 0;
288	+	cur.s = rord[ri][0];
289	+	if ((cur.d = rord[ri][1]) < 0 ||
290	+	(cur.price = pricerow(pm,cur.s)[cur.d]) >= best.price) {
291	+	if (pm->nrows > max2check) break;       /* sorted end */
292	+	continue;                       /* else skip this one */
293	+	}
294	+	cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
295	+	src_rem[cur.s] : dst_rem[cur.d];
296	+	/* don't just leave smidgen */
297	+	if (cur.amt > maxamt*1.02) cur.amt = maxamt;
298	+	dst_rem[cur.d] -= cur.amt;          /* add up opportunity costs */
299	+	for (i = pm->nrows; i--; )
300	+	if (i != cur.s)
301	+	cost_others += min_cost(src_rem[i], dst_rem, pm, i)
302	+	- src_cost[i];
303	+	dst_rem[cur.d] += cur.amt;          /* undo trial move */
304	+	cur.price += cost_others/cur.amt;   /* adjust effective price */
305	+	if (cur.price < best.price)         /* are we better than best? */
306	+	best = cur;
307	+	}
308	+	free(src_cost);                         /* clean up */
309	+	free(rord);
310	+	if ((best.s < 0) | (best.d < 0))        /* nothing left to move? */
311	+	return(.0);
312	+	/* else make the actual move */
313	+	mtx_coef(mig,best.s,best.d) += best.amt;
314	+	src_rem[best.s] -= best.amt;
315	+	dst_rem[best.d] -= best.amt;
316	+	return(best.amt);
317	+	memerr:
318	+	fprintf(stderr, "%s: Out of memory in migration_step()\n", progname);
319	+	exit(1);
320	+	}
321	+
322		/* Compute and insert migration along directed edge (may fork child) */
323		static MIGRATION *
324		create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
325		{
326	<	const double    end_thresh = 0.1/(from_rbf->nrbf*to_rbf->nrbf);
327	<	const double    check_thresh = 0.01;
295	<	const double    rel_thresh = 5e-6;
296	<	float           *pmtx;
326	>	const double    end_thresh = 5e-6;
327	>	PRICEMAT        pmtx;
328		MIGRATION       *newmig;
329		double          *src_rem, *dst_rem;
330		double          total_rem = 1., move_amt;
331	<	int             i;
331	>	int             i, j;
332		/* check if exists already */
333		for (newmig = from_rbf->ejl; newmig != NULL;
334		newmig = nextedge(from_rbf,newmig))
335		if (newmig->rbfv[1] == to_rbf)
336		return(NULL);
337		/* else allocate */
338	+	#ifdef DEBUG
339	+	fprintf(stderr, "Building path from (theta,phi) (%.1f,%.1f) ",
340	+	get_theta180(from_rbf->invec),
341	+	get_phi360(from_rbf->invec));
342	+	fprintf(stderr, "to (%.1f,%.1f) with %d x %d matrix\n",
343	+	get_theta180(to_rbf->invec),
344	+	get_phi360(to_rbf->invec),
345	+	from_rbf->nrbf, to_rbf->nrbf);
346	+	#endif
347		newmig = new_migration(from_rbf, to_rbf);
348		if (run_subprocess())
349		return(newmig);                 /* child continues */
350	<	pmtx = price_routes(from_rbf, to_rbf);
350	>	price_routes(&pmtx, from_rbf, to_rbf);
351		src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
352		dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
353		if ((src_rem == NULL) | (dst_rem == NULL)) {
#	Line 315 | Line 355 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
355		progname);
356		exit(1);
357		}
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
358		/* starting quantities */
359		memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
360		for (i = from_rbf->nrbf; i--; )
361		src_rem[i] = rbf_volume(&from_rbf->rbfa[i]) / from_rbf->vtotal;
362	<	for (i = to_rbf->nrbf; i--; )
363	<	dst_rem[i] = rbf_volume(&to_rbf->rbfa[i]) / to_rbf->vtotal;
362	>	for (j = to_rbf->nrbf; j--; )
363	>	dst_rem[j] = rbf_volume(&to_rbf->rbfa[j]) / to_rbf->vtotal;
364	>
365		do {                                    /* move a bit at a time */
366	<	move_amt = migration_step(newmig, src_rem, dst_rem, pmtx);
366	>	move_amt = migration_step(newmig, src_rem, dst_rem, &pmtx);
367		total_rem -= move_amt;
368	<	#ifdef DEBUG
369	<	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
368	>	} while ((total_rem > end_thresh) & (move_amt > 0));
369	>
370		for (i = from_rbf->nrbf; i--; ) {       /* normalize final matrix */
371	<	float       nf = rbf_volume(&from_rbf->rbfa[i]);
346	<	int         j;
371	>	double      nf = rbf_volume(&from_rbf->rbfa[i]);
372		if (nf <= FTINY) continue;
373		nf = from_rbf->vtotal / nf;
374		for (j = to_rbf->nrbf; j--; )
375	<	newmig->mtx[mtx_ndx(newmig,i,j)] *= nf;
375	>	mtx_coef(newmig,i,j) *= nf;     /* row now sums to 1.0 */
376		}
377		end_subprocess();                       /* exit here if subprocess */
378	<	free(pmtx);                             /* free working arrays */
378	>	free_routes(&pmtx);                     /* free working arrays */
379		free(src_rem);
380		free(dst_rem);
381		return(newmig);
#	Line 382 | Line 407 | overlaps_tri(const RBFNODE *bv0, const RBFNODE *bv1, c
407		return(vother[im_rev] != NULL);
408		}
409
410	<	/* Find context hull vertex to complete triangle (oriented call) */
410	>	/* Find convex hull vertex to complete triangle (oriented call) */
411		static RBFNODE *
412		find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rbf1)
413		{
#	Line 403 | Line 428 | find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rb
428		if (DOT(vp, vmid) <= FTINY)
429		continue;               /* wrong orientation */
430		area2 = .25*DOT(vp,vp);
431	<	VSUB(vp, rbf->invec, rbf0->invec);
431	>	VSUB(vp, rbf->invec, vmid);
432		dprod = -DOT(vp, vejn);
433		VSUM(vp, vp, vejn, dprod);      /* above guarantees non-zero */
434		dprod = DOT(vp, vmid) / VLEN(vp);
#	Line 459 | Line 484 | mesh_from_edge(MIGRATION *edge)
484		}
485		}
486		}
487	+
488	+	/* Add normal direction if missing */
489	+	static void
490	+	check_normal_incidence(void)
491	+	{
492	+	const int       saved_nprocs = nprocs;
493	+	RBFNODE         *near_rbf, *mir_rbf, *rbf;
494	+	double          bestd;
495	+	int             n, i, j;
496	+
497	+	if (dsf_list == NULL)
498	+	return;                         /* XXX should be error? */
499	+	near_rbf = dsf_list;
500	+	bestd = input_orient*near_rbf->invec[2];
501	+	if (single_plane_incident) {            /* ordered plane incidence? */
502	+	if (bestd >= 1.-2.*FTINY)
503	+	return;                 /* already have normal */
504	+	} else {
505	+	switch (inp_coverage) {
506	+	case INP_QUAD1:
507	+	case INP_QUAD2:
508	+	case INP_QUAD3:
509	+	case INP_QUAD4:
510	+	break;                  /* quadrilateral symmetry? */
511	+	default:
512	+	return;                 /* else we can interpolate */
513	+	}
514	+	for (rbf = near_rbf->next; rbf != NULL; rbf = rbf->next) {
515	+	const double    d = input_orient*rbf->invec[2];
516	+	if (d >= 1.-2.*FTINY)
517	+	return;         /* seems we have normal */
518	+	if (d > bestd) {
519	+	near_rbf = rbf;
520	+	bestd = d;
521	+	}
522	+	}
523	+	}
524	+	if (mig_list != NULL) {                 /* need to be called first */
525	+	fprintf(stderr, "%s: Late call to check_normal_incidence()\n",
526	+	progname);
527	+	exit(1);
528	+	}
529	+	#ifdef DEBUG
530	+	fprintf(stderr, "Interpolating normal incidence by mirroring (%.1f,%.1f)\n",
531	+	get_theta180(near_rbf->invec), get_phi360(near_rbf->invec));
532	+	#endif
533	+	/* mirror nearest incidence */
534	+	n = sizeof(RBFNODE) + sizeof(RBFVAL)*(near_rbf->nrbf-1);
535	+	mir_rbf = (RBFNODE *)malloc(n);
536	+	if (mir_rbf == NULL)
537	+	goto memerr;
538	+	memcpy(mir_rbf, near_rbf, n);
539	+	mir_rbf->ord = near_rbf->ord - 1;       /* not used, I think */
540	+	mir_rbf->next = NULL;
541	+	rev_rbf_symmetry(mir_rbf, MIRROR_X|MIRROR_Y);
542	+	nprocs = 1;                             /* compute migration matrix */
543	+	if (mig_list != create_migration(mir_rbf, near_rbf))
544	+	exit(1);                        /* XXX should never happen! */
545	+	n = 0;                                  /* count migrating particles */
546	+	for (i = 0; i < mtx_nrows(mig_list); i++)
547	+	for (j = 0; j < mtx_ncols(mig_list); j++)
548	+	n += (mtx_coef(mig_list,i,j) > FTINY);
549	+	rbf = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(n-1));
550	+	if (rbf == NULL)
551	+	goto memerr;
552	+	rbf->next = NULL; rbf->ejl = NULL;
553	+	rbf->invec[0] = rbf->invec[1] = 0; rbf->invec[2] = 1.;
554	+	rbf->nrbf = n;
555	+	rbf->vtotal = .5 + .5*mig_list->rbfv[1]->vtotal/mig_list->rbfv[0]->vtotal;
556	+	n = 0;                                  /* advect RBF lobes halfway */
557	+	for (i = 0; i < mtx_nrows(mig_list); i++) {
558	+	const RBFVAL        *rbf0i = &mig_list->rbfv[0]->rbfa[i];
559	+	const float         peak0 = rbf0i->peak;
560	+	const double        rad0 = R2ANG(rbf0i->crad);
561	+	FVECT               v0;
562	+	float               mv;
563	+	ovec_from_pos(v0, rbf0i->gx, rbf0i->gy);
564	+	for (j = 0; j < mtx_ncols(mig_list); j++)
565	+	if ((mv = mtx_coef(mig_list,i,j)) > FTINY) {
566	+	const RBFVAL    *rbf1j = &mig_list->rbfv[1]->rbfa[j];
567	+	double          rad2;
568	+	FVECT           v;
569	+	int             pos[2];
570	+	rad2 = R2ANG(rbf1j->crad);
571	+	rad2 = .5*(rad0*rad0 + rad2*rad2);
572	+	rbf->rbfa[n].peak = peak0 * mv * rbf->vtotal *
573	+	rad0*rad0/rad2;
574	+	rbf->rbfa[n].crad = ANG2R(sqrt(rad2));
575	+	ovec_from_pos(v, rbf1j->gx, rbf1j->gy);
576	+	geodesic(v, v0, v, .5, GEOD_REL);
577	+	pos_from_vec(pos, v);
578	+	rbf->rbfa[n].gx = pos[0];
579	+	rbf->rbfa[n].gy = pos[1];
580	+	++n;
581	+	}
582	+	}
583	+	rbf->vtotal *= mig_list->rbfv[0]->vtotal;
584	+	nprocs = saved_nprocs;                  /* final clean-up */
585	+	free(mir_rbf);
586	+	free(mig_list);
587	+	mig_list = near_rbf->ejl = NULL;
588	+	insert_dsf(rbf);                        /* insert interpolated normal */
589	+	return;
590	+	memerr:
591	+	fprintf(stderr, "%s: Out of memory in check_normal_incidence()\n",
592	+	progname);
593	+	exit(1);
594	+	}
595
596		/* Build our triangle mesh from recorded RBFs */
597		void
#	Line 467 | Line 600 | build_mesh(void)
600		double          best2 = M_PI*M_PI;
601		RBFNODE         *shrt_edj[2];
602		RBFNODE         *rbf0, *rbf1;
603	+	/* add normal if needed */
604	+	check_normal_incidence();
605		/* check if isotropic */
606		if (single_plane_incident) {
607		for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)

Diff Legend

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