ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2ttree.c
Revision: 2.27
Committed: Wed Mar 12 22:24:59 2014 UTC (10 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.26: +12 -2 lines
Log Message:
Changed progress bar to erase itself when done

File Contents

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