ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2klems.c
Revision: 2.8
Committed: Fri Aug 2 20:56:19 2013 UTC (10 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.7: +11 -1 lines
Log Message:
Added ability to use Dx, Dy and Dz instead of last 3 function parameters

File Contents

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