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.8 by greg, Thu Feb 24 20:14:26 2011 UTC vs.
Revision 3.37 by greg, Fri Jan 5 21:00:24 2018 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 (front exiting) */
132 < static int
133 < fo_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. - rx[0])*sq(cos(M_PI/180.*ab->lat[li].tmin)) +
155 >                rx[0]*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 (front exiting) */
165 < static int
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 | fo_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; }
212 <
213 < /* get projected solid angle for this angle basis index (universal) */
214 < static double
190 > /* Get projected solid angle for this angle basis index (universal) */
191 > double
192   io_getohm(int ndx, void *p)
193   {
194 +        static void     *last_p = NULL;
195          static int      last_li = -1;
196          static double   last_ohm;
197          ANGLE_BASIS     *ab = (ANGLE_BASIS *)p;
# Line 224 | Line 202 | io_getohm(int ndx, void *p)
202                  return -1.;
203          for (li = 0; ndx >= ab->lat[li].nphis; li++)
204                  ndx -= ab->lat[li].nphis;
205 <        if (li == last_li)                      /* cached latitude? */
205 >        if ((p == last_p) & (li == last_li))            /* cached latitude? */
206                  return last_ohm;
207 +        last_p = p;
208          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)));
209          theta = M_PI/180. * ab->lat[li].tmin;
210 +        theta1 = M_PI/180. * ab->lat[li+1].tmin;
211          return last_ohm = M_PI*(sq(cos(theta)) - sq(cos(theta1))) /
212                                  (double)ab->lat[li].nphis;
213   }
214  
215 < /* get vector for this angle basis index (back incident) */
216 < static int
217 < bi_getvec(FVECT v, int ndx, double randX, void *p)
215 > /* Get vector for this angle basis index (back incident) */
216 > int
217 > bi_getvec(FVECT v, double ndxr, void *p)
218   {
219 <        if (!fo_getvec(v, ndx, randX, p))
219 >        if (!fo_getvec(v, ndxr, p))
220                  return RC_FAIL;
221  
222          v[0] = -v[0];
# Line 249 | Line 226 | bi_getvec(FVECT v, int ndx, double randX, void *p)
226          return RC_GOOD;
227   }
228  
229 < /* get index corresponding to the vector (back incident) */
230 < static int
229 > /* Get index corresponding to the vector (back incident) */
230 > int
231   bi_getndx(const FVECT v, void *p)
232   {
233          FVECT  v2;
# Line 262 | Line 239 | bi_getndx(const FVECT v, void *p)
239          return fo_getndx(v2, p);
240   }
241  
242 < /* get vector for this angle basis index (back exiting) */
243 < static int
244 < bo_getvec(FVECT v, int ndx, double randX, void *p)
242 > /* Get vector for this angle basis index (back exiting) */
243 > int
244 > bo_getvec(FVECT v, double ndxr, void *p)
245   {
246 <        if (!fo_getvec(v, ndx, randX, p))
246 >        if (!fo_getvec(v, ndxr, p))
247                  return RC_FAIL;
248  
249          v[2] = -v[2];
# Line 274 | Line 251 | bo_getvec(FVECT v, int ndx, double randX, void *p)
251          return RC_GOOD;
252   }
253  
254 < /* get index corresponding to the vector (back exiting) */
255 < static int
254 > /* Get index corresponding to the vector (back exiting) */
255 > int
256   bo_getndx(const FVECT v, void *p)
257   {
258          FVECT  v2;
# Line 287 | Line 264 | bo_getndx(const FVECT v, void *p)
264          return fo_getndx(v2, p);
265   }
266  
267 < /* get vector for this angle basis index (front incident) */
268 < static int
269 < fi_getvec(FVECT v, int ndx, double randX, void *p)
267 > /* Get vector for this angle basis index (front incident) */
268 > int
269 > fi_getvec(FVECT v, double ndxr, void *p)
270   {
271 <        if (!fo_getvec(v, ndx, randX, p))
271 >        if (!fo_getvec(v, ndxr, p))
272                  return RC_FAIL;
273  
274          v[0] = -v[0];
# Line 300 | Line 277 | fi_getvec(FVECT v, int ndx, double randX, void *p)
277          return RC_GOOD;
278   }
279  
280 < /* get index corresponding to the vector (front incident) */
281 < static int
280 > /* Get index corresponding to the vector (front incident) */
281 > int
282   fi_getndx(const FVECT v, void *p)
283   {
284          FVECT  v2;
# Line 313 | Line 290 | fi_getndx(const FVECT v, void *p)
290          return fo_getndx(v2, p);
291   }
292  
293 + /* Get color or grayscale value for BSDF for the given direction pair */
294 + int
295 + mBSDF_color(float coef[], const SDMat *dp, int i, int o)
296 + {
297 +        C_COLOR cxy;
298 +
299 +        coef[0] = mBSDF_value(dp, i, o);
300 +        if (dp->chroma == NULL)
301 +                return 1;       /* grayscale */
302 +
303 +        c_decodeChroma(&cxy, mBSDF_chroma(dp,i,o));
304 +        c_toSharpRGB(&cxy, coef[0], coef);
305 +        coef[0] *= mtx_RGB_coef[0];
306 +        coef[1] *= mtx_RGB_coef[1];
307 +        coef[2] *= mtx_RGB_coef[2];
308 +        return 3;               /* RGB color */
309 + }
310 +
311   /* load custom BSDF angle basis */
312   static int
313   load_angle_basis(ezxml_t wab)
# Line 344 | Line 339 | load_angle_basis(ezxml_t wab)
339                                  ezxml_child(ezxml_child(wbb,
340                                          "ThetaBounds"), "UpperTheta")));
341                  if (!i)
342 <                        abase_list[nabases].lat[i].tmin =
348 <                                        -abase_list[nabases].lat[i+1].tmin;
342 >                        abase_list[nabases].lat[0].tmin = 0;
343                  else if (!fequal(atof(ezxml_txt(ezxml_child(ezxml_child(wbb,
344                                          "ThetaBounds"), "LowerTheta"))),
345                                  abase_list[nabases].lat[i].tmin)) {
346                          sprintf(SDerrorDetail, "Theta values disagree in '%s'",
347 <                                abname);
347 >                                        abname);
348                          return RC_DATERR;
349                  }
350                  abase_list[nabases].nangles +=
# Line 360 | Line 354 | load_angle_basis(ezxml_t wab)
354                                  (abase_list[nabases].lat[i].nphis == 1 &&
355                                  abase_list[nabases].lat[i].tmin > FTINY)) {
356                          sprintf(SDerrorDetail, "Illegal phi count in '%s'",
357 <                                abname);
357 >                                        abname);
358                          return RC_DATERR;
359                  }
360          }
# Line 406 | Line 400 | get_extrema(SDSpectralDF *df)
400  
401   /* load BSDF distribution for this wavelength */
402   static int
403 < load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc)
403 > load_bsdf_data(SDData *sd, ezxml_t wdb, int ct, int rowinc)
404   {
405          SDSpectralDF    *df;
406          SDMat           *dp;
407          char            *sdata;
414        int             tfront;
408          int             inbi, outbi;
409          int             i;
410                                          /* allocate BSDF component */
411          sdata = ezxml_txt(ezxml_child(wdb, "WavelengthDataDirection"));
412 +        if (!sdata)
413 +                return RC_FAIL;
414          /*
415           * Remember that front and back are reversed from WINDOW 6 orientations
416           */
417 <        if ((tfront = !strcasecmp(sdata, "Transmission Back")) ||
418 <                        (sd->tf == NULL &&
424 <                                !strcasecmp(sdata, "Transmission Front"))) {
425 <                if (sd->tf != NULL)
426 <                        SDfreeSpectralDF(sd->tf);
427 <                if ((sd->tf = SDnewSpectralDF(1)) == NULL)
417 >        if (!strcasecmp(sdata, "Transmission Front")) {
418 >                if (sd->tb == NULL && (sd->tb = SDnewSpectralDF(3)) == NULL)
419                          return RC_MEMERR;
420 +                df = sd->tb;
421 +        } else if (!strcasecmp(sdata, "Transmission Back")) {
422 +                if (sd->tf == NULL && (sd->tf = SDnewSpectralDF(3)) == NULL)
423 +                        return RC_MEMERR;
424                  df = sd->tf;
425          } else if (!strcasecmp(sdata, "Reflection Front")) {
426 <                if (sd->rb != NULL)     /* note back-front reversal */
432 <                        SDfreeSpectralDF(sd->rb);
433 <                if ((sd->rb = SDnewSpectralDF(1)) == NULL)
426 >                if (sd->rb == NULL && (sd->rb = SDnewSpectralDF(3)) == NULL)
427                          return RC_MEMERR;
428                  df = sd->rb;
429          } else if (!strcasecmp(sdata, "Reflection Back")) {
430 <                if (sd->rf != NULL)     /* note front-back reversal */
438 <                        SDfreeSpectralDF(sd->rf);
439 <                if ((sd->rf = SDnewSpectralDF(1)) == NULL)
430 >                if (sd->rf == NULL && (sd->rf = SDnewSpectralDF(3)) == NULL)
431                          return RC_MEMERR;
432                  df = sd->rf;
433          } else
434                  return RC_FAIL;
435 <        /* XXX should also check "ScatteringDataType" for consistency? */
435 >                                        /* free previous matrix if any */
436 >        if (df->comp[ct].dist != NULL) {
437 >                SDfreeMatrix(df->comp[ct].dist);
438 >                df->comp[ct].dist = NULL;
439 >        }
440                                          /* get angle bases */
441          sdata = ezxml_txt(ezxml_child(wdb,"ColumnAngleBasis"));
442          if (!sdata || !*sdata) {
# Line 476 | Line 471 | load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc)
471          dp->ib_priv = &abase_list[inbi];
472          dp->ob_priv = &abase_list[outbi];
473          if (df == sd->tf) {
474 <                if (tfront) {
475 <                        dp->ib_vec = &fi_getvec;
476 <                        dp->ib_ndx = &fi_getndx;
477 <                        dp->ob_vec = &bo_getvec;
478 <                        dp->ob_ndx = &bo_getndx;
479 <                } else {
480 <                        dp->ib_vec = &bi_getvec;
481 <                        dp->ib_ndx = &bi_getndx;
482 <                        dp->ob_vec = &fo_getvec;
488 <                        dp->ob_ndx = &fo_getndx;
489 <                }
474 >                dp->ib_vec = &fi_getvec;
475 >                dp->ib_ndx = &fi_getndx;
476 >                dp->ob_vec = &bo_getvec;
477 >                dp->ob_ndx = &bo_getndx;
478 >        } else if (df == sd->tb) {
479 >                dp->ib_vec = &bi_getvec;
480 >                dp->ib_ndx = &bi_getndx;
481 >                dp->ob_vec = &fo_getvec;
482 >                dp->ob_ndx = &fo_getndx;
483          } else if (df == sd->rf) {
484                  dp->ib_vec = &fi_getvec;
485                  dp->ib_ndx = &fi_getndx;
# Line 500 | Line 493 | load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc)
493          }
494          dp->ib_ohm = &io_getohm;
495          dp->ob_ohm = &io_getohm;
496 <        df->comp[0].cspec[0] = c_dfcolor; /* XXX monochrome for now */
497 <        df->comp[0].dist = dp;
505 <        df->comp[0].func = &SDhandleMtx;
496 >        df->comp[ct].dist = dp;
497 >        df->comp[ct].func = &SDhandleMtx;
498                                          /* read BSDF data */
499 <        sdata  = ezxml_txt(ezxml_child(wdb,"ScatteringData"));
499 >        sdata = ezxml_txt(ezxml_child(wdb, "ScatteringData"));
500          if (!sdata || !*sdata) {
501                  sprintf(SDerrorDetail, "Missing BSDF ScatteringData in '%s'",
502                                  sd->name);
503                  return RC_FORMERR;
504          }
505          for (i = 0; i < dp->ninc*dp->nout; i++) {
506 <                char  *sdnext = fskip(sdata);
506 >                char    *sdnext = fskip(sdata);
507 >                double  val;
508                  if (sdnext == NULL) {
509                          sprintf(SDerrorDetail,
510                                  "Bad/missing BSDF ScatteringData in '%s'",
511                                          sd->name);
512                          return RC_FORMERR;
513                  }
514 <                while (*sdnext && isspace(*sdnext))
514 >                while (isspace(*sdnext))
515                          sdnext++;
516                  if (*sdnext == ',') sdnext++;
517 +                if ((val = atof(sdata)) < 0)
518 +                        val = 0;        /* don't allow negative values */
519                  if (rowinc) {
520                          int     r = i/dp->nout;
521 <                        int     c = i - c*dp->nout;
522 <                        mBSDF_value(dp,r,c) = atof(sdata);
521 >                        int     c = i - r*dp->nout;
522 >                        mBSDF_value(dp,r,c) = val;
523                  } else
524 <                        dp->bsdf[i] = atof(sdata);
524 >                        dp->bsdf[i] = val;
525                  sdata = sdnext;
526          }
527 <        return get_extrema(df);
527 >        return (ct == mtx_Y) ? get_extrema(df) : RC_GOOD;
528   }
529  
530 < /* Subtract minimum (diffuse) scattering amount from BSDF */
530 > /* copy our RGB (x,y) primary chromaticities */
531 > static void
532 > copy_RGB_prims(C_COLOR cspec[])
533 > {
534 >        if (mtx_RGB_coef[1] < .001) {   /* need to initialize */
535 >                int     i = 3;
536 >                while (i--) {
537 >                        float   rgb[3];
538 >                        rgb[0] = rgb[1] = rgb[2] = .0f;
539 >                        rgb[i] = 1.f;
540 >                        mtx_RGB_coef[i] = c_fromSharpRGB(rgb, &mtx_RGB_prim[i]);
541 >                }
542 >        }
543 >        memcpy(cspec, mtx_RGB_prim, sizeof(mtx_RGB_prim));
544 > }
545 >
546 > /* encode chromaticity if XYZ -- reduce to one channel in any case */
547 > static SDSpectralDF *
548 > encode_chroma(SDSpectralDF *df)
549 > {
550 >        SDMat   *mpx, *mpy, *mpz;
551 >        int     n;
552 >
553 >        if (df == NULL || df->ncomp != 3)
554 >                return df;
555 >
556 >        mpy = (SDMat *)df->comp[mtx_Y].dist;
557 >        if (mpy == NULL) {
558 >                free(df);
559 >                return NULL;
560 >        }
561 >        mpx = (SDMat *)df->comp[mtx_X].dist;
562 >        mpz = (SDMat *)df->comp[mtx_Z].dist;
563 >        if (mpx == NULL || (mpx->ninc != mpy->ninc) | (mpx->nout != mpy->nout))
564 >                goto done;
565 >        if (mpz == NULL || (mpz->ninc != mpy->ninc) | (mpz->nout != mpy->nout))
566 >                goto done;
567 >        mpy->chroma = (C_CHROMA *)malloc(sizeof(C_CHROMA)*mpy->ninc*mpy->nout);
568 >        if (mpy->chroma == NULL)
569 >                goto done;              /* XXX punt */
570 >                                        /* encode chroma values */
571 >        for (n = mpy->ninc*mpy->nout; n--; ) {
572 >                const double    sum = mpx->bsdf[n] + mpy->bsdf[n] + mpz->bsdf[n];
573 >                C_COLOR         cxy;
574 >                if (sum > .0)
575 >                        c_cset(&cxy, mpx->bsdf[n]/sum, mpy->bsdf[n]/sum);
576 >                else
577 >                        c_cset(&cxy, 1./3., 1./3.);
578 >                mpy->chroma[n] = c_encodeChroma(&cxy);
579 >        }
580 > done:                                   /* free X & Z channels */
581 >        if (mpx != NULL) SDfreeMatrix(mpx);
582 >        if (mpz != NULL) SDfreeMatrix(mpz);
583 >        if (mpy->chroma == NULL)        /* grayscale after all? */
584 >                df->comp[0].cspec[0] = c_dfcolor;
585 >        else                            /* else copy RGB primaries */
586 >                copy_RGB_prims(df->comp[0].cspec);
587 >        df->ncomp = 1;                  /* return resized struct */
588 >        return (SDSpectralDF *)realloc(df, sizeof(SDSpectralDF));
589 > }
590 >
591 > /* subtract minimum (diffuse) scattering amount from BSDF */
592   static double
593 < subtract_min(SDMat *sm)
593 > subtract_min(C_COLOR *cs, SDMat *sm)
594   {
595 <        float   minv = sm->bsdf[0];
596 <        int     n = sm->ninc*sm->nout;
597 <        int     i;
595 >        const int       ncomp = 1 + 2*(sm->chroma != NULL);
596 >        float           min_coef[3], ymin, coef[3];
597 >        int             i, o, c;
598          
599 <        for (i = n; --i; )
600 <                if (sm->bsdf[i] < minv)
601 <                        minv = sm->bsdf[i];
602 <        for (i = n; i--; )
603 <                sm->bsdf[i] -= minv;
599 >        min_coef[0] = min_coef[1] = min_coef[2] = FHUGE;
600 >        for (i = 0; i < sm->ninc; i++)
601 >                for (o = 0; o < sm->nout; o++) {
602 >                        c = mBSDF_color(coef, sm, i, o);
603 >                        while (c--)
604 >                                if (coef[c] < min_coef[c])
605 >                                        min_coef[c] = coef[c];
606 >                }
607 >        ymin = 0;
608 >        for (c = ncomp; c--; )
609 >                ymin += min_coef[c];
610 >        if (ymin <= .01/M_PI)           /* not worth bothering about? */
611 >                return .0;
612 >        if (ncomp == 1) {               /* subtract grayscale minimum */
613 >                for (i = sm->ninc*sm->nout; i--; )
614 >                        sm->bsdf[i] -= ymin;
615 >                *cs = c_dfcolor;
616 >                return M_PI*ymin;
617 >        }
618 >                                        /* else subtract colored minimum */
619 >        for (i = 0; i < sm->ninc; i++)
620 >                for (o = 0; o < sm->nout; o++) {
621 >                        C_COLOR cxy;
622 >                        c = mBSDF_color(coef, sm, i, o);
623 >                        while (c--)
624 >                                coef[c] = (coef[c] - min_coef[c]) /
625 >                                                mtx_RGB_coef[c];
626 >                        if (c_fromSharpRGB(coef, &cxy) > 1e-5)
627 >                                mBSDF_chroma(sm,i,o) = c_encodeChroma(&cxy);
628 >                        mBSDF_value(sm,i,o) -= ymin;
629 >                }
630 >                                        /* return colored minimum */
631 >        for (i = 3; i--; )
632 >                coef[i] = min_coef[i]/mtx_RGB_coef[i];
633 >        c_fromSharpRGB(coef, cs);
634  
635 <        return minv*M_PI;               /* be sure to include multiplier */
635 >        return M_PI*ymin;
636   }
637  
638 < /* Extract and separate diffuse portion of BSDF */
639 < static void
638 > /* Extract and separate diffuse portion of BSDF & convert color */
639 > static SDSpectralDF *
640   extract_diffuse(SDValue *dv, SDSpectralDF *df)
641   {
556        int     n;
642  
643 +        df = encode_chroma(df);         /* reduce XYZ to Y + chroma */
644          if (df == NULL || df->ncomp <= 0) {
645                  dv->spec = c_dfcolor;
646                  dv->cieY = .0;
647 <                return;
647 >                return df;
648          }
649 <        dv->spec = df->comp[0].cspec[0];
650 <        dv->cieY = subtract_min((SDMat *)df->comp[0].dist);
651 <                                        /* in case of multiple components */
566 <        for (n = df->ncomp; --n; ) {
567 <                double  ymin = subtract_min((SDMat *)df->comp[n].dist);
568 <                c_cmix(&dv->spec, dv->cieY, &dv->spec, ymin, &df->comp[n].cspec[0]);
569 <                dv->cieY += ymin;
570 <        }
571 <        df->maxHemi -= dv->cieY;        /* adjust minimum hemispherical */
649 >                                        /* subtract minimum value */
650 >        dv->cieY = subtract_min(&dv->spec, (SDMat *)df->comp[0].dist);
651 >        df->maxHemi -= dv->cieY;        /* adjust maximum hemispherical */
652                                          /* make sure everything is set */
653          c_ccvt(&dv->spec, C_CSXY+C_CSSPEC);
654 +        return df;
655   }
656  
657   /* Load a BSDF matrix from an open XML file */
658   SDError
659   SDloadMtx(SDData *sd, ezxml_t wtl)
660   {
661 <        ezxml_t                 wld, wdb;
662 <        int                     rowIn;
663 <        struct BSDF_data        *dp;
664 <        char                    *txt;
665 <        int                     rval;
585 <        
661 >        ezxml_t         wld, wdb;
662 >        int             rowIn;
663 >        char            *txt;
664 >        int             rval;
665 >                                        /* basic checks and data ordering */
666          txt = ezxml_txt(ezxml_child(ezxml_child(wtl,
667                          "DataDefinition"), "IncidentDataStructure"));
668          if (txt == NULL || !*txt) {
# Line 601 | Line 681 | SDloadMtx(SDData *sd, ezxml_t wtl)
681                                  sd->name);
682                  return SDEsupport;
683          }
684 <                                /* get angle basis */
685 <        rval = load_angle_basis(ezxml_child(ezxml_child(wtl,
686 <                                "DataDefinition"), "AngleBasis"));
687 <        if (rval < 0)
688 <                return convert_errcode(rval);
689 <                                /* load BSDF components */
684 >                                        /* get angle bases */
685 >        for (wld = ezxml_child(ezxml_child(wtl, "DataDefinition"), "AngleBasis");
686 >                                wld != NULL; wld = wld->next) {
687 >                rval = load_angle_basis(wld);
688 >                if (rval < 0)
689 >                        return convert_errcode(rval);
690 >        }
691 >                                        /* load BSDF components */
692          for (wld = ezxml_child(wtl, "WavelengthData");
693                                  wld != NULL; wld = wld->next) {
694 <                if (strcasecmp(ezxml_txt(ezxml_child(wld,"Wavelength")),
695 <                                "Visible"))
696 <                        continue;       /* just visible for now */
694 >                const char      *cnm = ezxml_txt(ezxml_child(wld,"Wavelength"));
695 >                int             ct = -1;
696 >                if (!strcasecmp(cnm, "Visible"))
697 >                        ct = mtx_Y;
698 >                else if (!strcasecmp(cnm, "CIE-X"))
699 >                        ct = mtx_X;
700 >                else if (!strcasecmp(cnm, "CIE-Z"))
701 >                        ct = mtx_Z;
702 >                else
703 >                        continue;
704                  for (wdb = ezxml_child(wld, "WavelengthDataBlock");
705                                          wdb != NULL; wdb = wdb->next)
706 <                        if ((rval = load_bsdf_data(sd, wdb, rowIn)) < 0)
706 >                        if ((rval = load_bsdf_data(sd, wdb, ct, rowIn)) < 0)
707                                  return convert_errcode(rval);
708          }
709 <                                /* separate diffuse components */
710 <        extract_diffuse(&sd->rLambFront, sd->rf);
711 <        extract_diffuse(&sd->rLambBack, sd->rb);
712 <        extract_diffuse(&sd->tLamb, sd->tf);
713 <                                /* return success */
709 >                                        /* separate diffuse components */
710 >        sd->rf = extract_diffuse(&sd->rLambFront, sd->rf);
711 >        sd->rb = extract_diffuse(&sd->rLambBack, sd->rb);
712 >        if (sd->tf != NULL)
713 >                sd->tf = extract_diffuse(&sd->tLamb, sd->tf);
714 >        if (sd->tb != NULL)
715 >                sd->tb = extract_diffuse(&sd->tLamb, sd->tb);
716 >                                        /* return success */
717          return SDEnone;
718   }
719  
720   /* Get Matrix BSDF value */
721   static int
722   SDgetMtxBSDF(float coef[SDmaxCh], const FVECT outVec,
723 <                                const FVECT inVec, const void *dist)
723 >                                const FVECT inVec, SDComponent *sdc)
724   {
725 <        const SDMat     *dp = (const SDMat *)dist;
725 >        const SDMat     *dp;
726          int             i_ndx, o_ndx;
727 +                                        /* check arguments */
728 +        if ((coef == NULL) | (outVec == NULL) | (inVec == NULL) | (sdc == NULL)
729 +                        || (dp = (SDMat *)sdc->dist) == NULL)
730 +                return 0;
731                                          /* get angle indices */
732          i_ndx = mBSDF_incndx(dp, inVec);
733          o_ndx = mBSDF_outndx(dp, outVec);
# Line 642 | Line 738 | SDgetMtxBSDF(float coef[SDmaxCh], const FVECT outVec,
738          }
739          if ((i_ndx < 0) | (o_ndx < 0))
740                  return 0;               /* nothing from this component */
741 <        coef[0] = mBSDF_value(dp, i_ndx, o_ndx);
742 <        return 1;                       /* XXX monochrome for now */
741 >
742 >        return mBSDF_color(coef, dp, i_ndx, o_ndx);
743   }
744  
745 < /* Query solid angle for vector */
745 > /* Query solid angle for vector(s) */
746   static SDError
747 < SDqueryMtxProjSA(double *psa, const FVECT vec, int qflags, const void *dist)
747 > SDqueryMtxProjSA(double *psa, const FVECT v1, const RREAL *v2,
748 >                                        int qflags, SDComponent *sdc)
749   {
750 <        const SDMat     *dp = (const SDMat *)dist;
750 >        const SDMat     *dp;
751          double          inc_psa, out_psa;
752                                          /* check arguments */
753 <        if ((psa == NULL) | (vec == NULL) | (dp == NULL))
753 >        if ((psa == NULL) | (v1 == NULL) | (sdc == NULL) ||
754 >                        (dp = (SDMat *)sdc->dist) == NULL)
755                  return SDEargument;
756 +        if (v2 == NULL)
757 +                v2 = v1;
758                                          /* get projected solid angles */
759 <        inc_psa = mBSDF_incohm(dp, mBSDF_incndx(dp, vec));
760 <        out_psa = mBSDF_outohm(dp, mBSDF_outndx(dp, vec));
759 >        out_psa = mBSDF_outohm(dp, mBSDF_outndx(dp, v1));
760 >        inc_psa = mBSDF_incohm(dp, mBSDF_incndx(dp, v2));
761 >        if ((v1 != v2) & (out_psa <= 0) & (inc_psa <= 0)) {
762 >                inc_psa = mBSDF_outohm(dp, mBSDF_outndx(dp, v2));
763 >                out_psa = mBSDF_incohm(dp, mBSDF_incndx(dp, v1));
764 >        }
765  
766          switch (qflags) {               /* record based on flag settings */
663        case SDqueryVal:
664                psa[0] = .0;
665                /* fall through */
767          case SDqueryMax:
768                  if (inc_psa > psa[0])
769                          psa[0] = inc_psa;
# Line 670 | Line 771 | SDqueryMtxProjSA(double *psa, const FVECT vec, int qfl
771                          psa[0] = out_psa;
772                  break;
773          case SDqueryMin+SDqueryMax:
774 <                if (inc_psa > psa[0])
774 >                if (inc_psa > psa[1])
775                          psa[1] = inc_psa;
776 <                if (out_psa > psa[0])
776 >                if (out_psa > psa[1])
777                          psa[1] = out_psa;
778                  /* fall through */
779 +        case SDqueryVal:
780 +                if (qflags == SDqueryVal)
781 +                        psa[0] = M_PI;
782 +                /* fall through */
783          case SDqueryMin:
784 <                if ((inc_psa > .0) & (inc_psa < psa[0]))
784 >                if ((inc_psa > 0) & (inc_psa < psa[0]))
785                          psa[0] = inc_psa;
786 <                if ((out_psa > .0) & (out_psa < psa[0]))
786 >                if ((out_psa > 0) & (out_psa < psa[0]))
787                          psa[0] = out_psa;
788                  break;
789          }
790                                          /* make sure it's legal */
791 <        return (psa[0] <= .0) ? SDEinternal : SDEnone;
791 >        return (psa[0] <= 0) ? SDEinternal : SDEnone;
792   }
793  
794   /* Compute new cumulative distribution from BSDF */
# Line 721 | Line 826 | make_cdist(SDMatCDst *cd, const FVECT inVec, SDMat *dp
826   static const SDCDst *
827   SDgetMtxCDist(const FVECT inVec, SDComponent *sdc)
828   {
829 <        SDMat           *dp = (SDMat *)sdc->dist;
829 >        SDMat           *dp;
830          int             reverse;
831          SDMatCDst       myCD;
832          SDMatCDst       *cd, *cdlast;
833                                          /* check arguments */
834 <        if ((inVec == NULL) | (dp == NULL))
834 >        if ((inVec == NULL) | (sdc == NULL) ||
835 >                        (dp = (SDMat *)sdc->dist) == NULL)
836                  return NULL;
837          memset(&myCD, 0, sizeof(myCD));
838          myCD.indx = mBSDF_incndx(dp, inVec);
# Line 745 | Line 851 | SDgetMtxCDist(const FVECT inVec, SDComponent *sdc)
851                  reverse = 1;
852          }
853          cdlast = NULL;                  /* check for it in cache list */
854 <        for (cd = (SDMatCDst *)sdc->cdList;
855 <                                cd != NULL; cd = (SDMatCDst *)cd->next) {
854 >        /* PLACE MUTEX LOCK HERE FOR THREAD-SAFE */
855 >        for (cd = (SDMatCDst *)sdc->cdList; cd != NULL;
856 >                                        cdlast = cd, cd = cd->next)
857                  if (cd->indx == myCD.indx && (cd->calen == myCD.calen) &
858                                          (cd->ob_priv == myCD.ob_priv) &
859                                          (cd->ob_vec == myCD.ob_vec))
860                          break;
754                cdlast = cd;
755        }
861          if (cd == NULL) {               /* need to allocate new entry */
862                  cd = (SDMatCDst *)malloc(sizeof(SDMatCDst) +
863 <                                        myCD.calen*sizeof(myCD.carr[0]));
863 >                                        sizeof(myCD.carr[0])*myCD.calen);
864                  if (cd == NULL)
865                          return NULL;
866                  *cd = myCD;             /* compute cumulative distribution */
# Line 767 | Line 872 | SDgetMtxCDist(const FVECT inVec, SDComponent *sdc)
872          }
873          if (cdlast != NULL) {           /* move entry to head of cache list */
874                  cdlast->next = cd->next;
875 <                cd->next = sdc->cdList;
875 >                cd->next = (SDMatCDst *)sdc->cdList;
876                  sdc->cdList = (SDCDst *)cd;
877          }
878 +        /* END MUTEX LOCK */
879          return (SDCDst *)cd;            /* ready to go */
880   }
881  
882   /* Sample cumulative distribution */
883   static SDError
884 < SDsampMtxCDist(FVECT outVec, double randX, const SDCDst *cdp)
884 > SDsampMtxCDist(FVECT ioVec, double randX, const SDCDst *cdp)
885   {
886          const unsigned  maxval = ~0;
887          const SDMatCDst *mcd = (const SDMatCDst *)cdp;
888          const unsigned  target = randX*maxval;
889          int             i, iupper, ilower;
890                                          /* check arguments */
891 <        if ((outVec == NULL) | (mcd == NULL))
891 >        if ((ioVec == NULL) | (mcd == NULL))
892                  return SDEargument;
893                                          /* binary search to find index */
894          ilower = 0; iupper = mcd->calen;
895          while ((i = (iupper + ilower) >> 1) != ilower)
896 <                if ((long)target >= (long)mcd->carr[i])
896 >                if (target >= mcd->carr[i])
897                          ilower = i;
898                  else
899                          iupper = i;
# Line 795 | Line 901 | SDsampMtxCDist(FVECT outVec, double randX, const SDCDs
901          randX = (randX*maxval - mcd->carr[ilower]) /
902                          (double)(mcd->carr[iupper] - mcd->carr[ilower]);
903                                          /* convert index to vector */
904 <        if ((*mcd->ob_vec)(outVec, i, randX, mcd->ob_priv))
904 >        if ((*mcd->ob_vec)(ioVec, i+randX, mcd->ob_priv))
905                  return SDEnone;
906 <        strcpy(SDerrorDetail, "BSDF sampling fault");
906 >        strcpy(SDerrorDetail, "Matrix BSDF sampling fault");
907          return SDEinternal;
908   }
909  
910   /* Fixed resolution BSDF methods */
911 < SDFunc                  SDhandleMtx = {
911 > const SDFunc            SDhandleMtx = {
912                                  &SDgetMtxBSDF,
913                                  &SDqueryMtxProjSA,
914                                  &SDgetMtxCDist,

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines