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.19 by greg, Sat Mar 8 18:16:48 2014 UTC vs.
Revision 2.30 by greg, Thu Aug 21 10:33:48 2014 UTC

# Line 20 | Line 20 | static const char RCSid[] = "$Id$";
20   #include "bsdfrep.h"
21  
22   #ifndef NEIGH_FACT2
23 < #define NEIGH_FACT2     0.2     /* empirical neighborhood distance weight */
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) */
28   static int              nchild = 0;
29  
30 < 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 */
30 > /* Compute average DSF value at the given radius from central vector */
31 > static double
32 > eval_DSFsurround(const RBFNODE *rbf, const FVECT outvec, const double rad)
33 > {
34 >        const int       ninc = 12;
35 >        const double    phinc = 2.*M_PI/ninc;
36 >        double          sum = 0;
37 >        int             n = 0;
38 >        FVECT           tvec;
39 >        int             i;
40 >                                                /* compute initial vector */
41 >        if (output_orient*outvec[2] >= 1.-FTINY) {
42 >                tvec[0] = tvec[2] = 0;
43 >                tvec[1] = 1;
44 >        } else {
45 >                tvec[0] = tvec[1] = 0;
46 >                tvec[2] = 1;
47 >        }
48 >        geodesic(tvec, outvec, tvec, rad, GEOD_RAD);
49 >                                                /* average surrounding DSF */
50 >        for (i = 0; i < ninc; i++) {
51 >                if (i) spinvector(tvec, tvec, outvec, phinc);
52 >                if (tvec[2] > 0 ^ output_orient > 0)
53 >                        continue;
54 >                sum += eval_rbfrep(rbf, tvec) * output_orient*tvec[2];
55 >                ++n;
56 >        }
57 >        if (n < 2)                              /* should never happen! */
58 >                return(sum);
59 >        return(sum/(double)n);
60 > }
61  
62 < #define pricerow(p,i)   ((p)->price + (i)*(p)->ncols)
63 < #define psortrow(p,i)   ((p)->sord + (i)*(p)->ncols)
62 > /* Estimate single-lobe radius for DSF at the given outgoing angle */
63 > static double
64 > est_DSFrad(const RBFNODE *rbf, const FVECT outvec)
65 > {
66 >        const double    rad_epsilon = 0.03;
67 >        const double    DSFtarget = 0.60653066 * eval_rbfrep(rbf,outvec)
68 >                                                * output_orient*outvec[2];
69 >        double          inside_rad = rad_epsilon;
70 >        double          outside_rad = 0.5;
71 >        double          DSFinside = eval_DSFsurround(rbf, outvec, inside_rad);
72 >        double          DSFoutside = eval_DSFsurround(rbf, outvec, outside_rad);
73 > #define interp_rad      inside_rad + (outside_rad-inside_rad) * \
74 >                                (DSFtarget-DSFinside) / (DSFoutside-DSFinside)
75 >                                                /* interpolation search */
76 >        while (outside_rad-inside_rad > rad_epsilon) {
77 >                double  test_rad = interp_rad;
78 >                double  DSFtest = eval_DSFsurround(rbf, outvec, test_rad);
79 >                if (DSFtarget < DSFtest) {
80 >                        inside_rad = test_rad;
81 >                        DSFinside = DSFtest;
82 >                } else {
83 >                        outside_rad = test_rad;
84 >                        DSFoutside = DSFtest;
85 >                }
86 >        }
87 >        return(interp_rad);
88 > #undef interp_rad
89 > }
90  
91 + /* Compute average BSDF peak from current DSF's */
92 + static void
93 + comp_bsdf_spec(void)
94 + {
95 +        double          peak_sum = 0;
96 +        double          rad_sum = 0;
97 +        int             n = 0;
98 +        RBFNODE         *rbf;
99 +        FVECT           sdv;
100 +
101 +        if (dsf_list == NULL) {
102 +                bsdf_spec_peak = 0;
103 +                bsdf_spec_crad = 0;
104 +                return;
105 +        }
106 +        for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
107 +                sdv[0] = -rbf->invec[0];
108 +                sdv[1] = -rbf->invec[1];
109 +                sdv[2] = rbf->invec[2]*(2*(input_orient==output_orient) - 1);
110 +                peak_sum += eval_rbfrep(rbf, sdv);
111 +                rad_sum += est_DSFrad(rbf, sdv);
112 +                ++n;
113 +        }
114 +        bsdf_spec_peak = peak_sum/(double)n;
115 +        bsdf_spec_crad = ANG2R( rad_sum/(double)n );
116 + }
117 +
118   /* Create a new migration holder (sharing memory for multiprocessing) */
119   static MIGRATION *
120   new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
# Line 150 | Line 228 | compute_nDSFs(const RBFNODE *rbf0, const RBFNODE *rbf1
228  
229          for (x = GRIDRES; x--; )
230              for (y = GRIDRES; y--; ) {
231 <                ovec_from_pos(dv, x, y);
232 <                dsf_grid[x][y].val[0] = nf0 * eval_rbfrep(rbf0, dv);
233 <                dsf_grid[x][y].val[1] = nf1 * eval_rbfrep(rbf1, dv);
231 >                ovec_from_pos(dv, x, y);        /* cube root (brightness) */
232 >                dsf_grid[x][y].val[0] = pow(nf0*eval_rbfrep(rbf0, dv), .3333);
233 >                dsf_grid[x][y].val[1] = pow(nf1*eval_rbfrep(rbf1, dv), .3333);
234              }
235   }      
236  
# Line 165 | Line 243 | neighborhood_dist2(int x0, int y0, int x1, int y1)
243          double  d;
244          int     p[4];
245          int     i, j;
168
169        if ((x0 == x1) & (y0 == y1))
170                return(0.);
246                                                  /* check radius */
247          p[0] = x0; p[1] = y0; p[2] = x1; p[3] = y1;
248          for (i = 4; i--; ) {
# Line 183 | Line 258 | neighborhood_dist2(int x0, int y0, int x1, int y1)
258          return(sum2 / (4*rad*(rad+1) + 1));
259   }
260  
261 < /* Comparison routine needed for sorting price row */
262 < static int
263 < msrt_cmp(void *b, const void *p1, const void *p2)
261 > /* Compute distance between two RBF lobes */
262 > double
263 > lobe_distance(RBFVAL *rbf1, RBFVAL *rbf2)
264   {
265 <        PRICEMAT        *pm = (PRICEMAT *)b;
266 <        float           c1 = pm->prow[*(const short *)p1];
267 <        float           c2 = pm->prow[*(const short *)p2];
268 <
269 <        if (c1 > c2) return(1);
270 <        if (c1 < c2) return(-1);
271 <        return(0);
265 >        FVECT   vfrom, vto;
266 >        double  d, res;
267 >                                        /* quadratic cost function */
268 >        ovec_from_pos(vfrom, rbf1->gx, rbf1->gy);
269 >        ovec_from_pos(vto, rbf2->gx, rbf2->gy);
270 >        d = Acos(DOT(vfrom, vto));
271 >        res = d*d;
272 >        d = R2ANG(rbf2->crad) - R2ANG(rbf1->crad);
273 >        res += d*d;
274 >                                        /* neighborhood difference */
275 >        res += NEIGH_FACT2 * neighborhood_dist2( rbf1->gx, rbf1->gy,
276 >                                                rbf2->gx, rbf2->gy );
277 >        return(res);
278   }
279  
199 /* Compute (and allocate) migration price matrix for optimization */
200 static void
201 price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, const RBFNODE *to_rbf)
202 {
203        FVECT   *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
204        int     i, j;
280  
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);
210        pm->sord = (short *)malloc(sizeof(short) * pm->nrows*pm->ncols);
211        
212        if ((pm->price == NULL) | (pm->sord == NULL) | (vto == NULL)) {
213                fprintf(stderr, "%s: Out of memory in migration_costs()\n",
214                                progname);
215                exit(1);
216        }
217        for (j = to_rbf->nrbf; j--; )           /* save repetitive ops. */
218                ovec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);
219
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          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(srow, pm->ncols, sizeof(short), pm, &msrt_cmp);
240        }
241        free(vto);
242 }
243
244 /* Free price matrix */
245 static void
246 free_routes(PRICEMAT *pm)
247 {
248        free(pm->price); pm->price = NULL;
249        free(pm->sord); pm->sord = NULL;
250 }
251
252 /* Compute minimum (optimistic) cost for moving the given source material */
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;
260                                                /* move cheapest first */
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 * prow[d];
266                amt2move -= amt;
267        }
268        return(total_cost);
269 }
270
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, PRICEMAT *pm)
297 {
298        const int       max2check = 100;
299        const double    maxamt = 1./(double)pm->ncols;
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;  /* cost per amount moved */
306                double  amt;    /* amount we can move */
307        } cur, best;
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                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);
346                                                /* find best source & dest. */
347        best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
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            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                                                /* 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)
366                                        - src_cost[i];
367            dst_rem[cur.d] += cur.amt;          /* undo trial move */
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;
373        }
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 */
379        mtx_coef(mig,best.s,best.d) += best.amt;
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
281   /* Compute and insert migration along directed edge (may fork child) */
282   static MIGRATION *
283   create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
284   {
393        const double    end_thresh = 5e-6;
394        PRICEMAT        pmtx;
285          MIGRATION       *newmig;
396        double          *src_rem, *dst_rem;
397        double          total_rem = 1., move_amt;
286          int             i, j;
287                                                  /* check if exists already */
288          for (newmig = from_rbf->ejl; newmig != NULL;
# Line 414 | Line 302 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
302          newmig = new_migration(from_rbf, to_rbf);
303          if (run_subprocess())
304                  return(newmig);                 /* child continues */
417        price_routes(&pmtx, from_rbf, to_rbf);
418        src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
419        dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
420        if ((src_rem == NULL) | (dst_rem == NULL)) {
421                fprintf(stderr, "%s: Out of memory in create_migration()\n",
422                                progname);
423                exit(1);
424        }
425                                                /* starting quantities */
426        memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
427        for (i = from_rbf->nrbf; i--; )
428                src_rem[i] = rbf_volume(&from_rbf->rbfa[i]) / from_rbf->vtotal;
429        for (j = to_rbf->nrbf; j--; )
430                dst_rem[j] = rbf_volume(&to_rbf->rbfa[j]) / to_rbf->vtotal;
305  
306 <        do {                                    /* move a bit at a time */
307 <                move_amt = migration_step(newmig, src_rem, dst_rem, &pmtx);
308 <                total_rem -= move_amt;
435 <        } while ((total_rem > end_thresh) & (move_amt > 0));
306 >                                                /* compute transport plan */
307 >        compute_nDSFs(from_rbf, to_rbf);
308 >        plan_transport(newmig);
309  
310          for (i = from_rbf->nrbf; i--; ) {       /* normalize final matrix */
311              double      nf = rbf_volume(&from_rbf->rbfa[i]);
# Line 442 | Line 315 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
315                  mtx_coef(newmig,i,j) *= nf;     /* row now sums to 1.0 */
316          }
317          end_subprocess();                       /* exit here if subprocess */
445        free_routes(&pmtx);                     /* free working arrays */
446        free(src_rem);
447        free(dst_rem);
318          return(newmig);
319   }
320  
# Line 534 | Line 404 | mesh_from_edge(MIGRATION *edge)
404                                  ej1 = create_migration(tvert[0], edge->rbfv[1]);
405                          mesh_from_edge(ej0);
406                          mesh_from_edge(ej1);
407 +                        return;
408                  }
409 <        } else if (tvert[1] == NULL) {          /* grow mesh on left */
409 >        }
410 >        if (tvert[1] == NULL) {                 /* grow mesh on left */
411                  tvert[1] = find_chull_vert(edge->rbfv[1], edge->rbfv[0]);
412                  if (tvert[1] != NULL) {
413                          if (tvert[1]->ord > edge->rbfv[0]->ord)
# Line 556 | Line 428 | mesh_from_edge(MIGRATION *edge)
428   static void
429   check_normal_incidence(void)
430   {
431 <        static const FVECT      norm_vec = {.0, .0, 1.};
431 >        static FVECT            norm_vec = {.0, .0, 1.};
432          const int               saved_nprocs = nprocs;
433          RBFNODE                 *near_rbf, *mir_rbf, *rbf;
434          double                  bestd;
# Line 606 | Line 478 | check_normal_incidence(void)
478          memcpy(mir_rbf, near_rbf, n);
479          mir_rbf->ord = near_rbf->ord - 1;       /* not used, I think */
480          mir_rbf->next = NULL;
481 +        mir_rbf->ejl = NULL;
482          rev_rbf_symmetry(mir_rbf, MIRROR_X|MIRROR_Y);
483          nprocs = 1;                             /* compute migration matrix */
484 <        if (mig_list != create_migration(mir_rbf, near_rbf))
484 >        if (create_migration(mir_rbf, near_rbf) == NULL)
485                  exit(1);                        /* XXX should never happen! */
486 <                                                /* interpolate normal dist. */
487 <        rbf = e_advect_rbf(mig_list, norm_vec, 2*near_rbf->nrbf);
486 >        norm_vec[2] = input_orient;             /* interpolate normal dist. */
487 >        rbf = e_advect_rbf(mig_list, norm_vec, 0);
488          nprocs = saved_nprocs;                  /* final clean-up */
489          free(mir_rbf);
490          free(mig_list);
# Line 631 | Line 504 | build_mesh(void)
504          double          best2 = M_PI*M_PI;
505          RBFNODE         *shrt_edj[2];
506          RBFNODE         *rbf0, *rbf1;
507 +                                                /* average specular peak */
508 +        comp_bsdf_spec();
509                                                  /* add normal if needed */
510          check_normal_incidence();
511                                                  /* check if isotropic */

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines