ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/pabopto2xml.c
Revision: 2.12
Committed: Thu Sep 20 01:23:36 2012 UTC (11 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.11: +176 -5 lines
Log Message: 
Nearing completion -- needs prologue, epilogue and testing

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.12 static const char RCSid[] = "$Id: pabopto2xml.c,v 2.11 2012/09/19 22:03:37 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     /* percentage to cull (<0 to turn off) */
93     int pctcull = 90;
94     /* sampling order */
95     int samp_order = 0;
96    
97 greg 2.8 /* 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 greg 2.1
106 greg 2.3 /* Compute outgoing vector from grid position */
107     static void
108 greg 2.10 ovec_from_pos(FVECT vec, int xpos, int ypos)
109 greg 2.1 {
110 greg 2.3 double uv[2];
111     double r2;
112    
113     SDsquare2disk(uv, (1./GRIDRES)*(xpos+.5), (1./GRIDRES)*(ypos+.5));
114     /* uniform hemispherical projection */
115     r2 = uv[0]*uv[0] + uv[1]*uv[1];
116     vec[0] = vec[1] = sqrt(2. - r2);
117     vec[0] *= uv[0];
118     vec[1] *= uv[1];
119 greg 2.10 vec[2] = output_orient*(1. - r2);
120 greg 2.1 }
121    
122 greg 2.10 /* Compute grid position from normalized input/output vector */
123 greg 2.1 static void
124     pos_from_vec(int pos[2], const FVECT vec)
125     {
126     double sq[2]; /* uniform hemispherical projection */
127 greg 2.10 double norm = 1./sqrt(1. + fabs(vec[2]));
128 greg 2.1
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 greg 2.5 /* Evaluate RBF for DSF at the given normalized outgoing direction */
136 greg 2.1 static double
137 greg 2.10 eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
138 greg 2.1 {
139     double res = .0;
140     const RBFVAL *rbfp;
141     FVECT odir;
142     double sig2;
143     int n;
144    
145 greg 2.12 if (rp == NULL)
146     return(.0);
147 greg 2.1 rbfp = rp->rbfa;
148     for (n = rp->nrbf; n--; rbfp++) {
149 greg 2.10 ovec_from_pos(odir, rbfp->gx, rbfp->gy);
150 greg 2.2 sig2 = R2ANG(rbfp->crad);
151     sig2 = (DOT(odir,outvec) - 1.) / (sig2*sig2);
152 greg 2.1 if (sig2 > -19.)
153 greg 2.5 res += rbfp->peak * exp(sig2);
154 greg 2.1 }
155     return(res);
156     }
157    
158 greg 2.10 /* 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 greg 2.3 /* Count up filled nodes and build RBF representation from current grid */
180 greg 2.10 static RBFNODE *
181 greg 2.3 make_rbfrep(void)
182     {
183 greg 2.6 int niter = 16;
184 greg 2.12 int minrad = ANG2R(pow(2., 1.-samp_order));
185 greg 2.6 double lastVar, thisVar = 100.;
186 greg 2.3 int nn;
187 greg 2.10 RBFNODE *newnode;
188 greg 2.3 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 greg 2.6 nn += dsf_grid[i][j].nval;
194 greg 2.3 /* allocate RBF array */
195 greg 2.10 newnode = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(nn-1));
196 greg 2.3 if (newnode == NULL) {
197 greg 2.8 fputs("Out of memory in make_rbfrep()\n", stderr);
198 greg 2.3 exit(1);
199     }
200     newnode->next = NULL;
201 greg 2.7 newnode->ejl = NULL;
202 greg 2.3 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 greg 2.10 newnode->invec[2] = input_orient*sqrt(1. - newnode->invec[2]*newnode->invec[2]);
206 greg 2.8 newnode->vtotal = 0;
207 greg 2.3 newnode->nrbf = nn;
208     nn = 0; /* fill RBF array */
209     for (i = 0; i < GRIDRES; i++)
210     for (j = 0; j < GRIDRES; j++)
211 greg 2.5 if (dsf_grid[i][j].nval) {
212 greg 2.6 newnode->rbfa[nn].peak = dsf_grid[i][j].vsum;
213 greg 2.5 newnode->rbfa[nn].crad = RSCA*dsf_grid[i][j].crad + .5;
214 greg 2.3 newnode->rbfa[nn].gx = i;
215     newnode->rbfa[nn].gy = j;
216 greg 2.12 if (newnode->rbfa[nn].crad < minrad)
217     minrad = newnode->rbfa[nn].crad;
218 greg 2.3 ++nn;
219     }
220 greg 2.6 /* iterate to improve interpolation accuracy */
221     do {
222 greg 2.4 double dsum = .0, dsum2 = .0;
223 greg 2.3 nn = 0;
224     for (i = 0; i < GRIDRES; i++)
225     for (j = 0; j < GRIDRES; j++)
226 greg 2.5 if (dsf_grid[i][j].nval) {
227 greg 2.3 FVECT odir;
228 greg 2.6 double corr;
229 greg 2.10 ovec_from_pos(odir, i, j);
230 greg 2.6 newnode->rbfa[nn++].peak *= corr =
231 greg 2.5 dsf_grid[i][j].vsum /
232 greg 2.3 eval_rbfrep(newnode, odir);
233 greg 2.4 dsum += corr - 1.;
234     dsum2 += (corr-1.)*(corr-1.);
235 greg 2.3 }
236 greg 2.6 lastVar = thisVar;
237     thisVar = dsum2/(double)nn;
238 greg 2.10 #ifdef DEBUG
239 greg 2.4 fprintf(stderr, "Avg., RMS error: %.1f%% %.1f%%\n",
240     100.*dsum/(double)nn,
241 greg 2.6 100.*sqrt(thisVar));
242 greg 2.10 #endif
243 greg 2.6 } while (--niter > 0 && lastVar-thisVar > 0.02*lastVar);
244    
245 greg 2.8 nn = 0; /* compute sum for normalization */
246     while (nn < newnode->nrbf)
247     newnode->vtotal += rbf_volume(&newnode->rbfa[nn++]);
248    
249 greg 2.10 insert_dsf(newnode);
250 greg 2.12 /* adjust sampling resolution */
251     samp_order = log(2./R2ANG(minrad))/log(2.) + .5;
252    
253 greg 2.10 return(newnode);
254 greg 2.3 }
255    
256 greg 2.1 /* Load a set of measurements corresponding to a particular incident angle */
257     static int
258 greg 2.10 load_pabopto_meas(const char *fname)
259 greg 2.1 {
260     FILE *fp = fopen(fname, "r");
261     int inp_is_DSF = -1;
262 greg 2.10 double new_phi, theta_out, phi_out, val;
263 greg 2.1 char buf[2048];
264     int n, c;
265    
266     if (fp == NULL) {
267     fputs(fname, stderr);
268     fputs(": cannot open\n", stderr);
269     return(0);
270     }
271 greg 2.5 memset(dsf_grid, 0, sizeof(dsf_grid));
272 greg 2.10 #ifdef DEBUG
273     fprintf(stderr, "Loading measurement file '%s'...\n", fname);
274     #endif
275 greg 2.1 /* read header information */
276     while ((c = getc(fp)) == '#' || c == EOF) {
277     if (fgets(buf, sizeof(buf), fp) == NULL) {
278     fputs(fname, stderr);
279     fputs(": unexpected EOF\n", stderr);
280     fclose(fp);
281     return(0);
282     }
283     if (!strcmp(buf, "format: theta phi DSF\n")) {
284     inp_is_DSF = 1;
285     continue;
286     }
287     if (!strcmp(buf, "format: theta phi BSDF\n")) {
288     inp_is_DSF = 0;
289     continue;
290     }
291     if (sscanf(buf, "intheta %lf", &theta_in_deg) == 1)
292     continue;
293 greg 2.10 if (sscanf(buf, "inphi %lf", &new_phi) == 1)
294 greg 2.1 continue;
295     if (sscanf(buf, "incident_angle %lf %lf",
296 greg 2.10 &theta_in_deg, &new_phi) == 2)
297 greg 2.1 continue;
298     }
299     if (inp_is_DSF < 0) {
300     fputs(fname, stderr);
301     fputs(": unknown format\n", stderr);
302     fclose(fp);
303     return(0);
304     }
305 greg 2.10 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 greg 2.1 while (fscanf(fp, "%lf %lf %lf\n", &theta_out, &phi_out, &val) == 3) {
319     FVECT ovec;
320     int pos[2];
321    
322 greg 2.10 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 greg 2.1 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 greg 2.5 if (!inp_is_DSF)
335     val *= ovec[2]; /* convert from BSDF to DSF */
336 greg 2.1
337     pos_from_vec(pos, ovec);
338    
339 greg 2.5 dsf_grid[pos[0]][pos[1]].vsum += val;
340     dsf_grid[pos[0]][pos[1]].nval++;
341 greg 2.1 }
342     n = 0;
343     while ((c = getc(fp)) != EOF)
344     n += !isspace(c);
345     if (n)
346     fprintf(stderr,
347     "%s: warning: %d unexpected characters past EOD\n",
348     fname, n);
349     fclose(fp);
350     return(1);
351     }
352    
353     /* Compute radii for non-empty bins */
354     /* (distance to furthest empty bin for which non-empty bin is the closest) */
355     static void
356     compute_radii(void)
357     {
358 greg 2.4 unsigned int fill_grid[GRIDRES][GRIDRES];
359     unsigned short fill_cnt[GRIDRES][GRIDRES];
360 greg 2.2 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 greg 2.1 for (i = 0; i < GRIDRES; i++)
366     for (jn = 0; jn < GRIDRES; jn++) {
367     j = (i&1) ? jn : GRIDRES-1-jn;
368 greg 2.5 if (dsf_grid[i][j].nval) /* find empty grid pos. */
369 greg 2.1 continue;
370 greg 2.10 ovec_from_pos(ovec0, i, j);
371 greg 2.1 inear = jnear = -1; /* find nearest non-empty */
372 greg 2.2 lastang2 = M_PI*M_PI;
373 greg 2.1 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 greg 2.5 if (!dsf_grid[ii][jj].nval)
380 greg 2.1 continue;
381 greg 2.10 ovec_from_pos(ovec1, ii, jj);
382 greg 2.2 ang2 = 2. - 2.*DOT(ovec0,ovec1);
383     if (ang2 >= lastang2)
384 greg 2.1 continue;
385 greg 2.2 lastang2 = ang2;
386 greg 2.1 inear = ii; jnear = jj;
387     }
388     }
389 greg 2.2 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 greg 2.5 if (r > dsf_grid[inear][jnear].crad)
396     dsf_grid[inear][jnear].crad = r;
397 greg 2.2 /* next search radius */
398 greg 2.10 r = ang2*(2.*GRIDRES/M_PI) + 3;
399 greg 2.1 }
400 greg 2.4 /* blur radii over hemisphere */
401 greg 2.1 memset(fill_grid, 0, sizeof(fill_grid));
402 greg 2.4 memset(fill_cnt, 0, sizeof(fill_cnt));
403 greg 2.1 for (i = 0; i < GRIDRES; i++)
404     for (j = 0; j < GRIDRES; j++) {
405 greg 2.5 if (!dsf_grid[i][j].crad)
406 greg 2.4 continue; /* missing distance */
407 greg 2.5 r = R2ANG(dsf_grid[i][j].crad)*(2.*RSCA*GRIDRES/M_PI);
408 greg 2.1 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 greg 2.4 if ((ii-i)*(ii-i) + (jj-j)*(jj-j) > r*r)
415 greg 2.1 continue;
416 greg 2.5 fill_grid[ii][jj] += dsf_grid[i][j].crad;
417 greg 2.4 fill_cnt[ii][jj]++;
418 greg 2.1 }
419     }
420     }
421 greg 2.6 /* copy back blurred radii */
422 greg 2.1 for (i = 0; i < GRIDRES; i++)
423     for (j = 0; j < GRIDRES; j++)
424 greg 2.4 if (fill_cnt[i][j])
425 greg 2.5 dsf_grid[i][j].crad = fill_grid[i][j]/fill_cnt[i][j];
426 greg 2.1 }
427    
428 greg 2.6 /* Cull points for more uniform distribution, leave all nval 0 or 1 */
429 greg 2.1 static void
430     cull_values(void)
431     {
432 greg 2.2 FVECT ovec0, ovec1;
433     double maxang, maxang2;
434     int i, j, ii, jj, r;
435 greg 2.1 /* simple greedy algorithm */
436     for (i = 0; i < GRIDRES; i++)
437     for (j = 0; j < GRIDRES; j++) {
438 greg 2.5 if (!dsf_grid[i][j].nval)
439 greg 2.1 continue;
440 greg 2.5 if (!dsf_grid[i][j].crad)
441 greg 2.2 continue; /* shouldn't happen */
442 greg 2.10 ovec_from_pos(ovec0, i, j);
443 greg 2.5 maxang = 2.*R2ANG(dsf_grid[i][j].crad);
444 greg 2.2 if (maxang > ovec0[2]) /* clamp near horizon */
445     maxang = ovec0[2];
446     r = maxang*(2.*GRIDRES/M_PI) + 1;
447     maxang2 = maxang*maxang;
448 greg 2.1 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 greg 2.5 if (!dsf_grid[ii][jj].nval)
455 greg 2.1 continue;
456 greg 2.2 if ((ii == i) & (jj == j))
457     continue; /* don't get self-absorbed */
458 greg 2.10 ovec_from_pos(ovec1, ii, jj);
459 greg 2.2 if (2. - 2.*DOT(ovec0,ovec1) >= maxang2)
460 greg 2.1 continue;
461 greg 2.2 /* absorb sum */
462 greg 2.5 dsf_grid[i][j].vsum += dsf_grid[ii][jj].vsum;
463     dsf_grid[i][j].nval += dsf_grid[ii][jj].nval;
464 greg 2.2 /* keep value, though */
465 greg 2.6 dsf_grid[ii][jj].vsum /= (float)dsf_grid[ii][jj].nval;
466 greg 2.5 dsf_grid[ii][jj].nval = 0;
467 greg 2.1 }
468     }
469     }
470 greg 2.6 /* 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 greg 2.1 }
478    
479 greg 2.8 /* Compute (and allocate) migration price matrix for optimization */
480     static float *
481 greg 2.10 price_routes(const RBFNODE *from_rbf, const RBFNODE *to_rbf)
482 greg 2.8 {
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 greg 2.10 ovec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);
494 greg 2.8
495     for (i = from_rbf->nrbf; i--; ) {
496     const double from_ang = R2ANG(from_rbf->rbfa[i].crad);
497     FVECT vfrom;
498 greg 2.10 ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
499 greg 2.8 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 greg 2.9 const double maxamt = 2./(mtx_nrows(mig)*mtx_ncols(mig));
562 greg 2.8 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 greg 2.10 make_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
626 greg 2.8 {
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 greg 2.10 #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 greg 2.8 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 greg 2.10 /* 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 greg 2.8 {
737 greg 2.10 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 greg 2.8
790 greg 2.10 /* 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 greg 2.8 exit(1);
861     }
862 greg 2.10 /* 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 greg 2.12 { /* else use triangle mesh */
927 greg 2.10 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 greg 2.12 #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 greg 2.10 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 greg 2.11 /* 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 greg 2.10 /* 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 greg 2.11 float mbfact, mcfact;
1083     int n, i, j, k;
1084     FVECT v0, v1, v2;
1085 greg 2.10 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 greg 2.11 /* put in standard order */
1092     order_triangle(miga);
1093 greg 2.12 #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 greg 2.10 /* figure out position */
1102 greg 2.11 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 greg 2.12 #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 greg 2.10 rbf = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(n-1));
1125 greg 2.8 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 greg 2.10 rbf->nrbf = n;
1132 greg 2.11 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 greg 2.12 rbf->rbfa[n].crad = ANG2R(sqrt(srad2 + t*rad2k*rad2k));
1166 greg 2.11 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 greg 2.10 return(rbf);
1177 greg 2.8 }
1178 greg 2.1
1179 greg 2.12 /* 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 greg 2.1 #if 1
1286 greg 2.12 /* 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 greg 2.1 /* Test main produces a Radiance model from the given input file */
1323     int
1324     main(int argc, char *argv[])
1325     {
1326     char buf[128];
1327     FILE *pfp;
1328     double bsdf;
1329     FVECT dir;
1330     int i, j, n;
1331    
1332     if (argc != 2) {
1333     fprintf(stderr, "Usage: %s input.dat > output.rad\n", argv[0]);
1334     return(1);
1335     }
1336 greg 2.10 if (!load_pabopto_meas(argv[1]))
1337 greg 2.1 return(1);
1338    
1339     compute_radii();
1340     cull_values();
1341 greg 2.3 make_rbfrep();
1342     /* produce spheres at meas. */
1343     puts("void plastic yellow\n0\n0\n5 .6 .4 .01 .04 .08\n");
1344 greg 2.1 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 greg 2.5 if (dsf_grid[i][j].vsum > .0f) {
1349 greg 2.10 ovec_from_pos(dir, i, j);
1350 greg 2.5 bsdf = dsf_grid[i][j].vsum / dir[2];
1351     if (dsf_grid[i][j].nval) {
1352 greg 2.3 printf("pink cone c%04d\n0\n0\n8\n", ++n);
1353     printf("\t%.6g %.6g %.6g\n",
1354 greg 2.1 dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
1355 greg 2.3 printf("\t%.6g %.6g %.6g\n",
1356 greg 2.1 dir[0]*(bsdf+.005), dir[1]*(bsdf+.005),
1357     dir[2]*(bsdf+.005));
1358 greg 2.3 puts("\t.003\t0\n");
1359     } else {
1360 greg 2.10 ovec_from_pos(dir, i, j);
1361 greg 2.3 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 greg 2.1 }
1366     /* output continuous surface */
1367     puts("void trans tgreen\n0\n0\n7 .7 1 .7 .04 .04 .9 .9\n");
1368     fflush(stdout);
1369 greg 2.5 sprintf(buf, "gensurf tgreen bsdf - - - %d %d", GRIDRES-1, GRIDRES-1);
1370 greg 2.1 pfp = popen(buf, "w");
1371     if (pfp == NULL) {
1372     fputs(buf, stderr);
1373     fputs(": cannot start command\n", stderr);
1374     return(1);
1375     }
1376     for (i = 0; i < GRIDRES; i++)
1377     for (j = 0; j < GRIDRES; j++) {
1378 greg 2.10 ovec_from_pos(dir, i, j);
1379 greg 2.5 bsdf = eval_rbfrep(dsf_list, dir) / dir[2];
1380 greg 2.1 fprintf(pfp, "%.8e %.8e %.8e\n",
1381     dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
1382     }
1383     return(pclose(pfp)==0 ? 0 : 1);
1384     }
1385     #endif