ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdfrep.c
Revision: 2.5
Committed: Wed Nov 7 03:04:23 2012 UTC (11 years, 5 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.4: +17 -6 lines
Log Message:
Added GRIDRES to saved BSDF representation

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: bsdfrep.c,v 2.4 2012/10/23 05:10:42 greg Exp $";
3 #endif
4 /*
5 * Support BSDF representation as radial basis functions.
6 *
7 * G. Ward
8 */
9
10 #define _USE_MATH_DEFINES
11 #include <stdlib.h>
12 #include <string.h>
13 #include <math.h>
14 #include "rtio.h"
15 #include "resolu.h"
16 #include "bsdfrep.h"
17 /* active grid resolution */
18 int grid_res = GRIDRES;
19
20 /* coverage/symmetry using INP_QUAD? flags */
21 int inp_coverage = 0;
22 /* all incident angles in-plane so far? */
23 int single_plane_incident = -1;
24
25 /* input/output orientations */
26 int input_orient = 0;
27 int output_orient = 0;
28
29 /* processed incident DSF measurements */
30 RBFNODE *dsf_list = NULL;
31
32 /* RBF-linking matrices (edges) */
33 MIGRATION *mig_list = NULL;
34
35 /* current input direction */
36 double theta_in_deg, phi_in_deg;
37
38 /* Register new input direction */
39 int
40 new_input_direction(double new_theta, double new_phi)
41 {
42 if (!input_orient) /* check input orientation */
43 input_orient = 1 - 2*(new_theta > 90.);
44 else if (input_orient > 0 ^ new_theta < 90.) {
45 fprintf(stderr,
46 "%s: Cannot handle input angles on both sides of surface\n",
47 progname);
48 return(0);
49 }
50 /* normalize angle ranges */
51 while (new_theta < -180.)
52 new_theta += 360.;
53 while (new_theta > 180.)
54 new_theta -= 360.;
55 if (new_theta < 0) {
56 new_theta = -new_theta;
57 new_phi += 180.;
58 }
59 while (new_phi < 0)
60 new_phi += 360.;
61 while (new_phi >= 360.)
62 new_phi -= 360.;
63 if (single_plane_incident > 0) /* check input coverage */
64 single_plane_incident = (round(new_phi) == round(phi_in_deg));
65 else if (single_plane_incident < 0)
66 single_plane_incident = 1;
67 theta_in_deg = new_theta; /* assume it's OK */
68 phi_in_deg = new_phi;
69 if ((1. < new_phi) & (new_phi < 89.))
70 inp_coverage |= INP_QUAD1;
71 else if ((91. < new_phi) & (new_phi < 179.))
72 inp_coverage |= INP_QUAD2;
73 else if ((181. < new_phi) & (new_phi < 269.))
74 inp_coverage |= INP_QUAD3;
75 else if ((271. < new_phi) & (new_phi < 359.))
76 inp_coverage |= INP_QUAD4;
77 return(1);
78 }
79
80 /* Apply symmetry to the given vector based on distribution */
81 int
82 use_symmetry(FVECT vec)
83 {
84 double phi = get_phi360(vec);
85
86 switch (inp_coverage) {
87 case INP_QUAD1|INP_QUAD2|INP_QUAD3|INP_QUAD4:
88 break;
89 case INP_QUAD1|INP_QUAD2:
90 if ((-FTINY > phi) | (phi > 180.+FTINY))
91 goto mir_y;
92 break;
93 case INP_QUAD2|INP_QUAD3:
94 if ((90.-FTINY > phi) | (phi > 270.+FTINY))
95 goto mir_x;
96 break;
97 case INP_QUAD3|INP_QUAD4:
98 if ((180.-FTINY > phi) | (phi > 360.+FTINY))
99 goto mir_y;
100 break;
101 case INP_QUAD4|INP_QUAD1:
102 if ((270.-FTINY > phi) & (phi > 90.+FTINY))
103 goto mir_x;
104 break;
105 case INP_QUAD1:
106 if ((-FTINY > phi) | (phi > 90.+FTINY))
107 switch ((int)(phi*(1./90.))) {
108 case 1: goto mir_x;
109 case 2: goto mir_xy;
110 case 3: goto mir_y;
111 }
112 break;
113 case INP_QUAD2:
114 if ((90.-FTINY > phi) | (phi > 180.+FTINY))
115 switch ((int)(phi*(1./90.))) {
116 case 0: goto mir_x;
117 case 2: goto mir_y;
118 case 3: goto mir_xy;
119 }
120 break;
121 case INP_QUAD3:
122 if ((180.-FTINY > phi) | (phi > 270.+FTINY))
123 switch ((int)(phi*(1./90.))) {
124 case 0: goto mir_xy;
125 case 1: goto mir_y;
126 case 3: goto mir_x;
127 }
128 break;
129 case INP_QUAD4:
130 if ((270.-FTINY > phi) | (phi > 360.+FTINY))
131 switch ((int)(phi*(1./90.))) {
132 case 0: goto mir_y;
133 case 1: goto mir_xy;
134 case 2: goto mir_x;
135 }
136 break;
137 default:
138 fprintf(stderr, "%s: Illegal input coverage (%d)\n",
139 progname, inp_coverage);
140 exit(1);
141 }
142 return(0); /* in range */
143 mir_x:
144 vec[0] = -vec[0];
145 return(MIRROR_X);
146 mir_y:
147 vec[1] = -vec[1];
148 return(MIRROR_Y);
149 mir_xy:
150 vec[0] = -vec[0];
151 vec[1] = -vec[1];
152 return(MIRROR_X|MIRROR_Y);
153 }
154
155 /* Reverse symmetry based on what was done before */
156 void
157 rev_symmetry(FVECT vec, int sym)
158 {
159 if (sym & MIRROR_X)
160 vec[0] = -vec[0];
161 if (sym & MIRROR_Y)
162 vec[1] = -vec[1];
163 }
164
165 /* Reverse symmetry for an RBF distribution */
166 void
167 rev_rbf_symmetry(RBFNODE *rbf, int sym)
168 {
169 int n;
170
171 rev_symmetry(rbf->invec, sym);
172 if (sym & MIRROR_X)
173 for (n = rbf->nrbf; n-- > 0; )
174 rbf->rbfa[n].gx = grid_res-1 - rbf->rbfa[n].gx;
175 if (sym & MIRROR_Y)
176 for (n = rbf->nrbf; n-- > 0; )
177 rbf->rbfa[n].gy = grid_res-1 - rbf->rbfa[n].gy;
178 }
179
180 /* Compute volume associated with Gaussian lobe */
181 double
182 rbf_volume(const RBFVAL *rbfp)
183 {
184 double rad = R2ANG(rbfp->crad);
185
186 return((2.*M_PI) * rbfp->peak * rad*rad);
187 }
188
189 /* Compute outgoing vector from grid position */
190 void
191 ovec_from_pos(FVECT vec, int xpos, int ypos)
192 {
193 double uv[2];
194 double r2;
195
196 SDsquare2disk(uv, (1./grid_res)*(xpos+.5), (1./grid_res)*(ypos+.5));
197 /* uniform hemispherical projection */
198 r2 = uv[0]*uv[0] + uv[1]*uv[1];
199 vec[0] = vec[1] = sqrt(2. - r2);
200 vec[0] *= uv[0];
201 vec[1] *= uv[1];
202 vec[2] = output_orient*(1. - r2);
203 }
204
205 /* Compute grid position from normalized input/output vector */
206 void
207 pos_from_vec(int pos[2], const FVECT vec)
208 {
209 double sq[2]; /* uniform hemispherical projection */
210 double norm = 1./sqrt(1. + fabs(vec[2]));
211
212 SDdisk2square(sq, vec[0]*norm, vec[1]*norm);
213
214 pos[0] = (int)(sq[0]*grid_res);
215 pos[1] = (int)(sq[1]*grid_res);
216 }
217
218 /* Evaluate RBF for DSF at the given normalized outgoing direction */
219 double
220 eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
221 {
222 double res = .0;
223 const RBFVAL *rbfp;
224 FVECT odir;
225 double sig2;
226 int n;
227
228 if (rp == NULL)
229 return(.0);
230 rbfp = rp->rbfa;
231 for (n = rp->nrbf; n--; rbfp++) {
232 ovec_from_pos(odir, rbfp->gx, rbfp->gy);
233 sig2 = R2ANG(rbfp->crad);
234 sig2 = (DOT(odir,outvec) - 1.) / (sig2*sig2);
235 if (sig2 > -19.)
236 res += rbfp->peak * exp(sig2);
237 }
238 return(res);
239 }
240
241 /* Insert a new directional scattering function in our global list */
242 int
243 insert_dsf(RBFNODE *newrbf)
244 {
245 RBFNODE *rbf, *rbf_last;
246 int pos;
247 /* check for redundant meas. */
248 for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
249 if (DOT(rbf->invec, newrbf->invec) >= 1.-FTINY) {
250 fprintf(stderr,
251 "%s: Duplicate incident measurement (ignored)\n",
252 progname);
253 free(newrbf);
254 return(-1);
255 }
256 /* keep in ascending theta order */
257 for (rbf_last = NULL, rbf = dsf_list; rbf != NULL;
258 rbf_last = rbf, rbf = rbf->next)
259 if (single_plane_incident && input_orient*rbf->invec[2] <
260 input_orient*newrbf->invec[2])
261 break;
262 if (rbf_last == NULL) { /* insert new node in list */
263 newrbf->ord = 0;
264 newrbf->next = dsf_list;
265 dsf_list = newrbf;
266 } else {
267 newrbf->ord = rbf_last->ord + 1;
268 newrbf->next = rbf;
269 rbf_last->next = newrbf;
270 }
271 rbf_last = newrbf;
272 while (rbf != NULL) { /* update ordinal positions */
273 rbf->ord = rbf_last->ord + 1;
274 rbf_last = rbf;
275 rbf = rbf->next;
276 }
277 return(newrbf->ord);
278 }
279
280 /* Get the DSF indicated by its ordinal position */
281 RBFNODE *
282 get_dsf(int ord)
283 {
284 RBFNODE *rbf;
285
286 for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
287 if (rbf->ord == ord)
288 return(rbf);
289 return(NULL);
290 }
291
292 /* Get triangle surface orientation (unnormalized) */
293 void
294 tri_orient(FVECT vres, const FVECT v1, const FVECT v2, const FVECT v3)
295 {
296 FVECT v2minus1, v3minus2;
297
298 VSUB(v2minus1, v2, v1);
299 VSUB(v3minus2, v3, v2);
300 VCROSS(vres, v2minus1, v3minus2);
301 }
302
303 /* Determine if vertex order is reversed (inward normal) */
304 int
305 is_rev_tri(const FVECT v1, const FVECT v2, const FVECT v3)
306 {
307 FVECT tor;
308
309 tri_orient(tor, v1, v2, v3);
310
311 return(DOT(tor, v2) < 0.);
312 }
313
314 /* Find vertices completing triangles on either side of the given edge */
315 int
316 get_triangles(RBFNODE *rbfv[2], const MIGRATION *mig)
317 {
318 const MIGRATION *ej1, *ej2;
319 RBFNODE *tv;
320
321 rbfv[0] = rbfv[1] = NULL;
322 if (mig == NULL)
323 return(0);
324 for (ej1 = mig->rbfv[0]->ejl; ej1 != NULL;
325 ej1 = nextedge(mig->rbfv[0],ej1)) {
326 if (ej1 == mig)
327 continue;
328 tv = opp_rbf(mig->rbfv[0],ej1);
329 for (ej2 = tv->ejl; ej2 != NULL; ej2 = nextedge(tv,ej2))
330 if (opp_rbf(tv,ej2) == mig->rbfv[1]) {
331 rbfv[is_rev_tri(mig->rbfv[0]->invec,
332 mig->rbfv[1]->invec,
333 tv->invec)] = tv;
334 break;
335 }
336 }
337 return((rbfv[0] != NULL) + (rbfv[1] != NULL));
338 }
339
340 /* Clear our BSDF representation and free memory */
341 void
342 clear_bsdf_rep(void)
343 {
344 while (mig_list != NULL) {
345 MIGRATION *mig = mig_list;
346 mig_list = mig->next;
347 free(mig);
348 }
349 while (dsf_list != NULL) {
350 RBFNODE *rbf = dsf_list;
351 dsf_list = rbf->next;
352 free(rbf);
353 }
354 inp_coverage = 0;
355 single_plane_incident = -1;
356 input_orient = output_orient = 0;
357 grid_res = GRIDRES;
358 }
359
360 /* Write our BSDF mesh interpolant out to the given binary stream */
361 void
362 save_bsdf_rep(FILE *ofp)
363 {
364 RBFNODE *rbf;
365 MIGRATION *mig;
366 int i, n;
367 /* finish header */
368 fprintf(ofp, "SYMMETRY=%d\n", !single_plane_incident * inp_coverage);
369 fprintf(ofp, "IO_SIDES= %d %d\n", input_orient, output_orient);
370 fprintf(ofp, "GRIDRES=%d\n", grid_res);
371 fputformat(BSDFREP_FMT, ofp);
372 fputc('\n', ofp);
373 /* write each DSF */
374 for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
375 putint(rbf->ord, 4, ofp);
376 putflt(rbf->invec[0], ofp);
377 putflt(rbf->invec[1], ofp);
378 putflt(rbf->invec[2], ofp);
379 putflt(rbf->vtotal, ofp);
380 putint(rbf->nrbf, 4, ofp);
381 for (i = 0; i < rbf->nrbf; i++) {
382 putflt(rbf->rbfa[i].peak, ofp);
383 putint(rbf->rbfa[i].crad, 2, ofp);
384 putint(rbf->rbfa[i].gx, 1, ofp);
385 putint(rbf->rbfa[i].gy, 1, ofp);
386 }
387 }
388 putint(-1, 4, ofp); /* terminator */
389 /* write each migration matrix */
390 for (mig = mig_list; mig != NULL; mig = mig->next) {
391 int zerocnt = 0;
392 putint(mig->rbfv[0]->ord, 4, ofp);
393 putint(mig->rbfv[1]->ord, 4, ofp);
394 /* write out as sparse data */
395 n = mtx_nrows(mig) * mtx_ncols(mig);
396 for (i = 0; i < n; i++) {
397 if (zerocnt == 0xff) {
398 putint(0xff, 1, ofp); zerocnt = 0;
399 }
400 if (mig->mtx[i] != 0) {
401 putint(zerocnt, 1, ofp); zerocnt = 0;
402 putflt(mig->mtx[i], ofp);
403 } else
404 ++zerocnt;
405 }
406 putint(zerocnt, 1, ofp);
407 }
408 putint(-1, 4, ofp); /* terminator */
409 putint(-1, 4, ofp);
410 if (fflush(ofp) == EOF) {
411 fprintf(stderr, "%s: error writing BSDF interpolant\n",
412 progname);
413 exit(1);
414 }
415 }
416
417 /* Check header line for critical information */
418 static int
419 headline(char *s, void *p)
420 {
421 char fmt[32];
422
423 if (!strncmp(s, "SYMMETRY=", 9)) {
424 inp_coverage = atoi(s+9);
425 single_plane_incident = !inp_coverage;
426 return(0);
427 }
428 if (!strncmp(s, "IO_SIDES=", 9)) {
429 sscanf(s+9, "%d %d", &input_orient, &output_orient);
430 return(0);
431 }
432 if (!strncmp(s, "GRIDRES=", 8)) {
433 sscanf(s+8, "%d", &grid_res);
434 return(0);
435 }
436 if (formatval(fmt, s) && strcmp(fmt, BSDFREP_FMT))
437 return(-1);
438 return(0);
439 }
440
441 /* Read a BSDF mesh interpolant from the given binary stream */
442 int
443 load_bsdf_rep(FILE *ifp)
444 {
445 RBFNODE rbfh;
446 int from_ord, to_ord;
447 int i;
448
449 clear_bsdf_rep();
450 if (ifp == NULL)
451 return(0);
452 if (getheader(ifp, headline, NULL) < 0 || single_plane_incident < 0 |
453 !input_orient | !output_orient) {
454 fprintf(stderr, "%s: missing/bad format for BSDF interpolant\n",
455 progname);
456 return(0);
457 }
458 rbfh.next = NULL; /* read each DSF */
459 rbfh.ejl = NULL;
460 while ((rbfh.ord = getint(4, ifp)) >= 0) {
461 RBFNODE *newrbf;
462
463 rbfh.invec[0] = getflt(ifp);
464 rbfh.invec[1] = getflt(ifp);
465 rbfh.invec[2] = getflt(ifp);
466 rbfh.vtotal = getflt(ifp);
467 rbfh.nrbf = getint(4, ifp);
468 newrbf = (RBFNODE *)malloc(sizeof(RBFNODE) +
469 sizeof(RBFVAL)*(rbfh.nrbf-1));
470 if (newrbf == NULL)
471 goto memerr;
472 memcpy(newrbf, &rbfh, sizeof(RBFNODE));
473 for (i = 0; i < rbfh.nrbf; i++) {
474 newrbf->rbfa[i].peak = getflt(ifp);
475 newrbf->rbfa[i].crad = getint(2, ifp) & 0xffff;
476 newrbf->rbfa[i].gx = getint(1, ifp) & 0xff;
477 newrbf->rbfa[i].gy = getint(1, ifp) & 0xff;
478 }
479 if (feof(ifp))
480 goto badEOF;
481 /* insert in global list */
482 if (insert_dsf(newrbf) != rbfh.ord) {
483 fprintf(stderr, "%s: error adding DSF\n", progname);
484 return(0);
485 }
486 }
487 /* read each migration matrix */
488 while ((from_ord = getint(4, ifp)) >= 0 &&
489 (to_ord = getint(4, ifp)) >= 0) {
490 RBFNODE *from_rbf = get_dsf(from_ord);
491 RBFNODE *to_rbf = get_dsf(to_ord);
492 MIGRATION *newmig;
493 int n;
494
495 if ((from_rbf == NULL) | (to_rbf == NULL)) {
496 fprintf(stderr,
497 "%s: bad DSF reference in migration edge\n",
498 progname);
499 return(0);
500 }
501 n = from_rbf->nrbf * to_rbf->nrbf;
502 newmig = (MIGRATION *)malloc(sizeof(MIGRATION) +
503 sizeof(float)*(n-1));
504 if (newmig == NULL)
505 goto memerr;
506 newmig->rbfv[0] = from_rbf;
507 newmig->rbfv[1] = to_rbf;
508 memset(newmig->mtx, 0, sizeof(float)*n);
509 for (i = 0; ; ) { /* read sparse data */
510 int zc = getint(1, ifp) & 0xff;
511 if ((i += zc) >= n)
512 break;
513 if (zc == 0xff)
514 continue;
515 newmig->mtx[i++] = getflt(ifp);
516 }
517 if (feof(ifp))
518 goto badEOF;
519 /* insert in edge lists */
520 newmig->enxt[0] = from_rbf->ejl;
521 from_rbf->ejl = newmig;
522 newmig->enxt[1] = to_rbf->ejl;
523 to_rbf->ejl = newmig;
524 /* push onto global list */
525 newmig->next = mig_list;
526 mig_list = newmig;
527 }
528 return(1); /* success! */
529 memerr:
530 fprintf(stderr, "%s: Out of memory in load_bsdf_rep()\n", progname);
531 exit(1);
532 badEOF:
533 fprintf(stderr, "%s: Unexpected EOF in load_bsdf_rep()\n", progname);
534 return(0);
535 }