ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2ttree.c
Revision: 2.51
Committed: Fri May 15 03:02:32 2020 UTC (3 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.50: +2 -2 lines
Log Message: 
Fixed curly bracket in the wrong place (loop is WAAY too long!)

File Contents

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