src/cv/bsdf2klems.c

#ifndef lint
static const char RCSid[] = "$Id: bsdf2klems.c,v 2.18 2015/05/27 11:39:37 greg Exp $";
#endif
/*
 * Load measured BSDF interpolant and write out as XML file with Klems matrix.
 *
 *      G. Ward
 */

#define _USE_MATH_DEFINES
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "random.h"
#include "platform.h"
#include "calcomp.h"
#include "bsdfrep.h"
#include "bsdf_m.h"
                                /* assumed maximum # Klems patches */
#define MAXPATCHES      145
                                /* global argv[0] */
char                    *progname;
                                /* selected basis function name */
static const char       *kbasis = "LBNL/Klems Full";
                                /* number of BSDF samples per patch */
static int              npsamps = 256;
                                /* limit on number of RBF lobes */
static int              lobe_lim = 15000;
                                /* progress bar length */
static int              do_prog = 79;


/* Start new progress bar */
#define prog_start(s)   if (do_prog) fprintf(stderr, "%s: %s...\n", progname, s); else

/* Draw progress bar of the appropriate length */
static void
prog_show(double frac)
{
        static unsigned call_cnt = 0;
        static char     lastc[] = "-\\|/";
        char            pbar[256];
        int             nchars;

        if (do_prog <= 1) return;
        if (do_prog > sizeof(pbar)-2)
                do_prog = sizeof(pbar)-2;
        if (frac < 0) frac = 0;
        else if (frac >= 1) frac = .9999;
        nchars = do_prog*frac;
        pbar[0] = '\r';
        memset(pbar+1, '*', nchars);
        pbar[nchars+1] = lastc[call_cnt++ & 3];
        memset(pbar+2+nchars, '-', do_prog-nchars-1);
        pbar[do_prog+1] = '\0';
        fputs(pbar, stderr);
}

/* Finish progress bar */
static void
prog_done(void)
{
        int     n = do_prog;

        if (n <= 1) return;
        fputc('\r', stderr);
        while (n--)
                fputc(' ', stderr);
        fputc('\r', stderr);
}

/* Return angle basis corresponding to the given name */
static ANGLE_BASIS *
get_basis(const char *bn)
{
        int     n = nabases;
        
        while (n-- > 0)
                if (!strcasecmp(bn, abase_list[n].name))
                        return &abase_list[n];
        return NULL;
}

/* Output XML header to stdout */
static void
xml_header(int ac, char *av[])
{
        puts("<?xml version=\"1.0\" encoding=\"UTF-8\"?>");
        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\">");
        fputs("<!-- File produced by:", stdout);
        while (ac-- > 0) {
                fputc(' ', stdout);
                fputs(*av++, stdout);
        }
        puts(" -->");
}

/* Output XML prologue to stdout */
static void
xml_prologue(const SDData *sd)
{
        const char      *matn = (sd && sd->matn[0]) ? sd->matn :
                                bsdf_name[0] ? bsdf_name : "Unknown";
        const char      *makr = (sd && sd->makr[0]) ? sd->makr :
                                bsdf_manuf[0] ? bsdf_manuf : "Unknown";
        ANGLE_BASIS     *abp = get_basis(kbasis);
        int             i;

        if (abp == NULL) {
                fprintf(stderr, "%s: unknown angle basis '%s'\n", progname, kbasis);
                exit(1);
        }
        puts("<WindowElementType>System</WindowElementType>");
        puts("<FileType>BSDF</FileType>");
        puts("<Optical>");
        puts("<Layer>");
        puts("\t<Material>");
        printf("\t\t<Name>%s</Name>\n", matn);
        printf("\t\t<Manufacturer>%s</Manufacturer>\n", makr);
        if (sd && sd->dim[2] > .001) {
                printf("\t\t<Thickness unit=\"meter\">%.3f</Thickness>\n", sd->dim[2]);
                printf("\t\t<Width unit=\"meter\">%.3f</Width>\n", sd->dim[0]);
                printf("\t\t<Height unit=\"meter\">%.3f</Height>\n", sd->dim[1]);
        }
        puts("\t\t<DeviceType>Other</DeviceType>");
        puts("\t</Material>");
        if (sd && sd->mgf != NULL) {
                puts("\t<Geometry format=\"MGF\">");
                puts("\t\t<MGFblock unit=\"meter\">");
                fputs(sd->mgf, stdout);
                puts("</MGFblock>");
                puts("\t</Geometry>");
        }
        puts("\t<DataDefinition>");
        puts("\t\t<IncidentDataStructure>Columns</IncidentDataStructure>");
        puts("\t\t<AngleBasis>");
        printf("\t\t\t<AngleBasisName>%s</AngleBasisName>\n", kbasis);
        for (i = 0; abp->lat[i].nphis; i++) {
                puts("\t\t\t<AngleBasisBlock>");
                printf("\t\t\t<Theta>%g</Theta>\n", i ?
                                .5*(abp->lat[i].tmin + abp->lat[i+1].tmin) :
                                .0 );
                printf("\t\t\t<nPhis>%d</nPhis>\n", abp->lat[i].nphis);
                puts("\t\t\t<ThetaBounds>");
                printf("\t\t\t\t<LowerTheta>%g</LowerTheta>\n", abp->lat[i].tmin);
                printf("\t\t\t\t<UpperTheta>%g</UpperTheta>\n", abp->lat[i+1].tmin);
                puts("\t\t\t</ThetaBounds>");
                puts("\t\t\t</AngleBasisBlock>");
        }
        puts("\t\t</AngleBasis>");
        puts("\t</DataDefinition>");
}

