src/gen/mkillum4.c

#ifndef lint
static const char RCSid[] = "$Id: mkillum4.c,v 2.13 2008/11/10 19:08:18 greg Exp $";
#endif
/*
 * Routines for handling BSDF data within mkillum
 */

#include "mkillum.h"
#include "paths.h"
#include "ezxml.h"
#include <ctype.h>

#define MAXLATS         46              /* maximum number of latitudes */

/* BSDF angle specification */
typedef struct {
        char    name[64];               /* basis name */
        int     nangles;                /* total number of directions */
        struct {
                float   tmin;                   /* starting theta */
                short   nphis;                  /* number of phis (0 term) */
        }       lat[MAXLATS+1];         /* latitudes */
} ANGLE_BASIS;

#define MAXABASES       3               /* limit on defined bases */

static ANGLE_BASIS      abase_list[MAXABASES] = {
        {
                "LBNL/Klems Full", 145,
                { {-5., 1},
                {5., 8},
                {15., 16},
                {25., 20},
                {35., 24},
                {45., 24},
                {55., 24},
                {65., 16},
                {75., 12},
                {90., 0} }
        }, {
                "LBNL/Klems Half", 73,
                { {-6.5, 1},
                {6.5, 8},
                {19.5, 12},
                {32.5, 16},
                {46.5, 20},
                {61.5, 12},
                {76.5, 4},
                {90., 0} }
        }, {
                "LBNL/Klems Quarter", 41,
                { {-9., 1},
                {9., 8},
                {27., 12},
                {46., 12},
                {66., 8},
                {90., 0} }
        }
};

static int      nabases = 3;    /* current number of defined bases */


static int
ab_getvec(              /* get vector for this angle basis index */
        FVECT v,
        int ndx,
        void *p
)
{
        ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
        int     li;
        double  alt, azi, d;
        
        if ((ndx < 0) | (ndx >= ab->nangles))
                return(0);
        for (li = 0; ndx >= ab->lat[li].nphis; li++)
                ndx -= ab->lat[li].nphis;
        alt = PI/180.*0.5*(ab->lat[li].tmin + ab->lat[li+1].tmin);
        azi = 2.*PI*ndx/ab->lat[li].nphis;
        v[2] = d = cos(alt);
        d = sqrt(1. - d*d);     /* sin(alt) */
        v[0] = cos(azi)*d;
        v[1] = sin(azi)*d;
        return(1);
}


static int
ab_getndx(              /* get index corresponding to the given vector */
        FVECT v,
        void *p
)
{
        ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
        int     li, ndx;
        double  alt, azi, d;

        if ((v[2] < -1.0) | (v[2] > 1.0))
                return(-1);
        alt = 180.0/PI*acos(v[2]);
        azi = 180.0/PI*atan2(v[1], v[0]);
        if (azi < 0.0) azi += 360.0;
        for (li = 1; ab->lat[li].tmin <= alt; li++)
                if (!ab->lat[li].nphis)
                        return(-1);
        --li;
        ndx = (int)((1./360.)*azi*ab->lat[li].nphis + 0.5);
        if (ndx >= ab->lat[li].nphis) ndx = 0;
        while (li--)
                ndx += ab->lat[li].nphis;
        return(ndx);
}


static double
ab_getohm(              /* get solid angle for this angle basis index */
        int ndx,
        void *p
)
{
        ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
        int     li;
        double  tdia, pdia;
        
        if ((ndx < 0) | (ndx >= ab->nangles))
                return(0);
        for (li = 0; ndx >= ab->lat[li].nphis; li++)
                ndx -= ab->lat[li].nphis;
        if (ab->lat[li].nphis == 1) {           /* special case */
                if (ab->lat[li].tmin > FTINY)
                        error(USER, "unsupported BSDF coordinate system");
                tdia = PI/180. * ab->lat[li+1].tmin;
                return(PI*tdia*tdia);
        }
        tdia = PI/180.*(ab->lat[li+1].tmin - ab->lat[li].tmin);
        tdia *= sin(PI/180.*(ab->lat[li].tmin + ab->lat[li+1].tmin));
        pdia = 2.*PI/ab->lat[li].nphis;
        return(tdia*pdia);
}


static int
ab_getvecR(             /* get reverse vector for this angle basis index */
        FVECT v,
        int ndx,
        void *p
)
{
        if (!ab_getvec(v, ndx, p))
                return(0);

        v[0] = -v[0];
        v[1] = -v[1];
        v[2] = -v[2];

        return(1);
}


static int
ab_getndxR(             /* get index corresponding to the reverse vector */
        FVECT v,
        void *p
)
{
        FVECT  v2;
        
        v2[0] = -v[0];
        v2[1] = -v[1];
        v2[2] = -v[2];

        return ab_getndx(v2, p);
}


