ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdfmesh.c
Revision: 2.27
Committed: Wed Mar 26 02:52:31 2014 UTC (10 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.26: +22 -276 lines
Log Message:
Changed to optimized transport matrix computation

File Contents

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