ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdfmesh.c
Revision: 2.26
Committed: Wed Mar 26 00:11:30 2014 UTC (10 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.25: +40 -1 lines
Log Message:
Added routine to dump transport plan for debugging

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.26 static const char RCSid[] = "$Id: bsdfmesh.c,v 2.25 2014/03/24 06:07:46 greg Exp $";
3 greg 2.1 #endif
4     /*
5     * Create BSDF advection mesh from radial basis functions.
6     *
7     * G. Ward
8     */
9    
10     #ifndef _WIN32
11     #include <unistd.h>
12     #include <sys/wait.h>
13     #include <sys/mman.h>
14     #endif
15     #define _USE_MATH_DEFINES
16     #include <stdio.h>
17     #include <stdlib.h>
18     #include <string.h>
19     #include <math.h>
20     #include "bsdfrep.h"
21 greg 2.19
22     #ifndef NEIGH_FACT2
23 greg 2.21 #define NEIGH_FACT2 0.1 /* empirical neighborhood distance weight */
24 greg 2.19 #endif
25 greg 2.1 /* number of processes to run */
26     int nprocs = 1;
27     /* number of children (-1 in child) */
28     static int nchild = 0;
29    
30 greg 2.3 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 greg 2.10 float *prow; /* current price row */
35 greg 2.3 } PRICEMAT; /* sorted pricing matrix */
36    
37     #define pricerow(p,i) ((p)->price + (i)*(p)->ncols)
38     #define psortrow(p,i) ((p)->sord + (i)*(p)->ncols)
39    
40 greg 2.2 /* Create a new migration holder (sharing memory for multiprocessing) */
41     static MIGRATION *
42     new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
43     {
44     size_t memlen = sizeof(MIGRATION) +
45     sizeof(float)*(from_rbf->nrbf*to_rbf->nrbf - 1);
46     MIGRATION *newmig;
47     #ifdef _WIN32
48     if (nprocs > 1)
49     fprintf(stderr, "%s: warning - multiprocessing not supported\n",
50     progname);
51     nprocs = 1;
52     newmig = (MIGRATION *)malloc(memlen);
53     #else
54     if (nprocs <= 1) { /* single process? */
55     newmig = (MIGRATION *)malloc(memlen);
56     } else { /* else need to share memory */
57     newmig = (MIGRATION *)mmap(NULL, memlen, PROT_READ|PROT_WRITE,
58     MAP_ANON|MAP_SHARED, -1, 0);
59     if ((void *)newmig == MAP_FAILED)
60     newmig = NULL;
61     }
62     #endif
63     if (newmig == NULL) {
64     fprintf(stderr, "%s: cannot allocate new migration\n", progname);
65     exit(1);
66     }
67     newmig->rbfv[0] = from_rbf;
68     newmig->rbfv[1] = to_rbf;
69     /* insert in edge lists */
70     newmig->enxt[0] = from_rbf->ejl;
71     from_rbf->ejl = newmig;
72     newmig->enxt[1] = to_rbf->ejl;
73     to_rbf->ejl = newmig;
74     newmig->next = mig_list; /* push onto global list */
75     return(mig_list = newmig);
76     }
77    
78     #ifdef _WIN32
79     #define await_children(n) (void)(n)
80     #define run_subprocess() 0
81     #define end_subprocess() (void)0
82     #else
83    
84     /* Wait for the specified number of child processes to complete */
85     static void
86     await_children(int n)
87     {
88     int exit_status = 0;
89    
90     if (n > nchild)
91     n = nchild;
92     while (n-- > 0) {
93     int status;
94     if (wait(&status) < 0) {
95     fprintf(stderr, "%s: missing child(ren)!\n", progname);
96     nchild = 0;
97     break;
98     }
99     --nchild;
100     if (status) { /* something wrong */
101     if ((status = WEXITSTATUS(status)))
102     exit_status = status;
103     else
104     exit_status += !exit_status;
105     fprintf(stderr, "%s: subprocess died\n", progname);
106     n = nchild; /* wait for the rest */
107     }
108     }
109     if (exit_status)
110     exit(exit_status);
111     }
112    
113     /* Start child process if multiprocessing selected */
114     static pid_t
115     run_subprocess(void)
116     {
117     int status;
118     pid_t pid;
119    
120     if (nprocs <= 1) /* any children requested? */
121     return(0);
122     await_children(nchild + 1 - nprocs); /* free up child process */
123     if ((pid = fork())) {
124     if (pid < 0) {
125     fprintf(stderr, "%s: cannot fork subprocess\n",
126     progname);
127 greg 2.6 await_children(nchild);
128 greg 2.2 exit(1);
129     }
130     ++nchild; /* subprocess started */
131     return(pid);
132     }
133     nchild = -1;
134     return(0); /* put child to work */
135     }
136    
137     /* If we are in subprocess, call exit */
138     #define end_subprocess() if (nchild < 0) _exit(0); else
139    
140     #endif /* ! _WIN32 */
141    
142 greg 2.19 /* 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 greg 2.20 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 greg 2.19 }
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 greg 2.3 /* Comparison routine needed for sorting price row */
184     static int
185     msrt_cmp(void *b, const void *p1, const void *p2)
186     {
187     PRICEMAT *pm = (PRICEMAT *)b;
188 greg 2.10 float c1 = pm->prow[*(const short *)p1];
189     float c2 = pm->prow[*(const short *)p2];
190 greg 2.3
191     if (c1 > c2) return(1);
192     if (c1 < c2) return(-1);
193     return(0);
194     }
195    
196 greg 2.1 /* Compute (and allocate) migration price matrix for optimization */
197 greg 2.3 static void
198     price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, const RBFNODE *to_rbf)
199 greg 2.1 {
200     FVECT *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
201     int i, j;
202    
203 greg 2.19 compute_nDSFs(from_rbf, to_rbf);
204 greg 2.3 pm->nrows = from_rbf->nrbf;
205     pm->ncols = to_rbf->nrbf;
206     pm->price = (float *)malloc(sizeof(float) * pm->nrows*pm->ncols);
207     pm->sord = (short *)malloc(sizeof(short) * pm->nrows*pm->ncols);
208    
209     if ((pm->price == NULL) | (pm->sord == NULL) | (vto == NULL)) {
210 greg 2.1 fprintf(stderr, "%s: Out of memory in migration_costs()\n",
211     progname);
212     exit(1);
213     }
214     for (j = to_rbf->nrbf; j--; ) /* save repetitive ops. */
215     ovec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);
216    
217     for (i = from_rbf->nrbf; i--; ) {
218     const double from_ang = R2ANG(from_rbf->rbfa[i].crad);
219     FVECT vfrom;
220 greg 2.10 short *srow;
221 greg 2.1 ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
222 greg 2.10 pm->prow = pricerow(pm,i);
223     srow = psortrow(pm,i);
224 greg 2.3 for (j = to_rbf->nrbf; j--; ) {
225 greg 2.20 double d; /* quadratic cost function */
226 greg 2.18 d = Acos(DOT(vfrom, vto[j]));
227 greg 2.13 pm->prow[j] = d*d;
228     d = R2ANG(to_rbf->rbfa[j].crad) - from_ang;
229 greg 2.19 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 greg 2.10 srow[j] = j;
235 greg 2.3 }
236 greg 2.10 qsort_r(srow, pm->ncols, sizeof(short), pm, &msrt_cmp);
237 greg 2.1 }
238     free(vto);
239     }
240    
241 greg 2.3 /* Free price matrix */
242     static void
243     free_routes(PRICEMAT *pm)
244 greg 2.1 {
245 greg 2.3 free(pm->price); pm->price = NULL;
246     free(pm->sord); pm->sord = NULL;
247 greg 2.1 }
248    
249     /* Compute minimum (optimistic) cost for moving the given source material */
250     static double
251 greg 2.3 min_cost(double amt2move, const double *avail, const PRICEMAT *pm, int s)
252 greg 2.1 {
253 greg 2.11 const short *srow = psortrow(pm,s);
254     const float *prow = pricerow(pm,s);
255 greg 2.1 double total_cost = 0;
256 greg 2.3 int j;
257 greg 2.1 /* move cheapest first */
258 greg 2.11 for (j = 0; (j < pm->ncols) & (amt2move > FTINY); j++) {
259     int d = srow[j];
260 greg 2.1 double amt = (amt2move < avail[d]) ? amt2move : avail[d];
261    
262 greg 2.11 total_cost += amt * prow[d];
263 greg 2.1 amt2move -= amt;
264     }
265     return(total_cost);
266     }
267    
268 greg 2.17 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 greg 2.11 static int
275     rmovcmp(void *b, const void *p1, const void *p2)
276     {
277     PRICEMAT *pm = (PRICEMAT *)b;
278 greg 2.17 const ROWSENT *re1 = (const ROWSENT *)p1;
279     const ROWSENT *re2 = (const ROWSENT *)p2;
280     double price_diff;
281    
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 greg 2.11 if (price_diff > 0) return(1);
287     if (price_diff < 0) return(-1);
288     return(0);
289     }
290    
291     /* Take a step in migration by choosing reasonable bucket to transfer */
292 greg 2.1 static double
293 greg 2.11 migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, PRICEMAT *pm)
294 greg 2.1 {
295 greg 2.11 const int max2check = 100;
296 greg 2.4 const double maxamt = 1./(double)pm->ncols;
297 greg 2.12 const double minamt = maxamt*1e-4;
298 greg 2.5 double *src_cost;
299 greg 2.17 ROWSENT *rord;
300 greg 2.1 struct {
301     int s, d; /* source and destination */
302 greg 2.17 double price; /* cost per amount moved */
303 greg 2.1 double amt; /* amount we can move */
304     } cur, best;
305 greg 2.11 int r2check, i, ri;
306     /*
307     * Check cheapest available routes only -- a higher adjusted
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 greg 2.17 * 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 greg 2.11 */
315 greg 2.17 #define discount(qr) 1.0
316 greg 2.11 /* most promising row order */
317 greg 2.17 rord = (ROWSENT *)malloc(sizeof(ROWSENT)*pm->nrows);
318 greg 2.11 if (rord == NULL)
319     goto memerr;
320     for (ri = pm->nrows; ri--; ) {
321 greg 2.17 rord[ri].s = ri;
322     rord[ri].d = -1;
323     rord[ri].dc = 1.f;
324 greg 2.11 if (src_rem[ri] <= minamt) /* enough source material? */
325     continue;
326     for (i = 0; i < pm->ncols; i++)
327 greg 2.17 if (dst_rem[ rord[ri].d = psortrow(pm,ri)[i] ] > minamt)
328 greg 2.11 break;
329     if (i >= pm->ncols) { /* moved all we can? */
330     free(rord);
331     return(.0);
332     }
333 greg 2.17 rord[ri].dc = discount(src_rem[ri]);
334 greg 2.11 }
335     if (pm->nrows > max2check) /* sort if too many sources */
336 greg 2.17 qsort_r(rord, pm->nrows, sizeof(ROWSENT), pm, &rmovcmp);
337 greg 2.5 /* allocate cost array */
338     src_cost = (double *)malloc(sizeof(double)*pm->nrows);
339 greg 2.11 if (src_cost == NULL)
340     goto memerr;
341 greg 2.3 for (i = pm->nrows; i--; ) /* starting costs for diff. */
342     src_cost[i] = min_cost(src_rem[i], dst_rem, pm, i);
343 greg 2.1 /* find best source & dest. */
344     best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
345 greg 2.11 if ((r2check = pm->nrows) > max2check)
346     r2check = max2check; /* put a limit on search */
347     for (ri = 0; ri < r2check; ri++) { /* check each source row */
348 greg 2.1 double cost_others = 0;
349 greg 2.17 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 greg 2.11 if (pm->nrows > max2check) break; /* sorted end */
353     continue; /* else skip this one */
354     }
355 greg 2.1 cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
356     src_rem[cur.s] : dst_rem[cur.d];
357 greg 2.11 /* don't just leave smidgen */
358     if (cur.amt > maxamt*1.02) cur.amt = maxamt;
359     dst_rem[cur.d] -= cur.amt; /* add up opportunity costs */
360 greg 2.3 for (i = pm->nrows; i--; )
361 greg 2.1 if (i != cur.s)
362 greg 2.11 cost_others += min_cost(src_rem[i], dst_rem, pm, i)
363 greg 2.1 - src_cost[i];
364     dst_rem[cur.d] += cur.amt; /* undo trial move */
365 greg 2.17 /* discount effective price */
366     cur.price = ( pricerow(pm,cur.s)[cur.d] + cost_others/cur.amt ) *
367     rord[ri].dc;
368 greg 2.1 if (cur.price < best.price) /* are we better than best? */
369 greg 2.11 best = cur;
370 greg 2.1 }
371 greg 2.11 free(src_cost); /* clean up */
372     free(rord);
373 greg 2.5 if ((best.s < 0) | (best.d < 0)) /* nothing left to move? */
374 greg 2.1 return(.0);
375 greg 2.5 /* else make the actual move */
376 greg 2.2 mtx_coef(mig,best.s,best.d) += best.amt;
377 greg 2.1 src_rem[best.s] -= best.amt;
378     dst_rem[best.d] -= best.amt;
379     return(best.amt);
380 greg 2.11 memerr:
381     fprintf(stderr, "%s: Out of memory in migration_step()\n", progname);
382     exit(1);
383 greg 2.17 #undef discount
384 greg 2.1 }
385    
386 greg 2.26 #ifdef DUMP_MATRIX
387     /* Dump transport plan and corresponding price matrix to a text file */
388     static void
389     dump_matrix(const MIGRATION *me, const PRICEMAT *pm)
390     {
391     char fname[256];
392     FILE *fp;
393     int i, j;
394    
395     sprintf(fname, "edge_%d-%d.txt", me->rbfv[0]->ord, me->rbfv[1]->ord);
396     if ((fp = fopen(fname, "w")) == NULL)
397     return;
398     for (j = 0; j < 2; j++) {
399     fprintf(fp, "Available from %d source RBF lobes in node %d:\n",
400     me->rbfv[j]->nrbf, me->rbfv[j]->ord);
401     for (i = 0; i < me->rbfv[j]->nrbf; i++)
402     fprintf(fp, " %.4e", rbf_volume(&me->rbfv[j]->rbfa[i]) /
403     me->rbfv[j]->vtotal);
404     fputc('\n', fp);
405     }
406     fprintf(fp, "Price (quadratic distance metric) matrix:\n");
407     for (i = 0; i < pm->nrows; i++) {
408     for (j = 0; j < pm->ncols; j++)
409     fprintf(fp, " %.4e", pricerow(pm,i)[j]);
410     fputc('\n', fp);
411     }
412     fprintf(fp, "Solution matrix (transport plan):\n");
413     for (i = 0; i < mtx_nrows(me); i++) {
414     for (j = 0; j < mtx_ncols(me); j++)
415     fprintf(fp, " %.4e", mtx_coef(me,i,j));
416     fputc('\n', fp);
417     }
418     fclose(fp);
419     }
420     #endif
421    
422 greg 2.1 /* Compute and insert migration along directed edge (may fork child) */
423     static MIGRATION *
424     create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
425     {
426 greg 2.2 const double end_thresh = 5e-6;
427 greg 2.3 PRICEMAT pmtx;
428 greg 2.1 MIGRATION *newmig;
429     double *src_rem, *dst_rem;
430     double total_rem = 1., move_amt;
431 greg 2.6 int i, j;
432 greg 2.1 /* check if exists already */
433     for (newmig = from_rbf->ejl; newmig != NULL;
434     newmig = nextedge(from_rbf,newmig))
435     if (newmig->rbfv[1] == to_rbf)
436     return(NULL);
437     /* else allocate */
438 greg 2.7 #ifdef DEBUG
439 greg 2.14 fprintf(stderr, "Building path from (theta,phi) (%.1f,%.1f) ",
440 greg 2.7 get_theta180(from_rbf->invec),
441     get_phi360(from_rbf->invec));
442 greg 2.14 fprintf(stderr, "to (%.1f,%.1f) with %d x %d matrix\n",
443 greg 2.7 get_theta180(to_rbf->invec),
444     get_phi360(to_rbf->invec),
445     from_rbf->nrbf, to_rbf->nrbf);
446     #endif
447 greg 2.1 newmig = new_migration(from_rbf, to_rbf);
448     if (run_subprocess())
449     return(newmig); /* child continues */
450 greg 2.3 price_routes(&pmtx, from_rbf, to_rbf);
451 greg 2.1 src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
452     dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
453     if ((src_rem == NULL) | (dst_rem == NULL)) {
454     fprintf(stderr, "%s: Out of memory in create_migration()\n",
455     progname);
456     exit(1);
457     }
458     /* starting quantities */
459     memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
460     for (i = from_rbf->nrbf; i--; )
461     src_rem[i] = rbf_volume(&from_rbf->rbfa[i]) / from_rbf->vtotal;
462 greg 2.6 for (j = to_rbf->nrbf; j--; )
463     dst_rem[j] = rbf_volume(&to_rbf->rbfa[j]) / to_rbf->vtotal;
464    
465 greg 2.1 do { /* move a bit at a time */
466 greg 2.3 move_amt = migration_step(newmig, src_rem, dst_rem, &pmtx);
467 greg 2.1 total_rem -= move_amt;
468 greg 2.2 } while ((total_rem > end_thresh) & (move_amt > 0));
469 greg 2.6
470 greg 2.1 for (i = from_rbf->nrbf; i--; ) { /* normalize final matrix */
471 greg 2.6 double nf = rbf_volume(&from_rbf->rbfa[i]);
472 greg 2.1 if (nf <= FTINY) continue;
473     nf = from_rbf->vtotal / nf;
474     for (j = to_rbf->nrbf; j--; )
475 greg 2.6 mtx_coef(newmig,i,j) *= nf; /* row now sums to 1.0 */
476 greg 2.1 }
477 greg 2.26 #ifdef DUMP_MATRIX
478     dump_matrix(newmig, &pmtx);
479     #endif
480 greg 2.1 end_subprocess(); /* exit here if subprocess */
481 greg 2.3 free_routes(&pmtx); /* free working arrays */
482 greg 2.1 free(src_rem);
483     free(dst_rem);
484     return(newmig);
485     }
486    
487     /* Check if prospective vertex would create overlapping triangle */
488     static int
489     overlaps_tri(const RBFNODE *bv0, const RBFNODE *bv1, const RBFNODE *pv)
490     {
491     const MIGRATION *ej;
492     RBFNODE *vother[2];
493     int im_rev;
494     /* find shared edge in mesh */
495     for (ej = pv->ejl; ej != NULL; ej = nextedge(pv,ej)) {
496     const RBFNODE *tv = opp_rbf(pv,ej);
497     if (tv == bv0) {
498     im_rev = is_rev_tri(ej->rbfv[0]->invec,
499     ej->rbfv[1]->invec, bv1->invec);
500     break;
501     }
502     if (tv == bv1) {
503     im_rev = is_rev_tri(ej->rbfv[0]->invec,
504     ej->rbfv[1]->invec, bv0->invec);
505     break;
506     }
507     }
508     if (!get_triangles(vother, ej)) /* triangle on same side? */
509     return(0);
510     return(vother[im_rev] != NULL);
511     }
512    
513 greg 2.14 /* Find convex hull vertex to complete triangle (oriented call) */
514 greg 2.1 static RBFNODE *
515     find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rbf1)
516     {
517     FVECT vmid, vejn, vp;
518     RBFNODE *rbf, *rbfbest = NULL;
519     double dprod, area2, bestarea2 = FHUGE, bestdprod = -.5;
520    
521     VSUB(vejn, rbf1->invec, rbf0->invec);
522     VADD(vmid, rbf0->invec, rbf1->invec);
523     if (normalize(vejn) == 0 || normalize(vmid) == 0)
524     return(NULL);
525     /* XXX exhaustive search */
526     /* Find triangle with minimum rotation from perpendicular */
527     for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
528     if ((rbf == rbf0) | (rbf == rbf1))
529     continue;
530     tri_orient(vp, rbf0->invec, rbf1->invec, rbf->invec);
531     if (DOT(vp, vmid) <= FTINY)
532     continue; /* wrong orientation */
533     area2 = .25*DOT(vp,vp);
534 greg 2.14 VSUB(vp, rbf->invec, vmid);
535 greg 2.1 dprod = -DOT(vp, vejn);
536     VSUM(vp, vp, vejn, dprod); /* above guarantees non-zero */
537     dprod = DOT(vp, vmid) / VLEN(vp);
538     if (dprod <= bestdprod + FTINY*(1 - 2*(area2 < bestarea2)))
539     continue; /* found better already */
540     if (overlaps_tri(rbf0, rbf1, rbf))
541     continue; /* overlaps another triangle */
542     rbfbest = rbf;
543     bestdprod = dprod; /* new one to beat */
544     bestarea2 = area2;
545     }
546     return(rbfbest);
547     }
548    
549     /* Create new migration edge and grow mesh recursively around it */
550     static void
551     mesh_from_edge(MIGRATION *edge)
552     {
553     MIGRATION *ej0, *ej1;
554     RBFNODE *tvert[2];
555    
556     if (edge == NULL)
557     return;
558     /* triangle on either side? */
559     get_triangles(tvert, edge);
560     if (tvert[0] == NULL) { /* grow mesh on right */
561     tvert[0] = find_chull_vert(edge->rbfv[0], edge->rbfv[1]);
562     if (tvert[0] != NULL) {
563     if (tvert[0]->ord > edge->rbfv[0]->ord)
564     ej0 = create_migration(edge->rbfv[0], tvert[0]);
565     else
566     ej0 = create_migration(tvert[0], edge->rbfv[0]);
567     if (tvert[0]->ord > edge->rbfv[1]->ord)
568     ej1 = create_migration(edge->rbfv[1], tvert[0]);
569     else
570     ej1 = create_migration(tvert[0], edge->rbfv[1]);
571     mesh_from_edge(ej0);
572     mesh_from_edge(ej1);
573     }
574     } else if (tvert[1] == NULL) { /* grow mesh on left */
575     tvert[1] = find_chull_vert(edge->rbfv[1], edge->rbfv[0]);
576     if (tvert[1] != NULL) {
577     if (tvert[1]->ord > edge->rbfv[0]->ord)
578     ej0 = create_migration(edge->rbfv[0], tvert[1]);
579     else
580     ej0 = create_migration(tvert[1], edge->rbfv[0]);
581     if (tvert[1]->ord > edge->rbfv[1]->ord)
582     ej1 = create_migration(edge->rbfv[1], tvert[1]);
583     else
584     ej1 = create_migration(tvert[1], edge->rbfv[1]);
585     mesh_from_edge(ej0);
586     mesh_from_edge(ej1);
587     }
588     }
589     }
590 greg 2.15
591     /* Add normal direction if missing */
592     static void
593     check_normal_incidence(void)
594     {
595 greg 2.25 static FVECT norm_vec = {.0, .0, 1.};
596 greg 2.16 const int saved_nprocs = nprocs;
597     RBFNODE *near_rbf, *mir_rbf, *rbf;
598     double bestd;
599     int n;
600 greg 2.15
601     if (dsf_list == NULL)
602     return; /* XXX should be error? */
603     near_rbf = dsf_list;
604     bestd = input_orient*near_rbf->invec[2];
605     if (single_plane_incident) { /* ordered plane incidence? */
606     if (bestd >= 1.-2.*FTINY)
607     return; /* already have normal */
608     } else {
609     switch (inp_coverage) {
610     case INP_QUAD1:
611     case INP_QUAD2:
612     case INP_QUAD3:
613     case INP_QUAD4:
614     break; /* quadrilateral symmetry? */
615     default:
616     return; /* else we can interpolate */
617     }
618     for (rbf = near_rbf->next; rbf != NULL; rbf = rbf->next) {
619     const double d = input_orient*rbf->invec[2];
620     if (d >= 1.-2.*FTINY)
621     return; /* seems we have normal */
622     if (d > bestd) {
623     near_rbf = rbf;
624     bestd = d;
625     }
626     }
627     }
628     if (mig_list != NULL) { /* need to be called first */
629     fprintf(stderr, "%s: Late call to check_normal_incidence()\n",
630     progname);
631     exit(1);
632     }
633     #ifdef DEBUG
634     fprintf(stderr, "Interpolating normal incidence by mirroring (%.1f,%.1f)\n",
635     get_theta180(near_rbf->invec), get_phi360(near_rbf->invec));
636     #endif
637     /* mirror nearest incidence */
638     n = sizeof(RBFNODE) + sizeof(RBFVAL)*(near_rbf->nrbf-1);
639     mir_rbf = (RBFNODE *)malloc(n);
640     if (mir_rbf == NULL)
641     goto memerr;
642     memcpy(mir_rbf, near_rbf, n);
643     mir_rbf->ord = near_rbf->ord - 1; /* not used, I think */
644     mir_rbf->next = NULL;
645 greg 2.22 mir_rbf->ejl = NULL;
646 greg 2.15 rev_rbf_symmetry(mir_rbf, MIRROR_X|MIRROR_Y);
647     nprocs = 1; /* compute migration matrix */
648 greg 2.22 if (create_migration(mir_rbf, near_rbf) == NULL)
649 greg 2.15 exit(1); /* XXX should never happen! */
650 greg 2.25 norm_vec[2] = input_orient; /* interpolate normal dist. */
651 greg 2.16 rbf = e_advect_rbf(mig_list, norm_vec, 2*near_rbf->nrbf);
652 greg 2.15 nprocs = saved_nprocs; /* final clean-up */
653     free(mir_rbf);
654     free(mig_list);
655     mig_list = near_rbf->ejl = NULL;
656     insert_dsf(rbf); /* insert interpolated normal */
657     return;
658     memerr:
659     fprintf(stderr, "%s: Out of memory in check_normal_incidence()\n",
660     progname);
661     exit(1);
662     }
663 greg 2.1
664     /* Build our triangle mesh from recorded RBFs */
665     void
666     build_mesh(void)
667     {
668     double best2 = M_PI*M_PI;
669     RBFNODE *shrt_edj[2];
670     RBFNODE *rbf0, *rbf1;
671 greg 2.15 /* add normal if needed */
672     check_normal_incidence();
673 greg 2.1 /* check if isotropic */
674     if (single_plane_incident) {
675     for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
676     if (rbf0->next != NULL)
677     create_migration(rbf0, rbf0->next);
678     await_children(nchild);
679     return;
680     }
681     shrt_edj[0] = shrt_edj[1] = NULL; /* start w/ shortest edge */
682     for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
683     for (rbf1 = rbf0->next; rbf1 != NULL; rbf1 = rbf1->next) {
684     double dist2 = 2. - 2.*DOT(rbf0->invec,rbf1->invec);
685     if (dist2 < best2) {
686     shrt_edj[0] = rbf0;
687     shrt_edj[1] = rbf1;
688     best2 = dist2;
689     }
690     }
691     if (shrt_edj[0] == NULL) {
692     fprintf(stderr, "%s: Cannot find shortest edge\n", progname);
693     exit(1);
694     }
695     /* build mesh from this edge */
696     if (shrt_edj[0]->ord < shrt_edj[1]->ord)
697     mesh_from_edge(create_migration(shrt_edj[0], shrt_edj[1]));
698     else
699     mesh_from_edge(create_migration(shrt_edj[1], shrt_edj[0]));
700     /* complete migrations */
701     await_children(nchild);
702     }