ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2ttree.c
Revision: 2.31
Committed: Fri Jan 29 16:21:55 2016 UTC (8 years, 3 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.30: +3 -3 lines
Log Message:
Updated pabopto2bsdf to handle tristimulus color -- change in .sir format

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: bsdf2ttree.c,v 2.30 2015/05/05 22:16:49 greg Exp $";
3 #endif
4 /*
5 * Load measured BSDF interpolant and write out as XML file with tensor tree.
6 *
7 * G. Ward
8 */
9
10 #define _USE_MATH_DEFINES
11 #include <stdio.h>
12 #include <stdlib.h>
13 #include <math.h>
14 #include "random.h"
15 #include "platform.h"
16 #include "rtprocess.h"
17 #include "calcomp.h"
18 #include "bsdfrep.h"
19 /* global argv[0] */
20 char *progname;
21 /* percentage to cull (<0 to turn off) */
22 static double pctcull = 90.;
23 /* sampling order */
24 static int samp_order = 6;
25 /* super-sampling threshold */
26 const double ssamp_thresh = 0.35;
27 /* number of super-samples */
28 #ifndef NSSAMP
29 #define NSSAMP 100
30 #endif
31 /* limit on number of RBF lobes */
32 static int lobe_lim = 15000;
33 /* progress bar length */
34 static int do_prog = 79;
35
36
37 /* Start new progress bar */
38 #define prog_start(s) if (do_prog) fprintf(stderr, "%s: %s...\n", progname, s); else
39
40 /* Draw progress bar of the appropriate length */
41 static void
42 prog_show(double frac)
43 {
44 static unsigned call_cnt = 0;
45 static char lastc[] = "-\\|/";
46 char pbar[256];
47 int nchars;
48
49 if (do_prog <= 1) return;
50 if (do_prog > sizeof(pbar)-2)
51 do_prog = sizeof(pbar)-2;
52 if (frac < 0) frac = 0;
53 else if (frac >= 1) frac = .9999;
54 nchars = do_prog*frac;
55 pbar[0] = '\r';
56 memset(pbar+1, '*', nchars);
57 pbar[nchars+1] = lastc[call_cnt++ & 3];
58 memset(pbar+2+nchars, '-', do_prog-nchars-1);
59 pbar[do_prog+1] = '\0';
60 fputs(pbar, stderr);
61 }
62
63 /* Finish progress bar */
64 static void
65 prog_done(void)
66 {
67 int n = do_prog;
68
69 if (n <= 1) return;
70 fputc('\r', stderr);
71 while (n--)
72 fputc(' ', stderr);
73 fputc('\r', stderr);
74 }
75
76 /* Output XML prologue to stdout */
77 static void
78 xml_prologue(int ac, char *av[])
79 {
80 puts("<?xml version=\"1.0\" encoding=\"UTF-8\"?>");
81 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\">");
82 fputs("<!-- File produced by:", stdout);
83 while (ac-- > 0) {
84 fputc(' ', stdout);
85 fputs(*av++, stdout);
86 }
87 puts(" -->");
88 puts("<WindowElementType>System</WindowElementType>");
89 puts("<FileType>BSDF</FileType>");
90 puts("<Optical>");
91 puts("<Layer>");
92 puts("\t<Material>");
93 printf("\t\t<Name>%s</Name>\n", bsdf_name[0] ? bsdf_name : "Unknown");
94 printf("\t\t<Manufacturer>%s</Manufacturer>\n",
95 bsdf_manuf[0] ? bsdf_manuf : "Unknown");
96 puts("\t\t<DeviceType>Other</DeviceType>");
97 puts("\t</Material>");
98 puts("\t<DataDefinition>");
99 printf("\t\t<IncidentDataStructure>TensorTree%c</IncidentDataStructure>\n",
100 single_plane_incident ? '3' : '4');
101 puts("\t</DataDefinition>");
102 }
103
104 /* Output XML data prologue to stdout */
105 static void
106 data_prologue()
107 {
108 static const char *bsdf_type[4] = {
109 "Reflection Front",
110 "Transmission Front",
111 "Transmission Back",
112 "Reflection Back"
113 };
114
115 puts("\t<WavelengthData>");
116 puts("\t\t<LayerNumber>System</LayerNumber>");
117 puts("\t\t<Wavelength unit=\"Integral\">Visible</Wavelength>");
118 puts("\t\t<SourceSpectrum>CIE Illuminant D65 1nm.ssp</SourceSpectrum>");
119 puts("\t\t<DetectorSpectrum>ASTM E308 1931 Y.dsp</DetectorSpectrum>");
120 puts("\t\t<WavelengthDataBlock>");
121 printf("\t\t\t<WavelengthDataDirection>%s</WavelengthDataDirection>\n",
122 bsdf_type[(input_orient>0)<<1 | (output_orient>0)]);
123 puts("\t\t\t<AngleBasis>LBNL/Shirley-Chiu</AngleBasis>");
124 puts("\t\t\t<ScatteringDataType>BTDF</ScatteringDataType>");
125 puts("\t\t\t<ScatteringData>");
126 }
127
128 /* Output XML data epilogue to stdout */
129 static void
130 data_epilogue(void)
131 {
132 puts("\t\t\t</ScatteringData>");
133 puts("\t\t</WavelengthDataBlock>");
134 puts("\t</WavelengthData>");
135 }
136
137 /* Output XML epilogue to stdout */
138 static void
139 xml_epilogue(void)
140 {
141 puts("</Layer>");
142 puts("</Optical>");
143 puts("</WindowElement>");
144 }
145
146 /* Compute absolute relative difference */
147 static double
148 abs_diff(double v1, double v0)
149 {
150 if ((v0 < 0) | (v1 < 0))
151 return(.0);
152 v1 = (v1-v0)*2./(v0+v1+.0001);
153 if (v1 < 0)
154 return(-v1);
155 return(v1);
156 }
157
158 /* Interpolate and output isotropic BSDF data */
159 static void
160 eval_isotropic(char *funame)
161 {
162 const int sqres = 1<<samp_order;
163 FILE *ofp = NULL;
164 int assignD = 0;
165 char cmd[128];
166 int ix, ox, oy;
167 double iovec[6];
168 float bsdf;
169
170 data_prologue(); /* begin output */
171 if (pctcull >= 0) {
172 sprintf(cmd, "rttree_reduce -a -h -ff -r 3 -t %f -g %d",
173 pctcull, samp_order);
174 fflush(stdout);
175 ofp = popen(cmd, "w");
176 if (ofp == NULL) {
177 fprintf(stderr, "%s: cannot create pipe to rttree_reduce\n",
178 progname);
179 exit(1);
180 }
181 SET_FILE_BINARY(ofp);
182 #ifdef getc_unlocked /* avoid lock/unlock overhead */
183 flockfile(ofp);
184 #endif
185 } else
186 fputs("{\n", stdout);
187 /* need to assign Dx, Dy, Dz? */
188 if (funame != NULL)
189 assignD = (fundefined(funame) < 6);
190 /* run through directions */
191 for (ix = 0; ix < sqres/2; ix++) {
192 RBFNODE *rbf = NULL;
193 iovec[0] = 2.*(ix+.5)/sqres - 1.;
194 iovec[1] = .0;
195 iovec[2] = input_orient * sqrt(1. - iovec[0]*iovec[0]);
196 if (funame == NULL)
197 rbf = advect_rbf(iovec, lobe_lim);
198 for (ox = 0; ox < sqres; ox++) {
199 float last_bsdf = -1;
200 for (oy = 0; oy < sqres; oy++) {
201 SDsquare2disk(iovec+3, (ox+.5)/sqres, (oy+.5)/sqres);
202 iovec[5] = output_orient *
203 sqrt(1. - iovec[3]*iovec[3] - iovec[4]*iovec[4]);
204 if (funame == NULL)
205 bsdf = eval_rbfrep(rbf, iovec+3);
206 else {
207 if (assignD) {
208 varset("Dx", '=', -iovec[3]);
209 varset("Dy", '=', -iovec[4]);
210 varset("Dz", '=', -iovec[5]);
211 ++eclock;
212 }
213 bsdf = funvalue(funame, 6, iovec);
214 #if (NSSAMP > 0)
215 if (abs_diff(bsdf, last_bsdf) > ssamp_thresh) {
216 int ssi;
217 double ssa[3], ssvec[6], sum = 0;
218 /* super-sample voxel */
219 for (ssi = NSSAMP; ssi--; ) {
220 SDmultiSamp(ssa, 3, (ssi+frandom()) *
221 (1./NSSAMP));
222 ssvec[0] = 2.*(ix+ssa[0])/sqres - 1.;
223 ssvec[1] = .0;
224 ssvec[2] = input_orient *
225 sqrt(1. - ssvec[0]*ssvec[0]);
226 SDsquare2disk(ssvec+3, (ox+ssa[1])/sqres,
227 (oy+ssa[2])/sqres);
228 ssvec[5] = output_orient *
229 sqrt(1. - ssvec[3]*ssvec[3] -
230 ssvec[4]*ssvec[4]);
231 if (assignD) {
232 varset("Dx", '=', -ssvec[3]);
233 varset("Dy", '=', -ssvec[4]);
234 varset("Dz", '=', -ssvec[5]);
235 ++eclock;
236 }
237 sum += funvalue(funame, 6, ssvec);
238 }
239 bsdf = sum/NSSAMP;
240 }
241 #endif
242 }
243 if (pctcull >= 0)
244 fwrite(&bsdf, sizeof(bsdf), 1, ofp);
245 else
246 printf("\t%.3e\n", bsdf);
247 last_bsdf = bsdf;
248 }
249 }
250 if (rbf != NULL)
251 free(rbf);
252 prog_show((ix+1.)*(2./sqres));
253 }
254 if (pctcull >= 0) { /* finish output */
255 if (pclose(ofp)) {
256 fprintf(stderr, "%s: error running '%s'\n",
257 progname, cmd);
258 exit(1);
259 }
260 } else {
261 for (ix = sqres*sqres*sqres/2; ix--; )
262 fputs("\t0\n", stdout);
263 fputs("}\n", stdout);
264 }
265 data_epilogue();
266 prog_done();
267 }
268
269 /* Interpolate and output anisotropic BSDF data */
270 static void
271 eval_anisotropic(char *funame)
272 {
273 const int sqres = 1<<samp_order;
274 FILE *ofp = NULL;
275 int assignD = 0;
276 char cmd[128];
277 int ix, iy, ox, oy;
278 double iovec[6];
279 float bsdf;
280
281 data_prologue(); /* begin output */
282 if (pctcull >= 0) {
283 sprintf(cmd, "rttree_reduce%s -h -ff -r 4 -t %f -g %d",
284 (input_orient>0 ^ output_orient>0) ? "" : " -a",
285 pctcull, samp_order);
286 fflush(stdout);
287 ofp = popen(cmd, "w");
288 if (ofp == NULL) {
289 fprintf(stderr, "%s: cannot create pipe to rttree_reduce\n",
290 progname);
291 exit(1);
292 }
293 SET_FILE_BINARY(ofp);
294 #ifdef getc_unlocked /* avoid lock/unlock overhead */
295 flockfile(ofp);
296 #endif
297 } else
298 fputs("{\n", stdout);
299 /* need to assign Dx, Dy, Dz? */
300 if (funame != NULL)
301 assignD = (fundefined(funame) < 6);
302 /* run through directions */
303 for (ix = 0; ix < sqres; ix++)
304 for (iy = 0; iy < sqres; iy++) {
305 RBFNODE *rbf = NULL; /* Klems reversal */
306 SDsquare2disk(iovec, 1.-(ix+.5)/sqres, 1.-(iy+.5)/sqres);
307 iovec[2] = input_orient *
308 sqrt(1. - iovec[0]*iovec[0] - iovec[1]*iovec[1]);
309 if (funame == NULL)
310 rbf = advect_rbf(iovec, lobe_lim);
311 for (ox = 0; ox < sqres; ox++) {
312 float last_bsdf = -1;
313 for (oy = 0; oy < sqres; oy++) {
314 SDsquare2disk(iovec+3, (ox+.5)/sqres, (oy+.5)/sqres);
315 iovec[5] = output_orient *
316 sqrt(1. - iovec[3]*iovec[3] - iovec[4]*iovec[4]);
317 if (funame == NULL)
318 bsdf = eval_rbfrep(rbf, iovec+3);
319 else {
320 if (assignD) {
321 varset("Dx", '=', -iovec[3]);
322 varset("Dy", '=', -iovec[4]);
323 varset("Dz", '=', -iovec[5]);
324 ++eclock;
325 }
326 bsdf = funvalue(funame, 6, iovec);
327 #if (NSSAMP > 0)
328 if (abs_diff(bsdf, last_bsdf) > ssamp_thresh) {
329 int ssi;
330 double ssa[4], ssvec[6], sum = 0;
331 /* super-sample voxel */
332 for (ssi = NSSAMP; ssi--; ) {
333 SDmultiSamp(ssa, 4, (ssi+frandom()) *
334 (1./NSSAMP));
335 SDsquare2disk(ssvec, 1.-(ix+ssa[0])/sqres,
336 1.-(iy+ssa[1])/sqres);
337 ssvec[2] = input_orient *
338 sqrt(1. - ssvec[0]*ssvec[0] -
339 ssvec[1]*ssvec[1]);
340 SDsquare2disk(ssvec+3, (ox+ssa[2])/sqres,
341 (oy+ssa[3])/sqres);
342 ssvec[5] = output_orient *
343 sqrt(1. - ssvec[3]*ssvec[3] -
344 ssvec[4]*ssvec[4]);
345 if (assignD) {
346 varset("Dx", '=', -ssvec[3]);
347 varset("Dy", '=', -ssvec[4]);
348 varset("Dz", '=', -ssvec[5]);
349 ++eclock;
350 }
351 sum += funvalue(funame, 6, ssvec);
352 }
353 bsdf = sum/NSSAMP;
354 }
355 #endif
356 }
357 if (pctcull >= 0)
358 fwrite(&bsdf, sizeof(bsdf), 1, ofp);
359 else
360 printf("\t%.3e\n", bsdf);
361 last_bsdf = bsdf;
362 }
363 }
364 if (rbf != NULL)
365 free(rbf);
366 prog_show((ix*sqres+iy+1.)/(sqres*sqres));
367 }
368 if (pctcull >= 0) { /* finish output */
369 if (pclose(ofp)) {
370 fprintf(stderr, "%s: error running '%s'\n",
371 progname, cmd);
372 exit(1);
373 }
374 } else
375 fputs("}\n", stdout);
376 data_epilogue();
377 prog_done();
378 }
379
380 /* Read in BSDF and interpolate as tensor tree representation */
381 int
382 main(int argc, char *argv[])
383 {
384 int dofwd = 0, dobwd = 1;
385 int i, na;
386
387 progname = argv[0];
388 esupport |= E_VARIABLE|E_FUNCTION|E_RCONST;
389 esupport &= ~(E_INCHAN|E_OUTCHAN);
390 scompile("PI:3.14159265358979323846", NULL, 0);
391 biggerlib();
392 for (i = 1; i < argc-1 && (argv[i][0] == '-') | (argv[i][0] == '+'); i++)
393 switch (argv[i][1]) { /* get options */
394 case 'e':
395 scompile(argv[++i], NULL, 0);
396 break;
397 case 'f':
398 if (!argv[i][2])
399 fcompile(argv[++i]);
400 else
401 dofwd = (argv[i][0] == '+');
402 break;
403 case 'b':
404 dobwd = (argv[i][0] == '+');
405 break;
406 case 't':
407 switch (argv[i][2]) {
408 case '3':
409 single_plane_incident = 1;
410 break;
411 case '4':
412 single_plane_incident = 0;
413 break;
414 case '\0':
415 pctcull = atof(argv[++i]);
416 break;
417 default:
418 goto userr;
419 }
420 break;
421 case 'g':
422 samp_order = atoi(argv[++i]);
423 break;
424 case 'l':
425 lobe_lim = atoi(argv[++i]);
426 break;
427 case 'p':
428 do_prog = atoi(argv[i]+2);
429 break;
430 default:
431 goto userr;
432 }
433 if (single_plane_incident >= 0) { /* function-based BSDF? */
434 void (*evf)(char *s) = single_plane_incident ?
435 &eval_isotropic : &eval_anisotropic;
436 if (i != argc-1 || fundefined(argv[i]) < 3) {
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_prologue(argc, argv); /* start XML output */
445 if (dofwd) {
446 input_orient = -1;
447 output_orient = -1;
448 prog_start("Evaluating outside reflectance");
449 (*evf)(argv[i]);
450 output_orient = 1;
451 prog_start("Evaluating outside->inside transmission");
452 (*evf)(argv[i]);
453 }
454 if (dobwd) {
455 input_orient = 1;
456 output_orient = 1;
457 prog_start("Evaluating inside reflectance");
458 (*evf)(argv[i]);
459 output_orient = -1;
460 prog_start("Evaluating inside->outside transmission");
461 (*evf)(argv[i]);
462 }
463 xml_epilogue(); /* finish XML output & exit */
464 return(0);
465 }
466 if (i < argc) { /* open input files if given */
467 int nbsdf = 0;
468 for ( ; i < argc; i++) { /* interpolate each component */
469 char pbuf[256];
470 FILE *fpin = fopen(argv[i], "rb");
471 if (fpin == NULL) {
472 fprintf(stderr, "%s: cannot open BSDF interpolant '%s'\n",
473 progname, argv[i]);
474 return(1);
475 }
476 if (!load_bsdf_rep(fpin))
477 return(1);
478 fclose(fpin);
479 if (!nbsdf++) /* start XML on first dist. */
480 xml_prologue(argc, argv);
481 sprintf(pbuf, "Interpolating component '%s'", argv[i]);
482 prog_start(pbuf);
483 if (single_plane_incident)
484 eval_isotropic(NULL);
485 else
486 eval_anisotropic(NULL);
487 }
488 xml_epilogue(); /* finish XML output & exit */
489 return(0);
490 }
491 SET_FILE_BINARY(stdin); /* load from stdin */
492 if (!load_bsdf_rep(stdin))
493 return(1);
494 xml_prologue(argc, argv); /* start XML output */
495 prog_start("Interpolating from standard input");
496 if (single_plane_incident) /* resample dist. */
497 eval_isotropic(NULL);
498 else
499 eval_anisotropic(NULL);
500 xml_epilogue(); /* finish XML output & exit */
501 return(0);
502 userr:
503 fprintf(stderr,
504 "Usage: %s [-g Nlog2][-t pctcull][-l maxlobes] [bsdf.sir ..] > bsdf.xml\n",
505 progname);
506 fprintf(stderr,
507 " or: %s -t{3|4} [-g Nlog2][-t pctcull][{+|-}for[ward]][{+|-}b[ackward]][-e expr][-f file] bsdf_func > bsdf.xml\n",
508 progname);
509 return(1);
510 }