src/cv/bsdf2klems.c

#ifndef lint
static const char RCSid[] = "$Id: bsdf2klems.c,v 2.12 2013/11/21 20:09:06 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)
{
        char    pbar[256];
        int     nchars;

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

/* Finish progress bar */
#define prog_done()     if (do_prog) fputc('\n',stderr); else

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