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

Comparing ray/src/cv/bsdf2klems.c (file contents):
Revision 2.6 by greg, Sat Jun 29 21:03:25 2013 UTC vs.
Revision 2.38 by greg, Sat Jun 7 05:09:45 2025 UTC

# Line 8 | Line 8 | static const char RCSid[] = "$Id$";
8   */
9  
10   #define _USE_MATH_DEFINES
11 #include <stdio.h>
11   #include <stdlib.h>
13 #include <string.h>
12   #include <math.h>
13 + #include <ctype.h>
14   #include "random.h"
15 < #include "platform.h"
15 > #include "rtio.h"
16   #include "calcomp.h"
17   #include "bsdfrep.h"
18   #include "bsdf_m.h"
19 +                                /* tristimulus components */
20 + enum {CIE_X, CIE_Y, CIE_Z};
21                                  /* assumed maximum # Klems patches */
22   #define MAXPATCHES      145
22                                /* global argv[0] */
23 char                    *progname;
23                                  /* selected basis function name */
24 < static const char       *kbasis = "LBNL/Klems Full";
24 > static const char       klems_full[] = "LBNL/Klems Full";
25 > static const char       klems_half[] = "LBNL/Klems Half";
26 > static const char       klems_quarter[] = "LBNL/Klems Quarter";
27 > static const char       *kbasis = klems_full;
28                                  /* number of BSDF samples per patch */
29 < static int              npsamps = 256;
29 > static int              npsamps = 1024;
30 >                                /* limit on number of RBF lobes */
31 > static int              lobe_lim = 15000;
32 >                                /* progress bar length */
33 > static int              do_prog = 79;
34  
35 < /* Return angle basis corresponding to the given name */
36 < ANGLE_BASIS *
37 < get_basis(const char *bn)
32 < {
33 <        int     n = nabases;
34 <        
35 <        while (n-- > 0)
36 <                if (!strcasecmp(bn, abase_list[n].name))
37 <                        return &abase_list[n];
38 <        return NULL;
39 < }
35 > #define MAXCARG         512     /* wrapBSDF command */
36 > static char             *wrapBSDF[MAXCARG] = {"wrapBSDF", "-W", "-UU"};
37 > static int              wbsdfac = 3;
38  
39 < /* Output XML header to stdout */
39 > /* Add argument to wrapBSDF, allocating space if !isstatic */
40   static void
41 < xml_header(int ac, char *av[])
41 > add_wbsdf(const char *arg, int isstatic)
42   {
43 <        puts("<?xml version=\"1.0\" encoding=\"UTF-8\"?>");
44 <        puts("<WindowElement xmlns=\"http://windows.lbl.gov\" xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\" xsi:schemaLocation=\"http://windows.lbl.gov/BSDF-v1.4.xsd\">");
45 <        fputs("<!-- File produced by:", stdout);
46 <        while (ac-- > 0) {
47 <                fputc(' ', stdout);
48 <                fputs(*av++, stdout);
43 >        if (arg == NULL)
44 >                return;
45 >        if (wbsdfac >= MAXCARG-1) {
46 >                fputs(progname, stderr);
47 >                fputs(": too many command arguments to wrapBSDF\n", stderr);
48 >                exit(1);
49          }
50 <        puts(" -->");
50 >        if (!*arg)
51 >                arg = "";
52 >        else if (!isstatic)
53 >                arg = savqstr((char *)arg);
54 >
55 >        wrapBSDF[wbsdfac++] = (char *)arg;
56   }
57  
58 < /* Output XML prologue to stdout */
58 > /* Start new progress bar */
59 > #define prog_start(s)   if (do_prog) fprintf(stderr, "%s: %s...\n", progname, s); else
60 >
61 > /* Draw progress bar of the appropriate length */
62   static void
63 < xml_prologue(const SDData *sd)
63 > prog_show(double frac)
64   {
65 <        const char      *matn = (sd && sd->matn[0]) ? sd->matn : "Name";
66 <        const char      *makr = (sd && sd->makr[0]) ? sd->makr : "Manufacturer";
67 <        ANGLE_BASIS     *abp = get_basis(kbasis);
68 <        int             i;
65 >        static unsigned call_cnt = 0;
66 >        static char     lastc[] = "-\\|/";
67 >        char            pbar[256];
68 >        int             nchars;
69  
70 <        if (abp == NULL) {
71 <                fprintf(stderr, "%s: unknown angle basis '%s'\n", progname, kbasis);
72 <                exit(1);
73 <        }
74 <        puts("<WindowElementType>System</WindowElementType>");
75 <        puts("<FileType>BSDF</FileType>");
76 <        puts("<Optical>");
77 <        puts("<Layer>");
78 <        puts("\t<Material>");
79 <        printf("\t\t<Name>%s</Name>\n", matn);
80 <        printf("\t\t<Manufacturer>%s</Manufacturer>\n", makr);
81 <        if (sd && sd->dim[2] > .001) {
76 <                printf("\t\t<Thickness unit=\"meter\">%.3f</Thickness>\n", sd->dim[2]);
77 <                printf("\t\t<Width unit=\"meter\">%.3f</Width>\n", sd->dim[0]);
78 <                printf("\t\t<Height unit=\"meter\">%.3f</Height>\n", sd->dim[1]);
79 <        }
80 <        puts("\t\t<DeviceType>Other</DeviceType>");
81 <        puts("\t</Material>");
82 <        if (sd && sd->mgf != NULL) {
83 <                puts("\t<Geometry format=\"MGF\">");
84 <                puts("\t\t<MGFblock unit=\"meter\">");
85 <                fputs(sd->mgf, stdout);
86 <                puts("</MGFblock>");
87 <                puts("\t</Geometry>");
88 <        }
89 <        puts("\t<DataDefinition>");
90 <        puts("\t\t<IncidentDataStructure>Columns</IncidentDataStructure>");
91 <        puts("\t\t<AngleBasis>");
92 <        printf("\t\t\t<AngleBasisName>%s</AngleBasisName>\n", kbasis);
93 <        for (i = 0; abp->lat[i].nphis; i++) {
94 <                puts("\t\t\t<AngleBasisBlock>");
95 <                printf("\t\t\t<Theta>%g</Theta>\n", i ?
96 <                                .5*(abp->lat[i].tmin + abp->lat[i+1].tmin) :
97 <                                .0 );
98 <                printf("\t\t\t<nPhis>%d</nPhis>\n", abp->lat[i].nphis);
99 <                puts("\t\t\t<ThetaBounds>");
100 <                printf("\t\t\t\t<LowerTheta>%g</LowerTheta>\n", abp->lat[i].tmin);
101 <                printf("\t\t\t\t<UpperTheta>%g</UpperTheta>\n", abp->lat[i+1].tmin);
102 <                puts("\t\t\t</ThetaBounds>");
103 <                puts("\t\t\t</AngleBasisBlock>");
104 <        }
105 <        puts("\t\t</AngleBasis>");
106 <        puts("\t</DataDefinition>");
70 >        if (do_prog <= 1) return;
71 >        if (do_prog > sizeof(pbar)-2)
72 >                do_prog = sizeof(pbar)-2;
73 >        if (frac < 0) frac = 0;
74 >        else if (frac >= 1) frac = .9999;
75 >        nchars = do_prog*frac;
76 >        pbar[0] = '\r';
77 >        memset(pbar+1, '*', nchars);
78 >        pbar[nchars+1] = lastc[call_cnt++ & 3];
79 >        memset(pbar+2+nchars, '-', do_prog-nchars-1);
80 >        pbar[do_prog+1] = '\0';
81 >        fputs(pbar, stderr);
82   }
83  
84 < /* Output XML data prologue to stdout */
84 > /* Finish progress bar */
85   static void
86 < data_prologue()
86 > prog_done(void)
87   {
88 <        static const char       *bsdf_type[4] = {
114 <                                        "Reflection Front",
115 <                                        "Transmission Front",
116 <                                        "Transmission Back",
117 <                                        "Reflection Back"
118 <                                };
88 >        int     n = do_prog;
89  
90 <        puts("\t<WavelengthData>");
91 <        puts("\t\t<LayerNumber>System</LayerNumber>");
92 <        puts("\t\t<Wavelength unit=\"Integral\">Visible</Wavelength>");
93 <        puts("\t\t<SourceSpectrum>CIE Illuminant D65 1nm.ssp</SourceSpectrum>");
94 <        puts("\t\t<DetectorSpectrum>ASTM E308 1931 Y.dsp</DetectorSpectrum>");
125 <        puts("\t\t<WavelengthDataBlock>");
126 <        printf("\t\t\t<WavelengthDataDirection>%s</WavelengthDataDirection>\n",
127 <                        bsdf_type[(input_orient>0)<<1 | (output_orient>0)]);
128 <        printf("\t\t\t<ColumnAngleBasis>%s</ColumnAngleBasis>\n", kbasis);
129 <        printf("\t\t\t<RowAngleBasis>%s</RowAngleBasis>\n", kbasis);
130 <        puts("\t\t\t<ScatteringDataType>BTDF</ScatteringDataType>");
131 <        puts("\t\t\t<ScatteringData>");
90 >        if (n <= 1) return;
91 >        fputc('\r', stderr);
92 >        while (n--)
93 >                fputc(' ', stderr);
94 >        fputc('\r', stderr);
95   }
96  
97 < /* Output XML data epilogue to stdout */
98 < static void
99 < data_epilogue(void)
97 > /* Return angle basis corresponding to the given name */
98 > static ANGLE_BASIS *
99 > get_basis(const char *bn)
100   {
101 <        puts("\t\t\t</ScatteringData>");
102 <        puts("\t\t</WavelengthDataBlock>");
103 <        puts("\t</WavelengthData>");
101 >        int     n = nabases;
102 >        
103 >        while (n-- > 0)
104 >                if (!strcasecmp(bn, abase_list[n].name))
105 >                        return &abase_list[n];
106 >        return NULL;
107   }
108  
109 < /* Output XML epilogue to stdout */
110 < static void
111 < xml_epilogue(void)
109 > /* Copy geometry string to file for wrapBSDF */
110 > static char *
111 > save_geom(const char *mgf)
112   {
113 <        puts("</Layer>");
114 <        puts("</Optical>");
115 <        puts("</WindowElement>");
113 >        char    *tfname = mktemp(savqstr(TEMPLATE));
114 >        int     fd = open(tfname, O_CREAT|O_WRONLY, 0600);
115 >
116 >        if (fd < 0)
117 >                return(NULL);
118 >        write(fd, mgf, strlen(mgf));
119 >        close(fd);
120 >        add_wbsdf("-g", 1);
121 >        add_wbsdf(tfname, 1);
122 >        return(tfname);
123   }
124  
125 + /* Open XYZ component file for output and add appropriate arguments */
126 + static FILE *
127 + open_component_file(int c)
128 + {
129 +        static const char       sname[3][6] = {"CIE-X", "CIE-Y", "CIE-Z"};
130 +        static const char       cname[4][4] = {"-rf", "-tf", "-tb", "-rb"};
131 +        char                    *tfname = mktemp(savqstr(TEMPLATE));
132 +        FILE                    *fp = fopen(tfname, "w");
133 +
134 +        if (fp == NULL) {
135 +                fprintf(stderr, "%s: cannot open '%s' for writing\n",
136 +                                progname, tfname);
137 +                exit(1);
138 +        }
139 +        add_wbsdf("-s", 1); add_wbsdf(sname[c], 1);
140 +        add_wbsdf(cname[(input_orient>0)<<1 | (output_orient>0)], 1);
141 +        add_wbsdf(tfname, 1);
142 +        return(fp);
143 + }
144 +
145   /* Load and resample XML BSDF description using Klems basis */
146   static void
147   eval_bsdf(const char *fname)
148   {
149          ANGLE_BASIS     *abp = get_basis(kbasis);
150 +        FILE            *cfp[3];
151          SDData          bsd;
152          SDError         ec;
153          FVECT           vin, vout;
154 <        SDValue         sv;
155 <        double          sum;
154 >        SDValue         sdv;
155 >        double          sum, xsum, ysum;
156          int             i, j, n;
157  
158 +        initurand(npsamps);
159          SDclearBSDF(&bsd, fname);               /* load BSDF file */
160          if ((ec = SDloadFile(&bsd, fname)) != SDEnone)
161                  goto err;
162 <        xml_prologue(&bsd);                     /* pass geometry */
162 >        if (bsd.mgf != NULL)                    /* save geometry */
163 >                save_geom(bsd.mgf);
164 >        if (bsd.matn[0])                        /* save identifier(s) */
165 >                strcpy(bsdf_name, bsd.matn);
166 >        if (bsd.makr[0])
167 >                strcpy(bsdf_manuf, bsd.makr);
168 >        if (bsd.dim[2] > 0) {                   /* save dimension(s) */
169 >                char    buf[64];
170 >                if ((bsd.dim[0] > 0) & (bsd.dim[1] > 0))
171 >                        sprintf(buf, "w=%g;h=%g;t=%g",
172 >                                        bsd.dim[0], bsd.dim[1], bsd.dim[2]);
173 >                else
174 >                        sprintf(buf, "t=%g", bsd.dim[2]);
175 >                add_wbsdf("-f", 1);
176 >                add_wbsdf(buf, 0);
177 >        }
178                                                  /* front reflection */
179          if (bsd.rf != NULL || bsd.rLambFront.cieY > .002) {
180              input_orient = 1; output_orient = 1;
181 <            data_prologue();
181 >            cfp[CIE_Y] = open_component_file(CIE_Y);
182 >            if (bsd.rf != NULL && bsd.rf->comp[0].cspec[2].flags) {
183 >                rbf_colorimetry = RBCtristimulus;
184 >                cfp[CIE_X] = open_component_file(CIE_X);
185 >                cfp[CIE_Z] = open_component_file(CIE_Z);
186 >            } else
187 >                rbf_colorimetry = RBCphotopic;
188              for (j = 0; j < abp->nangles; j++) {
189                  for (i = 0; i < abp->nangles; i++) {
190                      sum = 0;                    /* average over patches */
191 +                    xsum = ysum = 0;
192                      for (n = npsamps; n-- > 0; ) {
193                          fo_getvec(vout, j+(n+frandom())/npsamps, abp);
194                          fi_getvec(vin, i+urand(n), abp);
195 <                        ec = SDevalBSDF(&sv, vout, vin, &bsd);
195 >                        ec = SDevalBSDF(&sdv, vin, vout, &bsd);
196                          if (ec != SDEnone)
197                                  goto err;
198 <                        sum += sv.cieY;
198 >                        sum += sdv.cieY;
199 >                        if (rbf_colorimetry == RBCtristimulus) {
200 >                                xsum += sdv.cieY * sdv.spec.cx;
201 >                                ysum += sdv.cieY * sdv.spec.cy;
202 >                        }
203                      }
204 <                    printf("\t%.3e\n", sum/npsamps);
204 >                    fprintf(cfp[CIE_Y], "\t%.3e\n", sum/npsamps);
205 >                    if (rbf_colorimetry == RBCtristimulus) {
206 >                        fprintf(cfp[CIE_X], "\t%.3e\n", xsum*sum/(npsamps*ysum));
207 >                        fprintf(cfp[CIE_Z], "\t%.3e\n",
208 >                                (sum - xsum - ysum)*sum/(npsamps*ysum));
209 >                    }
210                  }
211 <                putchar('\n');                  /* extra space between rows */
211 >                fputc('\n', cfp[CIE_Y]);        /* extra space between rows */
212 >                if (rbf_colorimetry == RBCtristimulus) {
213 >                        fputc('\n', cfp[CIE_X]);
214 >                        fputc('\n', cfp[CIE_Z]);
215 >                }
216              }
217 <            data_epilogue();
217 >            if (fclose(cfp[CIE_Y])) {
218 >                fprintf(stderr, "%s: error writing Y output\n", progname);
219 >                exit(1);
220 >            }
221 >            if (rbf_colorimetry == RBCtristimulus &&
222 >                        (fclose(cfp[CIE_X]) || fclose(cfp[CIE_Z]))) {
223 >                fprintf(stderr, "%s: error writing X/Z output\n", progname);
224 >                exit(1);
225 >            }
226          }
227                                                  /* back reflection */
228          if (bsd.rb != NULL || bsd.rLambBack.cieY > .002) {
229              input_orient = -1; output_orient = -1;
230 <            data_prologue();
230 >            cfp[CIE_Y] = open_component_file(CIE_Y);
231 >            if (bsd.rb != NULL && bsd.rb->comp[0].cspec[2].flags) {
232 >                rbf_colorimetry = RBCtristimulus;
233 >                cfp[CIE_X] = open_component_file(CIE_X);
234 >                cfp[CIE_Z] = open_component_file(CIE_Z);
235 >            } else
236 >                rbf_colorimetry = RBCphotopic;
237              for (j = 0; j < abp->nangles; j++) {
238                  for (i = 0; i < abp->nangles; i++) {
239                      sum = 0;                    /* average over patches */
240 +                    xsum = ysum = 0;
241                      for (n = npsamps; n-- > 0; ) {
242                          bo_getvec(vout, j+(n+frandom())/npsamps, abp);
243                          bi_getvec(vin, i+urand(n), abp);
244 <                        ec = SDevalBSDF(&sv, vout, vin, &bsd);
244 >                        ec = SDevalBSDF(&sdv, vin, vout, &bsd);
245                          if (ec != SDEnone)
246                                  goto err;
247 <                        sum += sv.cieY;
247 >                        sum += sdv.cieY;
248 >                        if (rbf_colorimetry == RBCtristimulus) {
249 >                                xsum += sdv.cieY * sdv.spec.cx;
250 >                                ysum += sdv.cieY * sdv.spec.cy;
251 >                        }
252                      }
253 <                    printf("\t%.3e\n", sum/npsamps);
253 >                    fprintf(cfp[CIE_Y], "\t%.3e\n", sum/npsamps);
254 >                    if (rbf_colorimetry == RBCtristimulus) {
255 >                        fprintf(cfp[CIE_X], "\t%.3e\n", xsum*sum/(npsamps*ysum));
256 >                        fprintf(cfp[CIE_Z], "\t%.3e\n",
257 >                                (sum - xsum - ysum)*sum/(npsamps*ysum));
258 >                    }
259                  }
260 <                putchar('\n');                  /* extra space between rows */
260 >                if (rbf_colorimetry == RBCtristimulus) {
261 >                        fputc('\n', cfp[CIE_X]);
262 >                        fputc('\n', cfp[CIE_Z]);
263 >                }
264              }
265 <            data_epilogue();
265 >            if (fclose(cfp[CIE_Y])) {
266 >                fprintf(stderr, "%s: error writing Y output\n", progname);
267 >                exit(1);
268 >            }
269 >            if (rbf_colorimetry == RBCtristimulus &&
270 >                        (fclose(cfp[CIE_X]) || fclose(cfp[CIE_Z]))) {
271 >                fprintf(stderr, "%s: error writing X/Z output\n", progname);
272 >                exit(1);
273 >            }
274          }
275                                                  /* front transmission */
276 <        if (bsd.tf != NULL || bsd.tLamb.cieY > .002) {
276 >        if (bsd.tf != NULL || bsd.tLambFront.cieY > .002) {
277              input_orient = 1; output_orient = -1;
278 <            data_prologue();
278 >            cfp[CIE_Y] = open_component_file(CIE_Y);
279 >            if (bsd.tf != NULL && bsd.tf->comp[0].cspec[2].flags) {
280 >                rbf_colorimetry = RBCtristimulus;
281 >                cfp[CIE_X] = open_component_file(CIE_X);
282 >                cfp[CIE_Z] = open_component_file(CIE_Z);
283 >            } else
284 >                rbf_colorimetry = RBCphotopic;
285              for (j = 0; j < abp->nangles; j++) {
286                  for (i = 0; i < abp->nangles; i++) {
287                      sum = 0;                    /* average over patches */
288 +                    xsum = ysum = 0;
289                      for (n = npsamps; n-- > 0; ) {
290                          bo_getvec(vout, j+(n+frandom())/npsamps, abp);
291                          fi_getvec(vin, i+urand(n), abp);
292 <                        ec = SDevalBSDF(&sv, vout, vin, &bsd);
292 >                        ec = SDevalBSDF(&sdv, vin, vout, &bsd);
293                          if (ec != SDEnone)
294                                  goto err;
295 <                        sum += sv.cieY;
295 >                        sum += sdv.cieY;
296 >                        if (rbf_colorimetry == RBCtristimulus) {
297 >                                xsum += sdv.cieY * sdv.spec.cx;
298 >                                ysum += sdv.cieY * sdv.spec.cy;
299 >                        }
300                      }
301 <                    printf("\t%.3e\n", sum/npsamps);
301 >                    fprintf(cfp[CIE_Y], "\t%.3e\n", sum/npsamps);
302 >                    if (rbf_colorimetry == RBCtristimulus) {
303 >                        fprintf(cfp[CIE_X], "\t%.3e\n", xsum*sum/(npsamps*ysum));
304 >                        fprintf(cfp[CIE_Z], "\t%.3e\n",
305 >                                (sum - xsum - ysum)*sum/(npsamps*ysum));
306 >                    }
307                  }
308 <                putchar('\n');                  /* extra space between rows */
308 >                if (rbf_colorimetry == RBCtristimulus) {
309 >                        fputc('\n', cfp[CIE_X]);
310 >                        fputc('\n', cfp[CIE_Z]);
311 >                }
312              }
313 <            data_epilogue();
313 >            if (fclose(cfp[CIE_Y])) {
314 >                fprintf(stderr, "%s: error writing Y output\n", progname);
315 >                exit(1);
316 >            }
317 >            if (rbf_colorimetry == RBCtristimulus &&
318 >                        (fclose(cfp[CIE_X]) || fclose(cfp[CIE_Z]))) {
319 >                fprintf(stderr, "%s: error writing X/Z output\n", progname);
320 >                exit(1);
321 >            }
322          }
323                                                  /* back transmission */
324 <        if (bsd.tb != NULL) {
324 >        if ((bsd.tb != NULL) | (bsd.tf != NULL)) {
325              input_orient = -1; output_orient = 1;
326 <            data_prologue();
326 >            cfp[CIE_Y] = open_component_file(CIE_Y);
327 >            if (bsd.tb != NULL)
328 >                rbf_colorimetry = bsd.tb->comp[0].cspec[2].flags
329 >                                        ? RBCtristimulus : RBCphotopic ;
330 >            if (rbf_colorimetry == RBCtristimulus) {
331 >                cfp[CIE_X] = open_component_file(CIE_X);
332 >                cfp[CIE_Z] = open_component_file(CIE_Z);
333 >            }
334              for (j = 0; j < abp->nangles; j++) {
335                  for (i = 0; i < abp->nangles; i++) {
336 <                    sum = 0;                    /* average over patches */
336 >                    sum = 0;            /* average over patches */
337 >                    xsum = ysum = 0;
338                      for (n = npsamps; n-- > 0; ) {
339                          fo_getvec(vout, j+(n+frandom())/npsamps, abp);
340                          bi_getvec(vin, i+urand(n), abp);
341 <                        ec = SDevalBSDF(&sv, vout, vin, &bsd);
341 >                        ec = SDevalBSDF(&sdv, vin, vout, &bsd);
342                          if (ec != SDEnone)
343                                  goto err;
344 <                        sum += sv.cieY;
344 >                        sum += sdv.cieY;
345 >                        if (rbf_colorimetry == RBCtristimulus) {
346 >                                xsum += sdv.cieY * sdv.spec.cx;
347 >                                ysum += sdv.cieY * sdv.spec.cy;
348 >                        }
349                      }
350 <                    printf("\t%.3e\n", sum/npsamps);
350 >                    fprintf(cfp[CIE_Y], "\t%.3e\n", sum/npsamps);
351 >                    if (rbf_colorimetry == RBCtristimulus) {
352 >                        fprintf(cfp[CIE_X], "\t%.3e\n", xsum*sum/(npsamps*ysum));
353 >                        fprintf(cfp[CIE_Z], "\t%.3e\n",
354 >                                (sum - xsum - ysum)*sum/(npsamps*ysum));
355 >                    }
356                  }
357 <                putchar('\n');                  /* extra space between rows */
357 >                if (rbf_colorimetry == RBCtristimulus) {
358 >                        fputc('\n', cfp[CIE_X]);
359 >                        fputc('\n', cfp[CIE_Z]);
360 >                }
361              }
362 <            data_epilogue();
362 >            if (fclose(cfp[CIE_Y])) {
363 >                fprintf(stderr, "%s: error writing Y output\n", progname);
364 >                exit(1);
365 >            }
366 >            if (rbf_colorimetry == RBCtristimulus &&
367 >                        (fclose(cfp[CIE_X]) || fclose(cfp[CIE_Z]))) {
368 >                fprintf(stderr, "%s: error writing X/Z output\n", progname);
369 >                exit(1);
370 >            }
371          }
372          SDfreeBSDF(&bsd);                       /* all done */
373          return;
# Line 261 | Line 381 | static void
381   eval_function(char *funame)
382   {
383          ANGLE_BASIS     *abp = get_basis(kbasis);
384 +        int             assignD = (fundefined(funame) < 6);
385 +        FILE            *ofp = open_component_file(CIE_Y);
386          double          iovec[6];
387          double          sum;
388          int             i, j, n;
389  
390          initurand(npsamps);
269        data_prologue();                        /* begin output */
391          for (j = 0; j < abp->nangles; j++) {    /* run through directions */
392              for (i = 0; i < abp->nangles; i++) {
393                  sum = 0;
# Line 281 | Line 402 | eval_function(char *funame)
402                      else
403                          bi_getvec(iovec, i+urand(n), abp);
404  
405 +                    if (assignD) {
406 +                        varset("Dx", '=', -iovec[3]);
407 +                        varset("Dy", '=', -iovec[4]);
408 +                        varset("Dz", '=', -iovec[5]);
409 +                        ++eclock;
410 +                    }
411                      sum += funvalue(funame, 6, iovec);
412                  }
413 <                printf("\t%.3e\n", sum/npsamps);
413 >                fprintf(ofp, "\t%.3e\n", sum/npsamps);
414              }
415 <            putchar('\n');
415 >            fputc('\n', ofp);
416 >            prog_show((j+1.)/abp->nangles);
417          }
418 <        data_epilogue();                        /* finish output */
418 >        prog_done();
419 >        if (fclose(ofp)) {
420 >                fprintf(stderr, "%s: error writing Y output\n", progname);
421 >                exit(1);
422 >        }
423   }
424  
425   /* Interpolate and output a radial basis function BSDF representation */
426   static void
427   eval_rbf(void)
428   {
429 <        ANGLE_BASIS     *abp = get_basis(kbasis);
430 <        float           bsdfarr[MAXPATCHES*MAXPATCHES];
431 <        FVECT           vin, vout;
432 <        RBFNODE         *rbf;
433 <        double          sum;
434 <        int             i, j, n;
435 <                                                /* sanity check */
436 <        if (abp->nangles > MAXPATCHES) {
437 <                fprintf(stderr, "%s: too many patches!\n", progname);
438 <                exit(1);
439 <        }
440 <        data_prologue();                        /* begin output */
441 <        for (i = 0; i < abp->nangles; i++) {
442 <            if (input_orient > 0)               /* use incident patch center */
443 <                fi_getvec(vin, i+.5*(i>0), abp);
444 <            else
445 <                bi_getvec(vin, i+.5*(i>0), abp);
429 >    ANGLE_BASIS *abp = get_basis(kbasis);
430 >    float       (*XZarr)[2] = NULL;
431 >    float       bsdfarr[MAXPATCHES*MAXPATCHES];
432 >    FILE        *cfp[3];
433 >    FVECT       vin, vout;
434 >    double      sum, xsum, ysum, normf;
435 >    int         i, j, ni, no, nisamps, nosamps;
436 >                                        /* sanity check */
437 >    if (abp->nangles > MAXPATCHES) {
438 >        fprintf(stderr, "%s: too many patches!\n", progname);
439 >        exit(1);
440 >    }
441 >    memset(bsdfarr, 0, sizeof(bsdfarr));
442 >    if (rbf_colorimetry == RBCtristimulus)
443 >        XZarr = (float (*)[2])calloc(abp->nangles*abp->nangles, 2*sizeof(float));
444 >    nosamps = (int)(pow((double)npsamps, 0.67) + .5);
445 >    nisamps = (npsamps + (nosamps>>1)) / nosamps;
446 >    normf = 1./(double)(nisamps*nosamps);
447 >    for (i = 0; i < abp->nangles; i++) {
448 >        for (ni = nisamps; ni--; ) {            /* sample over incident patch */
449 >            RBFNODE     *rbf;
450 >            if (input_orient > 0)               /* vary incident patch loc. */
451 >                fi_getvec(vin, i+urand(ni), abp);
452 >            else
453 >                bi_getvec(vin, i+urand(ni), abp);
454  
455 <            rbf = advect_rbf(vin);              /* compute radial basis func */
455 >            rbf = advect_rbf(vin, lobe_lim);    /* compute radial basis func */
456  
457 <            for (j = 0; j < abp->nangles; j++) {
458 <                sum = 0;                        /* sample over exiting patch */
459 <                for (n = npsamps; n--; ) {
457 >            for (j = 0; j < abp->nangles; j++) {
458 >                sum = 0;                        /* sample over exiting patch */
459 >                xsum = ysum = 0;
460 >                for (no = nosamps; no--; ) {
461 >                    SDValue     sdv;
462                      if (output_orient > 0)
463 <                        fo_getvec(vout, j+(n+frandom())/npsamps, abp);
463 >                        fo_getvec(vout, j+(no+frandom())/nosamps, abp);
464                      else
465 <                        bo_getvec(vout, j+(n+frandom())/npsamps, abp);
465 >                        bo_getvec(vout, j+(no+frandom())/nosamps, abp);
466  
467 <                    sum += eval_rbfrep(rbf, vout) / vout[2];
467 >                    eval_rbfcol(&sdv, rbf, vout);
468 >                    sum += sdv.cieY;
469 >                    if (rbf_colorimetry == RBCtristimulus) {
470 >                        xsum += sdv.cieY * sdv.spec.cx;
471 >                        ysum += sdv.cieY * sdv.spec.cy;
472 >                    }
473                  }
474 <                bsdfarr[j*abp->nangles + i] = sum*output_orient/npsamps;
474 >                no = j*abp->nangles + i;
475 >                bsdfarr[no] += sum * normf;
476 >                if (rbf_colorimetry == RBCtristimulus) {
477 >                    XZarr[no][0] += xsum*sum*normf/ysum;
478 >                    XZarr[no][1] += (sum - xsum - ysum)*sum*normf/ysum;
479 >                }
480              }
481 <            if (rbf != NULL)
481 >            if (rbf != NULL)
482                  free(rbf);
483          }
484 <        n = 0;                                  /* write out our matrix */
485 <        for (j = 0; j < abp->nangles; j++) {
486 <            for (i = 0; i < abp->nangles; i++)
487 <                printf("\t%.3e\n", bsdfarr[n++]);
488 <            putchar('\n');
484 >        prog_show((i+1.)/abp->nangles);
485 >    }
486 >                                        /* write out our matrix */
487 >    cfp[CIE_Y] = open_component_file(CIE_Y);
488 >    no = 0;
489 >    for (j = 0; j < abp->nangles; j++) {
490 >        for (i = 0; i < abp->nangles; i++, no++)
491 >        fprintf(cfp[CIE_Y], "\t%.3e\n", bsdfarr[no]);
492 >        fputc('\n', cfp[CIE_Y]);
493 >    }
494 >    prog_done();
495 >    if (fclose(cfp[CIE_Y])) {
496 >        fprintf(stderr, "%s: error writing Y output\n", progname);
497 >        exit(1);
498 >    }
499 >    if (XZarr == NULL)                  /* no color? */
500 >        return;
501 >    cfp[CIE_X] = open_component_file(CIE_X);
502 >    cfp[CIE_Z] = open_component_file(CIE_Z);
503 >    no = 0;
504 >    for (j = 0; j < abp->nangles; j++) {
505 >        for (i = 0; i < abp->nangles; i++, no++) {
506 >        fprintf(cfp[CIE_X], "\t%.3e\n", XZarr[no][0]);
507 >        fprintf(cfp[CIE_Z], "\t%.3e\n", XZarr[no][1]);
508 >        }
509 >        fputc('\n', cfp[CIE_X]);
510 >        fputc('\n', cfp[CIE_Z]);
511 >    }
512 >    free(XZarr);
513 >    if (fclose(cfp[CIE_X]) || fclose(cfp[CIE_Z])) {
514 >        fprintf(stderr, "%s: error writing X/Z output\n", progname);
515 >        exit(1);
516 >    }
517 > }
518 >
519 > #if defined(_WIN32) || defined(_WIN64)
520 > /* Execute wrapBSDF command (may never return) */
521 > static int
522 > wrap_up(void)
523 > {
524 >        char    cmd[32700];
525 >
526 >        if (bsdf_manuf[0]) {
527 >                add_wbsdf("-f", 1);
528 >                strcpy(cmd, "m=");
529 >                strcpy(cmd+2, bsdf_manuf);
530 >                add_wbsdf(cmd, 0);
531          }
532 <        data_epilogue();                        /* finish output */
532 >        if (bsdf_name[0]) {
533 >                add_wbsdf("-f", 1);
534 >                strcpy(cmd, "n=");
535 >                strcpy(cmd+2, bsdf_name);
536 >                add_wbsdf(cmd, 0);
537 >        }
538 >        if (!convert_commandline(cmd, sizeof(cmd), wrapBSDF)) {
539 >                fputs(progname, stderr);
540 >                fputs(": command line too long in wrap_up()\n", stderr);
541 >                return(1);
542 >        }
543 >        return(system(cmd));
544   }
545 + #else
546 + /* Execute wrapBSDF command (may never return) */
547 + static int
548 + wrap_up(void)
549 + {
550 +        char    buf[256];
551 +        char    *compath = getpath((char *)wrapBSDF[0], getenv("PATH"), X_OK);
552  
553 +        if (compath == NULL) {
554 +                fprintf(stderr, "%s: cannot locate %s\n", progname, wrapBSDF[0]);
555 +                return(1);
556 +        }
557 +        if (bsdf_manuf[0]) {
558 +                add_wbsdf("-f", 1);
559 +                strcpy(buf, "m=");
560 +                strcpy(buf+2, bsdf_manuf);
561 +                add_wbsdf(buf, 0);
562 +        }
563 +        if (bsdf_name[0]) {
564 +                add_wbsdf("-f", 1);
565 +                strcpy(buf, "n=");
566 +                strcpy(buf+2, bsdf_name);
567 +                add_wbsdf(buf, 0);
568 +        }
569 +        execv(compath, wrapBSDF);       /* successful call never returns */
570 +        perror(compath);
571 +        return(1);
572 + }
573 + #endif
574 +
575 + #define HEAD_BUFLEN     10240
576 + static char     head_buf[HEAD_BUFLEN];
577 + static int      cur_headlen = 0;
578 +
579 + /* Record header line as comment associated with this SIR input */
580 + static int
581 + record2header(char *s)
582 + {
583 +        int     len = strlen(s);
584 +
585 +        if (cur_headlen+len >= HEAD_BUFLEN-6)
586 +                return(0);
587 +                                        /* includes EOL */
588 +        strcpy(head_buf+cur_headlen, s);
589 +        cur_headlen += len;
590 +
591 + #if defined(_WIN32) || defined(_WIN64)
592 +        if (head_buf[cur_headlen-1] == '\n')
593 +                head_buf[cur_headlen-1] = '\t';
594 + #endif
595 +        return(1);
596 + }
597 +
598 + /* Finish off header for this file */
599 + static void
600 + done_header(void)
601 + {
602 +        while (cur_headlen > 0 && isspace(head_buf[cur_headlen-1]))
603 +                --cur_headlen;
604 +        head_buf[cur_headlen] = '\0';
605 +        if (!cur_headlen)
606 +                return;
607 +        add_wbsdf("-C", 1);
608 +        add_wbsdf(head_buf, 0);
609 +        head_buf[cur_headlen=0] = '\0';
610 + }
611 +
612   /* Read in BSDF and interpolate as Klems matrix representation */
613   int
614   main(int argc, char *argv[])
615   {
616          int     dofwd = 0, dobwd = 1;
617 +        char    buf[1024];
618          char    *cp;
619          int     i, na;
620 <
621 <        progname = argv[0];
620 >                                                /* set global progname */
621 >        fixargv0(argv[0]);
622          esupport |= E_VARIABLE|E_FUNCTION|E_RCONST;
623          esupport &= ~(E_INCHAN|E_OUTCHAN);
624          scompile("PI:3.14159265358979323846", NULL, 0);
# Line 363 | Line 635 | main(int argc, char *argv[])
635                          single_plane_incident = 0;
636                          break;
637                  case 'f':
638 <                        if (!argv[i][2]) {
639 <                                fcompile(argv[++i]);
640 <                                single_plane_incident = 0;
638 >                        if ((argv[i][0] == '-') & !argv[i][2]) {
639 >                                if (strchr(argv[++i], '=') != NULL) {
640 >                                        add_wbsdf("-f", 1);
641 >                                        add_wbsdf(argv[i], 1);
642 >                                } else {
643 >                                        char    *fpath = getpath(argv[i],
644 >                                                            getrlibpath(), 0);
645 >                                        if (fpath == NULL) {
646 >                                                fprintf(stderr,
647 >                                                "%s: cannot find file '%s'\n",
648 >                                                        argv[0], argv[i]);
649 >                                                return(1);
650 >                                        }
651 >                                        fcompile(fpath);
652 >                                        single_plane_incident = 0;
653 >                                }
654                          } else
655                                  dofwd = (argv[i][0] == '+');
656                          break;
# Line 373 | Line 658 | main(int argc, char *argv[])
658                          dobwd = (argv[i][0] == '+');
659                          break;
660                  case 'h':
661 <                        kbasis = "LBNL/Klems Half";
661 >                        kbasis = klems_half;
662 >                        add_wbsdf("-a", 1);
663 >                        add_wbsdf("kh", 1);
664                          break;
665                  case 'q':
666 <                        kbasis = "LBNL/Klems Quarter";
666 >                        kbasis = klems_quarter;
667 >                        add_wbsdf("-a", 1);
668 >                        add_wbsdf("kq", 1);
669                          break;
670 +                case 'l':
671 +                        lobe_lim = atoi(argv[++i]);
672 +                        break;
673 +                case 'p':
674 +                        do_prog = atoi(argv[i]+2);
675 +                        break;
676 +                case 'C':
677 +                        add_wbsdf(argv[i], 1);
678 +                        add_wbsdf(argv[++i], 1);
679 +                        break;
680                  default:
681                          goto userr;
682                  }
683 +        if (kbasis == klems_full) {             /* default (full) basis? */
684 +                add_wbsdf("-a", 1);
685 +                add_wbsdf("kf", 1);
686 +        }
687 +        strcpy(buf, "File produced by: ");
688 +        if (convert_commandline(buf+18, sizeof(buf)-18, argv) != NULL) {
689 +                add_wbsdf("-C", 1); add_wbsdf(buf, 0);
690 +        }
691          if (single_plane_incident >= 0) {       /* function-based BSDF? */
692 <                if (i != argc-1 || fundefined(argv[i]) != 6) {
692 >                if (i != argc-1 || fundefined(argv[i]) < 3) {
693                          fprintf(stderr,
694          "%s: need single function with 6 arguments: bsdf(ix,iy,iz,ox,oy,oz)\n",
695                                          progname);
696 +                        fprintf(stderr, "\tor 3 arguments using Dx,Dy,Dz: bsdf(ix,iy,iz)\n");
697                          goto userr;
698                  }
699 <                xml_header(argc, argv);                 /* start XML output */
700 <                xml_prologue(NULL);
699 >                doptimize(1);                   /* optimize definitions */
700 >                ++eclock;
701                  if (dofwd) {
702                          input_orient = -1;
703                          output_orient = -1;
704 <                        eval_function(argv[i]);         /* outside reflectance */
704 >                        prog_start("Evaluating outside reflectance");
705 >                        eval_function(argv[i]);
706                          output_orient = 1;
707 <                        eval_function(argv[i]);         /* outside -> inside */
707 >                        prog_start("Evaluating outside->inside transmission");
708 >                        eval_function(argv[i]);
709                  }
710                  if (dobwd) {
711                          input_orient = 1;
712                          output_orient = 1;
713 <                        eval_function(argv[i]);         /* inside reflectance */
713 >                        prog_start("Evaluating inside reflectance");
714 >                        eval_function(argv[i]);
715                          output_orient = -1;
716 <                        eval_function(argv[i]);         /* inside -> outside */
716 >                        prog_start("Evaluating inside->outside transmission");
717 >                        eval_function(argv[i]);
718                  }
719 <                xml_epilogue();                 /* finish XML output & exit */
408 <                return(0);
719 >                return(wrap_up());
720          }
721                                                  /* XML input? */
722          if (i == argc-1 && (cp = argv[i]+strlen(argv[i])-4) > argv[i] &&
723                                  !strcasecmp(cp, ".xml")) {
413                xml_header(argc, argv);         /* start XML output */
724                  eval_bsdf(argv[i]);             /* load & resample BSDF */
725 <                xml_epilogue();                 /* finish XML output & exit */
416 <                return(0);
725 >                return(wrap_up());
726          }
727          if (i < argc) {                         /* open input files if given */
728                  int     nbsdf = 0;
# Line 424 | Line 733 | main(int argc, char *argv[])
733                                                  progname, argv[i]);
734                                  return(1);
735                          }
736 +                        sprintf(buf, "%s:\n", argv[i]);
737 +                        record2header(buf);
738 +                        sir_headshare = &record2header;
739                          if (!load_bsdf_rep(fpin))
740                                  return(1);
741                          fclose(fpin);
742 <                        if (!nbsdf++) {         /* start XML on first dist. */
743 <                                xml_header(argc, argv);
744 <                                xml_prologue(NULL);
433 <                        }
742 >                        done_header();
743 >                        sprintf(buf, "Interpolating component '%s'", argv[i]);
744 >                        prog_start(buf);
745                          eval_rbf();
746                  }
747 <                xml_epilogue();                 /* finish XML output & exit */
437 <                return(0);
747 >                return(wrap_up());
748          }
749          SET_FILE_BINARY(stdin);                 /* load from stdin */
750 +        record2header("<stdin>:\n");
751 +        sir_headshare = &record2header;
752          if (!load_bsdf_rep(stdin))
753                  return(1);
754 <        xml_header(argc, argv);                 /* start XML output */
755 <        xml_prologue(NULL);
754 >        done_header();
755 >        prog_start("Interpolating from standard input");
756          eval_rbf();                             /* resample dist. */
757 <        xml_epilogue();                         /* finish XML output & exit */
446 <        return(0);
757 >        return(wrap_up());
758   userr:
759          fprintf(stderr,
760 <        "Usage: %s [-n spp][-h|-q][bsdf.sir ..] > bsdf.xml\n", progname);
760 >        "Usage: %s [-n spp][-h|-q][-l maxlobes] [bsdf.sir ..] > bsdf.xml\n", progname);
761          fprintf(stderr,
762          "   or: %s [-n spp][-h|-q] bsdf_in.xml > bsdf_out.xml\n", progname);
763          fprintf(stderr,

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines