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

Comparing ray/src/cv/bsdfmesh.c (file contents):
Revision 2.8 by greg, Wed Mar 20 01:00:22 2013 UTC vs.
Revision 2.23 by greg, Wed Mar 12 00:39:43 2014 UTC

# Line 18 | Line 18 | static const char RCSid[] = "$Id$";
18   #include <string.h>
19   #include <math.h>
20   #include "bsdfrep.h"
21 +
22 + #ifndef NEIGH_FACT2
23 + #define NEIGH_FACT2     0.1     /* empirical neighborhood distance weight */
24 + #endif
25                                  /* number of processes to run */
26   int                     nprocs = 1;
27                                  /* number of children (-1 in child) */
# Line 27 | Line 31 | typedef struct {
31          int             nrows, ncols;   /* array size (matches migration) */
32          float           *price;         /* migration prices */
33          short           *sord;          /* sort for each row, low to high */
34 +        float           *prow;          /* current price row */
35   } PRICEMAT;                     /* sorted pricing matrix */
36  
37   #define pricerow(p,i)   ((p)->price + (i)*(p)->ncols)
# Line 134 | Line 139 | run_subprocess(void)
139  
140   #endif  /* ! _WIN32 */
141  
142 + /* Compute normalized distribution scattering functions for comparison */
143 + static void
144 + compute_nDSFs(const RBFNODE *rbf0, const RBFNODE *rbf1)
145 + {
146 +        const double    nf0 = (GRIDRES*GRIDRES) / rbf0->vtotal;
147 +        const double    nf1 = (GRIDRES*GRIDRES) / rbf1->vtotal;
148 +        int             x, y;
149 +        FVECT           dv;
150 +
151 +        for (x = GRIDRES; x--; )
152 +            for (y = GRIDRES; y--; ) {
153 +                ovec_from_pos(dv, x, y);        /* cube root (brightness) */
154 +                dsf_grid[x][y].val[0] = pow(nf0*eval_rbfrep(rbf0, dv), .3333);
155 +                dsf_grid[x][y].val[1] = pow(nf1*eval_rbfrep(rbf1, dv), .3333);
156 +            }
157 + }      
158 +
159 + /* Compute neighborhood distance-squared (dissimilarity) */
160 + static double
161 + neighborhood_dist2(int x0, int y0, int x1, int y1)
162 + {
163 +        int     rad = GRIDRES>>5;
164 +        double  sum2 = 0.;
165 +        double  d;
166 +        int     p[4];
167 +        int     i, j;
168 +
169 +        if ((x0 == x1) & (y0 == y1))
170 +                return(0.);
171 +                                                /* check radius */
172 +        p[0] = x0; p[1] = y0; p[2] = x1; p[3] = y1;
173 +        for (i = 4; i--; ) {
174 +                if (p[i] < rad) rad = p[i];
175 +                if (GRIDRES-1-p[i] < rad) rad = GRIDRES-1-p[i];
176 +        }
177 +        for (i = -rad; i <= rad; i++)
178 +            for (j = -rad; j <= rad; j++) {
179 +                d = dsf_grid[x0+i][y0+j].val[0] -
180 +                        dsf_grid[x1+i][y1+j].val[1];
181 +                sum2 += d*d;
182 +            }
183 +        return(sum2 / (4*rad*(rad+1) + 1));
184 + }
185 +
186   /* Comparison routine needed for sorting price row */
187   static int
188   msrt_cmp(void *b, const void *p1, const void *p2)
189   {
190          PRICEMAT        *pm = (PRICEMAT *)b;
191 <        int             ri = ((const short *)p1 - pm->sord) / pm->ncols;
192 <        float           c1 = pricerow(pm,ri)[*(const short *)p1];
144 <        float           c2 = pricerow(pm,ri)[*(const short *)p2];
191 >        float           c1 = pm->prow[*(const short *)p1];
192 >        float           c2 = pm->prow[*(const short *)p2];
193  
194          if (c1 > c2) return(1);
195          if (c1 < c2) return(-1);
# Line 155 | Line 203 | price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, co
203          FVECT   *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
204          int     i, j;
205  
206 +        compute_nDSFs(from_rbf, to_rbf);
207          pm->nrows = from_rbf->nrbf;
208          pm->ncols = to_rbf->nrbf;
209          pm->price = (float *)malloc(sizeof(float) * pm->nrows*pm->ncols);
# Line 171 | Line 220 | price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, co
220          for (i = from_rbf->nrbf; i--; ) {
221              const double        from_ang = R2ANG(from_rbf->rbfa[i].crad);
222              FVECT               vfrom;
223 +            short               *srow;
224              ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
225 +            pm->prow = pricerow(pm,i);
226 +            srow = psortrow(pm,i);
227              for (j = to_rbf->nrbf; j--; ) {
228 <                double          dprod = DOT(vfrom, vto[j]);
229 <                pricerow(pm,i)[j] = ((dprod >= 1.) ? .0 : acos(dprod)) +
230 <                                fabs(R2ANG(to_rbf->rbfa[j].crad) - from_ang);
231 <                psortrow(pm,i)[j] = j;
228 >                double  d;                      /* quadratic cost function */
229 >                d = Acos(DOT(vfrom, vto[j]));
230 >                pm->prow[j] = d*d;
231 >                d = R2ANG(to_rbf->rbfa[j].crad) - from_ang;
232 >                pm->prow[j] += d*d;
233 >                                                /* neighborhood difference */
234 >                pm->prow[j] += NEIGH_FACT2 * neighborhood_dist2(
235 >                                from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy,
236 >                                to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy );
237 >                srow[j] = j;
238              }
239 <            qsort_r(psortrow(pm,i), pm->ncols, sizeof(short), pm, &msrt_cmp);
239 >            qsort_r(srow, pm->ncols, sizeof(short), pm, &msrt_cmp);
240          }
241          free(vto);
242   }
# Line 195 | Line 253 | free_routes(PRICEMAT *pm)
253   static double
254   min_cost(double amt2move, const double *avail, const PRICEMAT *pm, int s)
255   {
256 +        const short     *srow = psortrow(pm,s);
257 +        const float     *prow = pricerow(pm,s);
258          double          total_cost = 0;
259          int             j;
200
201        if (amt2move <= FTINY)                  /* pre-emptive check */
202                return(.0);
260                                                  /* move cheapest first */
261 <        for (j = 0; j < pm->ncols && amt2move > FTINY; j++) {
262 <                int     d = psortrow(pm,s)[j];
261 >        for (j = 0; (j < pm->ncols) & (amt2move > FTINY); j++) {
262 >                int     d = srow[j];
263                  double  amt = (amt2move < avail[d]) ? amt2move : avail[d];
264  
265 <                total_cost += amt * pricerow(pm,s)[d];
265 >                total_cost += amt * prow[d];
266                  amt2move -= amt;
267          }
268          return(total_cost);
269   }
270  
271 < /* Take a step in migration by choosing optimal bucket to transfer */
271 > typedef struct {
272 >        short   s, d;           /* source and destination */
273 >        float   dc;             /* discount to push inventory */
274 > } ROWSENT;              /* row sort entry */
275 >
276 > /* Compare entries by discounted moving price */
277 > static int
278 > rmovcmp(void *b, const void *p1, const void *p2)
279 > {
280 >        PRICEMAT        *pm = (PRICEMAT *)b;
281 >        const ROWSENT   *re1 = (const ROWSENT *)p1;
282 >        const ROWSENT   *re2 = (const ROWSENT *)p2;
283 >        double          price_diff;
284 >
285 >        if (re1->d < 0) return(re2->d >= 0);
286 >        if (re2->d < 0) return(-1);
287 >        price_diff = re1->dc*pricerow(pm,re1->s)[re1->d] -
288 >                        re2->dc*pricerow(pm,re2->s)[re2->d];
289 >        if (price_diff > 0) return(1);
290 >        if (price_diff < 0) return(-1);
291 >        return(0);
292 > }
293 >
294 > /* Take a step in migration by choosing reasonable bucket to transfer */
295   static double
296 < migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, const PRICEMAT *pm)
296 > migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, PRICEMAT *pm)
297   {
298 +        const int       max2check = 100;
299          const double    maxamt = 1./(double)pm->ncols;
300 <        const double    minamt = maxamt*5e-6;
300 >        const double    minamt = maxamt*1e-4;
301          double          *src_cost;
302 +        ROWSENT         *rord;
303          struct {
304                  int     s, d;   /* source and destination */
305 <                double  price;  /* price estimate per amount moved */
305 >                double  price;  /* cost per amount moved */
306                  double  amt;    /* amount we can move */
307          } cur, best;
308 <        int             i;
308 >        int             r2check, i, ri;
309 >        /*
310 >         * Check cheapest available routes only -- a higher adjusted
311 >         * destination price implies that another source is closer, so
312 >         * we can hold off considering more expensive options until
313 >         * some other (hopefully better) moves have been made.
314 >         * A discount based on source remaining is supposed to prioritize
315 >         * movement from large lobes, but it doesn't seem to do much,
316 >         * so we have it set to 1.0 at the moment.
317 >         */
318 > #define discount(qr)    1.0
319 >                                                /* most promising row order */
320 >        rord = (ROWSENT *)malloc(sizeof(ROWSENT)*pm->nrows);
321 >        if (rord == NULL)
322 >                goto memerr;
323 >        for (ri = pm->nrows; ri--; ) {
324 >            rord[ri].s = ri;
325 >            rord[ri].d = -1;
326 >            rord[ri].dc = 1.f;
327 >            if (src_rem[ri] <= minamt)          /* enough source material? */
328 >                    continue;
329 >            for (i = 0; i < pm->ncols; i++)
330 >                if (dst_rem[ rord[ri].d = psortrow(pm,ri)[i] ] > minamt)
331 >                        break;
332 >            if (i >= pm->ncols) {               /* moved all we can? */
333 >                free(rord);
334 >                return(.0);
335 >            }
336 >            rord[ri].dc = discount(src_rem[ri]);
337 >        }
338 >        if (pm->nrows > max2check)              /* sort if too many sources */
339 >                qsort_r(rord, pm->nrows, sizeof(ROWSENT), pm, &rmovcmp);
340                                                  /* allocate cost array */
341          src_cost = (double *)malloc(sizeof(double)*pm->nrows);
342 <        if (src_cost == NULL) {
343 <                fprintf(stderr, "%s: Out of memory in migration_step()\n",
231 <                                progname);
232 <                exit(1);
233 <        }
342 >        if (src_cost == NULL)
343 >                goto memerr;
344          for (i = pm->nrows; i--; )              /* starting costs for diff. */
345                  src_cost[i] = min_cost(src_rem[i], dst_rem, pm, i);
236
346                                                  /* find best source & dest. */
347          best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
348 <        for (cur.s = pm->nrows; cur.s--; ) {
348 >        if ((r2check = pm->nrows) > max2check)
349 >                r2check = max2check;            /* put a limit on search */
350 >        for (ri = 0; ri < r2check; ri++) {      /* check each source row */
351              double      cost_others = 0;
352 <
353 <            if (src_rem[cur.s] <= minamt)
354 <                    continue;
355 <                                                /* examine cheapest dest. */
356 <            for (i = 0; i < pm->ncols; i++)
357 <                if (dst_rem[ cur.d = psortrow(pm,cur.s)[i] ] > minamt)
247 <                        break;
248 <            if (i >= pm->ncols)
249 <                break;
250 <            if ((cur.price = pricerow(pm,cur.s)[cur.d]) >= best.price)
251 <                continue;                       /* no point checking further */
352 >            cur.s = rord[ri].s;
353 >            if ((cur.d = rord[ri].d) < 0 ||
354 >                        rord[ri].dc*pricerow(pm,cur.s)[cur.d] >= best.price) {
355 >                if (pm->nrows > max2check) break;       /* sorted end */
356 >                continue;                       /* else skip this one */
357 >            }
358              cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
359                                  src_rem[cur.s] : dst_rem[cur.d];
360 <            if (cur.amt > maxamt) cur.amt = maxamt;
361 <            dst_rem[cur.d] -= cur.amt;          /* add up differential costs */
360 >                                                /* don't just leave smidgen */
361 >            if (cur.amt > maxamt*1.02) cur.amt = maxamt;
362 >            dst_rem[cur.d] -= cur.amt;          /* add up opportunity costs */
363              for (i = pm->nrows; i--; )
364                  if (i != cur.s)
365 <                        cost_others += min_cost(src_rem[i], dst_rem, pm, i)
365 >                    cost_others += min_cost(src_rem[i], dst_rem, pm, i)
366                                          - src_cost[i];
367              dst_rem[cur.d] += cur.amt;          /* undo trial move */
368 <            cur.price += cost_others/cur.amt;   /* adjust effective price */
368 >                                                /* discount effective price */
369 >            cur.price = ( pricerow(pm,cur.s)[cur.d] + cost_others/cur.amt ) *
370 >                                        rord[ri].dc;
371              if (cur.price < best.price)         /* are we better than best? */
372 <                    best = cur;
372 >                best = cur;
373          }
374 <        free(src_cost);                         /* finish up */
375 <
374 >        free(src_cost);                         /* clean up */
375 >        free(rord);
376          if ((best.s < 0) | (best.d < 0))        /* nothing left to move? */
377                  return(.0);
378                                                  /* else make the actual move */
# Line 271 | Line 380 | migration_step(MIGRATION *mig, double *src_rem, double
380          src_rem[best.s] -= best.amt;
381          dst_rem[best.d] -= best.amt;
382          return(best.amt);
383 + memerr:
384 +        fprintf(stderr, "%s: Out of memory in migration_step()\n", progname);
385 +        exit(1);
386 + #undef discount
387   }
388  
389   /* Compute and insert migration along directed edge (may fork child) */
# Line 290 | Line 403 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
403                          return(NULL);
404                                                  /* else allocate */
405   #ifdef DEBUG
406 <        fprintf(stderr, "Building path from (theta,phi) (%.0f,%.0f) ",
406 >        fprintf(stderr, "Building path from (theta,phi) (%.1f,%.1f) ",
407                          get_theta180(from_rbf->invec),
408                          get_phi360(from_rbf->invec));
409 <        fprintf(stderr, "to (%.0f,%.0f) with %d x %d matrix\n",
409 >        fprintf(stderr, "to (%.1f,%.1f) with %d x %d matrix\n",
410                          get_theta180(to_rbf->invec),
411                          get_phi360(to_rbf->invec),
412                          from_rbf->nrbf, to_rbf->nrbf);
# Line 361 | Line 474 | overlaps_tri(const RBFNODE *bv0, const RBFNODE *bv1, c
474          return(vother[im_rev] != NULL);
475   }
476  
477 < /* Find context hull vertex to complete triangle (oriented call) */
477 > /* Find convex hull vertex to complete triangle (oriented call) */
478   static RBFNODE *
479   find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rbf1)
480   {
# Line 382 | Line 495 | find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rb
495                  if (DOT(vp, vmid) <= FTINY)
496                          continue;               /* wrong orientation */
497                  area2 = .25*DOT(vp,vp);
498 <                VSUB(vp, rbf->invec, rbf0->invec);
498 >                VSUB(vp, rbf->invec, vmid);
499                  dprod = -DOT(vp, vejn);
500                  VSUM(vp, vp, vejn, dprod);      /* above guarantees non-zero */
501                  dprod = DOT(vp, vmid) / VLEN(vp);
# Line 439 | Line 552 | mesh_from_edge(MIGRATION *edge)
552          }
553   }
554  
555 < /* Compute minimum BSDF from histogram and clear it */
555 > /* Add normal direction if missing */
556   static void
557 < comp_bsdf_min()
557 > check_normal_incidence(void)
558   {
559 <        int     cnt;
560 <        int     i, target;
559 >        static const FVECT      norm_vec = {.0, .0, 1.};
560 >        const int               saved_nprocs = nprocs;
561 >        RBFNODE                 *near_rbf, *mir_rbf, *rbf;
562 >        double                  bestd;
563 >        int                     n;
564  
449        cnt = 0;
450        for (i = HISTLEN; i--; )
451                cnt += bsdf_hist[i];
565  
566 <        target = cnt/100;                       /* ignore bottom 1% */
567 <        cnt = 0;
568 <        for (i = 0; cnt <= target; i++)
569 <                cnt += bsdf_hist[i];
570 <        bsdf_min = histval(i-1);
571 <        memset(bsdf_hist, 0, sizeof(bsdf_hist));
566 >        if (dsf_list == NULL)
567 >                return;                         /* XXX should be error? */
568 >        near_rbf = dsf_list;
569 >        bestd = input_orient*near_rbf->invec[2];
570 >        if (single_plane_incident) {            /* ordered plane incidence? */
571 >                if (bestd >= 1.-2.*FTINY)
572 >                        return;                 /* already have normal */
573 >        } else {
574 >                switch (inp_coverage) {
575 >                case INP_QUAD1:
576 >                case INP_QUAD2:
577 >                case INP_QUAD3:
578 >                case INP_QUAD4:
579 >                        break;                  /* quadrilateral symmetry? */
580 >                default:
581 >                        return;                 /* else we can interpolate */
582 >                }
583 >                for (rbf = near_rbf->next; rbf != NULL; rbf = rbf->next) {
584 >                        const double    d = input_orient*rbf->invec[2];
585 >                        if (d >= 1.-2.*FTINY)
586 >                                return;         /* seems we have normal */
587 >                        if (d > bestd) {
588 >                                near_rbf = rbf;
589 >                                bestd = d;
590 >                        }
591 >                }
592 >        }
593 >        if (mig_list != NULL) {                 /* need to be called first */
594 >                fprintf(stderr, "%s: Late call to check_normal_incidence()\n",
595 >                                progname);
596 >                exit(1);
597 >        }
598 > #ifdef DEBUG
599 >        fprintf(stderr, "Interpolating normal incidence by mirroring (%.1f,%.1f)\n",
600 >                        get_theta180(near_rbf->invec), get_phi360(near_rbf->invec));
601 > #endif
602 >                                                /* mirror nearest incidence */
603 >        n = sizeof(RBFNODE) + sizeof(RBFVAL)*(near_rbf->nrbf-1);
604 >        mir_rbf = (RBFNODE *)malloc(n);
605 >        if (mir_rbf == NULL)
606 >                goto memerr;
607 >        memcpy(mir_rbf, near_rbf, n);
608 >        mir_rbf->ord = near_rbf->ord - 1;       /* not used, I think */
609 >        mir_rbf->next = NULL;
610 >        mir_rbf->ejl = NULL;
611 >        rev_rbf_symmetry(mir_rbf, MIRROR_X|MIRROR_Y);
612 >        nprocs = 1;                             /* compute migration matrix */
613 >        if (create_migration(mir_rbf, near_rbf) == NULL)
614 >                exit(1);                        /* XXX should never happen! */
615 >                                                /* interpolate normal dist. */
616 >        rbf = e_advect_rbf(mig_list, norm_vec, 2*near_rbf->nrbf);
617 >        nprocs = saved_nprocs;                  /* final clean-up */
618 >        free(mir_rbf);
619 >        free(mig_list);
620 >        mig_list = near_rbf->ejl = NULL;
621 >        insert_dsf(rbf);                        /* insert interpolated normal */
622 >        return;
623 > memerr:
624 >        fprintf(stderr, "%s: Out of memory in check_normal_incidence()\n",
625 >                                progname);
626 >        exit(1);
627   }
628          
629   /* Build our triangle mesh from recorded RBFs */
# Line 465 | Line 633 | build_mesh(void)
633          double          best2 = M_PI*M_PI;
634          RBFNODE         *shrt_edj[2];
635          RBFNODE         *rbf0, *rbf1;
636 +                                                /* add normal if needed */
637 +        check_normal_incidence();
638                                                  /* check if isotropic */
639          if (single_plane_incident) {
640                  for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
# Line 492 | Line 662 | build_mesh(void)
662                  mesh_from_edge(create_migration(shrt_edj[0], shrt_edj[1]));
663          else
664                  mesh_from_edge(create_migration(shrt_edj[1], shrt_edj[0]));
495                                                /* compute minimum BSDF */
496        comp_bsdf_min();
665                                                  /* complete migrations */
666          await_children(nchild);
667   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines