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.2 by greg, Fri Sep 21 05:53:21 2007 UTC vs.
Revision 2.10 by greg, Thu Mar 27 01:40:30 2008 UTC

# Line 9 | Line 9 | static const char RCSid[] = "$Id$";
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[64];               /* 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 MAXABASES       3               /* limit on defined bases */
25 +
26 + ANGLE_BASIS     abase_list[MAXABASES] = {
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 + static int      nabases = 3;    /* current number of defined bases */
61 +
62 +
63 + static int
64 + ab_getvec(              /* get vector for this angle basis index */
65 +        FVECT v,
66 +        int ndx,
67 +        void *p
68 + )
69 + {
70 +        ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
71 +        int     li;
72 +        double  alt, azi, d;
73 +        
74 +        if ((ndx < 0) | (ndx >= ab->nangles))
75 +                return(0);
76 +        for (li = 0; ndx >= ab->lat[li].nphis; li++)
77 +                ndx -= ab->lat[li].nphis;
78 +        alt = PI/180.*0.5*(ab->lat[li].tmin + ab->lat[li+1].tmin);
79 +        azi = 2.*PI*ndx/ab->lat[li].nphis;
80 +        d = sin(alt);
81 +        v[0] = cos(azi)*d;
82 +        v[1] = sin(azi)*d;
83 +        v[2] = cos(alt);
84 +        return(1);
85 + }
86 +
87 +
88 + static int
89 + ab_getndx(              /* get index corresponding to the given vector */
90 +        FVECT v,
91 +        void *p
92 + )
93 + {
94 +        ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
95 +        int     li, ndx;
96 +        double  alt, azi, d;
97 +
98 +        if ((v[2] < -1.0) | (v[2] > 1.0))
99 +                return(-1);
100 +        alt = 180.0/PI*acos(v[2]);
101 +        azi = 180.0/PI*atan2(v[1], v[0]);
102 +        if (azi < 0.0) azi += 360.0;
103 +        for (li = 1; ab->lat[li].tmin <= alt; li++)
104 +                if (!ab->lat[li].nphis)
105 +                        return(-1);
106 +        --li;
107 +        ndx = (int)((1./360.)*azi*ab->lat[li].nphis + 0.5);
108 +        if (ndx >= ab->lat[li].nphis) ndx = 0;
109 +        while (li--)
110 +                ndx += ab->lat[li].nphis;
111 +        return(ndx);
112 + }
113 +
114 +
115 + static double
116 + ab_getohm(              /* get solid angle for this angle basis index */
117 +        int ndx,
118 +        void *p
119 + )
120 + {
121 +        ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
122 +        int     li;
123 +        double  tdia, pdia;
124 +        
125 +        if ((ndx < 0) | (ndx >= ab->nangles))
126 +                return(0);
127 +        for (li = 0; ndx >= ab->lat[li].nphis; li++)
128 +                ndx -= ab->lat[li].nphis;
129 +        if (ab->lat[li].nphis == 1) {           /* special case */
130 +                if (ab->lat[li].tmin > FTINY)
131 +                        error(USER, "unsupported BSDF coordinate system");
132 +                tdia = PI/180. * ab->lat[li+1].tmin;
133 +                return(PI*tdia*tdia);
134 +        }
135 +        tdia = PI/180.*(ab->lat[li+1].tmin - ab->lat[li].tmin);
136 +        tdia *= sin(PI/180.*(ab->lat[li].tmin + ab->lat[li+1].tmin));
137 +        pdia = 2.*PI/ab->lat[li].nphis;
138 +        return(tdia*pdia);
139 + }
140 +
141 +
142 + static int
143 + ab_getvecR(             /* get reverse vector for this angle basis index */
144 +        FVECT v,
145 +        int ndx,
146 +        void *p
147 + )
148 + {
149 +        if (!ab_getvec(v, ndx, p))
150 +                return(0);
151 +
152 +        v[0] = -v[0];
153 +        v[1] = -v[1];
154 +
155 +        return(1);
156 + }
157 +
158 +
159 + static int
160 + ab_getndxR(             /* get index corresponding to the reverse vector */
161 +        FVECT v,
162 +        void *p
163 + )
164 + {
165 +        FVECT  v2;
166 +        
167 +        v2[0] = -v[0];
168 +        v2[1] = -v[1];
169 +        v2[2] = v[2];
170 +
171 +        return ab_getndx(v2, p);
172 + }
173 +
174 +
175 + static void
176 + load_bsdf_data(         /* load BSDF distribution for this wavelength */
177 +        struct BSDF_data *dp,
178 +        ezxml_t wdb
179 + )
180 + {
181 +        char  *cbasis = ezxml_txt(ezxml_child(wdb,"ColumnAngleBasis"));
182 +        char  *rbasis = ezxml_txt(ezxml_child(wdb,"RowAngleBasis"));
183 +        char  *sdata;
184 +        int  i;
185 +        
186 +        if ((cbasis == NULL) | (rbasis == NULL)) {
187 +                error(WARNING, "missing column/row basis for BSDF");
188 +                return;
189 +        }
190 +        /* XXX need to add routines for loading in foreign bases */
191 +        for (i = nabases; i--; )
192 +                if (!strcmp(cbasis, abase_list[i].name)) {
193 +                        dp->ninc = abase_list[i].nangles;
194 +                        dp->ib_priv = (void *)&abase_list[i];
195 +                        dp->ib_vec = ab_getvec;
196 +                        dp->ib_ndx = ab_getndx;
197 +                        dp->ib_ohm = ab_getohm;
198 +                        break;
199 +                }
200 +        if (i < 0) {
201 +                sprintf(errmsg, "unsupported ColumnAngleBasis '%s'", cbasis);
202 +                error(WARNING, errmsg);
203 +                return;
204 +        }
205 +        for (i = nabases; i--; )
206 +                if (!strcmp(rbasis, abase_list[i].name)) {
207 +                        dp->nout = abase_list[i].nangles;
208 +                        dp->ob_priv = (void *)&abase_list[i];
209 +                        dp->ob_vec = ab_getvecR;
210 +                        dp->ob_ndx = ab_getndxR;
211 +                        dp->ob_ohm = ab_getohm;
212 +                        break;
213 +                }
214 +        if (i < 0) {
215 +                sprintf(errmsg, "unsupported RowAngleBasis '%s'", cbasis);
216 +                error(WARNING, errmsg);
217 +                return;
218 +        }
219 +                                /* read BSDF data */
220 +        sdata  = ezxml_txt(ezxml_child(wdb,"ScatteringData"));
221 +        if (sdata == NULL) {
222 +                error(WARNING, "missing BSDF ScatteringData");
223 +                return;
224 +        }
225 +        dp->bsdf = (float *)malloc(sizeof(float)*dp->ninc*dp->nout);
226 +        if (dp->bsdf == NULL)
227 +                error(SYSTEM, "out of memory in load_bsdf_data");
228 +        for (i = 0; i < dp->ninc*dp->nout; i++) {
229 +                char  *sdnext = fskip(sdata);
230 +                if (sdnext == NULL) {
231 +                        error(WARNING, "bad/missing BSDF ScatteringData");
232 +                        free(dp->bsdf); dp->bsdf = NULL;
233 +                        return;
234 +                }
235 +                while (*sdnext && isspace(*sdnext))
236 +                        sdnext++;
237 +                if (*sdnext == ',') sdnext++;
238 +                dp->bsdf[i] = atof(sdata);
239 +                sdata = sdnext;
240 +        }
241 +        while (isspace(*sdata))
242 +                sdata++;
243 +        if (*sdata) {
244 +                sprintf(errmsg, "%d extra characters after BSDF ScatteringData",
245 +                                strlen(sdata));
246 +                error(WARNING, errmsg);
247 +        }
248 + }
249 +
250 +
251   struct BSDF_data *
252   load_BSDF(              /* load BSDF data from file */
253          char *fname
254   )
255   {
256          char                    *path;
257 <        ezxml_t                 fl, wld, wdb;
257 >        ezxml_t                 fl, wtl, wld, wdb;
258          struct BSDF_data        *dp;
259          
260          path = getpath(fname, getrlibpath(), R_OK);
# Line 31 | Line 269 | load_BSDF(             /* load BSDF data from file */
269                  error(WARNING, errmsg);
270                  return(NULL);
271          }
272 <        dp = (struct BSDF_data *)malloc(sizeof(struct BSDF_data));
273 <        if (dp == NULL)
274 <                goto memerr;
275 <        for (wld = ezxml_child(fl, "WavelengthData");
276 <                                fl != NULL; fl = fl->next) {
277 <                if (strcmp(ezxml_child(wld, "Wavelength")->txt, "Visible"))
272 >        if (ezxml_error(fl)[0]) {
273 >                sprintf(errmsg, "BSDF \"%s\" %s", path, ezxml_error(fl));
274 >                error(WARNING, errmsg);
275 >                ezxml_free(fl);
276 >                return(NULL);
277 >        }
278 >        if (strcmp(ezxml_name(fl), "WindowElement")) {
279 >                sprintf(errmsg,
280 >                        "BSDF \"%s\": top level node not 'WindowElement'",
281 >                                path);
282 >                error(WARNING, errmsg);
283 >                ezxml_free(fl);
284 >                return(NULL);
285 >        }
286 >        wtl = ezxml_child(ezxml_child(fl, "Optical"), "Layer");
287 >        dp = (struct BSDF_data *)calloc(1, sizeof(struct BSDF_data));
288 >        for (wld = ezxml_child(wtl, "WavelengthData");
289 >                                wld != NULL; wld = wld->next) {
290 >                if (strcmp(ezxml_txt(ezxml_child(wld,"Wavelength")), "Visible"))
291                          continue;
292                  wdb = ezxml_child(wld, "WavelengthDataBlock");
293                  if (wdb == NULL) continue;
294 <                if (strcmp(ezxml_child(wdb, "WavelengthDataDirection")->txt,
294 >                if (strcmp(ezxml_txt(ezxml_child(wdb,"WavelengthDataDirection")),
295                                          "Transmission Front"))
296                          continue;
297 +                load_bsdf_data(dp, wdb);        /* load front BTDF */
298 +                break;                          /* ignore the rest */
299          }
300 <        /* etc... */
301 <        ezxml_free(fl);
300 >        ezxml_free(fl);                         /* done with XML file */
301 >        if (dp->bsdf == NULL) {
302 >                sprintf(errmsg, "bad/missing BTDF data in \"%s\"", path);
303 >                error(WARNING, errmsg);
304 >                free_BSDF(dp);
305 >                dp = NULL;
306 >        }
307          return(dp);
50 – memerr:
51 –        error(SYSTEM, "out of memory in load_BSDF");
52 –        return NULL;    /* pro forma return */
308   }
309  
310  
# Line 60 | Line 315 | free_BSDF(             /* free BSDF data structure */
315   {
316          if (b == NULL)
317                  return;
318 <        free(b->inc_dir);
319 <        free(b->inc_rad);
65 <        free(b->out_dir);
66 <        free(b->out_rad);
67 <        free(b->bsdf);
318 >        if (b->bsdf != NULL)
319 >                free(b->bsdf);
320          free(b);
321   }
322  
323  
324 < void
324 > int
325   r_BSDF_incvec(          /* compute random input vector at given location */
326          FVECT v,
327          struct BSDF_data *b,
# Line 82 | Line 334 | r_BSDF_incvec(         /* compute random input vector at give
334          double  rad;
335          int     j;
336          
337 <        getBSDF_incvec(v, b, i);
338 <        rad = getBSDF_incrad(b, i);
337 >        if (!getBSDF_incvec(v, b, i))
338 >                return(0);
339 >        rad = sqrt(getBSDF_incohm(b, i) / PI);
340          multisamp(pert, 3, rv);
341          for (j = 0; j < 3; j++)
342                  v[j] += rad*(2.*pert[j] - 1.);
343          if (xm != NULL)
344                  multv3(v, v, xm);
345 <        normalize(v);
345 >        return(normalize(v) != 0.0);
346   }
347  
348  
349 < void
349 > int
350   r_BSDF_outvec(          /* compute random output vector at given location */
351          FVECT v,
352          struct BSDF_data *b,
# Line 106 | Line 359 | r_BSDF_outvec(         /* compute random output vector at giv
359          double  rad;
360          int     j;
361          
362 <        getBSDF_outvec(v, b, o);
363 <        rad = getBSDF_outrad(b, o);
362 >        if (!getBSDF_outvec(v, b, o))
363 >                return(0);
364 >        rad = sqrt(getBSDF_outohm(b, o) / PI);
365          multisamp(pert, 3, rv);
366          for (j = 0; j < 3; j++)
367                  v[j] += rad*(2.*pert[j] - 1.);
368          if (xm != NULL)
369                  multv3(v, v, xm);
370 <        normalize(v);
370 >        return(normalize(v) != 0.0);
371   }
372  
373  
# Line 132 | Line 386 | addrot(                        /* compute rotation (x,y,z) => (xp,yp,zp) */
386          char    **xfp = xfarg;
387          double  theta;
388  
389 +        if (yp[2]*yp[2] + zp[2]*zp[2] < 2.*FTINY*FTINY) {
390 +                /* Special case for X' along Z-axis */
391 +                theta = -atan2(yp[0], yp[1]);
392 +                *xfp++ = "-ry";
393 +                *xfp++ = xp[2] < 0.0 ? "90" : "-90";
394 +                *xfp++ = "-rz";
395 +                sprintf(bufs[bn], "%f", theta*(180./PI));
396 +                *xfp++ = bufs[bn++];
397 +                return(xfp - xfarg);
398 +        }
399          theta = atan2(yp[2], zp[2]);
400          if (!FEQ(theta,0.0)) {
401                  *xfp++ = "-rx";
# Line 156 | Line 420 | addrot(                        /* compute rotation (x,y,z) => (xp,yp,zp) */
420  
421  
422   int
423 < getBSDF_xfm(            /* compute transform for the given surface */
423 > getBSDF_xfm(            /* compute BSDF orient. -> world orient. transform */
424          MAT4 xm,
425          FVECT nrm,
426          UpDir ud
# Line 187 | Line 451 | getBSDF_xfm(           /* compute transform for the given surfa
451                  updir[2] = 1.;
452                  break;
453          case UDunknown:
190 –                error(WARNING, "unspecified up direction");
454                  return(0);
455          }
456          fcross(xdest, updir, nrm);
# Line 211 | Line 474 | redistribute(          /* pass distarr ray sums through BSDF *
474          MAT4 xm
475   )
476   {
477 <        MAT4    mymat;
478 <        COLORV  *outarr;
479 <        float   *inpcoef;
477 >        int     nout = 0;
478 >        MAT4    mymat, inmat;
479 >        COLORV  *idist;
480          COLORV  *cp, *csum;
218 –        uint16  *distcnt;
481          FVECT   dv;
482 <        double  oom, wt;
483 <        int     i, j, o;
484 <        int     cnt;
485 <        COLOR   col;
486 <                                        /* allocate temporary memory */
487 <        outarr = (COLORV *)malloc(b->nout * sizeof(COLOR));
488 <        distcnt = (uint16 *)calloc(nalt*nazi, sizeof(uint16));
489 <        inpcoef = (float *)malloc(b->ninc * sizeof(float));
228 <        if ((outarr == NULL) | (distcnt == NULL) | (inpcoef == NULL))
482 >        double  wt;
483 >        int     i, j, k, o;
484 >        COLOR   col, cinc;
485 >                                        /* copy incoming distribution */
486 >        if (b->ninc != distsiz)
487 >                error(INTERNAL, "error 1 in redistribute");
488 >        idist = (COLORV *)malloc(sizeof(COLOR)*distsiz);
489 >        if (idist == NULL)
490                  error(SYSTEM, "out of memory in redistribute");
491 <                                        /* compose matrix */
491 >        memcpy(idist, distarr, sizeof(COLOR)*distsiz);
492 >                                        /* compose direction transform */
493          for (i = 3; i--; ) {
494                  mymat[i][0] = u[i];
495                  mymat[i][1] = v[i];
# Line 245 | Line 507 | redistribute(          /* pass distarr ray sums through BSDF *
507                  for (j = 3; j--; )
508                          mymat[j][i] *= wt;
509          }
510 <                                        /* pass through BSDF */
511 <        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);
510 >        if (!invmat4(inmat, mymat))     /* need inverse as well */
511 >                error(INTERNAL, "cannot invert BSDF transform");
512          newdist(nalt*nazi);             /* resample distribution */
513 <        for (o = b->nout; o--; ) {
514 <                getBSDF_outvec(dv, b, o);
513 >        for (i = b->ninc; i--; ) {
514 >                getBSDF_incvec(dv, b, i);       /* compute incident irrad. */
515                  multv3(dv, dv, mymat);
516 <                j = (.5 + atan2(dv[1],dv[0])*(.5/PI))*nazi + .5;
517 <                if (j >= nazi) j = 0;
518 <                i = (0.9999 - dv[2]*dv[2])*nalt;
519 <                csum = &distarr[3*(i*nazi + j)];
520 <                cp = &outarr[3*o];
521 <                addcolor(csum, cp);
522 <                ++distcnt[i*nazi + j];
523 <        }
524 <        free(outarr);
525 <                                        /* fill in missing bits */
526 <        for (i = nalt; i--; )
527 <            for (j = nazi; j--; ) {
528 <                int     ii, jj, alt, azi;
529 <                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)
516 >                if (dv[2] < 0.0) dv[2] = -dv[2];
517 >                wt = getBSDF_incohm(b, i);
518 >                wt *= dv[2];                    /* solid_angle*cosine(theta) */
519 >                cp = &idist[3*i];
520 >                copycolor(cinc, cp);
521 >                scalecolor(cinc, wt);
522 >                for (k = nalt; k--; )           /* loop over distribution */
523 >                    for (j = nazi; j--; ) {
524 >                        flatdir(dv, (k + .5)/nalt, (double)j/nazi);
525 >                        multv3(dv, dv, inmat);
526 >                                                /* evaluate BSDF @ outgoing */
527 >                        o = getBSDF_outndx(b, dv);
528 >                        if (o < 0) {
529 >                                nout++;
530                                  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];
531                          }
532 +                        wt = BSDF_value(b, i, o);
533 +                        copycolor(col, cinc);
534 +                        scalecolor(col, wt);
535 +                        csum = &distarr[3*(k*nazi + j)];
536 +                        addcolor(csum, col);    /* sum into distribution */
537                      }
538 <                wt = 1./cnt;
539 <                scalecolor(csum, wt);
540 <            }
541 <                                        /* finish averages */
542 <        for (i = nalt; i--; )
543 <            for (j = nazi; j--; ) {
544 <                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);
538 >        }
539 >        free(idist);                    /* free temp space */
540 >        if (nout) {
541 >                sprintf(errmsg, "missing %.1f%% of BSDF directions",
542 >                                100.*nout/(b->ninc*nalt*nazi));
543 >                error(WARNING, errmsg);
544 >        }
545   }

Diff Legend

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