ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2ttree.c
Revision: 2.21
Committed: Sat Nov 9 05:47:49 2013 UTC (10 years, 4 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.20: +4 -3 lines
Log Message:
Added tracking of BSDF sample name

File Contents

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