ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2ttree.c
Revision: 2.52
Committed: Mon May 18 20:08:57 2020 UTC (3 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R3
Changes since 2.51: +30 -35 lines
Log Message: 
Added -n and -s options to bsdf2ttree to control super-sampling

File Contents

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