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.2 by greg, Fri Aug 24 20:55:28 2012 UTC vs.
Revision 2.19 by greg, Thu Oct 18 04:28:20 2012 UTC

#	Line 8 | Line 8 | static const char RCSid[] = "$Id$";
8		*      G.Ward
9		*/
10
11	+	#ifndef _WIN32
12	+	#include <unistd.h>
13	+	#include <sys/wait.h>
14	+	#include <sys/mman.h>
15	+	#endif
16		#define _USE_MATH_DEFINES
17		#include <stdio.h>
18		#include <stdlib.h>
#	Line 16 | Line 21 | static const char RCSid[] = "$Id$";
21		#include <math.h>
22		#include "bsdf.h"
23
24	+	#define DEBUG           1
25	+
26		#ifndef GRIDRES
27	<	#define GRIDRES         200             /* max. grid resolution per side */
27	>	#define GRIDRES         200             /* grid resolution per side */
28		#endif
29
30	<	#define RSCA            3.              /* radius scaling factor (empirical) */
24	<	#define MSCA            .2              /* magnitude scaling (empirical) */
30	>	#define MAXSAMPORD      7               /* don't sample finer than this */
31
32	<	#define R2ANG(c)        (((c)+.5)*(M_PI/(1<<16)))
32	>	#define RSCA            2.7             /* radius scaling factor (empirical) */
33	>
34	>	/* convert to/from coded radians */
35		#define ANG2R(r)        (int)((r)*((1<<16)/M_PI))
36	+	#define R2ANG(c)        (((c)+.5)*(M_PI/(1<<16)))
37
38		typedef struct {
39	<	float           vsum;           /* BSDF sum */
39	>	float           vsum;           /* DSF sum */
40		unsigned short  nval;           /* number of values in sum */
41		unsigned short  crad;           /* radius (coded angle) */
42		} GRIDVAL;                      /* grid value */
43
44		typedef struct {
45	<	float           bsdf;           /* lobe value at peak */
45	>	float           peak;           /* lobe value at peak */
46		unsigned short  crad;           /* radius (coded angle) */
47		unsigned char   gx, gy;         /* grid position */
48		} RBFVAL;                       /* radial basis function value */
49
50	<	typedef struct s_rbflist {
51	<	struct s_rbflist        *next;          /* next in our RBF list */
50	>	struct s_rbfnode;               /* forward declaration of RBF struct */
51	>
52	>	typedef struct s_migration {
53	>	struct s_migration      *next;          /* next in global edge list */
54	>	struct s_rbfnode        *rbfv[2];       /* from,to vertex */
55	>	struct s_migration      *enxt[2];       /* next from,to sibling */
56	>	float                   mtx[1];         /* matrix (extends struct) */
57	>	} MIGRATION;                    /* migration link (winged edge structure) */
58	>
59	>	typedef struct s_rbfnode {
60	>	struct s_rbfnode        *next;          /* next in global RBF list */
61	>	MIGRATION               *ejl;           /* edge list for this vertex */
62		FVECT                   invec;          /* incident vector direction */
63	+	double                  vtotal;         /* volume for normalization */
64		int                     nrbf;           /* number of RBFs */
65		RBFVAL                  rbfa[1];        /* RBF array (extends struct) */
66	<	} RBFLIST;                      /* RBF representation of BSDF @ 1 incidence */
66	>	} RBFNODE;                      /* RBF representation of DSF @ 1 incidence */
67
68		/* our loaded grid for this incident angle */
69	<	static double   theta_in_deg, phi_in_deg;
70	<	static GRIDVAL  bsdf_grid[GRIDRES][GRIDRES];
69	>	static double           theta_in_deg, phi_in_deg;
70	>	static GRIDVAL          dsf_grid[GRIDRES][GRIDRES];
71
72	<	/* processed incident BSDF measurements */
73	<	static RBFLIST  *bsdf_list = NULL;
72	>	/* all incident angles in-plane so far? */
73	>	static int              single_plane_incident = -1;
74
75	<	/* Count up non-empty nodes and build RBF representation from current grid */
76	<	static RBFLIST *
77	<	make_rbfrep(void)
75	>	/* represented incident quadrants */
76	>	#define INP_QUAD1       1               /* 0-90 degree quadrant */
77	>	#define INP_QUAD2       2               /* 90-180 degree quadrant */
78	>	#define INP_QUAD3       4               /* 180-270 degree quadrant */
79	>	#define INP_QUAD4       8               /* 270-360 degree quadrant */
80	>
81	>	static int              inp_coverage = 0;
82	>
83	>	/* symmetry operations */
84	>	#define MIRROR_X        1               /* mirror(ed) x-coordinate */
85	>	#define MIRROR_Y        2               /* mirror(ed) y-coordinate */
86	>
87	>	/* input/output orientations */
88	>	static int              input_orient = 0;
89	>	static int              output_orient = 0;
90	>
91	>	/* processed incident DSF measurements */
92	>	static RBFNODE          *dsf_list = NULL;
93	>
94	>	/* RBF-linking matrices (edges) */
95	>	static MIGRATION        *mig_list = NULL;
96	>
97	>	/* migration edges drawn in raster fashion */
98	>	static MIGRATION        *mig_grid[GRIDRES][GRIDRES];
99	>
100	>	#define mtx_nrows(m)    ((m)->rbfv[0]->nrbf)
101	>	#define mtx_ncols(m)    ((m)->rbfv[1]->nrbf)
102	>	#define mtx_ndx(m,i,j)  ((i)*mtx_ncols(m) + (j))
103	>	#define is_src(rbf,m)   ((rbf) == (m)->rbfv[0])
104	>	#define is_dest(rbf,m)  ((rbf) == (m)->rbfv[1])
105	>	#define nextedge(rbf,m) (m)->enxt[is_dest(rbf,m)]
106	>	#define opp_rbf(rbf,m)  (m)->rbfv[is_src(rbf,m)]
107	>
108	>	#define round(v)        (int)((v) + .5 - ((v) < -.5))
109	>
110	>	char                    *progname;
111	>
112	>	/* percentage to cull (<0 to turn off) */
113	>	int                     pctcull = 90;
114	>	/* number of processes to run */
115	>	int                     nprocs = 1;
116	>
117	>	/* number of children (-1 in child) */
118	>	int                     nchild = 0;
119	>
120	>	/* sampling order (set by data density) */
121	>	int                     samp_order = 0;
122	>
123	>	/* get phi value in degrees, [0,360) range */
124	>	#define get_phi360(v)   ((180./M_PI)*atan2((v)[1],(v)[0]) + 180.)
125	>
126	>	/* Apply symmetry to the given vector based on distribution */
127	>	static int
128	>	use_symmetry(FVECT vec)
129		{
130	<	int     nn = 0;
131	<	RBFLIST *newnode;
132	<	int     i, j;
133	<	/* count non-empty bins */
134	<	for (i = 0; i < GRIDRES; i++)
135	<	for (j = 0; j < GRIDRES; j++)
136	<	nn += (bsdf_grid[i][j].nval > 0);
137	<	/* allocate RBF array */
138	<	newnode = (RBFLIST *)malloc(sizeof(RBFLIST) + sizeof(RBFVAL)*(nn-1));
139	<	if (newnode == NULL) {
140	<	fputs("Out of memory in make_rbfrep\n", stderr);
130	>	double  phi = get_phi360(vec);
131	>
132	>	switch (inp_coverage) {
133	>	case INP_QUAD1|INP_QUAD2|INP_QUAD3|INP_QUAD4:
134	>	break;
135	>	case INP_QUAD1|INP_QUAD2:
136	>	if ((-FTINY > phi) | (phi > 180.+FTINY))
137	>	goto mir_y;
138	>	break;
139	>	case INP_QUAD2|INP_QUAD3:
140	>	if ((90.-FTINY > phi) | (phi > 270.+FTINY))
141	>	goto mir_x;
142	>	break;
143	>	case INP_QUAD3|INP_QUAD4:
144	>	if ((180.-FTINY > phi) | (phi > 360.+FTINY))
145	>	goto mir_y;
146	>	break;
147	>	case INP_QUAD4|INP_QUAD1:
148	>	if ((270.-FTINY > phi) & (phi > 90.+FTINY))
149	>	goto mir_x;
150	>	break;
151	>	case INP_QUAD1:
152	>	if ((-FTINY > phi) | (phi > 90.+FTINY))
153	>	switch ((int)(phi*(1./90.))) {
154	>	case 1: goto mir_x;
155	>	case 2: goto mir_xy;
156	>	case 3: goto mir_y;
157	>	}
158	>	break;
159	>	case INP_QUAD2:
160	>	if ((90.-FTINY > phi) | (phi > 180.+FTINY))
161	>	switch ((int)(phi*(1./90.))) {
162	>	case 0: goto mir_x;
163	>	case 2: goto mir_y;
164	>	case 3: goto mir_xy;
165	>	}
166	>	break;
167	>	case INP_QUAD3:
168	>	if ((180.-FTINY > phi) | (phi > 270.+FTINY))
169	>	switch ((int)(phi*(1./90.))) {
170	>	case 0: goto mir_xy;
171	>	case 1: goto mir_y;
172	>	case 3: goto mir_x;
173	>	}
174	>	break;
175	>	case INP_QUAD4:
176	>	if ((270.-FTINY > phi) | (phi > 360.+FTINY))
177	>	switch ((int)(phi*(1./90.))) {
178	>	case 0: goto mir_y;
179	>	case 1: goto mir_xy;
180	>	case 2: goto mir_x;
181	>	}
182	>	break;
183	>	default:
184	>	fprintf(stderr, "%s: Illegal input coverage (%d)\n",
185	>	progname, inp_coverage);
186		exit(1);
187		}
188	<	newnode->invec[2] = sin(M_PI/180.*theta_in_deg);
189	<	newnode->invec[0] = cos(M_PI/180.*phi_in_deg)*newnode->invec[2];
190	<	newnode->invec[1] = sin(M_PI/180.*phi_in_deg)*newnode->invec[2];
191	<	newnode->invec[2] = sqrt(1. - newnode->invec[2]*newnode->invec[2]);
192	<	newnode->nrbf = nn;
193	<	nn = 0;                 /* fill RBF array */
194	<	for (i = 0; i < GRIDRES; i++)
195	<	for (j = 0; j < GRIDRES; j++)
196	<	if (bsdf_grid[i][j].nval) {
197	<	newnode->rbfa[nn].bsdf = MSCA*bsdf_grid[i][j].vsum /
198	<	(double)bsdf_grid[i][j].nval;
83	<	newnode->rbfa[nn].crad = RSCA*bsdf_grid[i][j].crad + .5;
84	<	newnode->rbfa[nn].gx = i;
85	<	newnode->rbfa[nn].gy = j;
86	<	++nn;
87	<	}
88	<	newnode->next = bsdf_list;
89	<	return(bsdf_list = newnode);
188	>	return(0);              /* in range */
189	>	mir_x:
190	>	vec[0] = -vec[0];
191	>	return(MIRROR_X);
192	>	mir_y:
193	>	vec[1] = -vec[1];
194	>	return(MIRROR_Y);
195	>	mir_xy:
196	>	vec[0] = -vec[0];
197	>	vec[1] = -vec[1];
198	>	return(MIRROR_X|MIRROR_Y);
199		}
200
201	<	/* Compute grid position from normalized outgoing vector */
201	>	/* Reverse symmetry based on what was done before */
202		static void
203	<	pos_from_vec(int pos[2], const FVECT vec)
203	>	rev_symmetry(FVECT vec, int sym)
204		{
205	<	double  sq[2];          /* uniform hemispherical projection */
206	<	double  norm = 1./sqrt(1. + vec[2]);
205	>	if (sym & MIRROR_X)
206	>	vec[0] = -vec[0];
207	>	if (sym & MIRROR_Y)
208	>	vec[1] = -vec[1];
209	>	}
210
211	<	SDdisk2square(sq, vec[0]*norm, vec[1]*norm);
211	>	/* Reverse symmetry for an RBF distribution */
212	>	static void
213	>	rev_rbf_symmetry(RBFNODE *rbf, int sym)
214	>	{
215	>	int     n;
216
217	<	pos[0] = (int)(sq[0]*GRIDRES);
218	<	pos[1] = (int)(sq[1]*GRIDRES);
217	>	rev_symmetry(rbf->invec, sym);
218	>	if (sym & MIRROR_X)
219	>	for (n = rbf->nrbf; n-- > 0; )
220	>	rbf->rbfa[n].gx = GRIDRES-1 - rbf->rbfa[n].gx;
221	>	if (sym & MIRROR_Y)
222	>	for (n = rbf->nrbf; n-- > 0; )
223	>	rbf->rbfa[n].gy = GRIDRES-1 - rbf->rbfa[n].gy;
224		}
225
226	+	/* Compute volume associated with Gaussian lobe */
227	+	static double
228	+	rbf_volume(const RBFVAL *rbfp)
229	+	{
230	+	double  rad = R2ANG(rbfp->crad);
231	+
232	+	return((2.*M_PI) * rbfp->peak * rad*rad);
233	+	}
234	+
235		/* Compute outgoing vector from grid position */
236		static void
237	<	vec_from_pos(FVECT vec, int xpos, int ypos)
237	>	ovec_from_pos(FVECT vec, int xpos, int ypos)
238		{
239		double  uv[2];
240		double  r2;
#	Line 115 | Line 245 | vec_from_pos(FVECT vec, int xpos, int ypos)
245		vec[0] = vec[1] = sqrt(2. - r2);
246		vec[0] *= uv[0];
247		vec[1] *= uv[1];
248	<	vec[2] = 1. - r2;
248	>	vec[2] = output_orient*(1. - r2);
249		}
250
251	<	/* Evaluate RBF for BSDF at the given normalized outgoing direction */
251	>	/* Compute grid position from normalized input/output vector */
252	>	static void
253	>	pos_from_vec(int pos[2], const FVECT vec)
254	>	{
255	>	double  sq[2];          /* uniform hemispherical projection */
256	>	double  norm = 1./sqrt(1. + fabs(vec[2]));
257	>
258	>	SDdisk2square(sq, vec[0]*norm, vec[1]*norm);
259	>
260	>	pos[0] = (int)(sq[0]*GRIDRES);
261	>	pos[1] = (int)(sq[1]*GRIDRES);
262	>	}
263	>
264	>	/* Evaluate RBF for DSF at the given normalized outgoing direction */
265		static double
266	<	eval_rbfrep(const RBFLIST *rp, const FVECT outvec)
266	>	eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
267		{
268		double          res = .0;
269		const RBFVAL    *rbfp;
#	Line 128 | Line 271 | eval_rbfrep(const RBFLIST *rp, const FVECT outvec)
271		double          sig2;
272		int             n;
273
274	+	if (rp == NULL)
275	+	return(.0);
276		rbfp = rp->rbfa;
277		for (n = rp->nrbf; n--; rbfp++) {
278	<	vec_from_pos(odir, rbfp->gx, rbfp->gy);
278	>	ovec_from_pos(odir, rbfp->gx, rbfp->gy);
279		sig2 = R2ANG(rbfp->crad);
280		sig2 = (DOT(odir,outvec) - 1.) / (sig2*sig2);
281		if (sig2 > -19.)
282	<	res += rbfp->bsdf * exp(sig2);
282	>	res += rbfp->peak * exp(sig2);
283		}
284		return(res);
285		}
286
287	+	/* Insert a new directional scattering function in our global list */
288	+	static void
289	+	insert_dsf(RBFNODE *newrbf)
290	+	{
291	+	RBFNODE         *rbf, *rbf_last;
292	+	/* check for redundant meas. */
293	+	for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
294	+	if (DOT(rbf->invec, newrbf->invec) >= 1.-FTINY) {
295	+	fprintf(stderr,
296	+	"%s: Duplicate incident measurement (ignored)\n",
297	+	progname);
298	+	free(newrbf);
299	+	return;
300	+	}
301	+	/* keep in ascending theta order */
302	+	for (rbf_last = NULL, rbf = dsf_list;
303	+	single_plane_incident & (rbf != NULL);
304	+	rbf_last = rbf, rbf = rbf->next)
305	+	if (input_orient*rbf->invec[2] < input_orient*newrbf->invec[2])
306	+	break;
307	+	if (rbf_last == NULL) {
308	+	newrbf->next = dsf_list;
309	+	dsf_list = newrbf;
310	+	return;
311	+	}
312	+	newrbf->next = rbf;
313	+	rbf_last->next = newrbf;
314	+	}
315	+
316	+	/* Count up filled nodes and build RBF representation from current grid */
317	+	static RBFNODE *
318	+	make_rbfrep(void)
319	+	{
320	+	int     niter = 16;
321	+	int     minrad = ANG2R(pow(2., 1.-samp_order));
322	+	double  lastVar, thisVar = 100.;
323	+	int     nn;
324	+	RBFNODE *newnode;
325	+	int     i, j;
326	+
327	+	nn = 0;                 /* count selected bins */
328	+	for (i = 0; i < GRIDRES; i++)
329	+	for (j = 0; j < GRIDRES; j++)
330	+	nn += dsf_grid[i][j].nval;
331	+	/* allocate RBF array */
332	+	newnode = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(nn-1));
333	+	if (newnode == NULL) {
334	+	fprintf(stderr, "%s: Out of memory in make_rbfrep()\n", progname);
335	+	exit(1);
336	+	}
337	+	newnode->next = NULL;
338	+	newnode->ejl = NULL;
339	+	newnode->invec[2] = sin(M_PI/180.*theta_in_deg);
340	+	newnode->invec[0] = cos(M_PI/180.*phi_in_deg)*newnode->invec[2];
341	+	newnode->invec[1] = sin(M_PI/180.*phi_in_deg)*newnode->invec[2];
342	+	newnode->invec[2] = input_orient*sqrt(1. - newnode->invec[2]*newnode->invec[2]);
343	+	newnode->vtotal = 0;
344	+	newnode->nrbf = nn;
345	+	nn = 0;                 /* fill RBF array */
346	+	for (i = 0; i < GRIDRES; i++)
347	+	for (j = 0; j < GRIDRES; j++)
348	+	if (dsf_grid[i][j].nval) {
349	+	newnode->rbfa[nn].peak = dsf_grid[i][j].vsum;
350	+	newnode->rbfa[nn].crad = RSCA*dsf_grid[i][j].crad + .5;
351	+	newnode->rbfa[nn].gx = i;
352	+	newnode->rbfa[nn].gy = j;
353	+	if (newnode->rbfa[nn].crad < minrad)
354	+	minrad = newnode->rbfa[nn].crad;
355	+	++nn;
356	+	}
357	+	/* iterate to improve interpolation accuracy */
358	+	do {
359	+	double  dsum = 0, dsum2 = 0;
360	+	nn = 0;
361	+	for (i = 0; i < GRIDRES; i++)
362	+	for (j = 0; j < GRIDRES; j++)
363	+	if (dsf_grid[i][j].nval) {
364	+	FVECT   odir;
365	+	double  corr;
366	+	ovec_from_pos(odir, i, j);
367	+	newnode->rbfa[nn++].peak *= corr =
368	+	dsf_grid[i][j].vsum /
369	+	eval_rbfrep(newnode, odir);
370	+	dsum += corr - 1.;
371	+	dsum2 += (corr-1.)*(corr-1.);
372	+	}
373	+	lastVar = thisVar;
374	+	thisVar = dsum2/(double)nn;
375	+	#ifdef DEBUG
376	+	fprintf(stderr, "Avg., RMS error: %.1f%%  %.1f%%\n",
377	+	100.*dsum/(double)nn,
378	+	100.*sqrt(thisVar));
379	+	#endif
380	+	} while (--niter > 0 && lastVar-thisVar > 0.02*lastVar);
381	+
382	+	nn = 0;                 /* compute sum for normalization */
383	+	while (nn < newnode->nrbf)
384	+	newnode->vtotal += rbf_volume(&newnode->rbfa[nn++]);
385	+
386	+	insert_dsf(newnode);
387	+	/* adjust sampling resolution */
388	+	samp_order = log(2./R2ANG(minrad))/M_LN2 + .5;
389	+	if (samp_order > MAXSAMPORD)
390	+	samp_order = MAXSAMPORD;
391	+
392	+	return(newnode);
393	+	}
394	+
395		/* Load a set of measurements corresponding to a particular incident angle */
396		static int
397	<	load_bsdf_meas(const char *fname)
397	>	load_pabopto_meas(const char *fname)
398		{
399		FILE    *fp = fopen(fname, "r");
400		int     inp_is_DSF = -1;
401	<	double  theta_out, phi_out, val;
401	>	double  new_phi, theta_out, phi_out, val;
402		char    buf[2048];
403		int     n, c;
404
#	Line 154 | Line 407 | load_bsdf_meas(const char *fname)
407		fputs(": cannot open\n", stderr);
408		return(0);
409		}
410	<	memset(bsdf_grid, 0, sizeof(bsdf_grid));
410	>	memset(dsf_grid, 0, sizeof(dsf_grid));
411	>	#ifdef DEBUG
412	>	fprintf(stderr, "Loading measurement file '%s'...\n", fname);
413	>	#endif
414		/* read header information */
415		while ((c = getc(fp)) == '#' || c == EOF) {
416		if (fgets(buf, sizeof(buf), fp) == NULL) {
#	Line 173 | Line 429 | load_bsdf_meas(const char *fname)
429		}
430		if (sscanf(buf, "intheta %lf", &theta_in_deg) == 1)
431		continue;
432	<	if (sscanf(buf, "inphi %lf", &phi_in_deg) == 1)
432	>	if (sscanf(buf, "inphi %lf", &new_phi) == 1)
433		continue;
434		if (sscanf(buf, "incident_angle %lf %lf",
435	<	&theta_in_deg, &phi_in_deg) == 2)
435	>	&theta_in_deg, &new_phi) == 2)
436		continue;
437		}
438		if (inp_is_DSF < 0) {
#	Line 185 | Line 441 | load_bsdf_meas(const char *fname)
441		fclose(fp);
442		return(0);
443		}
444	<	ungetc(c, fp);          /* read actual data */
444	>	if (!input_orient)              /* check input orientation */
445	>	input_orient = 1 - 2*(theta_in_deg > 90.);
446	>	else if (input_orient > 0 ^ theta_in_deg < 90.) {
447	>	fputs("Cannot handle input angles on both sides of surface\n",
448	>	stderr);
449	>	exit(1);
450	>	}
451	>	if (single_plane_incident > 0)  /* check input coverage */
452	>	single_plane_incident = (round(new_phi) == round(phi_in_deg));
453	>	else if (single_plane_incident < 0)
454	>	single_plane_incident = 1;
455	>	phi_in_deg = new_phi;
456	>	new_phi += 360.*(new_phi < -FTINY);
457	>	if ((1. < new_phi) & (new_phi < 89.))
458	>	inp_coverage |= INP_QUAD1;
459	>	else if ((91. < new_phi) & (new_phi < 179.))
460	>	inp_coverage |= INP_QUAD2;
461	>	else if ((181. < new_phi) & (new_phi < 269.))
462	>	inp_coverage |= INP_QUAD3;
463	>	else if ((271. < new_phi) & (new_phi < 359.))
464	>	inp_coverage |= INP_QUAD4;
465	>	ungetc(c, fp);                  /* read actual data */
466		while (fscanf(fp, "%lf %lf %lf\n", &theta_out, &phi_out, &val) == 3) {
467		FVECT   ovec;
468		int     pos[2];
469
470	+	if (!output_orient)     /* check output orientation */
471	+	output_orient = 1 - 2*(theta_out > 90.);
472	+	else if (output_orient > 0 ^ theta_out < 90.) {
473	+	fputs("Cannot handle output angles on both sides of surface\n",
474	+	stderr);
475	+	exit(1);
476	+	}
477		ovec[2] = sin(M_PI/180.*theta_out);
478		ovec[0] = cos(M_PI/180.*phi_out) * ovec[2];
479		ovec[1] = sin(M_PI/180.*phi_out) * ovec[2];
480		ovec[2] = sqrt(1. - ovec[2]*ovec[2]);
481
482	<	if (inp_is_DSF)
483	<	val /= ovec[2]; /* convert from DSF to BSDF */
482	>	if (!inp_is_DSF)
483	>	val *= ovec[2]; /* convert from BSDF to DSF */
484
485		pos_from_vec(pos, ovec);
486
487	<	bsdf_grid[pos[0]][pos[1]].vsum += val;
488	<	bsdf_grid[pos[0]][pos[1]].nval++;
487	>	dsf_grid[pos[0]][pos[1]].vsum += val;
488	>	dsf_grid[pos[0]][pos[1]].nval++;
489		}
490		n = 0;
491		while ((c = getc(fp)) != EOF)
#	Line 219 | Line 503 | load_bsdf_meas(const char *fname)
503		static void
504		compute_radii(void)
505		{
506	<	unsigned short  fill_grid[GRIDRES][GRIDRES];
506	>	unsigned int    fill_grid[GRIDRES][GRIDRES];
507	>	unsigned short  fill_cnt[GRIDRES][GRIDRES];
508		FVECT           ovec0, ovec1;
509		double          ang2, lastang2;
225	–	int             r2, lastr2;
510		int             r, i, j, jn, ii, jj, inear, jnear;
511
512		r = GRIDRES/2;                          /* proceed in zig-zag */
513		for (i = 0; i < GRIDRES; i++)
514		for (jn = 0; jn < GRIDRES; jn++) {
515		j = (i&1) ? jn : GRIDRES-1-jn;
516	<	if (bsdf_grid[i][j].nval)       /* find empty grid pos. */
516	>	if (dsf_grid[i][j].nval)        /* find empty grid pos. */
517		continue;
518	<	vec_from_pos(ovec0, i, j);
518	>	ovec_from_pos(ovec0, i, j);
519		inear = jnear = -1;             /* find nearest non-empty */
520		lastang2 = M_PI*M_PI;
521		for (ii = i-r; ii <= i+r; ii++) {
#	Line 240 | Line 524 | compute_radii(void)
524		for (jj = j-r; jj <= j+r; jj++) {
525		if (jj < 0) continue;
526		if (jj >= GRIDRES) break;
527	<	if (!bsdf_grid[ii][jj].nval)
527	>	if (!dsf_grid[ii][jj].nval)
528		continue;
529	<	vec_from_pos(ovec1, ii, jj);
529	>	ovec_from_pos(ovec1, ii, jj);
530		ang2 = 2. - 2.*DOT(ovec0,ovec1);
531		if (ang2 >= lastang2)
532		continue;
#	Line 251 | Line 535 | compute_radii(void)
535		}
536		}
537		if (inear < 0) {
538	<	fputs("Could not find non-empty neighbor!\n", stderr);
538	>	fprintf(stderr,
539	>	"%s: Could not find non-empty neighbor!\n",
540	>	progname);
541		exit(1);
542		}
543		ang2 = sqrt(lastang2);
544		r = ANG2R(ang2);                /* record if > previous */
545	<	if (r > bsdf_grid[inear][jnear].crad)
546	<	bsdf_grid[inear][jnear].crad = r;
545	>	if (r > dsf_grid[inear][jnear].crad)
546	>	dsf_grid[inear][jnear].crad = r;
547		/* next search radius */
548	<	r = ang2*(2.*GRIDRES/M_PI) + 1;
548	>	r = ang2*(2.*GRIDRES/M_PI) + 3;
549		}
550	<	/* fill in neighbors */
550	>	/* blur radii over hemisphere */
551		memset(fill_grid, 0, sizeof(fill_grid));
552	+	memset(fill_cnt, 0, sizeof(fill_cnt));
553		for (i = 0; i < GRIDRES; i++)
554		for (j = 0; j < GRIDRES; j++) {
555	<	if (!bsdf_grid[i][j].nval)
556	<	continue;               /* no value -- skip */
557	<	if (bsdf_grid[i][j].crad)
271	<	continue;               /* has distance already */
272	<	r = GRIDRES/20;
273	<	lastr2 = 2*r*r + 1;
555	>	if (!dsf_grid[i][j].crad)
556	>	continue;               /* missing distance */
557	>	r = R2ANG(dsf_grid[i][j].crad)*(2.*RSCA*GRIDRES/M_PI);
558		for (ii = i-r; ii <= i+r; ii++) {
559		if (ii < 0) continue;
560		if (ii >= GRIDRES) break;
561		for (jj = j-r; jj <= j+r; jj++) {
562		if (jj < 0) continue;
563		if (jj >= GRIDRES) break;
564	<	if (!bsdf_grid[ii][jj].crad)
564	>	if ((ii-i)*(ii-i) + (jj-j)*(jj-j) > r*r)
565		continue;
566	<	/* OK to use approx. closest */
567	<	r2 = (ii-i)*(ii-i) + (jj-j)*(jj-j);
284	<	if (r2 >= lastr2)
285	<	continue;
286	<	fill_grid[i][j] = bsdf_grid[ii][jj].crad;
287	<	lastr2 = r2;
566	>	fill_grid[ii][jj] += dsf_grid[i][j].crad;
567	>	fill_cnt[ii][jj]++;
568		}
569		}
570		}
571	<	/* copy back filled entries */
571	>	/* copy back blurred radii */
572		for (i = 0; i < GRIDRES; i++)
573		for (j = 0; j < GRIDRES; j++)
574	<	if (fill_grid[i][j])
575	<	bsdf_grid[i][j].crad = fill_grid[i][j];
574	>	if (fill_cnt[i][j])
575	>	dsf_grid[i][j].crad = fill_grid[i][j]/fill_cnt[i][j];
576		}
577
578	<	/* Cull points for more uniform distribution */
578	>	/* Cull points for more uniform distribution, leave all nval 0 or 1 */
579		static void
580		cull_values(void)
581		{
#	Line 305 | Line 585 | cull_values(void)
585		/* simple greedy algorithm */
586		for (i = 0; i < GRIDRES; i++)
587		for (j = 0; j < GRIDRES; j++) {
588	<	if (!bsdf_grid[i][j].nval)
588	>	if (!dsf_grid[i][j].nval)
589		continue;
590	<	if (!bsdf_grid[i][j].crad)
590	>	if (!dsf_grid[i][j].crad)
591		continue;               /* shouldn't happen */
592	<	vec_from_pos(ovec0, i, j);
593	<	maxang = 2.*R2ANG(bsdf_grid[i][j].crad);
592	>	ovec_from_pos(ovec0, i, j);
593	>	maxang = 2.*R2ANG(dsf_grid[i][j].crad);
594		if (maxang > ovec0[2])          /* clamp near horizon */
595		maxang = ovec0[2];
596		r = maxang*(2.*GRIDRES/M_PI) + 1;
#	Line 321 | Line 601 | cull_values(void)
601		for (jj = j-r; jj <= j+r; jj++) {
602		if (jj < 0) continue;
603		if (jj >= GRIDRES) break;
604	<	if (!bsdf_grid[ii][jj].nval)
604	>	if (!dsf_grid[ii][jj].nval)
605		continue;
606		if ((ii == i) & (jj == j))
607		continue;       /* don't get self-absorbed */
608	<	vec_from_pos(ovec1, ii, jj);
608	>	ovec_from_pos(ovec1, ii, jj);
609		if (2. - 2.*DOT(ovec0,ovec1) >= maxang2)
610		continue;
611		/* absorb sum */
612	<	bsdf_grid[i][j].vsum += bsdf_grid[ii][jj].vsum;
613	<	bsdf_grid[i][j].nval += bsdf_grid[ii][jj].nval;
612	>	dsf_grid[i][j].vsum += dsf_grid[ii][jj].vsum;
613	>	dsf_grid[i][j].nval += dsf_grid[ii][jj].nval;
614		/* keep value, though */
615	<	bsdf_grid[ii][jj].vsum /= (double)bsdf_grid[ii][jj].nval;
616	<	bsdf_grid[ii][jj].nval = 0;
615	>	dsf_grid[ii][jj].vsum /= (float)dsf_grid[ii][jj].nval;
616	>	dsf_grid[ii][jj].nval = 0;
617		}
618		}
619		}
620	+	/* final averaging pass */
621	+	for (i = 0; i < GRIDRES; i++)
622	+	for (j = 0; j < GRIDRES; j++)
623	+	if (dsf_grid[i][j].nval > 1) {
624	+	dsf_grid[i][j].vsum /= (float)dsf_grid[i][j].nval;
625	+	dsf_grid[i][j].nval = 1;
626	+	}
627		}
628
629	+	/* Compute (and allocate) migration price matrix for optimization */
630	+	static float *
631	+	price_routes(const RBFNODE *from_rbf, const RBFNODE *to_rbf)
632	+	{
633	+	float   *pmtx = (float *)malloc(sizeof(float) *
634	+	from_rbf->nrbf * to_rbf->nrbf);
635	+	FVECT   *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
636	+	int     i, j;
637
638	+	if ((pmtx == NULL) | (vto == NULL)) {
639	+	fprintf(stderr, "%s: Out of memory in migration_costs()\n",
640	+	progname);
641	+	exit(1);
642	+	}
643	+	for (j = to_rbf->nrbf; j--; )           /* save repetitive ops. */
644	+	ovec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);
645	+
646	+	for (i = from_rbf->nrbf; i--; ) {
647	+	const double        from_ang = R2ANG(from_rbf->rbfa[i].crad);
648	+	FVECT               vfrom;
649	+	ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
650	+	for (j = to_rbf->nrbf; j--; )
651	+	pmtx[i*to_rbf->nrbf + j] = acos(DOT(vfrom, vto[j])) +
652	+	fabs(R2ANG(to_rbf->rbfa[j].crad) - from_ang);
653	+	}
654	+	free(vto);
655	+	return(pmtx);
656	+	}
657	+
658	+	/* Comparison routine needed for sorting price row */
659	+	static const float      *price_arr;
660	+	static int
661	+	msrt_cmp(const void *p1, const void *p2)
662	+	{
663	+	float   c1 = price_arr[*(const int *)p1];
664	+	float   c2 = price_arr[*(const int *)p2];
665	+
666	+	if (c1 > c2) return(1);
667	+	if (c1 < c2) return(-1);
668	+	return(0);
669	+	}
670	+
671	+	/* Compute minimum (optimistic) cost for moving the given source material */
672	+	static double
673	+	min_cost(double amt2move, const double *avail, const float *price, int n)
674	+	{
675	+	static int      *price_sort = NULL;
676	+	static int      n_alloc = 0;
677	+	double          total_cost = 0;
678	+	int             i;
679	+
680	+	if (amt2move <= FTINY)                  /* pre-emptive check */
681	+	return(0.);
682	+	if (n > n_alloc) {                      /* (re)allocate sort array */
683	+	if (n_alloc) free(price_sort);
684	+	price_sort = (int *)malloc(sizeof(int)*n);
685	+	if (price_sort == NULL) {
686	+	fprintf(stderr, "%s: Out of memory in min_cost()\n",
687	+	progname);
688	+	exit(1);
689	+	}
690	+	n_alloc = n;
691	+	}
692	+	for (i = n; i--; )
693	+	price_sort[i] = i;
694	+	price_arr = price;
695	+	qsort(price_sort, n, sizeof(int), &msrt_cmp);
696	+	/* move cheapest first */
697	+	for (i = 0; i < n && amt2move > FTINY; i++) {
698	+	int     d = price_sort[i];
699	+	double  amt = (amt2move < avail[d]) ? amt2move : avail[d];
700	+
701	+	total_cost += amt * price[d];
702	+	amt2move -= amt;
703	+	}
704	+	return(total_cost);
705	+	}
706	+
707	+	/* Take a step in migration by choosing optimal bucket to transfer */
708	+	static double
709	+	migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, const float *pmtx)
710	+	{
711	+	const double    maxamt = .1;
712	+	const double    minamt = maxamt*.0001;
713	+	static double   *src_cost = NULL;
714	+	static int      n_alloc = 0;
715	+	struct {
716	+	int     s, d;   /* source and destination */
717	+	double  price;  /* price estimate per amount moved */
718	+	double  amt;    /* amount we can move */
719	+	} cur, best;
720	+	int             i;
721	+
722	+	if (mtx_nrows(mig) > n_alloc) {         /* allocate cost array */
723	+	if (n_alloc)
724	+	free(src_cost);
725	+	src_cost = (double *)malloc(sizeof(double)*mtx_nrows(mig));
726	+	if (src_cost == NULL) {
727	+	fprintf(stderr, "%s: Out of memory in migration_step()\n",
728	+	progname);
729	+	exit(1);
730	+	}
731	+	n_alloc = mtx_nrows(mig);
732	+	}
733	+	for (i = mtx_nrows(mig); i--; )         /* starting costs for diff. */
734	+	src_cost[i] = min_cost(src_rem[i], dst_rem,
735	+	pmtx+i*mtx_ncols(mig), mtx_ncols(mig));
736	+
737	+	/* find best source & dest. */
738	+	best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
739	+	for (cur.s = mtx_nrows(mig); cur.s--; ) {
740	+	const float *price = pmtx + cur.s*mtx_ncols(mig);
741	+	double      cost_others = 0;
742	+	if (src_rem[cur.s] < minamt)
743	+	continue;
744	+	cur.d = -1;                         /* examine cheapest dest. */
745	+	for (i = mtx_ncols(mig); i--; )
746	+	if (dst_rem[i] > minamt &&
747	+	(cur.d < 0 || price[i] < price[cur.d]))
748	+	cur.d = i;
749	+	if (cur.d < 0)
750	+	return(.0);
751	+	if ((cur.price = price[cur.d]) >= best.price)
752	+	continue;                   /* no point checking further */
753	+	cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
754	+	src_rem[cur.s] : dst_rem[cur.d];
755	+	if (cur.amt > maxamt) cur.amt = maxamt;
756	+	dst_rem[cur.d] -= cur.amt;          /* add up differential costs */
757	+	for (i = mtx_nrows(mig); i--; )
758	+	if (i != cur.s)
759	+	cost_others += min_cost(src_rem[i], dst_rem,
760	+	price, mtx_ncols(mig))
761	+	- src_cost[i];
762	+	dst_rem[cur.d] += cur.amt;          /* undo trial move */
763	+	cur.price += cost_others/cur.amt;   /* adjust effective price */
764	+	if (cur.price < best.price)         /* are we better than best? */
765	+	best = cur;
766	+	}
767	+	if ((best.s < 0) | (best.d < 0))
768	+	return(.0);
769	+	/* make the actual move */
770	+	mig->mtx[mtx_ndx(mig,best.s,best.d)] += best.amt;
771	+	src_rem[best.s] -= best.amt;
772	+	dst_rem[best.d] -= best.amt;
773	+	return(best.amt);
774	+	}
775	+
776	+	#ifdef DEBUG
777	+	static char *
778	+	thetaphi(const FVECT v)
779	+	{
780	+	static char     buf[128];
781	+	double          theta, phi;
782	+
783	+	theta = 180./M_PI*acos(v[2]);
784	+	phi = 180./M_PI*atan2(v[1],v[0]);
785	+	sprintf(buf, "(%.0f,%.0f)", theta, phi);
786	+
787	+	return(buf);
788	+	}
789	+	#endif
790	+
791	+	/* Create a new migration holder (sharing memory for multiprocessing) */
792	+	static MIGRATION *
793	+	new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
794	+	{
795	+	size_t          memlen = sizeof(MIGRATION) +
796	+	sizeof(float)*(from_rbf->nrbf*to_rbf->nrbf - 1);
797	+	MIGRATION       *newmig;
798	+	#ifdef _WIN32
799	+	newmig = (MIGRATION *)malloc(memlen);
800	+	#else
801	+	if (nprocs <= 1) {                      /* single process? */
802	+	newmig = (MIGRATION *)malloc(memlen);
803	+	} else {                                /* else need to share memory */
804	+	newmig = (MIGRATION *)mmap(NULL, memlen, PROT_READ|PROT_WRITE,
805	+	MAP_ANON|MAP_SHARED, -1, 0);
806	+	if ((void *)newmig == MAP_FAILED)
807	+	newmig = NULL;
808	+	}
809	+	#endif
810	+	if (newmig == NULL) {
811	+	fprintf(stderr, "%s: cannot allocate new migration\n", progname);
812	+	exit(1);
813	+	}
814	+	newmig->rbfv[0] = from_rbf;
815	+	newmig->rbfv[1] = to_rbf;
816	+	/* insert in edge lists */
817	+	newmig->enxt[0] = from_rbf->ejl;
818	+	from_rbf->ejl = newmig;
819	+	newmig->enxt[1] = to_rbf->ejl;
820	+	to_rbf->ejl = newmig;
821	+	newmig->next = mig_list;                /* push onto global list */
822	+	return(mig_list = newmig);
823	+	}
824	+
825	+	#ifdef _WIN32
826	+	#define await_children(n)       (void)(n)
827	+	#define run_subprocess()        0
828	+	#define end_subprocess()        (void)0
829	+	#else
830	+
831	+	/* Wait for the specified number of child processes to complete */
832	+	static void
833	+	await_children(int n)
834	+	{
835	+	int     exit_status = 0;
836	+
837	+	if (n > nchild)
838	+	n = nchild;
839	+	while (n-- > 0) {
840	+	int     status;
841	+	if (wait(&status) < 0) {
842	+	fprintf(stderr, "%s: missing child(ren)!\n", progname);
843	+	nchild = 0;
844	+	break;
845	+	}
846	+	--nchild;
847	+	if (status) {                   /* something wrong */
848	+	if ((status = WEXITSTATUS(status)))
849	+	exit_status = status;
850	+	else
851	+	exit_status += !exit_status;
852	+	fprintf(stderr, "%s: subprocess died\n", progname);
853	+	n = nchild;             /* wait for the rest */
854	+	}
855	+	}
856	+	if (exit_status)
857	+	exit(exit_status);
858	+	}
859	+
860	+	/* Start child process if multiprocessing selected */
861	+	static pid_t
862	+	run_subprocess(void)
863	+	{
864	+	int     status;
865	+	pid_t   pid;
866	+
867	+	if (nprocs <= 1)                        /* any children requested? */
868	+	return(0);
869	+	await_children(nchild + 1 - nprocs);    /* free up child process */
870	+	if ((pid = fork())) {
871	+	if (pid < 0) {
872	+	fprintf(stderr, "%s: cannot fork subprocess\n",
873	+	progname);
874	+	exit(1);
875	+	}
876	+	++nchild;                       /* subprocess started */
877	+	return(pid);
878	+	}
879	+	nchild = -1;
880	+	return(0);                              /* put child to work */
881	+	}
882	+
883	+	/* If we are in subprocess, call exit */
884	+	#define end_subprocess()        if (nchild < 0) _exit(0); else
885	+
886	+	#endif  /* ! _WIN32 */
887	+
888	+	/* Compute and insert migration along directed edge (may fork child) */
889	+	static MIGRATION *
890	+	create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
891	+	{
892	+	const double    end_thresh = 0.1/(from_rbf->nrbf*to_rbf->nrbf);
893	+	const double    check_thresh = 0.01;
894	+	const double    rel_thresh = 5e-6;
895	+	float           *pmtx;
896	+	MIGRATION       *newmig;
897	+	double          *src_rem, *dst_rem;
898	+	double          total_rem = 1., move_amt;
899	+	int             i;
900	+	/* check if exists already */
901	+	for (newmig = from_rbf->ejl; newmig != NULL;
902	+	newmig = nextedge(from_rbf,newmig))
903	+	if (newmig->rbfv[1] == to_rbf)
904	+	return(NULL);
905	+	/* else allocate */
906	+	newmig = new_migration(from_rbf, to_rbf);
907	+	if (run_subprocess())
908	+	return(newmig);                 /* child continues */
909	+	pmtx = price_routes(from_rbf, to_rbf);
910	+	src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
911	+	dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
912	+	if ((src_rem == NULL) | (dst_rem == NULL)) {
913	+	fprintf(stderr, "%s: Out of memory in create_migration()\n",
914	+	progname);
915	+	exit(1);
916	+	}
917	+	#ifdef DEBUG
918	+	fprintf(stderr, "Building path from (theta,phi) %s ",
919	+	thetaphi(from_rbf->invec));
920	+	fprintf(stderr, "to %s", thetaphi(to_rbf->invec));
921	+	/* if (nchild) */ fputc('\n', stderr);
922	+	#endif
923	+	/* starting quantities */
924	+	memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
925	+	for (i = from_rbf->nrbf; i--; )
926	+	src_rem[i] = rbf_volume(&from_rbf->rbfa[i]) / from_rbf->vtotal;
927	+	for (i = to_rbf->nrbf; i--; )
928	+	dst_rem[i] = rbf_volume(&to_rbf->rbfa[i]) / to_rbf->vtotal;
929	+	do {                                    /* move a bit at a time */
930	+	move_amt = migration_step(newmig, src_rem, dst_rem, pmtx);
931	+	total_rem -= move_amt;
932	+	#ifdef DEBUG
933	+	if (!nchild)
934	+	/* fputc('.', stderr); */
935	+	fprintf(stderr, "%.9f remaining...\r", total_rem);
936	+	#endif
937	+	} while (total_rem > end_thresh && (total_rem > check_thresh) |
938	+	(move_amt > rel_thresh*total_rem));
939	+	#ifdef DEBUG
940	+	if (!nchild) fputs("\ndone.\n", stderr);
941	+	else fprintf(stderr, "finished with %.9f remaining\n", total_rem);
942	+	#endif
943	+	for (i = from_rbf->nrbf; i--; ) {       /* normalize final matrix */
944	+	float       nf = rbf_volume(&from_rbf->rbfa[i]);
945	+	int         j;
946	+	if (nf <= FTINY) continue;
947	+	nf = from_rbf->vtotal / nf;
948	+	for (j = to_rbf->nrbf; j--; )
949	+	newmig->mtx[mtx_ndx(newmig,i,j)] *= nf;
950	+	}
951	+	end_subprocess();                       /* exit here if subprocess */
952	+	free(pmtx);                             /* free working arrays */
953	+	free(src_rem);
954	+	free(dst_rem);
955	+	return(newmig);
956	+	}
957	+
958	+	/* Get triangle surface orientation (unnormalized) */
959	+	static void
960	+	tri_orient(FVECT vres, const FVECT v1, const FVECT v2, const FVECT v3)
961	+	{
962	+	FVECT   v2minus1, v3minus2;
963	+
964	+	VSUB(v2minus1, v2, v1);
965	+	VSUB(v3minus2, v3, v2);
966	+	VCROSS(vres, v2minus1, v3minus2);
967	+	}
968	+
969	+	/* Determine if vertex order is reversed (inward normal) */
970	+	static int
971	+	is_rev_tri(const FVECT v1, const FVECT v2, const FVECT v3)
972	+	{
973	+	FVECT   tor;
974	+
975	+	tri_orient(tor, v1, v2, v3);
976	+
977	+	return(DOT(tor, v2) < 0.);
978	+	}
979	+
980	+	/* Find vertices completing triangles on either side of the given edge */
981	+	static int
982	+	get_triangles(RBFNODE *rbfv[2], const MIGRATION *mig)
983	+	{
984	+	const MIGRATION *ej, *ej2;
985	+	RBFNODE         *tv;
986	+
987	+	rbfv[0] = rbfv[1] = NULL;
988	+	if (mig == NULL)
989	+	return(0);
990	+	for (ej = mig->rbfv[0]->ejl; ej != NULL;
991	+	ej = nextedge(mig->rbfv[0],ej)) {
992	+	if (ej == mig)
993	+	continue;
994	+	tv = opp_rbf(mig->rbfv[0],ej);
995	+	for (ej2 = tv->ejl; ej2 != NULL; ej2 = nextedge(tv,ej2))
996	+	if (opp_rbf(tv,ej2) == mig->rbfv[1]) {
997	+	rbfv[is_rev_tri(mig->rbfv[0]->invec,
998	+	mig->rbfv[1]->invec,
999	+	tv->invec)] = tv;
1000	+	break;
1001	+	}
1002	+	}
1003	+	return((rbfv[0] != NULL) + (rbfv[1] != NULL));
1004	+	}
1005	+
1006	+	/* Check if prospective vertex would create overlapping triangle */
1007	+	static int
1008	+	overlaps_tri(const RBFNODE *bv0, const RBFNODE *bv1, const RBFNODE *pv)
1009	+	{
1010	+	const MIGRATION *ej;
1011	+	RBFNODE         *vother[2];
1012	+	int             im_rev;
1013	+	/* find shared edge in mesh */
1014	+	for (ej = pv->ejl; ej != NULL; ej = nextedge(pv,ej)) {
1015	+	const RBFNODE   *tv = opp_rbf(pv,ej);
1016	+	if (tv == bv0) {
1017	+	im_rev = is_rev_tri(ej->rbfv[0]->invec,
1018	+	ej->rbfv[1]->invec, bv1->invec);
1019	+	break;
1020	+	}
1021	+	if (tv == bv1) {
1022	+	im_rev = is_rev_tri(ej->rbfv[0]->invec,
1023	+	ej->rbfv[1]->invec, bv0->invec);
1024	+	break;
1025	+	}
1026	+	}
1027	+	if (!get_triangles(vother, ej))         /* triangle on same side? */
1028	+	return(0);
1029	+	return(vother[im_rev] != NULL);
1030	+	}
1031	+
1032	+	/* Find context hull vertex to complete triangle (oriented call) */
1033	+	static RBFNODE *
1034	+	find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rbf1)
1035	+	{
1036	+	FVECT   vmid, vejn, vp;
1037	+	RBFNODE *rbf, *rbfbest = NULL;
1038	+	double  dprod, area2, bestarea2 = FHUGE, bestdprod = -.5;
1039	+
1040	+	VSUB(vejn, rbf1->invec, rbf0->invec);
1041	+	VADD(vmid, rbf0->invec, rbf1->invec);
1042	+	if (normalize(vejn) == 0 || normalize(vmid) == 0)
1043	+	return(NULL);
1044	+	/* XXX exhaustive search */
1045	+	/* Find triangle with minimum rotation from perpendicular */
1046	+	for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
1047	+	if ((rbf == rbf0) | (rbf == rbf1))
1048	+	continue;
1049	+	tri_orient(vp, rbf0->invec, rbf1->invec, rbf->invec);
1050	+	if (DOT(vp, vmid) <= FTINY)
1051	+	continue;               /* wrong orientation */
1052	+	area2 = .25*DOT(vp,vp);
1053	+	VSUB(vp, rbf->invec, rbf0->invec);
1054	+	dprod = -DOT(vp, vejn);
1055	+	VSUM(vp, vp, vejn, dprod);      /* above guarantees non-zero */
1056	+	dprod = DOT(vp, vmid) / VLEN(vp);
1057	+	if (dprod <= bestdprod + FTINY*(1 - 2*(area2 < bestarea2)))
1058	+	continue;               /* found better already */
1059	+	if (overlaps_tri(rbf0, rbf1, rbf))
1060	+	continue;               /* overlaps another triangle */
1061	+	rbfbest = rbf;
1062	+	bestdprod = dprod;              /* new one to beat */
1063	+	bestarea2 = area2;
1064	+	}
1065	+	return(rbfbest);
1066	+	}
1067	+
1068	+	/* Create new migration edge and grow mesh recursively around it */
1069	+	static void
1070	+	mesh_from_edge(MIGRATION *edge)
1071	+	{
1072	+	MIGRATION       *ej0, *ej1;
1073	+	RBFNODE         *tvert[2];
1074	+
1075	+	if (edge == NULL)
1076	+	return;
1077	+	/* triangle on either side? */
1078	+	get_triangles(tvert, edge);
1079	+	if (tvert[0] == NULL) {                 /* grow mesh on right */
1080	+	tvert[0] = find_chull_vert(edge->rbfv[0], edge->rbfv[1]);
1081	+	if (tvert[0] != NULL) {
1082	+	if (tvert[0] > edge->rbfv[0])
1083	+	ej0 = create_migration(edge->rbfv[0], tvert[0]);
1084	+	else
1085	+	ej0 = create_migration(tvert[0], edge->rbfv[0]);
1086	+	if (tvert[0] > edge->rbfv[1])
1087	+	ej1 = create_migration(edge->rbfv[1], tvert[0]);
1088	+	else
1089	+	ej1 = create_migration(tvert[0], edge->rbfv[1]);
1090	+	mesh_from_edge(ej0);
1091	+	mesh_from_edge(ej1);
1092	+	}
1093	+	} else if (tvert[1] == NULL) {          /* grow mesh on left */
1094	+	tvert[1] = find_chull_vert(edge->rbfv[1], edge->rbfv[0]);
1095	+	if (tvert[1] != NULL) {
1096	+	if (tvert[1] > edge->rbfv[0])
1097	+	ej0 = create_migration(edge->rbfv[0], tvert[1]);
1098	+	else
1099	+	ej0 = create_migration(tvert[1], edge->rbfv[0]);
1100	+	if (tvert[1] > edge->rbfv[1])
1101	+	ej1 = create_migration(edge->rbfv[1], tvert[1]);
1102	+	else
1103	+	ej1 = create_migration(tvert[1], edge->rbfv[1]);
1104	+	mesh_from_edge(ej0);
1105	+	mesh_from_edge(ej1);
1106	+	}
1107	+	}
1108	+	}
1109	+
1110	+	#ifdef DEBUG
1111	+	#include "random.h"
1112	+	#include "bmpfile.h"
1113	+	/* Hash pointer to byte value (must return 0 for NULL) */
1114	+	static int
1115	+	byte_hash(const void *p)
1116	+	{
1117	+	size_t  h = (size_t)p;
1118	+	h ^= (size_t)p >> 8;
1119	+	h ^= (size_t)p >> 16;
1120	+	h ^= (size_t)p >> 24;
1121	+	return(h & 0xff);
1122	+	}
1123	+	/* Write out BMP image showing edges */
1124	+	static void
1125	+	write_edge_image(const char *fname)
1126	+	{
1127	+	BMPHeader       *hdr = BMPmappedHeader(GRIDRES, GRIDRES, 0, 256);
1128	+	BMPWriter       *wtr;
1129	+	int             i, j;
1130	+
1131	+	fprintf(stderr, "Writing incident mesh drawing to '%s'\n", fname);
1132	+	hdr->compr = BI_RLE8;
1133	+	for (i = 256; --i; ) {                  /* assign random color map */
1134	+	hdr->palette[i].r = random() & 0xff;
1135	+	hdr->palette[i].g = random() & 0xff;
1136	+	hdr->palette[i].b = random() & 0xff;
1137	+	/* reject dark colors */
1138	+	i += (hdr->palette[i].r + hdr->palette[i].g +
1139	+	hdr->palette[i].b < 128);
1140	+	}
1141	+	hdr->palette[0].r = hdr->palette[0].g = hdr->palette[0].b = 0;
1142	+	/* open output */
1143	+	wtr = BMPopenOutputFile(fname, hdr);
1144	+	if (wtr == NULL) {
1145	+	free(hdr);
1146	+	return;
1147	+	}
1148	+	for (i = 0; i < GRIDRES; i++) {         /* write scanlines */
1149	+	for (j = 0; j < GRIDRES; j++)
1150	+	wtr->scanline[j] = byte_hash(mig_grid[i][j]);
1151	+	if (BMPwriteScanline(wtr) != BIR_OK)
1152	+	break;
1153	+	}
1154	+	BMPcloseOutput(wtr);                    /* close & clean up */
1155	+	}
1156	+	#endif
1157	+
1158	+	/* Draw edge list into mig_grid array */
1159	+	static void
1160	+	draw_edges()
1161	+	{
1162	+	int             nnull = 0, ntot = 0;
1163	+	MIGRATION       *ej;
1164	+	int             p0[2], p1[2];
1165	+
1166	+	/* memset(mig_grid, 0, sizeof(mig_grid)); */
1167	+	for (ej = mig_list; ej != NULL; ej = ej->next) {
1168	+	++ntot;
1169	+	pos_from_vec(p0, ej->rbfv[0]->invec);
1170	+	pos_from_vec(p1, ej->rbfv[1]->invec);
1171	+	if ((p0[0] == p1[0]) & (p0[1] == p1[1])) {
1172	+	++nnull;
1173	+	mig_grid[p0[0]][p0[1]] = ej;
1174	+	continue;
1175	+	}
1176	+	if (abs(p1[0]-p0[0]) > abs(p1[1]-p0[1])) {
1177	+	const int       xstep = 2*(p1[0] > p0[0]) - 1;
1178	+	const double    ystep = (double)((p1[1]-p0[1])*xstep) /
1179	+	(double)(p1[0]-p0[0]);
1180	+	int             x;
1181	+	double          y;
1182	+	for (x = p0[0], y = p0[1]+.5; x != p1[0];
1183	+	x += xstep, y += ystep)
1184	+	mig_grid[x][(int)y] = ej;
1185	+	mig_grid[x][(int)y] = ej;
1186	+	} else {
1187	+	const int       ystep = 2*(p1[1] > p0[1]) - 1;
1188	+	const double    xstep = (double)((p1[0]-p0[0])*ystep) /
1189	+	(double)(p1[1]-p0[1]);
1190	+	int             y;
1191	+	double          x;
1192	+	for (y = p0[1], x = p0[0]+.5; y != p1[1];
1193	+	y += ystep, x += xstep)
1194	+	mig_grid[(int)x][y] = ej;
1195	+	mig_grid[(int)x][y] = ej;
1196	+	}
1197	+	}
1198	+	if (nnull)
1199	+	fprintf(stderr, "Warning: %d of %d edges are null\n",
1200	+	nnull, ntot);
1201	+	#ifdef DEBUG
1202	+	write_edge_image("bsdf_edges.bmp");
1203	+	#endif
1204	+	}
1205	+
1206	+	/* Build our triangle mesh from recorded RBFs */
1207	+	static void
1208	+	build_mesh()
1209	+	{
1210	+	double          best2 = M_PI*M_PI;
1211	+	RBFNODE         *shrt_edj[2];
1212	+	RBFNODE         *rbf0, *rbf1;
1213	+	/* check if isotropic */
1214	+	if (single_plane_incident) {
1215	+	for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
1216	+	if (rbf0->next != NULL)
1217	+	create_migration(rbf0, rbf0->next);
1218	+	await_children(nchild);
1219	+	return;
1220	+	}
1221	+	shrt_edj[0] = shrt_edj[1] = NULL;       /* start w/ shortest edge */
1222	+	for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
1223	+	for (rbf1 = rbf0->next; rbf1 != NULL; rbf1 = rbf1->next) {
1224	+	double  dist2 = 2. - 2.*DOT(rbf0->invec,rbf1->invec);
1225	+	if (dist2 < best2) {
1226	+	shrt_edj[0] = rbf0;
1227	+	shrt_edj[1] = rbf1;
1228	+	best2 = dist2;
1229	+	}
1230	+	}
1231	+	if (shrt_edj[0] == NULL) {
1232	+	fprintf(stderr, "%s: Cannot find shortest edge\n", progname);
1233	+	exit(1);
1234	+	}
1235	+	/* build mesh from this edge */
1236	+	if (shrt_edj[0] < shrt_edj[1])
1237	+	mesh_from_edge(create_migration(shrt_edj[0], shrt_edj[1]));
1238	+	else
1239	+	mesh_from_edge(create_migration(shrt_edj[1], shrt_edj[0]));
1240	+	/* draw edge list into grid */
1241	+	draw_edges();
1242	+	/* complete migrations */
1243	+	await_children(nchild);
1244	+	}
1245	+
1246	+	/* Identify enclosing triangle for this position (flood fill raster check) */
1247	+	static int
1248	+	identify_tri(MIGRATION *miga[3], unsigned char vmap[GRIDRES][(GRIDRES+7)/8],
1249	+	int px, int py)
1250	+	{
1251	+	const int       btest = 1<<(py&07);
1252	+
1253	+	if (vmap[px][py>>3] & btest)            /* already visited here? */
1254	+	return(1);
1255	+	/* else mark it */
1256	+	vmap[px][py>>3] |= btest;
1257	+
1258	+	if (mig_grid[px][py] != NULL) {         /* are we on an edge? */
1259	+	int     i;
1260	+	for (i = 0; i < 3; i++) {
1261	+	if (miga[i] == mig_grid[px][py])
1262	+	return(1);
1263	+	if (miga[i] != NULL)
1264	+	continue;
1265	+	miga[i] = mig_grid[px][py];
1266	+	return(1);
1267	+	}
1268	+	return(0);                      /* outside triangle! */
1269	+	}
1270	+	/* check neighbors (flood) */
1271	+	if (px > 0 && !identify_tri(miga, vmap, px-1, py))
1272	+	return(0);
1273	+	if (px < GRIDRES-1 && !identify_tri(miga, vmap, px+1, py))
1274	+	return(0);
1275	+	if (py > 0 && !identify_tri(miga, vmap, px, py-1))
1276	+	return(0);
1277	+	if (py < GRIDRES-1 && !identify_tri(miga, vmap, px, py+1))
1278	+	return(0);
1279	+	return(1);                              /* this neighborhood done */
1280	+	}
1281	+
1282	+	/* Insert vertex in ordered list */
1283	+	static void
1284	+	insert_vert(RBFNODE **vlist, RBFNODE *v)
1285	+	{
1286	+	int     i, j;
1287	+
1288	+	for (i = 0; vlist[i] != NULL; i++) {
1289	+	if (v == vlist[i])
1290	+	return;
1291	+	if (v < vlist[i])
1292	+	break;
1293	+	}
1294	+	for (j = i; vlist[j] != NULL; j++)
1295	+	;
1296	+	while (j > i) {
1297	+	vlist[j] = vlist[j-1];
1298	+	--j;
1299	+	}
1300	+	vlist[i] = v;
1301	+	}
1302	+
1303	+	/* Sort triangle edges in standard order */
1304	+	static int
1305	+	order_triangle(MIGRATION *miga[3])
1306	+	{
1307	+	RBFNODE         *vert[7];
1308	+	MIGRATION       *ord[3];
1309	+	int             i;
1310	+	/* order vertices, first */
1311	+	memset(vert, 0, sizeof(vert));
1312	+	for (i = 3; i--; ) {
1313	+	if (miga[i] == NULL)
1314	+	return(0);
1315	+	insert_vert(vert, miga[i]->rbfv[0]);
1316	+	insert_vert(vert, miga[i]->rbfv[1]);
1317	+	}
1318	+	/* should be just 3 vertices */
1319	+	if ((vert[3] == NULL) | (vert[4] != NULL))
1320	+	return(0);
1321	+	/* identify edge 0 */
1322	+	for (i = 3; i--; )
1323	+	if (miga[i]->rbfv[0] == vert[0] &&
1324	+	miga[i]->rbfv[1] == vert[1]) {
1325	+	ord[0] = miga[i];
1326	+	break;
1327	+	}
1328	+	if (i < 0)
1329	+	return(0);
1330	+	/* identify edge 1 */
1331	+	for (i = 3; i--; )
1332	+	if (miga[i]->rbfv[0] == vert[1] &&
1333	+	miga[i]->rbfv[1] == vert[2]) {
1334	+	ord[1] = miga[i];
1335	+	break;
1336	+	}
1337	+	if (i < 0)
1338	+	return(0);
1339	+	/* identify edge 2 */
1340	+	for (i = 3; i--; )
1341	+	if (miga[i]->rbfv[0] == vert[0] &&
1342	+	miga[i]->rbfv[1] == vert[2]) {
1343	+	ord[2] = miga[i];
1344	+	break;
1345	+	}
1346	+	if (i < 0)
1347	+	return(0);
1348	+	/* reassign order */
1349	+	miga[0] = ord[0]; miga[1] = ord[1]; miga[2] = ord[2];
1350	+	return(1);
1351	+	}
1352	+
1353	+	/* Find edge(s) for interpolating the given vector, applying symmetry */
1354	+	static int
1355	+	get_interp(MIGRATION *miga[3], FVECT invec)
1356	+	{
1357	+	miga[0] = miga[1] = miga[2] = NULL;
1358	+	if (single_plane_incident) {            /* isotropic BSDF? */
1359	+	RBFNODE *rbf;                   /* find edge we're on */
1360	+	for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
1361	+	if (input_orient*rbf->invec[2] < input_orient*invec[2])
1362	+	break;
1363	+	if (rbf->next != NULL &&
1364	+	input_orient*rbf->next->invec[2] <
1365	+	input_orient*invec[2]) {
1366	+	for (miga[0] = rbf->ejl; miga[0] != NULL;
1367	+	miga[0] = nextedge(rbf,miga[0]))
1368	+	if (opp_rbf(rbf,miga[0]) == rbf->next)
1369	+	return(0);
1370	+	break;
1371	+	}
1372	+	}
1373	+	return(-1);                     /* outside range! */
1374	+	}
1375	+	{                                       /* else use triangle mesh */
1376	+	const int       sym = use_symmetry(invec);
1377	+	unsigned char   floodmap[GRIDRES][(GRIDRES+7)/8];
1378	+	int             pstart[2];
1379	+	RBFNODE         *vother;
1380	+	MIGRATION       *ej;
1381	+	int             i;
1382	+
1383	+	pos_from_vec(pstart, invec);
1384	+	memset(floodmap, 0, sizeof(floodmap));
1385	+	/* call flooding function */
1386	+	if (!identify_tri(miga, floodmap, pstart[0], pstart[1]))
1387	+	return(-1);             /* outside mesh */
1388	+	if ((miga[0] == NULL) | (miga[2] == NULL))
1389	+	return(-1);             /* should never happen */
1390	+	if (miga[1] == NULL)
1391	+	return(sym);            /* on edge */
1392	+	/* verify triangle */
1393	+	if (!order_triangle(miga)) {
1394	+	#ifdef DEBUG
1395	+	fputs("Munged triangle in get_interp()\n", stderr);
1396	+	#endif
1397	+	vother = NULL;          /* find triangle from edge */
1398	+	for (i = 3; i--; ) {
1399	+	RBFNODE     *tpair[2];
1400	+	if (get_triangles(tpair, miga[i]) &&
1401	+	(vother = tpair[ is_rev_tri(
1402	+	miga[i]->rbfv[0]->invec,
1403	+	miga[i]->rbfv[1]->invec,
1404	+	invec) ]) != NULL)
1405	+	break;
1406	+	}
1407	+	if (vother == NULL) {   /* couldn't find 3rd vertex */
1408	+	#ifdef DEBUG
1409	+	fputs("No triangle in get_interp()\n", stderr);
1410	+	#endif
1411	+	return(-1);
1412	+	}
1413	+	/* reassign other two edges */
1414	+	for (ej = vother->ejl; ej != NULL;
1415	+	ej = nextedge(vother,ej)) {
1416	+	RBFNODE *vorig = opp_rbf(vother,ej);
1417	+	if (vorig == miga[i]->rbfv[0])
1418	+	miga[(i+1)%3] = ej;
1419	+	else if (vorig == miga[i]->rbfv[1])
1420	+	miga[(i+2)%3] = ej;
1421	+	}
1422	+	if (!order_triangle(miga)) {
1423	+	#ifdef DEBUG
1424	+	fputs("Bad triangle in get_interp()\n", stderr);
1425	+	#endif
1426	+	return(-1);
1427	+	}
1428	+	}
1429	+	return(sym);                    /* return in standard order */
1430	+	}
1431	+	}
1432	+
1433	+	/* Advect and allocate new RBF along edge */
1434	+	static RBFNODE *
1435	+	e_advect_rbf(const MIGRATION *mig, const FVECT invec)
1436	+	{
1437	+	RBFNODE         *rbf;
1438	+	int             n, i, j;
1439	+	double          t, full_dist;
1440	+	/* get relative position */
1441	+	t = acos(DOT(invec, mig->rbfv[0]->invec));
1442	+	if (t < M_PI/GRIDRES) {                 /* near first DSF */
1443	+	n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[0]->nrbf-1);
1444	+	rbf = (RBFNODE *)malloc(n);
1445	+	if (rbf == NULL)
1446	+	goto memerr;
1447	+	memcpy(rbf, mig->rbfv[0], n);   /* just duplicate */
1448	+	return(rbf);
1449	+	}
1450	+	full_dist = acos(DOT(mig->rbfv[0]->invec, mig->rbfv[1]->invec));
1451	+	if (t > full_dist-M_PI/GRIDRES) {       /* near second DSF */
1452	+	n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[1]->nrbf-1);
1453	+	rbf = (RBFNODE *)malloc(n);
1454	+	if (rbf == NULL)
1455	+	goto memerr;
1456	+	memcpy(rbf, mig->rbfv[1], n);   /* just duplicate */
1457	+	return(rbf);
1458	+	}
1459	+	t /= full_dist;
1460	+	n = 0;                                  /* count migrating particles */
1461	+	for (i = 0; i < mtx_nrows(mig); i++)
1462	+	for (j = 0; j < mtx_ncols(mig); j++)
1463	+	n += (mig->mtx[mtx_ndx(mig,i,j)] > FTINY);
1464	+	#ifdef DEBUG
1465	+	fprintf(stderr, "Input RBFs have %d, %d nodes -> output has %d\n",
1466	+	mig->rbfv[0]->nrbf, mig->rbfv[1]->nrbf, n);
1467	+	#endif
1468	+	rbf = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(n-1));
1469	+	if (rbf == NULL)
1470	+	goto memerr;
1471	+	rbf->next = NULL; rbf->ejl = NULL;
1472	+	VCOPY(rbf->invec, invec);
1473	+	rbf->nrbf = n;
1474	+	rbf->vtotal = 1.-t + t*mig->rbfv[1]->vtotal/mig->rbfv[0]->vtotal;
1475	+	n = 0;                                  /* advect RBF lobes */
1476	+	for (i = 0; i < mtx_nrows(mig); i++) {
1477	+	const RBFVAL        *rbf0i = &mig->rbfv[0]->rbfa[i];
1478	+	const float         peak0 = rbf0i->peak;
1479	+	const double        rad0 = R2ANG(rbf0i->crad);
1480	+	FVECT               v0;
1481	+	float               mv;
1482	+	ovec_from_pos(v0, rbf0i->gx, rbf0i->gy);
1483	+	for (j = 0; j < mtx_ncols(mig); j++)
1484	+	if ((mv = mig->mtx[mtx_ndx(mig,i,j)]) > FTINY) {
1485	+	const RBFVAL    *rbf1j = &mig->rbfv[1]->rbfa[j];
1486	+	double          rad1 = R2ANG(rbf1j->crad);
1487	+	FVECT           v;
1488	+	int             pos[2];
1489	+	rbf->rbfa[n].peak = peak0 * mv * rbf->vtotal;
1490	+	rbf->rbfa[n].crad = ANG2R(sqrt(rad0*rad0*(1.-t) +
1491	+	rad1*rad1*t));
1492	+	ovec_from_pos(v, rbf1j->gx, rbf1j->gy);
1493	+	geodesic(v, v0, v, t, GEOD_REL);
1494	+	pos_from_vec(pos, v);
1495	+	rbf->rbfa[n].gx = pos[0];
1496	+	rbf->rbfa[n].gy = pos[1];
1497	+	++n;
1498	+	}
1499	+	}
1500	+	rbf->vtotal *= mig->rbfv[0]->vtotal;    /* turn ratio into actual */
1501	+	return(rbf);
1502	+	memerr:
1503	+	fprintf(stderr, "%s: Out of memory in e_advect_rbf()\n", progname);
1504	+	exit(1);
1505	+	return(NULL);   /* pro forma return */
1506	+	}
1507	+
1508	+	/* Partially advect between recorded incident angles and allocate new RBF */
1509	+	static RBFNODE *
1510	+	advect_rbf(const FVECT invec)
1511	+	{
1512	+	FVECT           sivec;
1513	+	MIGRATION       *miga[3];
1514	+	RBFNODE         *rbf;
1515	+	int             sym;
1516	+	float           mbfact, mcfact;
1517	+	int             n, i, j, k;
1518	+	FVECT           v0, v1, v2;
1519	+	double          s, t;
1520	+
1521	+	VCOPY(sivec, invec);                    /* find triangle/edge */
1522	+	sym = get_interp(miga, sivec);
1523	+	if (sym < 0)                            /* can't interpolate? */
1524	+	return(NULL);
1525	+	if (miga[1] == NULL) {                  /* advect along edge? */
1526	+	rbf = e_advect_rbf(miga[0], sivec);
1527	+	rev_rbf_symmetry(rbf, sym);
1528	+	return(rbf);
1529	+	}
1530	+	#ifdef DEBUG
1531	+	if (miga[0]->rbfv[0] != miga[2]->rbfv[0] |
1532	+	miga[0]->rbfv[1] != miga[1]->rbfv[0] |
1533	+	miga[1]->rbfv[1] != miga[2]->rbfv[1]) {
1534	+	fprintf(stderr, "%s: Triangle vertex screw-up!\n", progname);
1535	+	exit(1);
1536	+	}
1537	+	#endif
1538	+	/* figure out position */
1539	+	fcross(v0, miga[2]->rbfv[0]->invec, miga[2]->rbfv[1]->invec);
1540	+	normalize(v0);
1541	+	fcross(v2, miga[1]->rbfv[0]->invec, miga[1]->rbfv[1]->invec);
1542	+	normalize(v2);
1543	+	fcross(v1, sivec, miga[1]->rbfv[1]->invec);
1544	+	normalize(v1);
1545	+	s = acos(DOT(v0,v1)) / acos(DOT(v0,v2));
1546	+	geodesic(v1, miga[0]->rbfv[0]->invec, miga[0]->rbfv[1]->invec,
1547	+	s, GEOD_REL);
1548	+	t = acos(DOT(v1,sivec)) / acos(DOT(v1,miga[1]->rbfv[1]->invec));
1549	+	n = 0;                                  /* count migrating particles */
1550	+	for (i = 0; i < mtx_nrows(miga[0]); i++)
1551	+	for (j = 0; j < mtx_ncols(miga[0]); j++)
1552	+	for (k = (miga[0]->mtx[mtx_ndx(miga[0],i,j)] > FTINY) *
1553	+	mtx_ncols(miga[2]); k--; )
1554	+	n += (miga[2]->mtx[mtx_ndx(miga[2],i,k)] > FTINY &&
1555	+	miga[1]->mtx[mtx_ndx(miga[1],j,k)] > FTINY);
1556	+	#ifdef DEBUG
1557	+	fprintf(stderr, "Input RBFs have %d, %d, %d nodes -> output has %d\n",
1558	+	miga[0]->rbfv[0]->nrbf, miga[0]->rbfv[1]->nrbf,
1559	+	miga[2]->rbfv[1]->nrbf, n);
1560	+	#endif
1561	+	rbf = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(n-1));
1562	+	if (rbf == NULL) {
1563	+	fprintf(stderr, "%s: Out of memory in advect_rbf()\n", progname);
1564	+	exit(1);
1565	+	}
1566	+	rbf->next = NULL; rbf->ejl = NULL;
1567	+	VCOPY(rbf->invec, sivec);
1568	+	rbf->nrbf = n;
1569	+	n = 0;                                  /* compute RBF lobes */
1570	+	mbfact = s * miga[0]->rbfv[1]->vtotal/miga[0]->rbfv[0]->vtotal *
1571	+	(1.-t + t*miga[1]->rbfv[1]->vtotal/miga[1]->rbfv[0]->vtotal);
1572	+	mcfact = (1.-s) *
1573	+	(1.-t + t*miga[2]->rbfv[1]->vtotal/miga[2]->rbfv[0]->vtotal);
1574	+	for (i = 0; i < mtx_nrows(miga[0]); i++) {
1575	+	const RBFVAL        *rbf0i = &miga[0]->rbfv[0]->rbfa[i];
1576	+	const float         w0i = rbf0i->peak;
1577	+	const double        rad0i = R2ANG(rbf0i->crad);
1578	+	ovec_from_pos(v0, rbf0i->gx, rbf0i->gy);
1579	+	for (j = 0; j < mtx_ncols(miga[0]); j++) {
1580	+	const float     ma = miga[0]->mtx[mtx_ndx(miga[0],i,j)];
1581	+	const RBFVAL    *rbf1j;
1582	+	double          rad1j, srad2;
1583	+	if (ma <= FTINY)
1584	+	continue;
1585	+	rbf1j = &miga[0]->rbfv[1]->rbfa[j];
1586	+	rad1j = R2ANG(rbf1j->crad);
1587	+	srad2 = (1.-s)*(1.-t)*rad0i*rad0i + s*(1.-t)*rad1j*rad1j;
1588	+	ovec_from_pos(v1, rbf1j->gx, rbf1j->gy);
1589	+	geodesic(v1, v0, v1, s, GEOD_REL);
1590	+	for (k = 0; k < mtx_ncols(miga[2]); k++) {
1591	+	float               mb = miga[1]->mtx[mtx_ndx(miga[1],j,k)];
1592	+	float               mc = miga[2]->mtx[mtx_ndx(miga[2],i,k)];
1593	+	const RBFVAL        *rbf2k;
1594	+	double              rad2k;
1595	+	FVECT               vout;
1596	+	int                 pos[2];
1597	+	if ((mb <= FTINY) | (mc <= FTINY))
1598	+	continue;
1599	+	rbf2k = &miga[2]->rbfv[1]->rbfa[k];
1600	+	rbf->rbfa[n].peak = w0i * ma * (mb*mbfact + mc*mcfact);
1601	+	rad2k = R2ANG(rbf2k->crad);
1602	+	rbf->rbfa[n].crad = ANG2R(sqrt(srad2 + t*rad2k*rad2k));
1603	+	ovec_from_pos(v2, rbf2k->gx, rbf2k->gy);
1604	+	geodesic(vout, v1, v2, t, GEOD_REL);
1605	+	pos_from_vec(pos, vout);
1606	+	rbf->rbfa[n].gx = pos[0];
1607	+	rbf->rbfa[n].gy = pos[1];
1608	+	++n;
1609	+	}
1610	+	}
1611	+	}
1612	+	rbf->vtotal = miga[0]->rbfv[0]->vtotal * (mbfact + mcfact);
1613	+	rev_rbf_symmetry(rbf, sym);
1614	+	return(rbf);
1615	+	}
1616	+
1617	+	/* Interpolate and output isotropic BSDF data */
1618	+	static void
1619	+	interp_isotropic()
1620	+	{
1621	+	const int       sqres = 1<<samp_order;
1622	+	FILE            *ofp = NULL;
1623	+	char            cmd[128];
1624	+	int             ix, ox, oy;
1625	+	FVECT           ivec, ovec;
1626	+	double          bsdf;
1627	+	#if DEBUG
1628	+	fprintf(stderr, "Writing isotropic order %d ", samp_order);
1629	+	if (pctcull >= 0) fprintf(stderr, "data with %d%% culling\n", pctcull);
1630	+	else fputs("raw data\n", stderr);
1631	+	#endif
1632	+	if (pctcull >= 0) {                     /* begin output */
1633	+	sprintf(cmd, "rttree_reduce -h -a -fd -r 3 -t %d -g %d",
1634	+	pctcull, samp_order);
1635	+	fflush(stdout);
1636	+	ofp = popen(cmd, "w");
1637	+	if (ofp == NULL) {
1638	+	fprintf(stderr, "%s: cannot create pipe to rttree_reduce\n",
1639	+	progname);
1640	+	exit(1);
1641	+	}
1642	+	} else
1643	+	fputs("{\n", stdout);
1644	+	/* run through directions */
1645	+	for (ix = 0; ix < sqres/2; ix++) {
1646	+	RBFNODE *rbf;
1647	+	SDsquare2disk(ivec, (ix+.5)/sqres, .5);
1648	+	ivec[2] = input_orient *
1649	+	sqrt(1. - ivec[0]*ivec[0] - ivec[1]*ivec[1]);
1650	+	rbf = advect_rbf(ivec);
1651	+	for (ox = 0; ox < sqres; ox++)
1652	+	for (oy = 0; oy < sqres; oy++) {
1653	+	SDsquare2disk(ovec, (ox+.5)/sqres, (oy+.5)/sqres);
1654	+	ovec[2] = output_orient *
1655	+	sqrt(1. - ovec[0]*ovec[0] - ovec[1]*ovec[1]);
1656	+	bsdf = eval_rbfrep(rbf, ovec) / fabs(ovec[2]);
1657	+	if (pctcull >= 0)
1658	+	fwrite(&bsdf, sizeof(bsdf), 1, ofp);
1659	+	else
1660	+	printf("\t%.3e\n", bsdf);
1661	+	}
1662	+	free(rbf);
1663	+	}
1664	+	if (pctcull >= 0) {                     /* finish output */
1665	+	if (pclose(ofp)) {
1666	+	fprintf(stderr, "%s: error running '%s'\n",
1667	+	progname, cmd);
1668	+	exit(1);
1669	+	}
1670	+	} else {
1671	+	for (ix = sqres*sqres*sqres/2; ix--; )
1672	+	fputs("\t0\n", stdout);
1673	+	fputs("}\n", stdout);
1674	+	}
1675	+	}
1676	+
1677	+	/* Interpolate and output anisotropic BSDF data */
1678	+	static void
1679	+	interp_anisotropic()
1680	+	{
1681	+	const int       sqres = 1<<samp_order;
1682	+	FILE            *ofp = NULL;
1683	+	char            cmd[128];
1684	+	int             ix, iy, ox, oy;
1685	+	FVECT           ivec, ovec;
1686	+	double          bsdf;
1687	+	#if DEBUG
1688	+	fprintf(stderr, "Writing anisotropic order %d ", samp_order);
1689	+	if (pctcull >= 0) fprintf(stderr, "data with %d%% culling\n", pctcull);
1690	+	else fputs("raw data\n", stderr);
1691	+	#endif
1692	+	if (pctcull >= 0) {                     /* begin output */
1693	+	sprintf(cmd, "rttree_reduce -h -a -fd -r 4 -t %d -g %d",
1694	+	pctcull, samp_order);
1695	+	fflush(stdout);
1696	+	ofp = popen(cmd, "w");
1697	+	if (ofp == NULL) {
1698	+	fprintf(stderr, "%s: cannot create pipe to rttree_reduce\n",
1699	+	progname);
1700	+	exit(1);
1701	+	}
1702	+	} else
1703	+	fputs("{\n", stdout);
1704	+	/* run through directions */
1705	+	for (ix = 0; ix < sqres; ix++)
1706	+	for (iy = 0; iy < sqres; iy++) {
1707	+	RBFNODE *rbf;
1708	+	SDsquare2disk(ivec, (ix+.5)/sqres, (iy+.5)/sqres);
1709	+	ivec[2] = input_orient *
1710	+	sqrt(1. - ivec[0]*ivec[0] - ivec[1]*ivec[1]);
1711	+	rbf = advect_rbf(ivec);
1712	+	for (ox = 0; ox < sqres; ox++)
1713	+	for (oy = 0; oy < sqres; oy++) {
1714	+	SDsquare2disk(ovec, (ox+.5)/sqres, (oy+.5)/sqres);
1715	+	ovec[2] = output_orient *
1716	+	sqrt(1. - ovec[0]*ovec[0] - ovec[1]*ovec[1]);
1717	+	bsdf = eval_rbfrep(rbf, ovec) / fabs(ovec[2]);
1718	+	if (pctcull >= 0)
1719	+	fwrite(&bsdf, sizeof(bsdf), 1, ofp);
1720	+	else
1721	+	printf("\t%.3e\n", bsdf);
1722	+	}
1723	+	free(rbf);
1724	+	}
1725	+	if (pctcull >= 0) {                     /* finish output */
1726	+	if (pclose(ofp)) {
1727	+	fprintf(stderr, "%s: error running '%s'\n",
1728	+	progname, cmd);
1729	+	exit(1);
1730	+	}
1731	+	} else
1732	+	fputs("}\n", stdout);
1733	+	}
1734	+
1735		#if 1
1736	+	/* Read in BSDF files and interpolate as tensor tree representation */
1737	+	int
1738	+	main(int argc, char *argv[])
1739	+	{
1740	+	RBFNODE         *rbf;
1741	+	double          bsdf;
1742	+	int             i;
1743	+
1744	+	progname = argv[0];                     /* get options */
1745	+	while (argc > 2 && argv[1][0] == '-') {
1746	+	switch (argv[1][1]) {
1747	+	case 'n':
1748	+	nprocs = atoi(argv[2]);
1749	+	break;
1750	+	case 't':
1751	+	pctcull = atoi(argv[2]);
1752	+	break;
1753	+	default:
1754	+	goto userr;
1755	+	}
1756	+	argv += 2; argc -= 2;
1757	+	}
1758	+	if (argc < 3)
1759	+	goto userr;
1760	+	#ifdef _WIN32
1761	+	if (nprocs > 1) {
1762	+	fprintf(stderr, "%s: multiprocessing not supported\n",
1763	+	progname);
1764	+	return(1);
1765	+	}
1766	+	#endif
1767	+	for (i = 1; i < argc; i++) {            /* compile measurements */
1768	+	if (!load_pabopto_meas(argv[i]))
1769	+	return(1);
1770	+	compute_radii();
1771	+	cull_values();
1772	+	make_rbfrep();
1773	+	}
1774	+	build_mesh();                           /* create interpolation */
1775	+	/* xml_prologue();                              /* start XML output */
1776	+	if (single_plane_incident)              /* resample dist. */
1777	+	interp_isotropic();
1778	+	else
1779	+	interp_anisotropic();
1780	+	/* xml_epilogue();                              /* finish XML output */
1781	+	return(0);
1782	+	userr:
1783	+	fprintf(stderr,
1784	+	"Usage: %s [-n nprocs][-t pctcull] meas1.dat meas2.dat .. > bsdf.xml\n",
1785	+	progname);
1786	+	return(1);
1787	+	}
1788	+	#else
1789		/* Test main produces a Radiance model from the given input file */
1790		int
1791		main(int argc, char *argv[])
#	Line 355 | Line 1800 | main(int argc, char *argv[])
1800		fprintf(stderr, "Usage: %s input.dat > output.rad\n", argv[0]);
1801		return(1);
1802		}
1803	<	if (!load_bsdf_meas(argv[1]))
1803	>	if (!load_pabopto_meas(argv[1]))
1804		return(1);
360	–	/* produce spheres at meas. */
361	–	puts("void plastic orange\n0\n0\n5 .6 .4 .01 .04 .08\n");
362	–	n = 0;
363	–	for (i = 0; i < GRIDRES; i++)
364	–	for (j = 0; j < GRIDRES; j++)
365	–	if (bsdf_grid[i][j].nval) {
366	–	double  bsdf = bsdf_grid[i][j].vsum /
367	–	(double)bsdf_grid[i][j].nval;
368	–	FVECT   dir;
1805
370	–	vec_from_pos(dir, i, j);
371	–	printf("orange sphere s%04d\n0\n0\n", ++n);
372	–	printf("4 %.6g %.6g %.6g .0015\n\n",
373	–	dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
374	–	}
1806		compute_radii();
1807		cull_values();
1808	<	/* highlight chosen values */
1808	>	make_rbfrep();
1809	>	/* produce spheres at meas. */
1810	>	puts("void plastic yellow\n0\n0\n5 .6 .4 .01 .04 .08\n");
1811		puts("void plastic pink\n0\n0\n5 .5 .05 .9 .04 .08\n");
1812		n = 0;
1813		for (i = 0; i < GRIDRES; i++)
1814		for (j = 0; j < GRIDRES; j++)
1815	<	if (bsdf_grid[i][j].nval) {
1816	<	bsdf = bsdf_grid[i][j].vsum /
1817	<	(double)bsdf_grid[i][j].nval;
1818	<	vec_from_pos(dir, i, j);
1819	<	printf("pink cone c%04d\n0\n0\n8\n", ++n);
1820	<	printf("\t%.6g %.6g %.6g\n",
1815	>	if (dsf_grid[i][j].vsum > .0f) {
1816	>	ovec_from_pos(dir, i, j);
1817	>	bsdf = dsf_grid[i][j].vsum / dir[2];
1818	>	if (dsf_grid[i][j].nval) {
1819	>	printf("pink cone c%04d\n0\n0\n8\n", ++n);
1820	>	printf("\t%.6g %.6g %.6g\n",
1821		dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
1822	<	printf("\t%.6g %.6g %.6g\n",
1822	>	printf("\t%.6g %.6g %.6g\n",
1823		dir[0]*(bsdf+.005), dir[1]*(bsdf+.005),
1824		dir[2]*(bsdf+.005));
1825	<	puts("\t.003\t0\n");
1825	>	puts("\t.003\t0\n");
1826	>	} else {
1827	>	ovec_from_pos(dir, i, j);
1828	>	printf("yellow sphere s%04d\n0\n0\n", ++n);
1829	>	printf("4 %.6g %.6g %.6g .0015\n\n",
1830	>	dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
1831	>	}
1832		}
1833		/* output continuous surface */
395	–	make_rbfrep();
1834		puts("void trans tgreen\n0\n0\n7 .7 1 .7 .04 .04 .9 .9\n");
1835		fflush(stdout);
1836	<	sprintf(buf, "gensurf tgreen bsdf - - - %d %d", GRIDRES, GRIDRES);
1836	>	sprintf(buf, "gensurf tgreen bsdf - - - %d %d", GRIDRES-1, GRIDRES-1);
1837		pfp = popen(buf, "w");
1838		if (pfp == NULL) {
1839		fputs(buf, stderr);
#	Line 404 | Line 1842 | main(int argc, char *argv[])
1842		}
1843		for (i = 0; i < GRIDRES; i++)
1844		for (j = 0; j < GRIDRES; j++) {
1845	<	vec_from_pos(dir, i, j);
1846	<	bsdf = eval_rbfrep(bsdf_list, dir);
1845	>	ovec_from_pos(dir, i, j);
1846	>	bsdf = eval_rbfrep(dsf_list, dir) / dir[2];
1847		fprintf(pfp, "%.8e %.8e %.8e\n",
1848		dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
1849		}

Diff Legend

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