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.2 by greg, Fri Aug 24 20:55:28 2012 UTC vs.
Revision 2.9 by greg, Thu Sep 6 00:07:43 2012 UTC

# Line 20 | Line 20 | static const char RCSid[] = "$Id$";
20   #define GRIDRES         200             /* max. grid resolution per side */
21   #endif
22  
23 < #define RSCA            3.              /* radius scaling factor (empirical) */
24 < #define MSCA            .2              /* magnitude scaling (empirical) */
23 > #define RSCA            2.7             /* radius scaling factor (empirical) */
24  
25 < #define R2ANG(c)        (((c)+.5)*(M_PI/(1<<16)))
25 >                                        /* convert to/from coded radians */
26   #define ANG2R(r)        (int)((r)*((1<<16)/M_PI))
27 + #define R2ANG(c)        (((c)+.5)*(M_PI/(1<<16)))
28  
29   typedef struct {
30 <        float           vsum;           /* BSDF sum */
30 >        float           vsum;           /* DSF sum */
31          unsigned short  nval;           /* number of values in sum */
32          unsigned short  crad;           /* radius (coded angle) */
33   } GRIDVAL;                      /* grid value */
34  
35   typedef struct {
36 <        float           bsdf;           /* lobe value at peak */
36 >        float           peak;           /* lobe value at peak */
37          unsigned short  crad;           /* radius (coded angle) */
38          unsigned char   gx, gy;         /* grid position */
39   } RBFVAL;                       /* radial basis function value */
40  
41 < typedef struct s_rbflist {
42 <        struct s_rbflist        *next;          /* next in our RBF list */
41 > struct s_rbfnode;               /* forward declaration of RBF struct */
42 >
43 > typedef struct s_migration {
44 >        struct s_migration      *next;          /* next in global edge list */
45 >        struct s_rbfnode        *rbfv[2];       /* from,to vertex */
46 >        struct s_migration      *enxt[2];       /* next from,to sibling */
47 >        float                   mtx[1];         /* matrix (extends struct) */
48 > } MIGRATION;                    /* migration link (winged edge structure) */
49 >
50 > typedef struct s_rbfnode {
51 >        struct s_rbfnode        *next;          /* next in global RBF list */
52 >        MIGRATION               *ejl;           /* edge list for this vertex */
53          FVECT                   invec;          /* incident vector direction */
54 +        double                  vtotal;         /* volume for normalization */
55          int                     nrbf;           /* number of RBFs */
56          RBFVAL                  rbfa[1];        /* RBF array (extends struct) */
57 < } RBFLIST;                      /* RBF representation of BSDF @ 1 incidence */
57 > } RBFLIST;                      /* RBF representation of DSF @ 1 incidence */
58  
59                                  /* our loaded grid for this incident angle */
60   static double   theta_in_deg, phi_in_deg;
61 < static GRIDVAL  bsdf_grid[GRIDRES][GRIDRES];
61 > static GRIDVAL  dsf_grid[GRIDRES][GRIDRES];
62  
63 <                                /* processed incident BSDF measurements */
64 < static RBFLIST  *bsdf_list = NULL;
63 >                                /* processed incident DSF measurements */
64 > static RBFLIST          *dsf_list = NULL;
65  
66 < /* Count up non-empty nodes and build RBF representation from current grid */
67 < static RBFLIST *
57 < make_rbfrep(void)
58 < {
59 <        int     nn = 0;
60 <        RBFLIST *newnode;
61 <        int     i, j;
62 <                                /* count non-empty bins */
63 <        for (i = 0; i < GRIDRES; i++)
64 <            for (j = 0; j < GRIDRES; j++)
65 <                nn += (bsdf_grid[i][j].nval > 0);
66 <                                /* allocate RBF array */
67 <        newnode = (RBFLIST *)malloc(sizeof(RBFLIST) + sizeof(RBFVAL)*(nn-1));
68 <        if (newnode == NULL) {
69 <                fputs("Out of memory in make_rbfrep\n", stderr);
70 <                exit(1);
71 <        }
72 <        newnode->invec[2] = sin(M_PI/180.*theta_in_deg);
73 <        newnode->invec[0] = cos(M_PI/180.*phi_in_deg)*newnode->invec[2];
74 <        newnode->invec[1] = sin(M_PI/180.*phi_in_deg)*newnode->invec[2];
75 <        newnode->invec[2] = sqrt(1. - newnode->invec[2]*newnode->invec[2]);
76 <        newnode->nrbf = nn;
77 <        nn = 0;                 /* fill RBF array */
78 <        for (i = 0; i < GRIDRES; i++)
79 <            for (j = 0; j < GRIDRES; j++)
80 <                if (bsdf_grid[i][j].nval) {
81 <                        newnode->rbfa[nn].bsdf = MSCA*bsdf_grid[i][j].vsum /
82 <                                                (double)bsdf_grid[i][j].nval;
83 <                        newnode->rbfa[nn].crad = RSCA*bsdf_grid[i][j].crad + .5;
84 <                        newnode->rbfa[nn].gx = i;
85 <                        newnode->rbfa[nn].gy = j;
86 <                        ++nn;
87 <                }
88 <        newnode->next = bsdf_list;
89 <        return(bsdf_list = newnode);
90 < }
66 >                                /* RBF-linking matrices (edges) */
67 > static MIGRATION        *mig_list = NULL;
68  
69 < /* Compute grid position from normalized outgoing vector */
70 < static void
71 < pos_from_vec(int pos[2], const FVECT vec)
69 > #define mtx_nrows(m)    ((m)->rbfv[0]->nrbf)
70 > #define mtx_ncols(m)    ((m)->rbfv[1]->nrbf)
71 > #define mtx_ndx(m,i,j)  ((i)*mtx_ncols(m) + (j))
72 > #define is_src(rbf,m)   ((rbf) == (m)->rbfv[0])
73 > #define is_dest(rbf,m)  ((rbf) == (m)->rbfv[1])
74 > #define nextedge(rbf,m) (m)->enxt[is_dest(rbf,m)]
75 >
76 > /* Compute volume associated with Gaussian lobe */
77 > static double
78 > rbf_volume(const RBFVAL *rbfp)
79   {
80 <        double  sq[2];          /* uniform hemispherical projection */
97 <        double  norm = 1./sqrt(1. + vec[2]);
80 >        double  rad = R2ANG(rbfp->crad);
81  
82 <        SDdisk2square(sq, vec[0]*norm, vec[1]*norm);
100 <
101 <        pos[0] = (int)(sq[0]*GRIDRES);
102 <        pos[1] = (int)(sq[1]*GRIDRES);
82 >        return((2.*M_PI) * rbfp->peak * rad*rad);
83   }
84  
85   /* Compute outgoing vector from grid position */
# Line 118 | Line 98 | vec_from_pos(FVECT vec, int xpos, int ypos)
98          vec[2] = 1. - r2;
99   }
100  
101 < /* Evaluate RBF for BSDF at the given normalized outgoing direction */
101 > /* Compute grid position from normalized outgoing vector */
102 > static void
103 > pos_from_vec(int pos[2], const FVECT vec)
104 > {
105 >        double  sq[2];          /* uniform hemispherical projection */
106 >        double  norm = 1./sqrt(1. + vec[2]);
107 >
108 >        SDdisk2square(sq, vec[0]*norm, vec[1]*norm);
109 >
110 >        pos[0] = (int)(sq[0]*GRIDRES);
111 >        pos[1] = (int)(sq[1]*GRIDRES);
112 > }
113 >
114 > /* Evaluate RBF for DSF at the given normalized outgoing direction */
115   static double
116   eval_rbfrep(const RBFLIST *rp, const FVECT outvec)
117   {
# Line 134 | Line 127 | eval_rbfrep(const RBFLIST *rp, const FVECT outvec)
127                  sig2 = R2ANG(rbfp->crad);
128                  sig2 = (DOT(odir,outvec) - 1.) / (sig2*sig2);
129                  if (sig2 > -19.)
130 <                        res += rbfp->bsdf * exp(sig2);
130 >                        res += rbfp->peak * exp(sig2);
131          }
132          return(res);
133   }
134  
135 + /* Count up filled nodes and build RBF representation from current grid */
136 + static RBFLIST *
137 + make_rbfrep(void)
138 + {
139 +        int     niter = 16;
140 +        double  lastVar, thisVar = 100.;
141 +        int     nn;
142 +        RBFLIST *newnode;
143 +        int     i, j;
144 +        
145 +        nn = 0;                 /* count selected bins */
146 +        for (i = 0; i < GRIDRES; i++)
147 +            for (j = 0; j < GRIDRES; j++)
148 +                nn += dsf_grid[i][j].nval;
149 +                                /* allocate RBF array */
150 +        newnode = (RBFLIST *)malloc(sizeof(RBFLIST) + sizeof(RBFVAL)*(nn-1));
151 +        if (newnode == NULL) {
152 +                fputs("Out of memory in make_rbfrep()\n", stderr);
153 +                exit(1);
154 +        }
155 +        newnode->next = NULL;
156 +        newnode->ejl = NULL;
157 +        newnode->invec[2] = sin(M_PI/180.*theta_in_deg);
158 +        newnode->invec[0] = cos(M_PI/180.*phi_in_deg)*newnode->invec[2];
159 +        newnode->invec[1] = sin(M_PI/180.*phi_in_deg)*newnode->invec[2];
160 +        newnode->invec[2] = sqrt(1. - newnode->invec[2]*newnode->invec[2]);
161 +        newnode->vtotal = 0;
162 +        newnode->nrbf = nn;
163 +        nn = 0;                 /* fill RBF array */
164 +        for (i = 0; i < GRIDRES; i++)
165 +            for (j = 0; j < GRIDRES; j++)
166 +                if (dsf_grid[i][j].nval) {
167 +                        newnode->rbfa[nn].peak = dsf_grid[i][j].vsum;
168 +                        newnode->rbfa[nn].crad = RSCA*dsf_grid[i][j].crad + .5;
169 +                        newnode->rbfa[nn].gx = i;
170 +                        newnode->rbfa[nn].gy = j;
171 +                        ++nn;
172 +                }
173 +                                /* iterate to improve interpolation accuracy */
174 +        do {
175 +                double  dsum = .0, dsum2 = .0;
176 +                nn = 0;
177 +                for (i = 0; i < GRIDRES; i++)
178 +                    for (j = 0; j < GRIDRES; j++)
179 +                        if (dsf_grid[i][j].nval) {
180 +                                FVECT   odir;
181 +                                double  corr;
182 +                                vec_from_pos(odir, i, j);
183 +                                newnode->rbfa[nn++].peak *= corr =
184 +                                        dsf_grid[i][j].vsum /
185 +                                                eval_rbfrep(newnode, odir);
186 +                                dsum += corr - 1.;
187 +                                dsum2 += (corr-1.)*(corr-1.);
188 +                        }
189 +                lastVar = thisVar;
190 +                thisVar = dsum2/(double)nn;
191 +                /*
192 +                fprintf(stderr, "Avg., RMS error: %.1f%%  %.1f%%\n",
193 +                                        100.*dsum/(double)nn,
194 +                                        100.*sqrt(thisVar));
195 +                */
196 +        } while (--niter > 0 && lastVar-thisVar > 0.02*lastVar);
197 +
198 +        nn = 0;                 /* compute sum for normalization */
199 +        while (nn < newnode->nrbf)
200 +                newnode->vtotal += rbf_volume(&newnode->rbfa[nn++]);
201 +
202 +        newnode->next = dsf_list;
203 +        return(dsf_list = newnode);
204 + }
205 +
206   /* Load a set of measurements corresponding to a particular incident angle */
207   static int
208   load_bsdf_meas(const char *fname)
# Line 154 | Line 218 | load_bsdf_meas(const char *fname)
218                  fputs(": cannot open\n", stderr);
219                  return(0);
220          }
221 <        memset(bsdf_grid, 0, sizeof(bsdf_grid));
221 >        memset(dsf_grid, 0, sizeof(dsf_grid));
222                                  /* read header information */
223          while ((c = getc(fp)) == '#' || c == EOF) {
224                  if (fgets(buf, sizeof(buf), fp) == NULL) {
# Line 195 | Line 259 | load_bsdf_meas(const char *fname)
259                  ovec[1] = sin(M_PI/180.*phi_out) * ovec[2];
260                  ovec[2] = sqrt(1. - ovec[2]*ovec[2]);
261  
262 <                if (inp_is_DSF)
263 <                        val /= ovec[2]; /* convert from DSF to BSDF */
262 >                if (!inp_is_DSF)
263 >                        val *= ovec[2]; /* convert from BSDF to DSF */
264  
265                  pos_from_vec(pos, ovec);
266  
267 <                bsdf_grid[pos[0]][pos[1]].vsum += val;
268 <                bsdf_grid[pos[0]][pos[1]].nval++;
267 >                dsf_grid[pos[0]][pos[1]].vsum += val;
268 >                dsf_grid[pos[0]][pos[1]].nval++;
269          }
270          n = 0;
271          while ((c = getc(fp)) != EOF)
# Line 219 | Line 283 | load_bsdf_meas(const char *fname)
283   static void
284   compute_radii(void)
285   {
286 <        unsigned short  fill_grid[GRIDRES][GRIDRES];
286 >        unsigned int    fill_grid[GRIDRES][GRIDRES];
287 >        unsigned short  fill_cnt[GRIDRES][GRIDRES];
288          FVECT           ovec0, ovec1;
289          double          ang2, lastang2;
225 –        int             r2, lastr2;
290          int             r, i, j, jn, ii, jj, inear, jnear;
291  
292          r = GRIDRES/2;                          /* proceed in zig-zag */
293          for (i = 0; i < GRIDRES; i++)
294              for (jn = 0; jn < GRIDRES; jn++) {
295                  j = (i&1) ? jn : GRIDRES-1-jn;
296 <                if (bsdf_grid[i][j].nval)       /* find empty grid pos. */
296 >                if (dsf_grid[i][j].nval)        /* find empty grid pos. */
297                          continue;
298                  vec_from_pos(ovec0, i, j);
299                  inear = jnear = -1;             /* find nearest non-empty */
# Line 240 | Line 304 | compute_radii(void)
304                      for (jj = j-r; jj <= j+r; jj++) {
305                          if (jj < 0) continue;
306                          if (jj >= GRIDRES) break;
307 <                        if (!bsdf_grid[ii][jj].nval)
307 >                        if (!dsf_grid[ii][jj].nval)
308                                  continue;
309                          vec_from_pos(ovec1, ii, jj);
310                          ang2 = 2. - 2.*DOT(ovec0,ovec1);
# Line 256 | Line 320 | compute_radii(void)
320                  }
321                  ang2 = sqrt(lastang2);
322                  r = ANG2R(ang2);                /* record if > previous */
323 <                if (r > bsdf_grid[inear][jnear].crad)
324 <                        bsdf_grid[inear][jnear].crad = r;
323 >                if (r > dsf_grid[inear][jnear].crad)
324 >                        dsf_grid[inear][jnear].crad = r;
325                                                  /* next search radius */
326                  r = ang2*(2.*GRIDRES/M_PI) + 1;
327              }
328 <                                                /* fill in neighbors */
328 >                                                /* blur radii over hemisphere */
329          memset(fill_grid, 0, sizeof(fill_grid));
330 +        memset(fill_cnt, 0, sizeof(fill_cnt));
331          for (i = 0; i < GRIDRES; i++)
332              for (j = 0; j < GRIDRES; j++) {
333 <                if (!bsdf_grid[i][j].nval)
334 <                        continue;               /* no value -- skip */
335 <                if (bsdf_grid[i][j].crad)
271 <                        continue;               /* has distance already */
272 <                r = GRIDRES/20;
273 <                lastr2 = 2*r*r + 1;
333 >                if (!dsf_grid[i][j].crad)
334 >                        continue;               /* missing distance */
335 >                r = R2ANG(dsf_grid[i][j].crad)*(2.*RSCA*GRIDRES/M_PI);
336                  for (ii = i-r; ii <= i+r; ii++) {
337                      if (ii < 0) continue;
338                      if (ii >= GRIDRES) break;
339                      for (jj = j-r; jj <= j+r; jj++) {
340                          if (jj < 0) continue;
341                          if (jj >= GRIDRES) break;
342 <                        if (!bsdf_grid[ii][jj].crad)
342 >                        if ((ii-i)*(ii-i) + (jj-j)*(jj-j) > r*r)
343                                  continue;
344 <                                                /* OK to use approx. closest */
345 <                        r2 = (ii-i)*(ii-i) + (jj-j)*(jj-j);
284 <                        if (r2 >= lastr2)
285 <                                continue;
286 <                        fill_grid[i][j] = bsdf_grid[ii][jj].crad;
287 <                        lastr2 = r2;
344 >                        fill_grid[ii][jj] += dsf_grid[i][j].crad;
345 >                        fill_cnt[ii][jj]++;
346                      }
347                  }
348              }
349 <                                                /* copy back filled entries */
349 >                                                /* copy back blurred radii */
350          for (i = 0; i < GRIDRES; i++)
351              for (j = 0; j < GRIDRES; j++)
352 <                if (fill_grid[i][j])
353 <                        bsdf_grid[i][j].crad = fill_grid[i][j];
352 >                if (fill_cnt[i][j])
353 >                        dsf_grid[i][j].crad = fill_grid[i][j]/fill_cnt[i][j];
354   }
355  
356 < /* Cull points for more uniform distribution */
356 > /* Cull points for more uniform distribution, leave all nval 0 or 1 */
357   static void
358   cull_values(void)
359   {
# Line 305 | Line 363 | cull_values(void)
363                                                  /* simple greedy algorithm */
364          for (i = 0; i < GRIDRES; i++)
365              for (j = 0; j < GRIDRES; j++) {
366 <                if (!bsdf_grid[i][j].nval)
366 >                if (!dsf_grid[i][j].nval)
367                          continue;
368 <                if (!bsdf_grid[i][j].crad)
368 >                if (!dsf_grid[i][j].crad)
369                          continue;               /* shouldn't happen */
370                  vec_from_pos(ovec0, i, j);
371 <                maxang = 2.*R2ANG(bsdf_grid[i][j].crad);
371 >                maxang = 2.*R2ANG(dsf_grid[i][j].crad);
372                  if (maxang > ovec0[2])          /* clamp near horizon */
373                          maxang = ovec0[2];
374                  r = maxang*(2.*GRIDRES/M_PI) + 1;
# Line 321 | Line 379 | cull_values(void)
379                      for (jj = j-r; jj <= j+r; jj++) {
380                          if (jj < 0) continue;
381                          if (jj >= GRIDRES) break;
382 <                        if (!bsdf_grid[ii][jj].nval)
382 >                        if (!dsf_grid[ii][jj].nval)
383                                  continue;
384                          if ((ii == i) & (jj == j))
385                                  continue;       /* don't get self-absorbed */
# Line 329 | Line 387 | cull_values(void)
387                          if (2. - 2.*DOT(ovec0,ovec1) >= maxang2)
388                                  continue;
389                                                  /* absorb sum */
390 <                        bsdf_grid[i][j].vsum += bsdf_grid[ii][jj].vsum;
391 <                        bsdf_grid[i][j].nval += bsdf_grid[ii][jj].nval;
390 >                        dsf_grid[i][j].vsum += dsf_grid[ii][jj].vsum;
391 >                        dsf_grid[i][j].nval += dsf_grid[ii][jj].nval;
392                                                  /* keep value, though */
393 <                        bsdf_grid[ii][jj].vsum /= (double)bsdf_grid[ii][jj].nval;
394 <                        bsdf_grid[ii][jj].nval = 0;
393 >                        dsf_grid[ii][jj].vsum /= (float)dsf_grid[ii][jj].nval;
394 >                        dsf_grid[ii][jj].nval = 0;
395                      }
396                  }
397              }
398 +                                                /* final averaging pass */
399 +        for (i = 0; i < GRIDRES; i++)
400 +            for (j = 0; j < GRIDRES; j++)
401 +                if (dsf_grid[i][j].nval > 1) {
402 +                        dsf_grid[i][j].vsum /= (float)dsf_grid[i][j].nval;
403 +                        dsf_grid[i][j].nval = 1;
404 +                }
405   }
406  
407 + /* Compute (and allocate) migration price matrix for optimization */
408 + static float *
409 + price_routes(const RBFLIST *from_rbf, const RBFLIST *to_rbf)
410 + {
411 +        float   *pmtx = (float *)malloc(sizeof(float) *
412 +                                        from_rbf->nrbf * to_rbf->nrbf);
413 +        FVECT   *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
414 +        int     i, j;
415  
416 +        if ((pmtx == NULL) | (vto == NULL)) {
417 +                fputs("Out of memory in migration_costs()\n", stderr);
418 +                exit(1);
419 +        }
420 +        for (j = to_rbf->nrbf; j--; )           /* save repetitive ops. */
421 +                vec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);
422 +
423 +        for (i = from_rbf->nrbf; i--; ) {
424 +            const double        from_ang = R2ANG(from_rbf->rbfa[i].crad);
425 +            FVECT               vfrom;
426 +            vec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
427 +            for (j = to_rbf->nrbf; j--; )
428 +                pmtx[i*to_rbf->nrbf + j] = acos(DOT(vfrom, vto[j])) +
429 +                                fabs(R2ANG(to_rbf->rbfa[j].crad) - from_ang);
430 +        }
431 +        free(vto);
432 +        return(pmtx);
433 + }
434 +
435 + /* Comparison routine needed for sorting price row */
436 + static const float      *price_arr;
437 + static int
438 + msrt_cmp(const void *p1, const void *p2)
439 + {
440 +        float   c1 = price_arr[*(const int *)p1];
441 +        float   c2 = price_arr[*(const int *)p2];
442 +
443 +        if (c1 > c2) return(1);
444 +        if (c1 < c2) return(-1);
445 +        return(0);
446 + }
447 +
448 + /* Compute minimum (optimistic) cost for moving the given source material */
449 + static double
450 + min_cost(double amt2move, const double *avail, const float *price, int n)
451 + {
452 +        static int      *price_sort = NULL;
453 +        static int      n_alloc = 0;
454 +        double          total_cost = 0;
455 +        int             i;
456 +
457 +        if (amt2move <= FTINY)                  /* pre-emptive check */
458 +                return(0.);
459 +        if (n > n_alloc) {                      /* (re)allocate sort array */
460 +                if (n_alloc) free(price_sort);
461 +                price_sort = (int *)malloc(sizeof(int)*n);
462 +                if (price_sort == NULL) {
463 +                        fputs("Out of memory in min_cost()\n", stderr);
464 +                        exit(1);
465 +                }
466 +                n_alloc = n;
467 +        }
468 +        for (i = n; i--; )
469 +                price_sort[i] = i;
470 +        price_arr = price;
471 +        qsort(price_sort, n, sizeof(int), &msrt_cmp);
472 +                                                /* move cheapest first */
473 +        for (i = 0; i < n && amt2move > FTINY; i++) {
474 +                int     d = price_sort[i];
475 +                double  amt = (amt2move < avail[d]) ? amt2move : avail[d];
476 +
477 +                total_cost += amt * price[d];
478 +                amt2move -= amt;
479 +        }
480 + if (amt2move > 1e-5) fprintf(stderr, "%g leftover!\n", amt2move);
481 +        return(total_cost);
482 + }
483 +
484 + /* Take a step in migration by choosing optimal bucket to transfer */
485 + static double
486 + migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, const float *pmtx)
487 + {
488 +        static double   *src_cost = NULL;
489 +        int             n_alloc = 0;
490 +        const double    maxamt = 2./(mtx_nrows(mig)*mtx_ncols(mig));
491 +        double          amt = 0;
492 +        struct {
493 +                int     s, d;   /* source and destination */
494 +                double  price;  /* price estimate per amount moved */
495 +                double  amt;    /* amount we can move */
496 +        } cur, best;
497 +        int             i;
498 +
499 +        if (mtx_nrows(mig) > n_alloc) {         /* allocate cost array */
500 +                if (n_alloc)
501 +                        free(src_cost);
502 +                src_cost = (double *)malloc(sizeof(double)*mtx_nrows(mig));
503 +                if (src_cost == NULL) {
504 +                        fputs("Out of memory in migration_step()\n", stderr);
505 +                        exit(1);
506 +                }
507 +                n_alloc = mtx_nrows(mig);
508 +        }
509 +        for (i = mtx_nrows(mig); i--; )         /* starting costs for diff. */
510 +                src_cost[i] = min_cost(src_rem[i], dst_rem,
511 +                                        pmtx+i*mtx_ncols(mig), mtx_ncols(mig));
512 +
513 +                                                /* find best source & dest. */
514 +        best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
515 +        for (cur.s = mtx_nrows(mig); cur.s--; ) {
516 +            const float *price = pmtx + cur.s*mtx_ncols(mig);
517 +            double      cost_others = 0;
518 +            if (src_rem[cur.s] <= FTINY)
519 +                    continue;
520 +            cur.d = -1;                         /* examine cheapest dest. */
521 +            for (i = mtx_ncols(mig); i--; )
522 +                if (dst_rem[i] > FTINY &&
523 +                                (cur.d < 0 || price[i] < price[cur.d]))
524 +                        cur.d = i;
525 +            if (cur.d < 0)
526 +                    return(.0);
527 +            if ((cur.price = price[cur.d]) >= best.price)
528 +                    continue;                   /* no point checking further */
529 +            cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
530 +                                src_rem[cur.s] : dst_rem[cur.d];
531 +            if (cur.amt > maxamt) cur.amt = maxamt;
532 +            dst_rem[cur.d] -= cur.amt;          /* add up differential costs */
533 +            for (i = mtx_nrows(mig); i--; ) {
534 +                if (i == cur.s) continue;
535 +                cost_others += min_cost(src_rem[i], dst_rem, price, mtx_ncols(mig))
536 +                                        - src_cost[i];
537 +            }
538 +            dst_rem[cur.d] += cur.amt;          /* undo trial move */
539 +            cur.price += cost_others/cur.amt;   /* adjust effective price */
540 +            if (cur.price < best.price)         /* are we better than best? */
541 +                    best = cur;
542 +        }
543 +        if ((best.s < 0) | (best.d < 0))
544 +                return(.0);
545 +                                                /* make the actual move */
546 +        mig->mtx[mtx_ndx(mig,best.s,best.d)] += best.amt;
547 +        src_rem[best.s] -= best.amt;
548 +        dst_rem[best.d] -= best.amt;
549 +        return(best.amt);
550 + }
551 +
552 + /* Compute (and insert) migration along directed edge */
553 + static MIGRATION *
554 + make_migration(RBFLIST *from_rbf, RBFLIST *to_rbf)
555 + {
556 +        const double    end_thresh = 0.02/(from_rbf->nrbf*to_rbf->nrbf);
557 +        float           *pmtx = price_routes(from_rbf, to_rbf);
558 +        MIGRATION       *newmig = (MIGRATION *)malloc(sizeof(MIGRATION) +
559 +                                                        sizeof(float) *
560 +                                        (from_rbf->nrbf*to_rbf->nrbf - 1));
561 +        double          *src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
562 +        double          *dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
563 +        double          total_rem = 1.;
564 +        int             i;
565 +
566 +        if ((newmig == NULL) | (src_rem == NULL) | (dst_rem == NULL)) {
567 +                fputs("Out of memory in make_migration()\n", stderr);
568 +                exit(1);
569 +        }
570 +        newmig->next = NULL;
571 +        newmig->rbfv[0] = from_rbf;
572 +        newmig->rbfv[1] = to_rbf;
573 +        newmig->enxt[0] = newmig->enxt[1] = NULL;
574 +        memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
575 +                                                /* starting quantities */
576 +        for (i = from_rbf->nrbf; i--; )
577 +                src_rem[i] = rbf_volume(&from_rbf->rbfa[i]) / from_rbf->vtotal;
578 +        for (i = to_rbf->nrbf; i--; )
579 +                dst_rem[i] = rbf_volume(&to_rbf->rbfa[i]) / to_rbf->vtotal;
580 +                                                /* move a bit at a time */
581 +        while (total_rem > end_thresh)
582 +                total_rem -= migration_step(newmig, src_rem, dst_rem, pmtx);
583 +
584 +        free(pmtx);                             /* free working arrays */
585 +        free(src_rem);
586 +        free(dst_rem);
587 +        for (i = from_rbf->nrbf; i--; ) {       /* normalize final matrix */
588 +            float       nf = rbf_volume(&from_rbf->rbfa[i]);
589 +            int         j;
590 +            if (nf <= FTINY) continue;
591 +            nf = from_rbf->vtotal / nf;
592 +            for (j = to_rbf->nrbf; j--; )
593 +                newmig->mtx[mtx_ndx(newmig,i,j)] *= nf;
594 +        }
595 +                                                /* insert in edge lists */
596 +        newmig->enxt[0] = from_rbf->ejl;
597 +        from_rbf->ejl = newmig;
598 +        newmig->enxt[1] = to_rbf->ejl;
599 +        to_rbf->ejl = newmig;
600 +        newmig->next = mig_list;
601 +        return(mig_list = newmig);
602 + }
603 +
604 + #if 0
605 + /* Partially advect between the given RBFs to a newly allocated one */
606 + static RBFLIST *
607 + advect_rbf(const RBFLIST *from_rbf, const RBFLIST *to_rbf,
608 +                        const float *mtx, const FVECT invec)
609 + {
610 +        RBFLIST         *rbf;
611 +
612 +        if (from_rbf->nrbf > to_rbf->nrbf) {
613 +                fputs("Internal error: source RBF won't fit destination\n",
614 +                                stderr);
615 +                exit(1);
616 +        }
617 +        rbf = (RBFLIST *)malloc(sizeof(RBFLIST) + sizeof(RBFVAL)*(to_rbf->nrbf-1));
618 +        if (rbf == NULL) {
619 +                fputs("Out of memory in advect_rbf()\n", stderr);
620 +                exit(1);
621 +        }
622 +        rbf->next = NULL; rbf->ejl = NULL;
623 +        VCOPY(rbf->invec, invec);
624 +        rbf->vtotal = 0;
625 +        rbf->nrbf = to_rbf->nrbf;
626 +        
627 +        return rbf;
628 + }
629 + #endif
630 +
631   #if 1
632   /* Test main produces a Radiance model from the given input file */
633   int
# Line 357 | Line 645 | main(int argc, char *argv[])
645          }
646          if (!load_bsdf_meas(argv[1]))
647                  return(1);
360 –                                                /* produce spheres at meas. */
361 –        puts("void plastic orange\n0\n0\n5 .6 .4 .01 .04 .08\n");
362 –        n = 0;
363 –        for (i = 0; i < GRIDRES; i++)
364 –            for (j = 0; j < GRIDRES; j++)
365 –                if (bsdf_grid[i][j].nval) {
366 –                        double  bsdf = bsdf_grid[i][j].vsum /
367 –                                        (double)bsdf_grid[i][j].nval;
368 –                        FVECT   dir;
648  
370 –                        vec_from_pos(dir, i, j);
371 –                        printf("orange sphere s%04d\n0\n0\n", ++n);
372 –                        printf("4 %.6g %.6g %.6g .0015\n\n",
373 –                                        dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
374 –                }
649          compute_radii();
650          cull_values();
651 <                                                /* highlight chosen values */
651 >        make_rbfrep();
652 >                                                /* produce spheres at meas. */
653 >        puts("void plastic yellow\n0\n0\n5 .6 .4 .01 .04 .08\n");
654          puts("void plastic pink\n0\n0\n5 .5 .05 .9 .04 .08\n");
655          n = 0;
656          for (i = 0; i < GRIDRES; i++)
657              for (j = 0; j < GRIDRES; j++)
658 <                if (bsdf_grid[i][j].nval) {
383 <                        bsdf = bsdf_grid[i][j].vsum /
384 <                                        (double)bsdf_grid[i][j].nval;
658 >                if (dsf_grid[i][j].vsum > .0f) {
659                          vec_from_pos(dir, i, j);
660 <                        printf("pink cone c%04d\n0\n0\n8\n", ++n);
661 <                        printf("\t%.6g %.6g %.6g\n",
660 >                        bsdf = dsf_grid[i][j].vsum / dir[2];
661 >                        if (dsf_grid[i][j].nval) {
662 >                                printf("pink cone c%04d\n0\n0\n8\n", ++n);
663 >                                printf("\t%.6g %.6g %.6g\n",
664                                          dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
665 <                        printf("\t%.6g %.6g %.6g\n",
665 >                                printf("\t%.6g %.6g %.6g\n",
666                                          dir[0]*(bsdf+.005), dir[1]*(bsdf+.005),
667                                          dir[2]*(bsdf+.005));
668 <                        puts("\t.003\t0\n");
668 >                                puts("\t.003\t0\n");
669 >                        } else {
670 >                                vec_from_pos(dir, i, j);
671 >                                printf("yellow sphere s%04d\n0\n0\n", ++n);
672 >                                printf("4 %.6g %.6g %.6g .0015\n\n",
673 >                                        dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
674 >                        }
675                  }
676                                                  /* output continuous surface */
395 –        make_rbfrep();
677          puts("void trans tgreen\n0\n0\n7 .7 1 .7 .04 .04 .9 .9\n");
678          fflush(stdout);
679 <        sprintf(buf, "gensurf tgreen bsdf - - - %d %d", GRIDRES, GRIDRES);
679 >        sprintf(buf, "gensurf tgreen bsdf - - - %d %d", GRIDRES-1, GRIDRES-1);
680          pfp = popen(buf, "w");
681          if (pfp == NULL) {
682                  fputs(buf, stderr);
# Line 405 | Line 686 | main(int argc, char *argv[])
686          for (i = 0; i < GRIDRES; i++)
687              for (j = 0; j < GRIDRES; j++) {
688                  vec_from_pos(dir, i, j);
689 <                bsdf = eval_rbfrep(bsdf_list, dir);
689 >                bsdf = eval_rbfrep(dsf_list, dir) / dir[2];
690                  fprintf(pfp, "%.8e %.8e %.8e\n",
691                                  dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
692              }

Diff Legend

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