ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdfmesh.c
Revision: 2.33
Committed: Fri Aug 22 05:38:44 2014 UTC (9 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R2P2, rad5R0, rad4R2P1
Changes since 2.32: +4 -4 lines
Log Message:
Set minimum cosine to 0.02 to avoid blowing-up values near grazing

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.33 static const char RCSid[] = "$Id: bsdfmesh.c,v 2.32 2014/08/21 19:31:23 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.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 greg 2.33 sum += eval_rbfrep(rbf, tvec) * COSF(tvec[2]);
55 greg 2.30 ++n;
56     }
57     if (n < 2) /* should never happen! */
58     return(sum);
59     return(sum/(double)n);
60     }
61    
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.03;
67 greg 2.33 const double DSFtarget = 0.60653066 * eval_rbfrep(rbf,outvec) *
68     COSF(outvec[2]);
69 greg 2.30 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 greg 2.32 /* Newton's method (sort of) */
76     do {
77 greg 2.30 double test_rad = interp_rad;
78 greg 2.32 double DSFtest;
79     if (test_rad >= outside_rad)
80     return(test_rad);
81     if (test_rad <= inside_rad)
82     return(test_rad*(test_rad>0));
83     DSFtest = eval_DSFsurround(rbf, outvec, test_rad);
84     if (DSFtest > DSFtarget) {
85 greg 2.30 inside_rad = test_rad;
86     DSFinside = DSFtest;
87     } else {
88     outside_rad = test_rad;
89     DSFoutside = DSFtest;
90     }
91 greg 2.32 if (DSFoutside >= DSFinside)
92     return(test_rad);
93     } while (outside_rad-inside_rad > rad_epsilon);
94 greg 2.30 return(interp_rad);
95     #undef interp_rad
96     }
97    
98     /* Compute average BSDF peak from current DSF's */
99     static void
100     comp_bsdf_spec(void)
101     {
102     double peak_sum = 0;
103     double rad_sum = 0;
104     int n = 0;
105     RBFNODE *rbf;
106     FVECT sdv;
107    
108     if (dsf_list == NULL) {
109     bsdf_spec_peak = 0;
110 greg 2.31 bsdf_spec_rad = 0;
111 greg 2.30 return;
112     }
113     for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
114     sdv[0] = -rbf->invec[0];
115     sdv[1] = -rbf->invec[1];
116     sdv[2] = rbf->invec[2]*(2*(input_orient==output_orient) - 1);
117     peak_sum += eval_rbfrep(rbf, sdv);
118     rad_sum += est_DSFrad(rbf, sdv);
119     ++n;
120     }
121     bsdf_spec_peak = peak_sum/(double)n;
122 greg 2.31 bsdf_spec_rad = rad_sum/(double)n;
123 greg 2.30 }
124    
125 greg 2.2 /* Create a new migration holder (sharing memory for multiprocessing) */
126     static MIGRATION *
127     new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
128     {
129     size_t memlen = sizeof(MIGRATION) +
130     sizeof(float)*(from_rbf->nrbf*to_rbf->nrbf - 1);
131     MIGRATION *newmig;
132     #ifdef _WIN32
133     if (nprocs > 1)
134     fprintf(stderr, "%s: warning - multiprocessing not supported\n",
135     progname);
136     nprocs = 1;
137     newmig = (MIGRATION *)malloc(memlen);
138     #else
139     if (nprocs <= 1) { /* single process? */
140     newmig = (MIGRATION *)malloc(memlen);
141     } else { /* else need to share memory */
142     newmig = (MIGRATION *)mmap(NULL, memlen, PROT_READ|PROT_WRITE,
143     MAP_ANON|MAP_SHARED, -1, 0);
144     if ((void *)newmig == MAP_FAILED)
145     newmig = NULL;
146     }
147     #endif
148     if (newmig == NULL) {
149     fprintf(stderr, "%s: cannot allocate new migration\n", progname);
150     exit(1);
151     }
152     newmig->rbfv[0] = from_rbf;
153     newmig->rbfv[1] = to_rbf;
154     /* insert in edge lists */
155     newmig->enxt[0] = from_rbf->ejl;
156     from_rbf->ejl = newmig;
157     newmig->enxt[1] = to_rbf->ejl;
158     to_rbf->ejl = newmig;
159     newmig->next = mig_list; /* push onto global list */
160     return(mig_list = newmig);
161     }
162    
163     #ifdef _WIN32
164     #define await_children(n) (void)(n)
165     #define run_subprocess() 0
166     #define end_subprocess() (void)0
167     #else
168    
169     /* Wait for the specified number of child processes to complete */
170     static void
171     await_children(int n)
172     {
173     int exit_status = 0;
174    
175     if (n > nchild)
176     n = nchild;
177     while (n-- > 0) {
178     int status;
179     if (wait(&status) < 0) {
180     fprintf(stderr, "%s: missing child(ren)!\n", progname);
181     nchild = 0;
182     break;
183     }
184     --nchild;
185     if (status) { /* something wrong */
186     if ((status = WEXITSTATUS(status)))
187     exit_status = status;
188     else
189     exit_status += !exit_status;
190     fprintf(stderr, "%s: subprocess died\n", progname);
191     n = nchild; /* wait for the rest */
192     }
193     }
194     if (exit_status)
195     exit(exit_status);
196     }
197    
198     /* Start child process if multiprocessing selected */
199     static pid_t
200     run_subprocess(void)
201     {
202     int status;
203     pid_t pid;
204    
205     if (nprocs <= 1) /* any children requested? */
206     return(0);
207     await_children(nchild + 1 - nprocs); /* free up child process */
208     if ((pid = fork())) {
209     if (pid < 0) {
210     fprintf(stderr, "%s: cannot fork subprocess\n",
211     progname);
212 greg 2.6 await_children(nchild);
213 greg 2.2 exit(1);
214     }
215     ++nchild; /* subprocess started */
216     return(pid);
217     }
218     nchild = -1;
219     return(0); /* put child to work */
220     }
221    
222     /* If we are in subprocess, call exit */
223     #define end_subprocess() if (nchild < 0) _exit(0); else
224    
225     #endif /* ! _WIN32 */
226    
227 greg 2.19 /* Compute normalized distribution scattering functions for comparison */
228     static void
229     compute_nDSFs(const RBFNODE *rbf0, const RBFNODE *rbf1)
230     {
231     const double nf0 = (GRIDRES*GRIDRES) / rbf0->vtotal;
232     const double nf1 = (GRIDRES*GRIDRES) / rbf1->vtotal;
233     int x, y;
234     FVECT dv;
235    
236     for (x = GRIDRES; x--; )
237     for (y = GRIDRES; y--; ) {
238 greg 2.20 ovec_from_pos(dv, x, y); /* cube root (brightness) */
239     dsf_grid[x][y].val[0] = pow(nf0*eval_rbfrep(rbf0, dv), .3333);
240     dsf_grid[x][y].val[1] = pow(nf1*eval_rbfrep(rbf1, dv), .3333);
241 greg 2.19 }
242     }
243    
244     /* Compute neighborhood distance-squared (dissimilarity) */
245     static double
246     neighborhood_dist2(int x0, int y0, int x1, int y1)
247     {
248     int rad = GRIDRES>>5;
249     double sum2 = 0.;
250     double d;
251     int p[4];
252     int i, j;
253     /* check radius */
254     p[0] = x0; p[1] = y0; p[2] = x1; p[3] = y1;
255     for (i = 4; i--; ) {
256     if (p[i] < rad) rad = p[i];
257     if (GRIDRES-1-p[i] < rad) rad = GRIDRES-1-p[i];
258     }
259     for (i = -rad; i <= rad; i++)
260     for (j = -rad; j <= rad; j++) {
261     d = dsf_grid[x0+i][y0+j].val[0] -
262     dsf_grid[x1+i][y1+j].val[1];
263     sum2 += d*d;
264     }
265     return(sum2 / (4*rad*(rad+1) + 1));
266     }
267    
268 greg 2.27 /* Compute distance between two RBF lobes */
269     double
270     lobe_distance(RBFVAL *rbf1, RBFVAL *rbf2)
271     {
272     FVECT vfrom, vto;
273     double d, res;
274     /* quadratic cost function */
275     ovec_from_pos(vfrom, rbf1->gx, rbf1->gy);
276     ovec_from_pos(vto, rbf2->gx, rbf2->gy);
277     d = Acos(DOT(vfrom, vto));
278     res = d*d;
279     d = R2ANG(rbf2->crad) - R2ANG(rbf1->crad);
280     res += d*d;
281     /* neighborhood difference */
282     res += NEIGH_FACT2 * neighborhood_dist2( rbf1->gx, rbf1->gy,
283     rbf2->gx, rbf2->gy );
284     return(res);
285 greg 2.1 }
286    
287 greg 2.26
288 greg 2.1 /* Compute and insert migration along directed edge (may fork child) */
289     static MIGRATION *
290     create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
291     {
292     MIGRATION *newmig;
293 greg 2.6 int i, j;
294 greg 2.1 /* check if exists already */
295     for (newmig = from_rbf->ejl; newmig != NULL;
296     newmig = nextedge(from_rbf,newmig))
297     if (newmig->rbfv[1] == to_rbf)
298     return(NULL);
299     /* else allocate */
300 greg 2.7 #ifdef DEBUG
301 greg 2.14 fprintf(stderr, "Building path from (theta,phi) (%.1f,%.1f) ",
302 greg 2.7 get_theta180(from_rbf->invec),
303     get_phi360(from_rbf->invec));
304 greg 2.14 fprintf(stderr, "to (%.1f,%.1f) with %d x %d matrix\n",
305 greg 2.7 get_theta180(to_rbf->invec),
306     get_phi360(to_rbf->invec),
307     from_rbf->nrbf, to_rbf->nrbf);
308     #endif
309 greg 2.1 newmig = new_migration(from_rbf, to_rbf);
310     if (run_subprocess())
311     return(newmig); /* child continues */
312 greg 2.27
313     /* compute transport plan */
314     compute_nDSFs(from_rbf, to_rbf);
315     plan_transport(newmig);
316 greg 2.6
317 greg 2.1 for (i = from_rbf->nrbf; i--; ) { /* normalize final matrix */
318 greg 2.6 double nf = rbf_volume(&from_rbf->rbfa[i]);
319 greg 2.1 if (nf <= FTINY) continue;
320     nf = from_rbf->vtotal / nf;
321     for (j = to_rbf->nrbf; j--; )
322 greg 2.6 mtx_coef(newmig,i,j) *= nf; /* row now sums to 1.0 */
323 greg 2.1 }
324     end_subprocess(); /* exit here if subprocess */
325     return(newmig);
326     }
327    
328     /* Check if prospective vertex would create overlapping triangle */
329     static int
330     overlaps_tri(const RBFNODE *bv0, const RBFNODE *bv1, const RBFNODE *pv)
331     {
332     const MIGRATION *ej;
333     RBFNODE *vother[2];
334     int im_rev;
335     /* find shared edge in mesh */
336     for (ej = pv->ejl; ej != NULL; ej = nextedge(pv,ej)) {
337     const RBFNODE *tv = opp_rbf(pv,ej);
338     if (tv == bv0) {
339     im_rev = is_rev_tri(ej->rbfv[0]->invec,
340     ej->rbfv[1]->invec, bv1->invec);
341     break;
342     }
343     if (tv == bv1) {
344     im_rev = is_rev_tri(ej->rbfv[0]->invec,
345     ej->rbfv[1]->invec, bv0->invec);
346     break;
347     }
348     }
349     if (!get_triangles(vother, ej)) /* triangle on same side? */
350     return(0);
351     return(vother[im_rev] != NULL);
352     }
353    
354 greg 2.14 /* Find convex hull vertex to complete triangle (oriented call) */
355 greg 2.1 static RBFNODE *
356     find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rbf1)
357     {
358     FVECT vmid, vejn, vp;
359     RBFNODE *rbf, *rbfbest = NULL;
360     double dprod, area2, bestarea2 = FHUGE, bestdprod = -.5;
361    
362     VSUB(vejn, rbf1->invec, rbf0->invec);
363     VADD(vmid, rbf0->invec, rbf1->invec);
364     if (normalize(vejn) == 0 || normalize(vmid) == 0)
365     return(NULL);
366     /* XXX exhaustive search */
367     /* Find triangle with minimum rotation from perpendicular */
368     for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
369     if ((rbf == rbf0) | (rbf == rbf1))
370     continue;
371     tri_orient(vp, rbf0->invec, rbf1->invec, rbf->invec);
372     if (DOT(vp, vmid) <= FTINY)
373     continue; /* wrong orientation */
374     area2 = .25*DOT(vp,vp);
375 greg 2.14 VSUB(vp, rbf->invec, vmid);
376 greg 2.1 dprod = -DOT(vp, vejn);
377     VSUM(vp, vp, vejn, dprod); /* above guarantees non-zero */
378     dprod = DOT(vp, vmid) / VLEN(vp);
379     if (dprod <= bestdprod + FTINY*(1 - 2*(area2 < bestarea2)))
380     continue; /* found better already */
381     if (overlaps_tri(rbf0, rbf1, rbf))
382     continue; /* overlaps another triangle */
383     rbfbest = rbf;
384     bestdprod = dprod; /* new one to beat */
385     bestarea2 = area2;
386     }
387     return(rbfbest);
388     }
389    
390     /* Create new migration edge and grow mesh recursively around it */
391     static void
392     mesh_from_edge(MIGRATION *edge)
393     {
394     MIGRATION *ej0, *ej1;
395     RBFNODE *tvert[2];
396    
397     if (edge == NULL)
398     return;
399     /* triangle on either side? */
400     get_triangles(tvert, edge);
401     if (tvert[0] == NULL) { /* grow mesh on right */
402     tvert[0] = find_chull_vert(edge->rbfv[0], edge->rbfv[1]);
403     if (tvert[0] != NULL) {
404     if (tvert[0]->ord > edge->rbfv[0]->ord)
405     ej0 = create_migration(edge->rbfv[0], tvert[0]);
406     else
407     ej0 = create_migration(tvert[0], edge->rbfv[0]);
408     if (tvert[0]->ord > edge->rbfv[1]->ord)
409     ej1 = create_migration(edge->rbfv[1], tvert[0]);
410     else
411     ej1 = create_migration(tvert[0], edge->rbfv[1]);
412     mesh_from_edge(ej0);
413     mesh_from_edge(ej1);
414 greg 2.28 return;
415 greg 2.1 }
416 greg 2.28 }
417     if (tvert[1] == NULL) { /* grow mesh on left */
418 greg 2.1 tvert[1] = find_chull_vert(edge->rbfv[1], edge->rbfv[0]);
419     if (tvert[1] != NULL) {
420     if (tvert[1]->ord > edge->rbfv[0]->ord)
421     ej0 = create_migration(edge->rbfv[0], tvert[1]);
422     else
423     ej0 = create_migration(tvert[1], edge->rbfv[0]);
424     if (tvert[1]->ord > edge->rbfv[1]->ord)
425     ej1 = create_migration(edge->rbfv[1], tvert[1]);
426     else
427     ej1 = create_migration(tvert[1], edge->rbfv[1]);
428     mesh_from_edge(ej0);
429     mesh_from_edge(ej1);
430     }
431     }
432     }
433 greg 2.15
434     /* Add normal direction if missing */
435     static void
436     check_normal_incidence(void)
437     {
438 greg 2.25 static FVECT norm_vec = {.0, .0, 1.};
439 greg 2.16 const int saved_nprocs = nprocs;
440     RBFNODE *near_rbf, *mir_rbf, *rbf;
441     double bestd;
442     int n;
443 greg 2.15
444     if (dsf_list == NULL)
445     return; /* XXX should be error? */
446     near_rbf = dsf_list;
447     bestd = input_orient*near_rbf->invec[2];
448     if (single_plane_incident) { /* ordered plane incidence? */
449     if (bestd >= 1.-2.*FTINY)
450     return; /* already have normal */
451     } else {
452     switch (inp_coverage) {
453     case INP_QUAD1:
454     case INP_QUAD2:
455     case INP_QUAD3:
456     case INP_QUAD4:
457     break; /* quadrilateral symmetry? */
458     default:
459     return; /* else we can interpolate */
460     }
461     for (rbf = near_rbf->next; rbf != NULL; rbf = rbf->next) {
462     const double d = input_orient*rbf->invec[2];
463     if (d >= 1.-2.*FTINY)
464     return; /* seems we have normal */
465     if (d > bestd) {
466     near_rbf = rbf;
467     bestd = d;
468     }
469     }
470     }
471     if (mig_list != NULL) { /* need to be called first */
472     fprintf(stderr, "%s: Late call to check_normal_incidence()\n",
473     progname);
474     exit(1);
475     }
476     #ifdef DEBUG
477     fprintf(stderr, "Interpolating normal incidence by mirroring (%.1f,%.1f)\n",
478     get_theta180(near_rbf->invec), get_phi360(near_rbf->invec));
479     #endif
480     /* mirror nearest incidence */
481     n = sizeof(RBFNODE) + sizeof(RBFVAL)*(near_rbf->nrbf-1);
482     mir_rbf = (RBFNODE *)malloc(n);
483     if (mir_rbf == NULL)
484     goto memerr;
485     memcpy(mir_rbf, near_rbf, n);
486     mir_rbf->ord = near_rbf->ord - 1; /* not used, I think */
487     mir_rbf->next = NULL;
488 greg 2.22 mir_rbf->ejl = NULL;
489 greg 2.15 rev_rbf_symmetry(mir_rbf, MIRROR_X|MIRROR_Y);
490     nprocs = 1; /* compute migration matrix */
491 greg 2.22 if (create_migration(mir_rbf, near_rbf) == NULL)
492 greg 2.15 exit(1); /* XXX should never happen! */
493 greg 2.25 norm_vec[2] = input_orient; /* interpolate normal dist. */
494 greg 2.29 rbf = e_advect_rbf(mig_list, norm_vec, 0);
495 greg 2.15 nprocs = saved_nprocs; /* final clean-up */
496     free(mir_rbf);
497     free(mig_list);
498     mig_list = near_rbf->ejl = NULL;
499     insert_dsf(rbf); /* insert interpolated normal */
500     return;
501     memerr:
502     fprintf(stderr, "%s: Out of memory in check_normal_incidence()\n",
503     progname);
504     exit(1);
505     }
506 greg 2.1
507     /* Build our triangle mesh from recorded RBFs */
508     void
509     build_mesh(void)
510     {
511     double best2 = M_PI*M_PI;
512     RBFNODE *shrt_edj[2];
513     RBFNODE *rbf0, *rbf1;
514 greg 2.30 /* average specular peak */
515     comp_bsdf_spec();
516 greg 2.15 /* add normal if needed */
517     check_normal_incidence();
518 greg 2.1 /* check if isotropic */
519     if (single_plane_incident) {
520     for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
521     if (rbf0->next != NULL)
522     create_migration(rbf0, rbf0->next);
523     await_children(nchild);
524     return;
525     }
526     shrt_edj[0] = shrt_edj[1] = NULL; /* start w/ shortest edge */
527     for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
528     for (rbf1 = rbf0->next; rbf1 != NULL; rbf1 = rbf1->next) {
529     double dist2 = 2. - 2.*DOT(rbf0->invec,rbf1->invec);
530     if (dist2 < best2) {
531     shrt_edj[0] = rbf0;
532     shrt_edj[1] = rbf1;
533     best2 = dist2;
534     }
535     }
536     if (shrt_edj[0] == NULL) {
537     fprintf(stderr, "%s: Cannot find shortest edge\n", progname);
538     exit(1);
539     }
540     /* build mesh from this edge */
541     if (shrt_edj[0]->ord < shrt_edj[1]->ord)
542     mesh_from_edge(create_migration(shrt_edj[0], shrt_edj[1]));
543     else
544     mesh_from_edge(create_migration(shrt_edj[1], shrt_edj[0]));
545     /* complete migrations */
546     await_children(nchild);
547     }