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.11 by greg, Thu Oct 24 16:11:37 2013 UTC vs.
Revision 2.24 by greg, Sat Mar 15 18:11:37 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 135 | 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 +                                                /* check radius */
169 +        p[0] = x0; p[1] = y0; p[2] = x1; p[3] = y1;
170 +        for (i = 4; i--; ) {
171 +                if (p[i] < rad) rad = p[i];
172 +                if (GRIDRES-1-p[i] < rad) rad = GRIDRES-1-p[i];
173 +        }
174 +        for (i = -rad; i <= rad; i++)
175 +            for (j = -rad; j <= rad; j++) {
176 +                d = dsf_grid[x0+i][y0+j].val[0] -
177 +                        dsf_grid[x1+i][y1+j].val[1];
178 +                sum2 += d*d;
179 +            }
180 +        return(sum2 / (4*rad*(rad+1) + 1));
181 + }
182 +
183   /* Comparison routine needed for sorting price row */
184   static int
185   msrt_cmp(void *b, const void *p1, const void *p2)
# Line 155 | Line 200 | price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, co
200          FVECT   *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
201          int     i, j;
202  
203 +        compute_nDSFs(from_rbf, to_rbf);
204          pm->nrows = from_rbf->nrbf;
205          pm->ncols = to_rbf->nrbf;
206          pm->price = (float *)malloc(sizeof(float) * pm->nrows*pm->ncols);
# Line 176 | Line 222 | price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, co
222              pm->prow = pricerow(pm,i);
223              srow = psortrow(pm,i);
224              for (j = to_rbf->nrbf; j--; ) {
225 <                double          dprod = DOT(vfrom, vto[j]);
226 <                pm->prow[j] = ((dprod >= 1.) ? .0 : acos(dprod)) +
227 <                                fabs(R2ANG(to_rbf->rbfa[j].crad) - from_ang);
225 >                double  d;                      /* quadratic cost function */
226 >                d = Acos(DOT(vfrom, vto[j]));
227 >                pm->prow[j] = d*d;
228 >                d = R2ANG(to_rbf->rbfa[j].crad) - from_ang;
229 >                pm->prow[j] += d*d;
230 >                                                /* neighborhood difference */
231 >                pm->prow[j] += NEIGH_FACT2 * neighborhood_dist2(
232 >                                from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy,
233 >                                to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy );
234                  srow[j] = j;
235              }
236              qsort_r(srow, pm->ncols, sizeof(short), pm, &msrt_cmp);
# Line 213 | Line 265 | min_cost(double amt2move, const double *avail, const P
265          return(total_cost);
266   }
267  
268 < /* Compare entries by moving price */
268 > typedef struct {
269 >        short   s, d;           /* source and destination */
270 >        float   dc;             /* discount to push inventory */
271 > } ROWSENT;              /* row sort entry */
272 >
273 > /* Compare entries by discounted moving price */
274   static int
275   rmovcmp(void *b, const void *p1, const void *p2)
276   {
277          PRICEMAT        *pm = (PRICEMAT *)b;
278 <        const short     *ij1 = (const short *)p1;
279 <        const short     *ij2 = (const short *)p2;
280 <        float           price_diff;
278 >        const ROWSENT   *re1 = (const ROWSENT *)p1;
279 >        const ROWSENT   *re2 = (const ROWSENT *)p2;
280 >        double          price_diff;
281  
282 <        if (ij1[1] < 0) return(ij2[1] >= 0);
283 <        if (ij2[1] < 0) return(-1);
284 <        price_diff = pricerow(pm,ij1[0])[ij1[1]] - pricerow(pm,ij2[0])[ij2[1]];
282 >        if (re1->d < 0) return(re2->d >= 0);
283 >        if (re2->d < 0) return(-1);
284 >        price_diff = re1->dc*pricerow(pm,re1->s)[re1->d] -
285 >                        re2->dc*pricerow(pm,re2->s)[re2->d];
286          if (price_diff > 0) return(1);
287          if (price_diff < 0) return(-1);
288          return(0);
# Line 236 | Line 294 | migration_step(MIGRATION *mig, double *src_rem, double
294   {
295          const int       max2check = 100;
296          const double    maxamt = 1./(double)pm->ncols;
297 <        const double    minamt = maxamt*5e-6;
297 >        const double    minamt = maxamt*1e-4;
298          double          *src_cost;
299 <        short           (*rord)[2];
299 >        ROWSENT         *rord;
300          struct {
301                  int     s, d;   /* source and destination */
302 <                double  price;  /* price estimate per amount moved */
302 >                double  price;  /* cost per amount moved */
303                  double  amt;    /* amount we can move */
304          } cur, best;
305          int             r2check, i, ri;
# Line 250 | Line 308 | migration_step(MIGRATION *mig, double *src_rem, double
308           * destination price implies that another source is closer, so
309           * we can hold off considering more expensive options until
310           * some other (hopefully better) moves have been made.
311 +         * A discount based on source remaining is supposed to prioritize
312 +         * movement from large lobes, but it doesn't seem to do much,
313 +         * so we have it set to 1.0 at the moment.
314           */
315 + #define discount(qr)    1.0
316                                                  /* most promising row order */
317 <        rord = (short (*)[2])malloc(sizeof(short)*2*pm->nrows);
317 >        rord = (ROWSENT *)malloc(sizeof(ROWSENT)*pm->nrows);
318          if (rord == NULL)
319                  goto memerr;
320          for (ri = pm->nrows; ri--; ) {
321 <            rord[ri][0] = ri;
322 <            rord[ri][1] = -1;
321 >            rord[ri].s = ri;
322 >            rord[ri].d = -1;
323 >            rord[ri].dc = 1.f;
324              if (src_rem[ri] <= minamt)          /* enough source material? */
325                      continue;
326              for (i = 0; i < pm->ncols; i++)
327 <                if (dst_rem[ rord[ri][1] = psortrow(pm,ri)[i] ] > minamt)
327 >                if (dst_rem[ rord[ri].d = psortrow(pm,ri)[i] ] > minamt)
328                          break;
329              if (i >= pm->ncols) {               /* moved all we can? */
330                  free(rord);
331                  return(.0);
332              }
333 +            rord[ri].dc = discount(src_rem[ri]);
334          }
335          if (pm->nrows > max2check)              /* sort if too many sources */
336 <                qsort_r(rord, pm->nrows, sizeof(short)*2, pm, &rmovcmp);
336 >                qsort_r(rord, pm->nrows, sizeof(ROWSENT), pm, &rmovcmp);
337                                                  /* allocate cost array */
338          src_cost = (double *)malloc(sizeof(double)*pm->nrows);
339          if (src_cost == NULL)
# Line 282 | Line 346 | migration_step(MIGRATION *mig, double *src_rem, double
346                  r2check = max2check;            /* put a limit on search */
347          for (ri = 0; ri < r2check; ri++) {      /* check each source row */
348              double      cost_others = 0;
349 <            cur.s = rord[ri][0];
350 <            if ((cur.d = rord[ri][1]) < 0 ||
351 <                        (cur.price = pricerow(pm,cur.s)[cur.d]) >= best.price) {
349 >            cur.s = rord[ri].s;
350 >            if ((cur.d = rord[ri].d) < 0 ||
351 >                        rord[ri].dc*pricerow(pm,cur.s)[cur.d] >= best.price) {
352                  if (pm->nrows > max2check) break;       /* sorted end */
353                  continue;                       /* else skip this one */
354              }
# Line 298 | Line 362 | migration_step(MIGRATION *mig, double *src_rem, double
362                      cost_others += min_cost(src_rem[i], dst_rem, pm, i)
363                                          - src_cost[i];
364              dst_rem[cur.d] += cur.amt;          /* undo trial move */
365 <            cur.price += cost_others/cur.amt;   /* adjust effective price */
365 >                                                /* discount effective price */
366 >            cur.price = ( pricerow(pm,cur.s)[cur.d] + cost_others/cur.amt ) *
367 >                                        rord[ri].dc;
368              if (cur.price < best.price)         /* are we better than best? */
369                  best = cur;
370          }
# Line 314 | Line 380 | migration_step(MIGRATION *mig, double *src_rem, double
380   memerr:
381          fprintf(stderr, "%s: Out of memory in migration_step()\n", progname);
382          exit(1);
383 + #undef discount
384   }
385  
386   /* Compute and insert migration along directed edge (may fork child) */
# Line 333 | Line 400 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
400                          return(NULL);
401                                                  /* else allocate */
402   #ifdef DEBUG
403 <        fprintf(stderr, "Building path from (theta,phi) (%.0f,%.0f) ",
403 >        fprintf(stderr, "Building path from (theta,phi) (%.1f,%.1f) ",
404                          get_theta180(from_rbf->invec),
405                          get_phi360(from_rbf->invec));
406 <        fprintf(stderr, "to (%.0f,%.0f) with %d x %d matrix\n",
406 >        fprintf(stderr, "to (%.1f,%.1f) with %d x %d matrix\n",
407                          get_theta180(to_rbf->invec),
408                          get_phi360(to_rbf->invec),
409                          from_rbf->nrbf, to_rbf->nrbf);
# Line 404 | Line 471 | overlaps_tri(const RBFNODE *bv0, const RBFNODE *bv1, c
471          return(vother[im_rev] != NULL);
472   }
473  
474 < /* Find context hull vertex to complete triangle (oriented call) */
474 > /* Find convex hull vertex to complete triangle (oriented call) */
475   static RBFNODE *
476   find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rbf1)
477   {
# Line 425 | Line 492 | find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rb
492                  if (DOT(vp, vmid) <= FTINY)
493                          continue;               /* wrong orientation */
494                  area2 = .25*DOT(vp,vp);
495 <                VSUB(vp, rbf->invec, rbf0->invec);
495 >                VSUB(vp, rbf->invec, vmid);
496                  dprod = -DOT(vp, vejn);
497                  VSUM(vp, vp, vejn, dprod);      /* above guarantees non-zero */
498                  dprod = DOT(vp, vmid) / VLEN(vp);
# Line 481 | Line 548 | mesh_from_edge(MIGRATION *edge)
548                  }
549          }
550   }
551 +
552 + /* Add normal direction if missing */
553 + static void
554 + check_normal_incidence(void)
555 + {
556 +        static const FVECT      norm_vec = {.0, .0, 1.};
557 +        const int               saved_nprocs = nprocs;
558 +        RBFNODE                 *near_rbf, *mir_rbf, *rbf;
559 +        double                  bestd;
560 +        int                     n;
561 +
562 +
563 +        if (dsf_list == NULL)
564 +                return;                         /* XXX should be error? */
565 +        near_rbf = dsf_list;
566 +        bestd = input_orient*near_rbf->invec[2];
567 +        if (single_plane_incident) {            /* ordered plane incidence? */
568 +                if (bestd >= 1.-2.*FTINY)
569 +                        return;                 /* already have normal */
570 +        } else {
571 +                switch (inp_coverage) {
572 +                case INP_QUAD1:
573 +                case INP_QUAD2:
574 +                case INP_QUAD3:
575 +                case INP_QUAD4:
576 +                        break;                  /* quadrilateral symmetry? */
577 +                default:
578 +                        return;                 /* else we can interpolate */
579 +                }
580 +                for (rbf = near_rbf->next; rbf != NULL; rbf = rbf->next) {
581 +                        const double    d = input_orient*rbf->invec[2];
582 +                        if (d >= 1.-2.*FTINY)
583 +                                return;         /* seems we have normal */
584 +                        if (d > bestd) {
585 +                                near_rbf = rbf;
586 +                                bestd = d;
587 +                        }
588 +                }
589 +        }
590 +        if (mig_list != NULL) {                 /* need to be called first */
591 +                fprintf(stderr, "%s: Late call to check_normal_incidence()\n",
592 +                                progname);
593 +                exit(1);
594 +        }
595 + #ifdef DEBUG
596 +        fprintf(stderr, "Interpolating normal incidence by mirroring (%.1f,%.1f)\n",
597 +                        get_theta180(near_rbf->invec), get_phi360(near_rbf->invec));
598 + #endif
599 +                                                /* mirror nearest incidence */
600 +        n = sizeof(RBFNODE) + sizeof(RBFVAL)*(near_rbf->nrbf-1);
601 +        mir_rbf = (RBFNODE *)malloc(n);
602 +        if (mir_rbf == NULL)
603 +                goto memerr;
604 +        memcpy(mir_rbf, near_rbf, n);
605 +        mir_rbf->ord = near_rbf->ord - 1;       /* not used, I think */
606 +        mir_rbf->next = NULL;
607 +        mir_rbf->ejl = NULL;
608 +        rev_rbf_symmetry(mir_rbf, MIRROR_X|MIRROR_Y);
609 +        nprocs = 1;                             /* compute migration matrix */
610 +        if (create_migration(mir_rbf, near_rbf) == NULL)
611 +                exit(1);                        /* XXX should never happen! */
612 +                                                /* interpolate normal dist. */
613 +        rbf = e_advect_rbf(mig_list, norm_vec, 2*near_rbf->nrbf);
614 +        nprocs = saved_nprocs;                  /* final clean-up */
615 +        free(mir_rbf);
616 +        free(mig_list);
617 +        mig_list = near_rbf->ejl = NULL;
618 +        insert_dsf(rbf);                        /* insert interpolated normal */
619 +        return;
620 + memerr:
621 +        fprintf(stderr, "%s: Out of memory in check_normal_incidence()\n",
622 +                                progname);
623 +        exit(1);
624 + }
625          
626   /* Build our triangle mesh from recorded RBFs */
627   void
# Line 489 | Line 630 | build_mesh(void)
630          double          best2 = M_PI*M_PI;
631          RBFNODE         *shrt_edj[2];
632          RBFNODE         *rbf0, *rbf1;
633 +                                                /* add normal if needed */
634 +        check_normal_incidence();
635                                                  /* check if isotropic */
636          if (single_plane_incident) {
637                  for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)

Diff Legend

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