1 |
#ifndef lint |
2 |
static const char RCSid[] = "$Id: mkillum4.c,v 2.1 2007/09/18 19:51:07 greg Exp $"; |
3 |
#endif |
4 |
/* |
5 |
* Routines for handling BSDF data within mkillum |
6 |
*/ |
7 |
|
8 |
#include "mkillum.h" |
9 |
#include "paths.h" |
10 |
#include "ezxml.h" |
11 |
|
12 |
|
13 |
struct BSDF_data * |
14 |
load_BSDF( /* load BSDF data from file */ |
15 |
char *fname |
16 |
) |
17 |
{ |
18 |
char *path; |
19 |
ezxml_t fl, wld, wdb; |
20 |
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 |
fl = ezxml_parse_file(path); |
29 |
if (fl == NULL) { |
30 |
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 |
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 |
/* etc... */ |
48 |
ezxml_free(fl); |
49 |
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 |
|
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(nalt*nazi, 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] = PI*wt*wt * dv[2]; /* solid_angle*cosine(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(nalt*nazi); /* 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*(i*nazi + 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*(i*nazi + 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 (ii*ii + jj*jj != 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*(alt*nazi + 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*(i*nazi + j)]; |
320 |
wt = 1./cnt; |
321 |
scalecolor(csum, wt); |
322 |
} |
323 |
free(distcnt); |
324 |
} |