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.14 by greg, Fri Nov 8 03:42:13 2013 UTC vs.
Revision 2.18 by greg, Thu Mar 6 00:40:37 2014 UTC

# Line 177 | Line 177 | price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, co
177              srow = psortrow(pm,i);
178              for (j = to_rbf->nrbf; j--; ) {
179                  double          d;              /* quadratic cost function */
180 <                d = DOT(vfrom, vto[j]);
181 <                d = (d >= 1.) ? .0 : acos(d);
180 >                d = Acos(DOT(vfrom, vto[j]));
181                  pm->prow[j] = d*d;
182                  d = R2ANG(to_rbf->rbfa[j].crad) - from_ang;
183                  pm->prow[j] += d*d;    
# Line 216 | Line 215 | min_cost(double amt2move, const double *avail, const P
215          return(total_cost);
216   }
217  
218 < /* Compare entries by moving price */
218 > typedef struct {
219 >        short   s, d;           /* source and destination */
220 >        float   dc;             /* discount to push inventory */
221 > } ROWSENT;              /* row sort entry */
222 >
223 > /* Compare entries by discounted moving price */
224   static int
225   rmovcmp(void *b, const void *p1, const void *p2)
226   {
227          PRICEMAT        *pm = (PRICEMAT *)b;
228 <        const short     *ij1 = (const short *)p1;
229 <        const short     *ij2 = (const short *)p2;
230 <        float           price_diff;
228 >        const ROWSENT   *re1 = (const ROWSENT *)p1;
229 >        const ROWSENT   *re2 = (const ROWSENT *)p2;
230 >        double          price_diff;
231  
232 <        if (ij1[1] < 0) return(ij2[1] >= 0);
233 <        if (ij2[1] < 0) return(-1);
234 <        price_diff = pricerow(pm,ij1[0])[ij1[1]] - pricerow(pm,ij2[0])[ij2[1]];
232 >        if (re1->d < 0) return(re2->d >= 0);
233 >        if (re2->d < 0) return(-1);
234 >        price_diff = re1->dc*pricerow(pm,re1->s)[re1->d] -
235 >                        re2->dc*pricerow(pm,re2->s)[re2->d];
236          if (price_diff > 0) return(1);
237          if (price_diff < 0) return(-1);
238          return(0);
# Line 241 | Line 246 | migration_step(MIGRATION *mig, double *src_rem, double
246          const double    maxamt = 1./(double)pm->ncols;
247          const double    minamt = maxamt*1e-4;
248          double          *src_cost;
249 <        short           (*rord)[2];
249 >        ROWSENT         *rord;
250          struct {
251                  int     s, d;   /* source and destination */
252 <                double  price;  /* price estimate per amount moved */
252 >                double  price;  /* cost per amount moved */
253                  double  amt;    /* amount we can move */
254          } cur, best;
255          int             r2check, i, ri;
# Line 253 | Line 258 | migration_step(MIGRATION *mig, double *src_rem, double
258           * destination price implies that another source is closer, so
259           * we can hold off considering more expensive options until
260           * some other (hopefully better) moves have been made.
261 +         * A discount based on source remaining is supposed to prioritize
262 +         * movement from large lobes, but it doesn't seem to do much,
263 +         * so we have it set to 1.0 at the moment.
264           */
265 + #define discount(qr)    1.0
266                                                  /* most promising row order */
267 <        rord = (short (*)[2])malloc(sizeof(short)*2*pm->nrows);
267 >        rord = (ROWSENT *)malloc(sizeof(ROWSENT)*pm->nrows);
268          if (rord == NULL)
269                  goto memerr;
270          for (ri = pm->nrows; ri--; ) {
271 <            rord[ri][0] = ri;
272 <            rord[ri][1] = -1;
271 >            rord[ri].s = ri;
272 >            rord[ri].d = -1;
273 >            rord[ri].dc = 1.f;
274              if (src_rem[ri] <= minamt)          /* enough source material? */
275                      continue;
276              for (i = 0; i < pm->ncols; i++)
277 <                if (dst_rem[ rord[ri][1] = psortrow(pm,ri)[i] ] > minamt)
277 >                if (dst_rem[ rord[ri].d = psortrow(pm,ri)[i] ] > minamt)
278                          break;
279              if (i >= pm->ncols) {               /* moved all we can? */
280                  free(rord);
281                  return(.0);
282              }
283 +            rord[ri].dc = discount(src_rem[ri]);
284          }
285          if (pm->nrows > max2check)              /* sort if too many sources */
286 <                qsort_r(rord, pm->nrows, sizeof(short)*2, pm, &rmovcmp);
286 >                qsort_r(rord, pm->nrows, sizeof(ROWSENT), pm, &rmovcmp);
287                                                  /* allocate cost array */
288          src_cost = (double *)malloc(sizeof(double)*pm->nrows);
289          if (src_cost == NULL)
# Line 285 | Line 296 | migration_step(MIGRATION *mig, double *src_rem, double
296                  r2check = max2check;            /* put a limit on search */
297          for (ri = 0; ri < r2check; ri++) {      /* check each source row */
298              double      cost_others = 0;
299 <            cur.s = rord[ri][0];
300 <            if ((cur.d = rord[ri][1]) < 0 ||
301 <                        (cur.price = pricerow(pm,cur.s)[cur.d]) >= best.price) {
299 >            cur.s = rord[ri].s;
300 >            if ((cur.d = rord[ri].d) < 0 ||
301 >                        rord[ri].dc*pricerow(pm,cur.s)[cur.d] >= best.price) {
302                  if (pm->nrows > max2check) break;       /* sorted end */
303                  continue;                       /* else skip this one */
304              }
# Line 301 | Line 312 | migration_step(MIGRATION *mig, double *src_rem, double
312                      cost_others += min_cost(src_rem[i], dst_rem, pm, i)
313                                          - src_cost[i];
314              dst_rem[cur.d] += cur.amt;          /* undo trial move */
315 <            cur.price += cost_others/cur.amt;   /* adjust effective price */
315 >                                                /* discount effective price */
316 >            cur.price = ( pricerow(pm,cur.s)[cur.d] + cost_others/cur.amt ) *
317 >                                        rord[ri].dc;
318              if (cur.price < best.price)         /* are we better than best? */
319                  best = cur;
320          }
# Line 317 | Line 330 | migration_step(MIGRATION *mig, double *src_rem, double
330   memerr:
331          fprintf(stderr, "%s: Out of memory in migration_step()\n", progname);
332          exit(1);
333 + #undef discount
334   }
335  
336   /* Compute and insert migration along directed edge (may fork child) */
# Line 484 | Line 498 | mesh_from_edge(MIGRATION *edge)
498                  }
499          }
500   }
501 +
502 + /* Add normal direction if missing */
503 + static void
504 + check_normal_incidence(void)
505 + {
506 +        static const FVECT      norm_vec = {.0, .0, 1.};
507 +        const int               saved_nprocs = nprocs;
508 +        RBFNODE                 *near_rbf, *mir_rbf, *rbf;
509 +        double                  bestd;
510 +        int                     n;
511 +
512 +        if (dsf_list == NULL)
513 +                return;                         /* XXX should be error? */
514 +        near_rbf = dsf_list;
515 +        bestd = input_orient*near_rbf->invec[2];
516 +        if (single_plane_incident) {            /* ordered plane incidence? */
517 +                if (bestd >= 1.-2.*FTINY)
518 +                        return;                 /* already have normal */
519 +        } else {
520 +                switch (inp_coverage) {
521 +                case INP_QUAD1:
522 +                case INP_QUAD2:
523 +                case INP_QUAD3:
524 +                case INP_QUAD4:
525 +                        break;                  /* quadrilateral symmetry? */
526 +                default:
527 +                        return;                 /* else we can interpolate */
528 +                }
529 +                for (rbf = near_rbf->next; rbf != NULL; rbf = rbf->next) {
530 +                        const double    d = input_orient*rbf->invec[2];
531 +                        if (d >= 1.-2.*FTINY)
532 +                                return;         /* seems we have normal */
533 +                        if (d > bestd) {
534 +                                near_rbf = rbf;
535 +                                bestd = d;
536 +                        }
537 +                }
538 +        }
539 +        if (mig_list != NULL) {                 /* need to be called first */
540 +                fprintf(stderr, "%s: Late call to check_normal_incidence()\n",
541 +                                progname);
542 +                exit(1);
543 +        }
544 + #ifdef DEBUG
545 +        fprintf(stderr, "Interpolating normal incidence by mirroring (%.1f,%.1f)\n",
546 +                        get_theta180(near_rbf->invec), get_phi360(near_rbf->invec));
547 + #endif
548 +                                                /* mirror nearest incidence */
549 +        n = sizeof(RBFNODE) + sizeof(RBFVAL)*(near_rbf->nrbf-1);
550 +        mir_rbf = (RBFNODE *)malloc(n);
551 +        if (mir_rbf == NULL)
552 +                goto memerr;
553 +        memcpy(mir_rbf, near_rbf, n);
554 +        mir_rbf->ord = near_rbf->ord - 1;       /* not used, I think */
555 +        mir_rbf->next = NULL;
556 +        rev_rbf_symmetry(mir_rbf, MIRROR_X|MIRROR_Y);
557 +        nprocs = 1;                             /* compute migration matrix */
558 +        if (mig_list != create_migration(mir_rbf, near_rbf))
559 +                exit(1);                        /* XXX should never happen! */
560 +                                                /* interpolate normal dist. */
561 +        rbf = e_advect_rbf(mig_list, norm_vec, 2*near_rbf->nrbf);
562 +        nprocs = saved_nprocs;                  /* final clean-up */
563 +        free(mir_rbf);
564 +        free(mig_list);
565 +        mig_list = near_rbf->ejl = NULL;
566 +        insert_dsf(rbf);                        /* insert interpolated normal */
567 +        return;
568 + memerr:
569 +        fprintf(stderr, "%s: Out of memory in check_normal_incidence()\n",
570 +                                progname);
571 +        exit(1);
572 + }
573          
574   /* Build our triangle mesh from recorded RBFs */
575   void
# Line 492 | Line 578 | build_mesh(void)
578          double          best2 = M_PI*M_PI;
579          RBFNODE         *shrt_edj[2];
580          RBFNODE         *rbf0, *rbf1;
581 +                                                /* add normal if needed */
582 +        check_normal_incidence();
583                                                  /* check if isotropic */
584          if (single_plane_incident) {
585                  for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)

Diff Legend

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