[ViewVC] Diff of: radiance/ray/src/cv/pabopto2xml.c

Comparing ray/src/cv/pabopto2xml.c (file contents):
Revision 2.8 by greg, Wed Sep 5 00:39:32 2012 UTC vs.
Revision 2.13 by greg, Fri Sep 21 05:17:22 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 54 \| Line 56 \| typedef struct s_rbfnode {
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		/* RBF-linking matrices (edges) */
76		static MIGRATION *mig_list = NULL;
77
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)
#	Line 84 \| Line 109 \| *rbf_volume(const RBFVAL rbfp)**
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 95 \| 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 113 \| 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 121 \| 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 132 \| 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	+	/* keep in ascending theta order */
168	+	for (rbf_last = NULL, rbf = dsf_list;
169	+	single_plane_incident & (rbf != NULL);
170	+	rbf_last = rbf, rbf = rbf->next)
171	+	if (input_orientrbf->invec[2] < input_orientnewrbf->invec[2])
172	+	break;
173	+	if (rbf_last == NULL) {
174	+	newrbf->next = dsf_list;
175	+	dsf_list = newrbf;
176	+	return;
177	+	}
178	+	newrbf->next = rbf;
179	+	rbf_last->next = newrbf;
180	+	}
181	+
182		/* Count up filled nodes and build RBF representation from current grid */
183	<	static RBFLIST *
183	>	static RBFNODE *
184		make_rbfrep(void)
185		{
186		int niter = 16;
187	+	int minrad = ANG2R(pow(2., 1.-samp_order));
188		double lastVar, thisVar = 100.;
189		int nn;
190	<	RBFLIST *newnode;
190	>	RBFNODE *newnode;
191		int i, j;
192
193		nn = 0; /* count selected bins */
#	Line 147 \| Line 195 \| make_rbfrep(void)
195		for (j = 0; j < GRIDRES; j++)
196		nn += dsf_grid[i][j].nval;
197		/* allocate RBF array */
198	<	newnode = (RBFLIST )malloc(sizeof(RBFLIST) + sizeof(RBFVAL)(nn-1));
198	>	newnode = (RBFNODE )malloc(sizeof(RBFNODE) + sizeof(RBFVAL)(nn-1));
199		if (newnode == NULL) {
200		fputs("Out of memory in make_rbfrep()\n", stderr);
201		exit(1);
#	Line 157 \| Line 205 \| make_rbfrep(void)
205		newnode->invec[2] = sin(M_PI/180.*theta_in_deg);
206		newnode->invec[0] = cos(M_PI/180.phi_in_deg)newnode->invec[2];
207		newnode->invec[1] = sin(M_PI/180.phi_in_deg)newnode->invec[2];
208	<	newnode->invec[2] = sqrt(1. - newnode->invec[2]*newnode->invec[2]);
208	>	newnode->invec[2] = input_orientsqrt(1. - newnode->invec[2]newnode->invec[2]);
209		newnode->vtotal = 0;
210		newnode->nrbf = nn;
211		nn = 0; /* fill RBF array */
#	Line 168 \| Line 216 \| make_rbfrep(void)
216		newnode->rbfa[nn].crad = RSCA*dsf_grid[i][j].crad + .5;
217		newnode->rbfa[nn].gx = i;
218		newnode->rbfa[nn].gy = j;
219	+	if (newnode->rbfa[nn].crad < minrad)
220	+	minrad = newnode->rbfa[nn].crad;
221		++nn;
222		}
223		/* iterate to improve interpolation accuracy */
224		do {
225	<	double dsum = .0, dsum2 = .0;
225	>	double dsum = 0, dsum2 = 0;
226		nn = 0;
227		for (i = 0; i < GRIDRES; i++)
228		for (j = 0; j < GRIDRES; j++)
229		if (dsf_grid[i][j].nval) {
230		FVECT odir;
231		double corr;
232	<	vec_from_pos(odir, i, j);
232	>	ovec_from_pos(odir, i, j);
233		newnode->rbfa[nn++].peak *= corr =
234		dsf_grid[i][j].vsum /
235		eval_rbfrep(newnode, odir);
#	Line 188 \| Line 238 \| make_rbfrep(void)
238		}
239		lastVar = thisVar;
240		thisVar = dsum2/(double)nn;
241	<	/*
241	>	#ifdef DEBUG
242		fprintf(stderr, "Avg., RMS error: %.1f%% %.1f%%\n",
243		100.*dsum/(double)nn,
244		100.*sqrt(thisVar));
245	<	*/
245	>	#endif
246		} while (--niter > 0 && lastVar-thisVar > 0.02*lastVar);
247
248		nn = 0; /* compute sum for normalization */
249		while (nn < newnode->nrbf)
250		newnode->vtotal += rbf_volume(&newnode->rbfa[nn++]);
251
252	<	newnode->next = dsf_list;
253	<	return(dsf_list = newnode);
252	>	insert_dsf(newnode);
253	>	/* adjust sampling resolution */
254	>	samp_order = log(2./R2ANG(minrad))/M_LN2 + .5;
255	>
256	>	return(newnode);
257		}
258
259		/* Load a set of measurements corresponding to a particular incident angle */
260		static int
261	<	load_bsdf_meas(const char *fname)
261	>	load_pabopto_meas(const char *fname)
262		{
263		FILE *fp = fopen(fname, "r");
264		int inp_is_DSF = -1;
265	<	double theta_out, phi_out, val;
265	>	double new_phi, theta_out, phi_out, val;
266		char buf[2048];
267		int n, c;
268
#	Line 219 \| Line 272 \| *load_bsdf_meas(const char fname)**
272		return(0);
273		}
274		memset(dsf_grid, 0, sizeof(dsf_grid));
275	+	#ifdef DEBUG
276	+	fprintf(stderr, "Loading measurement file '%s'...\n", fname);
277	+	#endif
278		/* read header information */
279		while ((c = getc(fp)) == '#' \|\| c == EOF) {
280		if (fgets(buf, sizeof(buf), fp) == NULL) {
#	Line 237 \| Line 293 \| *load_bsdf_meas(const char fname)**
293		}
294		if (sscanf(buf, "intheta %lf", &theta_in_deg) == 1)
295		continue;
296	<	if (sscanf(buf, "inphi %lf", &phi_in_deg) == 1)
296	>	if (sscanf(buf, "inphi %lf", &new_phi) == 1)
297		continue;
298		if (sscanf(buf, "incident_angle %lf %lf",
299	<	&theta_in_deg, &phi_in_deg) == 2)
299	>	&theta_in_deg, &new_phi) == 2)
300		continue;
301		}
302		if (inp_is_DSF < 0) {
#	Line 249 \| Line 305 \| *load_bsdf_meas(const char fname)**
305		fclose(fp);
306		return(0);
307		}
308	<	ungetc(c, fp); /* read actual data */
308	>	if (!input_orient) /* check input orientation */
309	>	input_orient = 1 - 2*(theta_in_deg > 90.);
310	>	else if (input_orient > 0 ^ theta_in_deg < 90.) {
311	>	fputs("Cannot handle input angles on both sides of surface\n",
312	>	stderr);
313	>	exit(1);
314	>	}
315	>	if (single_plane_incident > 0) /* check if still in plane */
316	>	single_plane_incident = (round(new_phi) == round(phi_in_deg));
317	>	else if (single_plane_incident < 0)
318	>	single_plane_incident = 1;
319	>	phi_in_deg = new_phi;
320	>	ungetc(c, fp); /* read actual data */
321		while (fscanf(fp, "%lf %lf %lf\n", &theta_out, &phi_out, &val) == 3) {
322		FVECT ovec;
323		int pos[2];
324
325	+	if (!output_orient) /* check output orientation */
326	+	output_orient = 1 - 2*(theta_out > 90.);
327	+	else if (output_orient > 0 ^ theta_out < 90.) {
328	+	fputs("Cannot handle output angles on both sides of surface\n",
329	+	stderr);
330	+	exit(1);
331	+	}
332		ovec[2] = sin(M_PI/180.*theta_out);
333		ovec[0] = cos(M_PI/180.phi_out) ovec[2];
334		ovec[1] = sin(M_PI/180.phi_out) ovec[2];
#	Line 295 \| Line 370 \| compute_radii(void)
370		j = (i&1) ? jn : GRIDRES-1-jn;
371		if (dsf_grid[i][j].nval) /* find empty grid pos. */
372		continue;
373	<	vec_from_pos(ovec0, i, j);
373	>	ovec_from_pos(ovec0, i, j);
374		inear = jnear = -1; /* find nearest non-empty */
375		lastang2 = M_PI*M_PI;
376		for (ii = i-r; ii <= i+r; ii++) {
#	Line 306 \| Line 381 \| compute_radii(void)
381		if (jj >= GRIDRES) break;
382		if (!dsf_grid[ii][jj].nval)
383		continue;
384	<	vec_from_pos(ovec1, ii, jj);
384	>	ovec_from_pos(ovec1, ii, jj);
385		ang2 = 2. - 2.*DOT(ovec0,ovec1);
386		if (ang2 >= lastang2)
387		continue;
#	Line 323 \| Line 398 \| compute_radii(void)
398		if (r > dsf_grid[inear][jnear].crad)
399		dsf_grid[inear][jnear].crad = r;
400		/* next search radius */
401	<	r = ang2(2.GRIDRES/M_PI) + 1;
401	>	r = ang2(2.GRIDRES/M_PI) + 3;
402		}
403		/* blur radii over hemisphere */
404		memset(fill_grid, 0, sizeof(fill_grid));
#	Line 367 \| Line 442 \| cull_values(void)
442		continue;
443		if (!dsf_grid[i][j].crad)
444		continue; /* shouldn't happen */
445	<	vec_from_pos(ovec0, i, j);
445	>	ovec_from_pos(ovec0, i, j);
446		maxang = 2.*R2ANG(dsf_grid[i][j].crad);
447		if (maxang > ovec0[2]) /* clamp near horizon */
448		maxang = ovec0[2];
#	Line 383 \| Line 458 \| cull_values(void)
458		continue;
459		if ((ii == i) & (jj == j))
460		continue; /* don't get self-absorbed */
461	<	vec_from_pos(ovec1, ii, jj);
461	>	ovec_from_pos(ovec1, ii, jj);
462		if (2. - 2.*DOT(ovec0,ovec1) >= maxang2)
463		continue;
464		/* absorb sum */
#	Line 406 \| Line 481 \| cull_values(void)
481
482		/* Compute (and allocate) migration price matrix for optimization */
483		static float *
484	<	price_routes(const RBFLIST from_rbf, const RBFLIST to_rbf)
484	>	price_routes(const RBFNODE from_rbf, const RBFNODE to_rbf)
485		{
486		float pmtx = (float )malloc(sizeof(float) *
487		from_rbf->nrbf * to_rbf->nrbf);
#	Line 418 \| Line 493 \| *price_routes(const RBFLIST from_rbf, const RBFLIST t*
493		exit(1);
494		}
495		for (j = to_rbf->nrbf; j--; ) /* save repetitive ops. */
496	<	vec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);
496	>	ovec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);
497
498		for (i = from_rbf->nrbf; i--; ) {
499		const double from_ang = R2ANG(from_rbf->rbfa[i].crad);
500		FVECT vfrom;
501	<	vec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
501	>	ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
502		for (j = to_rbf->nrbf; j--; )
503		pmtx[i*to_rbf->nrbf + j] = acos(DOT(vfrom, vto[j])) +
504		fabs(R2ANG(to_rbf->rbfa[j].crad) - from_ang);
#	Line 477 \| Line 552 \| *min_cost(double amt2move, const double avail, const f**
552		total_cost += amt * price[d];
553		amt2move -= amt;
554		}
480	–	if (amt2move > 1e-5) fprintf(stderr, "%g leftover!\n", amt2move);
555		return(total_cost);
556		}
557
#	Line 487 \| Line 561 \| *migration_step(MIGRATION mig, double src_rem, double*
561		{
562		static double *src_cost = NULL;
563		int n_alloc = 0;
564	<	const double maxamt = 0.5/(mtx_nrows(mig)*mtx_ncols(mig));
564	>	const double maxamt = 2./(mtx_nrows(mig)*mtx_ncols(mig));
565		double amt = 0;
566		struct {
567		int s, d; /* source and destination */
#	Line 551 \| Line 625 \| *migration_step(MIGRATION mig, double src_rem, double*
625
626		/* Compute (and insert) migration along directed edge */
627		static MIGRATION *
628	<	make_migration(RBFLIST from_rbf, RBFLIST to_rbf)
628	>	make_migration(RBFNODE from_rbf, RBFNODE to_rbf)
629		{
630		const double end_thresh = 0.02/(from_rbf->nrbf*to_rbf->nrbf);
631		float *pmtx = price_routes(from_rbf, to_rbf);
#	Line 567 \| Line 641 \| *make_migration(RBFLIST from_rbf, RBFLIST to_rbf)*
641		fputs("Out of memory in make_migration()\n", stderr);
642		exit(1);
643		}
644	+	#ifdef DEBUG
645	+	{
646	+	double theta, phi;
647	+	theta = acos(from_rbf->invec[2])*(180./M_PI);
648	+	phi = atan2(from_rbf->invec[1],from_rbf->invec[0])*(180./M_PI);
649	+	fprintf(stderr, "Building path from (theta,phi) (%d,%d) to ",
650	+	round(theta), round(phi));
651	+	theta = acos(to_rbf->invec[2])*(180./M_PI);
652	+	phi = atan2(to_rbf->invec[1],to_rbf->invec[0])*(180./M_PI);
653	+	fprintf(stderr, "(%d,%d)", round(theta), round(phi));
654	+	}
655	+	#endif
656		newmig->next = NULL;
657		newmig->rbfv[0] = from_rbf;
658		newmig->rbfv[1] = to_rbf;
#	Line 578 \| Line 664 \| *make_migration(RBFLIST from_rbf, RBFLIST to_rbf)*
664		for (i = to_rbf->nrbf; i--; )
665		dst_rem[i] = rbf_volume(&to_rbf->rbfa[i]) / to_rbf->vtotal;
666		/* move a bit at a time */
667	<	while (total_rem > end_thresh)
667	>	while (total_rem > end_thresh) {
668		total_rem -= migration_step(newmig, src_rem, dst_rem, pmtx);
669	+	#ifdef DEBUG
670	+	fputc('.', stderr);
671	+	/XXX/break;
672	+	#endif
673	+	}
674	+	#ifdef DEBUG
675	+	fputs("done.\n", stderr);
676	+	#endif
677
678		free(pmtx); /* free working arrays */
679		free(src_rem);
#	Line 601 \| Line 695 \| *make_migration(RBFLIST from_rbf, RBFLIST to_rbf)*
695		return(mig_list = newmig);
696		}
697
698	<	#if 0
699	<	/* Partially advect between the given RBFs to a newly allocated one */
700	<	static RBFLIST *
607	<	advect_rbf(const RBFLIST from_rbf, const RBFLIST to_rbf,
608	<	const float *mtx, const FVECT invec)
698	>	/* Get triangle surface orientation (unnormalized) */
699	>	static void
700	>	tri_orient(FVECT vres, const FVECT v1, const FVECT v2, const FVECT v3)
701		{
702	<	RBFLIST *rbf;
702	>	FVECT v2minus1, v3minus2;
703
704	<	if (from_rbf->nrbf > to_rbf->nrbf) {
705	<	fputs("Internal error: source RBF won't fit destination\n",
706	<	stderr);
704	>	VSUB(v2minus1, v2, v1);
705	>	VSUB(v3minus2, v3, v2);
706	>	VCROSS(vres, v2minus1, v3minus2);
707	>	}
708	>
709	>	/* Determine if vertex order is reversed (inward normal) */
710	>	static int
711	>	is_rev_tri(const FVECT v1, const FVECT v2, const FVECT v3)
712	>	{
713	>	FVECT tor;
714	>
715	>	tri_orient(tor, v1, v2, v3);
716	>
717	>	return(DOT(tor, v2) < 0.);
718	>	}
719	>
720	>	/* Find vertices completing triangles on either side of the given edge */
721	>	static int
722	>	get_triangles(RBFNODE rbfv[2], const MIGRATION mig)
723	>	{
724	>	const MIGRATION ej, ej2;
725	>	RBFNODE *tv;
726	>
727	>	rbfv[0] = rbfv[1] = NULL;
728	>	if (mig == NULL)
729	>	return(0);
730	>	for (ej = mig->rbfv[0]->ejl; ej != NULL;
731	>	ej = nextedge(mig->rbfv[0],ej)) {
732	>	if (ej == mig)
733	>	continue;
734	>	tv = opp_rbf(mig->rbfv[0],ej);
735	>	for (ej2 = tv->ejl; ej2 != NULL; ej2 = nextedge(tv,ej2))
736	>	if (opp_rbf(tv,ej2) == mig->rbfv[1]) {
737	>	rbfv[is_rev_tri(mig->rbfv[0]->invec,
738	>	mig->rbfv[1]->invec,
739	>	tv->invec)] = tv;
740	>	break;
741	>	}
742	>	}
743	>	return((rbfv[0] != NULL) + (rbfv[1] != NULL));
744	>	}
745	>
746	>	/* Check if prospective vertex would create overlapping triangle */
747	>	static int
748	>	overlaps_tri(const RBFNODE bv0, const RBFNODE bv1, const RBFNODE *pv)
749	>	{
750	>	const MIGRATION *ej;
751	>	RBFNODE *vother[2];
752	>	int im_rev;
753	>	/* find shared edge in mesh */
754	>	for (ej = pv->ejl; ej != NULL; ej = nextedge(pv,ej)) {
755	>	const RBFNODE *tv = opp_rbf(pv,ej);
756	>	if (tv == bv0) {
757	>	im_rev = is_rev_tri(ej->rbfv[0]->invec,
758	>	ej->rbfv[1]->invec, bv1->invec);
759	>	break;
760	>	}
761	>	if (tv == bv1) {
762	>	im_rev = is_rev_tri(ej->rbfv[0]->invec,
763	>	ej->rbfv[1]->invec, bv0->invec);
764	>	break;
765	>	}
766	>	}
767	>	if (!get_triangles(vother, ej))
768	>	return(0);
769	>	return(vother[im_rev] != NULL);
770	>	}
771	>
772	>	/* Find context hull vertex to complete triangle (oriented call) */
773	>	static RBFNODE *
774	>	find_chull_vert(const RBFNODE rbf0, const RBFNODE rbf1)
775	>	{
776	>	FVECT vmid, vor;
777	>	RBFNODE rbf, rbfbest = NULL;
778	>	double dprod2, area2, bestarea2 = FHUGE, bestdprod2 = 0.5;
779	>
780	>	VADD(vmid, rbf0->invec, rbf1->invec);
781	>	if (normalize(vmid) == 0)
782	>	return(NULL);
783	>	/* XXX exhaustive search */
784	>	for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
785	>	if ((rbf == rbf0) \| (rbf == rbf1))
786	>	continue;
787	>	tri_orient(vor, rbf0->invec, rbf1->invec, rbf->invec);
788	>	dprod2 = DOT(vor, vmid);
789	>	if (dprod2 <= FTINY)
790	>	continue; /* wrong orientation */
791	>	area2 = DOT(vor, vor);
792	>	dprod2 *= dprod2 / area2;
793	>	if (dprod2 > bestdprod2 +
794	>	FTINY(1 - 2(area2 < bestarea2)) &&
795	>	!overlaps_tri(rbf0, rbf1, rbf)) {
796	>	rbfbest = rbf;
797	>	bestdprod2 = dprod2;
798	>	bestarea2 = area2;
799	>	}
800	>	}
801	>	return(rbfbest);
802	>	}
803	>
804	>	/* Create new migration edge and grow mesh recursively around it */
805	>	static void
806	>	mesh_from_edge(MIGRATION *edge)
807	>	{
808	>	MIGRATION ej0, ej1;
809	>	RBFNODE *tvert[2];
810	>	/* triangle on either side? */
811	>	get_triangles(tvert, edge);
812	>	if (tvert[0] == NULL) { /* grow mesh on right */
813	>	tvert[0] = find_chull_vert(edge->rbfv[0], edge->rbfv[1]);
814	>	if (tvert[0] != NULL) {
815	>	if (tvert[0] > edge->rbfv[0])
816	>	ej0 = make_migration(edge->rbfv[0], tvert[0]);
817	>	else
818	>	ej0 = make_migration(tvert[0], edge->rbfv[0]);
819	>	if (tvert[0] > edge->rbfv[1])
820	>	ej1 = make_migration(edge->rbfv[1], tvert[0]);
821	>	else
822	>	ej1 = make_migration(tvert[0], edge->rbfv[1]);
823	>	mesh_from_edge(ej0);
824	>	mesh_from_edge(ej1);
825	>	}
826	>	}
827	>	if (tvert[1] == NULL) { /* grow mesh on left */
828	>	tvert[1] = find_chull_vert(edge->rbfv[1], edge->rbfv[0]);
829	>	if (tvert[1] != NULL) {
830	>	if (tvert[1] > edge->rbfv[0])
831	>	ej0 = make_migration(edge->rbfv[0], tvert[1]);
832	>	else
833	>	ej0 = make_migration(tvert[1], edge->rbfv[0]);
834	>	if (tvert[1] > edge->rbfv[1])
835	>	ej1 = make_migration(edge->rbfv[1], tvert[1]);
836	>	else
837	>	ej1 = make_migration(tvert[1], edge->rbfv[1]);
838	>	mesh_from_edge(ej0);
839	>	mesh_from_edge(ej1);
840	>	}
841	>	}
842	>	}
843	>
844	>	#ifdef DEBUG
845	>	#include "random.h"
846	>	#include "bmpfile.h"
847	>	/* Hash pointer to byte value */
848	>	static int
849	>	byte_hash(const void *p)
850	>	{
851	>	size_t h = (size_t)p;
852	>	h ^= (size_t)p >> 8;
853	>	h ^= (size_t)p >> 16;
854	>	h ^= (size_t)p >> 24;
855	>	return(h & 0xff);
856	>	}
857	>	/* Write out BMP image showing edges */
858	>	static void
859	>	write_edge_image(const char *fname)
860	>	{
861	>	BMPHeader *hdr = BMPmappedHeader(GRIDRES, GRIDRES, 0, 256);
862	>	BMPWriter *wtr;
863	>	int i, j;
864	>
865	>	fprintf(stderr, "Writing incident mesh drawing to '%s'\n", fname);
866	>	hdr->compr = BI_RLE8;
867	>	for (i = 256; --i; ) { /* assign random color map */
868	>	hdr->palette[i].r = random() & 0xff;
869	>	hdr->palette[i].r = random() & 0xff;
870	>	hdr->palette[i].r = random() & 0xff;
871	>	}
872	>	hdr->palette[0].r = hdr->palette[0].g = hdr->palette[0].b = 0;
873	>	/* open output */
874	>	wtr = BMPopenOutputFile(fname, hdr);
875	>	if (wtr == NULL) {
876	>	free(hdr);
877	>	return;
878	>	}
879	>	for (i = 0; i < GRIDRES; i++) { /* write scanlines */
880	>	for (j = 0; j < GRIDRES; j++)
881	>	wtr->scanline[j] = byte_hash(mig_grid[i][j]);
882	>	if (BMPwriteScanline(wtr) != BIR_OK)
883	>	break;
884	>	}
885	>	BMPcloseOutput(wtr); /* close & clean up */
886	>	}
887	>	#endif
888	>
889	>	/* Draw edge list into mig_grid array */
890	>	static void
891	>	draw_edges()
892	>	{
893	>	int nnull = 0, ntot = 0;
894	>	MIGRATION *ej;
895	>	int p0[2], p1[2];
896	>
897	>	/* memset(mig_grid, 0, sizeof(mig_grid)); */
898	>	for (ej = mig_list; ej != NULL; ej = ej->next) {
899	>	++ntot;
900	>	pos_from_vec(p0, ej->rbfv[0]->invec);
901	>	pos_from_vec(p1, ej->rbfv[1]->invec);
902	>	if ((p0[0] == p1[0]) & (p0[1] == p1[1])) {
903	>	++nnull;
904	>	mig_grid[p0[0]][p0[1]] = ej;
905	>	continue;
906	>	}
907	>	if (abs(p1[0]-p0[0]) > abs(p1[1]-p0[1])) {
908	>	const int xstep = 2*(p1[0] > p0[0]) - 1;
909	>	const double ystep = (double)((p1[1]-p0[1])*xstep) /
910	>	(double)(p1[0]-p0[0]);
911	>	int x;
912	>	double y;
913	>	for (x = p0[0], y = p0[1]+.5; x != p1[0];
914	>	x += xstep, y += ystep)
915	>	mig_grid[x][(int)y] = ej;
916	>	mig_grid[x][(int)y] = ej;
917	>	} else {
918	>	const int ystep = 2*(p1[1] > p0[1]) - 1;
919	>	const double xstep = (double)((p1[0]-p0[0])*ystep) /
920	>	(double)(p1[1]-p0[1]);
921	>	int y;
922	>	double x;
923	>	for (y = p0[1], x = p0[0]+.5; y != p1[1];
924	>	y += ystep, x += xstep)
925	>	mig_grid[(int)x][y] = ej;
926	>	mig_grid[(int)x][y] = ej;
927	>	}
928	>	}
929	>	if (nnull)
930	>	fprintf(stderr, "Warning: %d of %d edges are null\n",
931	>	nnull, ntot);
932	>	#ifdef DEBUG
933	>	write_edge_image("bsdf_edges.bmp");
934	>	#endif
935	>	}
936	>
937	>	/* Build our triangle mesh from recorded RBFs */
938	>	static void
939	>	build_mesh()
940	>	{
941	>	double best2 = M_PI*M_PI;
942	>	RBFNODE *shrt_edj[2];
943	>	RBFNODE rbf0, rbf1;
944	>	/* check if isotropic */
945	>	if (single_plane_incident) {
946	>	for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
947	>	if (rbf0->next != NULL)
948	>	make_migration(rbf0, rbf0->next);
949	>	return;
950	>	}
951	>	/* start w/ shortest edge */
952	>	shrt_edj[0] = shrt_edj[1] = NULL;
953	>	for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
954	>	for (rbf1 = rbf0->next; rbf1 != NULL; rbf1 = rbf1->next) {
955	>	double dist2 = 2. - 2.*DOT(rbf0->invec,rbf1->invec);
956	>	if (dist2 < best2) {
957	>	shrt_edj[0] = rbf0;
958	>	shrt_edj[1] = rbf1;
959	>	best2 = dist2;
960	>	}
961	>	}
962	>	if (shrt_edj[0] == NULL) {
963	>	fputs("Cannot find shortest edge\n", stderr);
964		exit(1);
965		}
966	<	rbf = (RBFLIST )malloc(sizeof(RBFLIST) + sizeof(RBFVAL)(to_rbf->nrbf-1));
966	>	/* build mesh from this edge */
967	>	if (shrt_edj[0] < shrt_edj[1])
968	>	mesh_from_edge(make_migration(shrt_edj[0], shrt_edj[1]));
969	>	else
970	>	mesh_from_edge(make_migration(shrt_edj[1], shrt_edj[0]));
971	>	/* draw edge list into grid */
972	>	draw_edges();
973	>	}
974	>
975	>	/* Identify enclosing triangle for this position (flood fill raster check) */
976	>	static int
977	>	identify_tri(MIGRATION *miga[3], unsigned char vmap[GRIDRES][(GRIDRES+7)/8],
978	>	int px, int py)
979	>	{
980	>	const int btest = 1<<(py&07);
981	>
982	>	if (vmap[px][py>>3] & btest) /* already visited here? */
983	>	return(1);
984	>	/* else mark it */
985	>	vmap[px][py>>3] \|= btest;
986	>
987	>	if (mig_grid[px][py] != NULL) { /* are we on an edge? */
988	>	int i;
989	>	for (i = 0; i < 3; i++) {
990	>	if (miga[i] == mig_grid[px][py])
991	>	return(1);
992	>	if (miga[i] != NULL)
993	>	continue;
994	>	miga[i] = mig_grid[px][py];
995	>	return(1);
996	>	}
997	>	return(0); /* outside triangle! */
998	>	}
999	>	/* check neighbors (flood) */
1000	>	if (px > 0 && !identify_tri(miga, vmap, px-1, py))
1001	>	return(0);
1002	>	if (px < GRIDRES-1 && !identify_tri(miga, vmap, px+1, py))
1003	>	return(0);
1004	>	if (py > 0 && !identify_tri(miga, vmap, px, py-1))
1005	>	return(0);
1006	>	if (py < GRIDRES-1 && !identify_tri(miga, vmap, px, py+1))
1007	>	return(0);
1008	>	return(1); /* this neighborhood done */
1009	>	}
1010	>
1011	>	/* Find edge(s) for interpolating the given incident vector */
1012	>	static int
1013	>	get_interp(MIGRATION *miga[3], const FVECT invec)
1014	>	{
1015	>	miga[0] = miga[1] = miga[2] = NULL;
1016	>	if (single_plane_incident) { /* isotropic BSDF? */
1017	>	RBFNODE rbf; / find edge we're on */
1018	>	for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
1019	>	if (input_orientrbf->invec[2] < input_orientinvec[2])
1020	>	break;
1021	>	if (rbf->next != NULL &&
1022	>	input_orient*rbf->next->invec[2] <
1023	>	input_orient*invec[2]) {
1024	>	for (miga[0] = rbf->ejl; miga[0] != NULL;
1025	>	miga[0] = nextedge(rbf,miga[0]))
1026	>	if (opp_rbf(rbf,miga[0]) == rbf->next)
1027	>	return(1);
1028	>	break;
1029	>	}
1030	>	}
1031	>	return(0); /* outside range! */
1032	>	}
1033	>	{ /* else use triangle mesh */
1034	>	unsigned char floodmap[GRIDRES][(GRIDRES+7)/8];
1035	>	int pstart[2];
1036	>
1037	>	pos_from_vec(pstart, invec);
1038	>	memset(floodmap, 0, sizeof(floodmap));
1039	>	/* call flooding function */
1040	>	if (!identify_tri(miga, floodmap, pstart[0], pstart[1]))
1041	>	return(0); /* outside mesh */
1042	>	if ((miga[0] == NULL) \| (miga[2] == NULL))
1043	>	return(0); /* should never happen */
1044	>	if (miga[1] == NULL)
1045	>	return(1); /* on edge */
1046	>	return(3); /* else in triangle */
1047	>	}
1048	>	}
1049	>
1050	>	/* Advect and allocate new RBF along edge */
1051	>	static RBFNODE *
1052	>	e_advect_rbf(const MIGRATION *mig, const FVECT invec)
1053	>	{
1054	>	RBFNODE *rbf;
1055	>	int n, i, j;
1056	>	double t, full_dist;
1057	>	/* get relative position */
1058	>	t = acos(DOT(invec, mig->rbfv[0]->invec));
1059	>	if (t < M_PI/GRIDRES) { /* near first DSF */
1060	>	n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[0]->nrbf-1);
1061	>	rbf = (RBFNODE *)malloc(n);
1062	>	if (rbf == NULL)
1063	>	goto memerr;
1064	>	memcpy(rbf, mig->rbfv[0], n); /* just duplicate */
1065	>	return(rbf);
1066	>	}
1067	>	full_dist = acos(DOT(mig->rbfv[0]->invec, mig->rbfv[1]->invec));
1068	>	if (t > full_dist-M_PI/GRIDRES) { /* near second DSF */
1069	>	n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[1]->nrbf-1);
1070	>	rbf = (RBFNODE *)malloc(n);
1071	>	if (rbf == NULL)
1072	>	goto memerr;
1073	>	memcpy(rbf, mig->rbfv[1], n); /* just duplicate */
1074	>	return(rbf);
1075	>	}
1076	>	t /= full_dist;
1077	>	n = 0; /* count migrating particles */
1078	>	for (i = 0; i < mtx_nrows(mig); i++)
1079	>	for (j = 0; j < mtx_ncols(mig); j++)
1080	>	n += (mig->mtx[mtx_ndx(mig,i,j)] > FTINY);
1081	>	#ifdef DEBUG
1082	>	fprintf(stderr, "Input RBFs have %d, %d nodes -> output has %d\n",
1083	>	mig->rbfv[0]->nrbf, mig->rbfv[1]->nrbf, n);
1084	>	#endif
1085	>	rbf = (RBFNODE )malloc(sizeof(RBFNODE) + sizeof(RBFVAL)(n-1));
1086	>	if (rbf == NULL)
1087	>	goto memerr;
1088	>	rbf->next = NULL; rbf->ejl = NULL;
1089	>	VCOPY(rbf->invec, invec);
1090	>	rbf->nrbf = n;
1091	>	rbf->vtotal = 1.-t + t*mig->rbfv[1]->vtotal/mig->rbfv[0]->vtotal;
1092	>	n = 0; /* advect RBF lobes */
1093	>	for (i = 0; i < mtx_nrows(mig); i++) {
1094	>	const RBFVAL *rbf0i = &mig->rbfv[0]->rbfa[i];
1095	>	const float peak0 = rbf0i->peak;
1096	>	const double rad0 = R2ANG(rbf0i->crad);
1097	>	FVECT v0;
1098	>	float mv;
1099	>	ovec_from_pos(v0, rbf0i->gx, rbf0i->gy);
1100	>	for (j = 0; j < mtx_ncols(mig); j++)
1101	>	if ((mv = mig->mtx[mtx_ndx(mig,i,j)]) > FTINY) {
1102	>	const RBFVAL *rbf1j = &mig->rbfv[1]->rbfa[j];
1103	>	double rad1 = R2ANG(rbf1j->crad);
1104	>	FVECT v;
1105	>	int pos[2];
1106	>	rbf->rbfa[n].peak = peak0 * mv * rbf->vtotal;
1107	>	rbf->rbfa[n].crad = ANG2R(sqrt(rad0rad0(1.-t) +
1108	>	rad1rad1t));
1109	>	ovec_from_pos(v, rbf1j->gx, rbf1j->gy);
1110	>	geodesic(v, v0, v, t, GEOD_REL);
1111	>	pos_from_vec(pos, v);
1112	>	rbf->rbfa[n].gx = pos[0];
1113	>	rbf->rbfa[n].gy = pos[1];
1114	>	++n;
1115	>	}
1116	>	}
1117	>	rbf->vtotal = mig->rbfv[0]->vtotal; / turn ratio into actual */
1118	>	return(rbf);
1119	>	memerr:
1120	>	fputs("Out of memory in e_advect_rbf()\n", stderr);
1121	>	exit(1);
1122	>	return(NULL); /* pro forma return */
1123	>	}
1124	>
1125	>	/* Insert vertex in ordered list */
1126	>	static void
1127	>	insert_vert(RBFNODE *vlist, RBFNODE v)
1128	>	{
1129	>	int i, j;
1130	>
1131	>	for (i = 0; vlist[i] != NULL; i++) {
1132	>	if (v == vlist[i])
1133	>	return;
1134	>	if (v < vlist[i])
1135	>	break;
1136	>	}
1137	>	for (j = i; vlist[j] != NULL; j++)
1138	>	;
1139	>	while (j > i) {
1140	>	vlist[j] = vlist[j-1];
1141	>	--j;
1142	>	}
1143	>	vlist[i] = v;
1144	>	}
1145	>
1146	>	/* Sort triangle edges in standard order */
1147	>	static void
1148	>	order_triangle(MIGRATION *miga[3])
1149	>	{
1150	>	RBFNODE *vert[4];
1151	>	MIGRATION *ord[3];
1152	>	int i;
1153	>	/* order vertices, first */
1154	>	memset(vert, 0, sizeof(vert));
1155	>	for (i = 0; i < 3; i++) {
1156	>	insert_vert(vert, miga[i]->rbfv[0]);
1157	>	insert_vert(vert, miga[i]->rbfv[1]);
1158	>	}
1159	>	/* identify edge 0 */
1160	>	for (i = 0; i < 3; i++)
1161	>	if (miga[i]->rbfv[0] == vert[0] &&
1162	>	miga[i]->rbfv[1] == vert[1]) {
1163	>	ord[0] = miga[i];
1164	>	break;
1165	>	}
1166	>	/* identify edge 1 */
1167	>	for (i = 0; i < 3; i++)
1168	>	if (miga[i]->rbfv[0] == vert[1] &&
1169	>	miga[i]->rbfv[1] == vert[2]) {
1170	>	ord[1] = miga[i];
1171	>	break;
1172	>	}
1173	>	/* identify edge 2 */
1174	>	for (i = 0; i < 3; i++)
1175	>	if (miga[i]->rbfv[0] == vert[0] &&
1176	>	miga[i]->rbfv[1] == vert[2]) {
1177	>	ord[2] = miga[i];
1178	>	break;
1179	>	}
1180	>	miga[0] = ord[0]; miga[1] = ord[1]; miga[2] = ord[2];
1181	>	}
1182	>
1183	>	/* Partially advect between recorded incident angles and allocate new RBF */
1184	>	static RBFNODE *
1185	>	advect_rbf(const FVECT invec)
1186	>	{
1187	>	MIGRATION *miga[3];
1188	>	RBFNODE *rbf;
1189	>	float mbfact, mcfact;
1190	>	int n, i, j, k;
1191	>	FVECT v0, v1, v2;
1192	>	double s, t;
1193	>
1194	>	if (!get_interp(miga, invec)) /* can't interpolate? */
1195	>	return(NULL);
1196	>	if (miga[1] == NULL) /* advect along edge? */
1197	>	return(e_advect_rbf(miga[0], invec));
1198	>	/* put in standard order */
1199	>	order_triangle(miga);
1200	>	#ifdef DEBUG
1201	>	if (miga[0]->rbfv[0] != miga[2]->rbfv[0] \|
1202	>	miga[0]->rbfv[1] != miga[1]->rbfv[0] \|
1203	>	miga[1]->rbfv[1] != miga[2]->rbfv[1]) {
1204	>	fputs("Triangle vertex screw-up!\n", stderr);
1205	>	exit(1);
1206	>	}
1207	>	#endif
1208	>	/* figure out position */
1209	>	fcross(v0, miga[2]->rbfv[0]->invec, miga[2]->rbfv[1]->invec);
1210	>	normalize(v0);
1211	>	fcross(v2, miga[1]->rbfv[0]->invec, miga[1]->rbfv[1]->invec);
1212	>	normalize(v2);
1213	>	fcross(v1, invec, miga[1]->rbfv[1]->invec);
1214	>	normalize(v1);
1215	>	s = acos(DOT(v0,v1)) / acos(DOT(v0,v2));
1216	>	geodesic(v1, miga[0]->rbfv[0]->invec, miga[0]->rbfv[1]->invec,
1217	>	s, GEOD_REL);
1218	>	t = acos(DOT(v1,invec)) / acos(DOT(v1,miga[1]->rbfv[1]->invec));
1219	>	n = 0; /* count migrating particles */
1220	>	for (i = 0; i < mtx_nrows(miga[0]); i++)
1221	>	for (j = 0; j < mtx_ncols(miga[0]); j++)
1222	>	for (k = (miga[0]->mtx[mtx_ndx(miga[0],i,j)] > FTINY) *
1223	>	mtx_ncols(miga[2]); k--; )
1224	>	n += (miga[2]->mtx[mtx_ndx(miga[2],i,k)] > FTINY &&
1225	>	miga[1]->mtx[mtx_ndx(miga[1],j,k)] > FTINY);
1226	>	#ifdef DEBUG
1227	>	fprintf(stderr, "Input RBFs have %d, %d, %d nodes -> output has %d\n",
1228	>	miga[0]->rbfv[0]->nrbf, miga[0]->rbfv[1]->nrbf,
1229	>	miga[2]->rbfv[1]->nrbf, n);
1230	>	#endif
1231	>	rbf = (RBFNODE )malloc(sizeof(RBFNODE) + sizeof(RBFVAL)(n-1));
1232		if (rbf == NULL) {
1233		fputs("Out of memory in advect_rbf()\n", stderr);
1234		exit(1);
1235		}
1236		rbf->next = NULL; rbf->ejl = NULL;
1237		VCOPY(rbf->invec, invec);
1238	<	rbf->vtotal = 0;
1239	<	rbf->nrbf = to_rbf->nrbf;
1240	<
1241	<	return rbf;
1238	>	rbf->nrbf = n;
1239	>	n = 0; /* compute RBF lobes */
1240	>	mbfact = s * miga[0]->rbfv[1]->vtotal/miga[0]->rbfv[0]->vtotal *
1241	>	(1.-t + t*miga[1]->rbfv[1]->vtotal/miga[1]->rbfv[0]->vtotal);
1242	>	mcfact = (1.-s) *
1243	>	(1.-t + t*miga[2]->rbfv[1]->vtotal/miga[2]->rbfv[0]->vtotal);
1244	>	for (i = 0; i < mtx_nrows(miga[0]); i++) {
1245	>	const RBFVAL *rbf0i = &miga[0]->rbfv[0]->rbfa[i];
1246	>	const float w0i = rbf0i->peak;
1247	>	const double rad0i = R2ANG(rbf0i->crad);
1248	>	ovec_from_pos(v0, rbf0i->gx, rbf0i->gy);
1249	>	for (j = 0; j < mtx_ncols(miga[0]); j++) {
1250	>	const float ma = miga[0]->mtx[mtx_ndx(miga[0],i,j)];
1251	>	const RBFVAL *rbf1j;
1252	>	double rad1j, srad2;
1253	>	if (ma <= FTINY)
1254	>	continue;
1255	>	rbf1j = &miga[0]->rbfv[1]->rbfa[j];
1256	>	rad1j = R2ANG(rbf1j->crad);
1257	>	srad2 = (1.-s)(1.-t)rad0irad0i + s(1.-t)rad1jrad1j;
1258	>	ovec_from_pos(v1, rbf1j->gx, rbf1j->gy);
1259	>	geodesic(v1, v0, v1, s, GEOD_REL);
1260	>	for (k = 0; k < mtx_ncols(miga[2]); k++) {
1261	>	float mb = miga[1]->mtx[mtx_ndx(miga[1],j,k)];
1262	>	float mc = miga[2]->mtx[mtx_ndx(miga[2],i,k)];
1263	>	const RBFVAL *rbf2k;
1264	>	double rad2k;
1265	>	FVECT vout;
1266	>	int pos[2];
1267	>	if ((mb <= FTINY) \| (mc <= FTINY))
1268	>	continue;
1269	>	rbf2k = &miga[2]->rbfv[1]->rbfa[k];
1270	>	rbf->rbfa[n].peak = w0i * ma * (mbmbfact + mcmcfact);
1271	>	rad2k = R2ANG(rbf2k->crad);
1272	>	rbf->rbfa[n].crad = ANG2R(sqrt(srad2 + trad2krad2k));
1273	>	ovec_from_pos(v2, rbf2k->gx, rbf2k->gy);
1274	>	geodesic(vout, v1, v2, t, GEOD_REL);
1275	>	pos_from_vec(pos, vout);
1276	>	rbf->rbfa[n].gx = pos[0];
1277	>	rbf->rbfa[n].gy = pos[1];
1278	>	++n;
1279	>	}
1280	>	}
1281	>	}
1282	>	rbf->vtotal = miga[0]->rbfv[0]->vtotal * (mbfact + mcfact);
1283	>	return(rbf);
1284		}
1285	+
1286	+	/* Interpolate and output isotropic BSDF data */
1287	+	static void
1288	+	interp_isotropic()
1289	+	{
1290	+	const int sqres = 1<<samp_order;
1291	+	FILE *ofp = NULL;
1292	+	char cmd[128];
1293	+	int ix, ox, oy;
1294	+	FVECT ivec, ovec;
1295	+	double bsdf;
1296	+	#if DEBUG
1297	+	fprintf(stderr, "Writing isotropic order %d ", samp_order);
1298	+	if (pctcull >= 0) fprintf(stderr, "data with %d%% culling\n", pctcull);
1299	+	else fputs("raw data\n", stderr);
1300		#endif
1301	+	if (pctcull >= 0) { /* begin output */
1302	+	sprintf(cmd, "rttree_reduce -h -a -fd -r 3 -t %d -g %d",
1303	+	pctcull, samp_order);
1304	+	fflush(stdout);
1305	+	ofp = popen(cmd, "w");
1306	+	if (ofp == NULL) {
1307	+	fputs("Cannot create pipe for rttree_reduce\n", stderr);
1308	+	exit(1);
1309	+	}
1310	+	} else
1311	+	fputs("{\n", stdout);
1312	+	/* run through directions */
1313	+	for (ix = 0; ix < sqres/2; ix++) {
1314	+	RBFNODE *rbf;
1315	+	SDsquare2disk(ivec, (ix+.5)/sqres, .5);
1316	+	ivec[2] = input_orient *
1317	+	sqrt(1. - ivec[0]ivec[0] - ivec[1]ivec[1]);
1318	+	rbf = advect_rbf(ivec);
1319	+	for (ox = 0; ox < sqres; ox++)
1320	+	for (oy = 0; oy < sqres; oy++) {
1321	+	SDsquare2disk(ovec, (ox+.5)/sqres, (oy+.5)/sqres);
1322	+	ovec[2] = output_orient *
1323	+	sqrt(1. - ovec[0]ovec[0] - ovec[1]ovec[1]);
1324	+	bsdf = eval_rbfrep(rbf, ovec) / fabs(ovec[2]);
1325	+	if (pctcull >= 0)
1326	+	fwrite(&bsdf, sizeof(bsdf), 1, ofp);
1327	+	else
1328	+	printf("\t%.3e\n", bsdf);
1329	+	}
1330	+	free(rbf);
1331	+	}
1332	+	if (pctcull >= 0) { /* finish output */
1333	+	if (pclose(ofp)) {
1334	+	fprintf(stderr, "Error running '%s'\n", cmd);
1335	+	exit(1);
1336	+	}
1337	+	} else {
1338	+	for (ix = sqressqressqres/2; ix--; )
1339	+	fputs("\t0\n", stdout);
1340	+	fputs("}\n", stdout);
1341	+	}
1342	+	}
1343
1344	+	/* Interpolate and output anisotropic BSDF data */
1345	+	static void
1346	+	interp_anisotropic()
1347	+	{
1348	+	const int sqres = 1<<samp_order;
1349	+	FILE *ofp = NULL;
1350	+	char cmd[128];
1351	+	int ix, iy, ox, oy;
1352	+	FVECT ivec, ovec;
1353	+	double bsdf;
1354	+	#if DEBUG
1355	+	fprintf(stderr, "Writing anisotropic order %d ", samp_order);
1356	+	if (pctcull >= 0) fprintf(stderr, "data with %d%% culling\n", pctcull);
1357	+	else fputs("raw data\n", stderr);
1358	+	#endif
1359	+	if (pctcull >= 0) { /* begin output */
1360	+	sprintf(cmd, "rttree_reduce -h -a -fd -r 4 -t %d -g %d",
1361	+	pctcull, samp_order);
1362	+	fflush(stdout);
1363	+	ofp = popen(cmd, "w");
1364	+	if (ofp == NULL) {
1365	+	fputs("Cannot create pipe for rttree_reduce\n", stderr);
1366	+	exit(1);
1367	+	}
1368	+	} else
1369	+	fputs("{\n", stdout);
1370	+	/* run through directions */
1371	+	for (ix = 0; ix < sqres; ix++)
1372	+	for (iy = 0; iy < sqres; iy++) {
1373	+	RBFNODE *rbf;
1374	+	SDsquare2disk(ivec, (ix+.5)/sqres, (iy+.5)/sqres);
1375	+	ivec[2] = input_orient *
1376	+	sqrt(1. - ivec[0]ivec[0] - ivec[1]ivec[1]);
1377	+	rbf = advect_rbf(ivec);
1378	+	for (ox = 0; ox < sqres; ox++)
1379	+	for (oy = 0; oy < sqres; oy++) {
1380	+	SDsquare2disk(ovec, (ox+.5)/sqres, (oy+.5)/sqres);
1381	+	ovec[2] = output_orient *
1382	+	sqrt(1. - ovec[0]ovec[0] - ovec[1]ovec[1]);
1383	+	bsdf = eval_rbfrep(rbf, ovec) / fabs(ovec[2]);
1384	+	if (pctcull >= 0)
1385	+	fwrite(&bsdf, sizeof(bsdf), 1, ofp);
1386	+	else
1387	+	printf("\t%.3e\n", bsdf);
1388	+	}
1389	+	free(rbf);
1390	+	}
1391	+	if (pctcull >= 0) { /* finish output */
1392	+	if (pclose(ofp)) {
1393	+	fprintf(stderr, "Error running '%s'\n", cmd);
1394	+	exit(1);
1395	+	}
1396	+	} else
1397	+	fputs("}\n", stdout);
1398	+	}
1399	+
1400		#if 1
1401	+	/* Read in BSDF files and interpolate as tensor tree representation */
1402	+	int
1403	+	main(int argc, char *argv[])
1404	+	{
1405	+	RBFNODE *rbf;
1406	+	double bsdf;
1407	+	int i;
1408	+
1409	+	progname = argv[0];
1410	+	if (argc > 2 && !strcmp(argv[1], "-t")) {
1411	+	pctcull = atoi(argv[2]);
1412	+	argv += 2; argc -= 2;
1413	+	}
1414	+	if (argc < 3) {
1415	+	fprintf(stderr,
1416	+	"Usage: %s [-t pctcull] meas1.dat meas2.dat .. > bsdf.xml\n",
1417	+	progname);
1418	+	return(1);
1419	+	}
1420	+	for (i = 1; i < argc; i++) { /* compile measurements */
1421	+	if (!load_pabopto_meas(argv[i]))
1422	+	return(1);
1423	+	compute_radii();
1424	+	cull_values();
1425	+	make_rbfrep();
1426	+	}
1427	+	build_mesh(); /* create interpolation */
1428	+	/* xml_prologue(); /* start XML output */
1429	+	if (single_plane_incident) /* resample dist. */
1430	+	interp_isotropic();
1431	+	else
1432	+	interp_anisotropic();
1433	+	/* xml_epilogue(); /* finish XML output */
1434	+	return(0);
1435	+	}
1436	+	#else
1437		/* Test main produces a Radiance model from the given input file */
1438		int
1439		main(int argc, char *argv[])
#	Line 643 \| Line 1448 \| *main(int argc, char argv[])**
1448		fprintf(stderr, "Usage: %s input.dat > output.rad\n", argv[0]);
1449		return(1);
1450		}
1451	<	if (!load_bsdf_meas(argv[1]))
1451	>	if (!load_pabopto_meas(argv[1]))
1452		return(1);
1453
1454		compute_radii();
#	Line 656 \| Line 1461 \| *main(int argc, char argv[])**
1461		for (i = 0; i < GRIDRES; i++)
1462		for (j = 0; j < GRIDRES; j++)
1463		if (dsf_grid[i][j].vsum > .0f) {
1464	<	vec_from_pos(dir, i, j);
1464	>	ovec_from_pos(dir, i, j);
1465		bsdf = dsf_grid[i][j].vsum / dir[2];
1466		if (dsf_grid[i][j].nval) {
1467		printf("pink cone c%04d\n0\n0\n8\n", ++n);
#	Line 667 \| Line 1472 \| *main(int argc, char argv[])**
1472		dir[2]*(bsdf+.005));
1473		puts("\t.003\t0\n");
1474		} else {
1475	<	vec_from_pos(dir, i, j);
1475	>	ovec_from_pos(dir, i, j);
1476		printf("yellow sphere s%04d\n0\n0\n", ++n);
1477		printf("4 %.6g %.6g %.6g .0015\n\n",
1478		dir[0]bsdf, dir[1]bsdf, dir[2]*bsdf);
#	Line 685 \| Line 1490 \| *main(int argc, char argv[])**
1490		}
1491		for (i = 0; i < GRIDRES; i++)
1492		for (j = 0; j < GRIDRES; j++) {
1493	<	vec_from_pos(dir, i, j);
1493	>	ovec_from_pos(dir, i, j);
1494		bsdf = eval_rbfrep(dsf_list, dir) / dir[2];
1495		fprintf(pfp, "%.8e %.8e %.8e\n",
1496		dir[0]bsdf, dir[1]bsdf, dir[2]*bsdf);

Diff Legend

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

Comparing ray/src/cv/pabopto2xml.c (file contents): Revision 2.8 by greg, Wed Sep 5 00:39:32 2012 UTC vs. Revision 2.13 by greg, Fri Sep 21 05:17:22 2012 UTC

Diff Legend

Comparing ray/src/cv/pabopto2xml.c (file contents):
Revision 2.8 by greg, Wed Sep 5 00:39:32 2012 UTC vs.
Revision 2.13 by greg, Fri Sep 21 05:17:22 2012 UTC