ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2ttree.c
Revision: 2.64
Committed: Sat Jun 7 05:09:45 2025 UTC (14 hours, 37 minutes ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 2.63: +1 -3 lines
Log Message: 
refactor: Put some declarations into "paths.h" and included in "platform.h"

File Contents

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