ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2ttree.c
Revision: 2.55
Committed: Wed Oct 28 18:54:21 2020 UTC (3 years, 6 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.54: +2 -2 lines
Log Message: 
fix(bsdf2klems, bsdf2ttree): increased command-line buffer to accommodate comments in Windows

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.55 static const char RCSid[] = "$Id: bsdf2ttree.c,v 2.54 2020/10/26 21:28:18 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 <stdlib.h>
12     #include <math.h>
13 greg 2.53 #include <ctype.h>
14 greg 2.23 #include "random.h"
15 greg 2.1 #include "platform.h"
16 greg 2.34 #include "paths.h"
17     #include "rtio.h"
18 greg 2.11 #include "calcomp.h"
19 greg 2.1 #include "bsdfrep.h"
20     /* global argv[0] */
21     char *progname;
22 greg 2.41 /* reciprocity averaging option */
23     static const char *recip = " -a";
24 greg 2.1 /* percentage to cull (<0 to turn off) */
25 greg 2.31 static double pctcull = 90.;
26 greg 2.1 /* sampling order */
27 greg 2.31 static int samp_order = 6;
28 greg 2.15 /* super-sampling threshold */
29 greg 2.48 static double ssamp_thresh = 0.35;
30 greg 2.15 /* number of super-samples */
31 greg 2.52 static int nssamp = 256;
32 greg 2.18 /* limit on number of RBF lobes */
33     static int lobe_lim = 15000;
34 greg 2.26 /* progress bar length */
35     static int do_prog = 79;
36    
37 greg 2.34 #define MAXCARG 512 /* wrapBSDF command */
38     static char *wrapBSDF[MAXCARG] = {"wrapBSDF", "-U"};
39     static int wbsdfac = 2;
40    
41 greg 2.39 /* Add argument to wrapBSDF, allocating space if !isstatic */
42 greg 2.34 static void
43     add_wbsdf(const char *arg, int isstatic)
44     {
45     if (arg == NULL)
46     return;
47     if (wbsdfac >= MAXCARG-1) {
48     fputs(progname, stderr);
49     fputs(": too many command arguments to wrapBSDF\n", stderr);
50     exit(1);
51     }
52     if (!*arg)
53     arg = "";
54     else if (!isstatic)
55     arg = savqstr((char *)arg);
56    
57     wrapBSDF[wbsdfac++] = (char *)arg;
58     }
59    
60     /* Create Yuv component file and add appropriate arguments */
61     static char *
62     create_component_file(int c)
63     {
64     static const char sname[3][6] = {"CIE-Y", "CIE-u", "CIE-v"};
65     static const char cname[4][4] = {"-rf", "-tf", "-tb", "-rb"};
66     char *tfname = mktemp(savqstr(TEMPLATE));
67    
68     add_wbsdf("-s", 1); add_wbsdf(sname[c], 1);
69     add_wbsdf(cname[(input_orient>0)<<1 | (output_orient>0)], 1);
70     add_wbsdf(tfname, 1);
71     return(tfname);
72     }
73 greg 2.26
74     /* Start new progress bar */
75     #define prog_start(s) if (do_prog) fprintf(stderr, "%s: %s...\n", progname, s); else
76    
77     /* Draw progress bar of the appropriate length */
78     static void
79     prog_show(double frac)
80     {
81 greg 2.30 static unsigned call_cnt = 0;
82     static char lastc[] = "-\\|/";
83     char pbar[256];
84     int nchars;
85 greg 2.26
86 greg 2.30 if (do_prog <= 1) return;
87 greg 2.26 if (do_prog > sizeof(pbar)-2)
88     do_prog = sizeof(pbar)-2;
89     if (frac < 0) frac = 0;
90 greg 2.30 else if (frac >= 1) frac = .9999;
91     nchars = do_prog*frac;
92 greg 2.26 pbar[0] = '\r';
93     memset(pbar+1, '*', nchars);
94 greg 2.30 pbar[nchars+1] = lastc[call_cnt++ & 3];
95     memset(pbar+2+nchars, '-', do_prog-nchars-1);
96 greg 2.26 pbar[do_prog+1] = '\0';
97     fputs(pbar, stderr);
98     }
99    
100     /* Finish progress bar */
101 greg 2.27 static void
102     prog_done(void)
103     {
104     int n = do_prog;
105    
106     if (n <= 1) return;
107     fputc('\r', stderr);
108     while (n--)
109     fputc(' ', stderr);
110     fputc('\r', stderr);
111     }
112 greg 2.1
113 greg 2.15 /* Compute absolute relative difference */
114     static double
115     abs_diff(double v1, double v0)
116     {
117     if ((v0 < 0) | (v1 < 0))
118     return(.0);
119     v1 = (v1-v0)*2./(v0+v1+.0001);
120     if (v1 < 0)
121     return(-v1);
122     return(v1);
123     }
124    
125 greg 2.1 /* Interpolate and output isotropic BSDF data */
126     static void
127 greg 2.11 eval_isotropic(char *funame)
128 greg 2.1 {
129     const int sqres = 1<<samp_order;
130 greg 2.42 const double sqfact = 1./(double)sqres;
131 greg 2.50 float *val_last = NULL;
132     float *val_next = NULL;
133     SDValue *sdv_next = NULL;
134 greg 2.34 FILE *ofp, *uvfp[2];
135 greg 2.17 int assignD = 0;
136 greg 2.1 char cmd[128];
137     int ix, ox, oy;
138 greg 2.11 double iovec[6];
139 greg 2.34 float bsdf, uv[2];
140 greg 2.15
141 greg 2.11 if (pctcull >= 0) {
142 greg 2.41 sprintf(cmd, "rttree_reduce%s -h -ff -r 3 -t %f -g %d > %s",
143     recip, pctcull, samp_order, create_component_file(0));
144 greg 2.1 ofp = popen(cmd, "w");
145     if (ofp == NULL) {
146     fprintf(stderr, "%s: cannot create pipe to rttree_reduce\n",
147     progname);
148     exit(1);
149     }
150     SET_FILE_BINARY(ofp);
151 greg 2.25 #ifdef getc_unlocked /* avoid lock/unlock overhead */
152     flockfile(ofp);
153     #endif
154 greg 2.34 if (rbf_colorimetry == RBCtristimulus) {
155     double uvcull = 100. - (100.-pctcull)*.25;
156 greg 2.41 sprintf(cmd, "rttree_reduce%s -h -ff -r 3 -t %f -g %d > %s",
157     recip, uvcull, samp_order, create_component_file(1));
158 greg 2.34 uvfp[0] = popen(cmd, "w");
159 greg 2.41 sprintf(cmd, "rttree_reduce%s -h -ff -r 3 -t %f -g %d > %s",
160     recip, uvcull, samp_order, create_component_file(2));
161 greg 2.34 uvfp[1] = popen(cmd, "w");
162     if ((uvfp[0] == NULL) | (uvfp[1] == NULL)) {
163     fprintf(stderr, "%s: cannot open pipes to uv output\n",
164     progname);
165     exit(1);
166     }
167     SET_FILE_BINARY(uvfp[0]); SET_FILE_BINARY(uvfp[1]);
168     #ifdef getc_unlocked
169     flockfile(uvfp[0]); flockfile(uvfp[1]);
170     #endif
171     }
172     } else {
173     ofp = fopen(create_component_file(0), "w");
174     if (ofp == NULL) {
175     fprintf(stderr, "%s: cannot create Y output file\n",
176     progname);
177     exit(1);
178     }
179     fputs("{\n", ofp);
180     if (rbf_colorimetry == RBCtristimulus) {
181     uvfp[0] = fopen(create_component_file(1), "w");
182     uvfp[1] = fopen(create_component_file(2), "w");
183     if ((uvfp[0] == NULL) | (uvfp[1] == NULL)) {
184     fprintf(stderr, "%s: cannot create uv output file(s)\n",
185     progname);
186     exit(1);
187     }
188     fputs("{\n", uvfp[0]);
189     fputs("{\n", uvfp[1]);
190     }
191     }
192 greg 2.35 if (funame != NULL) /* need to assign Dx, Dy, Dz? */
193 greg 2.17 assignD = (fundefined(funame) < 6);
194 greg 2.50 val_last = (float *)calloc(sqres, sizeof(float));
195     if (funame == NULL)
196     sdv_next = (SDValue *)malloc(sizeof(SDValue)*sqres);
197     else
198     val_next = (float *)malloc(sizeof(float)*sqres);
199 greg 2.1 /* run through directions */
200     for (ix = 0; ix < sqres/2; ix++) {
201 greg 2.45 const int zipsgn = (ix & 1)*2 - 1;
202     RBFNODE *rbf = NULL;
203 greg 2.42 iovec[0] = 2.*sqfact*(ix+.5) - 1.;
204 greg 2.45 iovec[1] = zipsgn*sqfact*.5;
205     iovec[2] = input_orient * sqrt(1. - iovec[0]*iovec[0]
206     - iovec[1]*iovec[1]);
207 greg 2.11 if (funame == NULL)
208 greg 2.18 rbf = advect_rbf(iovec, lobe_lim);
209 greg 2.50 /* presample first row */
210     for (oy = 0; oy < sqres; oy++) {
211     SDsquare2disk(iovec+3, .5*sqfact, (oy+.5)*sqfact);
212 greg 2.47 iovec[5] = output_orient *
213     sqrt(1. - iovec[3]*iovec[3] - iovec[4]*iovec[4]);
214     if (funame == NULL) {
215 greg 2.50 eval_rbfcol(&sdv_next[oy], rbf, iovec+3);
216 greg 2.47 } else {
217     if (assignD) {
218     varset("Dx", '=', -iovec[3]);
219     varset("Dy", '=', -iovec[4]);
220     varset("Dz", '=', -iovec[5]);
221     ++eclock;
222     }
223 greg 2.50 val_next[oy] = funvalue(funame, 6, iovec);
224 greg 2.47 }
225 greg 2.50 }
226     for (ox = 0; ox < sqres; ox++) {
227 greg 2.49 /*
228     * Super-sample when we detect a difference from before
229 greg 2.50 * or after in this row, above or below.
230 greg 2.49 */
231 greg 2.1 for (oy = 0; oy < sqres; oy++) {
232 greg 2.50 if (ox < sqres-1) { /* keeping one row ahead... */
233     SDsquare2disk(iovec+3, (ox+1.5)*sqfact, (oy+.5)*sqfact);
234 greg 2.47 iovec[5] = output_orient *
235 greg 2.11 sqrt(1. - iovec[3]*iovec[3] - iovec[4]*iovec[4]);
236 greg 2.47 }
237 greg 2.34 if (funame == NULL) {
238 greg 2.50 SDValue sdv = sdv_next[oy];
239     bsdf = sdv.cieY;
240     if (ox < sqres-1)
241     eval_rbfcol(&sdv_next[oy], rbf, iovec+3);
242     if (abs_diff(bsdf, sdv_next[oy].cieY) > ssamp_thresh ||
243     (ox && abs_diff(bsdf, val_last[oy]) > ssamp_thresh) ||
244     (oy && abs_diff(bsdf, val_last[oy-1]) > ssamp_thresh) ||
245     (oy < sqres-1 &&
246     abs_diff(bsdf, sdv_next[oy+1].cieY) > ssamp_thresh)) {
247 greg 2.42 int ssi;
248     double ssa[2], sum = 0, usum = 0, vsum = 0;
249     /* super-sample voxel */
250 greg 2.52 for (ssi = nssamp; ssi--; ) {
251     SDmultiSamp(ssa, 2, (ssi+frandom()) /
252     (double)nssamp);
253 greg 2.42 SDsquare2disk(iovec+3, (ox+ssa[0])*sqfact,
254     (oy+ssa[1])*sqfact);
255     iovec[5] = output_orient *
256     sqrt(1. - iovec[3]*iovec[3] - iovec[4]*iovec[4]);
257     eval_rbfcol(&sdv, rbf, iovec+3);
258     sum += sdv.cieY;
259     if (rbf_colorimetry == RBCtristimulus) {
260     sdv.cieY /=
261     -2.*sdv.spec.cx + 12.*sdv.spec.cy + 3.;
262     usum += 4.*sdv.spec.cx * sdv.cieY;
263     vsum += 9.*sdv.spec.cy * sdv.cieY;
264     }
265     }
266 greg 2.52 bsdf = sum / (double)nssamp;
267 greg 2.42 if (rbf_colorimetry == RBCtristimulus) {
268     uv[0] = usum / (sum+FTINY);
269     uv[1] = vsum / (sum+FTINY);
270     }
271     } else
272 greg 2.34 if (rbf_colorimetry == RBCtristimulus) {
273     uv[0] = uv[1] = 1. /
274     (-2.*sdv.spec.cx + 12.*sdv.spec.cy + 3.);
275     uv[0] *= 4.*sdv.spec.cx;
276     uv[1] *= 9.*sdv.spec.cy;
277     }
278     } else {
279 greg 2.50 bsdf = val_next[oy];
280     if (ox < sqres-1) {
281 greg 2.47 if (assignD) {
282     varset("Dx", '=', -iovec[3]);
283     varset("Dy", '=', -iovec[4]);
284     varset("Dz", '=', -iovec[5]);
285     ++eclock;
286     }
287 greg 2.50 val_next[oy] = funvalue(funame, 6, iovec);
288 greg 2.17 }
289 greg 2.50 if (abs_diff(bsdf, val_next[oy]) > ssamp_thresh ||
290     (ox && abs_diff(bsdf, val_last[oy]) > ssamp_thresh) ||
291     (oy && abs_diff(bsdf, val_last[oy-1]) > ssamp_thresh) ||
292     (oy < sqres-1 &&
293     abs_diff(bsdf, val_next[oy+1]) > ssamp_thresh)) {
294 greg 2.28 int ssi;
295 greg 2.45 double ssa[4], ssvec[6], sum = 0;
296 greg 2.28 /* super-sample voxel */
297 greg 2.52 for (ssi = nssamp; ssi--; ) {
298     SDmultiSamp(ssa, 4, (ssi+frandom()) /
299     (double)nssamp);
300 greg 2.42 ssvec[0] = 2.*sqfact*(ix+ssa[0]) - 1.;
301 greg 2.45 ssvec[1] = zipsgn*sqfact*ssa[1];
302     ssvec[2] = 1. - ssvec[0]*ssvec[0]
303     - ssvec[1]*ssvec[1];
304     if (ssvec[2] < .0) {
305     ssvec[1] = 0;
306     ssvec[2] = 1. - ssvec[0]*ssvec[0];
307     }
308     ssvec[2] = input_orient * sqrt(ssvec[2]);
309     SDsquare2disk(ssvec+3, (ox+ssa[2])*sqfact,
310     (oy+ssa[3])*sqfact);
311 greg 2.15 ssvec[5] = output_orient *
312     sqrt(1. - ssvec[3]*ssvec[3] -
313     ssvec[4]*ssvec[4]);
314 greg 2.17 if (assignD) {
315 greg 2.28 varset("Dx", '=', -ssvec[3]);
316     varset("Dy", '=', -ssvec[4]);
317     varset("Dz", '=', -ssvec[5]);
318 greg 2.17 ++eclock;
319     }
320 greg 2.16 sum += funvalue(funame, 6, ssvec);
321 greg 2.15 }
322 greg 2.52 bsdf = sum / (double)nssamp;
323 greg 2.15 }
324     }
325 greg 2.1 if (pctcull >= 0)
326 greg 2.37 putbinary(&bsdf, sizeof(bsdf), 1, ofp);
327 greg 2.1 else
328 greg 2.34 fprintf(ofp, "\t%.3e\n", bsdf);
329    
330     if (rbf_colorimetry == RBCtristimulus) {
331     if (pctcull >= 0) {
332 greg 2.37 putbinary(&uv[0], sizeof(*uv), 1, uvfp[0]);
333     putbinary(&uv[1], sizeof(*uv), 1, uvfp[1]);
334 greg 2.34 } else {
335     fprintf(uvfp[0], "\t%.3e\n", uv[0]);
336     fprintf(uvfp[1], "\t%.3e\n", uv[1]);
337     }
338     }
339 greg 2.50 if (val_last != NULL)
340     val_last[oy] = bsdf;
341 greg 2.1 }
342 greg 2.15 }
343 greg 2.3 if (rbf != NULL)
344     free(rbf);
345 greg 2.42 prog_show((ix+1.)*(2.*sqfact));
346 greg 2.1 }
347 greg 2.35 prog_done();
348 greg 2.50 if (val_last != NULL) {
349     free(val_last);
350     if (val_next != NULL) free(val_next);
351     if (sdv_next != NULL) free(sdv_next);
352 greg 2.47 }
353 greg 2.1 if (pctcull >= 0) { /* finish output */
354     if (pclose(ofp)) {
355 greg 2.34 fprintf(stderr, "%s: error running rttree_reduce on Y\n",
356     progname);
357     exit(1);
358     }
359     if (rbf_colorimetry == RBCtristimulus &&
360     (pclose(uvfp[0]) || pclose(uvfp[1]))) {
361     fprintf(stderr, "%s: error running rttree_reduce on uv\n",
362     progname);
363 greg 2.1 exit(1);
364     }
365     } else {
366     for (ix = sqres*sqres*sqres/2; ix--; )
367 greg 2.34 fputs("\t0\n", ofp);
368     fputs("}\n", ofp);
369     if (fclose(ofp)) {
370     fprintf(stderr, "%s: error writing Y file\n",
371     progname);
372     exit(1);
373     }
374     if (rbf_colorimetry == RBCtristimulus) {
375     for (ix = sqres*sqres*sqres/2; ix--; ) {
376     fputs("\t0\n", uvfp[0]);
377     fputs("\t0\n", uvfp[1]);
378     }
379     fputs("}\n", uvfp[0]);
380     fputs("}\n", uvfp[1]);
381     if (fclose(uvfp[0]) || fclose(uvfp[1])) {
382     fprintf(stderr, "%s: error writing uv file(s)\n",
383     progname);
384     exit(1);
385     }
386     }
387 greg 2.1 }
388     }
389    
390     /* Interpolate and output anisotropic BSDF data */
391     static void
392 greg 2.11 eval_anisotropic(char *funame)
393 greg 2.1 {
394     const int sqres = 1<<samp_order;
395 greg 2.42 const double sqfact = 1./(double)sqres;
396 greg 2.50 float *val_last = NULL;
397     float *val_next = NULL;
398     SDValue *sdv_next = NULL;
399 greg 2.34 FILE *ofp, *uvfp[2];
400 greg 2.17 int assignD = 0;
401 greg 2.1 char cmd[128];
402     int ix, iy, ox, oy;
403 greg 2.11 double iovec[6];
404 greg 2.34 float bsdf, uv[2];
405 greg 2.15
406 greg 2.11 if (pctcull >= 0) {
407 greg 2.35 const char *avgopt = (input_orient>0 ^ output_orient>0)
408 greg 2.41 ? "" : recip;
409 greg 2.34 sprintf(cmd, "rttree_reduce%s -h -ff -r 4 -t %f -g %d > %s",
410 greg 2.35 avgopt, pctcull, samp_order,
411     create_component_file(0));
412 greg 2.1 ofp = popen(cmd, "w");
413     if (ofp == NULL) {
414     fprintf(stderr, "%s: cannot create pipe to rttree_reduce\n",
415     progname);
416     exit(1);
417     }
418 greg 2.25 SET_FILE_BINARY(ofp);
419     #ifdef getc_unlocked /* avoid lock/unlock overhead */
420     flockfile(ofp);
421     #endif
422 greg 2.34 if (rbf_colorimetry == RBCtristimulus) {
423     double uvcull = 100. - (100.-pctcull)*.25;
424     sprintf(cmd, "rttree_reduce%s -h -ff -r 4 -t %f -g %d > %s",
425 greg 2.35 avgopt, uvcull, samp_order,
426     create_component_file(1));
427 greg 2.34 uvfp[0] = popen(cmd, "w");
428     sprintf(cmd, "rttree_reduce%s -h -ff -r 4 -t %f -g %d > %s",
429 greg 2.35 avgopt, uvcull, samp_order,
430     create_component_file(2));
431 greg 2.34 uvfp[1] = popen(cmd, "w");
432     if ((uvfp[0] == NULL) | (uvfp[1] == NULL)) {
433     fprintf(stderr, "%s: cannot open pipes to uv output\n",
434     progname);
435     exit(1);
436     }
437     SET_FILE_BINARY(uvfp[0]); SET_FILE_BINARY(uvfp[1]);
438     #ifdef getc_unlocked
439     flockfile(uvfp[0]); flockfile(uvfp[1]);
440     #endif
441     }
442     } else {
443     ofp = fopen(create_component_file(0), "w");
444     if (ofp == NULL) {
445     fprintf(stderr, "%s: cannot create Y output file\n",
446     progname);
447     exit(1);
448     }
449     fputs("{\n", ofp);
450     if (rbf_colorimetry == RBCtristimulus) {
451     uvfp[0] = fopen(create_component_file(1), "w");
452     uvfp[1] = fopen(create_component_file(2), "w");
453     if ((uvfp[0] == NULL) | (uvfp[1] == NULL)) {
454     fprintf(stderr, "%s: cannot create uv output file(s)\n",
455     progname);
456     exit(1);
457     }
458     fputs("{\n", uvfp[0]);
459     fputs("{\n", uvfp[1]);
460     }
461     }
462 greg 2.35 if (funame != NULL) /* need to assign Dx, Dy, Dz? */
463 greg 2.17 assignD = (fundefined(funame) < 6);
464 greg 2.50 val_last = (float *)calloc(sqres, sizeof(float));
465     if (funame == NULL)
466     sdv_next = (SDValue *)malloc(sizeof(SDValue)*sqres);
467     else
468     val_next = (float *)malloc(sizeof(float)*sqres);
469 greg 2.1 /* run through directions */
470     for (ix = 0; ix < sqres; ix++)
471     for (iy = 0; iy < sqres; iy++) {
472 greg 2.11 RBFNODE *rbf = NULL; /* Klems reversal */
473 greg 2.42 SDsquare2disk(iovec, 1.-(ix+.5)*sqfact, 1.-(iy+.5)*sqfact);
474 greg 2.11 iovec[2] = input_orient *
475     sqrt(1. - iovec[0]*iovec[0] - iovec[1]*iovec[1]);
476     if (funame == NULL)
477 greg 2.18 rbf = advect_rbf(iovec, lobe_lim);
478 greg 2.50 /* presample first row */
479     for (oy = 0; oy < sqres; oy++) {
480     SDsquare2disk(iovec+3, .5*sqfact, (oy+.5)*sqfact);
481 greg 2.47 iovec[5] = output_orient *
482     sqrt(1. - iovec[3]*iovec[3] - iovec[4]*iovec[4]);
483     if (funame == NULL) {
484 greg 2.50 eval_rbfcol(&sdv_next[oy], rbf, iovec+3);
485 greg 2.47 } else {
486     if (assignD) {
487     varset("Dx", '=', -iovec[3]);
488     varset("Dy", '=', -iovec[4]);
489     varset("Dz", '=', -iovec[5]);
490     ++eclock;
491     }
492 greg 2.50 val_next[oy] = funvalue(funame, 6, iovec);
493 greg 2.47 }
494 greg 2.50 }
495     for (ox = 0; ox < sqres; ox++) {
496 greg 2.49 /*
497     * Super-sample when we detect a difference from before
498 greg 2.50 * or after in this row, above or below.
499 greg 2.49 */
500 greg 2.1 for (oy = 0; oy < sqres; oy++) {
501 greg 2.50 if (ox < sqres-1) { /* keeping one row ahead... */
502     SDsquare2disk(iovec+3, (ox+1.5)*sqfact, (oy+.5)*sqfact);
503 greg 2.47 iovec[5] = output_orient *
504 greg 2.11 sqrt(1. - iovec[3]*iovec[3] - iovec[4]*iovec[4]);
505 greg 2.47 }
506 greg 2.34 if (funame == NULL) {
507 greg 2.50 SDValue sdv = sdv_next[oy];
508     bsdf = sdv.cieY;
509     if (ox < sqres-1)
510     eval_rbfcol(&sdv_next[oy], rbf, iovec+3);
511     if (abs_diff(bsdf, sdv_next[oy].cieY) > ssamp_thresh ||
512     (ox && abs_diff(bsdf, val_last[oy]) > ssamp_thresh) ||
513     (oy && abs_diff(bsdf, val_last[oy-1]) > ssamp_thresh) ||
514     (oy < sqres-1 &&
515     abs_diff(bsdf, sdv_next[oy+1].cieY) > ssamp_thresh)) {
516 greg 2.42 int ssi;
517     double ssa[2], sum = 0, usum = 0, vsum = 0;
518     /* super-sample voxel */
519 greg 2.52 for (ssi = nssamp; ssi--; ) {
520     SDmultiSamp(ssa, 2, (ssi+frandom()) /
521     (double)nssamp);
522 greg 2.42 SDsquare2disk(iovec+3, (ox+ssa[0])*sqfact,
523     (oy+ssa[1])*sqfact);
524     iovec[5] = output_orient *
525     sqrt(1. - iovec[3]*iovec[3] - iovec[4]*iovec[4]);
526     eval_rbfcol(&sdv, rbf, iovec+3);
527     sum += sdv.cieY;
528     if (rbf_colorimetry == RBCtristimulus) {
529     sdv.cieY /=
530     -2.*sdv.spec.cx + 12.*sdv.spec.cy + 3.;
531     usum += 4.*sdv.spec.cx * sdv.cieY;
532     vsum += 9.*sdv.spec.cy * sdv.cieY;
533     }
534     }
535 greg 2.52 bsdf = sum / (double)nssamp;
536 greg 2.42 if (rbf_colorimetry == RBCtristimulus) {
537     uv[0] = usum / (sum+FTINY);
538     uv[1] = vsum / (sum+FTINY);
539     }
540     } else
541 greg 2.34 if (rbf_colorimetry == RBCtristimulus) {
542     uv[0] = uv[1] = 1. /
543     (-2.*sdv.spec.cx + 12.*sdv.spec.cy + 3.);
544     uv[0] *= 4.*sdv.spec.cx;
545     uv[1] *= 9.*sdv.spec.cy;
546     }
547     } else {
548 greg 2.50 bsdf = val_next[oy];
549     if (ox < sqres-1) {
550 greg 2.47 if (assignD) {
551     varset("Dx", '=', -iovec[3]);
552     varset("Dy", '=', -iovec[4]);
553     varset("Dz", '=', -iovec[5]);
554     ++eclock;
555     }
556 greg 2.50 val_next[oy] = funvalue(funame, 6, iovec);
557 greg 2.17 }
558 greg 2.50 if (abs_diff(bsdf, val_next[oy]) > ssamp_thresh ||
559     (ox && abs_diff(bsdf, val_last[oy]) > ssamp_thresh) ||
560     (oy && abs_diff(bsdf, val_last[oy-1]) > ssamp_thresh) ||
561     (oy < sqres-1 &&
562     abs_diff(bsdf, val_next[oy+1]) > ssamp_thresh)) {
563 greg 2.28 int ssi;
564     double ssa[4], ssvec[6], sum = 0;
565     /* super-sample voxel */
566 greg 2.52 for (ssi = nssamp; ssi--; ) {
567     SDmultiSamp(ssa, 4, (ssi+frandom()) /
568     (double)nssamp);
569 greg 2.42 SDsquare2disk(ssvec, 1.-(ix+ssa[0])*sqfact,
570     1.-(iy+ssa[1])*sqfact);
571 greg 2.28 ssvec[2] = input_orient *
572 greg 2.15 sqrt(1. - ssvec[0]*ssvec[0] -
573     ssvec[1]*ssvec[1]);
574 greg 2.42 SDsquare2disk(ssvec+3, (ox+ssa[2])*sqfact,
575     (oy+ssa[3])*sqfact);
576 greg 2.15 ssvec[5] = output_orient *
577     sqrt(1. - ssvec[3]*ssvec[3] -
578     ssvec[4]*ssvec[4]);
579 greg 2.17 if (assignD) {
580 greg 2.28 varset("Dx", '=', -ssvec[3]);
581     varset("Dy", '=', -ssvec[4]);
582     varset("Dz", '=', -ssvec[5]);
583 greg 2.17 ++eclock;
584     }
585 greg 2.16 sum += funvalue(funame, 6, ssvec);
586 greg 2.15 }
587 greg 2.52 bsdf = sum / (double)nssamp;
588 greg 2.15 }
589     }
590 greg 2.1 if (pctcull >= 0)
591 greg 2.37 putbinary(&bsdf, sizeof(bsdf), 1, ofp);
592 greg 2.1 else
593 greg 2.34 fprintf(ofp, "\t%.3e\n", bsdf);
594    
595     if (rbf_colorimetry == RBCtristimulus) {
596     if (pctcull >= 0) {
597 greg 2.37 putbinary(&uv[0], sizeof(*uv), 1, uvfp[0]);
598     putbinary(&uv[1], sizeof(*uv), 1, uvfp[1]);
599 greg 2.34 } else {
600     fprintf(uvfp[0], "\t%.3e\n", uv[0]);
601     fprintf(uvfp[1], "\t%.3e\n", uv[1]);
602     }
603     }
604 greg 2.50 if (val_last != NULL)
605     val_last[oy] = bsdf;
606 greg 2.1 }
607 greg 2.51 }
608 greg 2.3 if (rbf != NULL)
609     free(rbf);
610 greg 2.26 prog_show((ix*sqres+iy+1.)/(sqres*sqres));
611 greg 2.1 }
612 greg 2.35 prog_done();
613 greg 2.50 if (val_last != NULL) {
614     free(val_last);
615     if (val_next != NULL) free(val_next);
616     if (sdv_next != NULL) free(sdv_next);
617 greg 2.47 }
618 greg 2.1 if (pctcull >= 0) { /* finish output */
619     if (pclose(ofp)) {
620 greg 2.34 fprintf(stderr, "%s: error running rttree_reduce on Y\n",
621     progname);
622     exit(1);
623     }
624     if (rbf_colorimetry == RBCtristimulus &&
625     (pclose(uvfp[0]) || pclose(uvfp[1]))) {
626     fprintf(stderr, "%s: error running rttree_reduce on uv\n",
627     progname);
628 greg 2.1 exit(1);
629     }
630 greg 2.34 } else {
631     fputs("}\n", ofp);
632     if (fclose(ofp)) {
633     fprintf(stderr, "%s: error writing Y file\n",
634     progname);
635     exit(1);
636     }
637     if (rbf_colorimetry == RBCtristimulus) {
638     fputs("}\n", uvfp[0]);
639     fputs("}\n", uvfp[1]);
640     if (fclose(uvfp[0]) || fclose(uvfp[1])) {
641     fprintf(stderr, "%s: error writing uv file(s)\n",
642     progname);
643     exit(1);
644     }
645     }
646     }
647 greg 2.5 }
648    
649 schorsch 2.36 #if defined(_WIN32) || defined(_WIN64)
650 greg 2.34 /* Execute wrapBSDF command (may never return) */
651     static int
652     wrap_up(void)
653     {
654 greg 2.55 char cmd[32700];
655 greg 2.34
656     if (bsdf_manuf[0]) {
657     add_wbsdf("-f", 1);
658     strcpy(cmd, "m=");
659     strcpy(cmd+2, bsdf_manuf);
660     add_wbsdf(cmd, 0);
661     }
662     if (bsdf_name[0]) {
663     add_wbsdf("-f", 1);
664     strcpy(cmd, "n=");
665     strcpy(cmd+2, bsdf_name);
666     add_wbsdf(cmd, 0);
667     }
668     if (!convert_commandline(cmd, sizeof(cmd), wrapBSDF)) {
669     fputs(progname, stderr);
670     fputs(": command line too long in wrap_up()\n", stderr);
671     return(1);
672     }
673     return(system(cmd));
674     }
675     #else
676     /* Execute wrapBSDF command (may never return) */
677     static int
678     wrap_up(void)
679     {
680     char buf[256];
681     char *compath = getpath((char *)wrapBSDF[0], getenv("PATH"), X_OK);
682    
683     if (compath == NULL) {
684     fprintf(stderr, "%s: cannot locate %s\n", progname, wrapBSDF[0]);
685     return(1);
686     }
687     if (bsdf_manuf[0]) {
688     add_wbsdf("-f", 1);
689     strcpy(buf, "m=");
690     strcpy(buf+2, bsdf_manuf);
691     add_wbsdf(buf, 0);
692     }
693     if (bsdf_name[0]) {
694     add_wbsdf("-f", 1);
695     strcpy(buf, "n=");
696     strcpy(buf+2, bsdf_name);
697     add_wbsdf(buf, 0);
698     }
699     execv(compath, wrapBSDF); /* successful call never returns */
700     perror(compath);
701     return(1);
702     }
703     #endif
704    
705 greg 2.53 #define HEAD_BUFLEN 8192
706     static char head_buf[HEAD_BUFLEN];
707     static int cur_headlen = 0;
708    
709     /* Record header line as comment associated with this SIR input */
710     static int
711     record2header(char *s)
712     {
713     int len = strlen(s);
714    
715     if (cur_headlen+len >= HEAD_BUFLEN-6)
716     return(0);
717     /* includes EOL */
718     strcpy(head_buf+cur_headlen, s);
719     cur_headlen += len;
720    
721 greg 2.54 #if defined(_WIN32) || defined(_WIN64)
722     if (head_buf[cur_headlen-1] == '\n')
723     head_buf[cur_headlen-1] = '\t';
724     #endif
725 greg 2.53 return(1);
726     }
727    
728     /* Finish off header for this file */
729     static void
730     done_header(void)
731     {
732     while (cur_headlen > 0 && isspace(head_buf[cur_headlen-1]))
733     --cur_headlen;
734     head_buf[cur_headlen] = '\0';
735     if (!cur_headlen)
736     return;
737     add_wbsdf("-C", 1);
738     add_wbsdf(head_buf, 0);
739     head_buf[cur_headlen=0] = '\0';
740     }
741    
742 greg 2.1 /* Read in BSDF and interpolate as tensor tree representation */
743     int
744     main(int argc, char *argv[])
745     {
746 greg 2.34 static char tfmt[2][4] = {"t4", "t3"};
747 greg 2.11 int dofwd = 0, dobwd = 1;
748 greg 2.34 char buf[2048];
749 greg 2.11 int i, na;
750 greg 2.1
751 greg 2.5 progname = argv[0];
752 greg 2.11 esupport |= E_VARIABLE|E_FUNCTION|E_RCONST;
753     esupport &= ~(E_INCHAN|E_OUTCHAN);
754     scompile("PI:3.14159265358979323846", NULL, 0);
755     biggerlib();
756 greg 2.41 for (i = 1; i < argc && (argv[i][0] == '-') | (argv[i][0] == '+'); i++)
757 greg 2.11 switch (argv[i][1]) { /* get options */
758     case 'e':
759     scompile(argv[++i], NULL, 0);
760 greg 2.46 if (single_plane_incident < 0)
761     single_plane_incident = 0;
762 greg 2.11 break;
763     case 'f':
764 greg 2.34 if (!argv[i][2]) {
765     if (strchr(argv[++i], '=') != NULL) {
766     add_wbsdf("-f", 1);
767     add_wbsdf(argv[i], 1);
768 greg 2.43 } else {
769     char *fpath = getpath(argv[i],
770     getrlibpath(), 0);
771     if (fpath == NULL) {
772     fprintf(stderr,
773     "%s: cannot find file '%s'\n",
774     argv[0], argv[i]);
775     return(1);
776     }
777     fcompile(fpath);
778 greg 2.45 if (single_plane_incident < 0)
779     single_plane_incident = 0;
780 greg 2.43 }
781 greg 2.34 } else
782 greg 2.11 dofwd = (argv[i][0] == '+');
783     break;
784 greg 2.41 case 'a':
785     recip = (argv[i][0] == '+') ? " -a" : "";
786     break;
787 greg 2.11 case 'b':
788     dobwd = (argv[i][0] == '+');
789     break;
790 greg 2.52 case 'n':
791     nssamp = atoi(argv[++i]);
792     if (nssamp <= 0)
793     goto userr;
794     break;
795     case 's':
796     ssamp_thresh = atof(argv[++i]);
797     if (ssamp_thresh <= FTINY)
798     goto userr;
799     break;
800 greg 2.1 case 't':
801 greg 2.11 switch (argv[i][2]) {
802     case '3':
803 greg 2.13 single_plane_incident = 1;
804 greg 2.11 break;
805     case '4':
806 greg 2.13 single_plane_incident = 0;
807 greg 2.11 break;
808     case '\0':
809     pctcull = atof(argv[++i]);
810     break;
811     default:
812     goto userr;
813     }
814 greg 2.1 break;
815     case 'g':
816 greg 2.5 samp_order = atoi(argv[++i]);
817 greg 2.1 break;
818 greg 2.18 case 'l':
819     lobe_lim = atoi(argv[++i]);
820     break;
821 greg 2.26 case 'p':
822     do_prog = atoi(argv[i]+2);
823     break;
824 greg 2.34 case 'W':
825     add_wbsdf(argv[i], 1);
826     break;
827     case 'u':
828     case 'C':
829     add_wbsdf(argv[i], 1);
830     add_wbsdf(argv[++i], 1);
831     break;
832 greg 2.1 default:
833     goto userr;
834     }
835 greg 2.34 strcpy(buf, "File produced by: ");
836     if (convert_commandline(buf+18, sizeof(buf)-18, argv) != NULL) {
837     add_wbsdf("-C", 1); add_wbsdf(buf, 0);
838     }
839 greg 2.11 if (single_plane_incident >= 0) { /* function-based BSDF? */
840     void (*evf)(char *s) = single_plane_incident ?
841 greg 2.34 eval_isotropic : eval_anisotropic;
842 greg 2.17 if (i != argc-1 || fundefined(argv[i]) < 3) {
843 greg 2.12 fprintf(stderr,
844     "%s: need single function with 6 arguments: bsdf(ix,iy,iz,ox,oy,oz)\n",
845     progname);
846 greg 2.24 fprintf(stderr, "\tor 3 arguments using Dx,Dy,Dz: bsdf(ix,iy,iz)\n");
847 greg 2.11 goto userr;
848 greg 2.12 }
849 greg 2.17 ++eclock;
850 greg 2.34 add_wbsdf("-a", 1);
851     add_wbsdf(tfmt[single_plane_incident], 1);
852 greg 2.11 if (dofwd) {
853     input_orient = -1;
854     output_orient = -1;
855 greg 2.26 prog_start("Evaluating outside reflectance");
856     (*evf)(argv[i]);
857 greg 2.11 output_orient = 1;
858 greg 2.26 prog_start("Evaluating outside->inside transmission");
859     (*evf)(argv[i]);
860 greg 2.11 }
861     if (dobwd) {
862     input_orient = 1;
863     output_orient = 1;
864 greg 2.26 prog_start("Evaluating inside reflectance");
865     (*evf)(argv[i]);
866 greg 2.11 output_orient = -1;
867 greg 2.26 prog_start("Evaluating inside->outside transmission");
868     (*evf)(argv[i]);
869 greg 2.11 }
870 greg 2.34 return(wrap_up());
871 greg 2.11 }
872     if (i < argc) { /* open input files if given */
873 greg 2.12 int nbsdf = 0;
874 greg 2.11 for ( ; i < argc; i++) { /* interpolate each component */
875 greg 2.26 char pbuf[256];
876 greg 2.11 FILE *fpin = fopen(argv[i], "rb");
877     if (fpin == NULL) {
878     fprintf(stderr, "%s: cannot open BSDF interpolant '%s'\n",
879     progname, argv[i]);
880     return(1);
881     }
882 greg 2.53 sprintf(pbuf, "%s:\n", argv[i]);
883     record2header(pbuf);
884     sir_headshare = &record2header;
885 greg 2.11 if (!load_bsdf_rep(fpin))
886     return(1);
887     fclose(fpin);
888 greg 2.53 done_header();
889 greg 2.26 sprintf(pbuf, "Interpolating component '%s'", argv[i]);
890     prog_start(pbuf);
891 greg 2.34 if (!nbsdf++) {
892     add_wbsdf("-a", 1);
893     add_wbsdf(tfmt[single_plane_incident], 1);
894     }
895 greg 2.11 if (single_plane_incident)
896     eval_isotropic(NULL);
897     else
898     eval_anisotropic(NULL);
899     }
900 greg 2.34 return(wrap_up());
901 greg 2.11 }
902     SET_FILE_BINARY(stdin); /* load from stdin */
903     if (!load_bsdf_rep(stdin))
904 greg 2.1 return(1);
905 greg 2.26 prog_start("Interpolating from standard input");
906 greg 2.34 add_wbsdf("-a", 1);
907     add_wbsdf(tfmt[single_plane_incident], 1);
908 greg 2.1 if (single_plane_incident) /* resample dist. */
909 greg 2.11 eval_isotropic(NULL);
910 greg 2.1 else
911 greg 2.11 eval_anisotropic(NULL);
912 greg 2.34
913     return(wrap_up());
914 greg 2.1 userr:
915     fprintf(stderr,
916 greg 2.52 "Usage: %s [{+|-}a][-g Nlog2][-t pctcull][-n nss][-s thresh][-l maxlobes] [bsdf.sir ..] > bsdf.xml\n",
917 greg 2.11 progname);
918     fprintf(stderr,
919 greg 2.52 " or: %s -t{3|4} [{+|-}a][-g Nlog2][-t pctcull][-n nss][-s thresh][{+|-}for[ward]][{+|-}b[ackward]][-e expr][-f file] bsdf_func > bsdf.xml\n",
920 greg 2.1 progname);
921     return(1);
922     }