static void
load_bsdf_data(         /* load BSDF distribution for this wavelength */
        struct BSDF_data *dp,
        ezxml_t wdb
)
{
        char  *cbasis = ezxml_txt(ezxml_child(wdb,"ColumnAngleBasis"));
        char  *rbasis = ezxml_txt(ezxml_child(wdb,"RowAngleBasis"));
        char  *sdata;
        int  i;
        
        if ((cbasis == NULL) | (rbasis == NULL)) {
                error(WARNING, "missing column/row basis for BSDF");
                return;
        }
        /* XXX need to add routines for loading in foreign bases */
        for (i = nabases; i--; )
                if (!strcmp(cbasis, abase_list[i].name)) {
                        dp->ninc = abase_list[i].nangles;
                        dp->ib_priv = (void *)&abase_list[i];
                        dp->ib_vec = ab_getvecR;
                        dp->ib_ndx = ab_getndxR;
                        dp->ib_ohm = ab_getohm;
                        break;
                }
        if (i < 0) {
                sprintf(errmsg, "unsupported ColumnAngleBasis '%s'", cbasis);
                error(WARNING, errmsg);
                return;
        }
        for (i = nabases; i--; )
                if (!strcmp(rbasis, abase_list[i].name)) {
                        dp->nout = abase_list[i].nangles;
                        dp->ob_priv = (void *)&abase_list[i];
                        dp->ob_vec = ab_getvec;
                        dp->ob_ndx = ab_getndx;
                        dp->ob_ohm = ab_getohm;
                        break;
                }
        if (i < 0) {
                sprintf(errmsg, "unsupported RowAngleBasis '%s'", cbasis);
                error(WARNING, errmsg);
                return;
        }
                                /* read BSDF data */
        sdata  = ezxml_txt(ezxml_child(wdb,"ScatteringData"));
        if (sdata == NULL) {
                error(WARNING, "missing BSDF ScatteringData");
                return;
        }
        dp->bsdf = (float *)malloc(sizeof(float)*dp->ninc*dp->nout);
        if (dp->bsdf == NULL)
                error(SYSTEM, "out of memory in load_bsdf_data");
        for (i = 0; i < dp->ninc*dp->nout; i++) {
                char  *sdnext = fskip(sdata);
                if (sdnext == NULL) {
                        error(WARNING, "bad/missing BSDF ScatteringData");
                        free(dp->bsdf); dp->bsdf = NULL;
                        return;
                }
                while (*sdnext && isspace(*sdnext))
                        sdnext++;
                if (*sdnext == ',') sdnext++;
                dp->bsdf[i] = atof(sdata);
                sdata = sdnext;
        }
        while (isspace(*sdata))
                sdata++;
        if (*sdata) {
                sprintf(errmsg, "%d extra characters after BSDF ScatteringData",
                                strlen(sdata));
                error(WARNING, errmsg);
        }
}


struct BSDF_data *
load_BSDF(              /* load BSDF data from file */
        char *fname
)
{
        char                    *path;
        ezxml_t                 fl, wtl, wld, wdb;
        struct BSDF_data        *dp;
        
        path = getpath(fname, getrlibpath(), R_OK);
        if (path == NULL) {
                sprintf(errmsg, "cannot find BSDF file \"%s\"", fname);
                error(WARNING, errmsg);
                return(NULL);
        }
        fl = ezxml_parse_file(path);
        if (fl == NULL) {
                sprintf(errmsg, "cannot open BSDF \"%s\"", path);
                error(WARNING, errmsg);
                return(NULL);
        }
        if (ezxml_error(fl)[0]) {
                sprintf(errmsg, "BSDF \"%s\" %s", path, ezxml_error(fl));
                error(WARNING, errmsg);
                ezxml_free(fl);
                return(NULL);
        }
        if (strcmp(ezxml_name(fl), "WindowElement")) {
                sprintf(errmsg,
                        "BSDF \"%s\": top level node not 'WindowElement'",
                                path);
                error(WARNING, errmsg);
                ezxml_free(fl);
                return(NULL);
        }
        wtl = ezxml_child(ezxml_child(fl, "Optical"), "Layer");
        dp = (struct BSDF_data *)calloc(1, sizeof(struct BSDF_data));
        for (wld = ezxml_child(wtl, "WavelengthData");
                                wld != NULL; wld = wld->next) {
                if (strcmp(ezxml_txt(ezxml_child(wld,"Wavelength")), "Visible"))
                        continue;
                wdb = ezxml_child(wld, "WavelengthDataBlock");
                if (wdb == NULL) continue;
                if (strcmp(ezxml_txt(ezxml_child(wdb,"WavelengthDataDirection")),
                                        "Transmission Front"))
                        continue;
                load_bsdf_data(dp, wdb);        /* load front BTDF */
                break;                          /* ignore the rest */
        }
        ezxml_free(fl);                         /* done with XML file */
        if (dp->bsdf == NULL) {
                sprintf(errmsg, "bad/missing BTDF data in \"%s\"", path);
                error(WARNING, errmsg);
                free_BSDF(dp);
                dp = NULL;
        }
        return(dp);
}


