ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2klems.c
Revision: 2.15
Committed: Thu Aug 21 10:33:48 2014 UTC (10 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R2P1
Changes since 2.14: +2 -2 lines
Log Message:
Added grazing angle extrapolation to BSDF interpolation

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: bsdf2klems.c,v 2.14 2014/03/12 22:24:59 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 <= 1) 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 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
70 /* Return angle basis corresponding to the given name */
71 static ANGLE_BASIS *
72 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 /* 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 /* Output XML prologue to stdout */
97 static void
98 xml_prologue(const SDData *sd)
99 {
100 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 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 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 puts("\t\t<DeviceType>Other</DeviceType>");
124 puts("\t</Material>");
125 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 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 printf("\t\t\t<nPhis>%d</nPhis>\n", abp->lat[i].nphis);
142 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 /* Load and resample XML BSDF description using Klems basis */
196 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 xml_prologue(&bsd); /* pass geometry */
211 /* 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 fi_getvec(vin, i+urand(n), abp);
221 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 bi_getvec(vin, i+urand(n), abp);
242 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 fi_getvec(vin, i+urand(n), abp);
263 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 if ((bsd.tb != NULL) | (bsd.tf != NULL)) {
276 input_orient = -1; output_orient = 1;
277 data_prologue();
278 for (j = 0; j < abp->nangles; j++) {
279 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 goto err;
287 sum += sv.cieY;
288 }
289 printf("\t%.3e\n", sum/npsamps);
290 }
291 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 /* Interpolate and output a BSDF function using Klems basis */
303 static void
304 eval_function(char *funame)
305 {
306 ANGLE_BASIS *abp = get_basis(kbasis);
307 int assignD = (fundefined(funame) < 6);
308 double iovec[6];
309 double sum;
310 int i, j, n;
311
312 initurand(npsamps);
313 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 fi_getvec(iovec, i+urand(n), abp);
325 else
326 bi_getvec(iovec, i+urand(n), abp);
327
328 if (assignD) {
329 varset("Dx", '=', -iovec[3]);
330 varset("Dy", '=', -iovec[4]);
331 varset("Dz", '=', -iovec[5]);
332 ++eclock;
333 }
334 sum += funvalue(funame, 6, iovec);
335 }
336 printf("\t%.3e\n", sum/npsamps);
337 }
338 putchar('\n');
339 prog_show((j+1.)/abp->nangles);
340 }
341 data_epilogue(); /* finish output */
342 prog_done();
343 }
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 float bsdfarr[MAXPATCHES*MAXPATCHES];
351 FVECT vin, vout;
352 RBFNODE *rbf;
353 double sum;
354 int i, j, n;
355 /* 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 rbf = advect_rbf(vin, lobe_lim); /* compute radial basis func */
368
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
377 sum += eval_rbfrep(rbf, vout);
378 }
379 fo_getvec(vout, j+.5, abp); /* use centered secant */
380 bsdfarr[j*abp->nangles + i] = sum / (double)npsamps;
381 }
382 if (rbf != NULL)
383 free(rbf);
384 prog_show((i+1.)/abp->nangles);
385 }
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 prog_done();
394 }
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 for (i = 1; i < argc && (argv[i][0] == '-') | (argv[i][0] == '+'); i++)
410 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 case 'l':
437 lobe_lim = atoi(argv[++i]);
438 break;
439 case 'p':
440 do_prog = atoi(argv[i]+2);
441 break;
442 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 fprintf(stderr, "\tor 3 arguments using Dx,Dy,Dz: bsdf(ix,iy,iz)\n");
451 goto userr;
452 }
453 ++eclock;
454 xml_header(argc, argv); /* start XML output */
455 xml_prologue(NULL);
456 if (dofwd) {
457 input_orient = -1;
458 output_orient = -1;
459 prog_start("Evaluating outside reflectance");
460 eval_function(argv[i]);
461 output_orient = 1;
462 prog_start("Evaluating outside->inside transmission");
463 eval_function(argv[i]);
464 }
465 if (dobwd) {
466 input_orient = 1;
467 output_orient = 1;
468 prog_start("Evaluating inside reflectance");
469 eval_function(argv[i]);
470 output_orient = -1;
471 prog_start("Evaluating inside->outside transmission");
472 eval_function(argv[i]);
473 }
474 xml_epilogue(); /* finish XML output & exit */
475 return(0);
476 }
477 /* XML input? */
478 if (i == argc-1 && (cp = argv[i]+strlen(argv[i])-4) > argv[i] &&
479 !strcasecmp(cp, ".xml")) {
480 xml_header(argc, argv); /* start XML output */
481 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 char pbuf[256];
489 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 if (!nbsdf++) { /* start XML on first dist. */
499 xml_header(argc, argv);
500 xml_prologue(NULL);
501 }
502 sprintf(pbuf, "Interpolating component '%s'", argv[i]);
503 prog_start(pbuf);
504 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 xml_header(argc, argv); /* start XML output */
513 xml_prologue(NULL);
514 prog_start("Interpolating from standard input");
515 eval_rbf(); /* resample dist. */
516 xml_epilogue(); /* finish XML output & exit */
517 return(0);
518 userr:
519 fprintf(stderr,
520 "Usage: %s [-n spp][-h|-q][-l maxlobes] [bsdf.sir ..] > bsdf.xml\n", progname);
521 fprintf(stderr,
522 " or: %s [-n spp][-h|-q] bsdf_in.xml > bsdf_out.xml\n", progname);
523 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 }