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.21 by greg, Sun Mar 9 01:51:48 2014 UTC vs.
Revision 2.30 by greg, Thu Aug 21 10:33:48 2014 UTC

#	Line 27 | Line 27 | int                    nprocs = 1;
27		/* number of children (-1 in child) */
28		static int              nchild = 0;
29
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 */
30	>	/* Compute average DSF value at the given radius from central vector */
31	>	static double
32	>	eval_DSFsurround(const RBFNODE *rbf, const FVECT outvec, const double rad)
33	>	{
34	>	const int       ninc = 12;
35	>	const double    phinc = 2.*M_PI/ninc;
36	>	double          sum = 0;
37	>	int             n = 0;
38	>	FVECT           tvec;
39	>	int             i;
40	>	/* compute initial vector */
41	>	if (output_orient*outvec[2] >= 1.-FTINY) {
42	>	tvec[0] = tvec[2] = 0;
43	>	tvec[1] = 1;
44	>	} else {
45	>	tvec[0] = tvec[1] = 0;
46	>	tvec[2] = 1;
47	>	}
48	>	geodesic(tvec, outvec, tvec, rad, GEOD_RAD);
49	>	/* average surrounding DSF */
50	>	for (i = 0; i < ninc; i++) {
51	>	if (i) spinvector(tvec, tvec, outvec, phinc);
52	>	if (tvec[2] > 0 ^ output_orient > 0)
53	>	continue;
54	>	sum += eval_rbfrep(rbf, tvec) * output_orient*tvec[2];
55	>	++n;
56	>	}
57	>	if (n < 2)                              /* should never happen! */
58	>	return(sum);
59	>	return(sum/(double)n);
60	>	}
61
62	<	#define pricerow(p,i)   ((p)->price + (i)*(p)->ncols)
63	<	#define psortrow(p,i)   ((p)->sord + (i)*(p)->ncols)
62	>	/* Estimate single-lobe radius for DSF at the given outgoing angle */
63	>	static double
64	>	est_DSFrad(const RBFNODE *rbf, const FVECT outvec)
65	>	{
66	>	const double    rad_epsilon = 0.03;
67	>	const double    DSFtarget = 0.60653066 * eval_rbfrep(rbf,outvec)
68	>	* output_orient*outvec[2];
69	>	double          inside_rad = rad_epsilon;
70	>	double          outside_rad = 0.5;
71	>	double          DSFinside = eval_DSFsurround(rbf, outvec, inside_rad);
72	>	double          DSFoutside = eval_DSFsurround(rbf, outvec, outside_rad);
73	>	#define interp_rad      inside_rad + (outside_rad-inside_rad) * \
74	>	(DSFtarget-DSFinside) / (DSFoutside-DSFinside)
75	>	/* interpolation search */
76	>	while (outside_rad-inside_rad > rad_epsilon) {
77	>	double  test_rad = interp_rad;
78	>	double  DSFtest = eval_DSFsurround(rbf, outvec, test_rad);
79	>	if (DSFtarget < DSFtest) {
80	>	inside_rad = test_rad;
81	>	DSFinside = DSFtest;
82	>	} else {
83	>	outside_rad = test_rad;
84	>	DSFoutside = DSFtest;
85	>	}
86	>	}
87	>	return(interp_rad);
88	>	#undef interp_rad
89	>	}
90
91	+	/* Compute average BSDF peak from current DSF's */
92	+	static void
93	+	comp_bsdf_spec(void)
94	+	{
95	+	double          peak_sum = 0;
96	+	double          rad_sum = 0;
97	+	int             n = 0;
98	+	RBFNODE         *rbf;
99	+	FVECT           sdv;
100	+
101	+	if (dsf_list == NULL) {
102	+	bsdf_spec_peak = 0;
103	+	bsdf_spec_crad = 0;
104	+	return;
105	+	}
106	+	for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
107	+	sdv[0] = -rbf->invec[0];
108	+	sdv[1] = -rbf->invec[1];
109	+	sdv[2] = rbf->invec[2]*(2*(input_orient==output_orient) - 1);
110	+	peak_sum += eval_rbfrep(rbf, sdv);
111	+	rad_sum += est_DSFrad(rbf, sdv);
112	+	++n;
113	+	}
114	+	bsdf_spec_peak = peak_sum/(double)n;
115	+	bsdf_spec_crad = ANG2R( rad_sum/(double)n );
116	+	}
117	+
118		/* Create a new migration holder (sharing memory for multiprocessing) */
119		static MIGRATION *
120		new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
#	Line 165 | Line 243 | neighborhood_dist2(int x0, int y0, int x1, int y1)
243		double  d;
244		int     p[4];
245		int     i, j;
168	–
169	–	if ((x0 == x1) & (y0 == y1))
170	–	return(0.);
246		/* check radius */
247		p[0] = x0; p[1] = y0; p[2] = x1; p[3] = y1;
248		for (i = 4; i--; ) {
#	Line 183 | Line 258 | neighborhood_dist2(int x0, int y0, int x1, int y1)
258		return(sum2 / (4*rad*(rad+1) + 1));
259		}
260
261	<	/* Comparison routine needed for sorting price row */
262	<	static int
263	<	msrt_cmp(void *b, const void *p1, const void *p2)
261	>	/* Compute distance between two RBF lobes */
262	>	double
263	>	lobe_distance(RBFVAL *rbf1, RBFVAL *rbf2)
264		{
265	<	PRICEMAT        *pm = (PRICEMAT *)b;
266	<	float           c1 = pm->prow[*(const short *)p1];
267	<	float           c2 = pm->prow[*(const short *)p2];
268	<
269	<	if (c1 > c2) return(1);
270	<	if (c1 < c2) return(-1);
271	<	return(0);
265	>	FVECT   vfrom, vto;
266	>	double  d, res;
267	>	/* quadratic cost function */
268	>	ovec_from_pos(vfrom, rbf1->gx, rbf1->gy);
269	>	ovec_from_pos(vto, rbf2->gx, rbf2->gy);
270	>	d = Acos(DOT(vfrom, vto));
271	>	res = d*d;
272	>	d = R2ANG(rbf2->crad) - R2ANG(rbf1->crad);
273	>	res += d*d;
274	>	/* neighborhood difference */
275	>	res += NEIGH_FACT2 * neighborhood_dist2( rbf1->gx, rbf1->gy,
276	>	rbf2->gx, rbf2->gy );
277	>	return(res);
278		}
279
199	–	/* Compute (and allocate) migration price matrix for optimization */
200	–	static void
201	–	price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, const RBFNODE *to_rbf)
202	–	{
203	–	FVECT   *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
204	–	int     i, j;
280
206	–	compute_nDSFs(from_rbf, to_rbf);
207	–	pm->nrows = from_rbf->nrbf;
208	–	pm->ncols = to_rbf->nrbf;
209	–	pm->price = (float *)malloc(sizeof(float) * pm->nrows*pm->ncols);
210	–	pm->sord = (short *)malloc(sizeof(short) * pm->nrows*pm->ncols);
211	–
212	–	if ((pm->price == NULL) | (pm->sord == NULL) | (vto == NULL)) {
213	–	fprintf(stderr, "%s: Out of memory in migration_costs()\n",
214	–	progname);
215	–	exit(1);
216	–	}
217	–	for (j = to_rbf->nrbf; j--; )           /* save repetitive ops. */
218	–	ovec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);
219	–
220	–	for (i = from_rbf->nrbf; i--; ) {
221	–	const double        from_ang = R2ANG(from_rbf->rbfa[i].crad);
222	–	FVECT               vfrom;
223	–	short               *srow;
224	–	ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
225	–	pm->prow = pricerow(pm,i);
226	–	srow = psortrow(pm,i);
227	–	for (j = to_rbf->nrbf; j--; ) {
228	–	double  d;                      /* quadratic cost function */
229	–	d = Acos(DOT(vfrom, vto[j]));
230	–	pm->prow[j] = d*d;
231	–	d = R2ANG(to_rbf->rbfa[j].crad) - from_ang;
232	–	pm->prow[j] += d*d;
233	–	/* neighborhood difference */
234	–	pm->prow[j] += NEIGH_FACT2 * neighborhood_dist2(
235	–	from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy,
236	–	to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy );
237	–	srow[j] = j;
238	–	}
239	–	qsort_r(srow, pm->ncols, sizeof(short), pm, &msrt_cmp);
240	–	}
241	–	free(vto);
242	–	}
243	–
244	–	/* Free price matrix */
245	–	static void
246	–	free_routes(PRICEMAT *pm)
247	–	{
248	–	free(pm->price); pm->price = NULL;
249	–	free(pm->sord); pm->sord = NULL;
250	–	}
251	–
252	–	/* Compute minimum (optimistic) cost for moving the given source material */
253	–	static double
254	–	min_cost(double amt2move, const double *avail, const PRICEMAT *pm, int s)
255	–	{
256	–	const short     *srow = psortrow(pm,s);
257	–	const float     *prow = pricerow(pm,s);
258	–	double          total_cost = 0;
259	–	int             j;
260	–	/* move cheapest first */
261	–	for (j = 0; (j < pm->ncols) & (amt2move > FTINY); j++) {
262	–	int     d = srow[j];
263	–	double  amt = (amt2move < avail[d]) ? amt2move : avail[d];
264	–
265	–	total_cost += amt * prow[d];
266	–	amt2move -= amt;
267	–	}
268	–	return(total_cost);
269	–	}
270	–
271	–	typedef struct {
272	–	short   s, d;           /* source and destination */
273	–	float   dc;             /* discount to push inventory */
274	–	} ROWSENT;              /* row sort entry */
275	–
276	–	/* Compare entries by discounted moving price */
277	–	static int
278	–	rmovcmp(void *b, const void *p1, const void *p2)
279	–	{
280	–	PRICEMAT        *pm = (PRICEMAT *)b;
281	–	const ROWSENT   *re1 = (const ROWSENT *)p1;
282	–	const ROWSENT   *re2 = (const ROWSENT *)p2;
283	–	double          price_diff;
284	–
285	–	if (re1->d < 0) return(re2->d >= 0);
286	–	if (re2->d < 0) return(-1);
287	–	price_diff = re1->dc*pricerow(pm,re1->s)[re1->d] -
288	–	re2->dc*pricerow(pm,re2->s)[re2->d];
289	–	if (price_diff > 0) return(1);
290	–	if (price_diff < 0) return(-1);
291	–	return(0);
292	–	}
293	–
294	–	/* Take a step in migration by choosing reasonable bucket to transfer */
295	–	static double
296	–	migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, PRICEMAT *pm)
297	–	{
298	–	const int       max2check = 100;
299	–	const double    maxamt = 1./(double)pm->ncols;
300	–	const double    minamt = maxamt*1e-4;
301	–	double          *src_cost;
302	–	ROWSENT         *rord;
303	–	struct {
304	–	int     s, d;   /* source and destination */
305	–	double  price;  /* cost per amount moved */
306	–	double  amt;    /* amount we can move */
307	–	} cur, best;
308	–	int             r2check, i, ri;
309	–	/*
310	–	* Check cheapest available routes only -- a higher adjusted
311	–	* destination price implies that another source is closer, so
312	–	* we can hold off considering more expensive options until
313	–	* some other (hopefully better) moves have been made.
314	–	* A discount based on source remaining is supposed to prioritize
315	–	* movement from large lobes, but it doesn't seem to do much,
316	–	* so we have it set to 1.0 at the moment.
317	–	*/
318	–	#define discount(qr)    1.0
319	–	/* most promising row order */
320	–	rord = (ROWSENT *)malloc(sizeof(ROWSENT)*pm->nrows);
321	–	if (rord == NULL)
322	–	goto memerr;
323	–	for (ri = pm->nrows; ri--; ) {
324	–	rord[ri].s = ri;
325	–	rord[ri].d = -1;
326	–	rord[ri].dc = 1.f;
327	–	if (src_rem[ri] <= minamt)          /* enough source material? */
328	–	continue;
329	–	for (i = 0; i < pm->ncols; i++)
330	–	if (dst_rem[ rord[ri].d = psortrow(pm,ri)[i] ] > minamt)
331	–	break;
332	–	if (i >= pm->ncols) {               /* moved all we can? */
333	–	free(rord);
334	–	return(.0);
335	–	}
336	–	rord[ri].dc = discount(src_rem[ri]);
337	–	}
338	–	if (pm->nrows > max2check)              /* sort if too many sources */
339	–	qsort_r(rord, pm->nrows, sizeof(ROWSENT), pm, &rmovcmp);
340	–	/* allocate cost array */
341	–	src_cost = (double *)malloc(sizeof(double)*pm->nrows);
342	–	if (src_cost == NULL)
343	–	goto memerr;
344	–	for (i = pm->nrows; i--; )              /* starting costs for diff. */
345	–	src_cost[i] = min_cost(src_rem[i], dst_rem, pm, i);
346	–	/* find best source & dest. */
347	–	best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
348	–	if ((r2check = pm->nrows) > max2check)
349	–	r2check = max2check;            /* put a limit on search */
350	–	for (ri = 0; ri < r2check; ri++) {      /* check each source row */
351	–	double      cost_others = 0;
352	–	cur.s = rord[ri].s;
353	–	if ((cur.d = rord[ri].d) < 0 ||
354	–	rord[ri].dc*pricerow(pm,cur.s)[cur.d] >= best.price) {
355	–	if (pm->nrows > max2check) break;       /* sorted end */
356	–	continue;                       /* else skip this one */
357	–	}
358	–	cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
359	–	src_rem[cur.s] : dst_rem[cur.d];
360	–	/* don't just leave smidgen */
361	–	if (cur.amt > maxamt*1.02) cur.amt = maxamt;
362	–	dst_rem[cur.d] -= cur.amt;          /* add up opportunity costs */
363	–	for (i = pm->nrows; i--; )
364	–	if (i != cur.s)
365	–	cost_others += min_cost(src_rem[i], dst_rem, pm, i)
366	–	- src_cost[i];
367	–	dst_rem[cur.d] += cur.amt;          /* undo trial move */
368	–	/* discount effective price */
369	–	cur.price = ( pricerow(pm,cur.s)[cur.d] + cost_others/cur.amt ) *
370	–	rord[ri].dc;
371	–	if (cur.price < best.price)         /* are we better than best? */
372	–	best = cur;
373	–	}
374	–	free(src_cost);                         /* clean up */
375	–	free(rord);
376	–	if ((best.s < 0) | (best.d < 0))        /* nothing left to move? */
377	–	return(.0);
378	–	/* else make the actual move */
379	–	mtx_coef(mig,best.s,best.d) += best.amt;
380	–	src_rem[best.s] -= best.amt;
381	–	dst_rem[best.d] -= best.amt;
382	–	return(best.amt);
383	–	memerr:
384	–	fprintf(stderr, "%s: Out of memory in migration_step()\n", progname);
385	–	exit(1);
386	–	#undef discount
387	–	}
388	–
281		/* Compute and insert migration along directed edge (may fork child) */
282		static MIGRATION *
283		create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
284		{
393	–	const double    end_thresh = 5e-6;
394	–	PRICEMAT        pmtx;
285		MIGRATION       *newmig;
396	–	double          *src_rem, *dst_rem;
397	–	double          total_rem = 1., move_amt;
286		int             i, j;
287		/* check if exists already */
288		for (newmig = from_rbf->ejl; newmig != NULL;
#	Line 414 | Line 302 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
302		newmig = new_migration(from_rbf, to_rbf);
303		if (run_subprocess())
304		return(newmig);                 /* child continues */
417	–	price_routes(&pmtx, from_rbf, to_rbf);
418	–	src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
419	–	dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
420	–	if ((src_rem == NULL) | (dst_rem == NULL)) {
421	–	fprintf(stderr, "%s: Out of memory in create_migration()\n",
422	–	progname);
423	–	exit(1);
424	–	}
425	–	/* starting quantities */
426	–	memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
427	–	for (i = from_rbf->nrbf; i--; )
428	–	src_rem[i] = rbf_volume(&from_rbf->rbfa[i]) / from_rbf->vtotal;
429	–	for (j = to_rbf->nrbf; j--; )
430	–	dst_rem[j] = rbf_volume(&to_rbf->rbfa[j]) / to_rbf->vtotal;
305
306	<	do {                                    /* move a bit at a time */
307	<	move_amt = migration_step(newmig, src_rem, dst_rem, &pmtx);
308	<	total_rem -= move_amt;
435	<	} while ((total_rem > end_thresh) & (move_amt > 0));
306	>	/* compute transport plan */
307	>	compute_nDSFs(from_rbf, to_rbf);
308	>	plan_transport(newmig);
309
310		for (i = from_rbf->nrbf; i--; ) {       /* normalize final matrix */
311		double      nf = rbf_volume(&from_rbf->rbfa[i]);
#	Line 442 | Line 315 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
315		mtx_coef(newmig,i,j) *= nf;     /* row now sums to 1.0 */
316		}
317		end_subprocess();                       /* exit here if subprocess */
445	–	free_routes(&pmtx);                     /* free working arrays */
446	–	free(src_rem);
447	–	free(dst_rem);
318		return(newmig);
319		}
320
#	Line 534 | Line 404 | mesh_from_edge(MIGRATION *edge)
404		ej1 = create_migration(tvert[0], edge->rbfv[1]);
405		mesh_from_edge(ej0);
406		mesh_from_edge(ej1);
407	+	return;
408		}
409	<	} else if (tvert[1] == NULL) {          /* grow mesh on left */
409	>	}
410	>	if (tvert[1] == NULL) {                 /* grow mesh on left */
411		tvert[1] = find_chull_vert(edge->rbfv[1], edge->rbfv[0]);
412		if (tvert[1] != NULL) {
413		if (tvert[1]->ord > edge->rbfv[0]->ord)
#	Line 556 | Line 428 | mesh_from_edge(MIGRATION *edge)
428		static void
429		check_normal_incidence(void)
430		{
431	<	static const FVECT      norm_vec = {.0, .0, 1.};
431	>	static FVECT            norm_vec = {.0, .0, 1.};
432		const int               saved_nprocs = nprocs;
433		RBFNODE                 *near_rbf, *mir_rbf, *rbf;
434		double                  bestd;
#	Line 606 | Line 478 | check_normal_incidence(void)
478		memcpy(mir_rbf, near_rbf, n);
479		mir_rbf->ord = near_rbf->ord - 1;       /* not used, I think */
480		mir_rbf->next = NULL;
481	+	mir_rbf->ejl = NULL;
482		rev_rbf_symmetry(mir_rbf, MIRROR_X|MIRROR_Y);
483		nprocs = 1;                             /* compute migration matrix */
484	<	if (mig_list != create_migration(mir_rbf, near_rbf))
484	>	if (create_migration(mir_rbf, near_rbf) == NULL)
485		exit(1);                        /* XXX should never happen! */
486	<	/* interpolate normal dist. */
487	<	rbf = e_advect_rbf(mig_list, norm_vec, 2*near_rbf->nrbf);
486	>	norm_vec[2] = input_orient;             /* interpolate normal dist. */
487	>	rbf = e_advect_rbf(mig_list, norm_vec, 0);
488		nprocs = saved_nprocs;                  /* final clean-up */
489		free(mir_rbf);
490		free(mig_list);
#	Line 631 | Line 504 | build_mesh(void)
504		double          best2 = M_PI*M_PI;
505		RBFNODE         *shrt_edj[2];
506		RBFNODE         *rbf0, *rbf1;
507	+	/* average specular peak */
508	+	comp_bsdf_spec();
509		/* add normal if needed */
510		check_normal_incidence();
511		/* check if isotropic */

Diff Legend

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