ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2klems.c
Revision: 2.18
Committed: Wed May 27 11:39:37 2015 UTC (8 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R0
Changes since 2.17: +2 -1 lines
Log Message:
Added missing initialization call

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: bsdf2klems.c,v 2.17 2015/05/05 22:16:49 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 initurand(npsamps);
211 SDclearBSDF(&bsd, fname); /* load BSDF file */
212 if ((ec = SDloadFile(&bsd, fname)) != SDEnone)
213 goto err;
214 xml_prologue(&bsd); /* pass geometry */
215 /* front reflection */
216 if (bsd.rf != NULL || bsd.rLambFront.cieY > .002) {
217 input_orient = 1; output_orient = 1;
218 data_prologue();
219 for (j = 0; j < abp->nangles; j++) {
220 for (i = 0; i < abp->nangles; i++) {
221 sum = 0; /* average over patches */
222 for (n = npsamps; n-- > 0; ) {
223 fo_getvec(vout, j+(n+frandom())/npsamps, abp);
224 fi_getvec(vin, i+urand(n), abp);
225 ec = SDevalBSDF(&sv, vout, vin, &bsd);
226 if (ec != SDEnone)
227 goto err;
228 sum += sv.cieY;
229 }
230 printf("\t%.3e\n", sum/npsamps);
231 }
232 putchar('\n'); /* extra space between rows */
233 }
234 data_epilogue();
235 }
236 /* back reflection */
237 if (bsd.rb != NULL || bsd.rLambBack.cieY > .002) {
238 input_orient = -1; output_orient = -1;
239 data_prologue();
240 for (j = 0; j < abp->nangles; j++) {
241 for (i = 0; i < abp->nangles; i++) {
242 sum = 0; /* average over patches */
243 for (n = npsamps; n-- > 0; ) {
244 bo_getvec(vout, j+(n+frandom())/npsamps, abp);
245 bi_getvec(vin, i+urand(n), abp);
246 ec = SDevalBSDF(&sv, vout, vin, &bsd);
247 if (ec != SDEnone)
248 goto err;
249 sum += sv.cieY;
250 }
251 printf("\t%.3e\n", sum/npsamps);
252 }
253 putchar('\n'); /* extra space between rows */
254 }
255 data_epilogue();
256 }
257 /* front transmission */
258 if (bsd.tf != NULL || bsd.tLamb.cieY > .002) {
259 input_orient = 1; output_orient = -1;
260 data_prologue();
261 for (j = 0; j < abp->nangles; j++) {
262 for (i = 0; i < abp->nangles; i++) {
263 sum = 0; /* average over patches */
264 for (n = npsamps; n-- > 0; ) {
265 bo_getvec(vout, j+(n+frandom())/npsamps, abp);
266 fi_getvec(vin, i+urand(n), abp);
267 ec = SDevalBSDF(&sv, vout, vin, &bsd);
268 if (ec != SDEnone)
269 goto err;
270 sum += sv.cieY;
271 }
272 printf("\t%.3e\n", sum/npsamps);
273 }
274 putchar('\n'); /* extra space between rows */
275 }
276 data_epilogue();
277 }
278 /* back transmission */
279 if ((bsd.tb != NULL) | (bsd.tf != NULL)) {
280 input_orient = -1; output_orient = 1;
281 data_prologue();
282 for (j = 0; j < abp->nangles; j++) {
283 for (i = 0; i < abp->nangles; i++) {
284 sum = 0; /* average over patches */
285 for (n = npsamps; n-- > 0; ) {
286 fo_getvec(vout, j+(n+frandom())/npsamps, abp);
287 bi_getvec(vin, i+urand(n), abp);
288 ec = SDevalBSDF(&sv, vout, vin, &bsd);
289 if (ec != SDEnone)
290 goto err;
291 sum += sv.cieY;
292 }
293 printf("\t%.3e\n", sum/npsamps);
294 }
295 putchar('\n'); /* extra space between rows */
296 }
297 data_epilogue();
298 }
299 SDfreeBSDF(&bsd); /* all done */
300 return;
301 err:
302 SDreportError(ec, stderr);
303 exit(1);
304 }
305
306 /* Interpolate and output a BSDF function using Klems basis */
307 static void
308 eval_function(char *funame)
309 {
310 ANGLE_BASIS *abp = get_basis(kbasis);
311 int assignD = (fundefined(funame) < 6);
312 double iovec[6];
313 double sum;
314 int i, j, n;
315
316 initurand(npsamps);
317 data_prologue(); /* begin output */
318 for (j = 0; j < abp->nangles; j++) { /* run through directions */
319 for (i = 0; i < abp->nangles; i++) {
320 sum = 0;
321 for (n = npsamps; n--; ) { /* average over patches */
322 if (output_orient > 0)
323 fo_getvec(iovec+3, j+(n+frandom())/npsamps, abp);
324 else
325 bo_getvec(iovec+3, j+(n+frandom())/npsamps, abp);
326
327 if (input_orient > 0)
328 fi_getvec(iovec, i+urand(n), abp);
329 else
330 bi_getvec(iovec, i+urand(n), abp);
331
332 if (assignD) {
333 varset("Dx", '=', -iovec[3]);
334 varset("Dy", '=', -iovec[4]);
335 varset("Dz", '=', -iovec[5]);
336 ++eclock;
337 }
338 sum += funvalue(funame, 6, iovec);
339 }
340 printf("\t%.3e\n", sum/npsamps);
341 }
342 putchar('\n');
343 prog_show((j+1.)/abp->nangles);
344 }
345 data_epilogue(); /* finish output */
346 prog_done();
347 }
348
349 /* Interpolate and output a radial basis function BSDF representation */
350 static void
351 eval_rbf(void)
352 {
353 ANGLE_BASIS *abp = get_basis(kbasis);
354 float bsdfarr[MAXPATCHES*MAXPATCHES];
355 FVECT vin, vout;
356 RBFNODE *rbf;
357 double sum;
358 int i, j, n;
359 /* sanity check */
360 if (abp->nangles > MAXPATCHES) {
361 fprintf(stderr, "%s: too many patches!\n", progname);
362 exit(1);
363 }
364 data_prologue(); /* begin output */
365 for (i = 0; i < abp->nangles; i++) {
366 if (input_orient > 0) /* use incident patch center */
367 fi_getvec(vin, i+.5*(i>0), abp);
368 else
369 bi_getvec(vin, i+.5*(i>0), abp);
370
371 rbf = advect_rbf(vin, lobe_lim); /* compute radial basis func */
372
373 for (j = 0; j < abp->nangles; j++) {
374 sum = 0; /* sample over exiting patch */
375 for (n = npsamps; n--; ) {
376 if (output_orient > 0)
377 fo_getvec(vout, j+(n+frandom())/npsamps, abp);
378 else
379 bo_getvec(vout, j+(n+frandom())/npsamps, abp);
380
381 sum += eval_rbfrep(rbf, vout);
382 }
383 bsdfarr[j*abp->nangles + i] = sum / (double)npsamps;
384 }
385 if (rbf != NULL)
386 free(rbf);
387 prog_show((i+1.)/abp->nangles);
388 }
389 n = 0; /* write out our matrix */
390 for (j = 0; j < abp->nangles; j++) {
391 for (i = 0; i < abp->nangles; i++)
392 printf("\t%.3e\n", bsdfarr[n++]);
393 putchar('\n');
394 }
395 data_epilogue(); /* finish output */
396 prog_done();
397 }
398
399 /* Read in BSDF and interpolate as Klems matrix representation */
400 int
401 main(int argc, char *argv[])
402 {
403 int dofwd = 0, dobwd = 1;
404 char *cp;
405 int i, na;
406
407 progname = argv[0];
408 esupport |= E_VARIABLE|E_FUNCTION|E_RCONST;
409 esupport &= ~(E_INCHAN|E_OUTCHAN);
410 scompile("PI:3.14159265358979323846", NULL, 0);
411 biggerlib();
412 for (i = 1; i < argc && (argv[i][0] == '-') | (argv[i][0] == '+'); i++)
413 switch (argv[i][1]) { /* get options */
414 case 'n':
415 npsamps = atoi(argv[++i]);
416 if (npsamps <= 0)
417 goto userr;
418 break;
419 case 'e':
420 scompile(argv[++i], NULL, 0);
421 single_plane_incident = 0;
422 break;
423 case 'f':
424 if (!argv[i][2]) {
425 fcompile(argv[++i]);
426 single_plane_incident = 0;
427 } else
428 dofwd = (argv[i][0] == '+');
429 break;
430 case 'b':
431 dobwd = (argv[i][0] == '+');
432 break;
433 case 'h':
434 kbasis = "LBNL/Klems Half";
435 break;
436 case 'q':
437 kbasis = "LBNL/Klems Quarter";
438 break;
439 case 'l':
440 lobe_lim = atoi(argv[++i]);
441 break;
442 case 'p':
443 do_prog = atoi(argv[i]+2);
444 break;
445 default:
446 goto userr;
447 }
448 if (single_plane_incident >= 0) { /* function-based BSDF? */
449 if (i != argc-1 || fundefined(argv[i]) != 6) {
450 fprintf(stderr,
451 "%s: need single function with 6 arguments: bsdf(ix,iy,iz,ox,oy,oz)\n",
452 progname);
453 fprintf(stderr, "\tor 3 arguments using Dx,Dy,Dz: bsdf(ix,iy,iz)\n");
454 goto userr;
455 }
456 ++eclock;
457 xml_header(argc, argv); /* start XML output */
458 xml_prologue(NULL);
459 if (dofwd) {
460 input_orient = -1;
461 output_orient = -1;
462 prog_start("Evaluating outside reflectance");
463 eval_function(argv[i]);
464 output_orient = 1;
465 prog_start("Evaluating outside->inside transmission");
466 eval_function(argv[i]);
467 }
468 if (dobwd) {
469 input_orient = 1;
470 output_orient = 1;
471 prog_start("Evaluating inside reflectance");
472 eval_function(argv[i]);
473 output_orient = -1;
474 prog_start("Evaluating inside->outside transmission");
475 eval_function(argv[i]);
476 }
477 xml_epilogue(); /* finish XML output & exit */
478 return(0);
479 }
480 /* XML input? */
481 if (i == argc-1 && (cp = argv[i]+strlen(argv[i])-4) > argv[i] &&
482 !strcasecmp(cp, ".xml")) {
483 xml_header(argc, argv); /* start XML output */
484 eval_bsdf(argv[i]); /* load & resample BSDF */
485 xml_epilogue(); /* finish XML output & exit */
486 return(0);
487 }
488 if (i < argc) { /* open input files if given */
489 int nbsdf = 0;
490 for ( ; i < argc; i++) { /* interpolate each component */
491 char pbuf[256];
492 FILE *fpin = fopen(argv[i], "rb");
493 if (fpin == NULL) {
494 fprintf(stderr, "%s: cannot open BSDF interpolant '%s'\n",
495 progname, argv[i]);
496 return(1);
497 }
498 if (!load_bsdf_rep(fpin))
499 return(1);
500 fclose(fpin);
501 if (!nbsdf++) { /* start XML on first dist. */
502 xml_header(argc, argv);
503 xml_prologue(NULL);
504 }
505 sprintf(pbuf, "Interpolating component '%s'", argv[i]);
506 prog_start(pbuf);
507 eval_rbf();
508 }
509 xml_epilogue(); /* finish XML output & exit */
510 return(0);
511 }
512 SET_FILE_BINARY(stdin); /* load from stdin */
513 if (!load_bsdf_rep(stdin))
514 return(1);
515 xml_header(argc, argv); /* start XML output */
516 xml_prologue(NULL);
517 prog_start("Interpolating from standard input");
518 eval_rbf(); /* resample dist. */
519 xml_epilogue(); /* finish XML output & exit */
520 return(0);
521 userr:
522 fprintf(stderr,
523 "Usage: %s [-n spp][-h|-q][-l maxlobes] [bsdf.sir ..] > bsdf.xml\n", progname);
524 fprintf(stderr,
525 " or: %s [-n spp][-h|-q] bsdf_in.xml > bsdf_out.xml\n", progname);
526 fprintf(stderr,
527 " or: %s [-n spp][-h|-q][{+|-}for[ward]][{+|-}b[ackward]][-e expr][-f file] bsdf_func > bsdf.xml\n",
528 progname);
529 return(1);
530 }