ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2klems.c
Revision: 2.3
Committed: Tue Apr 23 23:19:09 2013 UTC (10 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.2: +45 -22 lines
Log Message:
New bsdf2klems tool seems to pass basic tests

File Contents

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