ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdfrep.c
Revision: 2.2
Committed: Sat Oct 20 07:02:00 2012 UTC (11 years, 6 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.1: +53 -8 lines
Log Message:
Bug fixes and minor improvements

File Contents

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