ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/pabopto2xml.c
(Generate patch)

Comparing ray/src/cv/pabopto2xml.c (file contents):
Revision 2.5 by greg, Sat Aug 25 22:39:03 2012 UTC vs.
Revision 2.12 by greg, Thu Sep 20 01:23:36 2012 UTC

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

Diff Legend

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