ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/mkillum4.c
(Generate patch)

Comparing ray/src/gen/mkillum4.c (file contents):
Revision 2.1 by greg, Tue Sep 18 19:51:07 2007 UTC vs.
Revision 2.6 by greg, Sat Dec 8 01:43:09 2007 UTC

# Line 7 | Line 7 | static const char RCSid[] = "$Id$";
7  
8   #include "mkillum.h"
9   #include "paths.h"
10 + #include "ezxml.h"
11  
12 + #define MAXLATS         46              /* maximum number of latitudes */
13  
14 + /* BSDF angle specification (terminate with nphi = -1) */
15 + typedef struct {
16 +        char    name[32];               /* basis name */
17 +        int     nangles;                /* total number of directions */
18 +        struct {
19 +                float   tmin;                   /* starting theta */
20 +                short   nphis;                  /* number of phis (0 term) */
21 +        }       lat[MAXLATS+1];         /* latitudes */
22 + } ANGLE_BASIS;
23 +
24 + #define NABASES         3               /* number of predefined bases */
25 +
26 + ANGLE_BASIS     abase_list[NABASES] = {
27 +        {
28 +                "LBNL/Klems Full", 145,
29 +                { {-5., 1},
30 +                {5., 8},
31 +                {15., 16},
32 +                {25., 20},
33 +                {35., 24},
34 +                {45., 24},
35 +                {55., 24},
36 +                {65., 16},
37 +                {75., 12},
38 +                {90., 0} }
39 +        }, {
40 +                "LBNL/Klems Half", 73,
41 +                { {-6.5, 1},
42 +                {6.5, 8},
43 +                {19.5, 12},
44 +                {32.5, 16},
45 +                {46.5, 20},
46 +                {61.5, 12},
47 +                {76.5, 4},
48 +                {90., 0} }
49 +        }, {
50 +                "LBNL/Klems Quarter", 41,
51 +                { {-9., 1},
52 +                {9., 8},
53 +                {27., 12},
54 +                {46., 12},
55 +                {66., 8},
56 +                {90., 0} }
57 +        }
58 + };
59 +
60 +
61 + static int
62 + ab_getvec(              /* get vector for this angle basis index */
63 +        FVECT v,
64 +        int ndx,
65 +        void *p
66 + )
67 + {
68 +        ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
69 +        int     li;
70 +        double  alt, azi, d;
71 +        
72 +        if ((ndx < 0) | (ndx >= ab->nangles))
73 +                return(0);
74 +        for (li = 0; ndx >= ab->lat[li].nphis; li++)
75 +                ndx -= ab->lat[li].nphis;
76 +        alt = PI/180.*0.5*(ab->lat[li].tmin + ab->lat[li+1].tmin);
77 +        azi = 2.*PI*ndx/ab->lat[li].nphis;
78 +        d = sin(alt);
79 +        v[0] = cos(azi)*d;
80 +        v[1] = sin(azi)*d;
81 +        v[2] = cos(alt);
82 +        return(1);
83 + }
84 +
85 +
86 + static int
87 + ab_getndx(              /* get index corresponding to the given vector */
88 +        FVECT v,
89 +        void *p
90 + )
91 + {
92 +        ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
93 +        int     li, ndx;
94 +        double  alt, azi, d;
95 +
96 +        if ((v[2] < 0.0) | (v[2] > 1.0))
97 +                return(-1);
98 +        alt = 180.0/PI*acos(v[2]);
99 +        azi = 180.0/PI*atan2(v[1], v[0]);
100 +        if (azi < 0.0) azi += 360.0;
101 +        for (li = 1; ab->lat[li].tmin <= alt; li++)
102 +                if (!ab->lat[li].nphis)
103 +                        return(-1);
104 +        --li;
105 +        ndx = (int)((1./360.)*azi*ab->lat[li].nphis + 0.5);
106 +        if (ndx >= ab->lat[li].nphis) ndx = 0;
107 +        while (li--)
108 +                ndx += ab->lat[li].nphis;
109 +        return(ndx);
110 + }
111 +
112 +
113 + static double
114 + ab_getohm(              /* get solid angle for this angle basis index */
115 +        int ndx,
116 +        void *p
117 + )
118 + {
119 +        ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
120 +        int     li;
121 +        double  tdia, pdia;
122 +        
123 +        if ((ndx < 0) | (ndx >= ab->nangles))
124 +                return(0);
125 +        for (li = 0; ndx >= ab->lat[li].nphis; li++)
126 +                ndx -= ab->lat[li].nphis;
127 +        tdia = PI/180.*(ab->lat[li+1].tmin - ab->lat[li].tmin);
128 +        pdia = 2.*PI/ab->lat[li].nphis;
129 +        return(tdia*pdia);
130 + }
131 +
132 +
133 + static int
134 + ab_getvecR(             /* get reverse vector for this angle basis index */
135 +        FVECT v,
136 +        int ndx,
137 +        void *p
138 + )
139 + {
140 +        if (!ab_getvec(v, ndx, p))
141 +                return(0);
142 +
143 +        v[0] = -v[0];
144 +        v[1] = -v[1];
145 +
146 +        return(1);
147 + }
148 +
149 +
150 + static int
151 + ab_getndxR(             /* get index corresponding to the reverse vector */
152 +        FVECT v,
153 +        void *p
154 + )
155 + {
156 +        FVECT  v2;
157 +        
158 +        v2[0] = -v[0];
159 +        v2[1] = -v[1];
160 +        v2[2] = v[2];
161 +
162 +        return ab_getndx(v2, p);
163 + }
164 +
165 + static void
166 + load_bsdf_data(         /* load BSDF distribution for this wavelength */
167 +        struct BSDF_data *dp,
168 +        ezxml_t wdb
169 + )
170 + {
171 +        const char  *cbasis = ezxml_txt(ezxml_child(wdb,"ColumnAngleBasis"));
172 +        const char  *rbasis = ezxml_txt(ezxml_child(wdb,"RowAngleBasis"));
173 +        int  i;
174 +        
175 +        for (i = NABASES; i--; )
176 +                if (!strcmp(cbasis, abase_list[i].name)) {
177 +                        dp->ninc = abase_list[i].nangles;
178 +                        dp->ib_priv = (void *)&abase_list[i];
179 +                        dp->ib_vec = ab_getvec;
180 +                        dp->ib_ndx = ab_getndx;
181 +                        dp->ib_ohm = ab_getohm;
182 +                        break;
183 +                }
184 +        if (i < 0) {
185 +                sprintf(errmsg, "unsupported incident basis '%s'", cbasis);
186 +                error(INTERNAL, errmsg);
187 +        }
188 +        for (i = NABASES; i--; )
189 +                if (!strcmp(rbasis, abase_list[i].name)) {
190 +                        dp->nout = abase_list[i].nangles;
191 +                        dp->ob_priv = (void *)&abase_list[i];
192 +                        dp->ob_vec = ab_getvecR;
193 +                        dp->ob_ndx = ab_getndxR;
194 +                        dp->ob_ohm = ab_getohm;
195 +                        break;
196 +                }
197 +        if (i < 0) {
198 +                sprintf(errmsg, "unsupported exitant basis '%s'", cbasis);
199 +                error(INTERNAL, errmsg);
200 +        }
201 + }
202 +
203 +
204   struct BSDF_data *
205   load_BSDF(              /* load BSDF data from file */
206          char *fname
207   )
208   {
209          char                    *path;
210 <        FILE                    *fp;
210 >        ezxml_t                 fl, wld, wdb;
211          struct BSDF_data        *dp;
212          
213          path = getpath(fname, getrlibpath(), R_OK);
# Line 24 | Line 216 | load_BSDF(             /* load BSDF data from file */
216                  error(WARNING, errmsg);
217                  return(NULL);
218          }
219 <        if ((fp = fopen(path, "r")) == NULL) {
219 >        fl = ezxml_parse_file(path);
220 >        if (fl == NULL) {
221                  sprintf(errmsg, "cannot open BSDF \"%s\"", path);
222                  error(WARNING, errmsg);
223                  return(NULL);
224          }
225 <        dp = (struct BSDF_data *)malloc(sizeof(struct BSDF_data));
226 <        if (dp == NULL)
227 <                goto memerr;
228 <        /* etc... */
229 <        fclose(fp);
225 >        if (ezxml_error(fl)[0]) {
226 >                sprintf(errmsg, "BSDF \"%s\": %s", path, ezxml_error(fl));
227 >                error(WARNING, errmsg);
228 >                ezxml_free(fl);
229 >                return(NULL);
230 >        }
231 >        dp = (struct BSDF_data *)calloc(1, sizeof(struct BSDF_data));
232 >        for (wld = ezxml_child(fl, "WavelengthData");
233 >                                fl != NULL; fl = fl->next) {
234 >                if (strcmp(ezxml_txt(ezxml_child(wld,"Wavelength")), "Visible"))
235 >                        continue;
236 >                wdb = ezxml_child(wld, "WavelengthDataBlock");
237 >                if (wdb == NULL) continue;
238 >                if (strcmp(ezxml_txt(ezxml_child(wdb,"WavelengthDataDirection")),
239 >                                        "Transmission Front"))
240 >                        continue;
241 >                load_bsdf_data(dp, wdb);        /* load front BTDF */
242 >                break;                          /* ignore the rest */
243 >        }
244 >        ezxml_free(fl);                         /* done with XML file */
245 >        if (dp->bsdf == NULL) {
246 >                sprintf(errmsg, "bad/missing BTDF data in \"%s\"", path);
247 >                error(WARNING, errmsg);
248 >                free_BSDF(dp);
249 >                dp = NULL;
250 >        }
251          return(dp);
38 – memerr:
39 –        error(SYSTEM, "out of memory in load_BSDF");
40 –        return NULL;    /* pro forma return */
252   }
253  
254  
# Line 48 | Line 259 | free_BSDF(             /* free BSDF data structure */
259   {
260          if (b == NULL)
261                  return;
262 <        free(b->inc_dir);
263 <        free(b->inc_rad);
53 <        free(b->out_dir);
54 <        free(b->out_rad);
55 <        free(b->bsdf);
262 >        if (b->bsdf != NULL)
263 >                free(b->bsdf);
264          free(b);
265   }
266  
267  
268 < void
268 > int
269   r_BSDF_incvec(          /* compute random input vector at given location */
270          FVECT v,
271          struct BSDF_data *b,
# Line 70 | Line 278 | r_BSDF_incvec(         /* compute random input vector at give
278          double  rad;
279          int     j;
280          
281 <        getBSDF_incvec(v, b, i);
282 <        rad = getBSDF_incrad(b, i);
281 >        if (!getBSDF_incvec(v, b, i))
282 >                return(0);
283 >        rad = sqrt(getBSDF_incohm(b, i) / PI);
284          multisamp(pert, 3, rv);
285          for (j = 0; j < 3; j++)
286                  v[j] += rad*(2.*pert[j] - 1.);
287          if (xm != NULL)
288                  multv3(v, v, xm);
289          normalize(v);
290 +        return(1);
291   }
292  
293  
294 < void
294 > int
295   r_BSDF_outvec(          /* compute random output vector at given location */
296          FVECT v,
297          struct BSDF_data *b,
# Line 94 | Line 304 | r_BSDF_outvec(         /* compute random output vector at giv
304          double  rad;
305          int     j;
306          
307 <        getBSDF_outvec(v, b, o);
308 <        rad = getBSDF_outrad(b, o);
307 >        if (!getBSDF_outvec(v, b, o))
308 >                return(0);
309 >        rad = sqrt(getBSDF_outohm(b, o) / PI);
310          multisamp(pert, 3, rv);
311          for (j = 0; j < 3; j++)
312                  v[j] += rad*(2.*pert[j] - 1.);
313          if (xm != NULL)
314                  multv3(v, v, xm);
315          normalize(v);
316 +        return(1);
317 + }
318 +
319 +
320 + #define  FEQ(a,b)       ((a)-(b) <= 1e-7 && (b)-(a) <= 1e-7)
321 +
322 + static int
323 + addrot(                 /* compute rotation (x,y,z) => (xp,yp,zp) */
324 +        char *xfarg[],
325 +        FVECT xp,
326 +        FVECT yp,
327 +        FVECT zp
328 + )
329 + {
330 +        static char     bufs[3][16];
331 +        int     bn = 0;
332 +        char    **xfp = xfarg;
333 +        double  theta;
334 +
335 +        theta = atan2(yp[2], zp[2]);
336 +        if (!FEQ(theta,0.0)) {
337 +                *xfp++ = "-rx";
338 +                sprintf(bufs[bn], "%f", theta*(180./PI));
339 +                *xfp++ = bufs[bn++];
340 +        }
341 +        theta = asin(-xp[2]);
342 +        if (!FEQ(theta,0.0)) {
343 +                *xfp++ = "-ry";
344 +                sprintf(bufs[bn], " %f", theta*(180./PI));
345 +                *xfp++ = bufs[bn++];
346 +        }
347 +        theta = atan2(xp[1], xp[0]);
348 +        if (!FEQ(theta,0.0)) {
349 +                *xfp++ = "-rz";
350 +                sprintf(bufs[bn], "%f", theta*(180./PI));
351 +                *xfp++ = bufs[bn++];
352 +        }
353 +        *xfp = NULL;
354 +        return(xfp - xfarg);
355 + }
356 +
357 +
358 + int
359 + getBSDF_xfm(            /* compute BSDF orient. -> world orient. transform */
360 +        MAT4 xm,
361 +        FVECT nrm,
362 +        UpDir ud
363 + )
364 + {
365 +        char    *xfargs[7];
366 +        XF      myxf;
367 +        FVECT   updir, xdest, ydest;
368 +
369 +        updir[0] = updir[1] = updir[2] = 0.;
370 +        switch (ud) {
371 +        case UDzneg:
372 +                updir[2] = -1.;
373 +                break;
374 +        case UDyneg:
375 +                updir[1] = -1.;
376 +                break;
377 +        case UDxneg:
378 +                updir[0] = -1.;
379 +                break;
380 +        case UDxpos:
381 +                updir[0] = 1.;
382 +                break;
383 +        case UDypos:
384 +                updir[1] = 1.;
385 +                break;
386 +        case UDzpos:
387 +                updir[2] = 1.;
388 +                break;
389 +        case UDunknown:
390 +                error(WARNING, "unspecified up direction");
391 +                return(0);
392 +        }
393 +        fcross(xdest, updir, nrm);
394 +        if (normalize(xdest) == 0.0)
395 +                return(0);
396 +        fcross(ydest, nrm, xdest);
397 +        xf(&myxf, addrot(xfargs, xdest, ydest, nrm), xfargs);
398 +        copymat4(xm, myxf.xfm);
399 +        return(1);
400 + }
401 +
402 +
403 + void
404 + redistribute(           /* pass distarr ray sums through BSDF */
405 +        struct BSDF_data *b,
406 +        int nalt,
407 +        int nazi,
408 +        FVECT u,
409 +        FVECT v,
410 +        FVECT w,
411 +        MAT4 xm
412 + )
413 + {
414 +        int     nout = 0;
415 +        MAT4    mymat, inmat;
416 +        COLORV  *idist;
417 +        COLORV  *cp, *csum;
418 +        FVECT   dv;
419 +        double  wt;
420 +        int     i, j, k, o;
421 +        COLOR   col, cinc;
422 +                                        /* copy incoming distribution */
423 +        if (b->ninc != distsiz)
424 +                error(INTERNAL, "error 1 in redistribute");
425 +        idist = (COLORV *)malloc(sizeof(COLOR)*distsiz);
426 +        if (idist == NULL)
427 +                error(SYSTEM, "out of memory in redistribute");
428 +        memcpy(idist, distarr, sizeof(COLOR)*distsiz);
429 +                                        /* compose direction transform */
430 +        for (i = 3; i--; ) {
431 +                mymat[i][0] = u[i];
432 +                mymat[i][1] = v[i];
433 +                mymat[i][2] = w[i];
434 +                mymat[i][3] = 0.;
435 +        }
436 +        mymat[3][0] = mymat[3][1] = mymat[3][2] = 0.;
437 +        mymat[3][3] = 1.;
438 +        if (xm != NULL)
439 +                multmat4(mymat, xm, mymat);
440 +        for (i = 3; i--; ) {            /* make sure it's normalized */
441 +                wt = 1./sqrt(   mymat[0][i]*mymat[0][i] +
442 +                                mymat[1][i]*mymat[1][i] +
443 +                                mymat[2][i]*mymat[2][i] );
444 +                for (j = 3; j--; )
445 +                        mymat[j][i] *= wt;
446 +        }
447 +        if (!invmat4(inmat, mymat))     /* need inverse as well */
448 +                error(INTERNAL, "cannot invert BSDF transform");
449 +        newdist(nalt*nazi);             /* resample distribution */
450 +        for (i = b->ninc; i--; ) {
451 +                getBSDF_incvec(dv, b, i);       /* compute incident irrad. */
452 +                multv3(dv, dv, mymat);
453 +                if (dv[2] < 0.0) dv[2] = -dv[2];
454 +                wt = getBSDF_incohm(b, i);
455 +                wt *= dv[2];                    /* solid_angle*cosine(theta) */
456 +                cp = &idist[3*i];
457 +                copycolor(cinc, cp);
458 +                scalecolor(cinc, wt);
459 +                for (k = nalt; k--; )           /* loop over distribution */
460 +                    for (j = nazi; j--; ) {
461 +                        flatdir(dv, (k + .5)/nalt, (double)j/nazi);
462 +                        multv3(dv, dv, inmat);
463 +                                                /* evaluate BSDF @ outgoing */
464 +                        o = getBSDF_outndx(b, dv);
465 +                        if (o < 0) {
466 +                                nout++;
467 +                                continue;
468 +                        }
469 +                        wt = BSDF_visible(b, i, o);
470 +                        copycolor(col, cinc);
471 +                        scalecolor(col, wt);
472 +                        csum = &distarr[3*(k*nazi + j)];
473 +                        addcolor(csum, col);    /* sum into distribution */
474 +                    }
475 +        }
476 +        free(idist);                    /* free temp space */
477 +        if (nout) {
478 +                sprintf(errmsg, "missing %.1f%% of BSDF directions",
479 +                                100.*nout/(b->ninc*nalt*nazi));
480 +                error(WARNING, errmsg);
481 +        }
482   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines