src/gen/mkillum4.c

#ifndef lint
static const char RCSid[] = "$Id: mkillum4.c,v 2.17 2009/05/30 22:19:08 greg Exp $";
#endif
/*
 * Routines for handling BSDF data within mkillum
 */

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

#ifndef NBSDFSAMPS
#define NBSDFSAMPS      32              /* BSDF resampling count */
#endif
#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  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_arr;
        double          dom, hemi_total;
        int             nneg;
        int             i, o;

        if (dp == NULL || dp->bsdf == NULL)
                return(0);
        omega_arr = (double *)calloc(dp->nout, sizeof(double));
        if (omega_arr == NULL)
                error(SYSTEM, "out of memory in check_bsdf_data");
        hemi_total = .0;
        for (o = dp->nout; o--; ) {
                FVECT   v;
                dom = getBSDF_outohm(dp,o);
                if (dom <= .0) {
                        error(WARNING, "zero/negative solid angle");
                        continue;
                }
                if (!getBSDF_outvec(v,dp,o) || v[2] < -FTINY) {
                        error(WARNING, "illegal outgoing BSDF direction");
                        free(omega_arr);
                        return(0);
                }
                hemi_total += omega_arr[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;          /* normalize solid angles */
        for (o = dp->nout; o--; )
                omega_arr[o] *= dom;
        nneg = 0;
        for (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_arr[o];
                        else if (f < -FTINY)
                                ++nneg;
                }
                if (hemi_total > 1.02) {
                        sprintf(errmsg, "BSDF direction passes %.1f%% of light",
                                        100.*hemi_total);
                        error(WARNING, errmsg);
                }
        }
        free(omega_arr);
        if (nneg > 0) {
                sprintf(errmsg, "%d negative BSDF values", nneg);
                error(WARNING, errmsg);
                return(0);
        }
        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);
}


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