ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2ttree.c
Revision: 2.28
Committed: Mon Mar 24 03:50:28 2014 UTC (10 years ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R2
Changes since 2.27: +29 -21 lines
Log Message:
Fixed one minor and one major bug in super-sampling (major bug in transmission)

File Contents

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