ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2ttree.c
Revision: 2.61
Committed: Fri Feb 23 03:47:57 2024 UTC (3 months, 3 weeks ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 2.60: +2 -1 lines
Log Message:
perf: Added array index optimization to calcomp routines

File Contents

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