src/gen/mkillum4.c

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

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


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);
        }
        dp = (struct BSDF_data *)calloc(1, sizeof(struct BSDF_data));
        if (dp == NULL)
                goto memerr;
        for (wld = ezxml_child(fl, "WavelengthData");
                                fl != NULL; fl = fl->next) {
                if (strcmp(ezxml_child(wld, "Wavelength")->txt, "Visible"))
                        continue;
                wdb = ezxml_child(wld, "WavelengthDataBlock");
                if (wdb == NULL) continue;
                if (strcmp(ezxml_child(wdb, "WavelengthDataDirection")->txt,
                                        "Transmission Front"))
                        continue;
        }
        /* etc... */
        ezxml_free(fl);
        return(dp);
memerr:
        error(SYSTEM, "out of memory in load_BSDF");
        return NULL;    /* pro forma return */
}


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


void
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;
        
        getBSDF_incvec(v, b, i);
        rad = getBSDF_incrad(b, i);
        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);
}


void
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;
        
        getBSDF_outvec(v, b, o);
        rad = getBSDF_outrad(b, o);
        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);
}


#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 transform for the given surface */
        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
)
{
        MAT4    mymat, inmat;
        COLORV  *idist;
        COLORV  *cp, *csum;
        FVECT   dv;
        double  wt;
        int     i, j, k, h;
        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_incrad(b, i);
                wt *= wt*PI * 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 */
                        wt = BSDF_visible(b, i, getBSDF_outndx(b, dv));
                        copycolor(col, cinc);
                        scalecolor(col, wt);
                        csum = &distarr[3*(k*nazi + j)];
                        addcolor(csum, col);    /* sum into distribution */
                    }
        }
        free(idist);                    /* free temp space */
}
Revision:	2.5
Committed:	Wed Dec 5 20:07:34 2007 UTC (17 years, 10 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.4:	+38 -85 lines
Log Message:	Improved BSDF resampling
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.5	static const char RCSid[] = "$Id: mkillum4.c,v 2.4 2007/12/05 05:02:58 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.3	#include "random.h"
11	greg	2.2	#include "ezxml.h"
12	greg	2.1
13
14			struct BSDF_data *
15			load_BSDF( /* load BSDF data from file */
16			char *fname
17			)
18			{
19			char *path;
20	greg	2.2	ezxml_t fl, wld, wdb;
21	greg	2.1	struct BSDF_data *dp;
22
23			path = getpath(fname, getrlibpath(), R_OK);
24			if (path == NULL) {
25			sprintf(errmsg, "cannot find BSDF file \"%s\"", fname);
26			error(WARNING, errmsg);
27			return(NULL);
28			}
29	greg	2.2	fl = ezxml_parse_file(path);
30			if (fl == NULL) {
31	greg	2.1	sprintf(errmsg, "cannot open BSDF \"%s\"", path);
32			error(WARNING, errmsg);
33			return(NULL);
34			}
35	greg	2.5	dp = (struct BSDF_data *)calloc(1, sizeof(struct BSDF_data));
36	greg	2.1	if (dp == NULL)
37			goto memerr;
38	greg	2.2	for (wld = ezxml_child(fl, "WavelengthData");
39			fl != NULL; fl = fl->next) {
40			if (strcmp(ezxml_child(wld, "Wavelength")->txt, "Visible"))
41			continue;
42			wdb = ezxml_child(wld, "WavelengthDataBlock");
43			if (wdb == NULL) continue;
44			if (strcmp(ezxml_child(wdb, "WavelengthDataDirection")->txt,
45			"Transmission Front"))
46			continue;
47			}
48	greg	2.1	/* etc... */
49	greg	2.2	ezxml_free(fl);
50	greg	2.1	return(dp);
51			memerr:
52			error(SYSTEM, "out of memory in load_BSDF");
53			return NULL; /* pro forma return */
54			}
55
56
57			void
58			free_BSDF( /* free BSDF data structure */
59			struct BSDF_data *b
60			)
61			{
62			if (b == NULL)
63			return;
64			free(b->bsdf);
65			free(b);
66			}
67
68
69			void
70			r_BSDF_incvec( /* compute random input vector at given location */
71			FVECT v,
72			struct BSDF_data *b,
73			int i,
74			double rv,
75			MAT4 xm
76			)
77			{
78			FVECT pert;
79			double rad;
80			int j;
81
82			getBSDF_incvec(v, b, i);
83			rad = getBSDF_incrad(b, i);
84			multisamp(pert, 3, rv);
85			for (j = 0; j < 3; j++)
86			v[j] += rad(2.pert[j] - 1.);
87			if (xm != NULL)
88			multv3(v, v, xm);
89			normalize(v);
90			}
91
92
93			void
94			r_BSDF_outvec( /* compute random output vector at given location */
95			FVECT v,
96			struct BSDF_data *b,
97			int o,
98			double rv,
99			MAT4 xm
100			)
101			{
102			FVECT pert;
103			double rad;
104			int j;
105
106			getBSDF_outvec(v, b, o);
107			rad = getBSDF_outrad(b, o);
108			multisamp(pert, 3, rv);
109			for (j = 0; j < 3; j++)
110			v[j] += rad(2.pert[j] - 1.);
111			if (xm != NULL)
112			multv3(v, v, xm);
113			normalize(v);
114			}
115	greg	2.2
116
117			#define FEQ(a,b) ((a)-(b) <= 1e-7 && (b)-(a) <= 1e-7)
118
119			static int
120			addrot( /* compute rotation (x,y,z) => (xp,yp,zp) */
121			char *xfarg[],
122			FVECT xp,
123			FVECT yp,
124			FVECT zp
125			)
126			{
127			static char bufs[3][16];
128			int bn = 0;
129			char **xfp = xfarg;
130			double theta;
131
132			theta = atan2(yp[2], zp[2]);
133			if (!FEQ(theta,0.0)) {
134			*xfp++ = "-rx";
135			sprintf(bufs[bn], "%f", theta*(180./PI));
136			*xfp++ = bufs[bn++];
137			}
138			theta = asin(-xp[2]);
139			if (!FEQ(theta,0.0)) {
140			*xfp++ = "-ry";
141			sprintf(bufs[bn], " %f", theta*(180./PI));
142			*xfp++ = bufs[bn++];
143			}
144			theta = atan2(xp[1], xp[0]);
145			if (!FEQ(theta,0.0)) {
146			*xfp++ = "-rz";
147			sprintf(bufs[bn], "%f", theta*(180./PI));
148			*xfp++ = bufs[bn++];
149			}
150			*xfp = NULL;
151			return(xfp - xfarg);
152			}
153
154
155			int
156			getBSDF_xfm( /* compute transform for the given surface */
157			MAT4 xm,
158			FVECT nrm,
159			UpDir ud
160			)
161			{
162			char *xfargs[7];
163			XF myxf;
164			FVECT updir, xdest, ydest;
165
166			updir[0] = updir[1] = updir[2] = 0.;
167			switch (ud) {
168			case UDzneg:
169			updir[2] = -1.;
170			break;
171			case UDyneg:
172			updir[1] = -1.;
173			break;
174			case UDxneg:
175			updir[0] = -1.;
176			break;
177			case UDxpos:
178			updir[0] = 1.;
179			break;
180			case UDypos:
181			updir[1] = 1.;
182			break;
183			case UDzpos:
184			updir[2] = 1.;
185			break;
186			case UDunknown:
187			error(WARNING, "unspecified up direction");
188			return(0);
189			}
190			fcross(xdest, updir, nrm);
191			if (normalize(xdest) == 0.0)
192			return(0);
193			fcross(ydest, nrm, xdest);
194			xf(&myxf, addrot(xfargs, xdest, ydest, nrm), xfargs);
195			copymat4(xm, myxf.xfm);
196			return(1);
197			}
198
199
200			void
201			redistribute( /* pass distarr ray sums through BSDF */
202			struct BSDF_data *b,
203			int nalt,
204			int nazi,
205			FVECT u,
206			FVECT v,
207			FVECT w,
208			MAT4 xm
209			)
210			{
211	greg	2.5	MAT4 mymat, inmat;
212			COLORV *idist;
213	greg	2.2	COLORV cp, csum;
214			FVECT dv;
215	greg	2.5	double wt;
216			int i, j, k, h;
217			COLOR col, cinc;
218			/* copy incoming distribution */
219			if (b->ninc != distsiz)
220			error(INTERNAL, "error 1 in redistribute");
221			idist = (COLORV )malloc(sizeof(COLOR)distsiz);
222			if (idist == NULL)
223	greg	2.2	error(SYSTEM, "out of memory in redistribute");
224	greg	2.5	memcpy(idist, distarr, sizeof(COLOR)*distsiz);
225			/* compose direction transform */
226	greg	2.2	for (i = 3; i--; ) {
227			mymat[i][0] = u[i];
228			mymat[i][1] = v[i];
229			mymat[i][2] = w[i];
230			mymat[i][3] = 0.;
231			}
232			mymat[3][0] = mymat[3][1] = mymat[3][2] = 0.;
233			mymat[3][3] = 1.;
234			if (xm != NULL)
235			multmat4(mymat, xm, mymat);
236			for (i = 3; i--; ) { /* make sure it's normalized */
237	greg	2.5	wt = 1./sqrt( mymat[0][i]*mymat[0][i] +
238	greg	2.2	mymat[1][i]*mymat[1][i] +
239			mymat[2][i]*mymat[2][i] );
240			for (j = 3; j--; )
241	greg	2.5	mymat[j][i] *= wt;
242	greg	2.2	}
243	greg	2.5	if (!invmat4(inmat, mymat)) /* need inverse as well */
244			error(INTERNAL, "cannot invert BSDF transform");
245			newdist(naltnazi); / resample distribution */
246	greg	2.4	for (i = b->ninc; i--; ) {
247	greg	2.5	getBSDF_incvec(dv, b, i); /* compute incident irrad. */
248	greg	2.2	multv3(dv, dv, mymat);
249	greg	2.5	if (dv[2] < 0.0) dv[2] = -dv[2];
250			wt = getBSDF_incrad(b, i);
251			wt = wtPI * dv[2]; /* solid_anglecosine(theta) /
252			cp = &idist[3*i];
253			copycolor(cinc, cp);
254			scalecolor(cinc, wt);
255			for (k = nalt; k--; ) /* loop over distribution */
256			for (j = nazi; j--; ) {
257			flatdir(dv, (k + .5)/nalt, (double)j/nazi);
258			multv3(dv, dv, inmat);
259			/* evaluate BSDF @ outgoing */
260			wt = BSDF_visible(b, i, getBSDF_outndx(b, dv));
261			copycolor(col, cinc);
262			scalecolor(col, wt);
263			csum = &distarr[3(knazi + j)];
264			addcolor(csum, col); /* sum into distribution */
265			}
266	greg	2.2	}
267	greg	2.5	free(idist); /* free temp space */
268	greg	2.2	}