/* Output XML data prologue to stdout */
static void
data_prologue()
{
        static const char       *bsdf_type[4] = {
                                        "Reflection Front",
                                        "Transmission Front",
                                        "Transmission Back",
                                        "Reflection Back"
                                };

        puts("\t<WavelengthData>");
        puts("\t\t<LayerNumber>System</LayerNumber>");
        puts("\t\t<Wavelength unit=\"Integral\">Visible</Wavelength>");
        puts("\t\t<SourceSpectrum>CIE Illuminant D65 1nm.ssp</SourceSpectrum>");
        puts("\t\t<DetectorSpectrum>ASTM E308 1931 Y.dsp</DetectorSpectrum>");
        puts("\t\t<WavelengthDataBlock>");
        printf("\t\t\t<WavelengthDataDirection>%s</WavelengthDataDirection>\n",
                        bsdf_type[(input_orient>0)<<1 | (output_orient>0)]);
        printf("\t\t\t<ColumnAngleBasis>%s</ColumnAngleBasis>\n", kbasis);
        printf("\t\t\t<RowAngleBasis>%s</RowAngleBasis>\n", kbasis);
        puts("\t\t\t<ScatteringDataType>BTDF</ScatteringDataType>");
        puts("\t\t\t<ScatteringData>");
}

/* Output XML data epilogue to stdout */
static void
data_epilogue(void)
{
        puts("\t\t\t</ScatteringData>");
        puts("\t\t</WavelengthDataBlock>");
        puts("\t</WavelengthData>");
}

/* Output XML epilogue to stdout */
static void
xml_epilogue(void)
{
        puts("</Layer>");
        puts("</Optical>");
        puts("</WindowElement>");
}

/* Load and resample XML BSDF description using Klems basis */
static void
eval_bsdf(const char *fname)
{
        ANGLE_BASIS     *abp = get_basis(kbasis);
        SDData          bsd;
        SDError         ec;
        FVECT           vin, vout;
        SDValue         sv;
        double          sum;
        int             i, j, n;

        initurand(npsamps);
        SDclearBSDF(&bsd, fname);               /* load BSDF file */
        if ((ec = SDloadFile(&bsd, fname)) != SDEnone)
                goto err;
        xml_prologue(&bsd);                     /* pass geometry */
                                                /* front reflection */
        if (bsd.rf != NULL || bsd.rLambFront.cieY > .002) {
            input_orient = 1; output_orient = 1;
            data_prologue();
            for (j = 0; j < abp->nangles; j++) {
                for (i = 0; i < abp->nangles; i++) {
                    sum = 0;                    /* average over patches */
                    for (n = npsamps; n-- > 0; ) {
                        fo_getvec(vout, j+(n+frandom())/npsamps, abp);
                        fi_getvec(vin, i+urand(n), abp);
                        ec = SDevalBSDF(&sv, vout, vin, &bsd);
                        if (ec != SDEnone)
                                goto err;
                        sum += sv.cieY;
                    }
                    printf("\t%.3e\n", sum/npsamps);
                }
                putchar('\n');                  /* extra space between rows */
            }
            data_epilogue();
        }
                                                /* back reflection */
        if (bsd.rb != NULL || bsd.rLambBack.cieY > .002) {
            input_orient = -1; output_orient = -1;
            data_prologue();
            for (j = 0; j < abp->nangles; j++) {
                for (i = 0; i < abp->nangles; i++) {
                    sum = 0;                    /* average over patches */
                    for (n = npsamps; n-- > 0; ) {
                        bo_getvec(vout, j+(n+frandom())/npsamps, abp);
                        bi_getvec(vin, i+urand(n), abp);
                        ec = SDevalBSDF(&sv, vout, vin, &bsd);
                        if (ec != SDEnone)
                                goto err;
                        sum += sv.cieY;
                    }
                    printf("\t%.3e\n", sum/npsamps);
                }
                putchar('\n');                  /* extra space between rows */
            }
            data_epilogue();
        }
                                                /* front transmission */
        if (bsd.tf != NULL || bsd.tLamb.cieY > .002) {
            input_orient = 1; output_orient = -1;
            data_prologue();
            for (j = 0; j < abp->nangles; j++) {
                for (i = 0; i < abp->nangles; i++) {
                    sum = 0;                    /* average over patches */
                    for (n = npsamps; n-- > 0; ) {
                        bo_getvec(vout, j+(n+frandom())/npsamps, abp);
                        fi_getvec(vin, i+urand(n), abp);
                        ec = SDevalBSDF(&sv, vout, vin, &bsd);
                        if (ec != SDEnone)
                                goto err;
                        sum += sv.cieY;
                    }
                    printf("\t%.3e\n", sum/npsamps);
                }
                putchar('\n');                  /* extra space between rows */
            }
            data_epilogue();
        }
                                                /* back transmission */
        if ((bsd.tb != NULL) | (bsd.tf != NULL)) {
            input_orient = -1; output_orient = 1;
            data_prologue();
            for (j = 0; j < abp->nangles; j++) {
                for (i = 0; i < abp->nangles; i++) {
                    sum = 0;            /* average over patches */
                    for (n = npsamps; n-- > 0; ) {
                        fo_getvec(vout, j+(n+frandom())/npsamps, abp);
                        bi_getvec(vin, i+urand(n), abp);
                        ec = SDevalBSDF(&sv, vout, vin, &bsd);
                        if (ec != SDEnone)
                                goto err;
                        sum += sv.cieY;
                    }
                    printf("\t%.3e\n", sum/npsamps);
                }
                putchar('\n');                  /* extra space between rows */
            }
            data_epilogue();
        }
        SDfreeBSDF(&bsd);                       /* all done */
        return;
err:
        SDreportError(ec, stderr);
        exit(1);
}

