src/gen/mkillum4.c

#ifndef lint
static const char RCSid[] = "$Id: mkillum4.c,v 2.5 2007/12/05 20:07:34 greg Exp $";
#endif
/*
 * Routines for handling BSDF data within mkillum
 */

#include "mkillum.h"
#include "paths.h"
#include "ezxml.h"

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

/* BSDF angle specification (terminate with nphi = -1) */
typedef struct {
        char    name[32];               /* 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 NABASES         3               /* number of predefined bases */

ANGLE_BASIS     abase_list[NABASES] = {
        {
                "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
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;
        d = sin(alt);
        v[0] = cos(azi)*d;
        v[1] = sin(azi)*d;
        v[2] = cos(alt);
        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] < 0.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;
        tdia = PI/180.*(ab->lat[li+1].tmin - ab->lat[li].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];

        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
)
{
        const char  *cbasis = ezxml_txt(ezxml_child(wdb,"ColumnAngleBasis"));
        const char  *rbasis = ezxml_txt(ezxml_child(wdb,"RowAngleBasis"));
        int  i;
        
        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_getvec;
                        dp->ib_ndx = ab_getndx;
                        dp->ib_ohm = ab_getohm;
                        break;
                }
        if (i < 0) {
                sprintf(errmsg, "unsupported incident basis '%s'", cbasis);
                error(INTERNAL, errmsg);
        }
        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_getvecR;
                        dp->ob_ndx = ab_getndxR;
                        dp->ob_ohm = ab_getohm;
                        break;
                }
        if (i < 0) {
                sprintf(errmsg, "unsupported exitant basis '%s'", cbasis);
                error(INTERNAL, errmsg);
        }
}


struct BSDF_data *
load_BSDF(              /* load BSDF data from file */
        char *fname
)
{
        char                    *path;
        ezxml_t                 fl, 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);
        }
        dp = (struct BSDF_data *)calloc(1, sizeof(struct BSDF_data));
        for (wld = ezxml_child(fl, "WavelengthData");
                                fl != NULL; fl = fl->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);
        normalize(v);
        return(1);
}


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);
        normalize(v);
        return(1);
}


#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;

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