ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2klems.c
Revision: 2.14
Committed: Wed Mar 12 22:24:59 2014 UTC (10 years ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R2
Changes since 2.13: +13 -3 lines
Log Message:
Changed progress bar to erase itself when done

File Contents

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