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.2 by greg, Sat Oct 20 07:02:00 2012 UTC vs.
Revision 2.17 by greg, Wed Mar 5 22:47:16 2014 UTC

#	Line 23 | Line 23 | int                    nprocs = 1;
23		/* number of children (-1 in child) */
24		static int              nchild = 0;
25
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	+	#define pricerow(p,i)   ((p)->price + (i)*(p)->ncols)
34	+	#define psortrow(p,i)   ((p)->sord + (i)*(p)->ncols)
35	+
36		/* Create a new migration holder (sharing memory for multiprocessing) */
37		static MIGRATION *
38		new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
#	Line 110 | 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 124 | 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 float *
153	<	price_routes(const RBFNODE *from_rbf, const RBFNODE *to_rbf)
152	>	static void
153	>	price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, const RBFNODE *to_rbf)
154		{
131	–	float   *pmtx = (float *)malloc(sizeof(float) *
132	–	from_rbf->nrbf * to_rbf->nrbf);
155		FVECT   *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
156		int     i, j;
157
158	<	if ((pmtx == NULL) | (vto == NULL)) {
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);
#	Line 144 | Line 171 | price_routes(const RBFNODE *from_rbf, const RBFNODE *t
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	<	for (j = to_rbf->nrbf; j--; )
177	<	pmtx[i*to_rbf->nrbf + j] = acos(DOT(vfrom, vto[j])) +
178	<	fabs(R2ANG(to_rbf->rbfa[j].crad) - from_ang);
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);
153	–	return(pmtx);
190		}
191
192	<	/* Comparison routine needed for sorting price row */
193	<	static const float      *price_arr;
194	<	static int
159	<	msrt_cmp(const void *p1, const void *p2)
192	>	/* Free price matrix */
193	>	static void
194	>	free_routes(PRICEMAT *pm)
195		{
196	<	float   c1 = price_arr[*(const int *)p1];
197	<	float   c2 = price_arr[*(const int *)p2];
163	<
164	<	if (c1 > c2) return(1);
165	<	if (c1 < c2) return(-1);
166	<	return(0);
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 float *price, int n)
202	>	min_cost(double amt2move, const double *avail, const PRICEMAT *pm, int s)
203		{
204	<	static int      *price_sort = NULL;
205	<	static int      n_alloc = 0;
204	>	const short     *srow = psortrow(pm,s);
205	>	const float     *prow = pricerow(pm,s);
206		double          total_cost = 0;
207	<	int             i;
177	<
178	<	if (amt2move <= FTINY)                  /* pre-emptive check */
179	<	return(0.);
180	<	if (n > n_alloc) {                      /* (re)allocate sort array */
181	<	if (n_alloc) free(price_sort);
182	<	price_sort = (int *)malloc(sizeof(int)*n);
183	<	if (price_sort == NULL) {
184	<	fprintf(stderr, "%s: Out of memory in min_cost()\n",
185	<	progname);
186	<	exit(1);
187	<	}
188	<	n_alloc = n;
189	<	}
190	<	for (i = n; i--; )
191	<	price_sort[i] = i;
192	<	price_arr = price;
193	<	qsort(price_sort, n, sizeof(int), &msrt_cmp);
207	>	int             j;
208		/* move cheapest first */
209	<	for (i = 0; i < n && amt2move > FTINY; i++) {
210	<	int     d = price_sort[i];
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 * price[d];
213	>	total_cost += amt * prow[d];
214		amt2move -= amt;
215		}
216		return(total_cost);
217		}
218
219	<	/* Take a step in migration by choosing optimal bucket to transfer */
219	>	typedef struct {
220	>	short   s, d;           /* source and destination */
221	>	float   dc;             /* discount to push inventory */
222	>	} ROWSENT;              /* row sort entry */
223	>
224	>	/* Compare entries by discounted moving price */
225	>	static int
226	>	rmovcmp(void *b, const void *p1, const void *p2)
227	>	{
228	>	PRICEMAT        *pm = (PRICEMAT *)b;
229	>	const ROWSENT   *re1 = (const ROWSENT *)p1;
230	>	const ROWSENT   *re2 = (const ROWSENT *)p2;
231	>	double          price_diff;
232	>
233	>	if (re1->d < 0) return(re2->d >= 0);
234	>	if (re2->d < 0) return(-1);
235	>	price_diff = re1->dc*pricerow(pm,re1->s)[re1->d] -
236	>	re2->dc*pricerow(pm,re2->s)[re2->d];
237	>	if (price_diff > 0) return(1);
238	>	if (price_diff < 0) return(-1);
239	>	return(0);
240	>	}
241	>
242	>	/* Take a step in migration by choosing reasonable bucket to transfer */
243		static double
244	<	migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, const float *pmtx)
244	>	migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, PRICEMAT *pm)
245		{
246	<	const double    maxamt = .1;
247	<	const double    minamt = maxamt*5e-6;
248	<	static double   *src_cost = NULL;
249	<	static int      n_alloc = 0;
246	>	const int       max2check = 100;
247	>	const double    maxamt = 1./(double)pm->ncols;
248	>	const double    minamt = maxamt*1e-4;
249	>	double          *src_cost;
250	>	ROWSENT         *rord;
251		struct {
252		int     s, d;   /* source and destination */
253	<	double  price;  /* price estimate per amount moved */
253	>	double  price;  /* cost per amount moved */
254		double  amt;    /* amount we can move */
255		} cur, best;
256	<	int             i;
257	<
258	<	if (mtx_nrows(mig) > n_alloc) {         /* allocate cost array */
259	<	if (n_alloc)
260	<	free(src_cost);
261	<	src_cost = (double *)malloc(sizeof(double)*mtx_nrows(mig));
262	<	if (src_cost == NULL) {
263	<	fprintf(stderr, "%s: Out of memory in migration_step()\n",
264	<	progname);
265	<	exit(1);
266	<	}
267	<	n_alloc = mtx_nrows(mig);
256	>	int             r2check, i, ri;
257	>	/*
258	>	* Check cheapest available routes only -- a higher adjusted
259	>	* destination price implies that another source is closer, so
260	>	* we can hold off considering more expensive options until
261	>	* some other (hopefully better) moves have been made.
262	>	* A discount based on source remaining is supposed to prioritize
263	>	* movement from large lobes, but it doesn't seem to do much,
264	>	* so we have it set to 1.0 at the moment.
265	>	*/
266	>	#define discount(qr)    1.0
267	>	/* most promising row order */
268	>	rord = (ROWSENT *)malloc(sizeof(ROWSENT)*pm->nrows);
269	>	if (rord == NULL)
270	>	goto memerr;
271	>	for (ri = pm->nrows; ri--; ) {
272	>	rord[ri].s = ri;
273	>	rord[ri].d = -1;
274	>	rord[ri].dc = 1.f;
275	>	if (src_rem[ri] <= minamt)          /* enough source material? */
276	>	continue;
277	>	for (i = 0; i < pm->ncols; i++)
278	>	if (dst_rem[ rord[ri].d = psortrow(pm,ri)[i] ] > minamt)
279	>	break;
280	>	if (i >= pm->ncols) {               /* moved all we can? */
281	>	free(rord);
282	>	return(.0);
283	>	}
284	>	rord[ri].dc = discount(src_rem[ri]);
285		}
286	<	for (i = mtx_nrows(mig); i--; )         /* starting costs for diff. */
287	<	src_cost[i] = min_cost(src_rem[i], dst_rem,
288	<	pmtx+i*mtx_ncols(mig), mtx_ncols(mig));
289	<
286	>	if (pm->nrows > max2check)              /* sort if too many sources */
287	>	qsort_r(rord, pm->nrows, sizeof(ROWSENT), pm, &rmovcmp);
288	>	/* allocate cost array */
289	>	src_cost = (double *)malloc(sizeof(double)*pm->nrows);
290	>	if (src_cost == NULL)
291	>	goto memerr;
292	>	for (i = pm->nrows; i--; )              /* starting costs for diff. */
293	>	src_cost[i] = min_cost(src_rem[i], dst_rem, pm, i);
294		/* find best source & dest. */
295		best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
296	<	for (cur.s = mtx_nrows(mig); cur.s--; ) {
297	<	const float *price = pmtx + cur.s*mtx_ncols(mig);
296	>	if ((r2check = pm->nrows) > max2check)
297	>	r2check = max2check;            /* put a limit on search */
298	>	for (ri = 0; ri < r2check; ri++) {      /* check each source row */
299		double      cost_others = 0;
300	<	if (src_rem[cur.s] <= minamt)
301	<	continue;
302	<	cur.d = -1;                         /* examine cheapest dest. */
303	<	for (i = mtx_ncols(mig); i--; )
304	<	if (dst_rem[i] > minamt &&
305	<	(cur.d < 0 || price[i] < price[cur.d]))
246	<	cur.d = i;
247	<	if (cur.d < 0)
248	<	return(.0);
249	<	if ((cur.price = price[cur.d]) >= best.price)
250	<	continue;                   /* no point checking further */
300	>	cur.s = rord[ri].s;
301	>	if ((cur.d = rord[ri].d) < 0 ||
302	>	rord[ri].dc*pricerow(pm,cur.s)[cur.d] >= best.price) {
303	>	if (pm->nrows > max2check) break;       /* sorted end */
304	>	continue;                       /* else skip this one */
305	>	}
306		cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
307		src_rem[cur.s] : dst_rem[cur.d];
308	<	if (cur.amt > maxamt) cur.amt = maxamt;
309	<	dst_rem[cur.d] -= cur.amt;          /* add up differential costs */
310	<	for (i = mtx_nrows(mig); i--; )
308	>	/* don't just leave smidgen */
309	>	if (cur.amt > maxamt*1.02) cur.amt = maxamt;
310	>	dst_rem[cur.d] -= cur.amt;          /* add up opportunity costs */
311	>	for (i = pm->nrows; i--; )
312		if (i != cur.s)
313	<	cost_others += min_cost(src_rem[i], dst_rem,
258	<	price, mtx_ncols(mig))
313	>	cost_others += min_cost(src_rem[i], dst_rem, pm, i)
314		- src_cost[i];
315		dst_rem[cur.d] += cur.amt;          /* undo trial move */
316	<	cur.price += cost_others/cur.amt;   /* adjust effective price */
316	>	/* discount effective price */
317	>	cur.price = ( pricerow(pm,cur.s)[cur.d] + cost_others/cur.amt ) *
318	>	rord[ri].dc;
319		if (cur.price < best.price)         /* are we better than best? */
320	<	best = cur;
320	>	best = cur;
321		}
322	<	if ((best.s < 0) | (best.d < 0))
322	>	free(src_cost);                         /* clean up */
323	>	free(rord);
324	>	if ((best.s < 0) | (best.d < 0))        /* nothing left to move? */
325		return(.0);
326	<	/* make the actual move */
326	>	/* else make the actual move */
327		mtx_coef(mig,best.s,best.d) += best.amt;
328		src_rem[best.s] -= best.amt;
329		dst_rem[best.d] -= best.amt;
330		return(best.amt);
331	+	memerr:
332	+	fprintf(stderr, "%s: Out of memory in migration_step()\n", progname);
333	+	exit(1);
334	+	#undef discount
335		}
336
274	–	#ifdef DEBUG
275	–	static char *
276	–	thetaphi(const FVECT v)
277	–	{
278	–	static char     buf[128];
279	–	double          theta, phi;
280	–
281	–	theta = 180./M_PI*acos(v[2]);
282	–	phi = 180./M_PI*atan2(v[1],v[0]);
283	–	sprintf(buf, "(%.0f,%.0f)", theta, phi);
284	–
285	–	return(buf);
286	–	}
287	–	#endif
288	–
337		/* Compute and insert migration along directed edge (may fork child) */
338		static MIGRATION *
339		create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
340		{
341		const double    end_thresh = 5e-6;
342	<	float           *pmtx;
342	>	PRICEMAT        pmtx;
343		MIGRATION       *newmig;
344		double          *src_rem, *dst_rem;
345		double          total_rem = 1., move_amt;
346	<	int             i;
346	>	int             i, j;
347		/* check if exists already */
348		for (newmig = from_rbf->ejl; newmig != NULL;
349		newmig = nextedge(from_rbf,newmig))
350		if (newmig->rbfv[1] == to_rbf)
351		return(NULL);
352		/* else allocate */
353	+	#ifdef DEBUG
354	+	fprintf(stderr, "Building path from (theta,phi) (%.1f,%.1f) ",
355	+	get_theta180(from_rbf->invec),
356	+	get_phi360(from_rbf->invec));
357	+	fprintf(stderr, "to (%.1f,%.1f) with %d x %d matrix\n",
358	+	get_theta180(to_rbf->invec),
359	+	get_phi360(to_rbf->invec),
360	+	from_rbf->nrbf, to_rbf->nrbf);
361	+	#endif
362		newmig = new_migration(from_rbf, to_rbf);
363		if (run_subprocess())
364		return(newmig);                 /* child continues */
365	<	pmtx = price_routes(from_rbf, to_rbf);
365	>	price_routes(&pmtx, from_rbf, to_rbf);
366		src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
367		dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
368		if ((src_rem == NULL) | (dst_rem == NULL)) {
#	Line 313 | Line 370 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
370		progname);
371		exit(1);
372		}
316	–	#ifdef DEBUG
317	–	fprintf(stderr, "Building path from (theta,phi) %s ",
318	–	thetaphi(from_rbf->invec));
319	–	fprintf(stderr, "to %s with %d x %d matrix\n",
320	–	thetaphi(to_rbf->invec),
321	–	from_rbf->nrbf, to_rbf->nrbf);
322	–	#endif
373		/* starting quantities */
374		memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
375		for (i = from_rbf->nrbf; i--; )
376		src_rem[i] = rbf_volume(&from_rbf->rbfa[i]) / from_rbf->vtotal;
377	<	for (i = to_rbf->nrbf; i--; )
378	<	dst_rem[i] = rbf_volume(&to_rbf->rbfa[i]) / to_rbf->vtotal;
377	>	for (j = to_rbf->nrbf; j--; )
378	>	dst_rem[j] = rbf_volume(&to_rbf->rbfa[j]) / to_rbf->vtotal;
379	>
380		do {                                    /* move a bit at a time */
381	<	move_amt = migration_step(newmig, src_rem, dst_rem, pmtx);
381	>	move_amt = migration_step(newmig, src_rem, dst_rem, &pmtx);
382		total_rem -= move_amt;
332	–	#ifdef DEBUG
333	–	if (!nchild)
334	–	fprintf(stderr, "\r%.9f remaining...", total_rem);
335	–	#endif
383		} while ((total_rem > end_thresh) & (move_amt > 0));
384	<	#ifdef DEBUG
338	<	if (!nchild) fputs("done.\n", stderr);
339	<	else fprintf(stderr, "finished with %.9f remaining\n", total_rem);
340	<	#endif
384	>
385		for (i = from_rbf->nrbf; i--; ) {       /* normalize final matrix */
386	<	float       nf = rbf_volume(&from_rbf->rbfa[i]);
343	<	int         j;
386	>	double      nf = rbf_volume(&from_rbf->rbfa[i]);
387		if (nf <= FTINY) continue;
388		nf = from_rbf->vtotal / nf;
389		for (j = to_rbf->nrbf; j--; )
390	<	mtx_coef(newmig,i,j) *= nf;
390	>	mtx_coef(newmig,i,j) *= nf;     /* row now sums to 1.0 */
391		}
392		end_subprocess();                       /* exit here if subprocess */
393	<	free(pmtx);                             /* free working arrays */
393	>	free_routes(&pmtx);                     /* free working arrays */
394		free(src_rem);
395		free(dst_rem);
396		return(newmig);
#	Line 379 | Line 422 | overlaps_tri(const RBFNODE *bv0, const RBFNODE *bv1, c
422		return(vother[im_rev] != NULL);
423		}
424
425	<	/* Find context hull vertex to complete triangle (oriented call) */
425	>	/* Find convex hull vertex to complete triangle (oriented call) */
426		static RBFNODE *
427		find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rbf1)
428		{
#	Line 400 | Line 443 | find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rb
443		if (DOT(vp, vmid) <= FTINY)
444		continue;               /* wrong orientation */
445		area2 = .25*DOT(vp,vp);
446	<	VSUB(vp, rbf->invec, rbf0->invec);
446	>	VSUB(vp, rbf->invec, vmid);
447		dprod = -DOT(vp, vejn);
448		VSUM(vp, vp, vejn, dprod);      /* above guarantees non-zero */
449		dprod = DOT(vp, vmid) / VLEN(vp);
#	Line 456 | Line 499 | mesh_from_edge(MIGRATION *edge)
499		}
500		}
501		}
502	+
503	+	/* Add normal direction if missing */
504	+	static void
505	+	check_normal_incidence(void)
506	+	{
507	+	static const FVECT      norm_vec = {.0, .0, 1.};
508	+	const int               saved_nprocs = nprocs;
509	+	RBFNODE                 *near_rbf, *mir_rbf, *rbf;
510	+	double                  bestd;
511	+	int                     n;
512	+
513	+	if (dsf_list == NULL)
514	+	return;                         /* XXX should be error? */
515	+	near_rbf = dsf_list;
516	+	bestd = input_orient*near_rbf->invec[2];
517	+	if (single_plane_incident) {            /* ordered plane incidence? */
518	+	if (bestd >= 1.-2.*FTINY)
519	+	return;                 /* already have normal */
520	+	} else {
521	+	switch (inp_coverage) {
522	+	case INP_QUAD1:
523	+	case INP_QUAD2:
524	+	case INP_QUAD3:
525	+	case INP_QUAD4:
526	+	break;                  /* quadrilateral symmetry? */
527	+	default:
528	+	return;                 /* else we can interpolate */
529	+	}
530	+	for (rbf = near_rbf->next; rbf != NULL; rbf = rbf->next) {
531	+	const double    d = input_orient*rbf->invec[2];
532	+	if (d >= 1.-2.*FTINY)
533	+	return;         /* seems we have normal */
534	+	if (d > bestd) {
535	+	near_rbf = rbf;
536	+	bestd = d;
537	+	}
538	+	}
539	+	}
540	+	if (mig_list != NULL) {                 /* need to be called first */
541	+	fprintf(stderr, "%s: Late call to check_normal_incidence()\n",
542	+	progname);
543	+	exit(1);
544	+	}
545	+	#ifdef DEBUG
546	+	fprintf(stderr, "Interpolating normal incidence by mirroring (%.1f,%.1f)\n",
547	+	get_theta180(near_rbf->invec), get_phi360(near_rbf->invec));
548	+	#endif
549	+	/* mirror nearest incidence */
550	+	n = sizeof(RBFNODE) + sizeof(RBFVAL)*(near_rbf->nrbf-1);
551	+	mir_rbf = (RBFNODE *)malloc(n);
552	+	if (mir_rbf == NULL)
553	+	goto memerr;
554	+	memcpy(mir_rbf, near_rbf, n);
555	+	mir_rbf->ord = near_rbf->ord - 1;       /* not used, I think */
556	+	mir_rbf->next = NULL;
557	+	rev_rbf_symmetry(mir_rbf, MIRROR_X|MIRROR_Y);
558	+	nprocs = 1;                             /* compute migration matrix */
559	+	if (mig_list != create_migration(mir_rbf, near_rbf))
560	+	exit(1);                        /* XXX should never happen! */
561	+	/* interpolate normal dist. */
562	+	rbf = e_advect_rbf(mig_list, norm_vec, 2*near_rbf->nrbf);
563	+	nprocs = saved_nprocs;                  /* final clean-up */
564	+	free(mir_rbf);
565	+	free(mig_list);
566	+	mig_list = near_rbf->ejl = NULL;
567	+	insert_dsf(rbf);                        /* insert interpolated normal */
568	+	return;
569	+	memerr:
570	+	fprintf(stderr, "%s: Out of memory in check_normal_incidence()\n",
571	+	progname);
572	+	exit(1);
573	+	}
574
575		/* Build our triangle mesh from recorded RBFs */
576		void
#	Line 464 | Line 579 | build_mesh(void)
579		double          best2 = M_PI*M_PI;
580		RBFNODE         *shrt_edj[2];
581		RBFNODE         *rbf0, *rbf1;
582	+	/* add normal if needed */
583	+	check_normal_incidence();
584		/* check if isotropic */
585		if (single_plane_incident) {
586		for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)

Diff Legend

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