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

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: bsdfmesh.c,v 2.32 2014/08/21 19:31:23 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
22 #ifndef NEIGH_FACT2
23 #define NEIGH_FACT2 0.1 /* empirical neighborhood distance weight */
24 #endif
25 /* number of processes to run */
26 int nprocs = 1;
27 /* number of children (-1 in child) */
28 static int nchild = 0;
29
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 sum += eval_rbfrep(rbf, tvec) * COSF(tvec[2]);
55 ++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 const double DSFtarget = 0.60653066 * eval_rbfrep(rbf,outvec) *
68 COSF(outvec[2]);
69 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 /* Newton's method (sort of) */
76 do {
77 double test_rad = interp_rad;
78 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 inside_rad = test_rad;
86 DSFinside = DSFtest;
87 } else {
88 outside_rad = test_rad;
89 DSFoutside = DSFtest;
90 }
91 if (DSFoutside >= DSFinside)
92 return(test_rad);
93 } while (outside_rad-inside_rad > rad_epsilon);
94 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 bsdf_spec_rad = 0;
111 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 bsdf_spec_rad = rad_sum/(double)n;
123 }
124
125 /* 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 await_children(nchild);
213 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 /* 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 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 }
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 /* 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 }
286
287
288 /* 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 int i, j;
294 /* 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 #ifdef DEBUG
301 fprintf(stderr, "Building path from (theta,phi) (%.1f,%.1f) ",
302 get_theta180(from_rbf->invec),
303 get_phi360(from_rbf->invec));
304 fprintf(stderr, "to (%.1f,%.1f) with %d x %d matrix\n",
305 get_theta180(to_rbf->invec),
306 get_phi360(to_rbf->invec),
307 from_rbf->nrbf, to_rbf->nrbf);
308 #endif
309 newmig = new_migration(from_rbf, to_rbf);
310 if (run_subprocess())
311 return(newmig); /* child continues */
312
313 /* compute transport plan */
314 compute_nDSFs(from_rbf, to_rbf);
315 plan_transport(newmig);
316
317 for (i = from_rbf->nrbf; i--; ) { /* normalize final matrix */
318 double nf = rbf_volume(&from_rbf->rbfa[i]);
319 if (nf <= FTINY) continue;
320 nf = from_rbf->vtotal / nf;
321 for (j = to_rbf->nrbf; j--; )
322 mtx_coef(newmig,i,j) *= nf; /* row now sums to 1.0 */
323 }
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 /* Find convex hull vertex to complete triangle (oriented call) */
355 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 VSUB(vp, rbf->invec, vmid);
376 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 return;
415 }
416 }
417 if (tvert[1] == NULL) { /* grow mesh on left */
418 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
434 /* Add normal direction if missing */
435 static void
436 check_normal_incidence(void)
437 {
438 static FVECT norm_vec = {.0, .0, 1.};
439 const int saved_nprocs = nprocs;
440 RBFNODE *near_rbf, *mir_rbf, *rbf;
441 double bestd;
442 int n;
443
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 mir_rbf->ejl = NULL;
489 rev_rbf_symmetry(mir_rbf, MIRROR_X|MIRROR_Y);
490 nprocs = 1; /* compute migration matrix */
491 if (create_migration(mir_rbf, near_rbf) == NULL)
492 exit(1); /* XXX should never happen! */
493 norm_vec[2] = input_orient; /* interpolate normal dist. */
494 rbf = e_advect_rbf(mig_list, norm_vec, 0);
495 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
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 /* average specular peak */
515 comp_bsdf_spec();
516 /* add normal if needed */
517 check_normal_incidence();
518 /* 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 }