ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdfmesh.c
Revision: 2.16
Committed: Wed Feb 19 05:16:06 2014 UTC (10 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.15: +8 -44 lines
Log Message:
Eliminated redundant code added in last change to bsdfmesh.c

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.16 static const char RCSid[] = "$Id: bsdfmesh.c,v 2.15 2014/02/18 16:06:51 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     /* number of processes to run */
22     int nprocs = 1;
23     /* number of children (-1 in child) */
24     static int nchild = 0;
25    
26 greg 2.3 typedef struct {
27     int nrows, ncols; /* array size (matches migration) */
28     float *price; /* migration prices */
29     short *sord; /* sort for each row, low to high */
30 greg 2.10 float *prow; /* current price row */
31 greg 2.3 } PRICEMAT; /* sorted pricing matrix */
32    
33     #define pricerow(p,i) ((p)->price + (i)*(p)->ncols)
34     #define psortrow(p,i) ((p)->sord + (i)*(p)->ncols)
35    
36 greg 2.2 /* Create a new migration holder (sharing memory for multiprocessing) */
37     static MIGRATION *
38     new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
39     {
40     size_t memlen = sizeof(MIGRATION) +
41     sizeof(float)*(from_rbf->nrbf*to_rbf->nrbf - 1);
42     MIGRATION *newmig;
43     #ifdef _WIN32
44     if (nprocs > 1)
45     fprintf(stderr, "%s: warning - multiprocessing not supported\n",
46     progname);
47     nprocs = 1;
48     newmig = (MIGRATION *)malloc(memlen);
49     #else
50     if (nprocs <= 1) { /* single process? */
51     newmig = (MIGRATION *)malloc(memlen);
52     } else { /* else need to share memory */
53     newmig = (MIGRATION *)mmap(NULL, memlen, PROT_READ|PROT_WRITE,
54     MAP_ANON|MAP_SHARED, -1, 0);
55     if ((void *)newmig == MAP_FAILED)
56     newmig = NULL;
57     }
58     #endif
59     if (newmig == NULL) {
60     fprintf(stderr, "%s: cannot allocate new migration\n", progname);
61     exit(1);
62     }
63     newmig->rbfv[0] = from_rbf;
64     newmig->rbfv[1] = to_rbf;
65     /* insert in edge lists */
66     newmig->enxt[0] = from_rbf->ejl;
67     from_rbf->ejl = newmig;
68     newmig->enxt[1] = to_rbf->ejl;
69     to_rbf->ejl = newmig;
70     newmig->next = mig_list; /* push onto global list */
71     return(mig_list = newmig);
72     }
73    
74     #ifdef _WIN32
75     #define await_children(n) (void)(n)
76     #define run_subprocess() 0
77     #define end_subprocess() (void)0
78     #else
79    
80     /* Wait for the specified number of child processes to complete */
81     static void
82     await_children(int n)
83     {
84     int exit_status = 0;
85    
86     if (n > nchild)
87     n = nchild;
88     while (n-- > 0) {
89     int status;
90     if (wait(&status) < 0) {
91     fprintf(stderr, "%s: missing child(ren)!\n", progname);
92     nchild = 0;
93     break;
94     }
95     --nchild;
96     if (status) { /* something wrong */
97     if ((status = WEXITSTATUS(status)))
98     exit_status = status;
99     else
100     exit_status += !exit_status;
101     fprintf(stderr, "%s: subprocess died\n", progname);
102     n = nchild; /* wait for the rest */
103     }
104     }
105     if (exit_status)
106     exit(exit_status);
107     }
108    
109     /* Start child process if multiprocessing selected */
110     static pid_t
111     run_subprocess(void)
112     {
113     int status;
114     pid_t pid;
115    
116     if (nprocs <= 1) /* any children requested? */
117     return(0);
118     await_children(nchild + 1 - nprocs); /* free up child process */
119     if ((pid = fork())) {
120     if (pid < 0) {
121     fprintf(stderr, "%s: cannot fork subprocess\n",
122     progname);
123 greg 2.6 await_children(nchild);
124 greg 2.2 exit(1);
125     }
126     ++nchild; /* subprocess started */
127     return(pid);
128     }
129     nchild = -1;
130     return(0); /* put child to work */
131     }
132    
133     /* If we are in subprocess, call exit */
134     #define end_subprocess() if (nchild < 0) _exit(0); else
135    
136     #endif /* ! _WIN32 */
137    
138 greg 2.3 /* 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 greg 2.10 float c1 = pm->prow[*(const short *)p1];
144     float c2 = pm->prow[*(const short *)p2];
145 greg 2.3
146     if (c1 > c2) return(1);
147     if (c1 < c2) return(-1);
148     return(0);
149     }
150    
151 greg 2.1 /* Compute (and allocate) migration price matrix for optimization */
152 greg 2.3 static void
153     price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, const RBFNODE *to_rbf)
154 greg 2.1 {
155     FVECT *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
156     int i, j;
157    
158 greg 2.3 pm->nrows = from_rbf->nrbf;
159     pm->ncols = to_rbf->nrbf;
160     pm->price = (float *)malloc(sizeof(float) * pm->nrows*pm->ncols);
161     pm->sord = (short *)malloc(sizeof(short) * pm->nrows*pm->ncols);
162    
163     if ((pm->price == NULL) | (pm->sord == NULL) | (vto == NULL)) {
164 greg 2.1 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 greg 2.10 short *srow;
175 greg 2.1 ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
176 greg 2.10 pm->prow = pricerow(pm,i);
177     srow = psortrow(pm,i);
178 greg 2.3 for (j = to_rbf->nrbf; j--; ) {
179 greg 2.13 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 greg 2.10 srow[j] = j;
186 greg 2.3 }
187 greg 2.10 qsort_r(srow, pm->ncols, sizeof(short), pm, &msrt_cmp);
188 greg 2.1 }
189     free(vto);
190     }
191    
192 greg 2.3 /* Free price matrix */
193     static void
194     free_routes(PRICEMAT *pm)
195 greg 2.1 {
196 greg 2.3 free(pm->price); pm->price = NULL;
197     free(pm->sord); pm->sord = NULL;
198 greg 2.1 }
199    
200     /* Compute minimum (optimistic) cost for moving the given source material */
201     static double
202 greg 2.3 min_cost(double amt2move, const double *avail, const PRICEMAT *pm, int s)
203 greg 2.1 {
204 greg 2.11 const short *srow = psortrow(pm,s);
205     const float *prow = pricerow(pm,s);
206 greg 2.1 double total_cost = 0;
207 greg 2.3 int j;
208 greg 2.1 /* move cheapest first */
209 greg 2.11 for (j = 0; (j < pm->ncols) & (amt2move > FTINY); j++) {
210     int d = srow[j];
211 greg 2.1 double amt = (amt2move < avail[d]) ? amt2move : avail[d];
212    
213 greg 2.11 total_cost += amt * prow[d];
214 greg 2.1 amt2move -= amt;
215     }
216     return(total_cost);
217     }
218    
219 greg 2.11 /* Compare entries by moving price */
220     static int
221     rmovcmp(void *b, const void *p1, const void *p2)
222     {
223     PRICEMAT *pm = (PRICEMAT *)b;
224     const short *ij1 = (const short *)p1;
225     const short *ij2 = (const short *)p2;
226     float price_diff;
227    
228     if (ij1[1] < 0) return(ij2[1] >= 0);
229     if (ij2[1] < 0) return(-1);
230     price_diff = pricerow(pm,ij1[0])[ij1[1]] - pricerow(pm,ij2[0])[ij2[1]];
231     if (price_diff > 0) return(1);
232     if (price_diff < 0) return(-1);
233     return(0);
234     }
235    
236     /* Take a step in migration by choosing reasonable bucket to transfer */
237 greg 2.1 static double
238 greg 2.11 migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, PRICEMAT *pm)
239 greg 2.1 {
240 greg 2.11 const int max2check = 100;
241 greg 2.4 const double maxamt = 1./(double)pm->ncols;
242 greg 2.12 const double minamt = maxamt*1e-4;
243 greg 2.5 double *src_cost;
244 greg 2.11 short (*rord)[2];
245 greg 2.1 struct {
246     int s, d; /* source and destination */
247     double price; /* price estimate per amount moved */
248     double amt; /* amount we can move */
249     } cur, best;
250 greg 2.11 int r2check, i, ri;
251     /*
252     * Check cheapest available routes only -- a higher adjusted
253     * destination price implies that another source is closer, so
254     * we can hold off considering more expensive options until
255     * some other (hopefully better) moves have been made.
256     */
257     /* most promising row order */
258     rord = (short (*)[2])malloc(sizeof(short)*2*pm->nrows);
259     if (rord == NULL)
260     goto memerr;
261     for (ri = pm->nrows; ri--; ) {
262     rord[ri][0] = ri;
263     rord[ri][1] = -1;
264     if (src_rem[ri] <= minamt) /* enough source material? */
265     continue;
266     for (i = 0; i < pm->ncols; i++)
267     if (dst_rem[ rord[ri][1] = psortrow(pm,ri)[i] ] > minamt)
268     break;
269     if (i >= pm->ncols) { /* moved all we can? */
270     free(rord);
271     return(.0);
272     }
273     }
274     if (pm->nrows > max2check) /* sort if too many sources */
275     qsort_r(rord, pm->nrows, sizeof(short)*2, pm, &rmovcmp);
276 greg 2.5 /* allocate cost array */
277     src_cost = (double *)malloc(sizeof(double)*pm->nrows);
278 greg 2.11 if (src_cost == NULL)
279     goto memerr;
280 greg 2.3 for (i = pm->nrows; i--; ) /* starting costs for diff. */
281     src_cost[i] = min_cost(src_rem[i], dst_rem, pm, i);
282 greg 2.1 /* find best source & dest. */
283     best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
284 greg 2.11 if ((r2check = pm->nrows) > max2check)
285     r2check = max2check; /* put a limit on search */
286     for (ri = 0; ri < r2check; ri++) { /* check each source row */
287 greg 2.1 double cost_others = 0;
288 greg 2.11 cur.s = rord[ri][0];
289     if ((cur.d = rord[ri][1]) < 0 ||
290     (cur.price = pricerow(pm,cur.s)[cur.d]) >= best.price) {
291     if (pm->nrows > max2check) break; /* sorted end */
292     continue; /* else skip this one */
293     }
294 greg 2.1 cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
295     src_rem[cur.s] : dst_rem[cur.d];
296 greg 2.11 /* don't just leave smidgen */
297     if (cur.amt > maxamt*1.02) cur.amt = maxamt;
298     dst_rem[cur.d] -= cur.amt; /* add up opportunity costs */
299 greg 2.3 for (i = pm->nrows; i--; )
300 greg 2.1 if (i != cur.s)
301 greg 2.11 cost_others += min_cost(src_rem[i], dst_rem, pm, i)
302 greg 2.1 - src_cost[i];
303     dst_rem[cur.d] += cur.amt; /* undo trial move */
304     cur.price += cost_others/cur.amt; /* adjust effective price */
305     if (cur.price < best.price) /* are we better than best? */
306 greg 2.11 best = cur;
307 greg 2.1 }
308 greg 2.11 free(src_cost); /* clean up */
309     free(rord);
310 greg 2.5 if ((best.s < 0) | (best.d < 0)) /* nothing left to move? */
311 greg 2.1 return(.0);
312 greg 2.5 /* else make the actual move */
313 greg 2.2 mtx_coef(mig,best.s,best.d) += best.amt;
314 greg 2.1 src_rem[best.s] -= best.amt;
315     dst_rem[best.d] -= best.amt;
316     return(best.amt);
317 greg 2.11 memerr:
318     fprintf(stderr, "%s: Out of memory in migration_step()\n", progname);
319     exit(1);
320 greg 2.1 }
321    
322     /* Compute and insert migration along directed edge (may fork child) */
323     static MIGRATION *
324     create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
325     {
326 greg 2.2 const double end_thresh = 5e-6;
327 greg 2.3 PRICEMAT pmtx;
328 greg 2.1 MIGRATION *newmig;
329     double *src_rem, *dst_rem;
330     double total_rem = 1., move_amt;
331 greg 2.6 int i, j;
332 greg 2.1 /* check if exists already */
333     for (newmig = from_rbf->ejl; newmig != NULL;
334     newmig = nextedge(from_rbf,newmig))
335     if (newmig->rbfv[1] == to_rbf)
336     return(NULL);
337     /* else allocate */
338 greg 2.7 #ifdef DEBUG
339 greg 2.14 fprintf(stderr, "Building path from (theta,phi) (%.1f,%.1f) ",
340 greg 2.7 get_theta180(from_rbf->invec),
341     get_phi360(from_rbf->invec));
342 greg 2.14 fprintf(stderr, "to (%.1f,%.1f) with %d x %d matrix\n",
343 greg 2.7 get_theta180(to_rbf->invec),
344     get_phi360(to_rbf->invec),
345     from_rbf->nrbf, to_rbf->nrbf);
346     #endif
347 greg 2.1 newmig = new_migration(from_rbf, to_rbf);
348     if (run_subprocess())
349     return(newmig); /* child continues */
350 greg 2.3 price_routes(&pmtx, from_rbf, to_rbf);
351 greg 2.1 src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
352     dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
353     if ((src_rem == NULL) | (dst_rem == NULL)) {
354     fprintf(stderr, "%s: Out of memory in create_migration()\n",
355     progname);
356     exit(1);
357     }
358     /* starting quantities */
359     memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
360     for (i = from_rbf->nrbf; i--; )
361     src_rem[i] = rbf_volume(&from_rbf->rbfa[i]) / from_rbf->vtotal;
362 greg 2.6 for (j = to_rbf->nrbf; j--; )
363     dst_rem[j] = rbf_volume(&to_rbf->rbfa[j]) / to_rbf->vtotal;
364    
365 greg 2.1 do { /* move a bit at a time */
366 greg 2.3 move_amt = migration_step(newmig, src_rem, dst_rem, &pmtx);
367 greg 2.1 total_rem -= move_amt;
368 greg 2.2 } while ((total_rem > end_thresh) & (move_amt > 0));
369 greg 2.6
370 greg 2.1 for (i = from_rbf->nrbf; i--; ) { /* normalize final matrix */
371 greg 2.6 double nf = rbf_volume(&from_rbf->rbfa[i]);
372 greg 2.1 if (nf <= FTINY) continue;
373     nf = from_rbf->vtotal / nf;
374     for (j = to_rbf->nrbf; j--; )
375 greg 2.6 mtx_coef(newmig,i,j) *= nf; /* row now sums to 1.0 */
376 greg 2.1 }
377     end_subprocess(); /* exit here if subprocess */
378 greg 2.3 free_routes(&pmtx); /* free working arrays */
379 greg 2.1 free(src_rem);
380     free(dst_rem);
381     return(newmig);
382     }
383    
384     /* Check if prospective vertex would create overlapping triangle */
385     static int
386     overlaps_tri(const RBFNODE *bv0, const RBFNODE *bv1, const RBFNODE *pv)
387     {
388     const MIGRATION *ej;
389     RBFNODE *vother[2];
390     int im_rev;
391     /* find shared edge in mesh */
392     for (ej = pv->ejl; ej != NULL; ej = nextedge(pv,ej)) {
393     const RBFNODE *tv = opp_rbf(pv,ej);
394     if (tv == bv0) {
395     im_rev = is_rev_tri(ej->rbfv[0]->invec,
396     ej->rbfv[1]->invec, bv1->invec);
397     break;
398     }
399     if (tv == bv1) {
400     im_rev = is_rev_tri(ej->rbfv[0]->invec,
401     ej->rbfv[1]->invec, bv0->invec);
402     break;
403     }
404     }
405     if (!get_triangles(vother, ej)) /* triangle on same side? */
406     return(0);
407     return(vother[im_rev] != NULL);
408     }
409    
410 greg 2.14 /* Find convex hull vertex to complete triangle (oriented call) */
411 greg 2.1 static RBFNODE *
412     find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rbf1)
413     {
414     FVECT vmid, vejn, vp;
415     RBFNODE *rbf, *rbfbest = NULL;
416     double dprod, area2, bestarea2 = FHUGE, bestdprod = -.5;
417    
418     VSUB(vejn, rbf1->invec, rbf0->invec);
419     VADD(vmid, rbf0->invec, rbf1->invec);
420     if (normalize(vejn) == 0 || normalize(vmid) == 0)
421     return(NULL);
422     /* XXX exhaustive search */
423     /* Find triangle with minimum rotation from perpendicular */
424     for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
425     if ((rbf == rbf0) | (rbf == rbf1))
426     continue;
427     tri_orient(vp, rbf0->invec, rbf1->invec, rbf->invec);
428     if (DOT(vp, vmid) <= FTINY)
429     continue; /* wrong orientation */
430     area2 = .25*DOT(vp,vp);
431 greg 2.14 VSUB(vp, rbf->invec, vmid);
432 greg 2.1 dprod = -DOT(vp, vejn);
433     VSUM(vp, vp, vejn, dprod); /* above guarantees non-zero */
434     dprod = DOT(vp, vmid) / VLEN(vp);
435     if (dprod <= bestdprod + FTINY*(1 - 2*(area2 < bestarea2)))
436     continue; /* found better already */
437     if (overlaps_tri(rbf0, rbf1, rbf))
438     continue; /* overlaps another triangle */
439     rbfbest = rbf;
440     bestdprod = dprod; /* new one to beat */
441     bestarea2 = area2;
442     }
443     return(rbfbest);
444     }
445    
446     /* Create new migration edge and grow mesh recursively around it */
447     static void
448     mesh_from_edge(MIGRATION *edge)
449     {
450     MIGRATION *ej0, *ej1;
451     RBFNODE *tvert[2];
452    
453     if (edge == NULL)
454     return;
455     /* triangle on either side? */
456     get_triangles(tvert, edge);
457     if (tvert[0] == NULL) { /* grow mesh on right */
458     tvert[0] = find_chull_vert(edge->rbfv[0], edge->rbfv[1]);
459     if (tvert[0] != NULL) {
460     if (tvert[0]->ord > edge->rbfv[0]->ord)
461     ej0 = create_migration(edge->rbfv[0], tvert[0]);
462     else
463     ej0 = create_migration(tvert[0], edge->rbfv[0]);
464     if (tvert[0]->ord > edge->rbfv[1]->ord)
465     ej1 = create_migration(edge->rbfv[1], tvert[0]);
466     else
467     ej1 = create_migration(tvert[0], edge->rbfv[1]);
468     mesh_from_edge(ej0);
469     mesh_from_edge(ej1);
470     }
471     } else if (tvert[1] == NULL) { /* grow mesh on left */
472     tvert[1] = find_chull_vert(edge->rbfv[1], edge->rbfv[0]);
473     if (tvert[1] != NULL) {
474     if (tvert[1]->ord > edge->rbfv[0]->ord)
475     ej0 = create_migration(edge->rbfv[0], tvert[1]);
476     else
477     ej0 = create_migration(tvert[1], edge->rbfv[0]);
478     if (tvert[1]->ord > edge->rbfv[1]->ord)
479     ej1 = create_migration(edge->rbfv[1], tvert[1]);
480     else
481     ej1 = create_migration(tvert[1], edge->rbfv[1]);
482     mesh_from_edge(ej0);
483     mesh_from_edge(ej1);
484     }
485     }
486     }
487 greg 2.15
488     /* Add normal direction if missing */
489     static void
490     check_normal_incidence(void)
491     {
492 greg 2.16 static const FVECT norm_vec = {.0, .0, 1.};
493     const int saved_nprocs = nprocs;
494     RBFNODE *near_rbf, *mir_rbf, *rbf;
495     double bestd;
496     int n;
497 greg 2.15
498     if (dsf_list == NULL)
499     return; /* XXX should be error? */
500     near_rbf = dsf_list;
501     bestd = input_orient*near_rbf->invec[2];
502     if (single_plane_incident) { /* ordered plane incidence? */
503     if (bestd >= 1.-2.*FTINY)
504     return; /* already have normal */
505     } else {
506     switch (inp_coverage) {
507     case INP_QUAD1:
508     case INP_QUAD2:
509     case INP_QUAD3:
510     case INP_QUAD4:
511     break; /* quadrilateral symmetry? */
512     default:
513     return; /* else we can interpolate */
514     }
515     for (rbf = near_rbf->next; rbf != NULL; rbf = rbf->next) {
516     const double d = input_orient*rbf->invec[2];
517     if (d >= 1.-2.*FTINY)
518     return; /* seems we have normal */
519     if (d > bestd) {
520     near_rbf = rbf;
521     bestd = d;
522     }
523     }
524     }
525     if (mig_list != NULL) { /* need to be called first */
526     fprintf(stderr, "%s: Late call to check_normal_incidence()\n",
527     progname);
528     exit(1);
529     }
530     #ifdef DEBUG
531     fprintf(stderr, "Interpolating normal incidence by mirroring (%.1f,%.1f)\n",
532     get_theta180(near_rbf->invec), get_phi360(near_rbf->invec));
533     #endif
534     /* mirror nearest incidence */
535     n = sizeof(RBFNODE) + sizeof(RBFVAL)*(near_rbf->nrbf-1);
536     mir_rbf = (RBFNODE *)malloc(n);
537     if (mir_rbf == NULL)
538     goto memerr;
539     memcpy(mir_rbf, near_rbf, n);
540     mir_rbf->ord = near_rbf->ord - 1; /* not used, I think */
541     mir_rbf->next = NULL;
542     rev_rbf_symmetry(mir_rbf, MIRROR_X|MIRROR_Y);
543     nprocs = 1; /* compute migration matrix */
544     if (mig_list != create_migration(mir_rbf, near_rbf))
545     exit(1); /* XXX should never happen! */
546 greg 2.16 /* interpolate normal dist. */
547     rbf = e_advect_rbf(mig_list, norm_vec, 2*near_rbf->nrbf);
548 greg 2.15 nprocs = saved_nprocs; /* final clean-up */
549     free(mir_rbf);
550     free(mig_list);
551     mig_list = near_rbf->ejl = NULL;
552     insert_dsf(rbf); /* insert interpolated normal */
553     return;
554     memerr:
555     fprintf(stderr, "%s: Out of memory in check_normal_incidence()\n",
556     progname);
557     exit(1);
558     }
559 greg 2.1
560     /* Build our triangle mesh from recorded RBFs */
561     void
562     build_mesh(void)
563     {
564     double best2 = M_PI*M_PI;
565     RBFNODE *shrt_edj[2];
566     RBFNODE *rbf0, *rbf1;
567 greg 2.15 /* add normal if needed */
568     check_normal_incidence();
569 greg 2.1 /* check if isotropic */
570     if (single_plane_incident) {
571     for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
572     if (rbf0->next != NULL)
573     create_migration(rbf0, rbf0->next);
574     await_children(nchild);
575     return;
576     }
577     shrt_edj[0] = shrt_edj[1] = NULL; /* start w/ shortest edge */
578     for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
579     for (rbf1 = rbf0->next; rbf1 != NULL; rbf1 = rbf1->next) {
580     double dist2 = 2. - 2.*DOT(rbf0->invec,rbf1->invec);
581     if (dist2 < best2) {
582     shrt_edj[0] = rbf0;
583     shrt_edj[1] = rbf1;
584     best2 = dist2;
585     }
586     }
587     if (shrt_edj[0] == NULL) {
588     fprintf(stderr, "%s: Cannot find shortest edge\n", progname);
589     exit(1);
590     }
591     /* build mesh from this edge */
592     if (shrt_edj[0]->ord < shrt_edj[1]->ord)
593     mesh_from_edge(create_migration(shrt_edj[0], shrt_edj[1]));
594     else
595     mesh_from_edge(create_migration(shrt_edj[1], shrt_edj[0]));
596     /* complete migrations */
597     await_children(nchild);
598     }