ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/pabopto2xml.c
Revision: 2.15
Committed: Sun Sep 23 16:45:20 2012 UTC (11 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.14: +124 -69 lines
Log Message: 
Added checks for bad triangles in get_interp()

File Contents

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