ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2klems.c
Revision: 2.17
Committed: Tue May 5 22:16:49 2015 UTC (9 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.16: +9 -6 lines
Log Message:
Improved progress bar feedback

File Contents

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