ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/Development/ray/src/cv/bsdf2ttree.c
Revision: 2.49
Committed: Fri May 8 22:34:21 2020 UTC (5 years, 5 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.48: +17 -3 lines
Log Message: 
Forgot to clear the flags from previous row

File Contents

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