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

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.5 static const char RCSid[] = "$Id: bsdfrep.c,v 2.4 2012/10/23 05:10:42 greg Exp $";
3 greg 2.1 #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 greg 2.2 #include <string.h>
13 greg 2.1 #include <math.h>
14     #include "rtio.h"
15     #include "resolu.h"
16     #include "bsdfrep.h"
17 greg 2.5 /* active grid resolution */
18     int grid_res = GRIDRES;
19    
20 greg 2.4 /* coverage/symmetry using INP_QUAD? flags */
21 greg 2.1 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 greg 2.5 rbf->rbfa[n].gx = grid_res-1 - rbf->rbfa[n].gx;
175 greg 2.1 if (sym & MIRROR_Y)
176     for (n = rbf->nrbf; n-- > 0; )
177 greg 2.5 rbf->rbfa[n].gy = grid_res-1 - rbf->rbfa[n].gy;
178 greg 2.1 }
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 greg 2.5 SDsquare2disk(uv, (1./grid_res)*(xpos+.5), (1./grid_res)*(ypos+.5));
197 greg 2.1 /* 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 greg 2.5 pos[0] = (int)(sq[0]*grid_res);
215     pos[1] = (int)(sq[1]*grid_res);
216 greg 2.1 }
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 greg 2.3 if (rbf->ord == ord)
288 greg 2.1 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 greg 2.4 const MIGRATION *ej1, *ej2;
319 greg 2.1 RBFNODE *tv;
320    
321     rbfv[0] = rbfv[1] = NULL;
322     if (mig == NULL)
323     return(0);
324 greg 2.4 for (ej1 = mig->rbfv[0]->ejl; ej1 != NULL;
325     ej1 = nextedge(mig->rbfv[0],ej1)) {
326     if (ej1 == mig)
327 greg 2.1 continue;
328 greg 2.4 tv = opp_rbf(mig->rbfv[0],ej1);
329 greg 2.1 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 greg 2.4 /* 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 greg 2.5 grid_res = GRIDRES;
358 greg 2.4 }
359    
360 greg 2.1 /* 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 greg 2.2 fprintf(ofp, "SYMMETRY=%d\n", !single_plane_incident * inp_coverage);
369     fprintf(ofp, "IO_SIDES= %d %d\n", input_orient, output_orient);
370 greg 2.5 fprintf(ofp, "GRIDRES=%d\n", grid_res);
371 greg 2.1 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 greg 2.2 for (mig = mig_list; mig != NULL; mig = mig->next) {
391     int zerocnt = 0;
392 greg 2.1 putint(mig->rbfv[0]->ord, 4, ofp);
393     putint(mig->rbfv[1]->ord, 4, ofp);
394 greg 2.2 /* write out as sparse data */
395 greg 2.1 n = mtx_nrows(mig) * mtx_ncols(mig);
396 greg 2.2 for (i = 0; i < n; i++) {
397 greg 2.3 if (zerocnt == 0xff) {
398     putint(0xff, 1, ofp); zerocnt = 0;
399 greg 2.2 }
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 greg 2.1 }
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 greg 2.2 /* 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 greg 2.5 if (!strncmp(s, "GRIDRES=", 8)) {
433     sscanf(s+8, "%d", &grid_res);
434     return(0);
435     }
436 greg 2.2 if (formatval(fmt, s) && strcmp(fmt, BSDFREP_FMT))
437     return(-1);
438     return(0);
439     }
440    
441 greg 2.1 /* 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 greg 2.4
449     clear_bsdf_rep();
450 greg 2.5 if (ifp == NULL)
451     return(0);
452 greg 2.2 if (getheader(ifp, headline, NULL) < 0 || single_plane_incident < 0 |
453     !input_orient | !output_orient) {
454 greg 2.1 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 greg 2.3 rbfh.vtotal = getflt(ifp);
467 greg 2.1 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 greg 2.2 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 greg 2.3 if (zc == 0xff)
514     continue;
515 greg 2.2 newmig->mtx[i++] = getflt(ifp);
516     }
517 greg 2.1 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     }