ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2rad.c
(Generate patch)

Comparing ray/src/cv/bsdf2rad.c (file contents):
Revision 2.9 by greg, Mon Mar 24 06:07:46 2014 UTC vs.
Revision 2.41 by greg, Tue Jun 3 21:31:51 2025 UTC

# Line 5 | Line 5 | static const char RCSid[] = "$Id$";
5   *  Plot 3-D BSDF output based on scattering interpolant or XML representation
6   */
7  
8 #define _USE_MATH_DEFINES
9 #include <stdio.h>
10 #include <string.h>
8   #include <stdlib.h>
9 < #include <math.h>
10 < #include "rtprocess.h"
9 > #include "rtio.h"
10 > #include "paths.h"
11 > #include "rtmath.h"
12   #include "bsdfrep.h"
13  
14 < const float     colarr[6][3] = {
15 <                .7, 1., .7,
16 <                1., .7, .7,
17 <                .7, .7, 1.,
18 <                1., .5, 1.,
19 <                1., 1., .5,
20 <                .5, 1., 1.
23 <        };
14 > #ifndef NINCIDENT
15 > #define NINCIDENT       37              /* number of samples/hemisphere */
16 > #endif
17 > #ifndef GRIDSTEP
18 > #define GRIDSTEP        2               /* our grid step size */
19 > #endif
20 > #define SAMPRES         (GRIDRES/GRIDSTEP)
21  
22 < char    *progname;
22 > int     front_comp = 0;                 /* front component flags (SDsamp*) */
23 > int     back_comp = 0;                  /* back component flags */
24 > double  overall_min = 1./PI;            /* overall minimum BSDF value */
25 > double  min_log10;                      /* smallest log10 value for plotting */
26 > double  overall_max = .0;               /* overall maximum BSDF value */
27  
28 < /* Produce a Radiance model plotting the indicated incident direction(s) */
28 > char    ourTempDir[TEMPLEN+1] = "";     /* our temporary directory */
29 >
30 > const char      frpref[] = "rf";
31 > const char      ftpref[] = "tf";
32 > const char      brpref[] = "rb";
33 > const char      btpref[] = "tb";
34 > const char      dsuffix[] = ".txt";
35 >
36 > const char      sph_fmat[] = "fBSDFmat";
37 > const char      sph_bmat[] = "bBSDFmat";
38 > const double    sph_rad = 10.;
39 > const double    sph_xoffset = 15.;
40 >
41 > #define bsdf_rad        (sph_rad*.25)
42 > #define arrow_rad       (bsdf_rad*.015)
43 >
44 > #define set_minlog()    overall_min = (overall_min < 1e-5) ? 1e-5 : overall_min; \
45 >                                min_log10 = log10(overall_min) - .1
46 >
47 > /* Get Fibonacci sphere vector (0 to NINCIDENT-1) */
48 > static RREAL *
49 > get_ivector(FVECT iv, int i)
50 > {
51 >        const double    phistep = PI*(3. - 2.236067978);
52 >        double          r;
53 >
54 >        iv[2] = 1. - (i+.5)*(1./NINCIDENT);
55 >        r = sqrt(1. - iv[2]*iv[2]);
56 >        iv[0] = r * cos((i+1.)*phistep);
57 >        iv[1] = r * sin((i+1.)*phistep);
58 >
59 >        return(iv);
60 > }
61 >
62 > /* Convert incident vector into sphere position */
63 > static RREAL *
64 > cvt_sposition(FVECT sp, const FVECT iv, int inc_side)
65 > {
66 >        sp[0] = -iv[0]*sph_rad + inc_side*sph_xoffset;
67 >        sp[1] = -iv[1]*sph_rad;
68 >        sp[2] = iv[2]*sph_rad;
69 >
70 >        return(sp);
71 > }
72 >
73 > /* Get temporary file name */
74 > static char *
75 > tfile_name(const char *prefix, const char *suffix, int i)
76 > {
77 >        static char     buf[256];
78 >
79 >        if (!ourTempDir[0]) {           /* create temporary directory */
80 >                mktemp(strcpy(ourTempDir,TEMPLATE));
81 >                if (mkdir(ourTempDir, 0777) < 0) {
82 >                        perror("mkdir");
83 >                        exit(1);
84 >                }
85 >        }
86 >        if (!prefix) prefix = "T";
87 >        if (!suffix) suffix = "";
88 >        sprintf(buf, "%s/%s%03d%s", ourTempDir, prefix, i, suffix);
89 >        return(buf);
90 > }
91 >
92 > /* Remove temporary directory & contents */
93 > static void
94 > cleanup_tmp(void)
95 > {
96 >        char    buf[128];
97 >
98 >        if (!ourTempDir[0])
99 >                return;
100 > #if defined(_WIN32) || defined(_WIN64)
101 >        sprintf(buf, "RMDIR %s /S /Q", ourTempDir);
102 > #else
103 >        sprintf(buf, "rm -rf %s", ourTempDir);
104 > #endif
105 >        system(buf);
106 > }
107 >
108 > /* Run the specified command, returning 1 if OK */
109 > static int
110 > run_cmd(const char *cmd)
111 > {
112 >        fflush(stdout);
113 >        if (system(cmd)) {
114 >                fprintf(stderr, "%s: error running: %s\n", progname, cmd);
115 >                return(0);
116 >        }
117 >        return(1);
118 > }
119 >
120 > /* Plot surface points for the given BSDF incident angle */
121 > static int
122 > plotBSDF(const char *fname, const FVECT ivec, int dfl, const SDData *sd)
123 > {
124 >        FILE            *fp = fopen(fname, "w");
125 >        int             i, j;
126 >
127 >        if (fp == NULL) {
128 >                fprintf(stderr, "%s: cannot open '%s' for writing\n",
129 >                                progname, fname);
130 >                return(0);
131 >        }
132 >        if (ivec[2] > 0) {
133 >                input_orient = 1;
134 >                output_orient = dfl&SDsampR ? 1 : -1;
135 >        } else {
136 >                input_orient = -1;
137 >                output_orient = dfl&SDsampR ? -1 : 1;
138 >        }
139 >        for (i = SAMPRES; i--; )
140 >            for (j = 0; j < SAMPRES; j++) {
141 >                FVECT   ovec;
142 >                SDValue sval;
143 >                double  bsdf;
144 >                ovec_from_pos(ovec, i*GRIDSTEP, j*GRIDSTEP);
145 >                if (SDreportError(SDevalBSDF(&sval, ivec,
146 >                                                ovec, sd), stderr))
147 >                        return(0);
148 >                if (sval.cieY > overall_max)
149 >                        overall_max = sval.cieY;
150 >                bsdf = (sval.cieY < overall_min) ? overall_min : sval.cieY;
151 >                bsdf = log10(bsdf) - min_log10;
152 >                fprintf(fp, "%.5f %.5f %.5f\n",
153 >                                ovec[0]*bsdf, ovec[1]*bsdf, ovec[2]*bsdf);
154 >            }
155 >        if (fclose(fp) == EOF) {
156 >                fprintf(stderr, "%s: error writing data to '%s'\n",
157 >                                progname, fname);
158 >                return(0);
159 >        }
160 >        return(1);
161 > }
162 >
163 > /* Build BSDF values from loaded XML file */
164 > static int
165 > build_wBSDF(const SDData *sd)
166 > {
167 >        FVECT   ivec;
168 >        int     i;
169 >
170 >        if (front_comp & SDsampR)
171 >                for (i = 0; i < NINCIDENT; i++) {
172 >                        get_ivector(ivec, i);
173 >                        if (!plotBSDF(tfile_name(frpref, dsuffix, i),
174 >                                        ivec, SDsampR, sd))
175 >                                return(0);
176 >                }
177 >        if (front_comp & SDsampT)
178 >                for (i = 0; i < NINCIDENT; i++) {
179 >                        get_ivector(ivec, i);
180 >                        if (!plotBSDF(tfile_name(ftpref, dsuffix, i),
181 >                                        ivec, SDsampT, sd))
182 >                                return(0);
183 >                }
184 >        if (back_comp & SDsampR)
185 >                for (i = 0; i < NINCIDENT; i++) {
186 >                        get_ivector(ivec, i);
187 >                        ivec[0] = -ivec[0]; ivec[2] = -ivec[2];
188 >                        if (!plotBSDF(tfile_name(brpref, dsuffix, i),
189 >                                        ivec, SDsampR, sd))
190 >                                return(0);
191 >                }
192 >        if (back_comp & SDsampT)
193 >                for (i = 0; i < NINCIDENT; i++) {
194 >                        get_ivector(ivec, i);
195 >                        ivec[0] = -ivec[0]; ivec[2] = -ivec[2];
196 >                        if (!plotBSDF(tfile_name(btpref, dsuffix, i),
197 >                                        ivec, SDsampT, sd))
198 >                                return(0);
199 >                }
200 >        return(1);
201 > }
202 >
203 > /* Plot surface points using radial basis function */
204 > static int
205 > plotRBF(const char *fname, const RBFNODE *rbf)
206 > {
207 >        FILE            *fp = fopen(fname, "w");
208 >        int             i, j;
209 >
210 >        if (fp == NULL) {
211 >                fprintf(stderr, "%s: cannot open '%s' for writing\n",
212 >                                progname, fname);
213 >                return(0);
214 >        }
215 >        for (i = SAMPRES; i--; )
216 >            for (j = 0; j < SAMPRES; j++) {
217 >                FVECT   ovec;
218 >                double  bsdf;
219 >                ovec_from_pos(ovec, i*GRIDSTEP, j*GRIDSTEP);
220 >                bsdf = eval_rbfrep(rbf, ovec);
221 >                if (bsdf > overall_max)
222 >                        overall_max = bsdf;
223 >                else if (bsdf < overall_min)
224 >                        bsdf = overall_min;
225 >                bsdf = log10(bsdf) - min_log10;
226 >                fprintf(fp, "%.5f %.5f %.5f\n",
227 >                                ovec[0]*bsdf, ovec[1]*bsdf, ovec[2]*bsdf);
228 >            }
229 >        if (fclose(fp) == EOF) {
230 >                fprintf(stderr, "%s: error writing data to '%s'\n",
231 >                                progname, fname);
232 >                return(0);
233 >        }
234 >        return(1);
235 > }
236 >
237 > /* Build BSDF values from scattering interpolant representation */
238 > static int
239 > build_wRBF(void)
240 > {
241 >        const char      *pref;
242 >        int             i;
243 >
244 >        if (input_orient > 0) {
245 >                if (output_orient > 0)
246 >                        pref = frpref;
247 >                else
248 >                        pref = ftpref;
249 >        } else if (output_orient < 0)
250 >                pref = brpref;
251 >        else
252 >                pref = btpref;
253 >
254 >        for (i = 0; i < NINCIDENT; i++) {
255 >                FVECT   ivec;
256 >                RBFNODE *rbf;
257 >                get_ivector(ivec, i);
258 >                if (input_orient < 0) {
259 >                        ivec[0] = -ivec[0]; ivec[2] = -ivec[2];
260 >                }
261 >                rbf = advect_rbf(ivec, 15000);
262 >                if (!plotRBF(tfile_name(pref, dsuffix, i), rbf))
263 >                        return(0);
264 >                if (rbf) free(rbf);
265 >        }
266 >        return(1);                              /* next call frees */
267 > }
268 >
269 > /* Put out mirror arrow for the given incident vector */
270 > static void
271 > put_mirror_arrow(const FVECT origin, const FVECT nrm)
272 > {
273 >        const double    arrow_len = 1.2*bsdf_rad;
274 >        const double    tip_len = 0.2*bsdf_rad;
275 >        static int      cnt = 1;
276 >        FVECT           refl;
277 >        int             i;
278 >
279 >        refl[0] = 2.*nrm[2]*nrm[0];
280 >        refl[1] = 2.*nrm[2]*nrm[1];
281 >        refl[2] = 2.*nrm[2]*nrm[2] - 1.;
282 >
283 >        printf("\n# Mirror arrow #%d\n", cnt);
284 >        printf("\nshaft_mat cylinder inc_dir%d\n0\n0\n7", cnt);
285 >        printf("\n\t%f %f %f\n\t%f %f %f\n\t%f\n",
286 >                        origin[0], origin[1], origin[2]+arrow_len,
287 >                        origin[0], origin[1], origin[2],
288 >                        arrow_rad);
289 >        printf("\nshaft_mat cylinder mir_dir%d\n0\n0\n7", cnt);
290 >        printf("\n\t%f %f %f\n\t%f %f %f\n\t%f\n",
291 >                        origin[0], origin[1], origin[2],
292 >                        origin[0] + arrow_len*refl[0],
293 >                        origin[1] + arrow_len*refl[1],
294 >                        origin[2] + arrow_len*refl[2],
295 >                        arrow_rad);
296 >        printf("\ntip_mat cone mir_tip%d\n0\n0\n8", cnt);
297 >        printf("\n\t%f %f %f\n\t%f %f %f\n\t%f 0\n",
298 >                        origin[0] + (arrow_len-.5*tip_len)*refl[0],
299 >                        origin[1] + (arrow_len-.5*tip_len)*refl[1],
300 >                        origin[2] + (arrow_len-.5*tip_len)*refl[2],
301 >                        origin[0] + (arrow_len+.5*tip_len)*refl[0],
302 >                        origin[1] + (arrow_len+.5*tip_len)*refl[1],
303 >                        origin[2] + (arrow_len+.5*tip_len)*refl[2],
304 >                        2.*arrow_rad);
305 >        ++cnt;
306 > }
307 >
308 > /* Put out transmitted direction arrow for the given incident vector */
309 > static void
310 > put_trans_arrow(const FVECT origin)
311 > {
312 >        const double    arrow_len = 1.2*bsdf_rad;
313 >        const double    tip_len = 0.2*bsdf_rad;
314 >        static int      cnt = 1;
315 >        int             i;
316 >
317 >        printf("\n# Transmission arrow #%d\n", cnt);
318 >        printf("\nshaft_mat cylinder trans_dir%d\n0\n0\n7", cnt);
319 >        printf("\n\t%f %f %f\n\t%f %f %f\n\t%f\n",
320 >                        origin[0], origin[1], origin[2],
321 >                        origin[0], origin[1], origin[2]-arrow_len,
322 >                        arrow_rad);
323 >        printf("\ntip_mat cone trans_tip%d\n0\n0\n8", cnt);
324 >        printf("\n\t%f %f %f\n\t%f %f %f\n\t%f 0\n",
325 >                        origin[0], origin[1], origin[2]-arrow_len+.5*tip_len,
326 >                        origin[0], origin[1], origin[2]-arrow_len-.5*tip_len,
327 >                        2.*arrow_rad);
328 >        ++cnt;
329 > }
330 >
331 > /* Compute rotation (x,y,z) => (xp,yp,zp) */
332 > static int
333 > addrot(char *xf, const FVECT xp, const FVECT yp, const FVECT zp)
334 > {
335 >        int     n = 0;
336 >        double  theta;
337 >
338 >        if (yp[2]*yp[2] + zp[2]*zp[2] < 2.*FTINY*FTINY) {
339 >                /* Special case for X' along Z-axis */
340 >                theta = -atan2(yp[0], yp[1]);
341 >                sprintf(xf, " -ry %f -rz %f",
342 >                                xp[2] < 0.0 ? 90.0 : -90.0,
343 >                                theta*(180./PI));
344 >                return(4);
345 >        }
346 >        theta = atan2(yp[2], zp[2]);
347 >        if (!FABSEQ(theta,0.0)) {
348 >                sprintf(xf, " -rx %f", theta*(180./PI));
349 >                while (*xf) ++xf;
350 >                n += 2;
351 >        }
352 >        theta = Asin(-xp[2]);
353 >        if (!FABSEQ(theta,0.0)) {
354 >                sprintf(xf, " -ry %f", theta*(180./PI));
355 >                while (*xf) ++xf;
356 >                n += 2;
357 >        }
358 >        theta = atan2(xp[1], xp[0]);
359 >        if (!FABSEQ(theta,0.0)) {
360 >                sprintf(xf, " -rz %f", theta*(180./PI));
361 >                /* while (*xf) ++xf; */
362 >                n += 2;
363 >        }
364 >        return(n);
365 > }
366 >
367 > /* Put out BSDF surfaces */
368 > static int
369 > put_BSDFs(void)
370 > {
371 >        const double    scalef = bsdf_rad/(log10(overall_max) - min_log10);
372 >        FVECT           ivec, sorg, nrm, upv;
373 >        RREAL           vMtx[3][3];
374 >        char            *fname;
375 >        char            cmdbuf[256];
376 >        char            rotargs[64];
377 >        int             nrota;
378 >        int             i;
379 >
380 >        printf("\n# Gensurf output corresponding to %d incident directions\n",
381 >                        NINCIDENT);
382 >
383 >        printf("\nvoid glow tip_mat\n0\n0\n4 1 0 1 0\n");
384 >        printf("\nvoid mixfunc shaft_mat\n4 tip_mat void 0.25 .\n0\n0\n");
385 >
386 >        for (i = 0; i < NINCIDENT; i++) {
387 >                get_ivector(ivec, i);
388 >                nrm[0] = -ivec[0]; nrm[1] = -ivec[1]; nrm[2] = ivec[2];
389 >                upv[0] = nrm[0]*nrm[1]*(nrm[2] - 1.);
390 >                upv[1] = nrm[0]*nrm[0] + nrm[1]*nrm[1]*nrm[2];
391 >                upv[2] = -nrm[1]*(nrm[0]*nrm[0] + nrm[1]*nrm[1]);
392 >                if (SDcompXform(vMtx, nrm, upv) != SDEnone)
393 >                        continue;
394 >                nrota = addrot(rotargs, vMtx[0], vMtx[1], vMtx[2]);
395 >                if (front_comp) {
396 >                        cvt_sposition(sorg, ivec, 1);
397 >                        printf("\nvoid colorfunc scale_pat\n");
398 >                        printf("10 bsdf_red bsdf_grn bsdf_blu bsdf2rad.cal\n");
399 >                        printf("\t-s %f -t %f %f %f\n0\n0\n",
400 >                                        bsdf_rad, sorg[0], sorg[1], sorg[2]);
401 >                        printf("\nscale_pat glow scale_mat\n0\n0\n4 1 1 1 0\n");
402 >                }
403 >                if (front_comp & SDsampR) {
404 >                        put_mirror_arrow(sorg, nrm);
405 >                        fname = tfile_name(frpref, dsuffix, i);
406 >                        sprintf(cmdbuf,
407 >                "gensurf scale_mat %s%d %s %s %s %d %d | xform %s -s %f -t %f %f %f",
408 >                                        frpref, i, fname, fname, fname, SAMPRES-1, SAMPRES-1,
409 >                                        rotargs, scalef, sorg[0], sorg[1], sorg[2]);
410 >                        if (!run_cmd(cmdbuf))
411 >                                return(0);
412 >                }
413 >                if (front_comp & SDsampT) {
414 >                        put_trans_arrow(sorg);
415 >                        fname = tfile_name(ftpref, dsuffix, i);
416 >                        sprintf(cmdbuf,
417 >                "gensurf scale_mat %s%d %s %s %s %d %d | xform -I %s -s %f -t %f %f %f",
418 >                                        ftpref, i, fname, fname, fname, SAMPRES-1, SAMPRES-1,
419 >                                        rotargs, scalef, sorg[0], sorg[1], sorg[2]);
420 >                        if (!run_cmd(cmdbuf))
421 >                                return(0);
422 >                }
423 >                if (back_comp) {
424 >                        cvt_sposition(sorg, ivec, -1);
425 >                        printf("\nvoid colorfunc scale_pat\n");
426 >                        printf("10 bsdf_red bsdf_grn bsdf_blu bsdf2rad.cal\n");
427 >                        printf("\t-s %f -t %f %f %f\n0\n0\n",
428 >                                        bsdf_rad, sorg[0], sorg[1], sorg[2]);
429 >                        printf("\nscale_pat glow scale_mat\n0\n0\n4 1 1 1 0\n");
430 >                }
431 >                if (back_comp & SDsampR) {
432 >                        put_mirror_arrow(sorg, nrm);
433 >                        fname = tfile_name(brpref, dsuffix, i);
434 >                        sprintf(cmdbuf,
435 >                "gensurf scale_mat %s%d %s %s %s %d %d | xform -I -ry 180 %s -s %f -t %f %f %f",
436 >                                        brpref, i, fname, fname, fname, SAMPRES-1, SAMPRES-1,
437 >                                        rotargs, scalef, sorg[0], sorg[1], sorg[2]);
438 >                        if (!run_cmd(cmdbuf))
439 >                                return(0);
440 >                }
441 >                if (back_comp & SDsampT) {
442 >                        put_trans_arrow(sorg);
443 >                        fname = tfile_name(btpref, dsuffix, i);
444 >                        sprintf(cmdbuf,
445 >                "gensurf scale_mat %s%d %s %s %s %d %d | xform -ry 180 %s -s %f -t %f %f %f",
446 >                                        btpref, i, fname, fname, fname, SAMPRES-1, SAMPRES-1,
447 >                                        rotargs, scalef, sorg[0], sorg[1], sorg[2]);
448 >                        if (!run_cmd(cmdbuf))
449 >                                return(0);
450 >                }
451 >        }
452 >        return(1);
453 > }
454 >
455 > /* Put our hemisphere material */
456 > static void
457 > put_matBSDF(const char *XMLfile)
458 > {
459 >        const char      *curdir = "./";
460 >
461 >        if (!XMLfile) {                 /* simple material */
462 >                printf("\n# Simplified material because we have no XML input\n");
463 >                printf("\nvoid brightfunc latlong\n2 latlong bsdf2rad.cal\n0\n0\n");
464 >                if ((front_comp|back_comp) & SDsampT)
465 >                        printf("\nlatlong trans %s\n0\n0\n7 .75 .75 .75 0 .04 .5 .8\n",
466 >                                        sph_fmat);
467 >                else
468 >                        printf("\nlatlong plastic %s\n0\n0\n5 .5 .5 .5 0 0\n",
469 >                                        sph_fmat);
470 >                printf("\ninherit alias %s %s\n", sph_bmat, sph_fmat);
471 >                return;
472 >        }
473 >        switch (XMLfile[0]) {           /* avoid RAYPATH search */
474 >        case '.':
475 >        case '~':
476 >        CASEDIRSEP:
477 >                curdir = "";
478 >                break;
479 >        case '\0':
480 >                fprintf(stderr, "%s: empty file name in put_matBSDF\n", progname);
481 >                exit(1);
482 >                break;
483 >        }
484 >        printf("\n# Actual BSDF materials for rendering the hemispheres\n");
485 >        printf("\nvoid BSDF BSDF_f\n6 0 \"%s%s\" upx upy upz bsdf2rad.cal\n0\n0\n",
486 >                        curdir, XMLfile);
487 >        printf("\nvoid plastic black\n0\n0\n5 0 0 0 0 0\n");
488 >        printf("\nvoid mixfunc %s\n4 BSDF_f black latlong bsdf2rad.cal\n0\n0\n",
489 >                        sph_fmat);
490 >        printf("\nvoid BSDF BSDF_b\n8 0 \"%s%s\" upx upy upz bsdf2rad.cal -ry 180\n0\n0\n",
491 >                        curdir, XMLfile);
492 >        printf("\nvoid mixfunc %s\n4 BSDF_b black latlong bsdf2rad.cal\n0\n0\n",
493 >                        sph_bmat);
494 > }
495 >
496 > /* Put out overhead parallel light source */
497 > static void
498 > put_source(void)
499 > {
500 >        printf("\n# Overhead parallel light source\n");
501 >        printf("\nvoid light bright\n0\n0\n3 2500 2500 2500\n");
502 >        printf("\nbright source light\n0\n0\n4 0 0 1 2\n");
503 >        printf("\n# Material used for labels\n");
504 >        printf("\nvoid trans vellum\n0\n0\n7 1 1 1 0 0 .5 0\n");
505 > }
506 >
507 > /* Put out hemisphere(s) */
508 > static void
509 > put_hemispheres(void)
510 > {
511 >        const int       nsegs = 131;
512 >
513 >        printf("\n# Hemisphere(s) for showing BSDF appearance (if XML file)\n");
514 >        if (front_comp) {
515 >                printf(
516 > "\n!genrev %s Front \"R*sin(A*t)\" \"R*cos(A*t)\" %d -e \"R:%g;A:%f\" -s | xform -t %g 0 0\n",
517 >                                sph_fmat, nsegs, sph_rad, 0.5*PI, sph_xoffset);
518 >                printf("\nvoid brighttext front_text\n3 helvet.fnt . FRONT\n0\n");
519 >                printf("12\n\t%f %f 0\n\t%f 0 0\n\t0 %f 0\n\t.01 1 -.1\n",
520 >                                -.22*sph_rad + sph_xoffset, -1.4*sph_rad,
521 >                                .35/5.*sph_rad, -1.6*.35/5.*sph_rad);
522 >                printf("\nfront_text alias front_label_mat vellum\n");
523 >                printf("\nfront_label_mat polygon front_label\n0\n0\n12");
524 >                printf("\n\t%f %f 0\n\t%f %f 0\n\t%f %f 0\n\t%f %f 0\n",
525 >                                -.25*sph_rad + sph_xoffset, -1.3*sph_rad,
526 >                                -.25*sph_rad + sph_xoffset, (-1.4-1.6*.35/5.-.1)*sph_rad,
527 >                                .25*sph_rad + sph_xoffset, (-1.4-1.6*.35/5.-.1)*sph_rad,
528 >                                .25*sph_rad + sph_xoffset, -1.3*sph_rad );
529 >        }
530 >        if (back_comp) {
531 >                printf(
532 > "\n!genrev %s Back \"R*cos(A*t)\" \"R*sin(A*t)\" %d -e \"R:%g;A:%f\" -s | xform -t %g 0 0\n",
533 >                                sph_bmat, nsegs, sph_rad, 0.5*PI, -sph_xoffset);
534 >                printf("\nvoid brighttext back_text\n3 helvet.fnt . BACK\n0\n");
535 >                printf("12\n\t%f %f 0\n\t%f 0 0\n\t0 %f 0\n\t.01 1 -.1\n",
536 >                                -.22*sph_rad - sph_xoffset, -1.4*sph_rad,
537 >                                .35/4.*sph_rad, -1.6*.35/4.*sph_rad);
538 >                printf("\nback_text alias back_label_mat vellum\n");
539 >                printf("\nback_label_mat polygon back_label\n0\n0\n12");
540 >                printf("\n\t%f %f 0\n\t%f %f 0\n\t%f %f 0\n\t%f %f 0\n",
541 >                                -.25*sph_rad - sph_xoffset, -1.3*sph_rad,
542 >                                -.25*sph_rad - sph_xoffset, (-1.4-1.6*.35/4.-.1)*sph_rad,
543 >                                .25*sph_rad - sph_xoffset, (-1.4-1.6*.35/4.-.1)*sph_rad,
544 >                                .25*sph_rad - sph_xoffset, -1.3*sph_rad );
545 >        }
546 > }
547 >
548 > /* Put out falsecolor scale and name label */
549 > static void
550 > put_scale(void)
551 > {
552 >        const double    max_log10 = log10(overall_max);
553 >        const double    leg_width = 2.*.75*(fabs(sph_xoffset) - sph_rad);
554 >        const double    leg_height = 2.*sph_rad;
555 >        const int       text_lines = 6;
556 >        const int       text_digits = 8;
557 >        char            fmt[16];
558 >        int             i;
559 >
560 >        printf("\n# BSDF legend with falsecolor scale\n");
561 >        printf("\nvoid colorfunc lscale\n10 sca_red(Py) sca_grn(Py) sca_blu(Py)");
562 >        printf("\n\tbsdf2rad.cal -s %f -t 0 %f 0\n0\n0\n", leg_height, -.5*leg_height);
563 >        sprintf(fmt, "%%.%df", text_digits-3);
564 >        for (i = 0; i < text_lines; i++) {
565 >                char    vbuf[16];
566 >                sprintf(vbuf, fmt, pow(10., (i+.5)/text_lines*(max_log10-min_log10)+min_log10));
567 >                printf("\nlscale brighttext lscale\n");
568 >                printf("3 helvet.fnt . %s\n0\n12\n", vbuf);
569 >                printf("\t%f %f 0\n", -.45*leg_width, ((i+.9)/text_lines-.5)*leg_height);
570 >                printf("\t%f 0 0\n", .8*leg_width/strlen(vbuf));
571 >                printf("\t0 %f 0\n", -.9/text_lines*leg_height);
572 >                printf("\t.01 1 -.1\n");
573 >        }
574 >        printf("\nlscale alias legend_mat vellum\n");
575 >        printf("\nlegend_mat polygon legend\n0\n0\n12");
576 >        printf("\n\t%f %f 0\n\t%f %f 0\n\t%f %f 0\n\t%f %f 0\n",
577 >                        -.5*leg_width, .5*leg_height,
578 >                        -.5*leg_width, -.5*leg_height,
579 >                        .5*leg_width, -.5*leg_height,
580 >                        .5*leg_width, .5*leg_height);
581 >        printf("\nvoid brighttext BSDFtitle\n3 helvet.fnt . BSDF\n0\n12\n");
582 >        printf("\t%f %f 0\n", -.25*leg_width, .7*leg_height);
583 >        printf("\t%f 0 0\n", .4/4.*leg_width);
584 >        printf("\t0 %f 0\n", -.1*leg_height);
585 >        printf("\t.01 1 -.1\n");
586 >        printf("\nBSDFtitle alias title_mat vellum\n");
587 >        printf("\ntitle_mat polygon title\n0\n0\n12");
588 >        printf("\n\t%f %f 0\n\t%f %f 0\n\t%f %f 0\n\t%f %f 0\n",
589 >                        -.3*leg_width, .75*leg_height,
590 >                        -.3*leg_width, .55*leg_height,
591 >                        .3*leg_width, .55*leg_height,
592 >                        .3*leg_width, .75*leg_height);
593 >        if (!bsdf_name[0])
594 >                return;
595 >        printf("\nvoid brighttext BSDFname\n3 helvet.fnt . \"%s\"\n0\n12\n", bsdf_name);
596 >        printf("\t%f %f 0\n", -.95*leg_width, -.6*leg_height);
597 >        printf("\t%f 0 0\n", 1.8/strlen(bsdf_name)*leg_width);
598 >        printf("\t0 %f 0\n", -.1*leg_height);
599 >        printf("\t.01 1 -.1\n");
600 >        printf("\nBSDFname alias name_mat vellum\n");
601 >        printf("\nname_mat polygon name\n0\n0\n12");
602 >        printf("\n\t%f %f 0\n\t%f %f 0\n\t%f %f 0\n\t%f %f 0\n",
603 >                        -leg_width, -.55*leg_height,
604 >                        -leg_width, -.75*leg_height,
605 >                        leg_width, -.75*leg_height,
606 >                        leg_width, -.55*leg_height);
607 > }
608 >
609 > /* Convert MGF to Radiance in output */
610 > static void
611 > convert_mgf(const char *mgfdata)
612 > {
613 >        int     len = strlen(mgfdata);
614 >        char    mgfn[128];
615 >        char    radfn[128];
616 >        char    cmdbuf[256];
617 >        float   xmin, xmax, ymin, ymax, zmin, zmax;
618 >        double  max_dim;
619 >        int     fd;
620 >        FILE    *fp;
621 >
622 >        if (!len) return;
623 >        strcpy(mgfn, tfile_name("geom", ".mgf", 0));
624 >        fd = open(mgfn, O_WRONLY|O_CREAT, 0666);
625 >        if (fd < 0 || write(fd, mgfdata, len) != len) {
626 >                fprintf(stderr, "%s: cannot write file '%s'\n",
627 >                                progname, mgfn);
628 >                return;
629 >        }
630 >        close(fd);
631 >        strcpy(radfn, tfile_name("geom", ".rad", 0));
632 >        sprintf(cmdbuf, "mgf2rad %s > %s", mgfn, radfn);
633 >        if (!run_cmd(cmdbuf))
634 >                return;
635 >        sprintf(cmdbuf, "getbbox -w -h %s", radfn);
636 >        if ((fp = popen(cmdbuf, "r")) == NULL ||
637 >                        fscanf(fp, "%f %f %f %f %f %f",
638 >                                &xmin, &xmax, &ymin, &ymax, &zmin, &zmax) != 6
639 >                        || pclose(fp) != 0) {
640 >                fprintf(stderr, "%s: error reading from command: %s\n",
641 >                                progname, cmdbuf);
642 >                return;
643 >        }
644 >        max_dim = ymax - ymin;
645 >        if (xmax - xmin > max_dim)
646 >                max_dim = xmax - xmin;
647 >        if (front_comp) {
648 >                printf("\n# BSDF system geometry (front view)\n");
649 >                sprintf(cmdbuf, "xform -t %f %f %f -s %f -t %f %f 0 %s",
650 >                                -.5*(xmin+xmax), -.5*(ymin+ymax), -zmax,
651 >                                1.5*sph_rad/max_dim,
652 >                                sph_xoffset, -2.5*sph_rad,
653 >                                radfn);
654 >                if (!run_cmd(cmdbuf))
655 >                        return;
656 >        }
657 >        if (back_comp) {
658 >                printf("\n# BSDF system geometry (back view)\n");
659 >                sprintf(cmdbuf, "xform -t %f %f %f -s %f -ry 180 -t %f %f 0 %s",
660 >                                -.5*(xmin+xmax), -.5*(ymin+ymax), -zmin,
661 >                                1.5*sph_rad/max_dim,
662 >                                -sph_xoffset, -2.5*sph_rad,
663 >                                radfn);
664 >                if (!run_cmd(cmdbuf))
665 >                        return;
666 >        }
667 > }
668 >
669 > /* Check RBF input header line & get minimum BSDF value */
670 > static int
671 > rbf_headline(char *s, void *p)
672 > {
673 >        char    fmt[MAXFMTLEN];
674 >
675 >        if (formatval(fmt, s)) {
676 >                if (strcmp(fmt, BSDFREP_FMT))
677 >                        return(-1);
678 >                return(0);
679 >        }
680 >        if (!strncmp(s, "IO_SIDES=", 9)) {
681 >                sscanf(s+9, "%d %d", &input_orient, &output_orient);
682 >                if (input_orient == output_orient) {
683 >                        if (input_orient > 0)
684 >                                front_comp |= SDsampR;
685 >                        else
686 >                                back_comp |= SDsampR;
687 >                } else if (input_orient > 0)
688 >                        front_comp |= SDsampT;
689 >                else
690 >                        back_comp |= SDsampT;
691 >                return(0);
692 >        }
693 >        if (!strncmp(s, "BSDFMIN=", 8)) {
694 >                sscanf(s+8, "%lf", &bsdf_min);
695 >                if (bsdf_min < overall_min)
696 >                        overall_min = bsdf_min;
697 >                return(0);
698 >        }
699 >        return(0);
700 > }
701 >
702 > /* Produce a Radiance model plotting the given BSDF representation */
703   int
704   main(int argc, char *argv[])
705   {
31        int     showPeaks = 0;
32        int     doTrans = 0;
706          int     inpXML = -1;
707 <        RBFNODE *rbf = NULL;
35 <        FILE    *fp;
36 <        char    buf[128];
707 >        double  myLim[2];
708          SDData  myBSDF;
709 <        double  bsdf, min_log;
710 <        FVECT   idir, odir;
711 <        int     i, j, n;
709 >        int     a, n;
710 >                                                /* set global progname */
711 >        fixargv0(argv[0]);
712                                                  /* check arguments */
713 <        progname = argv[0];
714 <        if (argc > 1 && !strcmp(argv[1], "-p")) {
715 <                ++showPeaks;
716 <                ++argv; --argc;
713 >        a = 1;
714 >        myLim[0] = -1; myLim[1] = -2;           /* specified BSDF range? */
715 >        if (argc > a+3 && argv[a][0] == '-' && argv[a][1] == 'r') {
716 >                myLim[0] = atof(argv[++a]);
717 >                myLim[1] = atof(argv[++a]);
718 >                ++a;
719          }
720 <        if (argc > 1 && !strcmp(argv[1], "-t")) {
721 <                ++doTrans;
49 <                ++argv; --argc;
50 <        }
51 <        if (argc >= 4 && (n = strlen(argv[1])-4) > 0) {
52 <                if (!strcasecmp(argv[1]+n, ".xml"))
720 >        if (argc > a && (n = strlen(argv[a])-4) > 0) {
721 >                if (!strcasecmp(argv[a]+n, ".xml"))
722                          inpXML = 1;
723 <                else if (!strcasecmp(argv[1]+n, ".sir"))
723 >                else if (!strcasecmp(argv[a]+n, ".sir"))
724                          inpXML = 0;
725          }
726 <        if (inpXML < 0) {
727 <                fprintf(stderr, "Usage: %s [-p] bsdf.sir theta1 phi1 .. > output.rad\n", progname);
728 <                fprintf(stderr, "   Or: %s [-t] bsdf.xml theta1 phi1 .. > output.rad\n", progname);
726 >        if (inpXML < 0 || inpXML & (argc > a+1)) {
727 >                fprintf(stderr, "Usage: %s [-r min max] bsdf.xml > output.rad\n", progname);
728 >                fprintf(stderr, "   Or: %s [-r min max] hemi1.sir hemi2.sir .. > output.rad\n", progname);
729                  return(1);
730          }
731 <                                                /* load input */
731 >        fputs("# ", stdout);                    /* copy our command */
732 >        printargs(argc, argv, stdout);
733 >                                                /* evaluate BSDF */
734          if (inpXML) {
735 <                SDclearBSDF(&myBSDF, argv[1]);
736 <                if (SDreportError(SDloadFile(&myBSDF, argv[1]), stderr))
735 >                SDclearBSDF(&myBSDF, argv[a]);
736 >                if (SDreportError(SDloadFile(&myBSDF, argv[a]), stderr))
737                          return(1);
738 <                bsdf_min = 1./M_PI;
739 <                if (myBSDF.rf != NULL && myBSDF.rLambFront.cieY < bsdf_min*M_PI)
740 <                        bsdf_min = myBSDF.rLambFront.cieY/M_PI;
741 <                if (myBSDF.rb != NULL && myBSDF.rLambBack.cieY < bsdf_min*M_PI)
742 <                        bsdf_min = myBSDF.rLambBack.cieY/M_PI;
743 <                if ((myBSDF.tf != NULL) | (myBSDF.tb != NULL) &&
744 <                                myBSDF.tLamb.cieY < bsdf_min*M_PI)
74 <                        bsdf_min = myBSDF.tLamb.cieY/M_PI;
75 <                if (doTrans && (myBSDF.tf == NULL) & (myBSDF.tb == NULL)) {
76 <                        fprintf(stderr, "%s: no transmitted component in '%s'\n",
77 <                                        progname, argv[1]);
738 >                if (myBSDF.rf != NULL) front_comp |= SDsampR;
739 >                if (myBSDF.tf != NULL) front_comp |= SDsampT;
740 >                if (myBSDF.rb != NULL) back_comp |= SDsampR;
741 >                if (myBSDF.tb != NULL) back_comp |= SDsampT;
742 >                if (!front_comp & !back_comp) {
743 >                        fprintf(stderr, "%s: nothing to plot in '%s'\n",
744 >                                        progname, argv[a]);
745                          return(1);
746                  }
747 <        } else {
748 <                fp = fopen(argv[1], "rb");
749 <                if (fp == NULL) {
750 <                        fprintf(stderr, "%s: cannot open BSDF interpolant '%s'\n",
751 <                                        progname, argv[1]);
752 <                        return(1);
747 >                if (myLim[0] >= 0)
748 >                        overall_min = myLim[0];
749 >                else {
750 >                        if (front_comp & SDsampR && myBSDF.rLambFront.cieY < overall_min*PI)
751 >                                overall_min = myBSDF.rLambFront.cieY/PI;
752 >                        if (back_comp & SDsampR && myBSDF.rLambBack.cieY < overall_min*PI)
753 >                                overall_min = myBSDF.rLambBack.cieY/PI;
754 >                        if (front_comp & SDsampT && myBSDF.tLambFront.cieY < overall_min*PI)
755 >                                overall_min = myBSDF.tLambFront.cieY/PI;
756 >                        if (back_comp & SDsampT && myBSDF.tLambBack.cieY < overall_min*PI)
757 >                                overall_min = myBSDF.tLambBack.cieY/PI;
758                  }
759 <                if (!load_bsdf_rep(fp))
759 >                set_minlog();
760 >                if (!build_wBSDF(&myBSDF))
761                          return(1);
762 <                fclose(fp);
763 <        }
764 < #ifdef DEBUG
765 <        fprintf(stderr, "Minimum BSDF set to %.4f\n", bsdf_min);
766 < #endif
767 <        min_log = log(bsdf_min*.5 + 1e-5);
768 <                                                /* output BSDF rep. */
769 <        for (n = 0; (n < 6) & (2*n+3 < argc); n++) {
770 <                double  theta = (M_PI/180.)*atof(argv[2*n+2]);
771 <                double  phi = (M_PI/180.)*atof(argv[2*n+3]);
99 <                if (theta < -FTINY) {
100 <                        fprintf(stderr, "%s: theta values must be positive\n",
101 <                                        progname);
762 >                if (myBSDF.matn[0])
763 >                        strcpy(bsdf_name, myBSDF.matn);
764 >                else
765 >                        strcpy(bsdf_name, myBSDF.name);
766 >                strcpy(bsdf_manuf, myBSDF.makr);
767 >                put_matBSDF(argv[a]);
768 >        } else {
769 >                FILE    *fp[4];
770 >                if (argc > a+4) {
771 >                        fprintf(stderr, "%s: more than 4 hemispheres!\n", progname);
772                          return(1);
773                  }
774 <                if (inpXML) {
775 <                        input_orient = (theta <= M_PI/2.) ? 1 : -1;
776 <                        output_orient = doTrans ? -input_orient : input_orient;
774 >                for (n = a; n < argc; n++) {
775 >                        fp[n-a] = fopen(argv[n], "rb");
776 >                        if (fp[n-a] == NULL) {
777 >                                fprintf(stderr, "%s: cannot open BSDF interpolant '%s'\n",
778 >                                                progname, argv[n]);
779 >                                return(1);
780 >                        }
781 >                        if (getheader(fp[n-a], rbf_headline, NULL) < 0) {
782 >                                fprintf(stderr, "%s: bad BSDF interpolant '%s'\n",
783 >                                                progname, argv[n]);
784 >                                return(1);
785 >                        }
786                  }
787 <                idir[2] = sin(theta);
788 <                idir[0] = idir[2] * cos(phi);
789 <                idir[1] = idir[2] * sin(phi);
790 <                idir[2] = input_orient * sqrt(1. - idir[2]*idir[2]);
791 < #ifdef DEBUG
792 <                fprintf(stderr, "Computing BSDF for incident direction (%.1f,%.1f)\n",
793 <                                get_theta180(idir), get_phi360(idir));
794 < #endif
116 <                if (!inpXML)
117 <                        rbf = advect_rbf(idir, 15000);
118 < #ifdef DEBUG
119 <                if (inpXML)
120 <                        fprintf(stderr, "Hemispherical %s: %.3f\n",
121 <                                (output_orient > 0 ? "reflection" : "transmission"),
122 <                                SDdirectHemi(idir, SDsampSp|SDsampDf |
123 <                                                (output_orient > 0 ?
124 <                                                 SDsampR : SDsampT), &myBSDF));
125 <                else if (rbf == NULL)
126 <                        fputs("Empty RBF\n", stderr);
127 <                else
128 <                        fprintf(stderr, "Hemispherical %s: %.3f\n",
129 <                                (output_orient > 0 ? "reflection" : "transmission"),
130 <                                rbf->vtotal);
131 < #endif
132 <                printf("void trans tmat\n0\n0\n7 %f %f %f .04 .04 .9 1\n",
133 <                                colarr[n][0], colarr[n][1], colarr[n][2]);
134 <                if (showPeaks && rbf != NULL) {
135 <                        printf("void plastic pmat\n0\n0\n5 %f %f %f .04 .08\n",
136 <                                1.-colarr[n][0], 1.-colarr[n][1], 1.-colarr[n][2]);
137 <                        for (i = 0; i < rbf->nrbf; i++) {
138 <                                ovec_from_pos(odir, rbf->rbfa[i].gx, rbf->rbfa[i].gy);
139 <                                bsdf = eval_rbfrep(rbf, odir) / (output_orient*odir[2]);
140 <                                bsdf = log(bsdf + 1e-5) - min_log;
141 <                                printf("pmat sphere p%d\n0\n0\n4 %f %f %f %f\n",
142 <                                        i+1, odir[0]*bsdf, odir[1]*bsdf, odir[2]*bsdf,
143 <                                                .007*bsdf);
787 >                if (myLim[0] >= 0)
788 >                        overall_min = myLim[0];
789 >                set_minlog();
790 >                for (n = a; n < argc; n++) {
791 >                        if (fseek(fp[n-a], 0L, SEEK_SET) < 0) {
792 >                                fprintf(stderr, "%s: cannot seek on '%s'\n",
793 >                                                progname, argv[n]);
794 >                                return(1);
795                          }
796 +                        if (!load_bsdf_rep(fp[n-a]))
797 +                                return(1);
798 +                        fclose(fp[n-a]);
799 +                        if (!build_wRBF())
800 +                                return(1);
801                  }
802 <                fflush(stdout);
147 <                sprintf(buf, "gensurf tmat bsdf%d - - - %d %d", n+1,
148 <                                                GRIDRES-1, GRIDRES-1);
149 <                fp = popen(buf, "w");
150 <                if (fp == NULL) {
151 <                        fprintf(stderr, "%s: cannot open '| %s'\n", progname, buf);
152 <                        return(1);
153 <                }
154 <                for (i = 0; i < GRIDRES; i++)
155 <                    for (j = 0; j < GRIDRES; j++) {
156 <                        ovec_from_pos(odir, i, j);
157 <                        if (inpXML) {
158 <                                SDValue sval;
159 <                                if (SDreportError(SDevalBSDF(&sval, odir,
160 <                                                        idir, &myBSDF), stderr))
161 <                                        return(1);
162 <                                bsdf = sval.cieY;
163 <                        } else
164 <                                bsdf = eval_rbfrep(rbf, odir) /
165 <                                                (output_orient*odir[2]);
166 <                        bsdf = log(bsdf + 1e-5) - min_log;
167 <                        fprintf(fp, "%.8e %.8e %.8e\n",
168 <                                        odir[0]*bsdf, odir[1]*bsdf, odir[2]*bsdf);
169 <                    }
170 <                if (rbf != NULL)
171 <                        free(rbf);
172 <                if (pclose(fp))
173 <                        return(1);
802 >                put_matBSDF(NULL);
803          }
804 +        if (myLim[1] > myLim[0])        /* override maximum BSDF? */
805 +                overall_max = myLim[1];
806 +        put_source();                   /* before hemispheres & labels */
807 +        put_hemispheres();
808 +        put_scale();
809 +        if (inpXML && myBSDF.mgf)
810 +                convert_mgf(myBSDF.mgf);
811 +        if (!put_BSDFs())               /* most of the output happens here */
812 +                return(1);
813 +        cleanup_tmp();
814          return(0);
815   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines