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.17 by greg, Wed Mar 5 22:47:16 2014 UTC vs.
Revision 2.40 by greg, Tue Apr 23 14:30:36 2019 UTC

# Line 7 | Line 7 | static const char RCSid[] = "$Id$";
7   *      G. Ward
8   */
9  
10 < #ifndef _WIN32
10 > #if !defined(_WIN32) && !defined(_WIN64)
11   #include <unistd.h>
12   #include <sys/wait.h>
13   #include <sys/mman.h>
# 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) */
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) * COSF(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.01;
67 >        const double    DSFtarget = 0.60653066 * eval_rbfrep(rbf,outvec) *
68 >                                                        COSF(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 >                                                /* Newton's method (sort of) */
76 >        do {
77 >                double  test_rad = interp_rad;
78 >                double  DSFtest;
79 >                if ((test_rad >= outside_rad) | (test_rad <= inside_rad))
80 >                        test_rad = .5*(inside_rad + outside_rad);
81 >                DSFtest = eval_DSFsurround(rbf, outvec, test_rad);
82 >                if (DSFtest > DSFtarget) {
83 >                        inside_rad = test_rad;
84 >                        DSFinside = DSFtest;
85 >                } else {
86 >                        outside_rad = test_rad;
87 >                        DSFoutside = DSFtest;
88 >                }
89 >        } while (outside_rad-inside_rad > rad_epsilon);
90  
91 +        return(.5*(inside_rad + outside_rad));
92 + #undef interp_rad
93 + }
94 +
95 + static int
96 + dbl_cmp(const void *p1, const void *p2)
97 + {
98 +        double  d1 = *(const double *)p1;
99 +        double  d2 = *(const double *)p2;
100 +
101 +        if (d1 > d2) return(1);
102 +        if (d1 < d2) return(-1);
103 +        return(0);
104 + }
105 +
106 + /* Conservative estimate of average BSDF value from current DSF's */
107 + static void
108 + comp_bsdf_spec(void)
109 + {
110 +        double          vmod_sum = 0;
111 +        double          rad_sum = 0;
112 +        int             n = 0;
113 +        double          *cost_list = NULL;
114 +        double          max_cost = 1.;
115 +        RBFNODE         *rbf;
116 +        FVECT           sdv;
117 +                                                /* sort by incident altitude */
118 +        for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
119 +                n++;
120 +        if (n >= 10)
121 +                cost_list = (double *)malloc(sizeof(double)*n);
122 +        if (cost_list == NULL) {
123 +                bsdf_spec_val = 0;
124 +                bsdf_spec_rad = 0;
125 +                return;
126 +        }
127 +        n = 0;
128 +        for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
129 +                cost_list[n++] = rbf->invec[2]*input_orient;
130 +        qsort(cost_list, n, sizeof(double), dbl_cmp);
131 +        max_cost = cost_list[(n+3)/4];          /* accept 25% nearest grazing */
132 +        free(cost_list);
133 +        n = 0;
134 +        for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
135 +                double  this_rad, cosfact, vest;
136 +                if (rbf->invec[2]*input_orient > max_cost)
137 +                        continue;
138 +                sdv[0] = -rbf->invec[0];
139 +                sdv[1] = -rbf->invec[1];
140 +                sdv[2] = rbf->invec[2]*(2*(input_orient==output_orient) - 1);
141 +                cosfact = COSF(sdv[2]);
142 +                this_rad = est_DSFrad(rbf, sdv);
143 +                vest = eval_rbfrep(rbf, sdv) * cosfact *
144 +                                (2.*M_PI) * this_rad*this_rad;
145 +                if (vest > rbf->vtotal)         /* don't over-estimate energy */
146 +                        vest = rbf->vtotal;
147 +                vmod_sum += vest / cosfact;     /* remove cosine factor */
148 +                rad_sum += this_rad;
149 +                ++n;
150 +        }
151 +        bsdf_spec_rad = rad_sum/(double)n;
152 +        bsdf_spec_val = vmod_sum/(2.*M_PI*n*bsdf_spec_rad*bsdf_spec_rad);
153 + }
154 +
155   /* Create a new migration holder (sharing memory for multiprocessing) */
156   static MIGRATION *
157   new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
# Line 40 | Line 159 | new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
159          size_t          memlen = sizeof(MIGRATION) +
160                                  sizeof(float)*(from_rbf->nrbf*to_rbf->nrbf - 1);
161          MIGRATION       *newmig;
162 < #ifdef _WIN32
162 > #if defined(_WIN32) || defined(_WIN64)
163          if (nprocs > 1)
164                  fprintf(stderr, "%s: warning - multiprocessing not supported\n",
165                                  progname);
# Line 71 | Line 190 | new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
190          return(mig_list = newmig);
191   }
192  
193 < #ifdef _WIN32
193 > #if defined(_WIN32) || defined(_WIN64)
194   #define await_children(n)       (void)(n)
195   #define run_subprocess()        0
196   #define end_subprocess()        (void)0
# Line 135 | Line 254 | run_subprocess(void)
254  
255   #endif  /* ! _WIN32 */
256  
257 < /* Comparison routine needed for sorting price row */
139 < static int
140 < msrt_cmp(void *b, const void *p1, const void *p2)
141 < {
142 <        PRICEMAT        *pm = (PRICEMAT *)b;
143 <        float           c1 = pm->prow[*(const short *)p1];
144 <        float           c2 = pm->prow[*(const short *)p2];
145 <
146 <        if (c1 > c2) return(1);
147 <        if (c1 < c2) return(-1);
148 <        return(0);
149 < }
150 <
151 < /* Compute (and allocate) migration price matrix for optimization */
257 > /* Compute normalized distribution scattering functions for comparison */
258   static void
259 < price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, const RBFNODE *to_rbf)
259 > compute_nDSFs(const RBFNODE *rbf0, const RBFNODE *rbf1)
260   {
261 <        FVECT   *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
262 <        int     i, j;
261 >        const double    nf0 = (GRIDRES*GRIDRES) / rbf0->vtotal;
262 >        const double    nf1 = (GRIDRES*GRIDRES) / rbf1->vtotal;
263 >        int             x, y;
264 >        FVECT           dv;
265  
266 <        pm->nrows = from_rbf->nrbf;
267 <        pm->ncols = to_rbf->nrbf;
268 <        pm->price = (float *)malloc(sizeof(float) * pm->nrows*pm->ncols);
269 <        pm->sord = (short *)malloc(sizeof(short) * pm->nrows*pm->ncols);
270 <        
163 <        if ((pm->price == NULL) | (pm->sord == NULL) | (vto == NULL)) {
164 <                fprintf(stderr, "%s: Out of memory in migration_costs()\n",
165 <                                progname);
166 <                exit(1);
167 <        }
168 <        for (j = to_rbf->nrbf; j--; )           /* save repetitive ops. */
169 <                ovec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);
170 <
171 <        for (i = from_rbf->nrbf; i--; ) {
172 <            const double        from_ang = R2ANG(from_rbf->rbfa[i].crad);
173 <            FVECT               vfrom;
174 <            short               *srow;
175 <            ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
176 <            pm->prow = pricerow(pm,i);
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);
182 <                pm->prow[j] = d*d;
183 <                d = R2ANG(to_rbf->rbfa[j].crad) - from_ang;
184 <                pm->prow[j] += d*d;    
185 <                srow[j] = j;
266 >        for (x = GRIDRES; x--; )
267 >            for (y = GRIDRES; y--; ) {
268 >                ovec_from_pos(dv, x, y);        /* cube root (brightness) */
269 >                dsf_grid[x][y].val[0] = pow(nf0*eval_rbfrep(rbf0, dv), .3333);
270 >                dsf_grid[x][y].val[1] = pow(nf1*eval_rbfrep(rbf1, dv), .3333);
271              }
272 <            qsort_r(srow, pm->ncols, sizeof(short), pm, &msrt_cmp);
188 <        }
189 <        free(vto);
190 < }
272 > }      
273  
274 < /* Free price matrix */
193 < static void
194 < free_routes(PRICEMAT *pm)
195 < {
196 <        free(pm->price); pm->price = NULL;
197 <        free(pm->sord); pm->sord = NULL;
198 < }
199 <
200 < /* Compute minimum (optimistic) cost for moving the given source material */
274 > /* Compute neighborhood distance-squared (dissimilarity) */
275   static double
276 < min_cost(double amt2move, const double *avail, const PRICEMAT *pm, int s)
276 > neighborhood_dist2(int x0, int y0, int x1, int y1)
277   {
278 <        const short     *srow = psortrow(pm,s);
279 <        const float     *prow = pricerow(pm,s);
280 <        double          total_cost = 0;
281 <        int             j;
282 <                                                /* move cheapest first */
283 <        for (j = 0; (j < pm->ncols) & (amt2move > FTINY); j++) {
284 <                int     d = srow[j];
285 <                double  amt = (amt2move < avail[d]) ? amt2move : avail[d];
286 <
287 <                total_cost += amt * prow[d];
214 <                amt2move -= amt;
278 >        int     rad = GRIDRES>>5;
279 >        double  sum2 = 0.;
280 >        double  d;
281 >        int     p[4];
282 >        int     i, j;
283 >                                                /* check radius */
284 >        p[0] = x0; p[1] = y0; p[2] = x1; p[3] = y1;
285 >        for (i = 4; i--; ) {
286 >                if (p[i] < rad) rad = p[i];
287 >                if (GRIDRES-1-p[i] < rad) rad = GRIDRES-1-p[i];
288          }
289 <        return(total_cost);
289 >        for (i = -rad; i <= rad; i++)
290 >            for (j = -rad; j <= rad; j++) {
291 >                d = dsf_grid[x0+i][y0+j].val[0] -
292 >                        dsf_grid[x1+i][y1+j].val[1];
293 >                sum2 += d*d;
294 >            }
295 >        return(sum2 / (4*rad*(rad+1) + 1));
296   }
297  
298 < typedef struct {
299 <        short   s, d;           /* source and destination */
300 <        float   dc;             /* discount to push inventory */
222 < } ROWSENT;              /* row sort entry */
223 <
224 < /* Compare entries by discounted moving price */
225 < static int
226 < rmovcmp(void *b, const void *p1, const void *p2)
298 > /* Compute distance between two RBF lobes */
299 > double
300 > lobe_distance(RBFVAL *rbf1, RBFVAL *rbf2)
301   {
302 <        PRICEMAT        *pm = (PRICEMAT *)b;
303 <        const ROWSENT   *re1 = (const ROWSENT *)p1;
304 <        const ROWSENT   *re2 = (const ROWSENT *)p2;
305 <        double          price_diff;
306 <
307 <        if (re1->d < 0) return(re2->d >= 0);
308 <        if (re2->d < 0) return(-1);
309 <        price_diff = re1->dc*pricerow(pm,re1->s)[re1->d] -
310 <                        re2->dc*pricerow(pm,re2->s)[re2->d];
311 <        if (price_diff > 0) return(1);
312 <        if (price_diff < 0) return(-1);
313 <        return(0);
302 >        FVECT   vfrom, vto;
303 >        double  d, res;
304 >                                        /* quadratic cost function */
305 >        ovec_from_pos(vfrom, rbf1->gx, rbf1->gy);
306 >        ovec_from_pos(vto, rbf2->gx, rbf2->gy);
307 >        d = Acos(DOT(vfrom, vto));
308 >        res = d*d;
309 >        d = R2ANG(rbf2->crad) - R2ANG(rbf1->crad);
310 >        res += d*d;
311 >                                        /* neighborhood difference */
312 >        res += NEIGH_FACT2 * neighborhood_dist2( rbf1->gx, rbf1->gy,
313 >                                                rbf2->gx, rbf2->gy );
314 >        return(res);
315   }
316  
242 – /* Take a step in migration by choosing reasonable bucket to transfer */
243 – static double
244 – migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, PRICEMAT *pm)
245 – {
246 –        const int       max2check = 100;
247 –        const double    maxamt = 1./(double)pm->ncols;
248 –        const double    minamt = maxamt*1e-4;
249 –        double          *src_cost;
250 –        ROWSENT         *rord;
251 –        struct {
252 –                int     s, d;   /* source and destination */
253 –                double  price;  /* cost per amount moved */
254 –                double  amt;    /* amount we can move */
255 –        } cur, best;
256 –        int             r2check, i, ri;
257 –        /*
258 –         * Check cheapest available routes only -- a higher adjusted
259 –         * destination price implies that another source is closer, so
260 –         * we can hold off considering more expensive options until
261 –         * some other (hopefully better) moves have been made.
262 –         * A discount based on source remaining is supposed to prioritize
263 –         * movement from large lobes, but it doesn't seem to do much,
264 –         * so we have it set to 1.0 at the moment.
265 –         */
266 – #define discount(qr)    1.0
267 –                                                /* most promising row order */
268 –        rord = (ROWSENT *)malloc(sizeof(ROWSENT)*pm->nrows);
269 –        if (rord == NULL)
270 –                goto memerr;
271 –        for (ri = pm->nrows; ri--; ) {
272 –            rord[ri].s = ri;
273 –            rord[ri].d = -1;
274 –            rord[ri].dc = 1.f;
275 –            if (src_rem[ri] <= minamt)          /* enough source material? */
276 –                    continue;
277 –            for (i = 0; i < pm->ncols; i++)
278 –                if (dst_rem[ rord[ri].d = psortrow(pm,ri)[i] ] > minamt)
279 –                        break;
280 –            if (i >= pm->ncols) {               /* moved all we can? */
281 –                free(rord);
282 –                return(.0);
283 –            }
284 –            rord[ri].dc = discount(src_rem[ri]);
285 –        }
286 –        if (pm->nrows > max2check)              /* sort if too many sources */
287 –                qsort_r(rord, pm->nrows, sizeof(ROWSENT), pm, &rmovcmp);
288 –                                                /* allocate cost array */
289 –        src_cost = (double *)malloc(sizeof(double)*pm->nrows);
290 –        if (src_cost == NULL)
291 –                goto memerr;
292 –        for (i = pm->nrows; i--; )              /* starting costs for diff. */
293 –                src_cost[i] = min_cost(src_rem[i], dst_rem, pm, i);
294 –                                                /* find best source & dest. */
295 –        best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
296 –        if ((r2check = pm->nrows) > max2check)
297 –                r2check = max2check;            /* put a limit on search */
298 –        for (ri = 0; ri < r2check; ri++) {      /* check each source row */
299 –            double      cost_others = 0;
300 –            cur.s = rord[ri].s;
301 –            if ((cur.d = rord[ri].d) < 0 ||
302 –                        rord[ri].dc*pricerow(pm,cur.s)[cur.d] >= best.price) {
303 –                if (pm->nrows > max2check) break;       /* sorted end */
304 –                continue;                       /* else skip this one */
305 –            }
306 –            cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
307 –                                src_rem[cur.s] : dst_rem[cur.d];
308 –                                                /* don't just leave smidgen */
309 –            if (cur.amt > maxamt*1.02) cur.amt = maxamt;
310 –            dst_rem[cur.d] -= cur.amt;          /* add up opportunity costs */
311 –            for (i = pm->nrows; i--; )
312 –                if (i != cur.s)
313 –                    cost_others += min_cost(src_rem[i], dst_rem, pm, i)
314 –                                        - src_cost[i];
315 –            dst_rem[cur.d] += cur.amt;          /* undo trial move */
316 –                                                /* discount effective price */
317 –            cur.price = ( pricerow(pm,cur.s)[cur.d] + cost_others/cur.amt ) *
318 –                                        rord[ri].dc;
319 –            if (cur.price < best.price)         /* are we better than best? */
320 –                best = cur;
321 –        }
322 –        free(src_cost);                         /* clean up */
323 –        free(rord);
324 –        if ((best.s < 0) | (best.d < 0))        /* nothing left to move? */
325 –                return(.0);
326 –                                                /* else make the actual move */
327 –        mtx_coef(mig,best.s,best.d) += best.amt;
328 –        src_rem[best.s] -= best.amt;
329 –        dst_rem[best.d] -= best.amt;
330 –        return(best.amt);
331 – memerr:
332 –        fprintf(stderr, "%s: Out of memory in migration_step()\n", progname);
333 –        exit(1);
334 – #undef discount
335 – }
317  
318   /* Compute and insert migration along directed edge (may fork child) */
319   static MIGRATION *
320   create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
321   {
341 –        const double    end_thresh = 5e-6;
342 –        PRICEMAT        pmtx;
322          MIGRATION       *newmig;
344 –        double          *src_rem, *dst_rem;
345 –        double          total_rem = 1., move_amt;
323          int             i, j;
324                                                  /* check if exists already */
325          for (newmig = from_rbf->ejl; newmig != NULL;
# Line 362 | Line 339 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
339          newmig = new_migration(from_rbf, to_rbf);
340          if (run_subprocess())
341                  return(newmig);                 /* child continues */
365 –        price_routes(&pmtx, from_rbf, to_rbf);
366 –        src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
367 –        dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
368 –        if ((src_rem == NULL) | (dst_rem == NULL)) {
369 –                fprintf(stderr, "%s: Out of memory in create_migration()\n",
370 –                                progname);
371 –                exit(1);
372 –        }
373 –                                                /* starting quantities */
374 –        memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
375 –        for (i = from_rbf->nrbf; i--; )
376 –                src_rem[i] = rbf_volume(&from_rbf->rbfa[i]) / from_rbf->vtotal;
377 –        for (j = to_rbf->nrbf; j--; )
378 –                dst_rem[j] = rbf_volume(&to_rbf->rbfa[j]) / to_rbf->vtotal;
342  
343 <        do {                                    /* move a bit at a time */
344 <                move_amt = migration_step(newmig, src_rem, dst_rem, &pmtx);
345 <                total_rem -= move_amt;
383 <        } while ((total_rem > end_thresh) & (move_amt > 0));
343 >                                                /* compute transport plan */
344 >        compute_nDSFs(from_rbf, to_rbf);
345 >        plan_transport(newmig);
346  
347          for (i = from_rbf->nrbf; i--; ) {       /* normalize final matrix */
348              double      nf = rbf_volume(&from_rbf->rbfa[i]);
# Line 390 | Line 352 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
352                  mtx_coef(newmig,i,j) *= nf;     /* row now sums to 1.0 */
353          }
354          end_subprocess();                       /* exit here if subprocess */
393 –        free_routes(&pmtx);                     /* free working arrays */
394 –        free(src_rem);
395 –        free(dst_rem);
355          return(newmig);
356   }
357  
# Line 482 | Line 441 | mesh_from_edge(MIGRATION *edge)
441                                  ej1 = create_migration(tvert[0], edge->rbfv[1]);
442                          mesh_from_edge(ej0);
443                          mesh_from_edge(ej1);
444 +                        return;
445                  }
446 <        } else if (tvert[1] == NULL) {          /* grow mesh on left */
446 >        }
447 >        if (tvert[1] == NULL) {                 /* grow mesh on left */
448                  tvert[1] = find_chull_vert(edge->rbfv[1], edge->rbfv[0]);
449                  if (tvert[1] != NULL) {
450                          if (tvert[1]->ord > edge->rbfv[0]->ord)
# Line 504 | Line 465 | mesh_from_edge(MIGRATION *edge)
465   static void
466   check_normal_incidence(void)
467   {
468 <        static const FVECT      norm_vec = {.0, .0, 1.};
468 >        static FVECT            norm_vec = {.0, .0, 1.};
469          const int               saved_nprocs = nprocs;
470          RBFNODE                 *near_rbf, *mir_rbf, *rbf;
471          double                  bestd;
# Line 527 | Line 488 | check_normal_incidence(void)
488                  default:
489                          return;                 /* else we can interpolate */
490                  }
491 <                for (rbf = near_rbf->next; rbf != NULL; rbf = rbf->next) {
491 >                for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
492                          const double    d = input_orient*rbf->invec[2];
493                          if (d >= 1.-2.*FTINY)
494                                  return;         /* seems we have normal */
# Line 554 | Line 515 | check_normal_incidence(void)
515          memcpy(mir_rbf, near_rbf, n);
516          mir_rbf->ord = near_rbf->ord - 1;       /* not used, I think */
517          mir_rbf->next = NULL;
518 +        mir_rbf->ejl = NULL;
519          rev_rbf_symmetry(mir_rbf, MIRROR_X|MIRROR_Y);
520          nprocs = 1;                             /* compute migration matrix */
521 <        if (mig_list != create_migration(mir_rbf, near_rbf))
521 >        if (create_migration(mir_rbf, near_rbf) == NULL)
522                  exit(1);                        /* XXX should never happen! */
523 <                                                /* interpolate normal dist. */
524 <        rbf = e_advect_rbf(mig_list, norm_vec, 2*near_rbf->nrbf);
523 >        norm_vec[2] = input_orient;             /* interpolate normal dist. */
524 >        rbf = e_advect_rbf(mig_list, norm_vec, 0);
525          nprocs = saved_nprocs;                  /* final clean-up */
526          free(mir_rbf);
527          free(mig_list);
# Line 576 | Line 538 | memerr:
538   void
539   build_mesh(void)
540   {
541 +        int             nrbfs = 0, nmigs = 0;
542          double          best2 = M_PI*M_PI;
543          RBFNODE         *shrt_edj[2];
544          RBFNODE         *rbf0, *rbf1;
545 +        const MIGRATION *ej;
546 +                                                /* average specular peak */
547 +        comp_bsdf_spec();
548                                                  /* add normal if needed */
549          check_normal_incidence();
550                                                  /* check if isotropic */
# Line 590 | Line 556 | build_mesh(void)
556                  return;
557          }
558          shrt_edj[0] = shrt_edj[1] = NULL;       /* start w/ shortest edge */
559 <        for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
559 >        for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next) {
560              for (rbf1 = rbf0->next; rbf1 != NULL; rbf1 = rbf1->next) {
561                  double  dist2 = 2. - 2.*DOT(rbf0->invec,rbf1->invec);
562                  if (dist2 < best2) {
# Line 598 | Line 564 | build_mesh(void)
564                          shrt_edj[1] = rbf1;
565                          best2 = dist2;
566                  }
567 +            }
568 +            ++nrbfs;
569          }
570          if (shrt_edj[0] == NULL) {
571                  fprintf(stderr, "%s: Cannot find shortest edge\n", progname);
# Line 608 | Line 576 | build_mesh(void)
576                  mesh_from_edge(create_migration(shrt_edj[0], shrt_edj[1]));
577          else
578                  mesh_from_edge(create_migration(shrt_edj[1], shrt_edj[0]));
579 +                                                /* count up edges */
580 +        for (ej = mig_list; ej != NULL; ej = ej->next)
581 +                ++nmigs;
582 +        if (nmigs < nrbfs-1)                    /* did meshing fail? */
583 +                fprintf(stderr,
584 +            "%s: warning - %d incident directions but only %d interpolant(s)\n",
585 +                                progname, nrbfs, nmigs);
586                                                  /* complete migrations */
587          await_children(nchild);
588   }

Diff Legend

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