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

Comparing ray/src/common/bsdf.c (file contents):
Revision 2.1 by greg, Wed Jun 17 20:41:47 2009 UTC vs.
Revision 2.2 by greg, Fri Jun 19 06:49:25 2009 UTC

# Line 71 | Line 71 | ab_getvec(             /* get vector for this angle basis index *
71   {
72          ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
73          int     li;
74 <        double  alt, azi, d;
74 >        double  pol, azi, d;
75          
76          if ((ndx < 0) | (ndx >= ab->nangles))
77                  return(0);
78          for (li = 0; ndx >= ab->lat[li].nphis; li++)
79                  ndx -= ab->lat[li].nphis;
80 <        alt = PI/180.*0.5*(ab->lat[li].tmin + ab->lat[li+1].tmin);
80 >        pol = PI/180.*0.5*(ab->lat[li].tmin + ab->lat[li+1].tmin);
81          azi = 2.*PI*ndx/ab->lat[li].nphis;
82 <        v[2] = d = cos(alt);
83 <        d = sqrt(1. - d*d);     /* sin(alt) */
82 >        v[2] = d = cos(pol);
83 >        d = sqrt(1. - d*d);     /* sin(pol) */
84          v[0] = cos(azi)*d;
85          v[1] = sin(azi)*d;
86          return(1);
# Line 95 | Line 95 | ab_getndx(             /* get index corresponding to the given ve
95   {
96          ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
97          int     li, ndx;
98 <        double  alt, azi, d;
98 >        double  pol, azi, d;
99  
100          if ((v[2] < -1.0) | (v[2] > 1.0))
101                  return(-1);
102 <        alt = 180.0/PI*acos(v[2]);
102 >        pol = 180.0/PI*acos(v[2]);
103          azi = 180.0/PI*atan2(v[1], v[0]);
104          if (azi < 0.0) azi += 360.0;
105 <        for (li = 1; ab->lat[li].tmin <= alt; li++)
105 >        for (li = 1; ab->lat[li].tmin <= pol; li++)
106                  if (!ab->lat[li].nphis)
107                          return(-1);
108          --li;
# Line 254 | Line 254 | check_bsdf_data(       /* check that BSDF data is sane */
254          struct BSDF_data *dp
255   )
256   {
257 <        double *        omega_arr;
258 <        double          dom, hemi_total;
257 >        double          *omega_iarr, *omega_oarr;
258 >        double          dom, contrib, hemi_total;
259          int             nneg;
260 +        FVECT           v;
261          int             i, o;
262  
263          if (dp == NULL || dp->bsdf == NULL)
264                  return(0);
265 <        omega_arr = (double *)calloc(dp->nout, sizeof(double));
266 <        if (omega_arr == NULL)
265 >        omega_iarr = (double *)calloc(dp->ninc, sizeof(double));
266 >        omega_oarr = (double *)calloc(dp->nout, sizeof(double));
267 >        if ((omega_iarr == NULL) | (omega_oarr == NULL))
268                  error(SYSTEM, "out of memory in check_bsdf_data");
269 +                                        /* incoming projected solid angles */
270          hemi_total = .0;
271 +        for (i = dp->ninc; i--; ) {
272 +                dom = getBSDF_incohm(dp,i);
273 +                if (dom <= .0) {
274 +                        error(WARNING, "zero/negative incoming solid angle");
275 +                        continue;
276 +                }
277 +                if (!getBSDF_incvec(v,dp,i) || v[2] > FTINY) {
278 +                        error(WARNING, "illegal incoming BSDF direction");
279 +                        free(omega_iarr); free(omega_oarr);
280 +                        return(0);
281 +                }
282 +                hemi_total += omega_iarr[i] = dom * -v[2];
283 +        }
284 +        if ((hemi_total > 1.02*PI) | (hemi_total < 0.98*PI)) {
285 +                sprintf(errmsg, "incoming BSDF hemisphere off by %.1f%%",
286 +                                100.*(hemi_total/PI - 1.));
287 +                error(WARNING, errmsg);
288 +        }
289 +        dom = PI / hemi_total;          /* fix normalization */
290 +        for (i = dp->ninc; i--; )
291 +                omega_iarr[i] *= dom;
292 +                                        /* outgoing projected solid angles */
293 +        hemi_total = .0;
294          for (o = dp->nout; o--; ) {
269                FVECT   v;
295                  dom = getBSDF_outohm(dp,o);
296                  if (dom <= .0) {
297 <                        error(WARNING, "zero/negative solid angle");
297 >                        error(WARNING, "zero/negative outgoing solid angle");
298                          continue;
299                  }
300                  if (!getBSDF_outvec(v,dp,o) || v[2] < -FTINY) {
301                          error(WARNING, "illegal outgoing BSDF direction");
302 <                        free(omega_arr);
302 >                        free(omega_iarr); free(omega_oarr);
303                          return(0);
304                  }
305 <                hemi_total += omega_arr[o] = dom*v[2];
305 >                hemi_total += omega_oarr[o] = dom * v[2];
306          }
307          if ((hemi_total > 1.02*PI) | (hemi_total < 0.98*PI)) {
308                  sprintf(errmsg, "outgoing BSDF hemisphere off by %.1f%%",
309                                  100.*(hemi_total/PI - 1.));
310                  error(WARNING, errmsg);
311          }
312 <        dom = PI / hemi_total;          /* normalize solid angles */
312 >        dom = PI / hemi_total;          /* fix normalization */
313          for (o = dp->nout; o--; )
314 <                omega_arr[o] *= dom;
315 <        nneg = 0;
316 <        for (i = dp->ninc; i--; ) {
314 >                omega_oarr[o] *= dom;
315 >        nneg = 0;                       /* check outgoing totals */
316 >        for (i = 0; i < dp->ninc; i++) {
317                  hemi_total = .0;
318                  for (o = dp->nout; o--; ) {
319                          double  f = BSDF_value(dp,i,o);
320 <                        if (f > .0)
321 <                                hemi_total += f*omega_arr[o];
322 <                        else if (f < -FTINY)
323 <                                ++nneg;
320 >                        if (f >= .0)
321 >                                hemi_total += f*omega_oarr[o];
322 >                        else {
323 >                                nneg += (f < -FTINY);
324 >                                BSDF_value(dp,i,o) = .0f;
325 >                        }
326                  }
327                  if (hemi_total > 1.02) {
328 <                        sprintf(errmsg, "BSDF direction passes %.1f%% of light",
329 <                                        100.*hemi_total);
328 >                        sprintf(errmsg,
329 >                        "incoming BSDF direction %d passes %.1f%% of light",
330 >                                        i, 100.*hemi_total);
331                          error(WARNING, errmsg);
332                  }
333          }
334 <        free(omega_arr);
335 <        if (nneg > 0) {
308 <                sprintf(errmsg, "%d negative BSDF values", nneg);
334 >        if (nneg) {
335 >                sprintf(errmsg, "%d negative BSDF values (ignored)", nneg);
336                  error(WARNING, errmsg);
310                return(0);
337          }
338 +                                        /* reverse roles and check again */
339 +        for (o = 0; o < dp->nout; o++) {
340 +                hemi_total = .0;
341 +                for (i = dp->ninc; i--; )
342 +                        hemi_total += BSDF_value(dp,i,o) * omega_iarr[i];
343 +
344 +                if (hemi_total > 1.02) {
345 +                        sprintf(errmsg,
346 +                        "outgoing BSDF direction %d collects %.1f%% of light",
347 +                                        o, 100.*hemi_total);
348 +                        error(WARNING, errmsg);
349 +                }
350 +        }
351 +        free(omega_iarr); free(omega_oarr);
352          return(1);
353   }
354  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines