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

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines