ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2ttree.c
Revision: 2.26
Committed: Wed Mar 12 21:15:31 2014 UTC (10 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.25: +49 -5 lines
Log Message:
Disabled DEBUG statements in BSDF programs and added progress bars

File Contents

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