void
free_BSDF(              /* free BSDF data structure */
        struct BSDF_data *b
)
{
        if (b == NULL)
                return;
        if (b->bsdf != NULL)
                free(b->bsdf);
        free(b);
}


int
r_BSDF_incvec(          /* compute random input vector at given location */
        FVECT v,
        struct BSDF_data *b,
        int i,
        double rv,
        MAT4 xm
)
{
        FVECT   pert;
        double  rad;
        int     j;
        
        if (!getBSDF_incvec(v, b, i))
                return(0);
        rad = sqrt(getBSDF_incohm(b, i) / PI);
        multisamp(pert, 3, rv);
        for (j = 0; j < 3; j++)
                v[j] += rad*(2.*pert[j] - 1.);
        if (xm != NULL)
                multv3(v, v, xm);
        return(normalize(v) != 0.0);
}


int
r_BSDF_outvec(          /* compute random output vector at given location */
        FVECT v,
        struct BSDF_data *b,
        int o,
        double rv,
        MAT4 xm
)
{
        FVECT   pert;
        double  rad;
        int     j;
        
        if (!getBSDF_outvec(v, b, o))
                return(0);
        rad = sqrt(getBSDF_outohm(b, o) / PI);
        multisamp(pert, 3, rv);
        for (j = 0; j < 3; j++)
                v[j] += rad*(2.*pert[j] - 1.);
        if (xm != NULL)
                multv3(v, v, xm);
        return(normalize(v) != 0.0);
}


#define  FEQ(a,b)       ((a)-(b) <= 1e-7 && (b)-(a) <= 1e-7)

static int
addrot(                 /* compute rotation (x,y,z) => (xp,yp,zp) */
        char *xfarg[],
        FVECT xp,
        FVECT yp,
        FVECT zp
)
{
        static char     bufs[3][16];
        int     bn = 0;
        char    **xfp = xfarg;
        double  theta;

        if (yp[2]*yp[2] + zp[2]*zp[2] < 2.*FTINY*FTINY) {
                /* Special case for X' along Z-axis */
                theta = -atan2(yp[0], yp[1]);
                *xfp++ = "-ry";
                *xfp++ = xp[2] < 0.0 ? "90" : "-90";
                *xfp++ = "-rz";
                sprintf(bufs[bn], "%f", theta*(180./PI));
                *xfp++ = bufs[bn++];
                return(xfp - xfarg);
        }
        theta = atan2(yp[2], zp[2]);
        if (!FEQ(theta,0.0)) {
                *xfp++ = "-rx";
                sprintf(bufs[bn], "%f", theta*(180./PI));
                *xfp++ = bufs[bn++];
        }
        theta = asin(-xp[2]);
        if (!FEQ(theta,0.0)) {
                *xfp++ = "-ry";
                sprintf(bufs[bn], " %f", theta*(180./PI));
                *xfp++ = bufs[bn++];
        }
        theta = atan2(xp[1], xp[0]);
        if (!FEQ(theta,0.0)) {
                *xfp++ = "-rz";
                sprintf(bufs[bn], "%f", theta*(180./PI));
                *xfp++ = bufs[bn++];
        }
        *xfp = NULL;
        return(xfp - xfarg);
}


int
getBSDF_xfm(            /* compute BSDF orient. -> world orient. transform */
        MAT4 xm,
        FVECT nrm,
        UpDir ud
)
{
        char    *xfargs[7];
        XF      myxf;
        FVECT   updir, xdest, ydest;

        updir[0] = updir[1] = updir[2] = 0.;
        switch (ud) {
        case UDzneg:
                updir[2] = -1.;
                break;
        case UDyneg:
                updir[1] = -1.;
                break;
        case UDxneg:
                updir[0] = -1.;
                break;
        case UDxpos:
                updir[0] = 1.;
                break;
        case UDypos:
                updir[1] = 1.;
                break;
        case UDzpos:
                updir[2] = 1.;
                break;
        case UDunknown:
                return(0);
        }
        fcross(xdest, updir, nrm);
        if (normalize(xdest) == 0.0)
                return(0);
        fcross(ydest, nrm, xdest);
        xf(&myxf, addrot(xfargs, xdest, ydest, nrm), xfargs);
        copymat4(xm, myxf.xfm);
        return(1);
}


