ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2klems.c
Revision: 2.13
Committed: Wed Mar 12 21:15:31 2014 UTC (10 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.12: +51 -8 lines
Log Message:
Disabled DEBUG statements in BSDF programs and added progress bars

File Contents

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