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

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.2 static const char RCSid[] = "$Id: bsdfrep.c,v 2.1 2012/10/19 04:14:29 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     /* 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 greg 2.2 fprintf(ofp, "SYMMETRY=%d\n", !single_plane_incident * inp_coverage);
346     fprintf(ofp, "IO_SIDES= %d %d\n", input_orient, output_orient);
347 greg 2.1 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 greg 2.2 for (mig = mig_list; mig != NULL; mig = mig->next) {
367     int zerocnt = 0;
368 greg 2.1 putint(mig->rbfv[0]->ord, 4, ofp);
369     putint(mig->rbfv[1]->ord, 4, ofp);
370 greg 2.2 /* write out as sparse data */
371 greg 2.1 n = mtx_nrows(mig) * mtx_ncols(mig);
372 greg 2.2 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 greg 2.1 }
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 greg 2.2 /* 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 greg 2.1 /* 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 greg 2.2 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 greg 2.1 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 greg 2.2 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 greg 2.1 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     }