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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines