ViewVC Help

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

root/radiance/ray/src/cv/pabopto2xml.c

Comparing ray/src/cv/pabopto2xml.c (file contents):
Revision 2.7 by greg, Sun Sep 2 15:33:15 2012 UTC vs.
Revision 2.15 by greg, Sun Sep 23 16:45:20 2012 UTC

#	Line 16 | Line 16 | static const char RCSid[] = "$Id$";
16		#include <math.h>
17		#include "bsdf.h"
18
19	+	#define DEBUG           1
20	+
21		#ifndef GRIDRES
22	<	#define GRIDRES         200             /* max. grid resolution per side */
22	>	#define GRIDRES         200             /* grid resolution per side */
23		#endif
24
25		#define RSCA            2.7             /* radius scaling factor (empirical) */
#	Line 51 | Line 53 | typedef struct s_rbfnode {
53		struct s_rbfnode        *next;          /* next in global RBF list */
54		MIGRATION               *ejl;           /* edge list for this vertex */
55		FVECT                   invec;          /* incident vector direction */
56	+	double                  vtotal;         /* volume for normalization */
57		int                     nrbf;           /* number of RBFs */
58		RBFVAL                  rbfa[1];        /* RBF array (extends struct) */
59	<	} RBFLIST;                      /* RBF representation of DSF @ 1 incidence */
59	>	} RBFNODE;                      /* RBF representation of DSF @ 1 incidence */
60
61		/* our loaded grid for this incident angle */
62	<	static double   theta_in_deg, phi_in_deg;
63	<	static GRIDVAL  dsf_grid[GRIDRES][GRIDRES];
62	>	static double           theta_in_deg, phi_in_deg;
63	>	static GRIDVAL          dsf_grid[GRIDRES][GRIDRES];
64
65	+	/* all incident angles in-plane so far? */
66	+	static int              single_plane_incident = -1;
67	+
68	+	/* input/output orientations */
69	+	static int              input_orient = 0;
70	+	static int              output_orient = 0;
71	+
72		/* processed incident DSF measurements */
73	<	static RBFLIST          *dsf_list = NULL;
73	>	static RBFNODE          *dsf_list = NULL;
74
75	<	/* edge (linking) matrices */
75	>	/* RBF-linking matrices (edges) */
76		static MIGRATION        *mig_list = NULL;
77
78	<	#define mtxval(m,i,j)   (m)->mtx[(i)*(m)->rbfv[1]->nrbf+(j)]
79	<	#define nextedge(rbf,m) (m)->enxt[(rbf)==(m)->rbfv[1]]
78	>	/* migration edges drawn in raster fashion */
79	>	static MIGRATION        *mig_grid[GRIDRES][GRIDRES];
80
81	+	#define mtx_nrows(m)    ((m)->rbfv[0]->nrbf)
82	+	#define mtx_ncols(m)    ((m)->rbfv[1]->nrbf)
83	+	#define mtx_ndx(m,i,j)  ((i)*mtx_ncols(m) + (j))
84	+	#define is_src(rbf,m)   ((rbf) == (m)->rbfv[0])
85	+	#define is_dest(rbf,m)  ((rbf) == (m)->rbfv[1])
86	+	#define nextedge(rbf,m) (m)->enxt[is_dest(rbf,m)]
87	+	#define opp_rbf(rbf,m)  (m)->rbfv[is_src(rbf,m)]
88	+
89	+	#define round(v)        (int)((v) + .5 - ((v) < -.5))
90	+
91	+	char                    *progname;
92	+
93	+	#ifdef DEBUG                    /* percentage to cull (<0 to turn off) */
94	+	int                     pctcull = -1;
95	+	#else
96	+	int                     pctcull = 90;
97	+	#endif
98	+	/* sampling order (set by data density) */
99	+	int                     samp_order = 0;
100	+
101	+	/* Compute volume associated with Gaussian lobe */
102	+	static double
103	+	rbf_volume(const RBFVAL *rbfp)
104	+	{
105	+	double  rad = R2ANG(rbfp->crad);
106	+
107	+	return((2.*M_PI) * rbfp->peak * rad*rad);
108	+	}
109	+
110		/* Compute outgoing vector from grid position */
111		static void
112	<	vec_from_pos(FVECT vec, int xpos, int ypos)
112	>	ovec_from_pos(FVECT vec, int xpos, int ypos)
113		{
114		double  uv[2];
115		double  r2;
#	Line 81 | Line 120 | vec_from_pos(FVECT vec, int xpos, int ypos)
120		vec[0] = vec[1] = sqrt(2. - r2);
121		vec[0] *= uv[0];
122		vec[1] *= uv[1];
123	<	vec[2] = 1. - r2;
123	>	vec[2] = output_orient*(1. - r2);
124		}
125
126	<	/* Compute grid position from normalized outgoing vector */
126	>	/* Compute grid position from normalized input/output vector */
127		static void
128		pos_from_vec(int pos[2], const FVECT vec)
129		{
130		double  sq[2];          /* uniform hemispherical projection */
131	<	double  norm = 1./sqrt(1. + vec[2]);
131	>	double  norm = 1./sqrt(1. + fabs(vec[2]));
132
133		SDdisk2square(sq, vec[0]*norm, vec[1]*norm);
134
#	Line 99 | Line 138 | pos_from_vec(int pos[2], const FVECT vec)
138
139		/* Evaluate RBF for DSF at the given normalized outgoing direction */
140		static double
141	<	eval_rbfrep(const RBFLIST *rp, const FVECT outvec)
141	>	eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
142		{
143		double          res = .0;
144		const RBFVAL    *rbfp;
#	Line 107 | Line 146 | eval_rbfrep(const RBFLIST *rp, const FVECT outvec)
146		double          sig2;
147		int             n;
148
149	+	if (rp == NULL)
150	+	return(.0);
151		rbfp = rp->rbfa;
152		for (n = rp->nrbf; n--; rbfp++) {
153	<	vec_from_pos(odir, rbfp->gx, rbfp->gy);
153	>	ovec_from_pos(odir, rbfp->gx, rbfp->gy);
154		sig2 = R2ANG(rbfp->crad);
155		sig2 = (DOT(odir,outvec) - 1.) / (sig2*sig2);
156		if (sig2 > -19.)
#	Line 118 | Line 159 | eval_rbfrep(const RBFLIST *rp, const FVECT outvec)
159		return(res);
160		}
161
162	+	/* Insert a new directional scattering function in our global list */
163	+	static void
164	+	insert_dsf(RBFNODE *newrbf)
165	+	{
166	+	RBFNODE         *rbf, *rbf_last;
167	+	/* check for redundant meas. */
168	+	for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
169	+	if (DOT(rbf->invec, newrbf->invec) >= 1.-FTINY) {
170	+	fputs("Duplicate incident measurement (ignored)\n", stderr);
171	+	free(newrbf);
172	+	return;
173	+	}
174	+	/* keep in ascending theta order */
175	+	for (rbf_last = NULL, rbf = dsf_list;
176	+	single_plane_incident & (rbf != NULL);
177	+	rbf_last = rbf, rbf = rbf->next)
178	+	if (input_orient*rbf->invec[2] < input_orient*newrbf->invec[2])
179	+	break;
180	+	if (rbf_last == NULL) {
181	+	newrbf->next = dsf_list;
182	+	dsf_list = newrbf;
183	+	return;
184	+	}
185	+	newrbf->next = rbf;
186	+	rbf_last->next = newrbf;
187	+	}
188	+
189		/* Count up filled nodes and build RBF representation from current grid */
190	<	static RBFLIST *
190	>	static RBFNODE *
191		make_rbfrep(void)
192		{
193		int     niter = 16;
194	+	int     minrad = ANG2R(pow(2., 1.-samp_order));
195		double  lastVar, thisVar = 100.;
196		int     nn;
197	<	RBFLIST *newnode;
197	>	RBFNODE *newnode;
198		int     i, j;
199
200		nn = 0;                 /* count selected bins */
#	Line 133 | Line 202 | make_rbfrep(void)
202		for (j = 0; j < GRIDRES; j++)
203		nn += dsf_grid[i][j].nval;
204		/* allocate RBF array */
205	<	newnode = (RBFLIST *)malloc(sizeof(RBFLIST) + sizeof(RBFVAL)*(nn-1));
205	>	newnode = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(nn-1));
206		if (newnode == NULL) {
207	<	fputs("Out of memory in make_rbfrep\n", stderr);
207	>	fputs("Out of memory in make_rbfrep()\n", stderr);
208		exit(1);
209		}
210		newnode->next = NULL;
#	Line 143 | Line 212 | make_rbfrep(void)
212		newnode->invec[2] = sin(M_PI/180.*theta_in_deg);
213		newnode->invec[0] = cos(M_PI/180.*phi_in_deg)*newnode->invec[2];
214		newnode->invec[1] = sin(M_PI/180.*phi_in_deg)*newnode->invec[2];
215	<	newnode->invec[2] = sqrt(1. - newnode->invec[2]*newnode->invec[2]);
215	>	newnode->invec[2] = input_orient*sqrt(1. - newnode->invec[2]*newnode->invec[2]);
216	>	newnode->vtotal = 0;
217		newnode->nrbf = nn;
218		nn = 0;                 /* fill RBF array */
219		for (i = 0; i < GRIDRES; i++)
#	Line 153 | Line 223 | make_rbfrep(void)
223		newnode->rbfa[nn].crad = RSCA*dsf_grid[i][j].crad + .5;
224		newnode->rbfa[nn].gx = i;
225		newnode->rbfa[nn].gy = j;
226	+	if (newnode->rbfa[nn].crad < minrad)
227	+	minrad = newnode->rbfa[nn].crad;
228		++nn;
229		}
230		/* iterate to improve interpolation accuracy */
231		do {
232	<	double  dsum = .0, dsum2 = .0;
232	>	double  dsum = 0, dsum2 = 0;
233		nn = 0;
234		for (i = 0; i < GRIDRES; i++)
235		for (j = 0; j < GRIDRES; j++)
236		if (dsf_grid[i][j].nval) {
237		FVECT   odir;
238		double  corr;
239	<	vec_from_pos(odir, i, j);
239	>	ovec_from_pos(odir, i, j);
240		newnode->rbfa[nn++].peak *= corr =
241		dsf_grid[i][j].vsum /
242		eval_rbfrep(newnode, odir);
#	Line 173 | Line 245 | make_rbfrep(void)
245		}
246		lastVar = thisVar;
247		thisVar = dsum2/(double)nn;
248	<	/*
248	>	#ifdef DEBUG
249		fprintf(stderr, "Avg., RMS error: %.1f%%  %.1f%%\n",
250		100.*dsum/(double)nn,
251		100.*sqrt(thisVar));
252	<	*/
252	>	#endif
253		} while (--niter > 0 && lastVar-thisVar > 0.02*lastVar);
254
255	<	newnode->next = dsf_list;
256	<	return(dsf_list = newnode);
255	>	nn = 0;                 /* compute sum for normalization */
256	>	while (nn < newnode->nrbf)
257	>	newnode->vtotal += rbf_volume(&newnode->rbfa[nn++]);
258	>
259	>	insert_dsf(newnode);
260	>	/* adjust sampling resolution */
261	>	samp_order = log(2./R2ANG(minrad))/M_LN2 + .5;
262	>
263	>	return(newnode);
264		}
265
266		/* Load a set of measurements corresponding to a particular incident angle */
267		static int
268	<	load_bsdf_meas(const char *fname)
268	>	load_pabopto_meas(const char *fname)
269		{
270		FILE    *fp = fopen(fname, "r");
271		int     inp_is_DSF = -1;
272	<	double  theta_out, phi_out, val;
272	>	double  new_phi, theta_out, phi_out, val;
273		char    buf[2048];
274		int     n, c;
275
#	Line 200 | Line 279 | load_bsdf_meas(const char *fname)
279		return(0);
280		}
281		memset(dsf_grid, 0, sizeof(dsf_grid));
282	+	#ifdef DEBUG
283	+	fprintf(stderr, "Loading measurement file '%s'...\n", fname);
284	+	#endif
285		/* read header information */
286		while ((c = getc(fp)) == '#' || c == EOF) {
287		if (fgets(buf, sizeof(buf), fp) == NULL) {
#	Line 218 | Line 300 | load_bsdf_meas(const char *fname)
300		}
301		if (sscanf(buf, "intheta %lf", &theta_in_deg) == 1)
302		continue;
303	<	if (sscanf(buf, "inphi %lf", &phi_in_deg) == 1)
303	>	if (sscanf(buf, "inphi %lf", &new_phi) == 1)
304		continue;
305		if (sscanf(buf, "incident_angle %lf %lf",
306	<	&theta_in_deg, &phi_in_deg) == 2)
306	>	&theta_in_deg, &new_phi) == 2)
307		continue;
308		}
309		if (inp_is_DSF < 0) {
#	Line 230 | Line 312 | load_bsdf_meas(const char *fname)
312		fclose(fp);
313		return(0);
314		}
315	<	ungetc(c, fp);          /* read actual data */
315	>	if (!input_orient)              /* check input orientation */
316	>	input_orient = 1 - 2*(theta_in_deg > 90.);
317	>	else if (input_orient > 0 ^ theta_in_deg < 90.) {
318	>	fputs("Cannot handle input angles on both sides of surface\n",
319	>	stderr);
320	>	exit(1);
321	>	}
322	>	if (single_plane_incident > 0)  /* check if still in plane */
323	>	single_plane_incident = (round(new_phi) == round(phi_in_deg));
324	>	else if (single_plane_incident < 0)
325	>	single_plane_incident = 1;
326	>	phi_in_deg = new_phi;
327	>	ungetc(c, fp);                  /* read actual data */
328		while (fscanf(fp, "%lf %lf %lf\n", &theta_out, &phi_out, &val) == 3) {
329		FVECT   ovec;
330		int     pos[2];
331
332	+	if (!output_orient)     /* check output orientation */
333	+	output_orient = 1 - 2*(theta_out > 90.);
334	+	else if (output_orient > 0 ^ theta_out < 90.) {
335	+	fputs("Cannot handle output angles on both sides of surface\n",
336	+	stderr);
337	+	exit(1);
338	+	}
339		ovec[2] = sin(M_PI/180.*theta_out);
340		ovec[0] = cos(M_PI/180.*phi_out) * ovec[2];
341		ovec[1] = sin(M_PI/180.*phi_out) * ovec[2];
#	Line 276 | Line 377 | compute_radii(void)
377		j = (i&1) ? jn : GRIDRES-1-jn;
378		if (dsf_grid[i][j].nval)        /* find empty grid pos. */
379		continue;
380	<	vec_from_pos(ovec0, i, j);
380	>	ovec_from_pos(ovec0, i, j);
381		inear = jnear = -1;             /* find nearest non-empty */
382		lastang2 = M_PI*M_PI;
383		for (ii = i-r; ii <= i+r; ii++) {
#	Line 287 | Line 388 | compute_radii(void)
388		if (jj >= GRIDRES) break;
389		if (!dsf_grid[ii][jj].nval)
390		continue;
391	<	vec_from_pos(ovec1, ii, jj);
391	>	ovec_from_pos(ovec1, ii, jj);
392		ang2 = 2. - 2.*DOT(ovec0,ovec1);
393		if (ang2 >= lastang2)
394		continue;
#	Line 304 | Line 405 | compute_radii(void)
405		if (r > dsf_grid[inear][jnear].crad)
406		dsf_grid[inear][jnear].crad = r;
407		/* next search radius */
408	<	r = ang2*(2.*GRIDRES/M_PI) + 1;
408	>	r = ang2*(2.*GRIDRES/M_PI) + 3;
409		}
410		/* blur radii over hemisphere */
411		memset(fill_grid, 0, sizeof(fill_grid));
#	Line 348 | Line 449 | cull_values(void)
449		continue;
450		if (!dsf_grid[i][j].crad)
451		continue;               /* shouldn't happen */
452	<	vec_from_pos(ovec0, i, j);
452	>	ovec_from_pos(ovec0, i, j);
453		maxang = 2.*R2ANG(dsf_grid[i][j].crad);
454		if (maxang > ovec0[2])          /* clamp near horizon */
455		maxang = ovec0[2];
#	Line 364 | Line 465 | cull_values(void)
465		continue;
466		if ((ii == i) & (jj == j))
467		continue;       /* don't get self-absorbed */
468	<	vec_from_pos(ovec1, ii, jj);
468	>	ovec_from_pos(ovec1, ii, jj);
469		if (2. - 2.*DOT(ovec0,ovec1) >= maxang2)
470		continue;
471		/* absorb sum */
#	Line 385 | Line 486 | cull_values(void)
486		}
487		}
488
489	+	/* Compute (and allocate) migration price matrix for optimization */
490	+	static float *
491	+	price_routes(const RBFNODE *from_rbf, const RBFNODE *to_rbf)
492	+	{
493	+	float   *pmtx = (float *)malloc(sizeof(float) *
494	+	from_rbf->nrbf * to_rbf->nrbf);
495	+	FVECT   *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
496	+	int     i, j;
497
498	+	if ((pmtx == NULL) | (vto == NULL)) {
499	+	fputs("Out of memory in migration_costs()\n", stderr);
500	+	exit(1);
501	+	}
502	+	for (j = to_rbf->nrbf; j--; )           /* save repetitive ops. */
503	+	ovec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);
504	+
505	+	for (i = from_rbf->nrbf; i--; ) {
506	+	const double        from_ang = R2ANG(from_rbf->rbfa[i].crad);
507	+	FVECT               vfrom;
508	+	ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
509	+	for (j = to_rbf->nrbf; j--; )
510	+	pmtx[i*to_rbf->nrbf + j] = acos(DOT(vfrom, vto[j])) +
511	+	fabs(R2ANG(to_rbf->rbfa[j].crad) - from_ang);
512	+	}
513	+	free(vto);
514	+	return(pmtx);
515	+	}
516	+
517	+	/* Comparison routine needed for sorting price row */
518	+	static const float      *price_arr;
519	+	static int
520	+	msrt_cmp(const void *p1, const void *p2)
521	+	{
522	+	float   c1 = price_arr[*(const int *)p1];
523	+	float   c2 = price_arr[*(const int *)p2];
524	+
525	+	if (c1 > c2) return(1);
526	+	if (c1 < c2) return(-1);
527	+	return(0);
528	+	}
529	+
530	+	/* Compute minimum (optimistic) cost for moving the given source material */
531	+	static double
532	+	min_cost(double amt2move, const double *avail, const float *price, int n)
533	+	{
534	+	static int      *price_sort = NULL;
535	+	static int      n_alloc = 0;
536	+	double          total_cost = 0;
537	+	int             i;
538	+
539	+	if (amt2move <= FTINY)                  /* pre-emptive check */
540	+	return(0.);
541	+	if (n > n_alloc) {                      /* (re)allocate sort array */
542	+	if (n_alloc) free(price_sort);
543	+	price_sort = (int *)malloc(sizeof(int)*n);
544	+	if (price_sort == NULL) {
545	+	fputs("Out of memory in min_cost()\n", stderr);
546	+	exit(1);
547	+	}
548	+	n_alloc = n;
549	+	}
550	+	for (i = n; i--; )
551	+	price_sort[i] = i;
552	+	price_arr = price;
553	+	qsort(price_sort, n, sizeof(int), &msrt_cmp);
554	+	/* move cheapest first */
555	+	for (i = 0; i < n && amt2move > FTINY; i++) {
556	+	int     d = price_sort[i];
557	+	double  amt = (amt2move < avail[d]) ? amt2move : avail[d];
558	+
559	+	total_cost += amt * price[d];
560	+	amt2move -= amt;
561	+	}
562	+	return(total_cost);
563	+	}
564	+
565	+	/* Take a step in migration by choosing optimal bucket to transfer */
566	+	static double
567	+	migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, const float *pmtx)
568	+	{
569	+	static double   *src_cost = NULL;
570	+	int             n_alloc = 0;
571	+	const double    maxamt = .1; /* 2./(mtx_nrows(mig)*mtx_ncols(mig)); */
572	+	double          amt = 0;
573	+	struct {
574	+	int     s, d;   /* source and destination */
575	+	double  price;  /* price estimate per amount moved */
576	+	double  amt;    /* amount we can move */
577	+	} cur, best;
578	+	int             i;
579	+
580	+	if (mtx_nrows(mig) > n_alloc) {         /* allocate cost array */
581	+	if (n_alloc)
582	+	free(src_cost);
583	+	src_cost = (double *)malloc(sizeof(double)*mtx_nrows(mig));
584	+	if (src_cost == NULL) {
585	+	fputs("Out of memory in migration_step()\n", stderr);
586	+	exit(1);
587	+	}
588	+	n_alloc = mtx_nrows(mig);
589	+	}
590	+	for (i = mtx_nrows(mig); i--; )         /* starting costs for diff. */
591	+	src_cost[i] = min_cost(src_rem[i], dst_rem,
592	+	pmtx+i*mtx_ncols(mig), mtx_ncols(mig));
593	+
594	+	/* find best source & dest. */
595	+	best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
596	+	for (cur.s = mtx_nrows(mig); cur.s--; ) {
597	+	const float *price = pmtx + cur.s*mtx_ncols(mig);
598	+	double      cost_others = 0;
599	+	if (src_rem[cur.s] <= FTINY)
600	+	continue;
601	+	cur.d = -1;                         /* examine cheapest dest. */
602	+	for (i = mtx_ncols(mig); i--; )
603	+	if (dst_rem[i] > FTINY &&
604	+	(cur.d < 0 || price[i] < price[cur.d]))
605	+	cur.d = i;
606	+	if (cur.d < 0)
607	+	return(.0);
608	+	if ((cur.price = price[cur.d]) >= best.price)
609	+	continue;                   /* no point checking further */
610	+	cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
611	+	src_rem[cur.s] : dst_rem[cur.d];
612	+	if (cur.amt > maxamt) cur.amt = maxamt;
613	+	dst_rem[cur.d] -= cur.amt;          /* add up differential costs */
614	+	for (i = mtx_nrows(mig); i--; ) {
615	+	if (i == cur.s) continue;
616	+	cost_others += min_cost(src_rem[i], dst_rem, price, mtx_ncols(mig))
617	+	- src_cost[i];
618	+	}
619	+	dst_rem[cur.d] += cur.amt;          /* undo trial move */
620	+	cur.price += cost_others/cur.amt;   /* adjust effective price */
621	+	if (cur.price < best.price)         /* are we better than best? */
622	+	best = cur;
623	+	}
624	+	if ((best.s < 0) | (best.d < 0))
625	+	return(.0);
626	+	/* make the actual move */
627	+	mig->mtx[mtx_ndx(mig,best.s,best.d)] += best.amt;
628	+	src_rem[best.s] -= best.amt;
629	+	dst_rem[best.d] -= best.amt;
630	+	return(best.amt);
631	+	}
632	+
633	+	#ifdef DEBUG
634	+	static char *
635	+	thetaphi(const FVECT v)
636	+	{
637	+	static char     buf[128];
638	+	double          theta, phi;
639	+
640	+	theta = 180./M_PI*acos(v[2]);
641	+	phi = 180./M_PI*atan2(v[1],v[0]);
642	+	sprintf(buf, "(%.0f,%.0f)", theta, phi);
643	+
644	+	return(buf);
645	+	}
646	+	#endif
647	+
648	+	/* Compute (and insert) migration along directed edge */
649	+	static MIGRATION *
650	+	make_migration(RBFNODE *from_rbf, RBFNODE *to_rbf, int creat_only)
651	+	{
652	+	const double    end_thresh = 0.02/(from_rbf->nrbf*to_rbf->nrbf);
653	+	float           *pmtx;
654	+	MIGRATION       *newmig;
655	+	double          *src_rem, *dst_rem;
656	+	double          total_rem = 1.;
657	+	int             i;
658	+	/* check if exists already */
659	+	for (newmig = from_rbf->ejl; newmig != NULL;
660	+	newmig = nextedge(from_rbf,newmig))
661	+	if (newmig->rbfv[1] == to_rbf)
662	+	return(creat_only ? (MIGRATION *)NULL : newmig);
663	+	/* else allocate */
664	+	pmtx = price_routes(from_rbf, to_rbf);
665	+	newmig = (MIGRATION *)malloc(sizeof(MIGRATION) + sizeof(float) *
666	+	(from_rbf->nrbf*to_rbf->nrbf - 1));
667	+	src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
668	+	dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
669	+	if ((newmig == NULL) | (src_rem == NULL) | (dst_rem == NULL)) {
670	+	fputs("Out of memory in make_migration()\n", stderr);
671	+	exit(1);
672	+	}
673	+	#ifdef DEBUG
674	+	{
675	+	fprintf(stderr, "Building path from (theta,phi) %s ",
676	+	thetaphi(from_rbf->invec));
677	+	fprintf(stderr, "to %s", thetaphi(to_rbf->invec));
678	+	}
679	+	#endif
680	+	newmig->next = NULL;
681	+	newmig->rbfv[0] = from_rbf;
682	+	newmig->rbfv[1] = to_rbf;
683	+	newmig->enxt[0] = newmig->enxt[1] = NULL;
684	+	memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
685	+	/* starting quantities */
686	+	for (i = from_rbf->nrbf; i--; )
687	+	src_rem[i] = rbf_volume(&from_rbf->rbfa[i]) / from_rbf->vtotal;
688	+	for (i = to_rbf->nrbf; i--; )
689	+	dst_rem[i] = rbf_volume(&to_rbf->rbfa[i]) / to_rbf->vtotal;
690	+	/* move a bit at a time */
691	+	while (total_rem > end_thresh) {
692	+	total_rem -= migration_step(newmig, src_rem, dst_rem, pmtx);
693	+	#ifdef DEBUG
694	+	/* fputc('.', stderr); */
695	+	fprintf(stderr, "\n%.9f remaining...", total_rem);
696	+	#endif
697	+	}
698	+	#ifdef DEBUG
699	+	fputs("done.\n", stderr);
700	+	#endif
701	+
702	+	free(pmtx);                             /* free working arrays */
703	+	free(src_rem);
704	+	free(dst_rem);
705	+	for (i = from_rbf->nrbf; i--; ) {       /* normalize final matrix */
706	+	float       nf = rbf_volume(&from_rbf->rbfa[i]);
707	+	int         j;
708	+	if (nf <= FTINY) continue;
709	+	nf = from_rbf->vtotal / nf;
710	+	for (j = to_rbf->nrbf; j--; )
711	+	newmig->mtx[mtx_ndx(newmig,i,j)] *= nf;
712	+	}
713	+	/* insert in edge lists */
714	+	newmig->enxt[0] = from_rbf->ejl;
715	+	from_rbf->ejl = newmig;
716	+	newmig->enxt[1] = to_rbf->ejl;
717	+	to_rbf->ejl = newmig;
718	+	newmig->next = mig_list;
719	+	return(mig_list = newmig);
720	+	}
721	+
722	+	/* Get triangle surface orientation (unnormalized) */
723	+	static void
724	+	tri_orient(FVECT vres, const FVECT v1, const FVECT v2, const FVECT v3)
725	+	{
726	+	FVECT   v2minus1, v3minus2;
727	+
728	+	VSUB(v2minus1, v2, v1);
729	+	VSUB(v3minus2, v3, v2);
730	+	VCROSS(vres, v2minus1, v3minus2);
731	+	}
732	+
733	+	/* Determine if vertex order is reversed (inward normal) */
734	+	static int
735	+	is_rev_tri(const FVECT v1, const FVECT v2, const FVECT v3)
736	+	{
737	+	FVECT   tor;
738	+
739	+	tri_orient(tor, v1, v2, v3);
740	+
741	+	return(DOT(tor, v2) < 0.);
742	+	}
743	+
744	+	/* Find vertices completing triangles on either side of the given edge */
745	+	static int
746	+	get_triangles(RBFNODE *rbfv[2], const MIGRATION *mig)
747	+	{
748	+	const MIGRATION *ej, *ej2;
749	+	RBFNODE         *tv;
750	+
751	+	rbfv[0] = rbfv[1] = NULL;
752	+	if (mig == NULL)
753	+	return(0);
754	+	for (ej = mig->rbfv[0]->ejl; ej != NULL;
755	+	ej = nextedge(mig->rbfv[0],ej)) {
756	+	if (ej == mig)
757	+	continue;
758	+	tv = opp_rbf(mig->rbfv[0],ej);
759	+	for (ej2 = tv->ejl; ej2 != NULL; ej2 = nextedge(tv,ej2))
760	+	if (opp_rbf(tv,ej2) == mig->rbfv[1]) {
761	+	rbfv[is_rev_tri(mig->rbfv[0]->invec,
762	+	mig->rbfv[1]->invec,
763	+	tv->invec)] = tv;
764	+	break;
765	+	}
766	+	}
767	+	return((rbfv[0] != NULL) + (rbfv[1] != NULL));
768	+	}
769	+
770	+	/* Check if prospective vertex would create overlapping triangle */
771	+	static int
772	+	overlaps_tri(const RBFNODE *bv0, const RBFNODE *bv1, const RBFNODE *pv)
773	+	{
774	+	const MIGRATION *ej;
775	+	RBFNODE         *vother[2];
776	+	int             im_rev;
777	+	/* find shared edge in mesh */
778	+	for (ej = pv->ejl; ej != NULL; ej = nextedge(pv,ej)) {
779	+	const RBFNODE   *tv = opp_rbf(pv,ej);
780	+	if (tv == bv0) {
781	+	im_rev = is_rev_tri(ej->rbfv[0]->invec,
782	+	ej->rbfv[1]->invec, bv1->invec);
783	+	break;
784	+	}
785	+	if (tv == bv1) {
786	+	im_rev = is_rev_tri(ej->rbfv[0]->invec,
787	+	ej->rbfv[1]->invec, bv0->invec);
788	+	break;
789	+	}
790	+	}
791	+	if (!get_triangles(vother, ej))         /* triangle on same side? */
792	+	return(0);
793	+	return(vother[im_rev] != NULL);
794	+	}
795	+
796	+	/* Find context hull vertex to complete triangle (oriented call) */
797	+	static RBFNODE *
798	+	find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rbf1)
799	+	{
800	+	FVECT   vmid, vejn, vp;
801	+	RBFNODE *rbf, *rbfbest = NULL;
802	+	double  dprod, area2, bestarea2 = FHUGE, bestdprod = -.5;
803	+
804	+	VSUB(vejn, rbf1->invec, rbf0->invec);
805	+	VADD(vmid, rbf0->invec, rbf1->invec);
806	+	if (normalize(vejn) == 0 || normalize(vmid) == 0)
807	+	return(NULL);
808	+	/* XXX exhaustive search */
809	+	/* Find triangle with minimum rotation from perpendicular */
810	+	for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
811	+	if ((rbf == rbf0) | (rbf == rbf1))
812	+	continue;
813	+	tri_orient(vp, rbf0->invec, rbf1->invec, rbf->invec);
814	+	if (DOT(vp, vmid) <= FTINY)
815	+	continue;               /* wrong orientation */
816	+	area2 = .25*DOT(vp,vp);
817	+	VSUB(vp, rbf->invec, rbf0->invec);
818	+	dprod = -DOT(vp, vejn);
819	+	VSUM(vp, vp, vejn, dprod);      /* above guarantees non-zero */
820	+	dprod = DOT(vp, vmid) / VLEN(vp);
821	+	if (dprod <= bestdprod + FTINY*(1 - 2*(area2 < bestarea2)))
822	+	continue;               /* found better already */
823	+	if (overlaps_tri(rbf0, rbf1, rbf))
824	+	continue;               /* overlaps another triangle */
825	+	rbfbest = rbf;
826	+	bestdprod = dprod;              /* new one to beat */
827	+	bestarea2 = area2;
828	+	}
829	+	return(rbfbest);
830	+	}
831	+
832	+	/* Create new migration edge and grow mesh recursively around it */
833	+	static void
834	+	mesh_from_edge(MIGRATION *edge)
835	+	{
836	+	MIGRATION       *ej0, *ej1;
837	+	RBFNODE         *tvert[2];
838	+
839	+	if (edge == NULL)
840	+	return;
841	+	/* triangle on either side? */
842	+	get_triangles(tvert, edge);
843	+	if (tvert[0] == NULL) {                 /* grow mesh on right */
844	+	tvert[0] = find_chull_vert(edge->rbfv[0], edge->rbfv[1]);
845	+	if (tvert[0] != NULL) {
846	+	if (tvert[0] > edge->rbfv[0])
847	+	ej0 = make_migration(edge->rbfv[0], tvert[0], 1);
848	+	else
849	+	ej0 = make_migration(tvert[0], edge->rbfv[0], 1);
850	+	if (tvert[0] > edge->rbfv[1])
851	+	ej1 = make_migration(edge->rbfv[1], tvert[0], 1);
852	+	else
853	+	ej1 = make_migration(tvert[0], edge->rbfv[1], 1);
854	+	mesh_from_edge(ej0);
855	+	mesh_from_edge(ej1);
856	+	}
857	+	} else if (tvert[1] == NULL) {          /* grow mesh on left */
858	+	tvert[1] = find_chull_vert(edge->rbfv[1], edge->rbfv[0]);
859	+	if (tvert[1] != NULL) {
860	+	if (tvert[1] > edge->rbfv[0])
861	+	ej0 = make_migration(edge->rbfv[0], tvert[1], 1);
862	+	else
863	+	ej0 = make_migration(tvert[1], edge->rbfv[0], 1);
864	+	if (tvert[1] > edge->rbfv[1])
865	+	ej1 = make_migration(edge->rbfv[1], tvert[1], 1);
866	+	else
867	+	ej1 = make_migration(tvert[1], edge->rbfv[1], 1);
868	+	mesh_from_edge(ej0);
869	+	mesh_from_edge(ej1);
870	+	}
871	+	}
872	+	}
873	+
874	+	#ifdef DEBUG
875	+	#include "random.h"
876	+	#include "bmpfile.h"
877	+	/* Hash pointer to byte value (must return 0 for NULL) */
878	+	static int
879	+	byte_hash(const void *p)
880	+	{
881	+	size_t  h = (size_t)p;
882	+	h ^= (size_t)p >> 8;
883	+	h ^= (size_t)p >> 16;
884	+	h ^= (size_t)p >> 24;
885	+	return(h & 0xff);
886	+	}
887	+	/* Write out BMP image showing edges */
888	+	static void
889	+	write_edge_image(const char *fname)
890	+	{
891	+	BMPHeader       *hdr = BMPmappedHeader(GRIDRES, GRIDRES, 0, 256);
892	+	BMPWriter       *wtr;
893	+	int             i, j;
894	+
895	+	fprintf(stderr, "Writing incident mesh drawing to '%s'\n", fname);
896	+	hdr->compr = BI_RLE8;
897	+	for (i = 256; --i; ) {                  /* assign random color map */
898	+	hdr->palette[i].r = random() & 0xff;
899	+	hdr->palette[i].g = random() & 0xff;
900	+	hdr->palette[i].b = random() & 0xff;
901	+	/* reject dark colors */
902	+	i += (hdr->palette[i].r + hdr->palette[i].g +
903	+	hdr->palette[i].b < 128);
904	+	}
905	+	hdr->palette[0].r = hdr->palette[0].g = hdr->palette[0].b = 0;
906	+	/* open output */
907	+	wtr = BMPopenOutputFile(fname, hdr);
908	+	if (wtr == NULL) {
909	+	free(hdr);
910	+	return;
911	+	}
912	+	for (i = 0; i < GRIDRES; i++) {         /* write scanlines */
913	+	for (j = 0; j < GRIDRES; j++)
914	+	wtr->scanline[j] = byte_hash(mig_grid[i][j]);
915	+	if (BMPwriteScanline(wtr) != BIR_OK)
916	+	break;
917	+	}
918	+	BMPcloseOutput(wtr);                    /* close & clean up */
919	+	}
920	+	#endif
921	+
922	+	/* Draw edge list into mig_grid array */
923	+	static void
924	+	draw_edges()
925	+	{
926	+	int             nnull = 0, ntot = 0;
927	+	MIGRATION       *ej;
928	+	int             p0[2], p1[2];
929	+
930	+	/* memset(mig_grid, 0, sizeof(mig_grid)); */
931	+	for (ej = mig_list; ej != NULL; ej = ej->next) {
932	+	++ntot;
933	+	pos_from_vec(p0, ej->rbfv[0]->invec);
934	+	pos_from_vec(p1, ej->rbfv[1]->invec);
935	+	if ((p0[0] == p1[0]) & (p0[1] == p1[1])) {
936	+	++nnull;
937	+	mig_grid[p0[0]][p0[1]] = ej;
938	+	continue;
939	+	}
940	+	if (abs(p1[0]-p0[0]) > abs(p1[1]-p0[1])) {
941	+	const int       xstep = 2*(p1[0] > p0[0]) - 1;
942	+	const double    ystep = (double)((p1[1]-p0[1])*xstep) /
943	+	(double)(p1[0]-p0[0]);
944	+	int             x;
945	+	double          y;
946	+	for (x = p0[0], y = p0[1]+.5; x != p1[0];
947	+	x += xstep, y += ystep)
948	+	mig_grid[x][(int)y] = ej;
949	+	mig_grid[x][(int)y] = ej;
950	+	} else {
951	+	const int       ystep = 2*(p1[1] > p0[1]) - 1;
952	+	const double    xstep = (double)((p1[0]-p0[0])*ystep) /
953	+	(double)(p1[1]-p0[1]);
954	+	int             y;
955	+	double          x;
956	+	for (y = p0[1], x = p0[0]+.5; y != p1[1];
957	+	y += ystep, x += xstep)
958	+	mig_grid[(int)x][y] = ej;
959	+	mig_grid[(int)x][y] = ej;
960	+	}
961	+	}
962	+	if (nnull)
963	+	fprintf(stderr, "Warning: %d of %d edges are null\n",
964	+	nnull, ntot);
965	+	#ifdef DEBUG
966	+	write_edge_image("bsdf_edges.bmp");
967	+	#endif
968	+	}
969	+
970	+	/* Build our triangle mesh from recorded RBFs */
971	+	static void
972	+	build_mesh()
973	+	{
974	+	double          best2 = M_PI*M_PI;
975	+	RBFNODE         *shrt_edj[2];
976	+	RBFNODE         *rbf0, *rbf1;
977	+	/* check if isotropic */
978	+	if (single_plane_incident) {
979	+	for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
980	+	if (rbf0->next != NULL)
981	+	make_migration(rbf0, rbf0->next, 1);
982	+	return;
983	+	}
984	+	/* start w/ shortest edge */
985	+	shrt_edj[0] = shrt_edj[1] = NULL;
986	+	for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
987	+	for (rbf1 = rbf0->next; rbf1 != NULL; rbf1 = rbf1->next) {
988	+	double  dist2 = 2. - 2.*DOT(rbf0->invec,rbf1->invec);
989	+	if (dist2 < best2) {
990	+	shrt_edj[0] = rbf0;
991	+	shrt_edj[1] = rbf1;
992	+	best2 = dist2;
993	+	}
994	+	}
995	+	if (shrt_edj[0] == NULL) {
996	+	fputs("Cannot find shortest edge\n", stderr);
997	+	exit(1);
998	+	}
999	+	/* build mesh from this edge */
1000	+	if (shrt_edj[0] < shrt_edj[1])
1001	+	mesh_from_edge(make_migration(shrt_edj[0], shrt_edj[1], 0));
1002	+	else
1003	+	mesh_from_edge(make_migration(shrt_edj[1], shrt_edj[0], 0));
1004	+	/* draw edge list into grid */
1005	+	draw_edges();
1006	+	}
1007	+
1008	+	/* Identify enclosing triangle for this position (flood fill raster check) */
1009	+	static int
1010	+	identify_tri(MIGRATION *miga[3], unsigned char vmap[GRIDRES][(GRIDRES+7)/8],
1011	+	int px, int py)
1012	+	{
1013	+	const int       btest = 1<<(py&07);
1014	+
1015	+	if (vmap[px][py>>3] & btest)            /* already visited here? */
1016	+	return(1);
1017	+	/* else mark it */
1018	+	vmap[px][py>>3] |= btest;
1019	+
1020	+	if (mig_grid[px][py] != NULL) {         /* are we on an edge? */
1021	+	int     i;
1022	+	for (i = 0; i < 3; i++) {
1023	+	if (miga[i] == mig_grid[px][py])
1024	+	return(1);
1025	+	if (miga[i] != NULL)
1026	+	continue;
1027	+	miga[i] = mig_grid[px][py];
1028	+	return(1);
1029	+	}
1030	+	return(0);                      /* outside triangle! */
1031	+	}
1032	+	/* check neighbors (flood) */
1033	+	if (px > 0 && !identify_tri(miga, vmap, px-1, py))
1034	+	return(0);
1035	+	if (px < GRIDRES-1 && !identify_tri(miga, vmap, px+1, py))
1036	+	return(0);
1037	+	if (py > 0 && !identify_tri(miga, vmap, px, py-1))
1038	+	return(0);
1039	+	if (py < GRIDRES-1 && !identify_tri(miga, vmap, px, py+1))
1040	+	return(0);
1041	+	return(1);                              /* this neighborhood done */
1042	+	}
1043	+
1044	+	/* Insert vertex in ordered list */
1045	+	static void
1046	+	insert_vert(RBFNODE **vlist, RBFNODE *v)
1047	+	{
1048	+	int     i, j;
1049	+
1050	+	for (i = 0; vlist[i] != NULL; i++) {
1051	+	if (v == vlist[i])
1052	+	return;
1053	+	if (v < vlist[i])
1054	+	break;
1055	+	}
1056	+	for (j = i; vlist[j] != NULL; j++)
1057	+	;
1058	+	while (j > i) {
1059	+	vlist[j] = vlist[j-1];
1060	+	--j;
1061	+	}
1062	+	vlist[i] = v;
1063	+	}
1064	+
1065	+	/* Sort triangle edges in standard order */
1066	+	static int
1067	+	order_triangle(MIGRATION *miga[3])
1068	+	{
1069	+	RBFNODE         *vert[7];
1070	+	MIGRATION       *ord[3];
1071	+	int             i;
1072	+	/* order vertices, first */
1073	+	memset(vert, 0, sizeof(vert));
1074	+	for (i = 3; i--; ) {
1075	+	if (miga[i] == NULL)
1076	+	return(0);
1077	+	insert_vert(vert, miga[i]->rbfv[0]);
1078	+	insert_vert(vert, miga[i]->rbfv[1]);
1079	+	}
1080	+	/* should be just 3 vertices */
1081	+	if ((vert[3] == NULL) | (vert[4] != NULL))
1082	+	return(0);
1083	+	/* identify edge 0 */
1084	+	for (i = 3; i--; )
1085	+	if (miga[i]->rbfv[0] == vert[0] &&
1086	+	miga[i]->rbfv[1] == vert[1]) {
1087	+	ord[0] = miga[i];
1088	+	break;
1089	+	}
1090	+	if (i < 0)
1091	+	return(0);
1092	+	/* identify edge 1 */
1093	+	for (i = 3; i--; )
1094	+	if (miga[i]->rbfv[0] == vert[1] &&
1095	+	miga[i]->rbfv[1] == vert[2]) {
1096	+	ord[1] = miga[i];
1097	+	break;
1098	+	}
1099	+	if (i < 0)
1100	+	return(0);
1101	+	/* identify edge 2 */
1102	+	for (i = 3; i--; )
1103	+	if (miga[i]->rbfv[0] == vert[0] &&
1104	+	miga[i]->rbfv[1] == vert[2]) {
1105	+	ord[2] = miga[i];
1106	+	break;
1107	+	}
1108	+	if (i < 0)
1109	+	return(0);
1110	+	/* reassign order */
1111	+	miga[0] = ord[0]; miga[1] = ord[1]; miga[2] = ord[2];
1112	+	return(1);
1113	+	}
1114	+
1115	+	/* Find edge(s) for interpolating the given incident vector */
1116	+	static int
1117	+	get_interp(MIGRATION *miga[3], const FVECT invec)
1118	+	{
1119	+	miga[0] = miga[1] = miga[2] = NULL;
1120	+	if (single_plane_incident) {            /* isotropic BSDF? */
1121	+	RBFNODE *rbf;                   /* find edge we're on */
1122	+	for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
1123	+	if (input_orient*rbf->invec[2] < input_orient*invec[2])
1124	+	break;
1125	+	if (rbf->next != NULL &&
1126	+	input_orient*rbf->next->invec[2] <
1127	+	input_orient*invec[2]) {
1128	+	for (miga[0] = rbf->ejl; miga[0] != NULL;
1129	+	miga[0] = nextedge(rbf,miga[0]))
1130	+	if (opp_rbf(rbf,miga[0]) == rbf->next)
1131	+	return(1);
1132	+	break;
1133	+	}
1134	+	}
1135	+	return(0);                      /* outside range! */
1136	+	}
1137	+	{                                       /* else use triangle mesh */
1138	+	unsigned char   floodmap[GRIDRES][(GRIDRES+7)/8];
1139	+	int             pstart[2];
1140	+	RBFNODE         *vother;
1141	+	MIGRATION       *ej;
1142	+	int             i;
1143	+
1144	+	pos_from_vec(pstart, invec);
1145	+	memset(floodmap, 0, sizeof(floodmap));
1146	+	/* call flooding function */
1147	+	if (!identify_tri(miga, floodmap, pstart[0], pstart[1]))
1148	+	return(0);              /* outside mesh */
1149	+	if ((miga[0] == NULL) | (miga[2] == NULL))
1150	+	return(0);              /* should never happen */
1151	+	if (miga[1] == NULL)
1152	+	return(1);              /* on edge */
1153	+	/* verify triangle */
1154	+	if (!order_triangle(miga)) {
1155	+	#ifdef DEBUG
1156	+	fputs("Munged triangle in get_interp()\n", stderr);
1157	+	#endif
1158	+	vother = NULL;          /* find triangle from edge */
1159	+	for (i = 3; i--; ) {
1160	+	RBFNODE     *tpair[2];
1161	+	if (get_triangles(tpair, miga[i]) &&
1162	+	(vother = tpair[ is_rev_tri(
1163	+	miga[i]->rbfv[0]->invec,
1164	+	miga[i]->rbfv[1]->invec,
1165	+	invec) ]) != NULL)
1166	+	break;
1167	+	}
1168	+	if (vother == NULL) {   /* couldn't find 3rd vertex */
1169	+	#ifdef DEBUG
1170	+	fputs("No triangle in get_interp()\n", stderr);
1171	+	#endif
1172	+	return(0);
1173	+	}
1174	+	/* reassign other two edges */
1175	+	for (ej = vother->ejl; ej != NULL;
1176	+	ej = nextedge(vother,ej)) {
1177	+	RBFNODE *vorig = opp_rbf(vother,ej);
1178	+	if (vorig == miga[i]->rbfv[0])
1179	+	miga[(i+1)%3] = ej;
1180	+	else if (vorig == miga[i]->rbfv[1])
1181	+	miga[(i+2)%3] = ej;
1182	+	}
1183	+	if (!order_triangle(miga)) {
1184	+	#ifdef DEBUG
1185	+	fputs("Bad triangle in get_interp()\n", stderr);
1186	+	#endif
1187	+	return(0);
1188	+	}
1189	+	}
1190	+	return(3);                      /* return in standard order */
1191	+	}
1192	+	}
1193	+
1194	+	/* Advect and allocate new RBF along edge */
1195	+	static RBFNODE *
1196	+	e_advect_rbf(const MIGRATION *mig, const FVECT invec)
1197	+	{
1198	+	RBFNODE         *rbf;
1199	+	int             n, i, j;
1200	+	double          t, full_dist;
1201	+	/* get relative position */
1202	+	t = acos(DOT(invec, mig->rbfv[0]->invec));
1203	+	if (t < M_PI/GRIDRES) {                 /* near first DSF */
1204	+	n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[0]->nrbf-1);
1205	+	rbf = (RBFNODE *)malloc(n);
1206	+	if (rbf == NULL)
1207	+	goto memerr;
1208	+	memcpy(rbf, mig->rbfv[0], n);   /* just duplicate */
1209	+	return(rbf);
1210	+	}
1211	+	full_dist = acos(DOT(mig->rbfv[0]->invec, mig->rbfv[1]->invec));
1212	+	if (t > full_dist-M_PI/GRIDRES) {       /* near second DSF */
1213	+	n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[1]->nrbf-1);
1214	+	rbf = (RBFNODE *)malloc(n);
1215	+	if (rbf == NULL)
1216	+	goto memerr;
1217	+	memcpy(rbf, mig->rbfv[1], n);   /* just duplicate */
1218	+	return(rbf);
1219	+	}
1220	+	t /= full_dist;
1221	+	n = 0;                                  /* count migrating particles */
1222	+	for (i = 0; i < mtx_nrows(mig); i++)
1223	+	for (j = 0; j < mtx_ncols(mig); j++)
1224	+	n += (mig->mtx[mtx_ndx(mig,i,j)] > FTINY);
1225	+	#ifdef DEBUG
1226	+	fprintf(stderr, "Input RBFs have %d, %d nodes -> output has %d\n",
1227	+	mig->rbfv[0]->nrbf, mig->rbfv[1]->nrbf, n);
1228	+	#endif
1229	+	rbf = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(n-1));
1230	+	if (rbf == NULL)
1231	+	goto memerr;
1232	+	rbf->next = NULL; rbf->ejl = NULL;
1233	+	VCOPY(rbf->invec, invec);
1234	+	rbf->nrbf = n;
1235	+	rbf->vtotal = 1.-t + t*mig->rbfv[1]->vtotal/mig->rbfv[0]->vtotal;
1236	+	n = 0;                                  /* advect RBF lobes */
1237	+	for (i = 0; i < mtx_nrows(mig); i++) {
1238	+	const RBFVAL        *rbf0i = &mig->rbfv[0]->rbfa[i];
1239	+	const float         peak0 = rbf0i->peak;
1240	+	const double        rad0 = R2ANG(rbf0i->crad);
1241	+	FVECT               v0;
1242	+	float               mv;
1243	+	ovec_from_pos(v0, rbf0i->gx, rbf0i->gy);
1244	+	for (j = 0; j < mtx_ncols(mig); j++)
1245	+	if ((mv = mig->mtx[mtx_ndx(mig,i,j)]) > FTINY) {
1246	+	const RBFVAL    *rbf1j = &mig->rbfv[1]->rbfa[j];
1247	+	double          rad1 = R2ANG(rbf1j->crad);
1248	+	FVECT           v;
1249	+	int             pos[2];
1250	+	rbf->rbfa[n].peak = peak0 * mv * rbf->vtotal;
1251	+	rbf->rbfa[n].crad = ANG2R(sqrt(rad0*rad0*(1.-t) +
1252	+	rad1*rad1*t));
1253	+	ovec_from_pos(v, rbf1j->gx, rbf1j->gy);
1254	+	geodesic(v, v0, v, t, GEOD_REL);
1255	+	pos_from_vec(pos, v);
1256	+	rbf->rbfa[n].gx = pos[0];
1257	+	rbf->rbfa[n].gy = pos[1];
1258	+	++n;
1259	+	}
1260	+	}
1261	+	rbf->vtotal *= mig->rbfv[0]->vtotal;    /* turn ratio into actual */
1262	+	return(rbf);
1263	+	memerr:
1264	+	fputs("Out of memory in e_advect_rbf()\n", stderr);
1265	+	exit(1);
1266	+	return(NULL);   /* pro forma return */
1267	+	}
1268	+
1269	+	/* Partially advect between recorded incident angles and allocate new RBF */
1270	+	static RBFNODE *
1271	+	advect_rbf(const FVECT invec)
1272	+	{
1273	+	MIGRATION       *miga[3];
1274	+	RBFNODE         *rbf;
1275	+	float           mbfact, mcfact;
1276	+	int             n, i, j, k;
1277	+	FVECT           v0, v1, v2;
1278	+	double          s, t;
1279	+
1280	+	if (!get_interp(miga, invec))           /* can't interpolate? */
1281	+	return(NULL);
1282	+	if (miga[1] == NULL)                    /* advect along edge? */
1283	+	return(e_advect_rbf(miga[0], invec));
1284	+	#ifdef DEBUG
1285	+	if (miga[0]->rbfv[0] != miga[2]->rbfv[0] |
1286	+	miga[0]->rbfv[1] != miga[1]->rbfv[0] |
1287	+	miga[1]->rbfv[1] != miga[2]->rbfv[1]) {
1288	+	fputs("Triangle vertex screw-up!\n", stderr);
1289	+	exit(1);
1290	+	}
1291	+	#endif
1292	+	/* figure out position */
1293	+	fcross(v0, miga[2]->rbfv[0]->invec, miga[2]->rbfv[1]->invec);
1294	+	normalize(v0);
1295	+	fcross(v2, miga[1]->rbfv[0]->invec, miga[1]->rbfv[1]->invec);
1296	+	normalize(v2);
1297	+	fcross(v1, invec, miga[1]->rbfv[1]->invec);
1298	+	normalize(v1);
1299	+	s = acos(DOT(v0,v1)) / acos(DOT(v0,v2));
1300	+	geodesic(v1, miga[0]->rbfv[0]->invec, miga[0]->rbfv[1]->invec,
1301	+	s, GEOD_REL);
1302	+	t = acos(DOT(v1,invec)) / acos(DOT(v1,miga[1]->rbfv[1]->invec));
1303	+	n = 0;                                  /* count migrating particles */
1304	+	for (i = 0; i < mtx_nrows(miga[0]); i++)
1305	+	for (j = 0; j < mtx_ncols(miga[0]); j++)
1306	+	for (k = (miga[0]->mtx[mtx_ndx(miga[0],i,j)] > FTINY) *
1307	+	mtx_ncols(miga[2]); k--; )
1308	+	n += (miga[2]->mtx[mtx_ndx(miga[2],i,k)] > FTINY &&
1309	+	miga[1]->mtx[mtx_ndx(miga[1],j,k)] > FTINY);
1310	+	#ifdef DEBUG
1311	+	fprintf(stderr, "Input RBFs have %d, %d, %d nodes -> output has %d\n",
1312	+	miga[0]->rbfv[0]->nrbf, miga[0]->rbfv[1]->nrbf,
1313	+	miga[2]->rbfv[1]->nrbf, n);
1314	+	#endif
1315	+	rbf = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(n-1));
1316	+	if (rbf == NULL) {
1317	+	fputs("Out of memory in advect_rbf()\n", stderr);
1318	+	exit(1);
1319	+	}
1320	+	rbf->next = NULL; rbf->ejl = NULL;
1321	+	VCOPY(rbf->invec, invec);
1322	+	rbf->nrbf = n;
1323	+	n = 0;                                  /* compute RBF lobes */
1324	+	mbfact = s * miga[0]->rbfv[1]->vtotal/miga[0]->rbfv[0]->vtotal *
1325	+	(1.-t + t*miga[1]->rbfv[1]->vtotal/miga[1]->rbfv[0]->vtotal);
1326	+	mcfact = (1.-s) *
1327	+	(1.-t + t*miga[2]->rbfv[1]->vtotal/miga[2]->rbfv[0]->vtotal);
1328	+	for (i = 0; i < mtx_nrows(miga[0]); i++) {
1329	+	const RBFVAL        *rbf0i = &miga[0]->rbfv[0]->rbfa[i];
1330	+	const float         w0i = rbf0i->peak;
1331	+	const double        rad0i = R2ANG(rbf0i->crad);
1332	+	ovec_from_pos(v0, rbf0i->gx, rbf0i->gy);
1333	+	for (j = 0; j < mtx_ncols(miga[0]); j++) {
1334	+	const float     ma = miga[0]->mtx[mtx_ndx(miga[0],i,j)];
1335	+	const RBFVAL    *rbf1j;
1336	+	double          rad1j, srad2;
1337	+	if (ma <= FTINY)
1338	+	continue;
1339	+	rbf1j = &miga[0]->rbfv[1]->rbfa[j];
1340	+	rad1j = R2ANG(rbf1j->crad);
1341	+	srad2 = (1.-s)*(1.-t)*rad0i*rad0i + s*(1.-t)*rad1j*rad1j;
1342	+	ovec_from_pos(v1, rbf1j->gx, rbf1j->gy);
1343	+	geodesic(v1, v0, v1, s, GEOD_REL);
1344	+	for (k = 0; k < mtx_ncols(miga[2]); k++) {
1345	+	float               mb = miga[1]->mtx[mtx_ndx(miga[1],j,k)];
1346	+	float               mc = miga[2]->mtx[mtx_ndx(miga[2],i,k)];
1347	+	const RBFVAL        *rbf2k;
1348	+	double              rad2k;
1349	+	FVECT               vout;
1350	+	int                 pos[2];
1351	+	if ((mb <= FTINY) | (mc <= FTINY))
1352	+	continue;
1353	+	rbf2k = &miga[2]->rbfv[1]->rbfa[k];
1354	+	rbf->rbfa[n].peak = w0i * ma * (mb*mbfact + mc*mcfact);
1355	+	rad2k = R2ANG(rbf2k->crad);
1356	+	rbf->rbfa[n].crad = ANG2R(sqrt(srad2 + t*rad2k*rad2k));
1357	+	ovec_from_pos(v2, rbf2k->gx, rbf2k->gy);
1358	+	geodesic(vout, v1, v2, t, GEOD_REL);
1359	+	pos_from_vec(pos, vout);
1360	+	rbf->rbfa[n].gx = pos[0];
1361	+	rbf->rbfa[n].gy = pos[1];
1362	+	++n;
1363	+	}
1364	+	}
1365	+	}
1366	+	rbf->vtotal = miga[0]->rbfv[0]->vtotal * (mbfact + mcfact);
1367	+	return(rbf);
1368	+	}
1369	+
1370	+	/* Interpolate and output isotropic BSDF data */
1371	+	static void
1372	+	interp_isotropic()
1373	+	{
1374	+	const int       sqres = 1<<samp_order;
1375	+	FILE            *ofp = NULL;
1376	+	char            cmd[128];
1377	+	int             ix, ox, oy;
1378	+	FVECT           ivec, ovec;
1379	+	double          bsdf;
1380	+	#if DEBUG
1381	+	fprintf(stderr, "Writing isotropic order %d ", samp_order);
1382	+	if (pctcull >= 0) fprintf(stderr, "data with %d%% culling\n", pctcull);
1383	+	else fputs("raw data\n", stderr);
1384	+	#endif
1385	+	if (pctcull >= 0) {                     /* begin output */
1386	+	sprintf(cmd, "rttree_reduce -h -a -fd -r 3 -t %d -g %d",
1387	+	pctcull, samp_order);
1388	+	fflush(stdout);
1389	+	ofp = popen(cmd, "w");
1390	+	if (ofp == NULL) {
1391	+	fputs("Cannot create pipe for rttree_reduce\n", stderr);
1392	+	exit(1);
1393	+	}
1394	+	} else
1395	+	fputs("{\n", stdout);
1396	+	/* run through directions */
1397	+	for (ix = 0; ix < sqres/2; ix++) {
1398	+	RBFNODE *rbf;
1399	+	SDsquare2disk(ivec, (ix+.5)/sqres, .5);
1400	+	ivec[2] = input_orient *
1401	+	sqrt(1. - ivec[0]*ivec[0] - ivec[1]*ivec[1]);
1402	+	rbf = advect_rbf(ivec);
1403	+	for (ox = 0; ox < sqres; ox++)
1404	+	for (oy = 0; oy < sqres; oy++) {
1405	+	SDsquare2disk(ovec, (ox+.5)/sqres, (oy+.5)/sqres);
1406	+	ovec[2] = output_orient *
1407	+	sqrt(1. - ovec[0]*ovec[0] - ovec[1]*ovec[1]);
1408	+	bsdf = eval_rbfrep(rbf, ovec) / fabs(ovec[2]);
1409	+	if (pctcull >= 0)
1410	+	fwrite(&bsdf, sizeof(bsdf), 1, ofp);
1411	+	else
1412	+	printf("\t%.3e\n", bsdf);
1413	+	}
1414	+	free(rbf);
1415	+	}
1416	+	if (pctcull >= 0) {                     /* finish output */
1417	+	if (pclose(ofp)) {
1418	+	fprintf(stderr, "Error running '%s'\n", cmd);
1419	+	exit(1);
1420	+	}
1421	+	} else {
1422	+	for (ix = sqres*sqres*sqres/2; ix--; )
1423	+	fputs("\t0\n", stdout);
1424	+	fputs("}\n", stdout);
1425	+	}
1426	+	}
1427	+
1428	+	/* Interpolate and output anisotropic BSDF data */
1429	+	static void
1430	+	interp_anisotropic()
1431	+	{
1432	+	const int       sqres = 1<<samp_order;
1433	+	FILE            *ofp = NULL;
1434	+	char            cmd[128];
1435	+	int             ix, iy, ox, oy;
1436	+	FVECT           ivec, ovec;
1437	+	double          bsdf;
1438	+	#if DEBUG
1439	+	fprintf(stderr, "Writing anisotropic order %d ", samp_order);
1440	+	if (pctcull >= 0) fprintf(stderr, "data with %d%% culling\n", pctcull);
1441	+	else fputs("raw data\n", stderr);
1442	+	#endif
1443	+	if (pctcull >= 0) {                     /* begin output */
1444	+	sprintf(cmd, "rttree_reduce -h -a -fd -r 4 -t %d -g %d",
1445	+	pctcull, samp_order);
1446	+	fflush(stdout);
1447	+	ofp = popen(cmd, "w");
1448	+	if (ofp == NULL) {
1449	+	fputs("Cannot create pipe for rttree_reduce\n", stderr);
1450	+	exit(1);
1451	+	}
1452	+	} else
1453	+	fputs("{\n", stdout);
1454	+	/* run through directions */
1455	+	for (ix = 0; ix < sqres; ix++)
1456	+	for (iy = 0; iy < sqres; iy++) {
1457	+	RBFNODE *rbf;
1458	+	SDsquare2disk(ivec, (ix+.5)/sqres, (iy+.5)/sqres);
1459	+	ivec[2] = input_orient *
1460	+	sqrt(1. - ivec[0]*ivec[0] - ivec[1]*ivec[1]);
1461	+	rbf = advect_rbf(ivec);
1462	+	for (ox = 0; ox < sqres; ox++)
1463	+	for (oy = 0; oy < sqres; oy++) {
1464	+	SDsquare2disk(ovec, (ox+.5)/sqres, (oy+.5)/sqres);
1465	+	ovec[2] = output_orient *
1466	+	sqrt(1. - ovec[0]*ovec[0] - ovec[1]*ovec[1]);
1467	+	bsdf = eval_rbfrep(rbf, ovec) / fabs(ovec[2]);
1468	+	if (pctcull >= 0)
1469	+	fwrite(&bsdf, sizeof(bsdf), 1, ofp);
1470	+	else
1471	+	printf("\t%.3e\n", bsdf);
1472	+	}
1473	+	free(rbf);
1474	+	}
1475	+	if (pctcull >= 0) {                     /* finish output */
1476	+	if (pclose(ofp)) {
1477	+	fprintf(stderr, "Error running '%s'\n", cmd);
1478	+	exit(1);
1479	+	}
1480	+	} else
1481	+	fputs("}\n", stdout);
1482	+	}
1483	+
1484		#if 1
1485	+	/* Read in BSDF files and interpolate as tensor tree representation */
1486	+	int
1487	+	main(int argc, char *argv[])
1488	+	{
1489	+	RBFNODE         *rbf;
1490	+	double          bsdf;
1491	+	int             i;
1492	+
1493	+	progname = argv[0];
1494	+	if (argc > 2 && !strcmp(argv[1], "-t")) {
1495	+	pctcull = atoi(argv[2]);
1496	+	argv += 2; argc -= 2;
1497	+	}
1498	+	if (argc < 3) {
1499	+	fprintf(stderr,
1500	+	"Usage: %s [-t pctcull] meas1.dat meas2.dat .. > bsdf.xml\n",
1501	+	progname);
1502	+	return(1);
1503	+	}
1504	+	for (i = 1; i < argc; i++) {            /* compile measurements */
1505	+	if (!load_pabopto_meas(argv[i]))
1506	+	return(1);
1507	+	compute_radii();
1508	+	cull_values();
1509	+	make_rbfrep();
1510	+	}
1511	+	build_mesh();                           /* create interpolation */
1512	+	/* xml_prologue();                              /* start XML output */
1513	+	if (single_plane_incident)              /* resample dist. */
1514	+	interp_isotropic();
1515	+	else
1516	+	interp_anisotropic();
1517	+	/* xml_epilogue();                              /* finish XML output */
1518	+	return(0);
1519	+	}
1520	+	#else
1521		/* Test main produces a Radiance model from the given input file */
1522		int
1523		main(int argc, char *argv[])
#	Line 401 | Line 1532 | main(int argc, char *argv[])
1532		fprintf(stderr, "Usage: %s input.dat > output.rad\n", argv[0]);
1533		return(1);
1534		}
1535	<	if (!load_bsdf_meas(argv[1]))
1535	>	if (!load_pabopto_meas(argv[1]))
1536		return(1);
1537
1538		compute_radii();
#	Line 414 | Line 1545 | main(int argc, char *argv[])
1545		for (i = 0; i < GRIDRES; i++)
1546		for (j = 0; j < GRIDRES; j++)
1547		if (dsf_grid[i][j].vsum > .0f) {
1548	<	vec_from_pos(dir, i, j);
1548	>	ovec_from_pos(dir, i, j);
1549		bsdf = dsf_grid[i][j].vsum / dir[2];
1550		if (dsf_grid[i][j].nval) {
1551		printf("pink cone c%04d\n0\n0\n8\n", ++n);
#	Line 425 | Line 1556 | main(int argc, char *argv[])
1556		dir[2]*(bsdf+.005));
1557		puts("\t.003\t0\n");
1558		} else {
1559	<	vec_from_pos(dir, i, j);
1559	>	ovec_from_pos(dir, i, j);
1560		printf("yellow sphere s%04d\n0\n0\n", ++n);
1561		printf("4 %.6g %.6g %.6g .0015\n\n",
1562		dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
#	Line 443 | Line 1574 | main(int argc, char *argv[])
1574		}
1575		for (i = 0; i < GRIDRES; i++)
1576		for (j = 0; j < GRIDRES; j++) {
1577	<	vec_from_pos(dir, i, j);
1577	>	ovec_from_pos(dir, i, j);
1578		bsdf = eval_rbfrep(dsf_list, dir) / dir[2];
1579		fprintf(pfp, "%.8e %.8e %.8e\n",
1580		dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines