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

Comparing ray/src/common/bsdf_m.c (file contents):
Revision 3.4 by greg, Sat Feb 19 01:48:59 2011 UTC vs.
Revision 3.32 by greg, Wed Apr 8 02:41:02 2015 UTC

# Line 11 | Line 11 | static const char RCSid[] = "$Id$";
11   *
12   */
13  
14 + #define _USE_MATH_DEFINES
15   #include "rtio.h"
15 – #include <stdlib.h>
16   #include <math.h>
17   #include <ctype.h>
18   #include "ezxml.h"
# Line 28 | Line 28 | static const char RCSid[] = "$Id$";
28   #define RC_INTERR       (-4)
29   #define RC_MEMERR       (-5)
30  
31 < #define MAXLATS         46              /* maximum number of latitudes */
32 <
33 < /* BSDF angle specification */
34 < typedef struct {
35 <        char    name[64];               /* basis name */
36 <        int     nangles;                /* total number of directions */
37 <        struct {
38 <                float   tmin;                   /* starting theta */
39 <                int     nphis;                  /* number of phis (0 term) */
40 <        } lat[MAXLATS+1];               /* latitudes */
41 < } ANGLE_BASIS;
42 <
43 < #define MAXABASES       7               /* limit on defined bases */
44 <
45 < static ANGLE_BASIS      abase_list[MAXABASES] = {
31 > ANGLE_BASIS     abase_list[MAXABASES] = {
32          {
33                  "LBNL/Klems Full", 145,
34 <                { {-5., 1},
34 >                { {0., 1},
35                  {5., 8},
36                  {15., 16},
37                  {25., 20},
# Line 57 | Line 43 | static ANGLE_BASIS     abase_list[MAXABASES] = {
43                  {90., 0} }
44          }, {
45                  "LBNL/Klems Half", 73,
46 <                { {-6.5, 1},
46 >                { {0., 1},
47                  {6.5, 8},
48                  {19.5, 12},
49                  {32.5, 16},
# Line 67 | Line 53 | static ANGLE_BASIS     abase_list[MAXABASES] = {
53                  {90., 0} }
54          }, {
55                  "LBNL/Klems Quarter", 41,
56 <                { {-9., 1},
56 >                { {0., 1},
57                  {9., 8},
58                  {27., 12},
59                  {46., 12},
# Line 76 | Line 62 | static ANGLE_BASIS     abase_list[MAXABASES] = {
62          }
63   };
64  
65 < static int      nabases = 3;    /* current number of defined bases */
65 > int             nabases = 3;            /* current number of defined bases */
66  
67 + C_COLOR         mtx_RGB_prim[3];        /* our RGB primaries  */
68 + float           mtx_RGB_coef[3];        /* corresponding Y coefficients */
69 +
70 + enum {mtx_Y, mtx_X, mtx_Z};             /* matrix components (mtx_Y==0) */
71 +
72 + /* check if two real values are near enough to equal */
73   static int
74   fequal(double a, double b)
75   {
76 <        if (b != .0)
76 >        if (b != 0)
77                  a = a/b - 1.;
78          return (a <= 1e-6) & (a >= -1e-6);
79   }
80  
81 < /* returns the name of the given tag */
90 < #ifdef ezxml_name
91 < #undef ezxml_name
92 < static char *
93 < ezxml_name(ezxml_t xml)
94 < {
95 <        if (xml == NULL)
96 <                return NULL;
97 <        return xml->name;
98 < }
99 < #endif
100 <
101 < /* returns the given tag's character content or empty string if none */
102 < #ifdef ezxml_txt
103 < #undef ezxml_txt
104 < static char *
105 < ezxml_txt(ezxml_t xml)
106 < {
107 <        if (xml == NULL)
108 <                return "";
109 <        return xml->txt;
110 < }
111 < #endif
112 <
113 < /* Convert error to standard BSDF code */
81 > /* convert error to standard BSDF code */
82   static SDError
83   convert_errcode(int ec)
84   {
# Line 131 | Line 99 | convert_errcode(int ec)
99          return SDEunknown;
100   }
101  
102 < /* Allocate a BSDF matrix of the given size */
102 > /* allocate a BSDF matrix of the given size */
103   static SDMat *
104   SDnewMatrix(int ni, int no)
105   {
# Line 155 | Line 123 | SDnewMatrix(int ni, int no)
123   }
124  
125   /* Free a BSDF matrix */
126 < #define SDfreeMatrix            free
126 > void
127 > SDfreeMatrix(void *ptr)
128 > {
129 >        SDMat   *mp = (SDMat *)ptr;
130  
131 < /* get vector for this angle basis index */
132 < static int
133 < ab_getvec(FVECT v, int ndx, double randX, void *p)
131 >        if (mp->chroma != NULL) free(mp->chroma);
132 >        free(ptr);
133 > }
134 >
135 > /* compute square of real value */
136 > static double sq(double x) { return x*x; }
137 >
138 > /* Get vector for this angle basis index (front exiting) */
139 > int
140 > fo_getvec(FVECT v, double ndxr, void *p)
141   {
142 <        ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
142 >        ANGLE_BASIS     *ab = (ANGLE_BASIS *)p;
143 >        int             ndx = (int)ndxr;
144 >        double          randX = ndxr - ndx;
145          double          rx[2];
146          int             li;
147 <        double          pol, azi, d;
147 >        double          azi, d;
148          
149 <        if ((ndx < 0) | (ndx >= ab->nangles))
149 >        if ((ndxr < 0) | (ndx >= ab->nangles))
150                  return RC_FAIL;
151          for (li = 0; ndx >= ab->lat[li].nphis; li++)
152                  ndx -= ab->lat[li].nphis;
153          SDmultiSamp(rx, 2, randX);
154 <        pol = M_PI/180.*( (1.-rx[0])*ab->lat[li].tmin +
155 <                                rx[0]*ab->lat[li+1].tmin );
154 >        d = (1. - randX)*sq(cos(M_PI/180.*ab->lat[li].tmin)) +
155 >                randX*sq(cos(M_PI/180.*ab->lat[li+1].tmin));
156 >        v[2] = d = sqrt(d);     /* cos(pol) */
157          azi = 2.*M_PI*(ndx + rx[1] - .5)/ab->lat[li].nphis;
177 –        v[2] = d = cos(pol);
158          d = sqrt(1. - d*d);     /* sin(pol) */
159          v[0] = cos(azi)*d;
160          v[1] = sin(azi)*d;
161          return RC_GOOD;
162   }
163  
164 < /* get index corresponding to the given vector */
165 < static int
166 < ab_getndx(const FVECT v, void *p)
164 > /* Get index corresponding to the given vector (front exiting) */
165 > int
166 > fo_getndx(const FVECT v, void *p)
167   {
168 <        ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
168 >        ANGLE_BASIS     *ab = (ANGLE_BASIS *)p;
169          int     li, ndx;
170 <        double  pol, azi, d;
170 >        double  pol, azi;
171  
172          if (v == NULL)
173                  return -1;
174 <        if ((v[2] < .0) | (v[2] > 1.0))
174 >        if ((v[2] < 0) | (v[2] > 1.))
175                  return -1;
176 <        pol = 180.0/M_PI*acos(v[2]);
176 >        pol = 180.0/M_PI*Acos(v[2]);
177          azi = 180.0/M_PI*atan2(v[1], v[0]);
178          if (azi < 0.0) azi += 360.0;
179          for (li = 1; ab->lat[li].tmin <= pol; li++)
# Line 207 | Line 187 | ab_getndx(const FVECT v, void *p)
187          return ndx;
188   }
189  
190 < /* compute square of real value */
191 < static double sq(double x) { return x*x; }
192 <
213 < /* get projected solid angle for this angle basis index */
214 < static double
215 < ab_getohm(int ndx, void *p)
190 > /* Get projected solid angle for this angle basis index (universal) */
191 > double
192 > io_getohm(int ndx, void *p)
193   {
194          static int      last_li = -1;
195          static double   last_ohm;
# Line 227 | Line 204 | ab_getohm(int ndx, void *p)
204          if (li == last_li)                      /* cached latitude? */
205                  return last_ohm;
206          last_li = li;
230 –        theta1 = M_PI/180. * ab->lat[li+1].tmin;
231 –        if (ab->lat[li].nphis == 1)             /* special case */
232 –                return last_ohm = M_PI*(1. - sq(cos(theta1)));
207          theta = M_PI/180. * ab->lat[li].tmin;
208 +        theta1 = M_PI/180. * ab->lat[li+1].tmin;
209          return last_ohm = M_PI*(sq(cos(theta)) - sq(cos(theta1))) /
210                                  (double)ab->lat[li].nphis;
211   }
212  
213 < /* get reverse vector for this angle basis index */
214 < static int
215 < ab_getvecR(FVECT v, int ndx, double randX, void *p)
213 > /* Get vector for this angle basis index (back incident) */
214 > int
215 > bi_getvec(FVECT v, double ndxr, void *p)
216   {
217 <        int     na = (*(ANGLE_BASIS *)p).nangles;
243 <
244 <        if (!ab_getvec(v, ndx, randX, p))
217 >        if (!fo_getvec(v, ndxr, p))
218                  return RC_FAIL;
219  
220          v[0] = -v[0];
# Line 251 | Line 224 | ab_getvecR(FVECT v, int ndx, double randX, void *p)
224          return RC_GOOD;
225   }
226  
227 < /* get index corresponding to the reverse vector */
228 < static int
229 < ab_getndxR(const FVECT v, void *p)
227 > /* Get index corresponding to the vector (back incident) */
228 > int
229 > bi_getndx(const FVECT v, void *p)
230   {
231          FVECT  v2;
232          
# Line 261 | Line 234 | ab_getndxR(const FVECT v, void *p)
234          v2[1] = -v[1];
235          v2[2] = -v[2];
236  
237 <        return ab_getndx(v2, p);
237 >        return fo_getndx(v2, p);
238   }
239  
240 + /* Get vector for this angle basis index (back exiting) */
241 + int
242 + bo_getvec(FVECT v, double ndxr, void *p)
243 + {
244 +        if (!fo_getvec(v, ndxr, p))
245 +                return RC_FAIL;
246 +
247 +        v[2] = -v[2];
248 +
249 +        return RC_GOOD;
250 + }
251 +
252 + /* Get index corresponding to the vector (back exiting) */
253 + int
254 + bo_getndx(const FVECT v, void *p)
255 + {
256 +        FVECT  v2;
257 +        
258 +        v2[0] = v[0];
259 +        v2[1] = v[1];
260 +        v2[2] = -v[2];
261 +
262 +        return fo_getndx(v2, p);
263 + }
264 +
265 + /* Get vector for this angle basis index (front incident) */
266 + int
267 + fi_getvec(FVECT v, double ndxr, void *p)
268 + {
269 +        if (!fo_getvec(v, ndxr, p))
270 +                return RC_FAIL;
271 +
272 +        v[0] = -v[0];
273 +        v[1] = -v[1];
274 +
275 +        return RC_GOOD;
276 + }
277 +
278 + /* Get index corresponding to the vector (front incident) */
279 + int
280 + fi_getndx(const FVECT v, void *p)
281 + {
282 +        FVECT  v2;
283 +        
284 +        v2[0] = -v[0];
285 +        v2[1] = -v[1];
286 +        v2[2] = v[2];
287 +
288 +        return fo_getndx(v2, p);
289 + }
290 +
291 + /* Get color or grayscale value for BSDF for the given direction pair */
292 + int
293 + mBSDF_color(float coef[], const SDMat *dp, int i, int o)
294 + {
295 +        C_COLOR cxy;
296 +
297 +        coef[0] = mBSDF_value(dp, i, o);
298 +        if (dp->chroma == NULL)
299 +                return 1;       /* grayscale */
300 +
301 +        c_decodeChroma(&cxy, dp->chroma[o*dp->ninc + i]);
302 +        c_toSharpRGB(&cxy, coef[0], coef);
303 +        coef[0] *= mtx_RGB_coef[0];
304 +        coef[1] *= mtx_RGB_coef[1];
305 +        coef[2] *= mtx_RGB_coef[2];
306 +        return 3;               /* RGB color */
307 + }
308 +
309   /* load custom BSDF angle basis */
310   static int
311   load_angle_basis(ezxml_t wab)
# Line 295 | Line 337 | load_angle_basis(ezxml_t wab)
337                                  ezxml_child(ezxml_child(wbb,
338                                          "ThetaBounds"), "UpperTheta")));
339                  if (!i)
340 <                        abase_list[nabases].lat[i].tmin =
299 <                                        -abase_list[nabases].lat[i+1].tmin;
340 >                        abase_list[nabases].lat[0].tmin = 0;
341                  else if (!fequal(atof(ezxml_txt(ezxml_child(ezxml_child(wbb,
342                                          "ThetaBounds"), "LowerTheta"))),
343                                  abase_list[nabases].lat[i].tmin)) {
344                          sprintf(SDerrorDetail, "Theta values disagree in '%s'",
345 <                                abname);
345 >                                        abname);
346                          return RC_DATERR;
347                  }
348                  abase_list[nabases].nangles +=
# Line 311 | Line 352 | load_angle_basis(ezxml_t wab)
352                                  (abase_list[nabases].lat[i].nphis == 1 &&
353                                  abase_list[nabases].lat[i].tmin > FTINY)) {
354                          sprintf(SDerrorDetail, "Illegal phi count in '%s'",
355 <                                abname);
355 >                                        abname);
356                          return RC_DATERR;
357                  }
358          }
# Line 346 | Line 387 | get_extrema(SDSpectralDF *df)
387          }
388          free(ohma);
389                                          /* need incoming solid angles, too? */
390 <        if (dp->ninc < dp->nout || dp->ib_ohm != dp->ob_ohm ||
350 <                                dp->ib_priv != dp->ob_priv) {
390 >        if ((dp->ib_ohm != dp->ob_ohm) | (dp->ib_priv != dp->ob_priv)) {
391                  double  ohm;
392                  for (i = dp->ninc; i--; )
393                          if ((ohm = mBSDF_incohm(dp,i)) < df->minProjSA)
# Line 358 | Line 398 | get_extrema(SDSpectralDF *df)
398  
399   /* load BSDF distribution for this wavelength */
400   static int
401 < load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc)
401 > load_bsdf_data(SDData *sd, ezxml_t wdb, int ct, int rowinc)
402   {
403          SDSpectralDF    *df;
404          SDMat           *dp;
# Line 367 | Line 407 | load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc)
407          int             i;
408                                          /* allocate BSDF component */
409          sdata = ezxml_txt(ezxml_child(wdb, "WavelengthDataDirection"));
410 +        if (!sdata)
411 +                return RC_FAIL;
412 +        /*
413 +         * Remember that front and back are reversed from WINDOW 6 orientations
414 +         */
415          if (!strcasecmp(sdata, "Transmission Front")) {
416 <                if (sd->tf != NULL)
372 <                        SDfreeSpectralDF(sd->tf);
373 <                if ((sd->tf = SDnewSpectralDF(1)) == NULL)
416 >                if (sd->tb == NULL && (sd->tb = SDnewSpectralDF(3)) == NULL)
417                          return RC_MEMERR;
418 +                df = sd->tb;
419 +        } else if (!strcasecmp(sdata, "Transmission Back")) {
420 +                if (sd->tf == NULL && (sd->tf = SDnewSpectralDF(3)) == NULL)
421 +                        return RC_MEMERR;
422                  df = sd->tf;
423          } else if (!strcasecmp(sdata, "Reflection Front")) {
424 <                if (sd->rf != NULL)
378 <                        SDfreeSpectralDF(sd->rf);
379 <                if ((sd->rf = SDnewSpectralDF(1)) == NULL)
424 >                if (sd->rb == NULL && (sd->rb = SDnewSpectralDF(3)) == NULL)
425                          return RC_MEMERR;
426 <                df = sd->rf;
426 >                df = sd->rb;
427          } else if (!strcasecmp(sdata, "Reflection Back")) {
428 <                if (sd->rb != NULL)
384 <                        SDfreeSpectralDF(sd->rb);
385 <                if ((sd->rb = SDnewSpectralDF(1)) == NULL)
428 >                if (sd->rf == NULL && (sd->rf = SDnewSpectralDF(3)) == NULL)
429                          return RC_MEMERR;
430 <                df = sd->rb;
430 >                df = sd->rf;
431          } else
432                  return RC_FAIL;
433 <        /* XXX should also check "ScatteringDataType" for consistency? */
433 >                                        /* free previous matrix if any */
434 >        if (df->comp[ct].dist != NULL) {
435 >                SDfreeMatrix(df->comp[ct].dist);
436 >                df->comp[ct].dist = NULL;
437 >        }
438                                          /* get angle bases */
439          sdata = ezxml_txt(ezxml_child(wdb,"ColumnAngleBasis"));
440          if (!sdata || !*sdata) {
# Line 419 | Line 466 | load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc)
466          dp = SDnewMatrix(abase_list[inbi].nangles, abase_list[outbi].nangles);
467          if (dp == NULL)
468                  return RC_MEMERR;
469 <        dp->ib_priv = (void *)&abase_list[inbi];
470 <        dp->ob_priv = (void *)&abase_list[outbi];
469 >        dp->ib_priv = &abase_list[inbi];
470 >        dp->ob_priv = &abase_list[outbi];
471          if (df == sd->tf) {
472 <                dp->ib_vec = ab_getvecR;
473 <                dp->ib_ndx = ab_getndxR;
474 <                dp->ob_vec = ab_getvec;
475 <                dp->ob_ndx = ab_getndx;
472 >                dp->ib_vec = &fi_getvec;
473 >                dp->ib_ndx = &fi_getndx;
474 >                dp->ob_vec = &bo_getvec;
475 >                dp->ob_ndx = &bo_getndx;
476 >        } else if (df == sd->tb) {
477 >                dp->ib_vec = &bi_getvec;
478 >                dp->ib_ndx = &bi_getndx;
479 >                dp->ob_vec = &fo_getvec;
480 >                dp->ob_ndx = &fo_getndx;
481          } else if (df == sd->rf) {
482 <                dp->ib_vec = ab_getvec;
483 <                dp->ib_ndx = ab_getndx;
484 <                dp->ob_vec = ab_getvec;
485 <                dp->ob_ndx = ab_getndx;
482 >                dp->ib_vec = &fi_getvec;
483 >                dp->ib_ndx = &fi_getndx;
484 >                dp->ob_vec = &fo_getvec;
485 >                dp->ob_ndx = &fo_getndx;
486          } else /* df == sd->rb */ {
487 <                dp->ib_vec = ab_getvecR;
488 <                dp->ib_ndx = ab_getndxR;
489 <                dp->ob_vec = ab_getvecR;
490 <                dp->ob_ndx = ab_getndxR;
487 >                dp->ib_vec = &bi_getvec;
488 >                dp->ib_ndx = &bi_getndx;
489 >                dp->ob_vec = &bo_getvec;
490 >                dp->ob_ndx = &bo_getndx;
491          }
492 <        dp->ib_ohm = ab_getohm;
493 <        dp->ob_ohm = ab_getohm;
494 <        df->comp[0].cspec[0] = c_dfcolor; /* XXX monochrome for now */
495 <        df->comp[0].dist = dp;
444 <        df->comp[0].func = &SDhandleMtx;
492 >        dp->ib_ohm = &io_getohm;
493 >        dp->ob_ohm = &io_getohm;
494 >        df->comp[ct].dist = dp;
495 >        df->comp[ct].func = &SDhandleMtx;
496                                          /* read BSDF data */
497 <        sdata  = ezxml_txt(ezxml_child(wdb,"ScatteringData"));
497 >        sdata = ezxml_txt(ezxml_child(wdb, "ScatteringData"));
498          if (!sdata || !*sdata) {
499                  sprintf(SDerrorDetail, "Missing BSDF ScatteringData in '%s'",
500                                  sd->name);
501                  return RC_FORMERR;
502          }
503          for (i = 0; i < dp->ninc*dp->nout; i++) {
504 <                char  *sdnext = fskip(sdata);
504 >                char    *sdnext = fskip(sdata);
505 >                double  val;
506                  if (sdnext == NULL) {
507                          sprintf(SDerrorDetail,
508                                  "Bad/missing BSDF ScatteringData in '%s'",
509                                          sd->name);
510                          return RC_FORMERR;
511                  }
512 <                while (*sdnext && isspace(*sdnext))
512 >                while (isspace(*sdnext))
513                          sdnext++;
514                  if (*sdnext == ',') sdnext++;
515 +                if ((val = atof(sdata)) < 0)
516 +                        val = 0;        /* don't allow negative values */
517                  if (rowinc) {
518                          int     r = i/dp->nout;
519 <                        int     c = i - c*dp->nout;
520 <                        mBSDF_value(dp,r,c) = atof(sdata);
519 >                        int     c = i - r*dp->nout;
520 >                        mBSDF_value(dp,r,c) = val;
521                  } else
522 <                        dp->bsdf[i] = atof(sdata);
522 >                        dp->bsdf[i] = val;
523                  sdata = sdnext;
524          }
525 <        return get_extrema(df);
525 >        return (ct == mtx_Y) ? get_extrema(df) : RC_GOOD;
526   }
527  
528 < /* Subtract minimum (diffuse) scattering amount from BSDF */
528 > /* copy our RGB (x,y) primary chromaticities */
529 > static void
530 > copy_RGB_prims(C_COLOR cspec[])
531 > {
532 >        if (mtx_RGB_coef[1] < .001) {   /* need to initialize */
533 >                int     i = 3;
534 >                while (i--) {
535 >                        float   rgb[3];
536 >                        rgb[0] = rgb[1] = rgb[2] = .0f;
537 >                        rgb[i] = 1.f;
538 >                        mtx_RGB_coef[i] = c_fromSharpRGB(rgb, &mtx_RGB_prim[i]);
539 >                }
540 >        }
541 >        memcpy(cspec, mtx_RGB_prim, sizeof(mtx_RGB_prim));
542 > }
543 >
544 > /* encode chromaticity if XYZ -- reduce to one channel in any case */
545 > static SDSpectralDF *
546 > encode_chroma(SDSpectralDF *df)
547 > {
548 >        SDMat   *mpx, *mpy, *mpz;
549 >        int     n;
550 >
551 >        if (df == NULL || df->ncomp != 3)
552 >                return df;
553 >
554 >        mpy = (SDMat *)df->comp[mtx_Y].dist;
555 >        if (mpy == NULL) {
556 >                free(df);
557 >                return NULL;
558 >        }
559 >        mpx = (SDMat *)df->comp[mtx_X].dist;
560 >        mpz = (SDMat *)df->comp[mtx_Z].dist;
561 >        if (mpx == NULL || (mpx->ninc != mpy->ninc) | (mpx->nout != mpy->nout))
562 >                goto done;
563 >        if (mpz == NULL || (mpz->ninc != mpy->ninc) | (mpz->nout != mpy->nout))
564 >                goto done;
565 >        mpy->chroma = (C_CHROMA *)malloc(sizeof(C_CHROMA)*mpy->ninc*mpy->nout);
566 >        if (mpy->chroma == NULL)
567 >                goto done;              /* XXX punt */
568 >                                        /* encode chroma values */
569 >        for (n = mpy->ninc*mpy->nout; n--; ) {
570 >                const double    sum = mpx->bsdf[n] + mpy->bsdf[n] + mpz->bsdf[n];
571 >                C_COLOR         cxy;
572 >                if (sum > .0)
573 >                        c_cset(&cxy, mpx->bsdf[n]/sum, mpy->bsdf[n]/sum);
574 >                else
575 >                        c_cset(&cxy, 1./3., 1./3.);
576 >                mpy->chroma[n] = c_encodeChroma(&cxy);
577 >        }
578 > done:                                   /* free X & Z channels */
579 >        if (mpx != NULL) SDfreeMatrix(mpx);
580 >        if (mpz != NULL) SDfreeMatrix(mpz);
581 >        if (mpy->chroma == NULL)        /* grayscale after all? */
582 >                df->comp[0].cspec[0] = c_dfcolor;
583 >        else                            /* else copy RGB primaries */
584 >                copy_RGB_prims(df->comp[0].cspec);
585 >        df->ncomp = 1;                  /* return resized struct */
586 >        return (SDSpectralDF *)realloc(df, sizeof(SDSpectralDF));
587 > }
588 >
589 > /* subtract minimum (diffuse) scattering amount from BSDF */
590   static double
591 < subtract_min(SDMat *sm)
591 > subtract_min(C_COLOR *cs, SDMat *sm)
592   {
593 <        float   minv = sm->bsdf[0];
594 <        int     n = sm->ninc*sm->nout;
595 <        int     i;
593 >        const int       ncomp = 1 + 2*(sm->chroma != NULL);
594 >        float           min_coef[3], ymin, coef[3];
595 >        int             i, o, c;
596          
597 <        for (i = n; --i; )
598 <                if (sm->bsdf[i] < minv)
599 <                        minv = sm->bsdf[i];
600 <        for (i = n; i--; )
601 <                sm->bsdf[i] -= minv;
597 >        min_coef[0] = min_coef[1] = min_coef[2] = FHUGE;
598 >        for (i = 0; i < sm->ninc; i++)
599 >                for (o = 0; o < sm->nout; o++) {
600 >                        c = mBSDF_color(coef, sm, i, o);
601 >                        while (c--)
602 >                                if (coef[c] < min_coef[c])
603 >                                        min_coef[c] = coef[c];
604 >                }
605 >        ymin = 0;
606 >        for (c = ncomp; c--; )
607 >                ymin += min_coef[c];
608 >        if (ymin <= .01/M_PI)           /* not worth bothering about? */
609 >                return .0;
610 >        if (ncomp == 1) {               /* subtract grayscale minimum */
611 >                for (i = sm->ninc*sm->nout; i--; )
612 >                        sm->bsdf[i] -= ymin;
613 >                *cs = c_dfcolor;
614 >                return M_PI*ymin;
615 >        }
616 >                                        /* else subtract colored minimum */
617 >        for (i = 0; i < sm->ninc; i++)
618 >                for (o = 0; o < sm->nout; o++) {
619 >                        C_COLOR cxy;
620 >                        c = mBSDF_color(coef, sm, i, o);
621 >                        while (c--)
622 >                                coef[c] = (coef[c] - min_coef[c]) /
623 >                                                mtx_RGB_coef[c];
624 >                        if (c_fromSharpRGB(coef, &cxy) > 1e-5)
625 >                                sm->chroma[o*sm->ninc + i] = c_encodeChroma(&cxy);
626 >                        mBSDF_value(sm,i,o) -= ymin;
627 >                }
628 >                                        /* return colored minimum */
629 >        for (i = 3; i--; )
630 >                coef[i] = min_coef[i]/mtx_RGB_coef[i];
631 >        c_fromSharpRGB(coef, cs);
632  
633 <        return minv*M_PI;               /* be sure to include multiplier */
633 >        return M_PI*ymin;
634   }
635  
636 < /* Extract and separate diffuse portion of BSDF */
637 < static void
636 > /* Extract and separate diffuse portion of BSDF & convert color */
637 > static SDSpectralDF *
638   extract_diffuse(SDValue *dv, SDSpectralDF *df)
639   {
495 –        int     n;
640  
641 +        df = encode_chroma(df);         /* reduce XYZ to Y + chroma */
642          if (df == NULL || df->ncomp <= 0) {
643                  dv->spec = c_dfcolor;
644                  dv->cieY = .0;
645 <                return;
645 >                return df;
646          }
647 <        dv->spec = df->comp[0].cspec[0];
648 <        dv->cieY = subtract_min((SDMat *)df->comp[0].dist);
649 <                                        /* in case of multiple components */
505 <        for (n = df->ncomp; --n; ) {
506 <                double  ymin = subtract_min((SDMat *)df->comp[n].dist);
507 <                c_cmix(&dv->spec, dv->cieY, &dv->spec, ymin, &df->comp[n].cspec[0]);
508 <                dv->cieY += ymin;
509 <        }
510 <        df->maxHemi -= dv->cieY;        /* adjust minimum hemispherical */
647 >                                        /* subtract minimum value */
648 >        dv->cieY = subtract_min(&dv->spec, (SDMat *)df->comp[0].dist);
649 >        df->maxHemi -= dv->cieY;        /* adjust maximum hemispherical */
650                                          /* make sure everything is set */
651          c_ccvt(&dv->spec, C_CSXY+C_CSSPEC);
652 +        return df;
653   }
654  
655   /* Load a BSDF matrix from an open XML file */
656   SDError
657   SDloadMtx(SDData *sd, ezxml_t wtl)
658   {
659 <        ezxml_t                 wld, wdb;
660 <        int                     rowIn;
661 <        struct BSDF_data        *dp;
662 <        char                    *txt;
663 <        int                     rval;
524 <        
659 >        ezxml_t         wld, wdb;
660 >        int             rowIn;
661 >        char            *txt;
662 >        int             rval;
663 >                                        /* basic checks and data ordering */
664          txt = ezxml_txt(ezxml_child(ezxml_child(wtl,
665                          "DataDefinition"), "IncidentDataStructure"));
666          if (txt == NULL || !*txt) {
# Line 540 | Line 679 | SDloadMtx(SDData *sd, ezxml_t wtl)
679                                  sd->name);
680                  return SDEsupport;
681          }
682 <                                /* get angle basis */
683 <        rval = load_angle_basis(ezxml_child(ezxml_child(wtl,
684 <                                "DataDefinition"), "AngleBasis"));
685 <        if (rval < 0)
686 <                return convert_errcode(rval);
687 <                                /* load BSDF components */
682 >                                        /* get angle bases */
683 >        for (wld = ezxml_child(ezxml_child(wtl, "DataDefinition"), "AngleBasis");
684 >                                wld != NULL; wld = wld->next) {
685 >                rval = load_angle_basis(wld);
686 >                if (rval < 0)
687 >                        return convert_errcode(rval);
688 >        }
689 >                                        /* load BSDF components */
690          for (wld = ezxml_child(wtl, "WavelengthData");
691                                  wld != NULL; wld = wld->next) {
692 <                if (strcasecmp(ezxml_txt(ezxml_child(wld,"Wavelength")),
693 <                                "Visible"))
694 <                        continue;       /* just visible for now */
692 >                const char      *cnm = ezxml_txt(ezxml_child(wld,"Wavelength"));
693 >                int             ct = -1;
694 >                if (!strcasecmp(cnm, "Visible"))
695 >                        ct = mtx_Y;
696 >                else if (!strcasecmp(cnm, "CIE-X"))
697 >                        ct = mtx_X;
698 >                else if (!strcasecmp(cnm, "CIE-Z"))
699 >                        ct = mtx_Z;
700 >                else
701 >                        continue;
702                  for (wdb = ezxml_child(wld, "WavelengthDataBlock");
703                                          wdb != NULL; wdb = wdb->next)
704 <                        if ((rval = load_bsdf_data(sd, wdb, rowIn)) < 0)
704 >                        if ((rval = load_bsdf_data(sd, wdb, ct, rowIn)) < 0)
705                                  return convert_errcode(rval);
706          }
707 <                                /* separate diffuse components */
708 <        extract_diffuse(&sd->rLambFront, sd->rf);
709 <        extract_diffuse(&sd->rLambBack, sd->rb);
710 <        extract_diffuse(&sd->tLamb, sd->tf);
711 <                                /* return success */
707 >                                        /* separate diffuse components */
708 >        sd->rf = extract_diffuse(&sd->rLambFront, sd->rf);
709 >        sd->rb = extract_diffuse(&sd->rLambBack, sd->rb);
710 >        if (sd->tf != NULL)
711 >                sd->tf = extract_diffuse(&sd->tLamb, sd->tf);
712 >        if (sd->tb != NULL)
713 >                sd->tb = extract_diffuse(&sd->tLamb, sd->tb);
714 >                                        /* return success */
715          return SDEnone;
716   }
717  
718   /* Get Matrix BSDF value */
719   static int
720   SDgetMtxBSDF(float coef[SDmaxCh], const FVECT outVec,
721 <                                const FVECT inVec, const void *dist)
721 >                                const FVECT inVec, SDComponent *sdc)
722   {
723 <        const SDMat     *dp = (const SDMat *)dist;
723 >        const SDMat     *dp;
724          int             i_ndx, o_ndx;
725 +                                        /* check arguments */
726 +        if ((coef == NULL) | (outVec == NULL) | (inVec == NULL) | (sdc == NULL)
727 +                        || (dp = (SDMat *)sdc->dist) == NULL)
728 +                return 0;
729                                          /* get angle indices */
730          i_ndx = mBSDF_incndx(dp, inVec);
731          o_ndx = mBSDF_outndx(dp, outVec);
# Line 581 | Line 736 | SDgetMtxBSDF(float coef[SDmaxCh], const FVECT outVec,
736          }
737          if ((i_ndx < 0) | (o_ndx < 0))
738                  return 0;               /* nothing from this component */
739 <        coef[0] = mBSDF_value(dp, i_ndx, o_ndx);
740 <        return 1;                       /* XXX monochrome for now */
739 >
740 >        return mBSDF_color(coef, dp, i_ndx, o_ndx);
741   }
742  
743 < /* Query solid angle for vector */
743 > /* Query solid angle for vector(s) */
744   static SDError
745 < SDqueryMtxProjSA(double *psa, const FVECT vec, int qflags, const void *dist)
745 > SDqueryMtxProjSA(double *psa, const FVECT v1, const RREAL *v2,
746 >                                        int qflags, SDComponent *sdc)
747   {
748 <        const SDMat     *dp = (const SDMat *)dist;
749 <
750 <        if (!(qflags & SDqueryInc+SDqueryOut))
748 >        const SDMat     *dp;
749 >        double          inc_psa, out_psa;
750 >                                        /* check arguments */
751 >        if ((psa == NULL) | (v1 == NULL) | (sdc == NULL) ||
752 >                        (dp = (SDMat *)sdc->dist) == NULL)
753                  return SDEargument;
754 <        if (qflags & SDqueryInc) {
755 <                double  inc_psa = mBSDF_incohm(dp, mBSDF_incndx(dp, vec));
756 <                if (inc_psa < .0)
757 <                        return SDEinternal;
758 <                switch (qflags & SDqueryMin+SDqueryMax) {
759 <                case SDqueryMax:
760 <                        if (inc_psa > psa[0])
761 <                                psa[0] = inc_psa;
762 <                        break;
763 <                case SDqueryMin+SDqueryMax:
764 <                        if (inc_psa > psa[1])
765 <                                psa[1] = inc_psa;
766 <                        /* fall through */
609 <                case SDqueryMin:
610 <                        if (inc_psa < psa[0])
611 <                                psa[0] = inc_psa;
612 <                        break;
613 <                case 0:
754 >        if (v2 == NULL)
755 >                v2 = v1;
756 >                                        /* get projected solid angles */
757 >        out_psa = mBSDF_outohm(dp, mBSDF_outndx(dp, v1));
758 >        inc_psa = mBSDF_incohm(dp, mBSDF_incndx(dp, v2));
759 >        if ((v1 != v2) & (out_psa <= 0) & (inc_psa <= 0)) {
760 >                inc_psa = mBSDF_outohm(dp, mBSDF_outndx(dp, v2));
761 >                out_psa = mBSDF_incohm(dp, mBSDF_incndx(dp, v1));
762 >        }
763 >
764 >        switch (qflags) {               /* record based on flag settings */
765 >        case SDqueryMax:
766 >                if (inc_psa > psa[0])
767                          psa[0] = inc_psa;
768 <                        break;
769 <                }
768 >                if (out_psa > psa[0])
769 >                        psa[0] = out_psa;
770 >                break;
771 >        case SDqueryMin+SDqueryMax:
772 >                if (inc_psa > psa[1])
773 >                        psa[1] = inc_psa;
774 >                if (out_psa > psa[1])
775 >                        psa[1] = out_psa;
776 >                /* fall through */
777 >        case SDqueryVal:
778 >                if (qflags == SDqueryVal)
779 >                        psa[0] = M_PI;
780 >                /* fall through */
781 >        case SDqueryMin:
782 >                if ((inc_psa > 0) & (inc_psa < psa[0]))
783 >                        psa[0] = inc_psa;
784 >                if ((out_psa > 0) & (out_psa < psa[0]))
785 >                        psa[0] = out_psa;
786 >                break;
787          }
788 <        if (qflags & SDqueryOut) {
789 <                double  out_psa = mBSDF_outohm(dp, mBSDF_outndx(dp, vec));
620 <                if (out_psa < .0)
621 <                        return SDEinternal;
622 <                switch (qflags & SDqueryMin+SDqueryMax) {
623 <                case SDqueryMax:
624 <                        if (out_psa > psa[0])
625 <                                psa[0] = out_psa;
626 <                        break;
627 <                case SDqueryMin+SDqueryMax:
628 <                        if (out_psa > psa[1])
629 <                                psa[1] = out_psa;
630 <                        /* fall through */
631 <                case SDqueryMin:
632 <                        if (out_psa < psa[0])
633 <                                psa[0] = out_psa;
634 <                        break;
635 <                case 0:
636 <                        psa[(qflags&SDqueryInc)!=0] = out_psa;
637 <                        break;
638 <                }
639 <        }
640 <        return SDEnone;
788 >                                        /* make sure it's legal */
789 >        return (psa[0] <= 0) ? SDEinternal : SDEnone;
790   }
791  
792   /* Compute new cumulative distribution from BSDF */
# Line 675 | Line 824 | make_cdist(SDMatCDst *cd, const FVECT inVec, SDMat *dp
824   static const SDCDst *
825   SDgetMtxCDist(const FVECT inVec, SDComponent *sdc)
826   {
827 <        SDMat           *dp = (SDMat *)sdc->dist;
827 >        SDMat           *dp;
828          int             reverse;
829          SDMatCDst       myCD;
830          SDMatCDst       *cd, *cdlast;
831 <
832 <        if (dp == NULL)
831 >                                        /* check arguments */
832 >        if ((inVec == NULL) | (sdc == NULL) ||
833 >                        (dp = (SDMat *)sdc->dist) == NULL)
834                  return NULL;
835          memset(&myCD, 0, sizeof(myCD));
836          myCD.indx = mBSDF_incndx(dp, inVec);
# Line 699 | Line 849 | SDgetMtxCDist(const FVECT inVec, SDComponent *sdc)
849                  reverse = 1;
850          }
851          cdlast = NULL;                  /* check for it in cache list */
852 <        for (cd = (SDMatCDst *)sdc->cdList;
853 <                                cd != NULL; cd = (SDMatCDst *)cd->next) {
852 >        for (cd = (SDMatCDst *)sdc->cdList; cd != NULL;
853 >                                        cdlast = cd, cd = cd->next)
854                  if (cd->indx == myCD.indx && (cd->calen == myCD.calen) &
855                                          (cd->ob_priv == myCD.ob_priv) &
856                                          (cd->ob_vec == myCD.ob_vec))
857                          break;
708 –                cdlast = cd;
709 –        }
858          if (cd == NULL) {               /* need to allocate new entry */
859                  cd = (SDMatCDst *)malloc(sizeof(SDMatCDst) +
860 <                                        myCD.calen*sizeof(myCD.carr[0]));
860 >                                        sizeof(myCD.carr[0])*myCD.calen);
861                  if (cd == NULL)
862                          return NULL;
863                  *cd = myCD;             /* compute cumulative distribution */
# Line 721 | Line 869 | SDgetMtxCDist(const FVECT inVec, SDComponent *sdc)
869          }
870          if (cdlast != NULL) {           /* move entry to head of cache list */
871                  cdlast->next = cd->next;
872 <                cd->next = sdc->cdList;
872 >                cd->next = (SDMatCDst *)sdc->cdList;
873                  sdc->cdList = (SDCDst *)cd;
874          }
875          return (SDCDst *)cd;            /* ready to go */
# Line 729 | Line 877 | SDgetMtxCDist(const FVECT inVec, SDComponent *sdc)
877  
878   /* Sample cumulative distribution */
879   static SDError
880 < SDsampMtxCDist(FVECT outVec, double randX, const SDCDst *cdp)
880 > SDsampMtxCDist(FVECT ioVec, double randX, const SDCDst *cdp)
881   {
882          const unsigned  maxval = ~0;
883          const SDMatCDst *mcd = (const SDMatCDst *)cdp;
884          const unsigned  target = randX*maxval;
885          int             i, iupper, ilower;
886 +                                        /* check arguments */
887 +        if ((ioVec == NULL) | (mcd == NULL))
888 +                return SDEargument;
889                                          /* binary search to find index */
890          ilower = 0; iupper = mcd->calen;
891          while ((i = (iupper + ilower) >> 1) != ilower)
892 <                if ((long)target >= (long)mcd->carr[i])
892 >                if (target >= mcd->carr[i])
893                          ilower = i;
894                  else
895                          iupper = i;
# Line 746 | Line 897 | SDsampMtxCDist(FVECT outVec, double randX, const SDCDs
897          randX = (randX*maxval - mcd->carr[ilower]) /
898                          (double)(mcd->carr[iupper] - mcd->carr[ilower]);
899                                          /* convert index to vector */
900 <        if ((*mcd->ob_vec)(outVec, i, randX, mcd->ob_priv))
900 >        if ((*mcd->ob_vec)(ioVec, i+randX, mcd->ob_priv))
901                  return SDEnone;
902 <        strcpy(SDerrorDetail, "BSDF sampling fault");
902 >        strcpy(SDerrorDetail, "Matrix BSDF sampling fault");
903          return SDEinternal;
904   }
905  

Diff Legend

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