/* Interpolate and output a BSDF function using Klems basis */
static void
eval_function(char *funame)
{
        ANGLE_BASIS     *abp = get_basis(kbasis);
        int             assignD = (fundefined(funame) < 6);
        double          iovec[6];
        double          sum;
        int             i, j, n;

        initurand(npsamps);
        data_prologue();                        /* begin output */
        for (j = 0; j < abp->nangles; j++) {    /* run through directions */
            for (i = 0; i < abp->nangles; i++) {
                sum = 0;
                for (n = npsamps; n--; ) {      /* average over patches */
                    if (output_orient > 0)
                        fo_getvec(iovec+3, j+(n+frandom())/npsamps, abp);
                    else
                        bo_getvec(iovec+3, j+(n+frandom())/npsamps, abp);

                    if (input_orient > 0)
                        fi_getvec(iovec, i+urand(n), abp);
                    else
                        bi_getvec(iovec, i+urand(n), abp);

                    if (assignD) {
                        varset("Dx", '=', -iovec[3]);
                        varset("Dy", '=', -iovec[4]);
                        varset("Dz", '=', -iovec[5]);
                        ++eclock;
                    }
                    sum += funvalue(funame, 6, iovec);
                }
                printf("\t%.3e\n", sum/npsamps);
            }
            putchar('\n');
            prog_show((j+1.)/abp->nangles);
        }
        data_epilogue();                        /* finish output */
        prog_done();
}

/* Interpolate and output a radial basis function BSDF representation */
static void
eval_rbf(void)
{
        ANGLE_BASIS     *abp = get_basis(kbasis);
        float           bsdfarr[MAXPATCHES*MAXPATCHES];
        FVECT           vin, vout;
        RBFNODE         *rbf;
        double          sum;
        int             i, j, n;
                                                /* sanity check */
        if (abp->nangles > MAXPATCHES) {
                fprintf(stderr, "%s: too many patches!\n", progname);
                exit(1);
        }
        data_prologue();                        /* begin output */
        for (i = 0; i < abp->nangles; i++) {
            if (input_orient > 0)               /* use incident patch center */
                fi_getvec(vin, i+.5*(i>0), abp);
            else
                bi_getvec(vin, i+.5*(i>0), abp);

            rbf = advect_rbf(vin, lobe_lim);    /* compute radial basis func */

            for (j = 0; j < abp->nangles; j++) {
                sum = 0;                        /* sample over exiting patch */
                for (n = npsamps; n--; ) {
                    if (output_orient > 0)
                        fo_getvec(vout, j+(n+frandom())/npsamps, abp);
                    else
                        bo_getvec(vout, j+(n+frandom())/npsamps, abp);

                    sum += eval_rbfrep(rbf, vout);
                }
                bsdfarr[j*abp->nangles + i] = sum / (double)npsamps;
            }
            if (rbf != NULL)
                free(rbf);
            prog_show((i+1.)/abp->nangles);
        }
        n = 0;                                  /* write out our matrix */
        for (j = 0; j < abp->nangles; j++) {
            for (i = 0; i < abp->nangles; i++)
                printf("\t%.3e\n", bsdfarr[n++]);
            putchar('\n');
        }
        data_epilogue();                        /* finish output */
        prog_done();
}

