ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdfrep.c
Revision: 2.4
Committed: Tue Oct 23 05:10:42 2012 UTC (11 years, 6 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.3: +28 -17 lines
Log Message:
Hopeful bug fix in BSDF interpolation code

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.4 static const char RCSid[] = "$Id: bsdfrep.c,v 2.3 2012/10/20 17:01:26 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.4 /* coverage/symmetry using INP_QUAD? flags */
18 greg 2.1 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 greg 2.3 if (rbf->ord == ord)
285 greg 2.1 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 greg 2.4 const MIGRATION *ej1, *ej2;
316 greg 2.1 RBFNODE *tv;
317    
318     rbfv[0] = rbfv[1] = NULL;
319     if (mig == NULL)
320     return(0);
321 greg 2.4 for (ej1 = mig->rbfv[0]->ejl; ej1 != NULL;
322     ej1 = nextedge(mig->rbfv[0],ej1)) {
323     if (ej1 == mig)
324 greg 2.1 continue;
325 greg 2.4 tv = opp_rbf(mig->rbfv[0],ej1);
326 greg 2.1 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 greg 2.4 /* Clear our BSDF representation and free memory */
338     void
339     clear_bsdf_rep(void)
340     {
341     while (mig_list != NULL) {
342     MIGRATION *mig = mig_list;
343     mig_list = mig->next;
344     free(mig);
345     }
346     while (dsf_list != NULL) {
347     RBFNODE *rbf = dsf_list;
348     dsf_list = rbf->next;
349     free(rbf);
350     }
351     inp_coverage = 0;
352     single_plane_incident = -1;
353     input_orient = output_orient = 0;
354     }
355    
356 greg 2.1 /* Write our BSDF mesh interpolant out to the given binary stream */
357     void
358     save_bsdf_rep(FILE *ofp)
359     {
360     RBFNODE *rbf;
361     MIGRATION *mig;
362     int i, n;
363     /* finish header */
364 greg 2.2 fprintf(ofp, "SYMMETRY=%d\n", !single_plane_incident * inp_coverage);
365     fprintf(ofp, "IO_SIDES= %d %d\n", input_orient, output_orient);
366 greg 2.1 fputformat(BSDFREP_FMT, ofp);
367     fputc('\n', ofp);
368     /* write each DSF */
369     for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
370     putint(rbf->ord, 4, ofp);
371     putflt(rbf->invec[0], ofp);
372     putflt(rbf->invec[1], ofp);
373     putflt(rbf->invec[2], ofp);
374     putflt(rbf->vtotal, ofp);
375     putint(rbf->nrbf, 4, ofp);
376     for (i = 0; i < rbf->nrbf; i++) {
377     putflt(rbf->rbfa[i].peak, ofp);
378     putint(rbf->rbfa[i].crad, 2, ofp);
379     putint(rbf->rbfa[i].gx, 1, ofp);
380     putint(rbf->rbfa[i].gy, 1, ofp);
381     }
382     }
383     putint(-1, 4, ofp); /* terminator */
384     /* write each migration matrix */
385 greg 2.2 for (mig = mig_list; mig != NULL; mig = mig->next) {
386     int zerocnt = 0;
387 greg 2.1 putint(mig->rbfv[0]->ord, 4, ofp);
388     putint(mig->rbfv[1]->ord, 4, ofp);
389 greg 2.2 /* write out as sparse data */
390 greg 2.1 n = mtx_nrows(mig) * mtx_ncols(mig);
391 greg 2.2 for (i = 0; i < n; i++) {
392 greg 2.3 if (zerocnt == 0xff) {
393     putint(0xff, 1, ofp); zerocnt = 0;
394 greg 2.2 }
395     if (mig->mtx[i] != 0) {
396     putint(zerocnt, 1, ofp); zerocnt = 0;
397     putflt(mig->mtx[i], ofp);
398     } else
399     ++zerocnt;
400     }
401     putint(zerocnt, 1, ofp);
402 greg 2.1 }
403     putint(-1, 4, ofp); /* terminator */
404     putint(-1, 4, ofp);
405     if (fflush(ofp) == EOF) {
406     fprintf(stderr, "%s: error writing BSDF interpolant\n",
407     progname);
408     exit(1);
409     }
410     }
411    
412 greg 2.2 /* Check header line for critical information */
413     static int
414     headline(char *s, void *p)
415     {
416     char fmt[32];
417    
418     if (!strncmp(s, "SYMMETRY=", 9)) {
419     inp_coverage = atoi(s+9);
420     single_plane_incident = !inp_coverage;
421     return(0);
422     }
423     if (!strncmp(s, "IO_SIDES=", 9)) {
424     sscanf(s+9, "%d %d", &input_orient, &output_orient);
425     return(0);
426     }
427     if (formatval(fmt, s) && strcmp(fmt, BSDFREP_FMT))
428     return(-1);
429     return(0);
430     }
431    
432 greg 2.1 /* Read a BSDF mesh interpolant from the given binary stream */
433     int
434     load_bsdf_rep(FILE *ifp)
435     {
436     RBFNODE rbfh;
437     int from_ord, to_ord;
438     int i;
439 greg 2.4
440     clear_bsdf_rep();
441 greg 2.2 if (getheader(ifp, headline, NULL) < 0 || single_plane_incident < 0 |
442     !input_orient | !output_orient) {
443 greg 2.1 fprintf(stderr, "%s: missing/bad format for BSDF interpolant\n",
444     progname);
445     return(0);
446     }
447     rbfh.next = NULL; /* read each DSF */
448     rbfh.ejl = NULL;
449     while ((rbfh.ord = getint(4, ifp)) >= 0) {
450     RBFNODE *newrbf;
451    
452     rbfh.invec[0] = getflt(ifp);
453     rbfh.invec[1] = getflt(ifp);
454     rbfh.invec[2] = getflt(ifp);
455 greg 2.3 rbfh.vtotal = getflt(ifp);
456 greg 2.1 rbfh.nrbf = getint(4, ifp);
457     newrbf = (RBFNODE *)malloc(sizeof(RBFNODE) +
458     sizeof(RBFVAL)*(rbfh.nrbf-1));
459     if (newrbf == NULL)
460     goto memerr;
461     memcpy(newrbf, &rbfh, sizeof(RBFNODE));
462     for (i = 0; i < rbfh.nrbf; i++) {
463     newrbf->rbfa[i].peak = getflt(ifp);
464     newrbf->rbfa[i].crad = getint(2, ifp) & 0xffff;
465     newrbf->rbfa[i].gx = getint(1, ifp) & 0xff;
466     newrbf->rbfa[i].gy = getint(1, ifp) & 0xff;
467     }
468     if (feof(ifp))
469     goto badEOF;
470     /* insert in global list */
471     if (insert_dsf(newrbf) != rbfh.ord) {
472     fprintf(stderr, "%s: error adding DSF\n", progname);
473     return(0);
474     }
475     }
476     /* read each migration matrix */
477     while ((from_ord = getint(4, ifp)) >= 0 &&
478     (to_ord = getint(4, ifp)) >= 0) {
479     RBFNODE *from_rbf = get_dsf(from_ord);
480     RBFNODE *to_rbf = get_dsf(to_ord);
481     MIGRATION *newmig;
482     int n;
483    
484     if ((from_rbf == NULL) | (to_rbf == NULL)) {
485     fprintf(stderr,
486     "%s: bad DSF reference in migration edge\n",
487     progname);
488     return(0);
489     }
490     n = from_rbf->nrbf * to_rbf->nrbf;
491     newmig = (MIGRATION *)malloc(sizeof(MIGRATION) +
492     sizeof(float)*(n-1));
493     if (newmig == NULL)
494     goto memerr;
495     newmig->rbfv[0] = from_rbf;
496     newmig->rbfv[1] = to_rbf;
497 greg 2.2 memset(newmig->mtx, 0, sizeof(float)*n);
498     for (i = 0; ; ) { /* read sparse data */
499     int zc = getint(1, ifp) & 0xff;
500     if ((i += zc) >= n)
501     break;
502 greg 2.3 if (zc == 0xff)
503     continue;
504 greg 2.2 newmig->mtx[i++] = getflt(ifp);
505     }
506 greg 2.1 if (feof(ifp))
507     goto badEOF;
508     /* insert in edge lists */
509     newmig->enxt[0] = from_rbf->ejl;
510     from_rbf->ejl = newmig;
511     newmig->enxt[1] = to_rbf->ejl;
512     to_rbf->ejl = newmig;
513     /* push onto global list */
514     newmig->next = mig_list;
515     mig_list = newmig;
516     }
517     return(1); /* success! */
518     memerr:
519     fprintf(stderr, "%s: Out of memory in load_bsdf_rep()\n", progname);
520     exit(1);
521     badEOF:
522     fprintf(stderr, "%s: Unexpected EOF in load_bsdf_rep()\n", progname);
523     return(0);
524     }