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.1 by greg, Fri Aug 24 15:20:18 2012 UTC vs.
Revision 2.12 by greg, Thu Sep 20 01:23:36 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            1.9             /* radius scaling factor (empirical) */
24	<	#define MSCA            .12             /* magnitude scaling (empirical) */
25	>	#define RSCA            2.7             /* radius scaling factor (empirical) */
26
27	+	/* convert to/from coded radians */
28	+	#define ANG2R(r)        (int)((r)*((1<<16)/M_PI))
29	+	#define R2ANG(c)        (((c)+.5)*(M_PI/(1<<16)))
30	+
31		typedef struct {
32	<	float           vsum;           /* BSDF sum */
32	>	float           vsum;           /* DSF sum */
33		unsigned short  nval;           /* number of values in sum */
34	<	unsigned short  hrad2;          /* half radius squared */
34	>	unsigned short  crad;           /* radius (coded angle) */
35		} GRIDVAL;                      /* grid value */
36
37		typedef struct {
38	<	float           bsdf;           /* BSDF value at peak */
39	<	unsigned short  rad;            /* radius */
38	>	float           peak;           /* lobe value at peak */
39	>	unsigned short  crad;           /* radius (coded angle) */
40		unsigned char   gx, gy;         /* grid position */
41		} RBFVAL;                       /* radial basis function value */
42
43	<	typedef struct s_rbflist {
44	<	struct s_rbflist        *next;          /* next in our RBF list */
43	>	struct s_rbfnode;               /* forward declaration of RBF struct */
44	>
45	>	typedef struct s_migration {
46	>	struct s_migration      *next;          /* next in global edge list */
47	>	struct s_rbfnode        *rbfv[2];       /* from,to vertex */
48	>	struct s_migration      *enxt[2];       /* next from,to sibling */
49	>	float                   mtx[1];         /* matrix (extends struct) */
50	>	} MIGRATION;                    /* migration link (winged edge structure) */
51	>
52	>	typedef struct s_rbfnode {
53	>	struct s_rbfnode        *next;          /* next in global RBF list */
54	>	MIGRATION               *ejl;           /* edge list for this vertex */
55		FVECT                   invec;          /* incident vector direction */
56	+	double                  vtotal;         /* volume for normalization */
57		int                     nrbf;           /* number of RBFs */
58		RBFVAL                  rbfa[1];        /* RBF array (extends struct) */
59	<	} RBFLIST;                      /* RBF representation of BSDF @ 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  bsdf_grid[GRIDRES][GRIDRES];
62	>	static double           theta_in_deg, phi_in_deg;
63	>	static GRIDVAL          dsf_grid[GRIDRES][GRIDRES];
64
65	<	/* processed incident BSDF measurements */
66	<	static RBFLIST  *bsdf_list = NULL;
65	>	/* all incident angles in-plane so far? */
66	>	static int              single_plane_incident = -1;
67
68	<	/* Count up non-empty nodes and build RBF representation from current grid */
69	<	static RBFLIST *
70	<	make_rbfrep(void)
55	<	{
56	<	int     nn = 0;
57	<	RBFLIST *newnode;
58	<	int     i, j;
59	<	/* count non-empty bins */
60	<	for (i = 0; i < GRIDRES; i++)
61	<	for (j = 0; j < GRIDRES; j++)
62	<	nn += (bsdf_grid[i][j].nval > 0);
63	<	/* allocate RBF array */
64	<	newnode = (RBFLIST *)malloc(sizeof(RBFLIST) + sizeof(RBFVAL)*(nn-1));
65	<	if (newnode == NULL) {
66	<	fputs("Out of memory in make_rbfrep\n", stderr);
67	<	exit(1);
68	<	}
69	<	newnode->invec[2] = sin(M_PI/180.*theta_in_deg);
70	<	newnode->invec[0] = cos(M_PI/180.*phi_in_deg)*newnode->invec[2];
71	<	newnode->invec[1] = sin(M_PI/180.*phi_in_deg)*newnode->invec[2];
72	<	newnode->invec[2] = sqrt(1. - newnode->invec[2]*newnode->invec[2]);
73	<	newnode->nrbf = nn;
74	<	nn = 0;                 /* fill RBF array */
75	<	for (i = 0; i < GRIDRES; i++)
76	<	for (j = 0; j < GRIDRES; j++)
77	<	if (bsdf_grid[i][j].nval) {
78	<	newnode->rbfa[nn].bsdf = MSCA*bsdf_grid[i][j].vsum /
79	<	(double)bsdf_grid[i][j].nval;
80	<	newnode->rbfa[nn].rad =
81	<	(int)(2.*RSCA*sqrt((double)bsdf_grid[i][j].hrad2) + .5);
82	<	newnode->rbfa[nn].gx = i;
83	<	newnode->rbfa[nn].gy = j;
84	<	++nn;
85	<	}
86	<	newnode->next = bsdf_list;
87	<	return(bsdf_list = newnode);
88	<	}
68	>	/* input/output orientations */
69	>	static int              input_orient = 0;
70	>	static int              output_orient = 0;
71
72	<	/* Compute grid position from normalized outgoing vector */
73	<	static void
92	<	pos_from_vec(int pos[2], const FVECT vec)
93	<	{
94	<	double  sq[2];          /* uniform hemispherical projection */
95	<	double  norm = 1./sqrt(1. + vec[2]);
72	>	/* processed incident DSF measurements */
73	>	static RBFNODE          *dsf_list = NULL;
74
75	<	SDdisk2square(sq, vec[0]*norm, vec[1]*norm);
75	>	/* RBF-linking matrices (edges) */
76	>	static MIGRATION        *mig_list = NULL;
77
78	<	pos[0] = (int)(sq[0]*GRIDRES);
79	<	pos[1] = (int)(sq[1]*GRIDRES);
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	>	/* percentage to cull (<0 to turn off) */
93	>	int                     pctcull = 90;
94	>	/* sampling order */
95	>	int                     samp_order = 0;
96	>
97	>	/* Compute volume associated with Gaussian lobe */
98	>	static double
99	>	rbf_volume(const RBFVAL *rbfp)
100	>	{
101	>	double  rad = R2ANG(rbfp->crad);
102	>
103	>	return((2.*M_PI) * rbfp->peak * rad*rad);
104		}
105
106		/* Compute outgoing vector from grid position */
107		static void
108	<	vec_from_pos(FVECT vec, int xpos, int ypos)
108	>	ovec_from_pos(FVECT vec, int xpos, int ypos)
109		{
110		double  uv[2];
111		double  r2;
#	Line 113 | Line 116 | vec_from_pos(FVECT vec, int xpos, int ypos)
116		vec[0] = vec[1] = sqrt(2. - r2);
117		vec[0] *= uv[0];
118		vec[1] *= uv[1];
119	<	vec[2] = 1. - r2;
119	>	vec[2] = output_orient*(1. - r2);
120		}
121
122	<	/* Evaluate RBF at this grid position */
123	<	static double
124	<	eval_rbfrep2(const RBFLIST *rp, int xi, int yi)
122	>	/* Compute grid position from normalized input/output vector */
123	>	static void
124	>	pos_from_vec(int pos[2], const FVECT vec)
125		{
126	<	double          res = .0;
127	<	const RBFVAL    *rbfp;
125	<	double          sig2;
126	<	int             x2, y2;
127	<	int             n;
126	>	double  sq[2];          /* uniform hemispherical projection */
127	>	double  norm = 1./sqrt(1. + fabs(vec[2]));
128
129	<	rbfp = rp->rbfa;
130	<	for (n = rp->nrbf; n--; rbfp++) {
131	<	x2 = (signed)rbfp->gx - xi;
132	<	x2 *= x2;
133	<	y2 = (signed)rbfp->gy - yi;
134	<	y2 *= y2;
135	<	sig2 = -.5*(x2 + y2)/(double)(rbfp->rad*rbfp->rad);
136	<	if (sig2 > -19.)
137	<	res += rbfp->bsdf * exp(sig2);
138	<	}
139	<	return(res);
129	>	SDdisk2square(sq, vec[0]*norm, vec[1]*norm);
130	>
131	>	pos[0] = (int)(sq[0]*GRIDRES);
132	>	pos[1] = (int)(sq[1]*GRIDRES);
133		}
134
135	<	/* Evaluate RBF for BSDF at the given normalized outgoing direction */
135	>	/* Evaluate RBF for DSF at the given normalized outgoing direction */
136		static double
137	<	eval_rbfrep(const RBFLIST *rp, const FVECT outvec)
137	>	eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
138		{
139		double          res = .0;
140		const RBFVAL    *rbfp;
#	Line 149 | Line 142 | eval_rbfrep(const RBFLIST *rp, const FVECT outvec)
142		double          sig2;
143		int             n;
144
145	+	if (rp == NULL)
146	+	return(.0);
147		rbfp = rp->rbfa;
148		for (n = rp->nrbf; n--; rbfp++) {
149	<	vec_from_pos(odir, rbfp->gx, rbfp->gy);
150	<	sig2 = (DOT(odir, outvec) - 1.) /
151	<	((M_PI*M_PI/(double)(GRIDRES*GRIDRES)) *
157	<	rbfp->rad*rbfp->rad);
149	>	ovec_from_pos(odir, rbfp->gx, rbfp->gy);
150	>	sig2 = R2ANG(rbfp->crad);
151	>	sig2 = (DOT(odir,outvec) - 1.) / (sig2*sig2);
152		if (sig2 > -19.)
153	<	res += rbfp->bsdf * exp(sig2);
153	>	res += rbfp->peak * exp(sig2);
154		}
155		return(res);
156		}
157
158	+	/* Insert a new directional scattering function in our global list */
159	+	static void
160	+	insert_dsf(RBFNODE *newrbf)
161	+	{
162	+	RBFNODE         *rbf, *rbf_last;
163	+
164	+	/* keep in ascending theta order */
165	+	for (rbf_last = NULL, rbf = dsf_list;
166	+	single_plane_incident & (rbf != NULL);
167	+	rbf_last = rbf, rbf = rbf->next)
168	+	if (input_orient*rbf->invec[2] < input_orient*newrbf->invec[2])
169	+	break;
170	+	if (rbf_last == NULL) {
171	+	newrbf->next = dsf_list;
172	+	dsf_list = newrbf;
173	+	return;
174	+	}
175	+	newrbf->next = rbf;
176	+	rbf_last->next = newrbf;
177	+	}
178	+
179	+	/* Count up filled nodes and build RBF representation from current grid */
180	+	static RBFNODE *
181	+	make_rbfrep(void)
182	+	{
183	+	int     niter = 16;
184	+	int     minrad = ANG2R(pow(2., 1.-samp_order));
185	+	double  lastVar, thisVar = 100.;
186	+	int     nn;
187	+	RBFNODE *newnode;
188	+	int     i, j;
189	+
190	+	nn = 0;                 /* count selected bins */
191	+	for (i = 0; i < GRIDRES; i++)
192	+	for (j = 0; j < GRIDRES; j++)
193	+	nn += dsf_grid[i][j].nval;
194	+	/* allocate RBF array */
195	+	newnode = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(nn-1));
196	+	if (newnode == NULL) {
197	+	fputs("Out of memory in make_rbfrep()\n", stderr);
198	+	exit(1);
199	+	}
200	+	newnode->next = NULL;
201	+	newnode->ejl = NULL;
202	+	newnode->invec[2] = sin(M_PI/180.*theta_in_deg);
203	+	newnode->invec[0] = cos(M_PI/180.*phi_in_deg)*newnode->invec[2];
204	+	newnode->invec[1] = sin(M_PI/180.*phi_in_deg)*newnode->invec[2];
205	+	newnode->invec[2] = input_orient*sqrt(1. - newnode->invec[2]*newnode->invec[2]);
206	+	newnode->vtotal = 0;
207	+	newnode->nrbf = nn;
208	+	nn = 0;                 /* fill RBF array */
209	+	for (i = 0; i < GRIDRES; i++)
210	+	for (j = 0; j < GRIDRES; j++)
211	+	if (dsf_grid[i][j].nval) {
212	+	newnode->rbfa[nn].peak = dsf_grid[i][j].vsum;
213	+	newnode->rbfa[nn].crad = RSCA*dsf_grid[i][j].crad + .5;
214	+	newnode->rbfa[nn].gx = i;
215	+	newnode->rbfa[nn].gy = j;
216	+	if (newnode->rbfa[nn].crad < minrad)
217	+	minrad = newnode->rbfa[nn].crad;
218	+	++nn;
219	+	}
220	+	/* iterate to improve interpolation accuracy */
221	+	do {
222	+	double  dsum = .0, dsum2 = .0;
223	+	nn = 0;
224	+	for (i = 0; i < GRIDRES; i++)
225	+	for (j = 0; j < GRIDRES; j++)
226	+	if (dsf_grid[i][j].nval) {
227	+	FVECT   odir;
228	+	double  corr;
229	+	ovec_from_pos(odir, i, j);
230	+	newnode->rbfa[nn++].peak *= corr =
231	+	dsf_grid[i][j].vsum /
232	+	eval_rbfrep(newnode, odir);
233	+	dsum += corr - 1.;
234	+	dsum2 += (corr-1.)*(corr-1.);
235	+	}
236	+	lastVar = thisVar;
237	+	thisVar = dsum2/(double)nn;
238	+	#ifdef DEBUG
239	+	fprintf(stderr, "Avg., RMS error: %.1f%%  %.1f%%\n",
240	+	100.*dsum/(double)nn,
241	+	100.*sqrt(thisVar));
242	+	#endif
243	+	} while (--niter > 0 && lastVar-thisVar > 0.02*lastVar);
244	+
245	+	nn = 0;                 /* compute sum for normalization */
246	+	while (nn < newnode->nrbf)
247	+	newnode->vtotal += rbf_volume(&newnode->rbfa[nn++]);
248	+
249	+	insert_dsf(newnode);
250	+	/* adjust sampling resolution */
251	+	samp_order = log(2./R2ANG(minrad))/log(2.) + .5;
252	+
253	+	return(newnode);
254	+	}
255	+
256		/* Load a set of measurements corresponding to a particular incident angle */
257		static int
258	<	load_bsdf_meas(const char *fname)
258	>	load_pabopto_meas(const char *fname)
259		{
260		FILE    *fp = fopen(fname, "r");
261		int     inp_is_DSF = -1;
262	<	double  theta_out, phi_out, val;
262	>	double  new_phi, theta_out, phi_out, val;
263		char    buf[2048];
264		int     n, c;
265
#	Line 176 | Line 268 | load_bsdf_meas(const char *fname)
268		fputs(": cannot open\n", stderr);
269		return(0);
270		}
271	<	memset(bsdf_grid, 0, sizeof(bsdf_grid));
271	>	memset(dsf_grid, 0, sizeof(dsf_grid));
272	>	#ifdef DEBUG
273	>	fprintf(stderr, "Loading measurement file '%s'...\n", fname);
274	>	#endif
275		/* read header information */
276		while ((c = getc(fp)) == '#' || c == EOF) {
277		if (fgets(buf, sizeof(buf), fp) == NULL) {
#	Line 195 | Line 290 | load_bsdf_meas(const char *fname)
290		}
291		if (sscanf(buf, "intheta %lf", &theta_in_deg) == 1)
292		continue;
293	<	if (sscanf(buf, "inphi %lf", &phi_in_deg) == 1)
293	>	if (sscanf(buf, "inphi %lf", &new_phi) == 1)
294		continue;
295		if (sscanf(buf, "incident_angle %lf %lf",
296	<	&theta_in_deg, &phi_in_deg) == 2)
296	>	&theta_in_deg, &new_phi) == 2)
297		continue;
298		}
299		if (inp_is_DSF < 0) {
#	Line 207 | Line 302 | load_bsdf_meas(const char *fname)
302		fclose(fp);
303		return(0);
304		}
305	<	ungetc(c, fp);          /* read actual data */
305	>	if (!input_orient)              /* check input orientation */
306	>	input_orient = 1 - 2*(theta_in_deg > 90.);
307	>	else if (input_orient > 0 ^ theta_in_deg < 90.) {
308	>	fputs("Cannot handle input angles on both sides of surface\n",
309	>	stderr);
310	>	exit(1);
311	>	}
312	>	if (single_plane_incident > 0)  /* check if still in plane */
313	>	single_plane_incident = (round(new_phi) == round(phi_in_deg));
314	>	else if (single_plane_incident < 0)
315	>	single_plane_incident = 1;
316	>	phi_in_deg = new_phi;
317	>	ungetc(c, fp);                  /* read actual data */
318		while (fscanf(fp, "%lf %lf %lf\n", &theta_out, &phi_out, &val) == 3) {
319		FVECT   ovec;
320		int     pos[2];
321
322	+	if (!output_orient)     /* check output orientation */
323	+	output_orient = 1 - 2*(theta_out > 90.);
324	+	else if (output_orient > 0 ^ theta_out < 90.) {
325	+	fputs("Cannot handle output angles on both sides of surface\n",
326	+	stderr);
327	+	exit(1);
328	+	}
329		ovec[2] = sin(M_PI/180.*theta_out);
330		ovec[0] = cos(M_PI/180.*phi_out) * ovec[2];
331		ovec[1] = sin(M_PI/180.*phi_out) * ovec[2];
332		ovec[2] = sqrt(1. - ovec[2]*ovec[2]);
333
334	<	if (inp_is_DSF)
335	<	val /= ovec[2]; /* convert from DSF to BSDF */
334	>	if (!inp_is_DSF)
335	>	val *= ovec[2]; /* convert from BSDF to DSF */
336
337		pos_from_vec(pos, ovec);
338
339	<	bsdf_grid[pos[0]][pos[1]].vsum += val;
340	<	bsdf_grid[pos[0]][pos[1]].nval++;
339	>	dsf_grid[pos[0]][pos[1]].vsum += val;
340	>	dsf_grid[pos[0]][pos[1]].nval++;
341		}
342		n = 0;
343		while ((c = getc(fp)) != EOF)
#	Line 241 | Line 355 | load_bsdf_meas(const char *fname)
355		static void
356		compute_radii(void)
357		{
358	<	unsigned char   fill_grid[GRIDRES][GRIDRES];
359	<	int             r, r2, lastr2;
360	<	int             i, j, jn, ii, jj, inear, jnear;
361	<	/* proceed in zig-zag */
362	<	lastr2 = GRIDRES*GRIDRES;
358	>	unsigned int    fill_grid[GRIDRES][GRIDRES];
359	>	unsigned short  fill_cnt[GRIDRES][GRIDRES];
360	>	FVECT           ovec0, ovec1;
361	>	double          ang2, lastang2;
362	>	int             r, i, j, jn, ii, jj, inear, jnear;
363	>
364	>	r = GRIDRES/2;                          /* proceed in zig-zag */
365		for (i = 0; i < GRIDRES; i++)
366		for (jn = 0; jn < GRIDRES; jn++) {
367		j = (i&1) ? jn : GRIDRES-1-jn;
368	<	if (bsdf_grid[i][j].nval)       /* find empty grid pos. */
368	>	if (dsf_grid[i][j].nval)        /* find empty grid pos. */
369		continue;
370	<	r = (int)sqrt((double)lastr2) + 2;
370	>	ovec_from_pos(ovec0, i, j);
371		inear = jnear = -1;             /* find nearest non-empty */
372	<	lastr2 = 2*GRIDRES*GRIDRES;
372	>	lastang2 = M_PI*M_PI;
373		for (ii = i-r; ii <= i+r; ii++) {
374		if (ii < 0) continue;
375		if (ii >= GRIDRES) break;
376		for (jj = j-r; jj <= j+r; jj++) {
377		if (jj < 0) continue;
378		if (jj >= GRIDRES) break;
379	<	if (!bsdf_grid[ii][jj].nval)
379	>	if (!dsf_grid[ii][jj].nval)
380		continue;
381	<	r2 = (ii-i)*(ii-i) + (jj-j)*(jj-j);
382	<	if (r2 >= lastr2)
381	>	ovec_from_pos(ovec1, ii, jj);
382	>	ang2 = 2. - 2.*DOT(ovec0,ovec1);
383	>	if (ang2 >= lastang2)
384		continue;
385	<	lastr2 = r2;
385	>	lastang2 = ang2;
386		inear = ii; jnear = jj;
387		}
388		}
389	<	/* record if > previous */
390	<	if (bsdf_grid[inear][jnear].hrad2 < lastr2)
391	<	bsdf_grid[inear][jnear].hrad2 = lastr2;
389	>	if (inear < 0) {
390	>	fputs("Could not find non-empty neighbor!\n", stderr);
391	>	exit(1);
392	>	}
393	>	ang2 = sqrt(lastang2);
394	>	r = ANG2R(ang2);                /* record if > previous */
395	>	if (r > dsf_grid[inear][jnear].crad)
396	>	dsf_grid[inear][jnear].crad = r;
397	>	/* next search radius */
398	>	r = ang2*(2.*GRIDRES/M_PI) + 3;
399		}
400	<	/* fill in others */
400	>	/* blur radii over hemisphere */
401		memset(fill_grid, 0, sizeof(fill_grid));
402	+	memset(fill_cnt, 0, sizeof(fill_cnt));
403		for (i = 0; i < GRIDRES; i++)
404		for (j = 0; j < GRIDRES; j++) {
405	<	if (!bsdf_grid[i][j].nval)
406	<	continue;
407	<	if (bsdf_grid[i][j].hrad2)
283	<	continue;
284	<	r = GRIDRES/20;
285	<	lastr2 = 2*r*r;
405	>	if (!dsf_grid[i][j].crad)
406	>	continue;               /* missing distance */
407	>	r = R2ANG(dsf_grid[i][j].crad)*(2.*RSCA*GRIDRES/M_PI);
408		for (ii = i-r; ii <= i+r; ii++) {
409		if (ii < 0) continue;
410		if (ii >= GRIDRES) break;
411		for (jj = j-r; jj <= j+r; jj++) {
412		if (jj < 0) continue;
413		if (jj >= GRIDRES) break;
414	<	if (!bsdf_grid[ii][jj].hrad2)
414	>	if ((ii-i)*(ii-i) + (jj-j)*(jj-j) > r*r)
415		continue;
416	<	r2 = (ii-i)*(ii-i) + (jj-j)*(jj-j);
417	<	if (r2 >= lastr2)
296	<	continue;
297	<	fill_grid[i][j] = bsdf_grid[ii][jj].hrad2;
298	<	lastr2 = r2;
416	>	fill_grid[ii][jj] += dsf_grid[i][j].crad;
417	>	fill_cnt[ii][jj]++;
418		}
419		}
420		}
421	+	/* copy back blurred radii */
422		for (i = 0; i < GRIDRES; i++)
423		for (j = 0; j < GRIDRES; j++)
424	<	if (fill_grid[i][j])
425	<	bsdf_grid[i][j].hrad2 = fill_grid[i][j];
424	>	if (fill_cnt[i][j])
425	>	dsf_grid[i][j].crad = fill_grid[i][j]/fill_cnt[i][j];
426		}
427
428	<	/* Cull points for more uniform distribution */
428	>	/* Cull points for more uniform distribution, leave all nval 0 or 1 */
429		static void
430		cull_values(void)
431		{
432	<	int     i, j, ii, jj, r, r2;
432	>	FVECT   ovec0, ovec1;
433	>	double  maxang, maxang2;
434	>	int     i, j, ii, jj, r;
435		/* simple greedy algorithm */
436		for (i = 0; i < GRIDRES; i++)
437		for (j = 0; j < GRIDRES; j++) {
438	<	if (!bsdf_grid[i][j].nval)
438	>	if (!dsf_grid[i][j].nval)
439		continue;
440	<	if (!bsdf_grid[i][j].hrad2)
441	<	continue;
442	<	r = (int)(2.*sqrt((double)bsdf_grid[i][j].hrad2) + .9999);
440	>	if (!dsf_grid[i][j].crad)
441	>	continue;               /* shouldn't happen */
442	>	ovec_from_pos(ovec0, i, j);
443	>	maxang = 2.*R2ANG(dsf_grid[i][j].crad);
444	>	if (maxang > ovec0[2])          /* clamp near horizon */
445	>	maxang = ovec0[2];
446	>	r = maxang*(2.*GRIDRES/M_PI) + 1;
447	>	maxang2 = maxang*maxang;
448		for (ii = i-r; ii <= i+r; ii++) {
449		if (ii < 0) continue;
450		if (ii >= GRIDRES) break;
451		for (jj = j-r; jj <= j+r; jj++) {
452		if (jj < 0) continue;
453		if (jj >= GRIDRES) break;
454	<	if (!bsdf_grid[ii][jj].nval)
454	>	if (!dsf_grid[ii][jj].nval)
455		continue;
456	<	r2 = (ii-i)*(ii-i) + (jj-j)*(jj-j);
457	<	if (!r2 | (r2 > r*r))
456	>	if ((ii == i) & (jj == j))
457	>	continue;       /* don't get self-absorbed */
458	>	ovec_from_pos(ovec1, ii, jj);
459	>	if (2. - 2.*DOT(ovec0,ovec1) >= maxang2)
460		continue;
461	<	/* absorb victim's value */
462	<	bsdf_grid[i][j].vsum += bsdf_grid[ii][jj].vsum;
463	<	bsdf_grid[i][j].nval += bsdf_grid[ii][jj].nval;
464	<	memset(&bsdf_grid[ii][jj], 0, sizeof(GRIDVAL));
461	>	/* absorb sum */
462	>	dsf_grid[i][j].vsum += dsf_grid[ii][jj].vsum;
463	>	dsf_grid[i][j].nval += dsf_grid[ii][jj].nval;
464	>	/* keep value, though */
465	>	dsf_grid[ii][jj].vsum /= (float)dsf_grid[ii][jj].nval;
466	>	dsf_grid[ii][jj].nval = 0;
467		}
468		}
469		}
470	+	/* final averaging pass */
471	+	for (i = 0; i < GRIDRES; i++)
472	+	for (j = 0; j < GRIDRES; j++)
473	+	if (dsf_grid[i][j].nval > 1) {
474	+	dsf_grid[i][j].vsum /= (float)dsf_grid[i][j].nval;
475	+	dsf_grid[i][j].nval = 1;
476	+	}
477		}
478
479	+	/* Compute (and allocate) migration price matrix for optimization */
480	+	static float *
481	+	price_routes(const RBFNODE *from_rbf, const RBFNODE *to_rbf)
482	+	{
483	+	float   *pmtx = (float *)malloc(sizeof(float) *
484	+	from_rbf->nrbf * to_rbf->nrbf);
485	+	FVECT   *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
486	+	int     i, j;
487
488	+	if ((pmtx == NULL) | (vto == NULL)) {
489	+	fputs("Out of memory in migration_costs()\n", stderr);
490	+	exit(1);
491	+	}
492	+	for (j = to_rbf->nrbf; j--; )           /* save repetitive ops. */
493	+	ovec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);
494	+
495	+	for (i = from_rbf->nrbf; i--; ) {
496	+	const double        from_ang = R2ANG(from_rbf->rbfa[i].crad);
497	+	FVECT               vfrom;
498	+	ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
499	+	for (j = to_rbf->nrbf; j--; )
500	+	pmtx[i*to_rbf->nrbf + j] = acos(DOT(vfrom, vto[j])) +
501	+	fabs(R2ANG(to_rbf->rbfa[j].crad) - from_ang);
502	+	}
503	+	free(vto);
504	+	return(pmtx);
505	+	}
506	+
507	+	/* Comparison routine needed for sorting price row */
508	+	static const float      *price_arr;
509	+	static int
510	+	msrt_cmp(const void *p1, const void *p2)
511	+	{
512	+	float   c1 = price_arr[*(const int *)p1];
513	+	float   c2 = price_arr[*(const int *)p2];
514	+
515	+	if (c1 > c2) return(1);
516	+	if (c1 < c2) return(-1);
517	+	return(0);
518	+	}
519	+
520	+	/* Compute minimum (optimistic) cost for moving the given source material */
521	+	static double
522	+	min_cost(double amt2move, const double *avail, const float *price, int n)
523	+	{
524	+	static int      *price_sort = NULL;
525	+	static int      n_alloc = 0;
526	+	double          total_cost = 0;
527	+	int             i;
528	+
529	+	if (amt2move <= FTINY)                  /* pre-emptive check */
530	+	return(0.);
531	+	if (n > n_alloc) {                      /* (re)allocate sort array */
532	+	if (n_alloc) free(price_sort);
533	+	price_sort = (int *)malloc(sizeof(int)*n);
534	+	if (price_sort == NULL) {
535	+	fputs("Out of memory in min_cost()\n", stderr);
536	+	exit(1);
537	+	}
538	+	n_alloc = n;
539	+	}
540	+	for (i = n; i--; )
541	+	price_sort[i] = i;
542	+	price_arr = price;
543	+	qsort(price_sort, n, sizeof(int), &msrt_cmp);
544	+	/* move cheapest first */
545	+	for (i = 0; i < n && amt2move > FTINY; i++) {
546	+	int     d = price_sort[i];
547	+	double  amt = (amt2move < avail[d]) ? amt2move : avail[d];
548	+
549	+	total_cost += amt * price[d];
550	+	amt2move -= amt;
551	+	}
552	+	return(total_cost);
553	+	}
554	+
555	+	/* Take a step in migration by choosing optimal bucket to transfer */
556	+	static double
557	+	migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, const float *pmtx)
558	+	{
559	+	static double   *src_cost = NULL;
560	+	int             n_alloc = 0;
561	+	const double    maxamt = 2./(mtx_nrows(mig)*mtx_ncols(mig));
562	+	double          amt = 0;
563	+	struct {
564	+	int     s, d;   /* source and destination */
565	+	double  price;  /* price estimate per amount moved */
566	+	double  amt;    /* amount we can move */
567	+	} cur, best;
568	+	int             i;
569	+
570	+	if (mtx_nrows(mig) > n_alloc) {         /* allocate cost array */
571	+	if (n_alloc)
572	+	free(src_cost);
573	+	src_cost = (double *)malloc(sizeof(double)*mtx_nrows(mig));
574	+	if (src_cost == NULL) {
575	+	fputs("Out of memory in migration_step()\n", stderr);
576	+	exit(1);
577	+	}
578	+	n_alloc = mtx_nrows(mig);
579	+	}
580	+	for (i = mtx_nrows(mig); i--; )         /* starting costs for diff. */
581	+	src_cost[i] = min_cost(src_rem[i], dst_rem,
582	+	pmtx+i*mtx_ncols(mig), mtx_ncols(mig));
583	+
584	+	/* find best source & dest. */
585	+	best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
586	+	for (cur.s = mtx_nrows(mig); cur.s--; ) {
587	+	const float *price = pmtx + cur.s*mtx_ncols(mig);
588	+	double      cost_others = 0;
589	+	if (src_rem[cur.s] <= FTINY)
590	+	continue;
591	+	cur.d = -1;                         /* examine cheapest dest. */
592	+	for (i = mtx_ncols(mig); i--; )
593	+	if (dst_rem[i] > FTINY &&
594	+	(cur.d < 0 || price[i] < price[cur.d]))
595	+	cur.d = i;
596	+	if (cur.d < 0)
597	+	return(.0);
598	+	if ((cur.price = price[cur.d]) >= best.price)
599	+	continue;                   /* no point checking further */
600	+	cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
601	+	src_rem[cur.s] : dst_rem[cur.d];
602	+	if (cur.amt > maxamt) cur.amt = maxamt;
603	+	dst_rem[cur.d] -= cur.amt;          /* add up differential costs */
604	+	for (i = mtx_nrows(mig); i--; ) {
605	+	if (i == cur.s) continue;
606	+	cost_others += min_cost(src_rem[i], dst_rem, price, mtx_ncols(mig))
607	+	- src_cost[i];
608	+	}
609	+	dst_rem[cur.d] += cur.amt;          /* undo trial move */
610	+	cur.price += cost_others/cur.amt;   /* adjust effective price */
611	+	if (cur.price < best.price)         /* are we better than best? */
612	+	best = cur;
613	+	}
614	+	if ((best.s < 0) | (best.d < 0))
615	+	return(.0);
616	+	/* make the actual move */
617	+	mig->mtx[mtx_ndx(mig,best.s,best.d)] += best.amt;
618	+	src_rem[best.s] -= best.amt;
619	+	dst_rem[best.d] -= best.amt;
620	+	return(best.amt);
621	+	}
622	+
623	+	/* Compute (and insert) migration along directed edge */
624	+	static MIGRATION *
625	+	make_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
626	+	{
627	+	const double    end_thresh = 0.02/(from_rbf->nrbf*to_rbf->nrbf);
628	+	float           *pmtx = price_routes(from_rbf, to_rbf);
629	+	MIGRATION       *newmig = (MIGRATION *)malloc(sizeof(MIGRATION) +
630	+	sizeof(float) *
631	+	(from_rbf->nrbf*to_rbf->nrbf - 1));
632	+	double          *src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
633	+	double          *dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
634	+	double          total_rem = 1.;
635	+	int             i;
636	+
637	+	if ((newmig == NULL) | (src_rem == NULL) | (dst_rem == NULL)) {
638	+	fputs("Out of memory in make_migration()\n", stderr);
639	+	exit(1);
640	+	}
641	+	#ifdef DEBUG
642	+	{
643	+	double  theta, phi;
644	+	theta = acos(from_rbf->invec[2])*(180./M_PI);
645	+	phi = atan2(from_rbf->invec[1],from_rbf->invec[0])*(180./M_PI);
646	+	fprintf(stderr, "Building path from (theta,phi) (%d,%d) to ",
647	+	round(theta), round(phi));
648	+	theta = acos(to_rbf->invec[2])*(180./M_PI);
649	+	phi = atan2(to_rbf->invec[1],to_rbf->invec[0])*(180./M_PI);
650	+	fprintf(stderr, "(%d,%d)\n", round(theta), round(phi));
651	+	}
652	+	#endif
653	+	newmig->next = NULL;
654	+	newmig->rbfv[0] = from_rbf;
655	+	newmig->rbfv[1] = to_rbf;
656	+	newmig->enxt[0] = newmig->enxt[1] = NULL;
657	+	memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
658	+	/* starting quantities */
659	+	for (i = from_rbf->nrbf; i--; )
660	+	src_rem[i] = rbf_volume(&from_rbf->rbfa[i]) / from_rbf->vtotal;
661	+	for (i = to_rbf->nrbf; i--; )
662	+	dst_rem[i] = rbf_volume(&to_rbf->rbfa[i]) / to_rbf->vtotal;
663	+	/* move a bit at a time */
664	+	while (total_rem > end_thresh)
665	+	total_rem -= migration_step(newmig, src_rem, dst_rem, pmtx);
666	+
667	+	free(pmtx);                             /* free working arrays */
668	+	free(src_rem);
669	+	free(dst_rem);
670	+	for (i = from_rbf->nrbf; i--; ) {       /* normalize final matrix */
671	+	float       nf = rbf_volume(&from_rbf->rbfa[i]);
672	+	int         j;
673	+	if (nf <= FTINY) continue;
674	+	nf = from_rbf->vtotal / nf;
675	+	for (j = to_rbf->nrbf; j--; )
676	+	newmig->mtx[mtx_ndx(newmig,i,j)] *= nf;
677	+	}
678	+	/* insert in edge lists */
679	+	newmig->enxt[0] = from_rbf->ejl;
680	+	from_rbf->ejl = newmig;
681	+	newmig->enxt[1] = to_rbf->ejl;
682	+	to_rbf->ejl = newmig;
683	+	newmig->next = mig_list;
684	+	return(mig_list = newmig);
685	+	}
686	+
687	+	/* Get triangle surface orientation (unnormalized) */
688	+	static void
689	+	tri_orient(FVECT vres, const FVECT v1, const FVECT v2, const FVECT v3)
690	+	{
691	+	FVECT   v2minus1, v3minus2;
692	+
693	+	VSUB(v2minus1, v2, v1);
694	+	VSUB(v3minus2, v3, v2);
695	+	VCROSS(vres, v2minus1, v3minus2);
696	+	}
697	+
698	+	/* Determine if vertex order is reversed (inward normal) */
699	+	static int
700	+	is_rev_tri(const FVECT v1, const FVECT v2, const FVECT v3)
701	+	{
702	+	FVECT   tor;
703	+
704	+	tri_orient(tor, v1, v2, v3);
705	+
706	+	return(DOT(tor, v2) < 0.);
707	+	}
708	+
709	+	/* Find vertices completing triangles on either side of the given edge */
710	+	static int
711	+	get_triangles(RBFNODE *rbfv[2], const MIGRATION *mig)
712	+	{
713	+	const MIGRATION *ej, *ej2;
714	+	RBFNODE         *tv;
715	+
716	+	rbfv[0] = rbfv[1] = NULL;
717	+	for (ej = mig->rbfv[0]->ejl; ej != NULL;
718	+	ej = nextedge(mig->rbfv[0],ej)) {
719	+	if (ej == mig)
720	+	continue;
721	+	tv = opp_rbf(mig->rbfv[0],ej);
722	+	for (ej2 = tv->ejl; ej2 != NULL; ej2 = nextedge(tv,ej2))
723	+	if (opp_rbf(tv,ej2) == mig->rbfv[1]) {
724	+	rbfv[is_rev_tri(mig->rbfv[0]->invec,
725	+	mig->rbfv[1]->invec,
726	+	tv->invec)] = tv;
727	+	break;
728	+	}
729	+	}
730	+	return((rbfv[0] != NULL) + (rbfv[1] != NULL));
731	+	}
732	+
733	+	/* Find context hull vertex to complete triangle (oriented call) */
734	+	static RBFNODE *
735	+	find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rbf1)
736	+	{
737	+	FVECT   vmid, vor;
738	+	RBFNODE *rbf, *rbfbest = NULL;
739	+	double  dprod2, bestdprod2 = 0.5;
740	+
741	+	VADD(vmid, rbf0->invec, rbf1->invec);
742	+	if (normalize(vmid) == 0)
743	+	return(NULL);
744	+	/* XXX exhaustive search */
745	+	for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
746	+	if ((rbf == rbf0) | (rbf == rbf1))
747	+	continue;
748	+	tri_orient(vor, rbf0->invec, rbf1->invec, rbf->invec);
749	+	dprod2 = DOT(vor, vmid);
750	+	if (dprod2 <= FTINY)
751	+	continue;               /* wrong orientation */
752	+	dprod2 *= dprod2 / DOT(vor,vor);
753	+	if (dprod2 > bestdprod2) {      /* more convex than prev? */
754	+	rbfbest = rbf;
755	+	bestdprod2 = dprod2;
756	+	}
757	+	}
758	+	return(rbf);
759	+	}
760	+
761	+	/* Create new migration edge and grow mesh recursively around it */
762	+	static void
763	+	mesh_from_edge(RBFNODE *rbf0, RBFNODE *rbf1)
764	+	{
765	+	MIGRATION       *newej;
766	+	RBFNODE         *tvert[2];
767	+
768	+	if (rbf0 < rbf1)                        /* avoid migration loops */
769	+	newej = make_migration(rbf0, rbf1);
770	+	else
771	+	newej = make_migration(rbf1, rbf0);
772	+	/* triangle on either side? */
773	+	get_triangles(tvert, newej);
774	+	if (tvert[0] == NULL) {                 /* recurse on new right edge */
775	+	tvert[0] = find_chull_vert(newej->rbfv[0], newej->rbfv[1]);
776	+	if (tvert[0] != NULL) {
777	+	mesh_from_edge(rbf0, tvert[0]);
778	+	mesh_from_edge(rbf1, tvert[0]);
779	+	}
780	+	}
781	+	if (tvert[1] == NULL) {                 /* recurse on new left edge */
782	+	tvert[1] = find_chull_vert(newej->rbfv[1], newej->rbfv[0]);
783	+	if (tvert[1] != NULL) {
784	+	mesh_from_edge(rbf0, tvert[1]);
785	+	mesh_from_edge(rbf1, tvert[1]);
786	+	}
787	+	}
788	+	}
789	+
790	+	/* Draw edge list into mig_grid array */
791	+	static void
792	+	draw_edges()
793	+	{
794	+	int             nnull = 0, ntot = 0;
795	+	MIGRATION       *ej;
796	+	int             p0[2], p1[2];
797	+
798	+	/* memset(mig_grid, 0, sizeof(mig_grid)); */
799	+	for (ej = mig_list; ej != NULL; ej = ej->next) {
800	+	++ntot;
801	+	pos_from_vec(p0, ej->rbfv[0]->invec);
802	+	pos_from_vec(p1, ej->rbfv[1]->invec);
803	+	if ((p0[0] == p1[0]) & (p0[1] == p1[1])) {
804	+	++nnull;
805	+	mig_grid[p0[0]][p0[1]] = ej;
806	+	continue;
807	+	}
808	+	if (abs(p1[0]-p0[0]) > abs(p1[1]-p0[1])) {
809	+	const int       xstep = 2*(p1[0] > p0[0]) - 1;
810	+	const double    ystep = (double)((p1[1]-p0[1])*xstep) /
811	+	(double)(p1[0]-p0[0]);
812	+	int             x;
813	+	double          y;
814	+	for (x = p0[0], y = p0[1]+.5; x != p1[0];
815	+	x += xstep, y += ystep)
816	+	mig_grid[x][(int)y] = ej;
817	+	mig_grid[x][(int)y] = ej;
818	+	} else {
819	+	const int       ystep = 2*(p1[1] > p0[1]) - 1;
820	+	const double    xstep = (double)((p1[0]-p0[0])*ystep) /
821	+	(double)(p1[1]-p0[1]);
822	+	int             y;
823	+	double          x;
824	+	for (y = p0[1], x = p0[0]+.5; y != p1[1];
825	+	y += ystep, x += xstep)
826	+	mig_grid[(int)x][y] = ej;
827	+	mig_grid[(int)x][y] = ej;
828	+	}
829	+	}
830	+	if (nnull)
831	+	fprintf(stderr, "Warning: %d of %d edges are null\n",
832	+	nnull, ntot);
833	+	}
834	+
835	+	/* Build our triangle mesh from recorded RBFs */
836	+	static void
837	+	build_mesh()
838	+	{
839	+	double          best2 = M_PI*M_PI;
840	+	RBFNODE         *rbf, *rbf_near = NULL;
841	+	/* check if isotropic */
842	+	if (single_plane_incident) {
843	+	for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
844	+	if (rbf->next != NULL)
845	+	make_migration(rbf, rbf->next);
846	+	return;
847	+	}
848	+	/* find RBF nearest to head */
849	+	if (dsf_list == NULL)
850	+	return;
851	+	for (rbf = dsf_list->next; rbf != NULL; rbf = rbf->next) {
852	+	double  dist2 = 2. - 2.*DOT(dsf_list->invec,rbf->invec);
853	+	if (dist2 < best2) {
854	+	rbf_near = rbf;
855	+	best2 = dist2;
856	+	}
857	+	}
858	+	if (rbf_near == NULL) {
859	+	fputs("Cannot find nearest point for first edge\n", stderr);
860	+	exit(1);
861	+	}
862	+	/* build mesh from this edge */
863	+	mesh_from_edge(dsf_list, rbf_near);
864	+	/* draw edge list into grid */
865	+	draw_edges();
866	+	}
867	+
868	+	/* Identify enclosing triangle for this position (flood fill raster check) */
869	+	static int
870	+	identify_tri(MIGRATION *miga[3], unsigned char vmap[GRIDRES][(GRIDRES+7)/8],
871	+	int px, int py)
872	+	{
873	+	const int       btest = 1<<(py&07);
874	+
875	+	if (vmap[px][py>>3] & btest)            /* already visited here? */
876	+	return(1);
877	+	/* else mark it */
878	+	vmap[px][py>>3] |= btest;
879	+
880	+	if (mig_grid[px][py] != NULL) {         /* are we on an edge? */
881	+	int     i;
882	+	for (i = 0; i < 3; i++) {
883	+	if (miga[i] == mig_grid[px][py])
884	+	return(1);
885	+	if (miga[i] != NULL)
886	+	continue;
887	+	miga[i] = mig_grid[px][py];
888	+	return(1);
889	+	}
890	+	return(0);                      /* outside triangle! */
891	+	}
892	+	/* check neighbors (flood) */
893	+	if (px > 0 && !identify_tri(miga, vmap, px-1, py))
894	+	return(0);
895	+	if (px < GRIDRES-1 && !identify_tri(miga, vmap, px+1, py))
896	+	return(0);
897	+	if (py > 0 && !identify_tri(miga, vmap, px, py-1))
898	+	return(0);
899	+	if (py < GRIDRES-1 && !identify_tri(miga, vmap, px, py+1))
900	+	return(0);
901	+	return(1);                              /* this neighborhood done */
902	+	}
903	+
904	+	/* Find edge(s) for interpolating the given incident vector */
905	+	static int
906	+	get_interp(MIGRATION *miga[3], const FVECT invec)
907	+	{
908	+	miga[0] = miga[1] = miga[2] = NULL;
909	+	if (single_plane_incident) {            /* isotropic BSDF? */
910	+	RBFNODE *rbf;                   /* find edge we're on */
911	+	for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
912	+	if (input_orient*rbf->invec[2] < input_orient*invec[2])
913	+	break;
914	+	if (rbf->next != NULL &&
915	+	input_orient*rbf->next->invec[2] <
916	+	input_orient*invec[2]) {
917	+	for (miga[0] = rbf->ejl; miga[0] != NULL;
918	+	miga[0] = nextedge(rbf,miga[0]))
919	+	if (opp_rbf(rbf,miga[0]) == rbf->next)
920	+	return(1);
921	+	break;
922	+	}
923	+	}
924	+	return(0);                      /* outside range! */
925	+	}
926	+	{                                       /* else use triangle mesh */
927	+	unsigned char   floodmap[GRIDRES][(GRIDRES+7)/8];
928	+	int             pstart[2];
929	+
930	+	pos_from_vec(pstart, invec);
931	+	memset(floodmap, 0, sizeof(floodmap));
932	+	/* call flooding function */
933	+	if (!identify_tri(miga, floodmap, pstart[0], pstart[1]))
934	+	return(0);              /* outside mesh */
935	+	if ((miga[0] == NULL) | (miga[2] == NULL))
936	+	return(0);              /* should never happen */
937	+	if (miga[1] == NULL)
938	+	return(1);              /* on edge */
939	+	return(3);                      /* else in triangle */
940	+	}
941	+	}
942	+
943	+	/* Advect and allocate new RBF along edge */
944	+	static RBFNODE *
945	+	e_advect_rbf(const MIGRATION *mig, const FVECT invec)
946	+	{
947	+	RBFNODE         *rbf;
948	+	int             n, i, j;
949	+	double          t, full_dist;
950	+	/* get relative position */
951	+	t = acos(DOT(invec, mig->rbfv[0]->invec));
952	+	if (t < M_PI/GRIDRES) {                 /* near first DSF */
953	+	n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[0]->nrbf-1);
954	+	rbf = (RBFNODE *)malloc(n);
955	+	if (rbf == NULL)
956	+	goto memerr;
957	+	memcpy(rbf, mig->rbfv[0], n);   /* just duplicate */
958	+	return(rbf);
959	+	}
960	+	full_dist = acos(DOT(mig->rbfv[0]->invec, mig->rbfv[1]->invec));
961	+	if (t > full_dist-M_PI/GRIDRES) {       /* near second DSF */
962	+	n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[1]->nrbf-1);
963	+	rbf = (RBFNODE *)malloc(n);
964	+	if (rbf == NULL)
965	+	goto memerr;
966	+	memcpy(rbf, mig->rbfv[1], n);   /* just duplicate */
967	+	return(rbf);
968	+	}
969	+	t /= full_dist;
970	+	n = 0;                                  /* count migrating particles */
971	+	for (i = 0; i < mtx_nrows(mig); i++)
972	+	for (j = 0; j < mtx_ncols(mig); j++)
973	+	n += (mig->mtx[mtx_ndx(mig,i,j)] > FTINY);
974	+	#ifdef DEBUG
975	+	fprintf(stderr, "Input RBFs have %d, %d nodes -> output has %d\n",
976	+	mig->rbfv[0]->nrbf, mig->rbfv[1]->nrbf, n);
977	+	#endif
978	+	rbf = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(n-1));
979	+	if (rbf == NULL)
980	+	goto memerr;
981	+	rbf->next = NULL; rbf->ejl = NULL;
982	+	VCOPY(rbf->invec, invec);
983	+	rbf->nrbf = n;
984	+	rbf->vtotal = 1.-t + t*mig->rbfv[1]->vtotal/mig->rbfv[0]->vtotal;
985	+	n = 0;                                  /* advect RBF lobes */
986	+	for (i = 0; i < mtx_nrows(mig); i++) {
987	+	const RBFVAL        *rbf0i = &mig->rbfv[0]->rbfa[i];
988	+	const float         peak0 = rbf0i->peak;
989	+	const double        rad0 = R2ANG(rbf0i->crad);
990	+	FVECT               v0;
991	+	float               mv;
992	+	ovec_from_pos(v0, rbf0i->gx, rbf0i->gy);
993	+	for (j = 0; j < mtx_ncols(mig); j++)
994	+	if ((mv = mig->mtx[mtx_ndx(mig,i,j)]) > FTINY) {
995	+	const RBFVAL    *rbf1j = &mig->rbfv[1]->rbfa[j];
996	+	double          rad1 = R2ANG(rbf1j->crad);
997	+	FVECT           v;
998	+	int             pos[2];
999	+	rbf->rbfa[n].peak = peak0 * mv * rbf->vtotal;
1000	+	rbf->rbfa[n].crad = ANG2R(sqrt(rad0*rad0*(1.-t) +
1001	+	rad1*rad1*t));
1002	+	ovec_from_pos(v, rbf1j->gx, rbf1j->gy);
1003	+	geodesic(v, v0, v, t, GEOD_REL);
1004	+	pos_from_vec(pos, v);
1005	+	rbf->rbfa[n].gx = pos[0];
1006	+	rbf->rbfa[n].gy = pos[1];
1007	+	++n;
1008	+	}
1009	+	}
1010	+	rbf->vtotal *= mig->rbfv[0]->vtotal;    /* turn ratio into actual */
1011	+	return(rbf);
1012	+	memerr:
1013	+	fputs("Out of memory in e_advect_rbf()\n", stderr);
1014	+	exit(1);
1015	+	return(NULL);   /* pro forma return */
1016	+	}
1017	+
1018	+	/* Insert vertex in ordered list */
1019	+	static void
1020	+	insert_vert(RBFNODE **vlist, RBFNODE *v)
1021	+	{
1022	+	int     i, j;
1023	+
1024	+	for (i = 0; vlist[i] != NULL; i++) {
1025	+	if (v == vlist[i])
1026	+	return;
1027	+	if (v < vlist[i])
1028	+	break;
1029	+	}
1030	+	for (j = i; vlist[j] != NULL; j++)
1031	+	;
1032	+	while (j > i) {
1033	+	vlist[j] = vlist[j-1];
1034	+	--j;
1035	+	}
1036	+	vlist[i] = v;
1037	+	}
1038	+
1039	+	/* Sort triangle edges in standard order */
1040	+	static void
1041	+	order_triangle(MIGRATION *miga[3])
1042	+	{
1043	+	RBFNODE         *vert[4];
1044	+	MIGRATION       *ord[3];
1045	+	int             i;
1046	+	/* order vertices, first */
1047	+	memset(vert, 0, sizeof(vert));
1048	+	for (i = 0; i < 3; i++) {
1049	+	insert_vert(vert, miga[i]->rbfv[0]);
1050	+	insert_vert(vert, miga[i]->rbfv[1]);
1051	+	}
1052	+	/* identify edge 0 */
1053	+	for (i = 0; i < 3; i++)
1054	+	if (miga[i]->rbfv[0] == vert[0] &&
1055	+	miga[i]->rbfv[1] == vert[1]) {
1056	+	ord[0] = miga[i];
1057	+	break;
1058	+	}
1059	+	/* identify edge 1 */
1060	+	for (i = 0; i < 3; i++)
1061	+	if (miga[i]->rbfv[0] == vert[1] &&
1062	+	miga[i]->rbfv[1] == vert[2]) {
1063	+	ord[1] = miga[i];
1064	+	break;
1065	+	}
1066	+	/* identify edge 2 */
1067	+	for (i = 0; i < 3; i++)
1068	+	if (miga[i]->rbfv[0] == vert[0] &&
1069	+	miga[i]->rbfv[1] == vert[2]) {
1070	+	ord[2] = miga[i];
1071	+	break;
1072	+	}
1073	+	miga[0] = ord[0]; miga[1] = ord[1]; miga[2] = ord[2];
1074	+	}
1075	+
1076	+	/* Partially advect between recorded incident angles and allocate new RBF */
1077	+	static RBFNODE *
1078	+	advect_rbf(const FVECT invec)
1079	+	{
1080	+	MIGRATION       *miga[3];
1081	+	RBFNODE         *rbf;
1082	+	float           mbfact, mcfact;
1083	+	int             n, i, j, k;
1084	+	FVECT           v0, v1, v2;
1085	+	double          s, t;
1086	+
1087	+	if (!get_interp(miga, invec))           /* can't interpolate? */
1088	+	return(NULL);
1089	+	if (miga[1] == NULL)                    /* along edge? */
1090	+	return(e_advect_rbf(miga[0], invec));
1091	+	/* put in standard order */
1092	+	order_triangle(miga);
1093	+	#ifdef DEBUG
1094	+	if (miga[0]->rbfv[0] != miga[2]->rbfv[0] |
1095	+	miga[0]->rbfv[1] != miga[1]->rbfv[0] |
1096	+	miga[1]->rbfv[1] != miga[2]->rbfv[1]) {
1097	+	fputs("Triangle vertex screw-up!\n", stderr);
1098	+	exit(1);
1099	+	}
1100	+	#endif
1101	+	/* figure out position */
1102	+	fcross(v0, miga[2]->rbfv[0]->invec, miga[2]->rbfv[1]->invec);
1103	+	normalize(v0);
1104	+	fcross(v2, miga[1]->rbfv[0]->invec, miga[1]->rbfv[1]->invec);
1105	+	normalize(v2);
1106	+	fcross(v1, invec, miga[1]->rbfv[1]->invec);
1107	+	normalize(v1);
1108	+	s = acos(DOT(v0,v1)) / acos(DOT(v0,v2));
1109	+	geodesic(v1, miga[0]->rbfv[0]->invec, miga[0]->rbfv[1]->invec,
1110	+	s, GEOD_REL);
1111	+	t = acos(DOT(v1,invec)) / acos(DOT(v1,miga[1]->rbfv[1]->invec));
1112	+	n = 0;                                  /* count migrating particles */
1113	+	for (i = 0; i < mtx_nrows(miga[0]); i++)
1114	+	for (j = 0; j < mtx_ncols(miga[0]); j++)
1115	+	for (k = (miga[0]->mtx[mtx_ndx(miga[0],i,j)] > FTINY) *
1116	+	mtx_ncols(miga[2]); k--; )
1117	+	n += (miga[2]->mtx[mtx_ndx(miga[2],i,k)] > FTINY &&
1118	+	miga[1]->mtx[mtx_ndx(miga[1],j,k)] > FTINY);
1119	+	#ifdef DEBUG
1120	+	fprintf(stderr, "Input RBFs have %d, %d, %d nodes -> output has %d\n",
1121	+	miga[0]->rbfv[0]->nrbf, miga[0]->rbfv[1]->nrbf,
1122	+	miga[2]->rbfv[1]->nrbf, n);
1123	+	#endif
1124	+	rbf = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(n-1));
1125	+	if (rbf == NULL) {
1126	+	fputs("Out of memory in advect_rbf()\n", stderr);
1127	+	exit(1);
1128	+	}
1129	+	rbf->next = NULL; rbf->ejl = NULL;
1130	+	VCOPY(rbf->invec, invec);
1131	+	rbf->nrbf = n;
1132	+	n = 0;                                  /* compute RBF lobes */
1133	+	mbfact = s * miga[0]->rbfv[1]->vtotal/miga[0]->rbfv[0]->vtotal *
1134	+	(1.-t + t*miga[1]->rbfv[1]->vtotal/miga[1]->rbfv[0]->vtotal);
1135	+	mcfact = (1.-s) *
1136	+	(1.-t + t*miga[2]->rbfv[1]->vtotal/miga[2]->rbfv[0]->vtotal);
1137	+	for (i = 0; i < mtx_nrows(miga[0]); i++) {
1138	+	const RBFVAL        *rbf0i = &miga[0]->rbfv[0]->rbfa[i];
1139	+	const float         w0i = rbf0i->peak;
1140	+	const double        rad0i = R2ANG(rbf0i->crad);
1141	+	ovec_from_pos(v0, rbf0i->gx, rbf0i->gy);
1142	+	for (j = 0; j < mtx_ncols(miga[0]); j++) {
1143	+	const float     ma = miga[0]->mtx[mtx_ndx(miga[0],i,j)];
1144	+	const RBFVAL    *rbf1j;
1145	+	double          rad1j, srad2;
1146	+	if (ma <= FTINY)
1147	+	continue;
1148	+	rbf1j = &miga[0]->rbfv[1]->rbfa[j];
1149	+	rad1j = R2ANG(rbf1j->crad);
1150	+	srad2 = (1.-s)*(1.-t)*rad0i*rad0i + s*(1.-t)*rad1j*rad1j;
1151	+	ovec_from_pos(v1, rbf1j->gx, rbf1j->gy);
1152	+	geodesic(v1, v0, v1, s, GEOD_REL);
1153	+	for (k = 0; k < mtx_ncols(miga[2]); k++) {
1154	+	float               mb = miga[1]->mtx[mtx_ndx(miga[1],j,k)];
1155	+	float               mc = miga[2]->mtx[mtx_ndx(miga[2],i,k)];
1156	+	const RBFVAL        *rbf2k;
1157	+	double              rad2k;
1158	+	FVECT               vout;
1159	+	int                 pos[2];
1160	+	if ((mb <= FTINY) | (mc <= FTINY))
1161	+	continue;
1162	+	rbf2k = &miga[2]->rbfv[1]->rbfa[k];
1163	+	rbf->rbfa[n].peak = w0i * ma * (mb*mbfact + mc*mcfact);
1164	+	rad2k = R2ANG(rbf2k->crad);
1165	+	rbf->rbfa[n].crad = ANG2R(sqrt(srad2 + t*rad2k*rad2k));
1166	+	ovec_from_pos(v2, rbf2k->gx, rbf2k->gy);
1167	+	geodesic(vout, v1, v2, t, GEOD_REL);
1168	+	pos_from_vec(pos, vout);
1169	+	rbf->rbfa[n].gx = pos[0];
1170	+	rbf->rbfa[n].gy = pos[1];
1171	+	++n;
1172	+	}
1173	+	}
1174	+	}
1175	+	rbf->vtotal = miga[0]->rbfv[0]->vtotal * (mbfact + mcfact);
1176	+	return(rbf);
1177	+	}
1178	+
1179	+	/* Interpolate and output isotropic BSDF data */
1180	+	static void
1181	+	interp_isotropic()
1182	+	{
1183	+	const int       sqres = 1<<samp_order;
1184	+	FILE            *ofp = NULL;
1185	+	char            cmd[128];
1186	+	int             ix, ox, oy;
1187	+	FVECT           ivec, ovec;
1188	+	double          bsdf;
1189	+
1190	+	if (pctcull >= 0) {                     /* begin output */
1191	+	sprintf(cmd, "rttree_reduce -h -a -fd -r 3 -t %d -g %d",
1192	+	pctcull, samp_order);
1193	+	fflush(stdout);
1194	+	ofp = popen(cmd, "w");
1195	+	if (ofp == NULL) {
1196	+	fputs("Cannot create pipe for rttree_reduce\n", stderr);
1197	+	exit(1);
1198	+	}
1199	+	} else
1200	+	fputs("{\n", stdout);
1201	+	/* run through directions */
1202	+	for (ix = 0; ix < sqres/2; ix++) {
1203	+	RBFNODE *rbf;
1204	+	SDsquare2disk(ivec, (ix+.5)/sqres, .5);
1205	+	ivec[2] = input_orient *
1206	+	sqrt(1. - ivec[0]*ivec[0] - ivec[1]*ivec[1]);
1207	+	rbf = advect_rbf(ivec);
1208	+	for (ox = 0; ox < sqres; ox++)
1209	+	for (oy = 0; oy < sqres; oy++) {
1210	+	SDsquare2disk(ovec, (ox+.5)/sqres, (oy+.5)/sqres);
1211	+	ovec[2] = output_orient *
1212	+	sqrt(1. - ovec[0]*ovec[0] - ovec[1]*ovec[1]);
1213	+	bsdf = eval_rbfrep(rbf, ovec) / fabs(ovec[2]);
1214	+	if (pctcull >= 0)
1215	+	fwrite(&bsdf, sizeof(bsdf), 1, ofp);
1216	+	else
1217	+	printf("\t%.3e\n", bsdf);
1218	+	}
1219	+	free(rbf);
1220	+	}
1221	+	if (pctcull >= 0) {                     /* finish output */
1222	+	if (pclose(ofp)) {
1223	+	fprintf(stderr, "Error running '%s'\n", cmd);
1224	+	exit(1);
1225	+	}
1226	+	} else {
1227	+	for (ix = sqres*sqres*sqres/2; ix--; )
1228	+	fputs("\t0\n", stdout);
1229	+	fputs("}\n", stdout);
1230	+	}
1231	+	}
1232	+
1233	+	/* Interpolate and output anisotropic BSDF data */
1234	+	static void
1235	+	interp_anisotropic()
1236	+	{
1237	+	const int       sqres = 1<<samp_order;
1238	+	FILE            *ofp = NULL;
1239	+	char            cmd[128];
1240	+	int             ix, iy, ox, oy;
1241	+	FVECT           ivec, ovec;
1242	+	double          bsdf;
1243	+
1244	+	if (pctcull >= 0) {                     /* begin output */
1245	+	sprintf(cmd, "rttree_reduce -h -a -fd -r 4 -t %d -g %d",
1246	+	pctcull, samp_order);
1247	+	fflush(stdout);
1248	+	ofp = popen(cmd, "w");
1249	+	if (ofp == NULL) {
1250	+	fputs("Cannot create pipe for rttree_reduce\n", stderr);
1251	+	exit(1);
1252	+	}
1253	+	} else
1254	+	fputs("{\n", stdout);
1255	+	/* run through directions */
1256	+	for (ix = 0; ix < sqres; ix++)
1257	+	for (iy = 0; iy < sqres; iy++) {
1258	+	RBFNODE *rbf;
1259	+	SDsquare2disk(ivec, (ix+.5)/sqres, (iy+.5)/sqres);
1260	+	ivec[2] = input_orient *
1261	+	sqrt(1. - ivec[0]*ivec[0] - ivec[1]*ivec[1]);
1262	+	rbf = advect_rbf(ivec);
1263	+	for (ox = 0; ox < sqres; ox++)
1264	+	for (oy = 0; oy < sqres; oy++) {
1265	+	SDsquare2disk(ovec, (ox+.5)/sqres, (oy+.5)/sqres);
1266	+	ovec[2] = output_orient *
1267	+	sqrt(1. - ovec[0]*ovec[0] - ovec[1]*ovec[1]);
1268	+	bsdf = eval_rbfrep(rbf, ovec) / fabs(ovec[2]);
1269	+	if (pctcull >= 0)
1270	+	fwrite(&bsdf, sizeof(bsdf), 1, ofp);
1271	+	else
1272	+	printf("\t%.3e\n", bsdf);
1273	+	}
1274	+	free(rbf);
1275	+	}
1276	+	if (pctcull >= 0) {                     /* finish output */
1277	+	if (pclose(ofp)) {
1278	+	fprintf(stderr, "Error running '%s'\n", cmd);
1279	+	exit(1);
1280	+	}
1281	+	} else
1282	+	fputs("}\n", stdout);
1283	+	}
1284	+
1285		#if 1
1286	+	/* Read in BSDF files and interpolate as tensor tree representation */
1287	+	int
1288	+	main(int argc, char *argv[])
1289	+	{
1290	+	RBFNODE         *rbf;
1291	+	double          bsdf;
1292	+	int             i;
1293	+
1294	+	progname = argv[0];
1295	+	if (argc > 2 && !strcmp(argv[1], "-t")) {
1296	+	pctcull = atoi(argv[2]);
1297	+	argv += 2; argc -= 2;
1298	+	}
1299	+	if (argc < 3) {
1300	+	fprintf(stderr,
1301	+	"Usage: %s [-t pctcull] meas1.dat meas2.dat .. > bsdf.xml\n",
1302	+	progname);
1303	+	return(1);
1304	+	}
1305	+	for (i = 1; i < argc; i++) {            /* compile measurements */
1306	+	if (!load_pabopto_meas(argv[i]))
1307	+	return(1);
1308	+	compute_radii();
1309	+	cull_values();
1310	+	make_rbfrep();
1311	+	}
1312	+	build_mesh();                           /* create interpolation */
1313	+	/* xml_prologue();                              /* start XML output */
1314	+	if (single_plane_incident)              /* resample dist. */
1315	+	interp_isotropic();
1316	+	else
1317	+	interp_anisotropic();
1318	+	/* xml_epilogue();                              /* finish XML output */
1319	+	return(0);
1320	+	}
1321	+	#else
1322		/* Test main produces a Radiance model from the given input file */
1323		int
1324		main(int argc, char *argv[])
#	Line 354 | Line 1333 | main(int argc, char *argv[])
1333		fprintf(stderr, "Usage: %s input.dat > output.rad\n", argv[0]);
1334		return(1);
1335		}
1336	<	if (!load_bsdf_meas(argv[1]))
1336	>	if (!load_pabopto_meas(argv[1]))
1337		return(1);
359	–	/* produce spheres at meas. */
360	–	puts("void plastic orange\n0\n0\n5 .6 .4 .01 .04 .08\n");
361	–	n = 0;
362	–	for (i = 0; i < GRIDRES; i++)
363	–	for (j = 0; j < GRIDRES; j++)
364	–	if (bsdf_grid[i][j].nval) {
365	–	double  bsdf = bsdf_grid[i][j].vsum /
366	–	(double)bsdf_grid[i][j].nval;
367	–	FVECT   dir;
1338
369	–	vec_from_pos(dir, i, j);
370	–	printf("orange sphere s%04d\n0\n0\n", ++n);
371	–	printf("4 %.6g %.6g %.6g .0015\n\n",
372	–	dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
373	–	}
1339		compute_radii();
1340		cull_values();
1341	<	/* highlight chosen values */
1341	>	make_rbfrep();
1342	>	/* produce spheres at meas. */
1343	>	puts("void plastic yellow\n0\n0\n5 .6 .4 .01 .04 .08\n");
1344		puts("void plastic pink\n0\n0\n5 .5 .05 .9 .04 .08\n");
1345		n = 0;
1346		for (i = 0; i < GRIDRES; i++)
1347		for (j = 0; j < GRIDRES; j++)
1348	<	if (bsdf_grid[i][j].nval) {
1349	<	bsdf = bsdf_grid[i][j].vsum /
1350	<	(double)bsdf_grid[i][j].nval;
1351	<	vec_from_pos(dir, i, j);
1352	<	printf("pink cone c%04d\n0\n0\n8\n", ++n);
1353	<	printf("\t%.6g %.6g %.6g\n",
1348	>	if (dsf_grid[i][j].vsum > .0f) {
1349	>	ovec_from_pos(dir, i, j);
1350	>	bsdf = dsf_grid[i][j].vsum / dir[2];
1351	>	if (dsf_grid[i][j].nval) {
1352	>	printf("pink cone c%04d\n0\n0\n8\n", ++n);
1353	>	printf("\t%.6g %.6g %.6g\n",
1354		dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
1355	<	printf("\t%.6g %.6g %.6g\n",
1355	>	printf("\t%.6g %.6g %.6g\n",
1356		dir[0]*(bsdf+.005), dir[1]*(bsdf+.005),
1357		dir[2]*(bsdf+.005));
1358	<	puts("\t.003\t0\n");
1358	>	puts("\t.003\t0\n");
1359	>	} else {
1360	>	ovec_from_pos(dir, i, j);
1361	>	printf("yellow sphere s%04d\n0\n0\n", ++n);
1362	>	printf("4 %.6g %.6g %.6g .0015\n\n",
1363	>	dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
1364	>	}
1365		}
1366		/* output continuous surface */
394	–	make_rbfrep();
1367		puts("void trans tgreen\n0\n0\n7 .7 1 .7 .04 .04 .9 .9\n");
1368		fflush(stdout);
1369	<	sprintf(buf, "gensurf tgreen bsdf - - - %d %d", GRIDRES, GRIDRES);
1369	>	sprintf(buf, "gensurf tgreen bsdf - - - %d %d", GRIDRES-1, GRIDRES-1);
1370		pfp = popen(buf, "w");
1371		if (pfp == NULL) {
1372		fputs(buf, stderr);
#	Line 403 | Line 1375 | main(int argc, char *argv[])
1375		}
1376		for (i = 0; i < GRIDRES; i++)
1377		for (j = 0; j < GRIDRES; j++) {
1378	<	vec_from_pos(dir, i, j);
1379	<	bsdf = eval_rbfrep(bsdf_list, dir);
1378	>	ovec_from_pos(dir, i, j);
1379	>	bsdf = eval_rbfrep(dsf_list, dir) / dir[2];
1380		fprintf(pfp, "%.8e %.8e %.8e\n",
1381		dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
1382		}

Diff Legend

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