ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdfmesh.c
Revision: 2.7
Committed: Sat Nov 10 19:47:42 2012 UTC (11 years, 4 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.6: +12 -11 lines
Log Message:
Minor bug fixes

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: bsdfmesh.c,v 2.6 2012/11/09 02:16:29 greg Exp $";
3 #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 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 } PRICEMAT; /* sorted pricing matrix */
31
32 #define pricerow(p,i) ((p)->price + (i)*(p)->ncols)
33 #define psortrow(p,i) ((p)->sord + (i)*(p)->ncols)
34
35 /* Create a new migration holder (sharing memory for multiprocessing) */
36 static MIGRATION *
37 new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
38 {
39 size_t memlen = sizeof(MIGRATION) +
40 sizeof(float)*(from_rbf->nrbf*to_rbf->nrbf - 1);
41 MIGRATION *newmig;
42 #ifdef _WIN32
43 if (nprocs > 1)
44 fprintf(stderr, "%s: warning - multiprocessing not supported\n",
45 progname);
46 nprocs = 1;
47 newmig = (MIGRATION *)malloc(memlen);
48 #else
49 if (nprocs <= 1) { /* single process? */
50 newmig = (MIGRATION *)malloc(memlen);
51 } else { /* else need to share memory */
52 newmig = (MIGRATION *)mmap(NULL, memlen, PROT_READ|PROT_WRITE,
53 MAP_ANON|MAP_SHARED, -1, 0);
54 if ((void *)newmig == MAP_FAILED)
55 newmig = NULL;
56 }
57 #endif
58 if (newmig == NULL) {
59 fprintf(stderr, "%s: cannot allocate new migration\n", progname);
60 exit(1);
61 }
62 newmig->rbfv[0] = from_rbf;
63 newmig->rbfv[1] = to_rbf;
64 /* insert in edge lists */
65 newmig->enxt[0] = from_rbf->ejl;
66 from_rbf->ejl = newmig;
67 newmig->enxt[1] = to_rbf->ejl;
68 to_rbf->ejl = newmig;
69 newmig->next = mig_list; /* push onto global list */
70 return(mig_list = newmig);
71 }
72
73 #ifdef _WIN32
74 #define await_children(n) (void)(n)
75 #define run_subprocess() 0
76 #define end_subprocess() (void)0
77 #else
78
79 /* Wait for the specified number of child processes to complete */
80 static void
81 await_children(int n)
82 {
83 int exit_status = 0;
84
85 if (n > nchild)
86 n = nchild;
87 while (n-- > 0) {
88 int status;
89 if (wait(&status) < 0) {
90 fprintf(stderr, "%s: missing child(ren)!\n", progname);
91 nchild = 0;
92 break;
93 }
94 --nchild;
95 if (status) { /* something wrong */
96 if ((status = WEXITSTATUS(status)))
97 exit_status = status;
98 else
99 exit_status += !exit_status;
100 fprintf(stderr, "%s: subprocess died\n", progname);
101 n = nchild; /* wait for the rest */
102 }
103 }
104 if (exit_status)
105 exit(exit_status);
106 }
107
108 /* Start child process if multiprocessing selected */
109 static pid_t
110 run_subprocess(void)
111 {
112 int status;
113 pid_t pid;
114
115 if (nprocs <= 1) /* any children requested? */
116 return(0);
117 await_children(nchild + 1 - nprocs); /* free up child process */
118 if ((pid = fork())) {
119 if (pid < 0) {
120 fprintf(stderr, "%s: cannot fork subprocess\n",
121 progname);
122 await_children(nchild);
123 exit(1);
124 }
125 ++nchild; /* subprocess started */
126 return(pid);
127 }
128 nchild = -1;
129 return(0); /* put child to work */
130 }
131
132 /* If we are in subprocess, call exit */
133 #define end_subprocess() if (nchild < 0) _exit(0); else
134
135 #endif /* ! _WIN32 */
136
137 /* Comparison routine needed for sorting price row */
138 static int
139 msrt_cmp(void *b, const void *p1, const void *p2)
140 {
141 PRICEMAT *pm = (PRICEMAT *)b;
142 int ri = ((const short *)p1 - pm->sord) / pm->ncols;
143 float c1 = pricerow(pm,ri)[*(const short *)p1];
144 float c2 = pricerow(pm,ri)[*(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 */
152 static void
153 price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, const RBFNODE *to_rbf)
154 {
155 FVECT *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
156 int i, j;
157
158 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 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 ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
175 for (j = to_rbf->nrbf; j--; ) {
176 double dprod = DOT(vfrom, vto[j]);
177 pricerow(pm,i)[j] = ((dprod >= 1.) ? .0 : acos(dprod)) +
178 fabs(R2ANG(to_rbf->rbfa[j].crad) - from_ang);
179 psortrow(pm,i)[j] = j;
180 }
181 qsort_r(psortrow(pm,i), pm->ncols, sizeof(short), pm, &msrt_cmp);
182 }
183 free(vto);
184 }
185
186 /* Free price matrix */
187 static void
188 free_routes(PRICEMAT *pm)
189 {
190 free(pm->price); pm->price = NULL;
191 free(pm->sord); pm->sord = NULL;
192 }
193
194 /* Compute minimum (optimistic) cost for moving the given source material */
195 static double
196 min_cost(double amt2move, const double *avail, const PRICEMAT *pm, int s)
197 {
198 double total_cost = 0;
199 int j;
200
201 if (amt2move <= FTINY) /* pre-emptive check */
202 return(.0);
203 /* move cheapest first */
204 for (j = 0; j < pm->ncols && amt2move > FTINY; j++) {
205 int d = psortrow(pm,s)[j];
206 double amt = (amt2move < avail[d]) ? amt2move : avail[d];
207
208 total_cost += amt * pricerow(pm,s)[d];
209 amt2move -= amt;
210 }
211 return(total_cost);
212 }
213
214 /* Take a step in migration by choosing optimal bucket to transfer */
215 static double
216 migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, const PRICEMAT *pm)
217 {
218 const double maxamt = 1./(double)pm->ncols;
219 const double minamt = maxamt*5e-6;
220 double *src_cost;
221 struct {
222 int s, d; /* source and destination */
223 double price; /* price estimate per amount moved */
224 double amt; /* amount we can move */
225 } cur, best;
226 int i;
227 /* allocate cost array */
228 src_cost = (double *)malloc(sizeof(double)*pm->nrows);
229 if (src_cost == NULL) {
230 fprintf(stderr, "%s: Out of memory in migration_step()\n",
231 progname);
232 exit(1);
233 }
234 for (i = pm->nrows; i--; ) /* starting costs for diff. */
235 src_cost[i] = min_cost(src_rem[i], dst_rem, pm, i);
236
237 /* find best source & dest. */
238 best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
239 for (cur.s = pm->nrows; cur.s--; ) {
240 double cost_others = 0;
241
242 if (src_rem[cur.s] <= minamt)
243 continue;
244 /* examine cheapest dest. */
245 for (i = 0; i < pm->ncols; i++)
246 if (dst_rem[ cur.d = psortrow(pm,cur.s)[i] ] > minamt)
247 break;
248 if (i >= pm->ncols)
249 break;
250 if ((cur.price = pricerow(pm,cur.s)[cur.d]) >= best.price)
251 continue; /* no point checking further */
252 cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
253 src_rem[cur.s] : dst_rem[cur.d];
254 if (cur.amt > maxamt) cur.amt = maxamt;
255 dst_rem[cur.d] -= cur.amt; /* add up differential costs */
256 for (i = pm->nrows; i--; )
257 if (i != cur.s)
258 cost_others += min_cost(src_rem[i], dst_rem, pm, i)
259 - src_cost[i];
260 dst_rem[cur.d] += cur.amt; /* undo trial move */
261 cur.price += cost_others/cur.amt; /* adjust effective price */
262 if (cur.price < best.price) /* are we better than best? */
263 best = cur;
264 }
265 free(src_cost); /* finish up */
266
267 if ((best.s < 0) | (best.d < 0)) /* nothing left to move? */
268 return(.0);
269 /* else make the actual move */
270 mtx_coef(mig,best.s,best.d) += best.amt;
271 src_rem[best.s] -= best.amt;
272 dst_rem[best.d] -= best.amt;
273 return(best.amt);
274 }
275
276 /* Compute and insert migration along directed edge (may fork child) */
277 static MIGRATION *
278 create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
279 {
280 const double end_thresh = 5e-6;
281 PRICEMAT pmtx;
282 MIGRATION *newmig;
283 double *src_rem, *dst_rem;
284 double total_rem = 1., move_amt;
285 int i, j;
286 /* check if exists already */
287 for (newmig = from_rbf->ejl; newmig != NULL;
288 newmig = nextedge(from_rbf,newmig))
289 if (newmig->rbfv[1] == to_rbf)
290 return(NULL);
291 /* else allocate */
292 #ifdef DEBUG
293 fprintf(stderr, "Building path from (theta,phi) (%.0f,%.0f) ",
294 get_theta180(from_rbf->invec),
295 get_phi360(from_rbf->invec));
296 fprintf(stderr, "to (%.0f,%.0f) with %d x %d matrix\n",
297 get_theta180(to_rbf->invec),
298 get_phi360(to_rbf->invec),
299 from_rbf->nrbf, to_rbf->nrbf);
300 #endif
301 newmig = new_migration(from_rbf, to_rbf);
302 if (run_subprocess())
303 return(newmig); /* child continues */
304 price_routes(&pmtx, from_rbf, to_rbf);
305 src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
306 dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
307 if ((src_rem == NULL) | (dst_rem == NULL)) {
308 fprintf(stderr, "%s: Out of memory in create_migration()\n",
309 progname);
310 exit(1);
311 }
312 /* starting quantities */
313 memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
314 for (i = from_rbf->nrbf; i--; )
315 src_rem[i] = rbf_volume(&from_rbf->rbfa[i]) / from_rbf->vtotal;
316 for (j = to_rbf->nrbf; j--; )
317 dst_rem[j] = rbf_volume(&to_rbf->rbfa[j]) / to_rbf->vtotal;
318
319 do { /* move a bit at a time */
320 move_amt = migration_step(newmig, src_rem, dst_rem, &pmtx);
321 total_rem -= move_amt;
322 } while ((total_rem > end_thresh) & (move_amt > 0));
323
324 for (i = from_rbf->nrbf; i--; ) { /* normalize final matrix */
325 double nf = rbf_volume(&from_rbf->rbfa[i]);
326 if (nf <= FTINY) continue;
327 nf = from_rbf->vtotal / nf;
328 for (j = to_rbf->nrbf; j--; )
329 mtx_coef(newmig,i,j) *= nf; /* row now sums to 1.0 */
330 }
331 end_subprocess(); /* exit here if subprocess */
332 free_routes(&pmtx); /* free working arrays */
333 free(src_rem);
334 free(dst_rem);
335 return(newmig);
336 }
337
338 /* Check if prospective vertex would create overlapping triangle */
339 static int
340 overlaps_tri(const RBFNODE *bv0, const RBFNODE *bv1, const RBFNODE *pv)
341 {
342 const MIGRATION *ej;
343 RBFNODE *vother[2];
344 int im_rev;
345 /* find shared edge in mesh */
346 for (ej = pv->ejl; ej != NULL; ej = nextedge(pv,ej)) {
347 const RBFNODE *tv = opp_rbf(pv,ej);
348 if (tv == bv0) {
349 im_rev = is_rev_tri(ej->rbfv[0]->invec,
350 ej->rbfv[1]->invec, bv1->invec);
351 break;
352 }
353 if (tv == bv1) {
354 im_rev = is_rev_tri(ej->rbfv[0]->invec,
355 ej->rbfv[1]->invec, bv0->invec);
356 break;
357 }
358 }
359 if (!get_triangles(vother, ej)) /* triangle on same side? */
360 return(0);
361 return(vother[im_rev] != NULL);
362 }
363
364 /* Find context hull vertex to complete triangle (oriented call) */
365 static RBFNODE *
366 find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rbf1)
367 {
368 FVECT vmid, vejn, vp;
369 RBFNODE *rbf, *rbfbest = NULL;
370 double dprod, area2, bestarea2 = FHUGE, bestdprod = -.5;
371
372 VSUB(vejn, rbf1->invec, rbf0->invec);
373 VADD(vmid, rbf0->invec, rbf1->invec);
374 if (normalize(vejn) == 0 || normalize(vmid) == 0)
375 return(NULL);
376 /* XXX exhaustive search */
377 /* Find triangle with minimum rotation from perpendicular */
378 for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
379 if ((rbf == rbf0) | (rbf == rbf1))
380 continue;
381 tri_orient(vp, rbf0->invec, rbf1->invec, rbf->invec);
382 if (DOT(vp, vmid) <= FTINY)
383 continue; /* wrong orientation */
384 area2 = .25*DOT(vp,vp);
385 VSUB(vp, rbf->invec, rbf0->invec);
386 dprod = -DOT(vp, vejn);
387 VSUM(vp, vp, vejn, dprod); /* above guarantees non-zero */
388 dprod = DOT(vp, vmid) / VLEN(vp);
389 if (dprod <= bestdprod + FTINY*(1 - 2*(area2 < bestarea2)))
390 continue; /* found better already */
391 if (overlaps_tri(rbf0, rbf1, rbf))
392 continue; /* overlaps another triangle */
393 rbfbest = rbf;
394 bestdprod = dprod; /* new one to beat */
395 bestarea2 = area2;
396 }
397 return(rbfbest);
398 }
399
400 /* Create new migration edge and grow mesh recursively around it */
401 static void
402 mesh_from_edge(MIGRATION *edge)
403 {
404 MIGRATION *ej0, *ej1;
405 RBFNODE *tvert[2];
406
407 if (edge == NULL)
408 return;
409 /* triangle on either side? */
410 get_triangles(tvert, edge);
411 if (tvert[0] == NULL) { /* grow mesh on right */
412 tvert[0] = find_chull_vert(edge->rbfv[0], edge->rbfv[1]);
413 if (tvert[0] != NULL) {
414 if (tvert[0]->ord > edge->rbfv[0]->ord)
415 ej0 = create_migration(edge->rbfv[0], tvert[0]);
416 else
417 ej0 = create_migration(tvert[0], edge->rbfv[0]);
418 if (tvert[0]->ord > edge->rbfv[1]->ord)
419 ej1 = create_migration(edge->rbfv[1], tvert[0]);
420 else
421 ej1 = create_migration(tvert[0], edge->rbfv[1]);
422 mesh_from_edge(ej0);
423 mesh_from_edge(ej1);
424 }
425 } else if (tvert[1] == NULL) { /* grow mesh on left */
426 tvert[1] = find_chull_vert(edge->rbfv[1], edge->rbfv[0]);
427 if (tvert[1] != NULL) {
428 if (tvert[1]->ord > edge->rbfv[0]->ord)
429 ej0 = create_migration(edge->rbfv[0], tvert[1]);
430 else
431 ej0 = create_migration(tvert[1], edge->rbfv[0]);
432 if (tvert[1]->ord > edge->rbfv[1]->ord)
433 ej1 = create_migration(edge->rbfv[1], tvert[1]);
434 else
435 ej1 = create_migration(tvert[1], edge->rbfv[1]);
436 mesh_from_edge(ej0);
437 mesh_from_edge(ej1);
438 }
439 }
440 }
441
442 /* Build our triangle mesh from recorded RBFs */
443 void
444 build_mesh(void)
445 {
446 double best2 = M_PI*M_PI;
447 RBFNODE *shrt_edj[2];
448 RBFNODE *rbf0, *rbf1;
449 /* check if isotropic */
450 if (single_plane_incident) {
451 for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
452 if (rbf0->next != NULL)
453 create_migration(rbf0, rbf0->next);
454 await_children(nchild);
455 return;
456 }
457 shrt_edj[0] = shrt_edj[1] = NULL; /* start w/ shortest edge */
458 for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
459 for (rbf1 = rbf0->next; rbf1 != NULL; rbf1 = rbf1->next) {
460 double dist2 = 2. - 2.*DOT(rbf0->invec,rbf1->invec);
461 if (dist2 < best2) {
462 shrt_edj[0] = rbf0;
463 shrt_edj[1] = rbf1;
464 best2 = dist2;
465 }
466 }
467 if (shrt_edj[0] == NULL) {
468 fprintf(stderr, "%s: Cannot find shortest edge\n", progname);
469 exit(1);
470 }
471 /* build mesh from this edge */
472 if (shrt_edj[0]->ord < shrt_edj[1]->ord)
473 mesh_from_edge(create_migration(shrt_edj[0], shrt_edj[1]));
474 else
475 mesh_from_edge(create_migration(shrt_edj[1], shrt_edj[0]));
476 /* complete migrations */
477 await_children(nchild);
478 }