ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdf2klems.c
Revision: 2.2
Committed: Tue Apr 23 04:40:23 2013 UTC (11 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.1: +47 -11 lines
Log Message:
Untested implementation of Klems conversion from radial basis function

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.2 static const char RCSid[] = "$Id: bsdf2klems.c,v 2.1 2013/04/21 23:01:14 greg Exp $";
3 greg 2.1 #endif
4     /*
5     * Load measured BSDF interpolant and write out as XML file with Klems matrix.
6     *
7     * G. Ward
8     */
9    
10     #define _USE_MATH_DEFINES
11     #include <stdio.h>
12     #include <stdlib.h>
13     #include <string.h>
14     #include <math.h>
15     #include "random.h"
16     #include "platform.h"
17     #include "calcomp.h"
18     #include "bsdfrep.h"
19     #include "bsdf_m.h"
20     /* global argv[0] */
21     char *progname;
22     /* selected basis function name */
23     static const char *kbasis = "LBNL/Klems Full";
24     /* number of BSDF samples per patch */
25     static int npsamps = 256;
26    
27     /* Return angle basis corresponding to the given name */
28     ANGLE_BASIS *
29     get_basis(const char *bn)
30     {
31     int n = nabases;
32    
33     while (n-- > 0)
34     if (!strcasecmp(bn, abase_list[n].name))
35     return &abase_list[n];
36     return NULL;
37     }
38    
39     /* Output XML prologue to stdout */
40     static void
41     xml_prologue(int ac, char *av[])
42     {
43     ANGLE_BASIS *abp = get_basis(kbasis);
44     int i;
45    
46     if (abp == NULL) {
47     fprintf(stderr, "%s: unknown angle basis '%s'\n", progname, kbasis);
48     exit(1);
49     }
50     puts("<?xml version=\"1.0\" encoding=\"UTF-8\"?>");
51     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\">");
52     fputs("<!-- File produced by:", stdout);
53     while (ac-- > 0) {
54     fputc(' ', stdout);
55     fputs(*av++, stdout);
56     }
57     puts(" -->");
58     puts("<WindowElementType>System</WindowElementType>");
59     puts("<FileType>BSDF</FileType>");
60     puts("<Optical>");
61     puts("<Layer>");
62     puts("\t<Material>");
63     puts("\t\t<Name>Name</Name>");
64     puts("\t\t<Manufacturer>Manufacturer</Manufacturer>");
65     puts("\t\t<DeviceType>Other</DeviceType>");
66     puts("\t</Material>");
67     puts("\t<DataDefinition>");
68     puts("\t\t<IncidentDataStructure>Columns</IncidentDataStructure>");
69     puts("\t\t<AngleBasis>");
70     printf("\t\t\t<AngleBasisName>%s</AngleBasisName>\n", kbasis);
71     for (i = 0; abp->lat[i].nphis; i++) {
72     puts("\t\t\t<AngleBasisBlock>");
73     printf("\t\t\t<Theta>%g</Theta>\n", i ?
74     .5*(abp->lat[i].tmin + abp->lat[i+1].tmin) :
75     .0 );
76     printf("\t\t\t<nPhis>%d</nPhis>", abp->lat[i].nphis);
77     puts("\t\t\t<ThetaBounds>");
78     printf("\t\t\t\t<LowerTheta>%g</LowerTheta>\n", abp->lat[i].tmin);
79     printf("\t\t\t\t<UpperTheta>%g</UpperTheta>\n", abp->lat[i+1].tmin);
80     puts("\t\t\t</ThetaBounds>");
81     puts("\t\t\t</AngleBasisBlock>");
82     }
83     puts("\t\t</AngleBasis>");
84     puts("\t</DataDefinition>");
85     }
86    
87     /* Output XML data prologue to stdout */
88     static void
89     data_prologue()
90     {
91     static const char *bsdf_type[4] = {
92     "Reflection Front",
93     "Transmission Front",
94     "Transmission Back",
95     "Reflection Back"
96     };
97    
98     puts("\t<WavelengthData>");
99     puts("\t\t<LayerNumber>System</LayerNumber>");
100     puts("\t\t<Wavelength unit=\"Integral\">Visible</Wavelength>");
101     puts("\t\t<SourceSpectrum>CIE Illuminant D65 1nm.ssp</SourceSpectrum>");
102     puts("\t\t<DetectorSpectrum>ASTM E308 1931 Y.dsp</DetectorSpectrum>");
103     puts("\t\t<WavelengthDataBlock>");
104     printf("\t\t\t<WavelengthDataDirection>%s</WavelengthDataDirection>\n",
105     bsdf_type[(input_orient>0)<<1 | (output_orient>0)]);
106     printf("\t\t\t<ColumnAngleBasis>%s</ColumnAngleBasis>\n", kbasis);
107     printf("\t\t\t<RowAngleBasis>%s</RowAngleBasis>\n", kbasis);
108     puts("\t\t\t<ScatteringDataType>BTDF</ScatteringDataType>");
109     puts("\t\t\t<ScatteringData>");
110     }
111    
112     /* Output XML data epilogue to stdout */
113     static void
114     data_epilogue(void)
115     {
116     puts("\t\t\t</ScatteringData>");
117     puts("\t\t</WavelengthDataBlock>");
118     puts("\t</WavelengthData>");
119     }
120    
121     /* Output XML epilogue to stdout */
122     static void
123     xml_epilogue(void)
124     {
125     puts("</Layer>");
126     puts("</Optical>");
127     puts("</WindowElement>");
128     }
129    
130 greg 2.2 /* Load and resample XML BSDF description using Klems basis */
131 greg 2.1 static void
132     eval_bsdf(const char *fname)
133     {
134     ANGLE_BASIS *abp = get_basis(kbasis);
135     SDData bsd;
136     SDError ec;
137     FVECT vin, vout;
138     SDValue sv;
139     double sum;
140     int i, j, n;
141    
142     SDclearBSDF(&bsd, fname); /* load BSDF file */
143     if ((ec = SDloadFile(&bsd, fname)) != SDEnone)
144     goto err;
145     /* front reflection */
146     if (bsd.rf != NULL || bsd.rLambFront.cieY > .002) {
147     input_orient = 1; output_orient = 1;
148     data_prologue();
149     for (j = 0; j < abp->nangles; j++) {
150     for (i = 0; i < abp->nangles; i++) {
151     sum = 0; /* average over patches */
152     for (n = npsamps; n-- > 0; ) {
153     fo_getvec(vout, j+(n+frandom())/npsamps, abp);
154     fi_getvec(vin, i+(n+frandom())/npsamps, abp);
155     ec = SDevalBSDF(&sv, vout, vin, &bsd);
156     if (ec != SDEnone)
157     goto err;
158     sum += sv.cieY;
159     }
160     printf("\t%.3e\n", sum/npsamps);
161     }
162     putchar('\n'); /* extra space between rows */
163     }
164     data_epilogue();
165     }
166     /* back reflection */
167     if (bsd.rb != NULL || bsd.rLambBack.cieY > .002) {
168     input_orient = -1; output_orient = -1;
169     data_prologue();
170     for (j = 0; j < abp->nangles; j++) {
171     for (i = 0; i < abp->nangles; i++) {
172     sum = 0; /* average over patches */
173     for (n = npsamps; n-- > 0; ) {
174     bo_getvec(vout, j+(n+frandom())/npsamps, abp);
175     bi_getvec(vin, i+(n+frandom())/npsamps, abp);
176     ec = SDevalBSDF(&sv, vout, vin, &bsd);
177     if (ec != SDEnone)
178     goto err;
179     sum += sv.cieY;
180     }
181     printf("\t%.3e\n", sum/npsamps);
182     }
183     putchar('\n'); /* extra space between rows */
184     }
185     data_epilogue();
186     }
187     /* front transmission */
188     if (bsd.tf != NULL || bsd.tLamb.cieY > .002) {
189     input_orient = 1; output_orient = -1;
190     data_prologue();
191     for (j = 0; j < abp->nangles; j++) {
192     for (i = 0; i < abp->nangles; i++) {
193     sum = 0; /* average over patches */
194     for (n = npsamps; n-- > 0; ) {
195     bo_getvec(vout, j+(n+frandom())/npsamps, abp);
196     fi_getvec(vin, i+(n+frandom())/npsamps, abp);
197     ec = SDevalBSDF(&sv, vout, vin, &bsd);
198     if (ec != SDEnone)
199     goto err;
200     sum += sv.cieY;
201     }
202     printf("\t%.3e\n", sum/npsamps);
203     }
204     putchar('\n'); /* extra space between rows */
205     }
206     data_epilogue();
207     }
208     /* back transmission */
209     if (bsd.tb != NULL) {
210     input_orient = -1; output_orient = 1;
211     data_prologue();
212     for (j = 0; j < abp->nangles; j++) {
213     for (i = 0; i < abp->nangles; i++) {
214     sum = 0; /* average over patches */
215     for (n = npsamps; n-- > 0; ) {
216     fo_getvec(vout, j+(n+frandom())/npsamps, abp);
217     bi_getvec(vin, i+(n+frandom())/npsamps, abp);
218     ec = SDevalBSDF(&sv, vout, vin, &bsd);
219     if (ec != SDEnone)
220     goto err;
221     sum += sv.cieY;
222     }
223     printf("\t%.3e\n", sum/npsamps);
224     }
225     putchar('\n'); /* extra space between rows */
226     }
227     data_epilogue();
228     }
229     SDfreeBSDF(&bsd); /* all done */
230     return;
231     err:
232     SDreportError(ec, stderr);
233     exit(1);
234     }
235    
236 greg 2.2 /* Interpolate and output a BSDF function using Klems basis */
237 greg 2.1 static void
238     eval_function(char *funame)
239     {
240     ANGLE_BASIS *abp = get_basis(kbasis);
241     double iovec[6];
242     double sum;
243     int i, j, n;
244    
245     data_prologue(); /* begin output */
246     for (j = 0; j < abp->nangles; j++) { /* run through directions */
247     for (i = 0; i < abp->nangles; i++) {
248     sum = 0;
249     for (n = npsamps; n--; ) { /* average over patches */
250     if (output_orient > 0)
251     fo_getvec(iovec+3, j+(n+frandom())/npsamps, abp);
252     else
253     bo_getvec(iovec+3, j+(n+frandom())/npsamps, abp);
254    
255     if (input_orient > 0)
256 greg 2.2 fi_getvec(iovec, i+(n+frandom())/npsamps, abp);
257 greg 2.1 else
258 greg 2.2 bi_getvec(iovec, i+(n+frandom())/npsamps, abp);
259 greg 2.1
260     sum += funvalue(funame, 6, iovec);
261     }
262     printf("\t%.3e\n", sum/npsamps);
263     }
264     putchar('\n');
265     }
266     data_epilogue(); /* finish output */
267     }
268    
269     /* Interpolate and output a radial basis function BSDF representation */
270     static void
271     eval_rbf(void)
272     {
273 greg 2.2 #define MAXPATCHES 145
274 greg 2.1 ANGLE_BASIS *abp = get_basis(kbasis);
275 greg 2.2 float bsdfarr[MAXPATCHES*MAXPATCHES];
276     FVECT vin, vout;
277     RBFNODE *rbf;
278 greg 2.1 double sum;
279     int i, j, n;
280 greg 2.2 /* sanity check */
281     if (abp->nangles > MAXPATCHES) {
282     fprintf(stderr, "%s: too many patches!\n", progname);
283     exit(1);
284     }
285     data_prologue(); /* begin output */
286     for (i = 0; i < abp->nangles; i++) {
287     if (input_orient > 0) /* use incident patch center */
288     fi_getvec(vin, i+.5*(i>0), abp);
289     else
290     bi_getvec(vin, i+.5*(i>0), abp);
291    
292     rbf = advect_rbf(vin); /* compute radial basis func */
293    
294     for (j = 0; j < abp->nangles; j++) {
295     sum = 0; /* sample over exiting patch */
296     for (n = npsamps; n--; ) {
297     if (output_orient > 0)
298     fo_getvec(vout, j+(n+frandom())/npsamps, abp);
299     else
300     bo_getvec(vout, j+(n+frandom())/npsamps, abp);
301 greg 2.1
302 greg 2.2 sum += eval_rbfrep(rbf, vout) / vout[2];
303     }
304     bsdfarr[j*abp->nangles + i] = sum*output_orient/npsamps;
305     }
306     }
307     n = 0; /* write out our matrix */
308     for (j = 0; j < abp->nangles; j++) {
309     for (i = 0; i < abp->nangles; i++)
310     printf("\t%.3e\n", bsdfarr[n++]);
311     putchar('\n');
312     }
313     data_epilogue(); /* finish output */
314     #undef MAXPATCHES
315 greg 2.1 }
316    
317     /* Read in BSDF and interpolate as Klems matrix representation */
318     int
319     main(int argc, char *argv[])
320     {
321     int dofwd = 0, dobwd = 1;
322     char *cp;
323     int i, na;
324    
325     progname = argv[0];
326     esupport |= E_VARIABLE|E_FUNCTION|E_RCONST;
327     esupport &= ~(E_INCHAN|E_OUTCHAN);
328     scompile("PI:3.14159265358979323846", NULL, 0);
329     biggerlib();
330 greg 2.2 for (i = 1; i < argc && (argv[i][0] == '-') | (argv[i][0] == '+'); i++)
331 greg 2.1 switch (argv[i][1]) { /* get options */
332     case 'n':
333     npsamps = atoi(argv[++i]);
334     if (npsamps <= 0)
335     goto userr;
336     break;
337     case 'e':
338     scompile(argv[++i], NULL, 0);
339     single_plane_incident = 0;
340     break;
341     case 'f':
342     if (!argv[i][2]) {
343     fcompile(argv[++i]);
344     single_plane_incident = 0;
345     } else
346     dofwd = (argv[i][0] == '+');
347     break;
348     case 'b':
349     dobwd = (argv[i][0] == '+');
350     break;
351     case 'h':
352     kbasis = "LBNL/Klems Half";
353     break;
354     case 'q':
355     kbasis = "LBNL/Klems Quarter";
356     break;
357     default:
358     goto userr;
359     }
360     if (single_plane_incident >= 0) { /* function-based BSDF? */
361     if (i != argc-1 || fundefined(argv[i]) != 6) {
362     fprintf(stderr,
363     "%s: need single function with 6 arguments: bsdf(ix,iy,iz,ox,oy,oz)\n",
364     progname);
365     goto userr;
366     }
367     xml_prologue(argc, argv); /* start XML output */
368     if (dofwd) {
369     input_orient = -1;
370     output_orient = -1;
371     eval_function(argv[i]); /* outside reflectance */
372     output_orient = 1;
373     eval_function(argv[i]); /* outside -> inside */
374     }
375     if (dobwd) {
376     input_orient = 1;
377     output_orient = 1;
378     eval_function(argv[i]); /* inside reflectance */
379     output_orient = -1;
380     eval_function(argv[i]); /* inside -> outside */
381     }
382     xml_epilogue(); /* finish XML output & exit */
383     return(0);
384     }
385 greg 2.2 /* XML input? */
386     if (i == argc-1 && (cp = argv[i]+strlen(argv[i])-4) > argv[i] &&
387     !strcasecmp(cp, ".xml")) {
388 greg 2.1 xml_prologue(argc, argv); /* start XML output */
389     eval_bsdf(argv[i]); /* load & resample BSDF */
390     xml_epilogue(); /* finish XML output & exit */
391     return(0);
392     }
393     if (i < argc) { /* open input files if given */
394     int nbsdf = 0;
395     for ( ; i < argc; i++) { /* interpolate each component */
396     FILE *fpin = fopen(argv[i], "rb");
397     if (fpin == NULL) {
398     fprintf(stderr, "%s: cannot open BSDF interpolant '%s'\n",
399     progname, argv[i]);
400     return(1);
401     }
402     if (!load_bsdf_rep(fpin))
403     return(1);
404     fclose(fpin);
405     if (!nbsdf++) /* start XML on first dist. */
406     xml_prologue(argc, argv);
407     eval_rbf();
408     }
409     xml_epilogue(); /* finish XML output & exit */
410     return(0);
411     }
412     SET_FILE_BINARY(stdin); /* load from stdin */
413     if (!load_bsdf_rep(stdin))
414     return(1);
415     xml_prologue(argc, argv); /* start XML output */
416     eval_rbf(); /* resample dist. */
417     xml_epilogue(); /* finish XML output & exit */
418     return(0);
419     userr:
420     fprintf(stderr,
421     "Usage: %s [-n spp][-h|-q][bsdf.sir ..] > bsdf.xml\n",
422     progname);
423     fprintf(stderr,
424     " or: %s [-n spp][-h|-q] bsdf_in.xml > bsdf_out.xml\n",
425     progname);
426     fprintf(stderr,
427     " or: %s [-n spp][-h|-q][{+|-}for[ward]][{+|-}b[ackward]][-e expr][-f file] bsdf_func > bsdf.xml\n",
428     progname);
429     return(1);
430     }