ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2ttree.c
Revision: 2.15
Committed: Wed May 15 03:22:31 2013 UTC (10 years, 10 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.14: +75 -21 lines
Log Message:
Added super-sampling at peaks

File Contents

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