ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2klems.c
Revision: 2.8
Committed: Fri Aug 2 20:56:19 2013 UTC (10 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.7: +11 -1 lines
Log Message:
Added ability to use Dx, Dy and Dz instead of last 3 function parameters

File Contents

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