src/cv/bsdf2klems.c

#ifndef lint
static const char RCSid[] = "$Id: bsdf2klems.c,v 2.16 2014/10/01 23:32:42 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;

        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]) != 6) {
                        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.17
Committed:	Tue May 5 22:16:49 2015 UTC (10 years ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.16:	+9 -6 lines
Log Message:	Improved progress bar feedback
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.17	static const char RCSid[] = "$Id: bsdf2klems.c,v 2.16 2014/10/01 23:32:42 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			SDclearBSDF(&bsd, fname); /* load BSDF file */
211			if ((ec = SDloadFile(&bsd, fname)) != SDEnone)
212			goto err;
213	greg	2.3	xml_prologue(&bsd); /* pass geometry */
214	greg	2.1	/* front reflection */
215			if (bsd.rf != NULL \|\| bsd.rLambFront.cieY > .002) {
216			input_orient = 1; output_orient = 1;
217			data_prologue();
218			for (j = 0; j < abp->nangles; j++) {
219			for (i = 0; i < abp->nangles; i++) {
220			sum = 0; /* average over patches */
221			for (n = npsamps; n-- > 0; ) {
222			fo_getvec(vout, j+(n+frandom())/npsamps, abp);
223	greg	2.5	fi_getvec(vin, i+urand(n), abp);
224	greg	2.1	ec = SDevalBSDF(&sv, vout, vin, &bsd);
225			if (ec != SDEnone)
226			goto err;
227			sum += sv.cieY;
228			}
229			printf("\t%.3e\n", sum/npsamps);
230			}
231			putchar('\n'); /* extra space between rows */
232			}
233			data_epilogue();
234			}
235			/* back reflection */
236			if (bsd.rb != NULL \|\| bsd.rLambBack.cieY > .002) {
237			input_orient = -1; output_orient = -1;
238			data_prologue();
239			for (j = 0; j < abp->nangles; j++) {
240			for (i = 0; i < abp->nangles; i++) {
241			sum = 0; /* average over patches */
242			for (n = npsamps; n-- > 0; ) {
243			bo_getvec(vout, j+(n+frandom())/npsamps, abp);
244	greg	2.5	bi_getvec(vin, i+urand(n), abp);
245	greg	2.1	ec = SDevalBSDF(&sv, vout, vin, &bsd);
246			if (ec != SDEnone)
247			goto err;
248			sum += sv.cieY;
249			}
250			printf("\t%.3e\n", sum/npsamps);
251			}
252			putchar('\n'); /* extra space between rows */
253			}
254			data_epilogue();
255			}
256			/* front transmission */
257			if (bsd.tf != NULL \|\| bsd.tLamb.cieY > .002) {
258			input_orient = 1; output_orient = -1;
259			data_prologue();
260			for (j = 0; j < abp->nangles; j++) {
261			for (i = 0; i < abp->nangles; i++) {
262			sum = 0; /* average over patches */
263			for (n = npsamps; n-- > 0; ) {
264			bo_getvec(vout, j+(n+frandom())/npsamps, abp);
265	greg	2.5	fi_getvec(vin, i+urand(n), abp);
266	greg	2.1	ec = SDevalBSDF(&sv, vout, vin, &bsd);
267			if (ec != SDEnone)
268			goto err;
269			sum += sv.cieY;
270			}
271			printf("\t%.3e\n", sum/npsamps);
272			}
273			putchar('\n'); /* extra space between rows */
274			}
275			data_epilogue();
276			}
277			/* back transmission */
278	greg	2.9	if ((bsd.tb != NULL) \| (bsd.tf != NULL)) {
279	greg	2.1	input_orient = -1; output_orient = 1;
280			data_prologue();
281			for (j = 0; j < abp->nangles; j++) {
282	greg	2.9	for (i = 0; i < abp->nangles; i++) {
283			sum = 0; /* average over patches */
284			for (n = npsamps; n-- > 0; ) {
285			fo_getvec(vout, j+(n+frandom())/npsamps, abp);
286			bi_getvec(vin, i+urand(n), abp);
287			ec = SDevalBSDF(&sv, vout, vin, &bsd);
288			if (ec != SDEnone)
289	greg	2.1	goto err;
290	greg	2.9	sum += sv.cieY;
291	greg	2.1	}
292	greg	2.9	printf("\t%.3e\n", sum/npsamps);
293			}
294	greg	2.1	putchar('\n'); /* extra space between rows */
295			}
296			data_epilogue();
297			}
298			SDfreeBSDF(&bsd); /* all done */
299			return;
300			err:
301			SDreportError(ec, stderr);
302			exit(1);
303			}
304
305	greg	2.2	/* Interpolate and output a BSDF function using Klems basis */
306	greg	2.1	static void
307			eval_function(char *funame)
308			{
309			ANGLE_BASIS *abp = get_basis(kbasis);
310	greg	2.8	int assignD = (fundefined(funame) < 6);
311	greg	2.1	double iovec[6];
312			double sum;
313			int i, j, n;
314
315	greg	2.4	initurand(npsamps);
316	greg	2.1	data_prologue(); /* begin output */
317			for (j = 0; j < abp->nangles; j++) { /* run through directions */
318			for (i = 0; i < abp->nangles; i++) {
319			sum = 0;
320			for (n = npsamps; n--; ) { /* average over patches */
321			if (output_orient > 0)
322			fo_getvec(iovec+3, j+(n+frandom())/npsamps, abp);
323			else
324			bo_getvec(iovec+3, j+(n+frandom())/npsamps, abp);
325
326			if (input_orient > 0)
327	greg	2.4	fi_getvec(iovec, i+urand(n), abp);
328	greg	2.1	else
329	greg	2.4	bi_getvec(iovec, i+urand(n), abp);
330	greg	2.1
331	greg	2.8	if (assignD) {
332			varset("Dx", '=', -iovec[3]);
333			varset("Dy", '=', -iovec[4]);
334			varset("Dz", '=', -iovec[5]);
335			++eclock;
336			}
337	greg	2.1	sum += funvalue(funame, 6, iovec);
338			}
339			printf("\t%.3e\n", sum/npsamps);
340			}
341			putchar('\n');
342	greg	2.13	prog_show((j+1.)/abp->nangles);
343	greg	2.1	}
344			data_epilogue(); /* finish output */
345	greg	2.13	prog_done();
346	greg	2.1	}
347
348			/* Interpolate and output a radial basis function BSDF representation */
349			static void
350			eval_rbf(void)
351			{
352			ANGLE_BASIS *abp = get_basis(kbasis);
353	greg	2.2	float bsdfarr[MAXPATCHES*MAXPATCHES];
354			FVECT vin, vout;
355			RBFNODE *rbf;
356	greg	2.1	double sum;
357			int i, j, n;
358	greg	2.2	/* sanity check */
359			if (abp->nangles > MAXPATCHES) {
360			fprintf(stderr, "%s: too many patches!\n", progname);
361			exit(1);
362			}
363			data_prologue(); /* begin output */
364			for (i = 0; i < abp->nangles; i++) {
365			if (input_orient > 0) /* use incident patch center */
366			fi_getvec(vin, i+.5*(i>0), abp);
367			else
368			bi_getvec(vin, i+.5*(i>0), abp);
369
370	greg	2.10	rbf = advect_rbf(vin, lobe_lim); /* compute radial basis func */
371	greg	2.2
372			for (j = 0; j < abp->nangles; j++) {
373			sum = 0; /* sample over exiting patch */
374			for (n = npsamps; n--; ) {
375			if (output_orient > 0)
376			fo_getvec(vout, j+(n+frandom())/npsamps, abp);
377			else
378			bo_getvec(vout, j+(n+frandom())/npsamps, abp);
379	greg	2.1
380	greg	2.12	sum += eval_rbfrep(rbf, vout);
381	greg	2.2	}
382	greg	2.15	bsdfarr[j*abp->nangles + i] = sum / (double)npsamps;
383	greg	2.2	}
384	greg	2.6	if (rbf != NULL)
385			free(rbf);
386	greg	2.13	prog_show((i+1.)/abp->nangles);
387	greg	2.2	}
388			n = 0; /* write out our matrix */
389			for (j = 0; j < abp->nangles; j++) {
390			for (i = 0; i < abp->nangles; i++)
391			printf("\t%.3e\n", bsdfarr[n++]);
392			putchar('\n');
393			}
394			data_epilogue(); /* finish output */
395	greg	2.13	prog_done();
396	greg	2.1	}
397
398			/* Read in BSDF and interpolate as Klems matrix representation */
399			int
400			main(int argc, char *argv[])
401			{
402			int dofwd = 0, dobwd = 1;
403			char *cp;
404			int i, na;
405
406			progname = argv[0];
407			esupport \|= E_VARIABLE\|E_FUNCTION\|E_RCONST;
408			esupport &= ~(E_INCHAN\|E_OUTCHAN);
409			scompile("PI:3.14159265358979323846", NULL, 0);
410			biggerlib();
411	greg	2.2	for (i = 1; i < argc && (argv[i][0] == '-') \| (argv[i][0] == '+'); i++)
412	greg	2.1	switch (argv[i][1]) { /* get options */
413			case 'n':
414			npsamps = atoi(argv[++i]);
415			if (npsamps <= 0)
416			goto userr;
417			break;
418			case 'e':
419			scompile(argv[++i], NULL, 0);
420			single_plane_incident = 0;
421			break;
422			case 'f':
423			if (!argv[i][2]) {
424			fcompile(argv[++i]);
425			single_plane_incident = 0;
426			} else
427			dofwd = (argv[i][0] == '+');
428			break;
429			case 'b':
430			dobwd = (argv[i][0] == '+');
431			break;
432			case 'h':
433			kbasis = "LBNL/Klems Half";
434			break;
435			case 'q':
436			kbasis = "LBNL/Klems Quarter";
437			break;
438	greg	2.10	case 'l':
439			lobe_lim = atoi(argv[++i]);
440			break;
441	greg	2.13	case 'p':
442			do_prog = atoi(argv[i]+2);
443			break;
444	greg	2.1	default:
445			goto userr;
446			}
447			if (single_plane_incident >= 0) { /* function-based BSDF? */
448			if (i != argc-1 \|\| fundefined(argv[i]) != 6) {
449			fprintf(stderr,
450			"%s: need single function with 6 arguments: bsdf(ix,iy,iz,ox,oy,oz)\n",
451			progname);
452	greg	2.13	fprintf(stderr, "\tor 3 arguments using Dx,Dy,Dz: bsdf(ix,iy,iz)\n");
453	greg	2.1	goto userr;
454			}
455	greg	2.8	++eclock;
456	greg	2.3	xml_header(argc, argv); /* start XML output */
457			xml_prologue(NULL);
458	greg	2.1	if (dofwd) {
459			input_orient = -1;
460			output_orient = -1;
461	greg	2.13	prog_start("Evaluating outside reflectance");
462			eval_function(argv[i]);
463	greg	2.1	output_orient = 1;
464	greg	2.13	prog_start("Evaluating outside->inside transmission");
465			eval_function(argv[i]);
466	greg	2.1	}
467			if (dobwd) {
468			input_orient = 1;
469			output_orient = 1;
470	greg	2.13	prog_start("Evaluating inside reflectance");
471			eval_function(argv[i]);
472	greg	2.1	output_orient = -1;
473	greg	2.13	prog_start("Evaluating inside->outside transmission");
474			eval_function(argv[i]);
475	greg	2.1	}
476			xml_epilogue(); /* finish XML output & exit */
477			return(0);
478			}
479	greg	2.2	/* XML input? */
480			if (i == argc-1 && (cp = argv[i]+strlen(argv[i])-4) > argv[i] &&
481			!strcasecmp(cp, ".xml")) {
482	greg	2.3	xml_header(argc, argv); /* start XML output */
483	greg	2.1	eval_bsdf(argv[i]); /* load & resample BSDF */
484			xml_epilogue(); /* finish XML output & exit */
485			return(0);
486			}
487			if (i < argc) { /* open input files if given */
488			int nbsdf = 0;
489			for ( ; i < argc; i++) { /* interpolate each component */
490	greg	2.13	char pbuf[256];
491	greg	2.1	FILE *fpin = fopen(argv[i], "rb");
492			if (fpin == NULL) {
493			fprintf(stderr, "%s: cannot open BSDF interpolant '%s'\n",
494			progname, argv[i]);
495			return(1);
496			}
497			if (!load_bsdf_rep(fpin))
498			return(1);
499			fclose(fpin);
500	greg	2.3	if (!nbsdf++) { /* start XML on first dist. */
501			xml_header(argc, argv);
502			xml_prologue(NULL);
503			}
504	greg	2.13	sprintf(pbuf, "Interpolating component '%s'", argv[i]);
505			prog_start(pbuf);
506	greg	2.1	eval_rbf();
507			}
508			xml_epilogue(); /* finish XML output & exit */
509			return(0);
510			}
511			SET_FILE_BINARY(stdin); /* load from stdin */
512			if (!load_bsdf_rep(stdin))
513			return(1);
514	greg	2.3	xml_header(argc, argv); /* start XML output */
515			xml_prologue(NULL);
516	greg	2.13	prog_start("Interpolating from standard input");
517	greg	2.1	eval_rbf(); /* resample dist. */
518			xml_epilogue(); /* finish XML output & exit */
519			return(0);
520			userr:
521			fprintf(stderr,
522	greg	2.10	"Usage: %s [-n spp][-h\|-q][-l maxlobes] [bsdf.sir ..] > bsdf.xml\n", progname);
523	greg	2.1	fprintf(stderr,
524	greg	2.3	" or: %s [-n spp][-h\|-q] bsdf_in.xml > bsdf_out.xml\n", progname);
525	greg	2.1	fprintf(stderr,
526			" or: %s [-n spp][-h\|-q][{+\|-}for[ward]][{+\|-}b[ackward]][-e expr][-f file] bsdf_func > bsdf.xml\n",
527			progname);
528			return(1);
529			}