src/gen/mkillum4.c

#ifndef lint
static const char RCSid[] = "$Id: mkillum4.c,v 2.2 2007/09/21 05:53:21 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;
        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.3
Committed:	Wed Nov 21 18:51:05 2007 UTC (17 years, 5 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.2:	+2 -1 lines
Log Message:	Fixed implicit declarations
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.3	static const char RCSid[] = "$Id: mkillum4.c,v 2.2 2007/09/21 05:53:21 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			float *inpcoef;
218			COLORV cp, csum;
219			uint16 *distcnt;
220			FVECT dv;
221			double oom, wt;
222			int i, j, o;
223			int cnt;
224			COLOR col;
225			/* allocate temporary memory */
226			outarr = (COLORV )malloc(b->nout sizeof(COLOR));
227			distcnt = (uint16 )calloc(naltnazi, sizeof(uint16));
228			inpcoef = (float )malloc(b->ninc sizeof(float));
229			if ((outarr == NULL) \| (distcnt == NULL) \| (inpcoef == NULL))
230			error(SYSTEM, "out of memory in redistribute");
231			/* compose matrix */
232			for (i = 3; i--; ) {
233			mymat[i][0] = u[i];
234			mymat[i][1] = v[i];
235			mymat[i][2] = w[i];
236			mymat[i][3] = 0.;
237			}
238			mymat[3][0] = mymat[3][1] = mymat[3][2] = 0.;
239			mymat[3][3] = 1.;
240			if (xm != NULL)
241			multmat4(mymat, xm, mymat);
242			for (i = 3; i--; ) { /* make sure it's normalized */
243			wt = 1./sqrt( mymat[0][i]*mymat[0][i] +
244			mymat[1][i]*mymat[1][i] +
245			mymat[2][i]*mymat[2][i] );
246			for (j = 3; j--; )
247			mymat[j][i] *= wt;
248			}
249			/* pass through BSDF */
250			for (i = b->ninc; i--; ) { /* get input coefficients */
251			getBSDF_incvec(dv, b, i);
252			multv3(dv, dv, mymat);
253			wt = getBSDF_incrad(b, i);
254			inpcoef[i] = PIwtwt * dv[2]; /* solid_anglecosine(theta) /
255			}
256			for (o = b->nout; o--; ) {
257			csum = &outarr[3*o];
258			setcolor(csum, 0., 0., 0.);
259			oom = getBSDF_outrad(b, o);
260			oom = oom PI;
261			for (i = b->ninc; i--; ) {
262			wt = BSDF_data(b,i,o) * inpcoef[i] / oom;
263			cp = &distarr[3*i];
264			copycolor(col, cp);
265			scalecolor(col, wt);
266			addcolor(csum, col);
267			}
268			wt = 1./b->ninc;
269			scalecolor(csum, wt);
270			}
271			free(inpcoef);
272			newdist(naltnazi); / resample distribution */
273			for (o = b->nout; o--; ) {
274			getBSDF_outvec(dv, b, o);
275			multv3(dv, dv, mymat);
276			j = (.5 + atan2(dv[1],dv[0])(.5/PI))nazi + .5;
277			if (j >= nazi) j = 0;
278			i = (0.9999 - dv[2]dv[2])nalt;
279			csum = &distarr[3(inazi + j)];
280			cp = &outarr[3*o];
281			addcolor(csum, cp);
282			++distcnt[i*nazi + j];
283			}
284			free(outarr);
285			/* fill in missing bits */
286			for (i = nalt; i--; )
287			for (j = nazi; j--; ) {
288			int ii, jj, alt, azi;
289			if (distcnt[i*nazi + j])
290			continue;
291			csum = &distarr[3(inazi + j)];
292			setcolor(csum, 0., 0., 0.);
293			cnt = 0;
294			for (o = 0; !cnt; o++)
295			for (ii = -o; ii <= o; ii++) {
296			alt = i + ii;
297			if (alt < 0) continue;
298			if (alt >= nalt) break;
299			for (jj = -o; jj <= o; jj++) {
300			if (iiii + jjjj != o*o)
301			continue;
302			azi = j + jj;
303			if (azi >= nazi) azi -= nazi;
304			else if (azi < 0) azi += nazi;
305			if (!distcnt[alt*nazi + azi])
306			continue;
307			cp = &distarr[3(altnazi + azi)];
308			addcolor(csum, cp);
309			cnt += distcnt[alt*nazi + azi];
310			}
311			}
312			wt = 1./cnt;
313			scalecolor(csum, wt);
314			}
315			/* finish averages */
316			for (i = nalt; i--; )
317			for (j = nazi; j--; ) {
318			if ((cnt = distcnt[i*nazi + j]) <= 1)
319			continue;
320			csum = &distarr[3(inazi + j)];
321			wt = 1./cnt;
322			scalecolor(csum, wt);
323			}
324			free(distcnt);
325			}