src/gen/mkillum4.c

#ifndef lint
static const char RCSid[] = "$Id: mkillum4.c,v 2.3 2007/11/21 18:51:05 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 *)malloc(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->inc_dir);
        free(b->inc_rad);
        free(b->out_dir);
        free(b->out_rad);
        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;
        COLORV  *outarr;
        COLORV  *cp, *csum;
        uint16  *distcnt;
        FVECT   dv;
        double  wt0, wt1;
        int     i, j, o;
        int     cnt;
        COLOR   col;
                                        /* allocate temporary memory */
        outarr = (COLORV *)calloc(b->nout, sizeof(COLOR));
        distcnt = (uint16 *)calloc(nalt*nazi, sizeof(uint16));
        if ((outarr == NULL) | (distcnt == NULL))
                error(SYSTEM, "out of memory in redistribute");
                                        /* compose matrix */
        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 */
                wt0 = 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] *= wt0;
        }
                                        /* pass through BSDF */
        for (i = b->ninc; i--; ) {
                getBSDF_incvec(dv, b, i);
                multv3(dv, dv, mymat);
                wt0 = getBSDF_incrad(b, i);
                wt0 *= PI*wt0 * dv[2];          /* solid_angle*cosine(theta) */
                for (o = b->nout; o--; ) {
                        cp = &distarr[3*i];
                        csum = &outarr[3*o];
                        wt1 = wt0 * BSDF_data(b,i,o);
                        copycolor(col, cp);
                        scalecolor(col, wt1);
                        addcolor(csum, col);
                }
        }
        newdist(nalt*nazi);             /* resample distribution */
        for (o = b->nout; o--; ) {
                getBSDF_outvec(dv, b, o);
                multv3(dv, dv, mymat);
                j = (.5 + atan2(dv[1],dv[0])*(.5/PI))*nazi + .5;
                if (j >= nazi) j = 0;
                i = (0.9999 - dv[2]*dv[2])*nalt;
                csum = &distarr[3*(i*nazi + j)];
                cp = &outarr[3*o];
                addcolor(csum, cp);
                ++distcnt[i*nazi + j];
        }
        free(outarr);
                                        /* fill in missing bits */
        for (i = nalt; i--; )
            for (j = nazi; j--; ) {
                int     ii, jj, alt, azi;
                if (distcnt[i*nazi + j])
                        continue;
                csum = &distarr[3*(i*nazi + j)];
                setcolor(csum, 0., 0., 0.);
                cnt = 0;
                for (o = 0; !cnt; o++)
                    for (ii = -o; ii <= o; ii++) {
                        alt = i + ii;
                        if (alt < 0) continue;
                        if (alt >= nalt) break;
                        for (jj = -o; jj <= o; jj++) {
                            if (ii*ii + jj*jj != o*o)
                                continue;
                            azi = j + jj;
                            if (azi >= nazi) azi -= nazi;
                            else if (azi < 0) azi += nazi;
                            if (!distcnt[alt*nazi + azi])
                                continue;
                            cp = &distarr[3*(alt*nazi + azi)];
                            addcolor(csum, cp);
                            cnt += distcnt[alt*nazi + azi];
                        }
                    }
                wt0 = 1./cnt;
                scalecolor(csum, wt0);
            }
                                        /* finish averages */
        for (i = nalt; i--; )
            for (j = nazi; j--; ) {
                if ((cnt = distcnt[i*nazi + j]) <= 1)
                        continue;
                csum = &distarr[3*(i*nazi + j)];
                wt0 = 1./cnt;
                scalecolor(csum, wt0);
            }
        free(distcnt);
}
Revision:	2.4
Committed:	Wed Dec 5 05:02:58 2007 UTC (16 years, 6 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.3:	+17 -27 lines
Log Message:	Fixed resampling coefficients
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.4	static const char RCSid[] = "$Id: mkillum4.c,v 2.3 2007/11/21 18:51:05 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			dp = (struct BSDF_data *)malloc(sizeof(struct BSDF_data));
36			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->inc_dir);
65			free(b->inc_rad);
66			free(b->out_dir);
67			free(b->out_rad);
68			free(b->bsdf);
69			free(b);
70			}
71
72
73			void
74			r_BSDF_incvec( /* compute random input vector at given location */
75			FVECT v,
76			struct BSDF_data *b,
77			int i,
78			double rv,
79			MAT4 xm
80			)
81			{
82			FVECT pert;
83			double rad;
84			int j;
85
86			getBSDF_incvec(v, b, i);
87			rad = getBSDF_incrad(b, i);
88			multisamp(pert, 3, rv);
89			for (j = 0; j < 3; j++)
90			v[j] += rad(2.pert[j] - 1.);
91			if (xm != NULL)
92			multv3(v, v, xm);
93			normalize(v);
94			}
95
96
97			void
98			r_BSDF_outvec( /* compute random output vector at given location */
99			FVECT v,
100			struct BSDF_data *b,
101			int o,
102			double rv,
103			MAT4 xm
104			)
105			{
106			FVECT pert;
107			double rad;
108			int j;
109
110			getBSDF_outvec(v, b, o);
111			rad = getBSDF_outrad(b, o);
112			multisamp(pert, 3, rv);
113			for (j = 0; j < 3; j++)
114			v[j] += rad(2.pert[j] - 1.);
115			if (xm != NULL)
116			multv3(v, v, xm);
117			normalize(v);
118			}
119	greg	2.2
120
121			#define FEQ(a,b) ((a)-(b) <= 1e-7 && (b)-(a) <= 1e-7)
122
123			static int
124			addrot( /* compute rotation (x,y,z) => (xp,yp,zp) */
125			char *xfarg[],
126			FVECT xp,
127			FVECT yp,
128			FVECT zp
129			)
130			{
131			static char bufs[3][16];
132			int bn = 0;
133			char **xfp = xfarg;
134			double theta;
135
136			theta = atan2(yp[2], zp[2]);
137			if (!FEQ(theta,0.0)) {
138			*xfp++ = "-rx";
139			sprintf(bufs[bn], "%f", theta*(180./PI));
140			*xfp++ = bufs[bn++];
141			}
142			theta = asin(-xp[2]);
143			if (!FEQ(theta,0.0)) {
144			*xfp++ = "-ry";
145			sprintf(bufs[bn], " %f", theta*(180./PI));
146			*xfp++ = bufs[bn++];
147			}
148			theta = atan2(xp[1], xp[0]);
149			if (!FEQ(theta,0.0)) {
150			*xfp++ = "-rz";
151			sprintf(bufs[bn], "%f", theta*(180./PI));
152			*xfp++ = bufs[bn++];
153			}
154			*xfp = NULL;
155			return(xfp - xfarg);
156			}
157
158
159			int
160			getBSDF_xfm( /* compute transform for the given surface */
161			MAT4 xm,
162			FVECT nrm,
163			UpDir ud
164			)
165			{
166			char *xfargs[7];
167			XF myxf;
168			FVECT updir, xdest, ydest;
169
170			updir[0] = updir[1] = updir[2] = 0.;
171			switch (ud) {
172			case UDzneg:
173			updir[2] = -1.;
174			break;
175			case UDyneg:
176			updir[1] = -1.;
177			break;
178			case UDxneg:
179			updir[0] = -1.;
180			break;
181			case UDxpos:
182			updir[0] = 1.;
183			break;
184			case UDypos:
185			updir[1] = 1.;
186			break;
187			case UDzpos:
188			updir[2] = 1.;
189			break;
190			case UDunknown:
191			error(WARNING, "unspecified up direction");
192			return(0);
193			}
194			fcross(xdest, updir, nrm);
195			if (normalize(xdest) == 0.0)
196			return(0);
197			fcross(ydest, nrm, xdest);
198			xf(&myxf, addrot(xfargs, xdest, ydest, nrm), xfargs);
199			copymat4(xm, myxf.xfm);
200			return(1);
201			}
202
203
204			void
205			redistribute( /* pass distarr ray sums through BSDF */
206			struct BSDF_data *b,
207			int nalt,
208			int nazi,
209			FVECT u,
210			FVECT v,
211			FVECT w,
212			MAT4 xm
213			)
214			{
215			MAT4 mymat;
216			COLORV *outarr;
217			COLORV cp, csum;
218			uint16 *distcnt;
219			FVECT dv;
220	greg	2.4	double wt0, wt1;
221	greg	2.2	int i, j, o;
222			int cnt;
223			COLOR col;
224			/* allocate temporary memory */
225	greg	2.4	outarr = (COLORV *)calloc(b->nout, sizeof(COLOR));
226	greg	2.2	distcnt = (uint16 )calloc(naltnazi, sizeof(uint16));
227	greg	2.4	if ((outarr == NULL) \| (distcnt == NULL))
228	greg	2.2	error(SYSTEM, "out of memory in redistribute");
229			/* compose matrix */
230			for (i = 3; i--; ) {
231			mymat[i][0] = u[i];
232			mymat[i][1] = v[i];
233			mymat[i][2] = w[i];
234			mymat[i][3] = 0.;
235			}
236			mymat[3][0] = mymat[3][1] = mymat[3][2] = 0.;
237			mymat[3][3] = 1.;
238			if (xm != NULL)
239			multmat4(mymat, xm, mymat);
240			for (i = 3; i--; ) { /* make sure it's normalized */
241	greg	2.4	wt0 = 1./sqrt( mymat[0][i]*mymat[0][i] +
242	greg	2.2	mymat[1][i]*mymat[1][i] +
243			mymat[2][i]*mymat[2][i] );
244			for (j = 3; j--; )
245	greg	2.4	mymat[j][i] *= wt0;
246	greg	2.2	}
247			/* pass through BSDF */
248	greg	2.4	for (i = b->ninc; i--; ) {
249	greg	2.2	getBSDF_incvec(dv, b, i);
250			multv3(dv, dv, mymat);
251	greg	2.4	wt0 = getBSDF_incrad(b, i);
252			wt0 = PIwt0 * dv[2]; /* solid_anglecosine(theta) /
253			for (o = b->nout; o--; ) {
254	greg	2.2	cp = &distarr[3*i];
255	greg	2.4	csum = &outarr[3*o];
256			wt1 = wt0 * BSDF_data(b,i,o);
257	greg	2.2	copycolor(col, cp);
258	greg	2.4	scalecolor(col, wt1);
259	greg	2.2	addcolor(csum, col);
260			}
261			}
262			newdist(naltnazi); / resample distribution */
263			for (o = b->nout; o--; ) {
264			getBSDF_outvec(dv, b, o);
265			multv3(dv, dv, mymat);
266			j = (.5 + atan2(dv[1],dv[0])(.5/PI))nazi + .5;
267			if (j >= nazi) j = 0;
268			i = (0.9999 - dv[2]dv[2])nalt;
269			csum = &distarr[3(inazi + j)];
270			cp = &outarr[3*o];
271			addcolor(csum, cp);
272			++distcnt[i*nazi + j];
273			}
274			free(outarr);
275			/* fill in missing bits */
276			for (i = nalt; i--; )
277			for (j = nazi; j--; ) {
278			int ii, jj, alt, azi;
279			if (distcnt[i*nazi + j])
280			continue;
281			csum = &distarr[3(inazi + j)];
282			setcolor(csum, 0., 0., 0.);
283			cnt = 0;
284			for (o = 0; !cnt; o++)
285			for (ii = -o; ii <= o; ii++) {
286			alt = i + ii;
287			if (alt < 0) continue;
288			if (alt >= nalt) break;
289			for (jj = -o; jj <= o; jj++) {
290			if (iiii + jjjj != o*o)
291			continue;
292			azi = j + jj;
293			if (azi >= nazi) azi -= nazi;
294			else if (azi < 0) azi += nazi;
295			if (!distcnt[alt*nazi + azi])
296			continue;
297			cp = &distarr[3(altnazi + azi)];
298			addcolor(csum, cp);
299			cnt += distcnt[alt*nazi + azi];
300			}
301			}
302	greg	2.4	wt0 = 1./cnt;
303			scalecolor(csum, wt0);
304	greg	2.2	}
305			/* finish averages */
306			for (i = nalt; i--; )
307			for (j = nazi; j--; ) {
308			if ((cnt = distcnt[i*nazi + j]) <= 1)
309			continue;
310			csum = &distarr[3(inazi + j)];
311	greg	2.4	wt0 = 1./cnt;
312			scalecolor(csum, wt0);
313	greg	2.2	}
314			free(distcnt);
315			}