ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2klems.c
Revision: 2.12
Committed: Thu Nov 21 20:09:06 2013 UTC (10 years, 4 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.11: +4 -3 lines
Log Message:
Changed secant correction to patch centers for better normalization

File Contents

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