void
redistribute(           /* pass distarr ray sums through BSDF */
        struct BSDF_data *b,
        int nalt,
        int nazi,
        FVECT u,
        FVECT v,
        FVECT w,
        MAT4 xm
)
{
        int     nout = 0;
        MAT4    mymat, inmat;
        COLORV  *idist;
        COLORV  *cp, *csum;
        FVECT   dv;
        double  wt;
        int     i, j, k, o;
        COLOR   col, cinc;
                                        /* copy incoming distribution */
        if (b->ninc > distsiz)
                error(INTERNAL, "error 1 in redistribute");
        idist = (COLORV *)malloc(sizeof(COLOR)*b->ninc);
        if (idist == NULL)
                error(SYSTEM, "out of memory in redistribute");
        memcpy(idist, distarr, sizeof(COLOR)*b->ninc);
                                        /* compose direction transform */
        for (i = 3; i--; ) {
                mymat[i][0] = u[i];
                mymat[i][1] = v[i];
                mymat[i][2] = w[i];
                mymat[i][3] = 0.;
        }
        mymat[3][0] = mymat[3][1] = mymat[3][2] = 0.;
        mymat[3][3] = 1.;
        if (xm != NULL)
                multmat4(mymat, xm, mymat);
        for (i = 3; i--; ) {            /* make sure it's normalized */
                wt = 1./sqrt(   mymat[0][i]*mymat[0][i] +
                                mymat[1][i]*mymat[1][i] +
                                mymat[2][i]*mymat[2][i] );
                for (j = 3; j--; )
                        mymat[j][i] *= wt;
        }
        if (!invmat4(inmat, mymat))     /* need inverse as well */
                error(INTERNAL, "cannot invert BSDF transform");
        newdist(nalt*nazi);             /* resample distribution */
        for (i = b->ninc; i--; ) {
                getBSDF_incvec(dv, b, i);       /* compute incident irrad. */
                multv3(dv, dv, mymat);
                if (dv[2] < 0.0) dv[2] = -dv[2];
                wt = getBSDF_incohm(b, i);
                wt *= dv[2];                    /* solid_angle*cosine(theta) */
                cp = &idist[3*i];
                copycolor(cinc, cp);
                scalecolor(cinc, wt);
                for (k = nalt; k--; )           /* loop over distribution */
                    for (j = nazi; j--; ) {
                        flatdir(dv, (k + .5)/nalt, (double)j/nazi);
                        multv3(dv, dv, inmat);
                                                /* evaluate BSDF @ outgoing */
                        o = getBSDF_outndx(b, dv);
                        if (o < 0) {
                                nout++;
                                continue;
                        }
                        wt = BSDF_value(b, i, o);
                        copycolor(col, cinc);
                        scalecolor(col, wt);
                        csum = &distarr[3*(k*nazi + j)];
                        addcolor(csum, col);    /* sum into distribution */
                    }
        }
        free(idist);                    /* free temp space */
        if (nout) {
                sprintf(errmsg, "missing %.1f%% of BSDF directions",
                                100.*nout/(b->ninc*nalt*nazi));
                error(WARNING, errmsg);
        }
}
Revision:	2.14
Committed:	Sun Jan 18 17:41:38 2009 UTC (16 years, 3 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.13:	+4 -4 lines
Log Message:	Bug fix for "error 1 in redistribute"
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.14	static const char RCSid[] = "$Id: mkillum4.c,v 2.13 2008/11/10 19:08:18 greg Exp $";
3	greg	2.1	#endif
4			/*
5			* Routines for handling BSDF data within mkillum
6			*/
7
8			#include "mkillum.h"
9			#include "paths.h"
10	greg	2.2	#include "ezxml.h"
11	greg	2.12	#include <ctype.h>
12	greg	2.1
13	greg	2.6	#define MAXLATS 46 /* maximum number of latitudes */
14
15	greg	2.11	/* BSDF angle specification */
16	greg	2.6	typedef struct {
17	greg	2.7	char name[64]; /* basis name */
18	greg	2.6	int nangles; /* total number of directions */
19			struct {
20			float tmin; /* starting theta */
21			short nphis; /* number of phis (0 term) */
22			} lat[MAXLATS+1]; /* latitudes */
23			} ANGLE_BASIS;
24
25	greg	2.7	#define MAXABASES 3 /* limit on defined bases */
26	greg	2.6
27	greg	2.13	static ANGLE_BASIS abase_list[MAXABASES] = {
28	greg	2.6	{
29			"LBNL/Klems Full", 145,
30			{ {-5., 1},
31			{5., 8},
32			{15., 16},
33			{25., 20},
34			{35., 24},
35			{45., 24},
36			{55., 24},
37			{65., 16},
38			{75., 12},
39			{90., 0} }
40			}, {
41			"LBNL/Klems Half", 73,
42			{ {-6.5, 1},
43			{6.5, 8},
44			{19.5, 12},
45			{32.5, 16},
46			{46.5, 20},
47			{61.5, 12},
48			{76.5, 4},
49			{90., 0} }
50			}, {
51			"LBNL/Klems Quarter", 41,
52			{ {-9., 1},
53			{9., 8},
54			{27., 12},
55			{46., 12},
56			{66., 8},
57			{90., 0} }
58			}
59			};
60
61	greg	2.7	static int nabases = 3; /* current number of defined bases */
62
63	greg	2.6
64			static int
65			ab_getvec( /* get vector for this angle basis index */
66			FVECT v,
67			int ndx,
68			void *p
69			)
70			{
71			ANGLE_BASIS ab = (ANGLE_BASIS )p;
72			int li;
73			double alt, azi, d;
74
75			if ((ndx < 0) \| (ndx >= ab->nangles))
76			return(0);
77			for (li = 0; ndx >= ab->lat[li].nphis; li++)
78			ndx -= ab->lat[li].nphis;
79			alt = PI/180.0.5(ab->lat[li].tmin + ab->lat[li+1].tmin);
80			azi = 2.PIndx/ab->lat[li].nphis;
81	greg	2.11	v[2] = d = cos(alt);
82			d = sqrt(1. - dd); / sin(alt) */
83	greg	2.6	v[0] = cos(azi)*d;
84			v[1] = sin(azi)*d;
85			return(1);
86			}
87
88
89			static int
90			ab_getndx( /* get index corresponding to the given vector */
91			FVECT v,
92			void *p
93			)
94			{
95			ANGLE_BASIS ab = (ANGLE_BASIS )p;
96			int li, ndx;
97			double alt, azi, d;
98
99	greg	2.7	if ((v[2] < -1.0) \| (v[2] > 1.0))
100	greg	2.6	return(-1);
101			alt = 180.0/PI*acos(v[2]);
102			azi = 180.0/PI*atan2(v[1], v[0]);
103			if (azi < 0.0) azi += 360.0;
104			for (li = 1; ab->lat[li].tmin <= alt; li++)
105			if (!ab->lat[li].nphis)
106			return(-1);
107			--li;
108			ndx = (int)((1./360.)aziab->lat[li].nphis + 0.5);
109			if (ndx >= ab->lat[li].nphis) ndx = 0;
110			while (li--)
111			ndx += ab->lat[li].nphis;
112			return(ndx);
113			}
114
115
116			static double
117			ab_getohm( /* get solid angle for this angle basis index */
118			int ndx,
119			void *p
120			)
121			{
122			ANGLE_BASIS ab = (ANGLE_BASIS )p;
123			int li;
124			double tdia, pdia;
125
126			if ((ndx < 0) \| (ndx >= ab->nangles))
127			return(0);
128			for (li = 0; ndx >= ab->lat[li].nphis; li++)
129			ndx -= ab->lat[li].nphis;
130	greg	2.7	if (ab->lat[li].nphis == 1) { /* special case */
131			if (ab->lat[li].tmin > FTINY)
132			error(USER, "unsupported BSDF coordinate system");
133			tdia = PI/180. * ab->lat[li+1].tmin;
134			return(PItdiatdia);
135			}
136	greg	2.6	tdia = PI/180.*(ab->lat[li+1].tmin - ab->lat[li].tmin);
137	greg	2.7	tdia = sin(PI/180.(ab->lat[li].tmin + ab->lat[li+1].tmin));
138	greg	2.6	pdia = 2.*PI/ab->lat[li].nphis;
139			return(tdia*pdia);
140			}
141
142
143			static int
144			ab_getvecR( /* get reverse vector for this angle basis index */
145			FVECT v,
146			int ndx,
147			void *p
148			)
149			{
150			if (!ab_getvec(v, ndx, p))
151			return(0);
152
153			v[0] = -v[0];
154			v[1] = -v[1];
155	greg	2.11	v[2] = -v[2];
156	greg	2.6
157			return(1);
158			}
159
160
161			static int
162			ab_getndxR( /* get index corresponding to the reverse vector */
163			FVECT v,
164			void *p
165			)
166			{
167			FVECT v2;
168
169			v2[0] = -v[0];
170			v2[1] = -v[1];
171	greg	2.11	v2[2] = -v[2];
172	greg	2.6
173			return ab_getndx(v2, p);
174			}
175
176	greg	2.7
177	greg	2.6	static void
178			load_bsdf_data( /* load BSDF distribution for this wavelength */
179			struct BSDF_data *dp,
180			ezxml_t wdb
181			)
182			{
183	greg	2.7	char *cbasis = ezxml_txt(ezxml_child(wdb,"ColumnAngleBasis"));
184			char *rbasis = ezxml_txt(ezxml_child(wdb,"RowAngleBasis"));
185			char *sdata;
186	greg	2.6	int i;
187
188	greg	2.7	if ((cbasis == NULL) \| (rbasis == NULL)) {
189			error(WARNING, "missing column/row basis for BSDF");
190			return;
191			}
192			/* XXX need to add routines for loading in foreign bases */
193			for (i = nabases; i--; )
194	greg	2.6	if (!strcmp(cbasis, abase_list[i].name)) {
195			dp->ninc = abase_list[i].nangles;
196			dp->ib_priv = (void *)&abase_list[i];
197	greg	2.11	dp->ib_vec = ab_getvecR;
198			dp->ib_ndx = ab_getndxR;
199	greg	2.6	dp->ib_ohm = ab_getohm;
200			break;
201			}
202			if (i < 0) {
203	greg	2.7	sprintf(errmsg, "unsupported ColumnAngleBasis '%s'", cbasis);
204			error(WARNING, errmsg);
205			return;
206	greg	2.6	}
207	greg	2.7	for (i = nabases; i--; )
208	greg	2.6	if (!strcmp(rbasis, abase_list[i].name)) {
209			dp->nout = abase_list[i].nangles;
210			dp->ob_priv = (void *)&abase_list[i];
211	greg	2.11	dp->ob_vec = ab_getvec;
212			dp->ob_ndx = ab_getndx;
213	greg	2.6	dp->ob_ohm = ab_getohm;
214			break;
215			}
216			if (i < 0) {
217	greg	2.7	sprintf(errmsg, "unsupported RowAngleBasis '%s'", cbasis);
218			error(WARNING, errmsg);
219			return;
220			}
221			/* read BSDF data */
222			sdata = ezxml_txt(ezxml_child(wdb,"ScatteringData"));
223			if (sdata == NULL) {
224			error(WARNING, "missing BSDF ScatteringData");
225			return;
226			}
227			dp->bsdf = (float )malloc(sizeof(float)dp->ninc*dp->nout);
228			if (dp->bsdf == NULL)
229			error(SYSTEM, "out of memory in load_bsdf_data");
230			for (i = 0; i < dp->ninc*dp->nout; i++) {
231			char *sdnext = fskip(sdata);
232			if (sdnext == NULL) {
233			error(WARNING, "bad/missing BSDF ScatteringData");
234			free(dp->bsdf); dp->bsdf = NULL;
235			return;
236			}
237			while (sdnext && isspace(sdnext))
238			sdnext++;
239			if (*sdnext == ',') sdnext++;
240			dp->bsdf[i] = atof(sdata);
241			sdata = sdnext;
242			}
243			while (isspace(*sdata))
244			sdata++;
245			if (*sdata) {
246			sprintf(errmsg, "%d extra characters after BSDF ScatteringData",
247			strlen(sdata));
248			error(WARNING, errmsg);
249	greg	2.6	}
250			}
251
252	greg	2.1
253			struct BSDF_data *
254			load_BSDF( /* load BSDF data from file */
255			char *fname
256			)
257			{
258			char *path;
259	greg	2.10	ezxml_t fl, wtl, wld, wdb;
260	greg	2.1	struct BSDF_data *dp;
261
262			path = getpath(fname, getrlibpath(), R_OK);
263			if (path == NULL) {
264			sprintf(errmsg, "cannot find BSDF file \"%s\"", fname);
265			error(WARNING, errmsg);
266			return(NULL);
267			}
268	greg	2.2	fl = ezxml_parse_file(path);
269			if (fl == NULL) {
270	greg	2.1	sprintf(errmsg, "cannot open BSDF \"%s\"", path);
271			error(WARNING, errmsg);
272			return(NULL);
273			}
274	greg	2.6	if (ezxml_error(fl)[0]) {
275	greg	2.10	sprintf(errmsg, "BSDF \"%s\" %s", path, ezxml_error(fl));
276	greg	2.6	error(WARNING, errmsg);
277			ezxml_free(fl);
278			return(NULL);
279			}
280	greg	2.10	if (strcmp(ezxml_name(fl), "WindowElement")) {
281			sprintf(errmsg,
282			"BSDF \"%s\": top level node not 'WindowElement'",
283			path);
284			error(WARNING, errmsg);
285			ezxml_free(fl);
286			return(NULL);
287			}
288			wtl = ezxml_child(ezxml_child(fl, "Optical"), "Layer");
289	greg	2.5	dp = (struct BSDF_data *)calloc(1, sizeof(struct BSDF_data));
290	greg	2.10	for (wld = ezxml_child(wtl, "WavelengthData");
291	greg	2.8	wld != NULL; wld = wld->next) {
292	greg	2.6	if (strcmp(ezxml_txt(ezxml_child(wld,"Wavelength")), "Visible"))
293	greg	2.2	continue;
294			wdb = ezxml_child(wld, "WavelengthDataBlock");
295			if (wdb == NULL) continue;
296	greg	2.6	if (strcmp(ezxml_txt(ezxml_child(wdb,"WavelengthDataDirection")),
297	greg	2.2	"Transmission Front"))
298			continue;
299	greg	2.6	load_bsdf_data(dp, wdb); /* load front BTDF */
300			break; /* ignore the rest */
301			}
302			ezxml_free(fl); /* done with XML file */
303			if (dp->bsdf == NULL) {
304			sprintf(errmsg, "bad/missing BTDF data in \"%s\"", path);
305			error(WARNING, errmsg);
306			free_BSDF(dp);
307			dp = NULL;
308	greg	2.2	}
309	greg	2.1	return(dp);
310			}
311
312
313			void
314			free_BSDF( /* free BSDF data structure */
315			struct BSDF_data *b
316			)
317			{
318			if (b == NULL)
319			return;
320	greg	2.6	if (b->bsdf != NULL)
321			free(b->bsdf);
322	greg	2.1	free(b);
323			}
324
325
326	greg	2.6	int
327	greg	2.1	r_BSDF_incvec( /* compute random input vector at given location */
328			FVECT v,
329			struct BSDF_data *b,
330			int i,
331			double rv,
332			MAT4 xm
333			)
334			{
335			FVECT pert;
336			double rad;
337			int j;
338
339	greg	2.6	if (!getBSDF_incvec(v, b, i))
340			return(0);
341			rad = sqrt(getBSDF_incohm(b, i) / PI);
342	greg	2.1	multisamp(pert, 3, rv);
343			for (j = 0; j < 3; j++)
344			v[j] += rad(2.pert[j] - 1.);
345			if (xm != NULL)
346			multv3(v, v, xm);
347	greg	2.7	return(normalize(v) != 0.0);
348	greg	2.1	}
349
350
351	greg	2.6	int
352	greg	2.1	r_BSDF_outvec( /* compute random output vector at given location */
353			FVECT v,
354			struct BSDF_data *b,
355			int o,
356			double rv,
357			MAT4 xm
358			)
359			{
360			FVECT pert;
361			double rad;
362			int j;
363
364	greg	2.6	if (!getBSDF_outvec(v, b, o))
365			return(0);
366			rad = sqrt(getBSDF_outohm(b, o) / PI);
367	greg	2.1	multisamp(pert, 3, rv);
368			for (j = 0; j < 3; j++)
369			v[j] += rad(2.pert[j] - 1.);
370			if (xm != NULL)
371			multv3(v, v, xm);
372	greg	2.7	return(normalize(v) != 0.0);
373	greg	2.1	}
374	greg	2.2
375
376			#define FEQ(a,b) ((a)-(b) <= 1e-7 && (b)-(a) <= 1e-7)
377
378			static int
379			addrot( /* compute rotation (x,y,z) => (xp,yp,zp) */
380			char *xfarg[],
381			FVECT xp,
382			FVECT yp,
383			FVECT zp
384			)
385			{
386			static char bufs[3][16];
387			int bn = 0;
388			char **xfp = xfarg;
389			double theta;
390
391	greg	2.9	if (yp[2]yp[2] + zp[2]zp[2] < 2.FTINYFTINY) {
392			/* Special case for X' along Z-axis */
393			theta = -atan2(yp[0], yp[1]);
394			*xfp++ = "-ry";
395			*xfp++ = xp[2] < 0.0 ? "90" : "-90";
396			*xfp++ = "-rz";
397			sprintf(bufs[bn], "%f", theta*(180./PI));
398			*xfp++ = bufs[bn++];
399			return(xfp - xfarg);
400			}
401	greg	2.2	theta = atan2(yp[2], zp[2]);
402			if (!FEQ(theta,0.0)) {
403			*xfp++ = "-rx";
404			sprintf(bufs[bn], "%f", theta*(180./PI));
405			*xfp++ = bufs[bn++];
406			}
407			theta = asin(-xp[2]);
408			if (!FEQ(theta,0.0)) {
409			*xfp++ = "-ry";
410			sprintf(bufs[bn], " %f", theta*(180./PI));
411			*xfp++ = bufs[bn++];
412			}
413			theta = atan2(xp[1], xp[0]);
414			if (!FEQ(theta,0.0)) {
415			*xfp++ = "-rz";
416			sprintf(bufs[bn], "%f", theta*(180./PI));
417			*xfp++ = bufs[bn++];
418			}
419			*xfp = NULL;
420			return(xfp - xfarg);
421			}
422
423
424			int
425	greg	2.6	getBSDF_xfm( /* compute BSDF orient. -> world orient. transform */
426	greg	2.2	MAT4 xm,
427			FVECT nrm,
428			UpDir ud
429			)
430			{
431			char *xfargs[7];
432			XF myxf;
433			FVECT updir, xdest, ydest;
434
435			updir[0] = updir[1] = updir[2] = 0.;
436			switch (ud) {
437			case UDzneg:
438			updir[2] = -1.;
439			break;
440			case UDyneg:
441			updir[1] = -1.;
442			break;
443			case UDxneg:
444			updir[0] = -1.;
445			break;
446			case UDxpos:
447			updir[0] = 1.;
448			break;
449			case UDypos:
450			updir[1] = 1.;
451			break;
452			case UDzpos:
453			updir[2] = 1.;
454			break;
455			case UDunknown:
456			return(0);
457			}
458			fcross(xdest, updir, nrm);
459			if (normalize(xdest) == 0.0)
460			return(0);
461			fcross(ydest, nrm, xdest);
462			xf(&myxf, addrot(xfargs, xdest, ydest, nrm), xfargs);
463			copymat4(xm, myxf.xfm);
464			return(1);
465			}
466
467
468			void
469			redistribute( /* pass distarr ray sums through BSDF */
470			struct BSDF_data *b,
471			int nalt,
472			int nazi,
473			FVECT u,
474			FVECT v,
475			FVECT w,
476			MAT4 xm
477			)
478			{
479	greg	2.6	int nout = 0;
480	greg	2.5	MAT4 mymat, inmat;
481			COLORV *idist;
482	greg	2.2	COLORV cp, csum;
483			FVECT dv;
484	greg	2.5	double wt;
485	greg	2.6	int i, j, k, o;
486	greg	2.5	COLOR col, cinc;
487			/* copy incoming distribution */
488	greg	2.14	if (b->ninc > distsiz)
489	greg	2.5	error(INTERNAL, "error 1 in redistribute");
490	greg	2.14	idist = (COLORV )malloc(sizeof(COLOR)b->ninc);
491	greg	2.5	if (idist == NULL)
492	greg	2.2	error(SYSTEM, "out of memory in redistribute");
493	greg	2.14	memcpy(idist, distarr, sizeof(COLOR)*b->ninc);
494	greg	2.5	/* compose direction transform */
495	greg	2.2	for (i = 3; i--; ) {
496			mymat[i][0] = u[i];
497			mymat[i][1] = v[i];
498			mymat[i][2] = w[i];
499			mymat[i][3] = 0.;
500			}
501			mymat[3][0] = mymat[3][1] = mymat[3][2] = 0.;
502			mymat[3][3] = 1.;
503			if (xm != NULL)
504			multmat4(mymat, xm, mymat);
505			for (i = 3; i--; ) { /* make sure it's normalized */
506	greg	2.5	wt = 1./sqrt( mymat[0][i]*mymat[0][i] +
507	greg	2.2	mymat[1][i]*mymat[1][i] +
508			mymat[2][i]*mymat[2][i] );
509			for (j = 3; j--; )
510	greg	2.5	mymat[j][i] *= wt;
511	greg	2.2	}
512	greg	2.5	if (!invmat4(inmat, mymat)) /* need inverse as well */
513			error(INTERNAL, "cannot invert BSDF transform");
514			newdist(naltnazi); / resample distribution */
515	greg	2.4	for (i = b->ninc; i--; ) {
516	greg	2.5	getBSDF_incvec(dv, b, i); /* compute incident irrad. */
517	greg	2.2	multv3(dv, dv, mymat);
518	greg	2.5	if (dv[2] < 0.0) dv[2] = -dv[2];
519	greg	2.6	wt = getBSDF_incohm(b, i);
520			wt = dv[2]; / solid_anglecosine(theta) /
521	greg	2.5	cp = &idist[3*i];
522			copycolor(cinc, cp);
523			scalecolor(cinc, wt);
524			for (k = nalt; k--; ) /* loop over distribution */
525			for (j = nazi; j--; ) {
526			flatdir(dv, (k + .5)/nalt, (double)j/nazi);
527			multv3(dv, dv, inmat);
528			/* evaluate BSDF @ outgoing */
529	greg	2.6	o = getBSDF_outndx(b, dv);
530			if (o < 0) {
531			nout++;
532			continue;
533			}
534	greg	2.7	wt = BSDF_value(b, i, o);
535	greg	2.5	copycolor(col, cinc);
536			scalecolor(col, wt);
537			csum = &distarr[3(knazi + j)];
538			addcolor(csum, col); /* sum into distribution */
539			}
540	greg	2.2	}
541	greg	2.5	free(idist); /* free temp space */
542	greg	2.6	if (nout) {
543			sprintf(errmsg, "missing %.1f%% of BSDF directions",
544			100.nout/(b->nincnalt*nazi));
545			error(WARNING, errmsg);
546			}
547	greg	2.2	}