ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2ttree.c
Revision: 2.63
Committed: Tue Jun 3 21:31:51 2025 UTC (3 days, 20 hours ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.62: +3 -5 lines
Log Message: 
refactor: More consistent use of global char * progname and fixargv0()

File Contents

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