ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2klems.c
Revision: 2.10
Committed: Thu Sep 26 17:05:00 2013 UTC (10 years, 6 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.9: +8 -3 lines
Log Message:
Added -l option to limit maximum number of RBF lobes for interpolation

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.10 static const char RCSid[] = "$Id: bsdf2klems.c,v 2.9 2013/08/11 14:32:34 greg Exp $";
3 greg 2.1 #endif
4     /*
5     * Load measured BSDF interpolant and write out as XML file with Klems matrix.
6     *
7     * G. Ward
8     */
9    
10     #define _USE_MATH_DEFINES
11     #include <stdio.h>
12     #include <stdlib.h>
13     #include <string.h>
14     #include <math.h>
15     #include "random.h"
16     #include "platform.h"
17     #include "calcomp.h"
18     #include "bsdfrep.h"
19     #include "bsdf_m.h"
20 greg 2.4 /* assumed maximum # Klems patches */
21     #define MAXPATCHES 145
22 greg 2.1 /* global argv[0] */
23     char *progname;
24     /* selected basis function name */
25     static const char *kbasis = "LBNL/Klems Full";
26     /* number of BSDF samples per patch */
27     static int npsamps = 256;
28 greg 2.10 /* limit on number of RBF lobes */
29     static int lobe_lim = 15000;
30 greg 2.1
31     /* Return angle basis corresponding to the given name */
32     ANGLE_BASIS *
33     get_basis(const char *bn)
34     {
35     int n = nabases;
36    
37     while (n-- > 0)
38     if (!strcasecmp(bn, abase_list[n].name))
39     return &abase_list[n];
40     return NULL;
41     }
42    
43 greg 2.3 /* Output XML header to stdout */
44     static void
45     xml_header(int ac, char *av[])
46     {
47     puts("<?xml version=\"1.0\" encoding=\"UTF-8\"?>");
48     puts("<WindowElement xmlns=\"http://windows.lbl.gov\" xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\" xsi:schemaLocation=\"http://windows.lbl.gov/BSDF-v1.4.xsd\">");
49     fputs("<!-- File produced by:", stdout);
50     while (ac-- > 0) {
51     fputc(' ', stdout);
52     fputs(*av++, stdout);
53     }
54     puts(" -->");
55     }
56    
57 greg 2.1 /* Output XML prologue to stdout */
58     static void
59 greg 2.3 xml_prologue(const SDData *sd)
60 greg 2.1 {
61 greg 2.3 const char *matn = (sd && sd->matn[0]) ? sd->matn : "Name";
62     const char *makr = (sd && sd->makr[0]) ? sd->makr : "Manufacturer";
63 greg 2.1 ANGLE_BASIS *abp = get_basis(kbasis);
64     int i;
65    
66     if (abp == NULL) {
67     fprintf(stderr, "%s: unknown angle basis '%s'\n", progname, kbasis);
68     exit(1);
69     }
70     puts("<WindowElementType>System</WindowElementType>");
71     puts("<FileType>BSDF</FileType>");
72     puts("<Optical>");
73     puts("<Layer>");
74     puts("\t<Material>");
75 greg 2.3 printf("\t\t<Name>%s</Name>\n", matn);
76     printf("\t\t<Manufacturer>%s</Manufacturer>\n", makr);
77     if (sd && sd->dim[2] > .001) {
78     printf("\t\t<Thickness unit=\"meter\">%.3f</Thickness>\n", sd->dim[2]);
79     printf("\t\t<Width unit=\"meter\">%.3f</Width>\n", sd->dim[0]);
80     printf("\t\t<Height unit=\"meter\">%.3f</Height>\n", sd->dim[1]);
81     }
82 greg 2.1 puts("\t\t<DeviceType>Other</DeviceType>");
83     puts("\t</Material>");
84 greg 2.3 if (sd && sd->mgf != NULL) {
85     puts("\t<Geometry format=\"MGF\">");
86     puts("\t\t<MGFblock unit=\"meter\">");
87     fputs(sd->mgf, stdout);
88     puts("</MGFblock>");
89     puts("\t</Geometry>");
90     }
91 greg 2.1 puts("\t<DataDefinition>");
92     puts("\t\t<IncidentDataStructure>Columns</IncidentDataStructure>");
93     puts("\t\t<AngleBasis>");
94     printf("\t\t\t<AngleBasisName>%s</AngleBasisName>\n", kbasis);
95     for (i = 0; abp->lat[i].nphis; i++) {
96     puts("\t\t\t<AngleBasisBlock>");
97     printf("\t\t\t<Theta>%g</Theta>\n", i ?
98     .5*(abp->lat[i].tmin + abp->lat[i+1].tmin) :
99     .0 );
100 greg 2.3 printf("\t\t\t<nPhis>%d</nPhis>\n", abp->lat[i].nphis);
101 greg 2.1 puts("\t\t\t<ThetaBounds>");
102     printf("\t\t\t\t<LowerTheta>%g</LowerTheta>\n", abp->lat[i].tmin);
103     printf("\t\t\t\t<UpperTheta>%g</UpperTheta>\n", abp->lat[i+1].tmin);
104     puts("\t\t\t</ThetaBounds>");
105     puts("\t\t\t</AngleBasisBlock>");
106     }
107     puts("\t\t</AngleBasis>");
108     puts("\t</DataDefinition>");
109     }
110    
111     /* Output XML data prologue to stdout */
112     static void
113     data_prologue()
114     {
115     static const char *bsdf_type[4] = {
116     "Reflection Front",
117     "Transmission Front",
118     "Transmission Back",
119     "Reflection Back"
120     };
121    
122     puts("\t<WavelengthData>");
123     puts("\t\t<LayerNumber>System</LayerNumber>");
124     puts("\t\t<Wavelength unit=\"Integral\">Visible</Wavelength>");
125     puts("\t\t<SourceSpectrum>CIE Illuminant D65 1nm.ssp</SourceSpectrum>");
126     puts("\t\t<DetectorSpectrum>ASTM E308 1931 Y.dsp</DetectorSpectrum>");
127     puts("\t\t<WavelengthDataBlock>");
128     printf("\t\t\t<WavelengthDataDirection>%s</WavelengthDataDirection>\n",
129     bsdf_type[(input_orient>0)<<1 | (output_orient>0)]);
130     printf("\t\t\t<ColumnAngleBasis>%s</ColumnAngleBasis>\n", kbasis);
131     printf("\t\t\t<RowAngleBasis>%s</RowAngleBasis>\n", kbasis);
132     puts("\t\t\t<ScatteringDataType>BTDF</ScatteringDataType>");
133     puts("\t\t\t<ScatteringData>");
134     }
135    
136     /* Output XML data epilogue to stdout */
137     static void
138     data_epilogue(void)
139     {
140     puts("\t\t\t</ScatteringData>");
141     puts("\t\t</WavelengthDataBlock>");
142     puts("\t</WavelengthData>");
143     }
144    
145     /* Output XML epilogue to stdout */
146     static void
147     xml_epilogue(void)
148     {
149     puts("</Layer>");
150     puts("</Optical>");
151     puts("</WindowElement>");
152     }
153    
154 greg 2.2 /* Load and resample XML BSDF description using Klems basis */
155 greg 2.1 static void
156     eval_bsdf(const char *fname)
157     {
158     ANGLE_BASIS *abp = get_basis(kbasis);
159     SDData bsd;
160     SDError ec;
161     FVECT vin, vout;
162     SDValue sv;
163     double sum;
164     int i, j, n;
165    
166     SDclearBSDF(&bsd, fname); /* load BSDF file */
167     if ((ec = SDloadFile(&bsd, fname)) != SDEnone)
168     goto err;
169 greg 2.3 xml_prologue(&bsd); /* pass geometry */
170 greg 2.1 /* front reflection */
171     if (bsd.rf != NULL || bsd.rLambFront.cieY > .002) {
172     input_orient = 1; output_orient = 1;
173     data_prologue();
174     for (j = 0; j < abp->nangles; j++) {
175     for (i = 0; i < abp->nangles; i++) {
176     sum = 0; /* average over patches */
177     for (n = npsamps; n-- > 0; ) {
178     fo_getvec(vout, j+(n+frandom())/npsamps, abp);
179 greg 2.5 fi_getvec(vin, i+urand(n), abp);
180 greg 2.1 ec = SDevalBSDF(&sv, vout, vin, &bsd);
181     if (ec != SDEnone)
182     goto err;
183     sum += sv.cieY;
184     }
185     printf("\t%.3e\n", sum/npsamps);
186     }
187     putchar('\n'); /* extra space between rows */
188     }
189     data_epilogue();
190     }
191     /* back reflection */
192     if (bsd.rb != NULL || bsd.rLambBack.cieY > .002) {
193     input_orient = -1; output_orient = -1;
194     data_prologue();
195     for (j = 0; j < abp->nangles; j++) {
196     for (i = 0; i < abp->nangles; i++) {
197     sum = 0; /* average over patches */
198     for (n = npsamps; n-- > 0; ) {
199     bo_getvec(vout, j+(n+frandom())/npsamps, abp);
200 greg 2.5 bi_getvec(vin, i+urand(n), abp);
201 greg 2.1 ec = SDevalBSDF(&sv, vout, vin, &bsd);
202     if (ec != SDEnone)
203     goto err;
204     sum += sv.cieY;
205     }
206     printf("\t%.3e\n", sum/npsamps);
207     }
208     putchar('\n'); /* extra space between rows */
209     }
210     data_epilogue();
211     }
212     /* front transmission */
213     if (bsd.tf != NULL || bsd.tLamb.cieY > .002) {
214     input_orient = 1; output_orient = -1;
215     data_prologue();
216     for (j = 0; j < abp->nangles; j++) {
217     for (i = 0; i < abp->nangles; i++) {
218     sum = 0; /* average over patches */
219     for (n = npsamps; n-- > 0; ) {
220     bo_getvec(vout, j+(n+frandom())/npsamps, abp);
221 greg 2.5 fi_getvec(vin, i+urand(n), abp);
222 greg 2.1 ec = SDevalBSDF(&sv, vout, vin, &bsd);
223     if (ec != SDEnone)
224     goto err;
225     sum += sv.cieY;
226     }
227     printf("\t%.3e\n", sum/npsamps);
228     }
229     putchar('\n'); /* extra space between rows */
230     }
231     data_epilogue();
232     }
233     /* back transmission */
234 greg 2.9 if ((bsd.tb != NULL) | (bsd.tf != NULL)) {
235 greg 2.1 input_orient = -1; output_orient = 1;
236     data_prologue();
237     for (j = 0; j < abp->nangles; j++) {
238 greg 2.9 for (i = 0; i < abp->nangles; i++) {
239     sum = 0; /* average over patches */
240     for (n = npsamps; n-- > 0; ) {
241     fo_getvec(vout, j+(n+frandom())/npsamps, abp);
242     bi_getvec(vin, i+urand(n), abp);
243     ec = SDevalBSDF(&sv, vout, vin, &bsd);
244     if (ec != SDEnone)
245 greg 2.1 goto err;
246 greg 2.9 sum += sv.cieY;
247 greg 2.1 }
248 greg 2.9 printf("\t%.3e\n", sum/npsamps);
249     }
250 greg 2.1 putchar('\n'); /* extra space between rows */
251     }
252     data_epilogue();
253     }
254     SDfreeBSDF(&bsd); /* all done */
255     return;
256     err:
257     SDreportError(ec, stderr);
258     exit(1);
259     }
260    
261 greg 2.2 /* Interpolate and output a BSDF function using Klems basis */
262 greg 2.1 static void
263     eval_function(char *funame)
264     {
265     ANGLE_BASIS *abp = get_basis(kbasis);
266 greg 2.8 int assignD = (fundefined(funame) < 6);
267 greg 2.1 double iovec[6];
268     double sum;
269     int i, j, n;
270    
271 greg 2.4 initurand(npsamps);
272 greg 2.1 data_prologue(); /* begin output */
273     for (j = 0; j < abp->nangles; j++) { /* run through directions */
274     for (i = 0; i < abp->nangles; i++) {
275     sum = 0;
276     for (n = npsamps; n--; ) { /* average over patches */
277     if (output_orient > 0)
278     fo_getvec(iovec+3, j+(n+frandom())/npsamps, abp);
279     else
280     bo_getvec(iovec+3, j+(n+frandom())/npsamps, abp);
281    
282     if (input_orient > 0)
283 greg 2.4 fi_getvec(iovec, i+urand(n), abp);
284 greg 2.1 else
285 greg 2.4 bi_getvec(iovec, i+urand(n), abp);
286 greg 2.1
287 greg 2.8 if (assignD) {
288     varset("Dx", '=', -iovec[3]);
289     varset("Dy", '=', -iovec[4]);
290     varset("Dz", '=', -iovec[5]);
291     ++eclock;
292     }
293 greg 2.1 sum += funvalue(funame, 6, iovec);
294     }
295     printf("\t%.3e\n", sum/npsamps);
296     }
297     putchar('\n');
298     }
299     data_epilogue(); /* finish output */
300     }
301    
302     /* Interpolate and output a radial basis function BSDF representation */
303     static void
304     eval_rbf(void)
305     {
306     ANGLE_BASIS *abp = get_basis(kbasis);
307 greg 2.2 float bsdfarr[MAXPATCHES*MAXPATCHES];
308     FVECT vin, vout;
309     RBFNODE *rbf;
310 greg 2.1 double sum;
311     int i, j, n;
312 greg 2.2 /* sanity check */
313     if (abp->nangles > MAXPATCHES) {
314     fprintf(stderr, "%s: too many patches!\n", progname);
315     exit(1);
316     }
317     data_prologue(); /* begin output */
318     for (i = 0; i < abp->nangles; i++) {
319     if (input_orient > 0) /* use incident patch center */
320     fi_getvec(vin, i+.5*(i>0), abp);
321     else
322     bi_getvec(vin, i+.5*(i>0), abp);
323    
324 greg 2.10 rbf = advect_rbf(vin, lobe_lim); /* compute radial basis func */
325 greg 2.2
326     for (j = 0; j < abp->nangles; j++) {
327     sum = 0; /* sample over exiting patch */
328     for (n = npsamps; n--; ) {
329     if (output_orient > 0)
330     fo_getvec(vout, j+(n+frandom())/npsamps, abp);
331     else
332     bo_getvec(vout, j+(n+frandom())/npsamps, abp);
333 greg 2.1
334 greg 2.2 sum += eval_rbfrep(rbf, vout) / vout[2];
335     }
336     bsdfarr[j*abp->nangles + i] = sum*output_orient/npsamps;
337     }
338 greg 2.6 if (rbf != NULL)
339     free(rbf);
340 greg 2.2 }
341     n = 0; /* write out our matrix */
342     for (j = 0; j < abp->nangles; j++) {
343     for (i = 0; i < abp->nangles; i++)
344     printf("\t%.3e\n", bsdfarr[n++]);
345     putchar('\n');
346     }
347     data_epilogue(); /* finish output */
348 greg 2.1 }
349    
350     /* Read in BSDF and interpolate as Klems matrix representation */
351     int
352     main(int argc, char *argv[])
353     {
354     int dofwd = 0, dobwd = 1;
355     char *cp;
356     int i, na;
357    
358     progname = argv[0];
359     esupport |= E_VARIABLE|E_FUNCTION|E_RCONST;
360     esupport &= ~(E_INCHAN|E_OUTCHAN);
361     scompile("PI:3.14159265358979323846", NULL, 0);
362     biggerlib();
363 greg 2.2 for (i = 1; i < argc && (argv[i][0] == '-') | (argv[i][0] == '+'); i++)
364 greg 2.1 switch (argv[i][1]) { /* get options */
365     case 'n':
366     npsamps = atoi(argv[++i]);
367     if (npsamps <= 0)
368     goto userr;
369     break;
370     case 'e':
371     scompile(argv[++i], NULL, 0);
372     single_plane_incident = 0;
373     break;
374     case 'f':
375     if (!argv[i][2]) {
376     fcompile(argv[++i]);
377     single_plane_incident = 0;
378     } else
379     dofwd = (argv[i][0] == '+');
380     break;
381     case 'b':
382     dobwd = (argv[i][0] == '+');
383     break;
384     case 'h':
385     kbasis = "LBNL/Klems Half";
386     break;
387     case 'q':
388     kbasis = "LBNL/Klems Quarter";
389     break;
390 greg 2.10 case 'l':
391     lobe_lim = atoi(argv[++i]);
392     break;
393 greg 2.1 default:
394     goto userr;
395     }
396     if (single_plane_incident >= 0) { /* function-based BSDF? */
397     if (i != argc-1 || fundefined(argv[i]) != 6) {
398     fprintf(stderr,
399     "%s: need single function with 6 arguments: bsdf(ix,iy,iz,ox,oy,oz)\n",
400     progname);
401 greg 2.8 fprintf(stderr, "\tor 3 arguments using Dx,Dy,Dz: bsdf(ix,iy,iz)\n",
402     progname);
403 greg 2.1 goto userr;
404     }
405 greg 2.8 ++eclock;
406 greg 2.3 xml_header(argc, argv); /* start XML output */
407     xml_prologue(NULL);
408 greg 2.1 if (dofwd) {
409     input_orient = -1;
410     output_orient = -1;
411     eval_function(argv[i]); /* outside reflectance */
412     output_orient = 1;
413     eval_function(argv[i]); /* outside -> inside */
414     }
415     if (dobwd) {
416     input_orient = 1;
417     output_orient = 1;
418     eval_function(argv[i]); /* inside reflectance */
419     output_orient = -1;
420     eval_function(argv[i]); /* inside -> outside */
421     }
422     xml_epilogue(); /* finish XML output & exit */
423     return(0);
424     }
425 greg 2.2 /* XML input? */
426     if (i == argc-1 && (cp = argv[i]+strlen(argv[i])-4) > argv[i] &&
427     !strcasecmp(cp, ".xml")) {
428 greg 2.3 xml_header(argc, argv); /* start XML output */
429 greg 2.1 eval_bsdf(argv[i]); /* load & resample BSDF */
430     xml_epilogue(); /* finish XML output & exit */
431     return(0);
432     }
433     if (i < argc) { /* open input files if given */
434     int nbsdf = 0;
435     for ( ; i < argc; i++) { /* interpolate each component */
436     FILE *fpin = fopen(argv[i], "rb");
437     if (fpin == NULL) {
438     fprintf(stderr, "%s: cannot open BSDF interpolant '%s'\n",
439     progname, argv[i]);
440     return(1);
441     }
442     if (!load_bsdf_rep(fpin))
443     return(1);
444     fclose(fpin);
445 greg 2.3 if (!nbsdf++) { /* start XML on first dist. */
446     xml_header(argc, argv);
447     xml_prologue(NULL);
448     }
449 greg 2.1 eval_rbf();
450     }
451     xml_epilogue(); /* finish XML output & exit */
452     return(0);
453     }
454     SET_FILE_BINARY(stdin); /* load from stdin */
455     if (!load_bsdf_rep(stdin))
456     return(1);
457 greg 2.3 xml_header(argc, argv); /* start XML output */
458     xml_prologue(NULL);
459 greg 2.1 eval_rbf(); /* resample dist. */
460     xml_epilogue(); /* finish XML output & exit */
461     return(0);
462     userr:
463     fprintf(stderr,
464 greg 2.10 "Usage: %s [-n spp][-h|-q][-l maxlobes] [bsdf.sir ..] > bsdf.xml\n", progname);
465 greg 2.1 fprintf(stderr,
466 greg 2.3 " or: %s [-n spp][-h|-q] bsdf_in.xml > bsdf_out.xml\n", progname);
467 greg 2.1 fprintf(stderr,
468     " or: %s [-n spp][-h|-q][{+|-}for[ward]][{+|-}b[ackward]][-e expr][-f file] bsdf_func > bsdf.xml\n",
469     progname);
470     return(1);
471     }