src/gen/mkillum4.c

#ifndef lint
static const char RCSid[] = "$Id: mkillum4.c,v 2.1 2007/09/18 19:51:07 greg Exp $";
#endif
/*
 * Routines for handling BSDF data within mkillum
 */

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