src/common/bsdf.c

#ifndef lint
static const char RCSid[] = "$Id: bsdf.c,v 2.1 2009/06/17 20:41:47 greg Exp $";
#endif
/*
 * Routines for handling BSDF data
 */

#include "standard.h"
#include "bsdf.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  pol, azi, d;
        
        if ((ndx < 0) | (ndx >= ab->nangles))
                return(0);
        for (li = 0; ndx >= ab->lat[li].nphis; li++)
                ndx -= ab->lat[li].nphis;
        pol = 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(pol);
        d = sqrt(1. - d*d);     /* sin(pol) */
        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  pol, azi, d;

        if ((v[2] < -1.0) | (v[2] > 1.0))
                return(-1);
        pol = 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 <= pol; 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  theta, theta1;
        
        if ((ndx < 0) | (ndx >= ab->nangles))
                return(0);
        for (li = 0; ndx >= ab->lat[li].nphis; li++)
                ndx -= ab->lat[li].nphis;
        theta1 = PI/180. * ab->lat[li+1].tmin;
        if (ab->lat[li].nphis == 1) {           /* special case */
                if (ab->lat[li].tmin > FTINY)
                        error(USER, "unsupported BSDF coordinate system");
                return(2.*PI*(1. - cos(theta1)));
        }
        theta = PI/180. * ab->lat[li].tmin;
        return(2.*PI*(cos(theta) - cos(theta1))/(double)ab->lat[li].nphis);
}


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);
        }
}


static int
check_bsdf_data(        /* check that BSDF data is sane */
        struct BSDF_data *dp
)
{
        double          *omega_iarr, *omega_oarr;
        double          dom, contrib, hemi_total;
        int             nneg;
        FVECT           v;
        int             i, o;

        if (dp == NULL || dp->bsdf == NULL)
                return(0);
        omega_iarr = (double *)calloc(dp->ninc, sizeof(double));
        omega_oarr = (double *)calloc(dp->nout, sizeof(double));
        if ((omega_iarr == NULL) | (omega_oarr == NULL))
                error(SYSTEM, "out of memory in check_bsdf_data");
                                        /* incoming projected solid angles */
        hemi_total = .0;
        for (i = dp->ninc; i--; ) {
                dom = getBSDF_incohm(dp,i);
                if (dom <= .0) {
                        error(WARNING, "zero/negative incoming solid angle");
                        continue;
                }
                if (!getBSDF_incvec(v,dp,i) || v[2] > FTINY) {
                        error(WARNING, "illegal incoming BSDF direction");
                        free(omega_iarr); free(omega_oarr);
                        return(0);
                }
                hemi_total += omega_iarr[i] = dom * -v[2];
        }
        if ((hemi_total > 1.02*PI) | (hemi_total < 0.98*PI)) {
                sprintf(errmsg, "incoming BSDF hemisphere off by %.1f%%",
                                100.*(hemi_total/PI - 1.));
                error(WARNING, errmsg);
        }
        dom = PI / hemi_total;          /* fix normalization */
        for (i = dp->ninc; i--; )
                omega_iarr[i] *= dom;
                                        /* outgoing projected solid angles */
        hemi_total = .0;
        for (o = dp->nout; o--; ) {
                dom = getBSDF_outohm(dp,o);
                if (dom <= .0) {
                        error(WARNING, "zero/negative outgoing solid angle");
                        continue;
                }
                if (!getBSDF_outvec(v,dp,o) || v[2] < -FTINY) {
                        error(WARNING, "illegal outgoing BSDF direction");
                        free(omega_iarr); free(omega_oarr);
                        return(0);
                }
                hemi_total += omega_oarr[o] = dom * v[2];
        }
        if ((hemi_total > 1.02*PI) | (hemi_total < 0.98*PI)) {
                sprintf(errmsg, "outgoing BSDF hemisphere off by %.1f%%",
                                100.*(hemi_total/PI - 1.));
                error(WARNING, errmsg);
        }
        dom = PI / hemi_total;          /* fix normalization */
        for (o = dp->nout; o--; )
                omega_oarr[o] *= dom;
        nneg = 0;                       /* check outgoing totals */
        for (i = 0; i < dp->ninc; i++) {
                hemi_total = .0;
                for (o = dp->nout; o--; ) {
                        double  f = BSDF_value(dp,i,o);
                        if (f >= .0)
                                hemi_total += f*omega_oarr[o];
                        else {
                                nneg += (f < -FTINY);
                                BSDF_value(dp,i,o) = .0f;
                        }
                }
                if (hemi_total > 1.02) {
                        sprintf(errmsg,
                        "incoming BSDF direction %d passes %.1f%% of light",
                                        i, 100.*hemi_total);
                        error(WARNING, errmsg);
                }
        }
        if (nneg) {
                sprintf(errmsg, "%d negative BSDF values (ignored)", nneg);
                error(WARNING, errmsg);
        }
                                        /* reverse roles and check again */
        for (o = 0; o < dp->nout; o++) {
                hemi_total = .0;
                for (i = dp->ninc; i--; )
                        hemi_total += BSDF_value(dp,i,o) * omega_iarr[i];

                if (hemi_total > 1.02) {
                        sprintf(errmsg,
                        "outgoing BSDF direction %d collects %.1f%% of light",
                                        o, 100.*hemi_total);
                        error(WARNING, errmsg);
                }
        }
        free(omega_iarr); free(omega_oarr);
        return(1);
}

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