ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2ttree.c
Revision: 2.29
Committed: Thu Aug 21 10:33:48 2014 UTC (10 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R2P2, rad4R2P1
Changes since 2.28: +3 -5 lines
Log Message:
Added grazing angle extrapolation to BSDF interpolation

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: bsdf2ttree.c,v 2.28 2014/03/24 03:50:28 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 else {
204 if (assignD) {
205 varset("Dx", '=', -iovec[3]);
206 varset("Dy", '=', -iovec[4]);
207 varset("Dz", '=', -iovec[5]);
208 ++eclock;
209 }
210 bsdf = funvalue(funame, 6, iovec);
211 #if (NSSAMP > 0)
212 if (abs_diff(bsdf, last_bsdf) > ssamp_thresh) {
213 int ssi;
214 double ssa[3], ssvec[6], sum = 0;
215 /* super-sample voxel */
216 for (ssi = NSSAMP; ssi--; ) {
217 SDmultiSamp(ssa, 3, (ssi+frandom()) *
218 (1./NSSAMP));
219 ssvec[0] = 2.*(ix+ssa[0])/sqres - 1.;
220 ssvec[1] = .0;
221 ssvec[2] = input_orient *
222 sqrt(1. - ssvec[0]*ssvec[0]);
223 SDsquare2disk(ssvec+3, (ox+ssa[1])/sqres,
224 (oy+ssa[2])/sqres);
225 ssvec[5] = output_orient *
226 sqrt(1. - ssvec[3]*ssvec[3] -
227 ssvec[4]*ssvec[4]);
228 if (assignD) {
229 varset("Dx", '=', -ssvec[3]);
230 varset("Dy", '=', -ssvec[4]);
231 varset("Dz", '=', -ssvec[5]);
232 ++eclock;
233 }
234 sum += funvalue(funame, 6, ssvec);
235 }
236 bsdf = sum/NSSAMP;
237 }
238 #endif
239 }
240 if (pctcull >= 0)
241 fwrite(&bsdf, sizeof(bsdf), 1, ofp);
242 else
243 printf("\t%.3e\n", bsdf);
244 last_bsdf = bsdf;
245 }
246 }
247 if (rbf != NULL)
248 free(rbf);
249 prog_show((ix+1.)*(2./sqres));
250 }
251 if (pctcull >= 0) { /* finish output */
252 if (pclose(ofp)) {
253 fprintf(stderr, "%s: error running '%s'\n",
254 progname, cmd);
255 exit(1);
256 }
257 } else {
258 for (ix = sqres*sqres*sqres/2; ix--; )
259 fputs("\t0\n", stdout);
260 fputs("}\n", stdout);
261 }
262 data_epilogue();
263 prog_done();
264 }
265
266 /* Interpolate and output anisotropic BSDF data */
267 static void
268 eval_anisotropic(char *funame)
269 {
270 const int sqres = 1<<samp_order;
271 FILE *ofp = NULL;
272 int assignD = 0;
273 char cmd[128];
274 int ix, iy, ox, oy;
275 double iovec[6];
276 float bsdf;
277
278 data_prologue(); /* begin output */
279 if (pctcull >= 0) {
280 sprintf(cmd, "rttree_reduce%s -h -ff -r 4 -t %f -g %d",
281 (input_orient>0 ^ output_orient>0) ? "" : " -a",
282 pctcull, samp_order);
283 fflush(stdout);
284 ofp = popen(cmd, "w");
285 if (ofp == NULL) {
286 fprintf(stderr, "%s: cannot create pipe to rttree_reduce\n",
287 progname);
288 exit(1);
289 }
290 SET_FILE_BINARY(ofp);
291 #ifdef getc_unlocked /* avoid lock/unlock overhead */
292 flockfile(ofp);
293 #endif
294 } else
295 fputs("{\n", stdout);
296 /* need to assign Dx, Dy, Dz? */
297 if (funame != NULL)
298 assignD = (fundefined(funame) < 6);
299 /* run through directions */
300 for (ix = 0; ix < sqres; ix++)
301 for (iy = 0; iy < sqres; iy++) {
302 RBFNODE *rbf = NULL; /* Klems reversal */
303 SDsquare2disk(iovec, 1.-(ix+.5)/sqres, 1.-(iy+.5)/sqres);
304 iovec[2] = input_orient *
305 sqrt(1. - iovec[0]*iovec[0] - iovec[1]*iovec[1]);
306 if (funame == NULL)
307 rbf = advect_rbf(iovec, lobe_lim);
308 for (ox = 0; ox < sqres; ox++) {
309 float last_bsdf = -1;
310 for (oy = 0; oy < sqres; oy++) {
311 SDsquare2disk(iovec+3, (ox+.5)/sqres, (oy+.5)/sqres);
312 iovec[5] = output_orient *
313 sqrt(1. - iovec[3]*iovec[3] - iovec[4]*iovec[4]);
314 if (funame == NULL)
315 bsdf = eval_rbfrep(rbf, iovec+3);
316 else {
317 if (assignD) {
318 varset("Dx", '=', -iovec[3]);
319 varset("Dy", '=', -iovec[4]);
320 varset("Dz", '=', -iovec[5]);
321 ++eclock;
322 }
323 bsdf = funvalue(funame, 6, iovec);
324 #if (NSSAMP > 0)
325 if (abs_diff(bsdf, last_bsdf) > ssamp_thresh) {
326 int ssi;
327 double ssa[4], ssvec[6], sum = 0;
328 /* super-sample voxel */
329 for (ssi = NSSAMP; ssi--; ) {
330 SDmultiSamp(ssa, 4, (ssi+frandom()) *
331 (1./NSSAMP));
332 SDsquare2disk(ssvec, 1.-(ix+ssa[0])/sqres,
333 1.-(iy+ssa[1])/sqres);
334 ssvec[2] = input_orient *
335 sqrt(1. - ssvec[0]*ssvec[0] -
336 ssvec[1]*ssvec[1]);
337 SDsquare2disk(ssvec+3, (ox+ssa[2])/sqres,
338 (oy+ssa[3])/sqres);
339 ssvec[5] = output_orient *
340 sqrt(1. - ssvec[3]*ssvec[3] -
341 ssvec[4]*ssvec[4]);
342 if (assignD) {
343 varset("Dx", '=', -ssvec[3]);
344 varset("Dy", '=', -ssvec[4]);
345 varset("Dz", '=', -ssvec[5]);
346 ++eclock;
347 }
348 sum += funvalue(funame, 6, ssvec);
349 }
350 bsdf = sum/NSSAMP;
351 }
352 #endif
353 }
354 if (pctcull >= 0)
355 fwrite(&bsdf, sizeof(bsdf), 1, ofp);
356 else
357 printf("\t%.3e\n", bsdf);
358 last_bsdf = bsdf;
359 }
360 }
361 if (rbf != NULL)
362 free(rbf);
363 prog_show((ix*sqres+iy+1.)/(sqres*sqres));
364 }
365 if (pctcull >= 0) { /* finish output */
366 if (pclose(ofp)) {
367 fprintf(stderr, "%s: error running '%s'\n",
368 progname, cmd);
369 exit(1);
370 }
371 } else
372 fputs("}\n", stdout);
373 data_epilogue();
374 prog_done();
375 }
376
377 /* Read in BSDF and interpolate as tensor tree representation */
378 int
379 main(int argc, char *argv[])
380 {
381 int dofwd = 0, dobwd = 1;
382 int i, na;
383
384 progname = argv[0];
385 esupport |= E_VARIABLE|E_FUNCTION|E_RCONST;
386 esupport &= ~(E_INCHAN|E_OUTCHAN);
387 scompile("PI:3.14159265358979323846", NULL, 0);
388 biggerlib();
389 for (i = 1; i < argc-1 && (argv[i][0] == '-') | (argv[i][0] == '+'); i++)
390 switch (argv[i][1]) { /* get options */
391 case 'e':
392 scompile(argv[++i], NULL, 0);
393 break;
394 case 'f':
395 if (!argv[i][2])
396 fcompile(argv[++i]);
397 else
398 dofwd = (argv[i][0] == '+');
399 break;
400 case 'b':
401 dobwd = (argv[i][0] == '+');
402 break;
403 case 't':
404 switch (argv[i][2]) {
405 case '3':
406 single_plane_incident = 1;
407 break;
408 case '4':
409 single_plane_incident = 0;
410 break;
411 case '\0':
412 pctcull = atof(argv[++i]);
413 break;
414 default:
415 goto userr;
416 }
417 break;
418 case 'g':
419 samp_order = atoi(argv[++i]);
420 break;
421 case 'l':
422 lobe_lim = atoi(argv[++i]);
423 break;
424 case 'p':
425 do_prog = atoi(argv[i]+2);
426 break;
427 default:
428 goto userr;
429 }
430 if (single_plane_incident >= 0) { /* function-based BSDF? */
431 void (*evf)(char *s) = single_plane_incident ?
432 &eval_isotropic : &eval_anisotropic;
433 if (i != argc-1 || fundefined(argv[i]) < 3) {
434 fprintf(stderr,
435 "%s: need single function with 6 arguments: bsdf(ix,iy,iz,ox,oy,oz)\n",
436 progname);
437 fprintf(stderr, "\tor 3 arguments using Dx,Dy,Dz: bsdf(ix,iy,iz)\n");
438 goto userr;
439 }
440 ++eclock;
441 xml_prologue(argc, argv); /* start XML output */
442 if (dofwd) {
443 input_orient = -1;
444 output_orient = -1;
445 prog_start("Evaluating outside reflectance");
446 (*evf)(argv[i]);
447 output_orient = 1;
448 prog_start("Evaluating outside->inside transmission");
449 (*evf)(argv[i]);
450 }
451 if (dobwd) {
452 input_orient = 1;
453 output_orient = 1;
454 prog_start("Evaluating inside reflectance");
455 (*evf)(argv[i]);
456 output_orient = -1;
457 prog_start("Evaluating inside->outside transmission");
458 (*evf)(argv[i]);
459 }
460 xml_epilogue(); /* finish XML output & exit */
461 return(0);
462 }
463 if (i < argc) { /* open input files if given */
464 int nbsdf = 0;
465 for ( ; i < argc; i++) { /* interpolate each component */
466 char pbuf[256];
467 FILE *fpin = fopen(argv[i], "rb");
468 if (fpin == NULL) {
469 fprintf(stderr, "%s: cannot open BSDF interpolant '%s'\n",
470 progname, argv[i]);
471 return(1);
472 }
473 if (!load_bsdf_rep(fpin))
474 return(1);
475 fclose(fpin);
476 if (!nbsdf++) /* start XML on first dist. */
477 xml_prologue(argc, argv);
478 sprintf(pbuf, "Interpolating component '%s'", argv[i]);
479 prog_start(pbuf);
480 if (single_plane_incident)
481 eval_isotropic(NULL);
482 else
483 eval_anisotropic(NULL);
484 }
485 xml_epilogue(); /* finish XML output & exit */
486 return(0);
487 }
488 SET_FILE_BINARY(stdin); /* load from stdin */
489 if (!load_bsdf_rep(stdin))
490 return(1);
491 xml_prologue(argc, argv); /* start XML output */
492 prog_start("Interpolating from standard input");
493 if (single_plane_incident) /* resample dist. */
494 eval_isotropic(NULL);
495 else
496 eval_anisotropic(NULL);
497 xml_epilogue(); /* finish XML output & exit */
498 return(0);
499 userr:
500 fprintf(stderr,
501 "Usage: %s [-g Nlog2][-t pctcull][-l maxlobes] [bsdf.sir ..] > bsdf.xml\n",
502 progname);
503 fprintf(stderr,
504 " or: %s -t{3|4} [-g Nlog2][-t pctcull][{+|-}for[ward]][{+|-}b[ackward]][-e expr][-f file] bsdf_func > bsdf.xml\n",
505 progname);
506 return(1);
507 }