/* Read in BSDF and interpolate as Klems matrix representation */
int
main(int argc, char *argv[])
{
        int     dofwd = 0, dobwd = 1;
        char    *cp;
        int     i, na;

        progname = argv[0];
        esupport |= E_VARIABLE|E_FUNCTION|E_RCONST;
        esupport &= ~(E_INCHAN|E_OUTCHAN);
        scompile("PI:3.14159265358979323846", NULL, 0);
        biggerlib();
        for (i = 1; i < argc && (argv[i][0] == '-') | (argv[i][0] == '+'); i++)
                switch (argv[i][1]) {           /* get options */
                case 'n':
                        npsamps = atoi(argv[++i]);
                        if (npsamps <= 0)
                                goto userr;
                        break;
                case 'e':
                        scompile(argv[++i], NULL, 0);
                        single_plane_incident = 0;
                        break;
                case 'f':
                        if (!argv[i][2]) {
                                fcompile(argv[++i]);
                                single_plane_incident = 0;
                        } else
                                dofwd = (argv[i][0] == '+');
                        break;
                case 'b':
                        dobwd = (argv[i][0] == '+');
                        break;
                case 'h':
                        kbasis = "LBNL/Klems Half";
                        break;
                case 'q':
                        kbasis = "LBNL/Klems Quarter";
                        break;
                case 'l':
                        lobe_lim = atoi(argv[++i]);
                        break;
                case 'p':
                        do_prog = atoi(argv[i]+2);
                        break;
                default:
                        goto userr;
                }
        if (single_plane_incident >= 0) {       /* function-based BSDF? */
                if (i != argc-1 || fundefined(argv[i]) < 3) {
                        fprintf(stderr,
        "%s: need single function with 6 arguments: bsdf(ix,iy,iz,ox,oy,oz)\n",
                                        progname);
                        fprintf(stderr, "\tor 3 arguments using Dx,Dy,Dz: bsdf(ix,iy,iz)\n");
                        goto userr;
                }
                ++eclock;
                xml_header(argc, argv);                 /* start XML output */
                xml_prologue(NULL);
                if (dofwd) {
                        input_orient = -1;
                        output_orient = -1;
                        prog_start("Evaluating outside reflectance");
                        eval_function(argv[i]);
                        output_orient = 1;
                        prog_start("Evaluating outside->inside transmission");
                        eval_function(argv[i]);
                }
                if (dobwd) {
                        input_orient = 1;
                        output_orient = 1;
                        prog_start("Evaluating inside reflectance");
                        eval_function(argv[i]);
                        output_orient = -1;
                        prog_start("Evaluating inside->outside transmission");
                        eval_function(argv[i]);
                }
                xml_epilogue();                 /* finish XML output & exit */
                return(0);
        }
                                                /* XML input? */
        if (i == argc-1 && (cp = argv[i]+strlen(argv[i])-4) > argv[i] &&
                                !strcasecmp(cp, ".xml")) {
                xml_header(argc, argv);         /* start XML output */
                eval_bsdf(argv[i]);             /* load & resample BSDF */
                xml_epilogue();                 /* finish XML output & exit */
                return(0);
        }
        if (i < argc) {                         /* open input files if given */
                int     nbsdf = 0;
                for ( ; i < argc; i++) {        /* interpolate each component */
                        char    pbuf[256];
                        FILE    *fpin = fopen(argv[i], "rb");
                        if (fpin == NULL) {
                                fprintf(stderr, "%s: cannot open BSDF interpolant '%s'\n",
                                                progname, argv[i]);
                                return(1);
                        }
                        if (!load_bsdf_rep(fpin))
                                return(1);
                        fclose(fpin);
                        if (!nbsdf++) {         /* start XML on first dist. */
                                xml_header(argc, argv);
                                xml_prologue(NULL);
                        }
                        sprintf(pbuf, "Interpolating component '%s'", argv[i]);
                        prog_start(pbuf);
                        eval_rbf();
                }
                xml_epilogue();                 /* finish XML output & exit */
                return(0);
        }
        SET_FILE_BINARY(stdin);                 /* load from stdin */
        if (!load_bsdf_rep(stdin))
                return(1);
        xml_header(argc, argv);                 /* start XML output */
        xml_prologue(NULL);
        prog_start("Interpolating from standard input");
        eval_rbf();                             /* resample dist. */
        xml_epilogue();                         /* finish XML output & exit */
        return(0);
userr:
        fprintf(stderr,
        "Usage: %s [-n spp][-h|-q][-l maxlobes] [bsdf.sir ..] > bsdf.xml\n", progname);
        fprintf(stderr,
        "   or: %s [-n spp][-h|-q] bsdf_in.xml > bsdf_out.xml\n", progname);
        fprintf(stderr,
        "   or: %s [-n spp][-h|-q][{+|-}for[ward]][{+|-}b[ackward]][-e expr][-f file] bsdf_func > bsdf.xml\n",
                                progname);
        return(1);
}
Revision:	2.19
Committed:	Sun Jan 24 14:28:51 2016 UTC (8 years, 3 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.18:	+2 -2 lines
Log Message:	Fixed check for number of arguments to allow 3-argument BSDFs
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.19	static const char RCSid[] = "$Id: bsdf2klems.c,v 2.18 2015/05/27 11:39:37 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	greg	2.4	/* assumed maximum # Klems patches */
21			#define MAXPATCHES 145
22	greg	2.1	/* global argv[0] */
23			char *progname;
24			/* selected basis function name */
25			static const char *kbasis = "LBNL/Klems Full";
26			/* number of BSDF samples per patch */
27			static int npsamps = 256;
28	greg	2.10	/* limit on number of RBF lobes */
29			static int lobe_lim = 15000;
30	greg	2.13	/* progress bar length */
31			static int do_prog = 79;
32
33
34			/* Start new progress bar */
35			#define prog_start(s) if (do_prog) fprintf(stderr, "%s: %s...\n", progname, s); else
36
37			/* Draw progress bar of the appropriate length */
38			static void
39			prog_show(double frac)
40			{
41	greg	2.17	static unsigned call_cnt = 0;
42			static char lastc[] = "-\\\|/";
43			char pbar[256];
44			int nchars;
45	greg	2.13
46	greg	2.14	if (do_prog <= 1) return;
47	greg	2.13	if (do_prog > sizeof(pbar)-2)
48			do_prog = sizeof(pbar)-2;
49			if (frac < 0) frac = 0;
50	greg	2.17	else if (frac >= 1) frac = .9999;
51			nchars = do_prog*frac;
52	greg	2.13	pbar[0] = '\r';
53			memset(pbar+1, '*', nchars);
54	greg	2.17	pbar[nchars+1] = lastc[call_cnt++ & 3];
55			memset(pbar+2+nchars, '-', do_prog-nchars-1);
56	greg	2.13	pbar[do_prog+1] = '\0';
57			fputs(pbar, stderr);
58			}
59
60			/* Finish progress bar */
61	greg	2.14	static void
62			prog_done(void)
63			{
64			int n = do_prog;
65
66			if (n <= 1) return;
67			fputc('\r', stderr);
68			while (n--)
69			fputc(' ', stderr);
70			fputc('\r', stderr);
71			}
72	greg	2.1
73			/* Return angle basis corresponding to the given name */
74	greg	2.13	static ANGLE_BASIS *
75	greg	2.1	get_basis(const char *bn)
76			{
77			int n = nabases;
78
79			while (n-- > 0)
80			if (!strcasecmp(bn, abase_list[n].name))
81			return &abase_list[n];
82			return NULL;
83			}
84
85	greg	2.3	/* Output XML header to stdout */
86			static void
87			xml_header(int ac, char *av[])
88			{
89			puts("<?xml version=\"1.0\" encoding=\"UTF-8\"?>");
90			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\">");
91			fputs("<!-- File produced by:", stdout);
92			while (ac-- > 0) {
93			fputc(' ', stdout);
94			fputs(*av++, stdout);
95			}
96			puts(" -->");
97			}
98
99	greg	2.1	/* Output XML prologue to stdout */
100			static void
101	greg	2.3	xml_prologue(const SDData *sd)
102	greg	2.1	{
103	greg	2.11	const char *matn = (sd && sd->matn[0]) ? sd->matn :
104			bsdf_name[0] ? bsdf_name : "Unknown";
105			const char *makr = (sd && sd->makr[0]) ? sd->makr :
106			bsdf_manuf[0] ? bsdf_manuf : "Unknown";
107	greg	2.1	ANGLE_BASIS *abp = get_basis(kbasis);
108			int i;
109
110			if (abp == NULL) {
111			fprintf(stderr, "%s: unknown angle basis '%s'\n", progname, kbasis);
112			exit(1);
113			}
114			puts("<WindowElementType>System</WindowElementType>");
115			puts("<FileType>BSDF</FileType>");
116			puts("<Optical>");
117			puts("<Layer>");
118			puts("\t<Material>");
119	greg	2.3	printf("\t\t<Name>%s</Name>\n", matn);
120			printf("\t\t<Manufacturer>%s</Manufacturer>\n", makr);
121			if (sd && sd->dim[2] > .001) {
122			printf("\t\t<Thickness unit=\"meter\">%.3f</Thickness>\n", sd->dim[2]);
123			printf("\t\t<Width unit=\"meter\">%.3f</Width>\n", sd->dim[0]);
124			printf("\t\t<Height unit=\"meter\">%.3f</Height>\n", sd->dim[1]);
125			}
126	greg	2.1	puts("\t\t<DeviceType>Other</DeviceType>");
127			puts("\t</Material>");
128	greg	2.3	if (sd && sd->mgf != NULL) {
129			puts("\t<Geometry format=\"MGF\">");
130			puts("\t\t<MGFblock unit=\"meter\">");
131			fputs(sd->mgf, stdout);
132			puts("</MGFblock>");
133			puts("\t</Geometry>");
134			}
135	greg	2.1	puts("\t<DataDefinition>");
136			puts("\t\t<IncidentDataStructure>Columns</IncidentDataStructure>");
137			puts("\t\t<AngleBasis>");
138			printf("\t\t\t<AngleBasisName>%s</AngleBasisName>\n", kbasis);
139			for (i = 0; abp->lat[i].nphis; i++) {
140			puts("\t\t\t<AngleBasisBlock>");
141			printf("\t\t\t<Theta>%g</Theta>\n", i ?
142			.5*(abp->lat[i].tmin + abp->lat[i+1].tmin) :
143			.0 );
144	greg	2.3	printf("\t\t\t<nPhis>%d</nPhis>\n", abp->lat[i].nphis);
145	greg	2.1	puts("\t\t\t<ThetaBounds>");
146			printf("\t\t\t\t<LowerTheta>%g</LowerTheta>\n", abp->lat[i].tmin);
147			printf("\t\t\t\t<UpperTheta>%g</UpperTheta>\n", abp->lat[i+1].tmin);
148			puts("\t\t\t</ThetaBounds>");
149			puts("\t\t\t</AngleBasisBlock>");
150			}
151			puts("\t\t</AngleBasis>");
152			puts("\t</DataDefinition>");
153			}
154
155			/* Output XML data prologue to stdout */
156			static void
157			data_prologue()
158			{
159			static const char *bsdf_type[4] = {
160			"Reflection Front",
161			"Transmission Front",
162			"Transmission Back",
163			"Reflection Back"
164			};
165
166			puts("\t<WavelengthData>");
167			puts("\t\t<LayerNumber>System</LayerNumber>");
168			puts("\t\t<Wavelength unit=\"Integral\">Visible</Wavelength>");
169			puts("\t\t<SourceSpectrum>CIE Illuminant D65 1nm.ssp</SourceSpectrum>");
170			puts("\t\t<DetectorSpectrum>ASTM E308 1931 Y.dsp</DetectorSpectrum>");
171			puts("\t\t<WavelengthDataBlock>");
172			printf("\t\t\t<WavelengthDataDirection>%s</WavelengthDataDirection>\n",
173			bsdf_type[(input_orient>0)<<1 \| (output_orient>0)]);
174			printf("\t\t\t<ColumnAngleBasis>%s</ColumnAngleBasis>\n", kbasis);
175			printf("\t\t\t<RowAngleBasis>%s</RowAngleBasis>\n", kbasis);
176			puts("\t\t\t<ScatteringDataType>BTDF</ScatteringDataType>");
177			puts("\t\t\t<ScatteringData>");
178			}
179
180			/* Output XML data epilogue to stdout */
181			static void
182			data_epilogue(void)
183			{
184			puts("\t\t\t</ScatteringData>");
185			puts("\t\t</WavelengthDataBlock>");
186			puts("\t</WavelengthData>");
187			}
188
189			/* Output XML epilogue to stdout */
190			static void
191			xml_epilogue(void)
192			{
193			puts("</Layer>");
194			puts("</Optical>");
195			puts("</WindowElement>");
196			}
197
198	greg	2.2	/* Load and resample XML BSDF description using Klems basis */
199	greg	2.1	static void
200			eval_bsdf(const char *fname)
201			{
202			ANGLE_BASIS *abp = get_basis(kbasis);
203			SDData bsd;
204			SDError ec;
205			FVECT vin, vout;
206			SDValue sv;
207			double sum;
208			int i, j, n;
209
210	greg	2.18	initurand(npsamps);
211	greg	2.1	SDclearBSDF(&bsd, fname); /* load BSDF file */
212			if ((ec = SDloadFile(&bsd, fname)) != SDEnone)
213			goto err;
214	greg	2.3	xml_prologue(&bsd); /* pass geometry */
215	greg	2.1	/* front reflection */
216			if (bsd.rf != NULL \|\| bsd.rLambFront.cieY > .002) {
217			input_orient = 1; output_orient = 1;
218			data_prologue();
219			for (j = 0; j < abp->nangles; j++) {
220			for (i = 0; i < abp->nangles; i++) {
221			sum = 0; /* average over patches */
222			for (n = npsamps; n-- > 0; ) {
223			fo_getvec(vout, j+(n+frandom())/npsamps, abp);
224	greg	2.5	fi_getvec(vin, i+urand(n), abp);
225	greg	2.1	ec = SDevalBSDF(&sv, vout, vin, &bsd);
226			if (ec != SDEnone)
227			goto err;
228			sum += sv.cieY;
229			}
230			printf("\t%.3e\n", sum/npsamps);
231			}
232			putchar('\n'); /* extra space between rows */
233			}
234			data_epilogue();
235			}
236			/* back reflection */
237			if (bsd.rb != NULL \|\| bsd.rLambBack.cieY > .002) {
238			input_orient = -1; output_orient = -1;
239			data_prologue();
240			for (j = 0; j < abp->nangles; j++) {
241			for (i = 0; i < abp->nangles; i++) {
242			sum = 0; /* average over patches */
243			for (n = npsamps; n-- > 0; ) {
244			bo_getvec(vout, j+(n+frandom())/npsamps, abp);
245	greg	2.5	bi_getvec(vin, i+urand(n), abp);
246	greg	2.1	ec = SDevalBSDF(&sv, vout, vin, &bsd);
247			if (ec != SDEnone)
248			goto err;
249			sum += sv.cieY;
250			}
251			printf("\t%.3e\n", sum/npsamps);
252			}
253			putchar('\n'); /* extra space between rows */
254			}
255			data_epilogue();
256			}
257			/* front transmission */
258			if (bsd.tf != NULL \|\| bsd.tLamb.cieY > .002) {
259			input_orient = 1; output_orient = -1;
260			data_prologue();
261			for (j = 0; j < abp->nangles; j++) {
262			for (i = 0; i < abp->nangles; i++) {
263			sum = 0; /* average over patches */
264			for (n = npsamps; n-- > 0; ) {
265			bo_getvec(vout, j+(n+frandom())/npsamps, abp);
266	greg	2.5	fi_getvec(vin, i+urand(n), abp);
267	greg	2.1	ec = SDevalBSDF(&sv, vout, vin, &bsd);
268			if (ec != SDEnone)
269			goto err;
270			sum += sv.cieY;
271			}
272			printf("\t%.3e\n", sum/npsamps);
273			}
274			putchar('\n'); /* extra space between rows */
275			}
276			data_epilogue();
277			}
278			/* back transmission */
279	greg	2.9	if ((bsd.tb != NULL) \| (bsd.tf != NULL)) {
280	greg	2.1	input_orient = -1; output_orient = 1;
281			data_prologue();
282			for (j = 0; j < abp->nangles; j++) {
283	greg	2.9	for (i = 0; i < abp->nangles; i++) {
284			sum = 0; /* average over patches */
285			for (n = npsamps; n-- > 0; ) {
286			fo_getvec(vout, j+(n+frandom())/npsamps, abp);
287			bi_getvec(vin, i+urand(n), abp);
288			ec = SDevalBSDF(&sv, vout, vin, &bsd);
289			if (ec != SDEnone)
290	greg	2.1	goto err;
291	greg	2.9	sum += sv.cieY;
292	greg	2.1	}
293	greg	2.9	printf("\t%.3e\n", sum/npsamps);
294			}
295	greg	2.1	putchar('\n'); /* extra space between rows */
296			}
297			data_epilogue();
298			}
299			SDfreeBSDF(&bsd); /* all done */
300			return;
301			err:
302			SDreportError(ec, stderr);
303			exit(1);
304			}
305
306	greg	2.2	/* Interpolate and output a BSDF function using Klems basis */
307	greg	2.1	static void
308			eval_function(char *funame)
309			{
310			ANGLE_BASIS *abp = get_basis(kbasis);
311	greg	2.8	int assignD = (fundefined(funame) < 6);
312	greg	2.1	double iovec[6];
313			double sum;
314			int i, j, n;
315
316	greg	2.4	initurand(npsamps);
317	greg	2.1	data_prologue(); /* begin output */
318			for (j = 0; j < abp->nangles; j++) { /* run through directions */
319			for (i = 0; i < abp->nangles; i++) {
320			sum = 0;
321			for (n = npsamps; n--; ) { /* average over patches */
322			if (output_orient > 0)
323			fo_getvec(iovec+3, j+(n+frandom())/npsamps, abp);
324			else
325			bo_getvec(iovec+3, j+(n+frandom())/npsamps, abp);
326
327			if (input_orient > 0)
328	greg	2.4	fi_getvec(iovec, i+urand(n), abp);
329	greg	2.1	else
330	greg	2.4	bi_getvec(iovec, i+urand(n), abp);
331	greg	2.1
332	greg	2.8	if (assignD) {
333			varset("Dx", '=', -iovec[3]);
334			varset("Dy", '=', -iovec[4]);
335			varset("Dz", '=', -iovec[5]);
336			++eclock;
337			}
338	greg	2.1	sum += funvalue(funame, 6, iovec);
339			}
340			printf("\t%.3e\n", sum/npsamps);
341			}
342			putchar('\n');
343	greg	2.13	prog_show((j+1.)/abp->nangles);
344	greg	2.1	}
345			data_epilogue(); /* finish output */
346	greg	2.13	prog_done();
347	greg	2.1	}
348
349			/* Interpolate and output a radial basis function BSDF representation */
350			static void
351			eval_rbf(void)
352			{
353			ANGLE_BASIS *abp = get_basis(kbasis);
354	greg	2.2	float bsdfarr[MAXPATCHES*MAXPATCHES];
355			FVECT vin, vout;
356			RBFNODE *rbf;
357	greg	2.1	double sum;
358			int i, j, n;
359	greg	2.2	/* sanity check */
360			if (abp->nangles > MAXPATCHES) {
361			fprintf(stderr, "%s: too many patches!\n", progname);
362			exit(1);
363			}
364			data_prologue(); /* begin output */
365			for (i = 0; i < abp->nangles; i++) {
366			if (input_orient > 0) /* use incident patch center */
367			fi_getvec(vin, i+.5*(i>0), abp);
368			else
369			bi_getvec(vin, i+.5*(i>0), abp);
370
371	greg	2.10	rbf = advect_rbf(vin, lobe_lim); /* compute radial basis func */
372	greg	2.2
373			for (j = 0; j < abp->nangles; j++) {
374			sum = 0; /* sample over exiting patch */
375			for (n = npsamps; n--; ) {
376			if (output_orient > 0)
377			fo_getvec(vout, j+(n+frandom())/npsamps, abp);
378			else
379			bo_getvec(vout, j+(n+frandom())/npsamps, abp);
380	greg	2.1
381	greg	2.12	sum += eval_rbfrep(rbf, vout);
382	greg	2.2	}
383	greg	2.15	bsdfarr[j*abp->nangles + i] = sum / (double)npsamps;
384	greg	2.2	}
385	greg	2.6	if (rbf != NULL)
386			free(rbf);
387	greg	2.13	prog_show((i+1.)/abp->nangles);
388	greg	2.2	}
389			n = 0; /* write out our matrix */
390			for (j = 0; j < abp->nangles; j++) {
391			for (i = 0; i < abp->nangles; i++)
392			printf("\t%.3e\n", bsdfarr[n++]);
393			putchar('\n');
394			}
395			data_epilogue(); /* finish output */
396	greg	2.13	prog_done();
397	greg	2.1	}
398
399			/* Read in BSDF and interpolate as Klems matrix representation */
400			int
401			main(int argc, char *argv[])
402			{
403			int dofwd = 0, dobwd = 1;
404			char *cp;
405			int i, na;
406
407			progname = argv[0];
408			esupport \|= E_VARIABLE\|E_FUNCTION\|E_RCONST;
409			esupport &= ~(E_INCHAN\|E_OUTCHAN);
410			scompile("PI:3.14159265358979323846", NULL, 0);
411			biggerlib();
412	greg	2.2	for (i = 1; i < argc && (argv[i][0] == '-') \| (argv[i][0] == '+'); i++)
413	greg	2.1	switch (argv[i][1]) { /* get options */
414			case 'n':
415			npsamps = atoi(argv[++i]);
416			if (npsamps <= 0)
417			goto userr;
418			break;
419			case 'e':
420			scompile(argv[++i], NULL, 0);
421			single_plane_incident = 0;
422			break;
423			case 'f':
424			if (!argv[i][2]) {
425			fcompile(argv[++i]);
426			single_plane_incident = 0;
427			} else
428			dofwd = (argv[i][0] == '+');
429			break;
430			case 'b':
431			dobwd = (argv[i][0] == '+');
432			break;
433			case 'h':
434			kbasis = "LBNL/Klems Half";
435			break;
436			case 'q':
437			kbasis = "LBNL/Klems Quarter";
438			break;
439	greg	2.10	case 'l':
440			lobe_lim = atoi(argv[++i]);
441			break;
442	greg	2.13	case 'p':
443			do_prog = atoi(argv[i]+2);
444			break;
445	greg	2.1	default:
446			goto userr;
447			}
448			if (single_plane_incident >= 0) { /* function-based BSDF? */
449	greg	2.19	if (i != argc-1 \|\| fundefined(argv[i]) < 3) {
450	greg	2.1	fprintf(stderr,
451			"%s: need single function with 6 arguments: bsdf(ix,iy,iz,ox,oy,oz)\n",
452			progname);
453	greg	2.13	fprintf(stderr, "\tor 3 arguments using Dx,Dy,Dz: bsdf(ix,iy,iz)\n");
454	greg	2.1	goto userr;
455			}
456	greg	2.8	++eclock;
457	greg	2.3	xml_header(argc, argv); /* start XML output */
458			xml_prologue(NULL);
459	greg	2.1	if (dofwd) {
460			input_orient = -1;
461			output_orient = -1;
462	greg	2.13	prog_start("Evaluating outside reflectance");
463			eval_function(argv[i]);
464	greg	2.1	output_orient = 1;
465	greg	2.13	prog_start("Evaluating outside->inside transmission");
466			eval_function(argv[i]);
467	greg	2.1	}
468			if (dobwd) {
469			input_orient = 1;
470			output_orient = 1;
471	greg	2.13	prog_start("Evaluating inside reflectance");
472			eval_function(argv[i]);
473	greg	2.1	output_orient = -1;
474	greg	2.13	prog_start("Evaluating inside->outside transmission");
475			eval_function(argv[i]);
476	greg	2.1	}
477			xml_epilogue(); /* finish XML output & exit */
478			return(0);
479			}
480	greg	2.2	/* XML input? */
481			if (i == argc-1 && (cp = argv[i]+strlen(argv[i])-4) > argv[i] &&
482			!strcasecmp(cp, ".xml")) {
483	greg	2.3	xml_header(argc, argv); /* start XML output */
484	greg	2.1	eval_bsdf(argv[i]); /* load & resample BSDF */
485			xml_epilogue(); /* finish XML output & exit */
486			return(0);
487			}
488			if (i < argc) { /* open input files if given */
489			int nbsdf = 0;
490			for ( ; i < argc; i++) { /* interpolate each component */
491	greg	2.13	char pbuf[256];
492	greg	2.1	FILE *fpin = fopen(argv[i], "rb");
493			if (fpin == NULL) {
494			fprintf(stderr, "%s: cannot open BSDF interpolant '%s'\n",
495			progname, argv[i]);
496			return(1);
497			}
498			if (!load_bsdf_rep(fpin))
499			return(1);
500			fclose(fpin);
501	greg	2.3	if (!nbsdf++) { /* start XML on first dist. */
502			xml_header(argc, argv);
503			xml_prologue(NULL);
504			}
505	greg	2.13	sprintf(pbuf, "Interpolating component '%s'", argv[i]);
506			prog_start(pbuf);
507	greg	2.1	eval_rbf();
508			}
509			xml_epilogue(); /* finish XML output & exit */
510			return(0);
511			}
512			SET_FILE_BINARY(stdin); /* load from stdin */
513			if (!load_bsdf_rep(stdin))
514			return(1);
515	greg	2.3	xml_header(argc, argv); /* start XML output */
516			xml_prologue(NULL);
517	greg	2.13	prog_start("Interpolating from standard input");
518	greg	2.1	eval_rbf(); /* resample dist. */
519			xml_epilogue(); /* finish XML output & exit */
520			return(0);
521			userr:
522			fprintf(stderr,
523	greg	2.10	"Usage: %s [-n spp][-h\|-q][-l maxlobes] [bsdf.sir ..] > bsdf.xml\n", progname);
524	greg	2.1	fprintf(stderr,
525	greg	2.3	" or: %s [-n spp][-h\|-q] bsdf_in.xml > bsdf_out.xml\n", progname);
526	greg	2.1	fprintf(stderr,
527			" or: %s [-n spp][-h\|-q][{+\|-}for[ward]][{+\|-}b[ackward]][-e expr][-f file] bsdf_func > bsdf.xml\n",
528			progname);
529			return(1);
530			}