ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdfmesh.c
Revision: 2.8
Committed: Wed Mar 20 01:00:22 2013 UTC (11 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.7: +22 -1 lines
Log Message:
Added minimum BSDF to output where we have no data.

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.8 static const char RCSid[] = "$Id: bsdfmesh.c,v 2.7 2012/11/10 19:47:42 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     } 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 greg 2.2 /* 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 greg 2.6 await_children(nchild);
123 greg 2.2 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 greg 2.3 /* 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 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     ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
175 greg 2.3 for (j = to_rbf->nrbf; j--; ) {
176 greg 2.7 double dprod = DOT(vfrom, vto[j]);
177     pricerow(pm,i)[j] = ((dprod >= 1.) ? .0 : acos(dprod)) +
178 greg 2.1 fabs(R2ANG(to_rbf->rbfa[j].crad) - from_ang);
179 greg 2.3 psortrow(pm,i)[j] = j;
180     }
181     qsort_r(psortrow(pm,i), pm->ncols, sizeof(short), pm, &msrt_cmp);
182 greg 2.1 }
183     free(vto);
184     }
185    
186 greg 2.3 /* Free price matrix */
187     static void
188     free_routes(PRICEMAT *pm)
189 greg 2.1 {
190 greg 2.3 free(pm->price); pm->price = NULL;
191     free(pm->sord); pm->sord = NULL;
192 greg 2.1 }
193    
194     /* Compute minimum (optimistic) cost for moving the given source material */
195     static double
196 greg 2.3 min_cost(double amt2move, const double *avail, const PRICEMAT *pm, int s)
197 greg 2.1 {
198     double total_cost = 0;
199 greg 2.3 int j;
200 greg 2.1
201     if (amt2move <= FTINY) /* pre-emptive check */
202 greg 2.5 return(.0);
203 greg 2.1 /* move cheapest first */
204 greg 2.3 for (j = 0; j < pm->ncols && amt2move > FTINY; j++) {
205     int d = psortrow(pm,s)[j];
206 greg 2.1 double amt = (amt2move < avail[d]) ? amt2move : avail[d];
207    
208 greg 2.3 total_cost += amt * pricerow(pm,s)[d];
209 greg 2.1 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 greg 2.3 migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, const PRICEMAT *pm)
217 greg 2.1 {
218 greg 2.4 const double maxamt = 1./(double)pm->ncols;
219 greg 2.2 const double minamt = maxamt*5e-6;
220 greg 2.5 double *src_cost;
221 greg 2.1 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 greg 2.5 /* 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 greg 2.1 }
234 greg 2.3 for (i = pm->nrows; i--; ) /* starting costs for diff. */
235     src_cost[i] = min_cost(src_rem[i], dst_rem, pm, i);
236 greg 2.1
237     /* find best source & dest. */
238     best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
239 greg 2.3 for (cur.s = pm->nrows; cur.s--; ) {
240 greg 2.1 double cost_others = 0;
241 greg 2.5
242 greg 2.2 if (src_rem[cur.s] <= minamt)
243 greg 2.1 continue;
244 greg 2.4 /* examine cheapest dest. */
245     for (i = 0; i < pm->ncols; i++)
246 greg 2.5 if (dst_rem[ cur.d = psortrow(pm,cur.s)[i] ] > minamt)
247 greg 2.4 break;
248     if (i >= pm->ncols)
249 greg 2.5 break;
250 greg 2.4 if ((cur.price = pricerow(pm,cur.s)[cur.d]) >= best.price)
251 greg 2.5 continue; /* no point checking further */
252 greg 2.1 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 greg 2.3 for (i = pm->nrows; i--; )
257 greg 2.1 if (i != cur.s)
258 greg 2.3 cost_others += min_cost(src_rem[i], dst_rem, pm, i)
259 greg 2.1 - 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 greg 2.5 free(src_cost); /* finish up */
266    
267     if ((best.s < 0) | (best.d < 0)) /* nothing left to move? */
268 greg 2.1 return(.0);
269 greg 2.5 /* else make the actual move */
270 greg 2.2 mtx_coef(mig,best.s,best.d) += best.amt;
271 greg 2.1 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 greg 2.2 const double end_thresh = 5e-6;
281 greg 2.3 PRICEMAT pmtx;
282 greg 2.1 MIGRATION *newmig;
283     double *src_rem, *dst_rem;
284     double total_rem = 1., move_amt;
285 greg 2.6 int i, j;
286 greg 2.1 /* 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 greg 2.7 #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 greg 2.1 newmig = new_migration(from_rbf, to_rbf);
302     if (run_subprocess())
303     return(newmig); /* child continues */
304 greg 2.3 price_routes(&pmtx, from_rbf, to_rbf);
305 greg 2.1 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 greg 2.6 for (j = to_rbf->nrbf; j--; )
317     dst_rem[j] = rbf_volume(&to_rbf->rbfa[j]) / to_rbf->vtotal;
318    
319 greg 2.1 do { /* move a bit at a time */
320 greg 2.3 move_amt = migration_step(newmig, src_rem, dst_rem, &pmtx);
321 greg 2.1 total_rem -= move_amt;
322 greg 2.2 } while ((total_rem > end_thresh) & (move_amt > 0));
323 greg 2.6
324 greg 2.1 for (i = from_rbf->nrbf; i--; ) { /* normalize final matrix */
325 greg 2.6 double nf = rbf_volume(&from_rbf->rbfa[i]);
326 greg 2.1 if (nf <= FTINY) continue;
327     nf = from_rbf->vtotal / nf;
328     for (j = to_rbf->nrbf; j--; )
329 greg 2.6 mtx_coef(newmig,i,j) *= nf; /* row now sums to 1.0 */
330 greg 2.1 }
331     end_subprocess(); /* exit here if subprocess */
332 greg 2.3 free_routes(&pmtx); /* free working arrays */
333 greg 2.1 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 greg 2.8
442     /* Compute minimum BSDF from histogram and clear it */
443     static void
444     comp_bsdf_min()
445     {
446     int cnt;
447     int i, target;
448    
449     cnt = 0;
450     for (i = HISTLEN; i--; )
451     cnt += bsdf_hist[i];
452    
453     target = cnt/100; /* ignore bottom 1% */
454     cnt = 0;
455     for (i = 0; cnt <= target; i++)
456     cnt += bsdf_hist[i];
457     bsdf_min = histval(i-1);
458     memset(bsdf_hist, 0, sizeof(bsdf_hist));
459     }
460 greg 2.1
461     /* Build our triangle mesh from recorded RBFs */
462     void
463     build_mesh(void)
464     {
465     double best2 = M_PI*M_PI;
466     RBFNODE *shrt_edj[2];
467     RBFNODE *rbf0, *rbf1;
468     /* check if isotropic */
469     if (single_plane_incident) {
470     for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
471     if (rbf0->next != NULL)
472     create_migration(rbf0, rbf0->next);
473     await_children(nchild);
474     return;
475     }
476     shrt_edj[0] = shrt_edj[1] = NULL; /* start w/ shortest edge */
477     for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
478     for (rbf1 = rbf0->next; rbf1 != NULL; rbf1 = rbf1->next) {
479     double dist2 = 2. - 2.*DOT(rbf0->invec,rbf1->invec);
480     if (dist2 < best2) {
481     shrt_edj[0] = rbf0;
482     shrt_edj[1] = rbf1;
483     best2 = dist2;
484     }
485     }
486     if (shrt_edj[0] == NULL) {
487     fprintf(stderr, "%s: Cannot find shortest edge\n", progname);
488     exit(1);
489     }
490     /* build mesh from this edge */
491     if (shrt_edj[0]->ord < shrt_edj[1]->ord)
492     mesh_from_edge(create_migration(shrt_edj[0], shrt_edj[1]));
493     else
494     mesh_from_edge(create_migration(shrt_edj[1], shrt_edj[0]));
495 greg 2.8 /* compute minimum BSDF */
496     comp_bsdf_min();
497 greg 2.1 /* complete migrations */
498     await_children(nchild);
499     }