--- ray/src/common/bsdf_m.c 2011/02/18 00:40:25 3.1 +++ ray/src/common/bsdf_m.c 2015/04/08 02:41:02 3.32 @@ -1,3 +1,6 @@ +#ifndef lint +static const char RCSid[] = "$Id: bsdf_m.c,v 3.32 2015/04/08 02:41:02 greg Exp $"; +#endif /* * bsdf_m.c * @@ -8,19 +11,14 @@ * */ -#include -#include +#define _USE_MATH_DEFINES +#include "rtio.h" #include -#include #include #include "ezxml.h" #include "bsdf.h" #include "bsdf_m.h" -#ifndef FTINY -#define FTINY 1e-6 -#endif - /* Function return codes */ #define RC_GOOD 1 #define RC_FAIL 0 @@ -30,24 +28,10 @@ #define RC_INTERR (-4) #define RC_MEMERR (-5) -#define MAXLATS 46 /* maximum number of latitudes */ - -/* BSDF angle specification */ -typedef struct { - char name[64]; /* basis name */ - int nangles; /* total number of directions */ - struct { - float tmin; /* starting theta */ - int nphis; /* number of phis (0 term) */ - } lat[MAXLATS+1]; /* latitudes */ -} ANGLE_BASIS; - -#define MAXABASES 7 /* limit on defined bases */ - -static ANGLE_BASIS abase_list[MAXABASES] = { +ANGLE_BASIS abase_list[MAXABASES] = { { "LBNL/Klems Full", 145, - { {-5., 1}, + { {0., 1}, {5., 8}, {15., 16}, {25., 20}, @@ -59,7 +43,7 @@ static ANGLE_BASIS abase_list[MAXABASES] = { {90., 0} } }, { "LBNL/Klems Half", 73, - { {-6.5, 1}, + { {0., 1}, {6.5, 8}, {19.5, 12}, {32.5, 16}, @@ -69,7 +53,7 @@ static ANGLE_BASIS abase_list[MAXABASES] = { {90., 0} } }, { "LBNL/Klems Quarter", 41, - { {-9., 1}, + { {0., 1}, {9., 8}, {27., 12}, {46., 12}, @@ -78,41 +62,23 @@ static ANGLE_BASIS abase_list[MAXABASES] = { } }; -static int nabases = 3; /* current number of defined bases */ +int nabases = 3; /* current number of defined bases */ +C_COLOR mtx_RGB_prim[3]; /* our RGB primaries */ +float mtx_RGB_coef[3]; /* corresponding Y coefficients */ + +enum {mtx_Y, mtx_X, mtx_Z}; /* matrix components (mtx_Y==0) */ + +/* check if two real values are near enough to equal */ static int fequal(double a, double b) { - if (b != .0) + if (b != 0) a = a/b - 1.; return (a <= 1e-6) & (a >= -1e-6); } -/* returns the name of the given tag */ -#ifdef ezxml_name -#undef ezxml_name -static char * -ezxml_name(ezxml_t xml) -{ - if (xml == NULL) - return NULL; - return xml->name; -} -#endif - -/* returns the given tag's character content or empty string if none */ -#ifdef ezxml_txt -#undef ezxml_txt -static char * -ezxml_txt(ezxml_t xml) -{ - if (xml == NULL) - return ""; - return xml->txt; -} -#endif - -/* Convert error to standard BSDF code */ +/* convert error to standard BSDF code */ static SDError convert_errcode(int ec) { @@ -133,7 +99,7 @@ convert_errcode(int ec) return SDEunknown; } -/* Allocate a BSDF matrix of the given size */ +/* allocate a BSDF matrix of the given size */ static SDMat * SDnewMatrix(int ni, int no) { @@ -157,45 +123,57 @@ SDnewMatrix(int ni, int no) } /* Free a BSDF matrix */ -#define SDfreeMatrix free +void +SDfreeMatrix(void *ptr) +{ + SDMat *mp = (SDMat *)ptr; -/* get vector for this angle basis index */ -static int -ab_getvec(FVECT v, int ndx, double randX, void *p) + if (mp->chroma != NULL) free(mp->chroma); + free(ptr); +} + +/* compute square of real value */ +static double sq(double x) { return x*x; } + +/* Get vector for this angle basis index (front exiting) */ +int +fo_getvec(FVECT v, double ndxr, void *p) { - ANGLE_BASIS *ab = (ANGLE_BASIS *)p; + ANGLE_BASIS *ab = (ANGLE_BASIS *)p; + int ndx = (int)ndxr; + double randX = ndxr - ndx; double rx[2]; int li; - double pol, azi, d; + double azi, d; - if ((ndx < 0) | (ndx >= ab->nangles)) + if ((ndxr < 0) | (ndx >= ab->nangles)) return RC_FAIL; for (li = 0; ndx >= ab->lat[li].nphis; li++) ndx -= ab->lat[li].nphis; SDmultiSamp(rx, 2, randX); - pol = M_PI/180.*( (1.-rx[0])*ab->lat[li].tmin + - rx[0]*ab->lat[li+1].tmin ); + d = (1. - randX)*sq(cos(M_PI/180.*ab->lat[li].tmin)) + + randX*sq(cos(M_PI/180.*ab->lat[li+1].tmin)); + v[2] = d = sqrt(d); /* cos(pol) */ azi = 2.*M_PI*(ndx + rx[1] - .5)/ab->lat[li].nphis; - v[2] = d = cos(pol); d = sqrt(1. - d*d); /* sin(pol) */ v[0] = cos(azi)*d; v[1] = sin(azi)*d; return RC_GOOD; } -/* get index corresponding to the given vector */ -static int -ab_getndx(const FVECT v, void *p) +/* Get index corresponding to the given vector (front exiting) */ +int +fo_getndx(const FVECT v, void *p) { - ANGLE_BASIS *ab = (ANGLE_BASIS *)p; + ANGLE_BASIS *ab = (ANGLE_BASIS *)p; int li, ndx; - double pol, azi, d; + double pol, azi; if (v == NULL) return -1; - if ((v[2] < .0) | (v[2] > 1.0)) + if ((v[2] < 0) | (v[2] > 1.)) return -1; - pol = 180.0/M_PI*acos(v[2]); + pol = 180.0/M_PI*Acos(v[2]); azi = 180.0/M_PI*atan2(v[1], v[0]); if (azi < 0.0) azi += 360.0; for (li = 1; ab->lat[li].tmin <= pol; li++) @@ -209,12 +187,9 @@ ab_getndx(const FVECT v, void *p) return ndx; } -/* compute square of real value */ -static double sq(double x) { return x*x; } - -/* get projected solid angle for this angle basis index */ -static double -ab_getohm(int ndx, void *p) +/* Get projected solid angle for this angle basis index (universal) */ +double +io_getohm(int ndx, void *p) { static int last_li = -1; static double last_ohm; @@ -229,21 +204,17 @@ ab_getohm(int ndx, void *p) if (li == last_li) /* cached latitude? */ return last_ohm; last_li = li; - theta1 = M_PI/180. * ab->lat[li+1].tmin; - if (ab->lat[li].nphis == 1) /* special case */ - return last_ohm = M_PI*(1. - sq(cos(theta1))); theta = M_PI/180. * ab->lat[li].tmin; + theta1 = M_PI/180. * ab->lat[li+1].tmin; return last_ohm = M_PI*(sq(cos(theta)) - sq(cos(theta1))) / (double)ab->lat[li].nphis; } -/* get reverse vector for this angle basis index */ -static int -ab_getvecR(FVECT v, int ndx, double randX, void *p) +/* Get vector for this angle basis index (back incident) */ +int +bi_getvec(FVECT v, double ndxr, void *p) { - int na = (*(ANGLE_BASIS *)p).nangles; - - if (!ab_getvec(v, ndx, randX, p)) + if (!fo_getvec(v, ndxr, p)) return RC_FAIL; v[0] = -v[0]; @@ -253,9 +224,9 @@ ab_getvecR(FVECT v, int ndx, double randX, void *p) return RC_GOOD; } -/* get index corresponding to the reverse vector */ -static int -ab_getndxR(const FVECT v, void *p) +/* Get index corresponding to the vector (back incident) */ +int +bi_getndx(const FVECT v, void *p) { FVECT v2; @@ -263,9 +234,78 @@ ab_getndxR(const FVECT v, void *p) v2[1] = -v[1]; v2[2] = -v[2]; - return ab_getndx(v2, p); + return fo_getndx(v2, p); } +/* Get vector for this angle basis index (back exiting) */ +int +bo_getvec(FVECT v, double ndxr, void *p) +{ + if (!fo_getvec(v, ndxr, p)) + return RC_FAIL; + + v[2] = -v[2]; + + return RC_GOOD; +} + +/* Get index corresponding to the vector (back exiting) */ +int +bo_getndx(const FVECT v, void *p) +{ + FVECT v2; + + v2[0] = v[0]; + v2[1] = v[1]; + v2[2] = -v[2]; + + return fo_getndx(v2, p); +} + +/* Get vector for this angle basis index (front incident) */ +int +fi_getvec(FVECT v, double ndxr, void *p) +{ + if (!fo_getvec(v, ndxr, p)) + return RC_FAIL; + + v[0] = -v[0]; + v[1] = -v[1]; + + return RC_GOOD; +} + +/* Get index corresponding to the vector (front incident) */ +int +fi_getndx(const FVECT v, void *p) +{ + FVECT v2; + + v2[0] = -v[0]; + v2[1] = -v[1]; + v2[2] = v[2]; + + return fo_getndx(v2, p); +} + +/* Get color or grayscale value for BSDF for the given direction pair */ +int +mBSDF_color(float coef[], const SDMat *dp, int i, int o) +{ + C_COLOR cxy; + + coef[0] = mBSDF_value(dp, i, o); + if (dp->chroma == NULL) + return 1; /* grayscale */ + + c_decodeChroma(&cxy, dp->chroma[o*dp->ninc + i]); + c_toSharpRGB(&cxy, coef[0], coef); + coef[0] *= mtx_RGB_coef[0]; + coef[1] *= mtx_RGB_coef[1]; + coef[2] *= mtx_RGB_coef[2]; + return 3; /* RGB color */ +} + /* load custom BSDF angle basis */ static int load_angle_basis(ezxml_t wab) @@ -297,13 +337,12 @@ load_angle_basis(ezxml_t wab) ezxml_child(ezxml_child(wbb, "ThetaBounds"), "UpperTheta"))); if (!i) - abase_list[nabases].lat[i].tmin = - -abase_list[nabases].lat[i+1].tmin; + abase_list[nabases].lat[0].tmin = 0; else if (!fequal(atof(ezxml_txt(ezxml_child(ezxml_child(wbb, "ThetaBounds"), "LowerTheta"))), abase_list[nabases].lat[i].tmin)) { sprintf(SDerrorDetail, "Theta values disagree in '%s'", - abname); + abname); return RC_DATERR; } abase_list[nabases].nangles += @@ -313,7 +352,7 @@ load_angle_basis(ezxml_t wab) (abase_list[nabases].lat[i].nphis == 1 && abase_list[nabases].lat[i].tmin > FTINY)) { sprintf(SDerrorDetail, "Illegal phi count in '%s'", - abname); + abname); return RC_DATERR; } } @@ -348,8 +387,7 @@ get_extrema(SDSpectralDF *df) } free(ohma); /* need incoming solid angles, too? */ - if (dp->ninc < dp->nout || dp->ib_ohm != dp->ob_ohm || - dp->ib_priv != dp->ob_priv) { + if ((dp->ib_ohm != dp->ob_ohm) | (dp->ib_priv != dp->ob_priv)) { double ohm; for (i = dp->ninc; i--; ) if ((ohm = mBSDF_incohm(dp,i)) < df->minProjSA) @@ -358,86 +396,45 @@ get_extrema(SDSpectralDF *df) return (df->maxHemi <= 1.01); } -/* skip integer in string */ -static char * -i_skip(char *s) -{ - while (isspace(*s)) - s++; - if (*s == '-' || *s == '+') - s++; - if (!isdigit(*s)) - return(NULL); - do - s++; - while (isdigit(*s)); - return(s); -} - -/* skip float in string */ -static char * -f_skip(char *s) -{ - register char *cp; - - while (isspace(*s)) - s++; - if (*s == '-' || *s == '+') - s++; - cp = s; - while (isdigit(*cp)) - cp++; - if (*cp == '.') { - cp++; s++; - while (isdigit(*cp)) - cp++; - } - if (cp == s) - return(NULL); - if (*cp == 'e' || *cp == 'E') - return(i_skip(cp+1)); - return(cp); -} - /* load BSDF distribution for this wavelength */ static int -load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc) +load_bsdf_data(SDData *sd, ezxml_t wdb, int ct, int rowinc) { - char *cbasis = ezxml_txt(ezxml_child(wdb,"ColumnAngleBasis")); - char *rbasis = ezxml_txt(ezxml_child(wdb,"RowAngleBasis")); SDSpectralDF *df; SDMat *dp; char *sdata; int inbi, outbi; int i; - - if ((!cbasis || !*cbasis) | (!rbasis || !*rbasis)) { - sprintf(SDerrorDetail, "Missing column/row basis for BSDF '%s'", - sd->name); - return RC_FORMERR; - } /* allocate BSDF component */ sdata = ezxml_txt(ezxml_child(wdb, "WavelengthDataDirection")); + if (!sdata) + return RC_FAIL; + /* + * Remember that front and back are reversed from WINDOW 6 orientations + */ if (!strcasecmp(sdata, "Transmission Front")) { - if (sd->tf != NULL) - SDfreeSpectralDF(sd->tf); - if ((sd->tf = SDnewSpectralDF(1)) == NULL) + if (sd->tb == NULL && (sd->tb = SDnewSpectralDF(3)) == NULL) return RC_MEMERR; + df = sd->tb; + } else if (!strcasecmp(sdata, "Transmission Back")) { + if (sd->tf == NULL && (sd->tf = SDnewSpectralDF(3)) == NULL) + return RC_MEMERR; df = sd->tf; } else if (!strcasecmp(sdata, "Reflection Front")) { - if (sd->rf != NULL) - SDfreeSpectralDF(sd->rf); - if ((sd->rf = SDnewSpectralDF(1)) == NULL) + if (sd->rb == NULL && (sd->rb = SDnewSpectralDF(3)) == NULL) return RC_MEMERR; - df = sd->rf; + df = sd->rb; } else if (!strcasecmp(sdata, "Reflection Back")) { - if (sd->rb != NULL) - SDfreeSpectralDF(sd->rb); - if ((sd->rb = SDnewSpectralDF(1)) == NULL) + if (sd->rf == NULL && (sd->rf = SDnewSpectralDF(3)) == NULL) return RC_MEMERR; - df = sd->rb; + df = sd->rf; } else return RC_FAIL; + /* free previous matrix if any */ + if (df->comp[ct].dist != NULL) { + SDfreeMatrix(df->comp[ct].dist); + df->comp[ct].dist = NULL; + } /* get angle bases */ sdata = ezxml_txt(ezxml_child(wdb,"ColumnAngleBasis")); if (!sdata || !*sdata) { @@ -446,11 +443,10 @@ load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc) return RC_FORMERR; } for (inbi = nabases; inbi--; ) - if (!strcasecmp(cbasis, abase_list[inbi].name)) + if (!strcasecmp(sdata, abase_list[inbi].name)) break; if (inbi < 0) { - sprintf(SDerrorDetail, "Undefined ColumnAngleBasis '%s'", - cbasis); + sprintf(SDerrorDetail, "Undefined ColumnAngleBasis '%s'", sdata); return RC_FORMERR; } sdata = ezxml_txt(ezxml_child(wdb,"RowAngleBasis")); @@ -460,128 +456,219 @@ load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc) return RC_FORMERR; } for (outbi = nabases; outbi--; ) - if (!strcasecmp(rbasis, abase_list[outbi].name)) + if (!strcasecmp(sdata, abase_list[outbi].name)) break; if (outbi < 0) { - sprintf(SDerrorDetail, "Undefined RowAngleBasis '%s'", cbasis); + sprintf(SDerrorDetail, "Undefined RowAngleBasis '%s'", sdata); return RC_FORMERR; } /* allocate BSDF matrix */ dp = SDnewMatrix(abase_list[inbi].nangles, abase_list[outbi].nangles); if (dp == NULL) return RC_MEMERR; - dp->ib_priv = (void *)&abase_list[inbi]; - dp->ob_priv = (void *)&abase_list[outbi]; + dp->ib_priv = &abase_list[inbi]; + dp->ob_priv = &abase_list[outbi]; if (df == sd->tf) { - dp->ib_vec = ab_getvecR; - dp->ib_ndx = ab_getndxR; - dp->ob_vec = ab_getvec; - dp->ob_ndx = ab_getndx; + dp->ib_vec = &fi_getvec; + dp->ib_ndx = &fi_getndx; + dp->ob_vec = &bo_getvec; + dp->ob_ndx = &bo_getndx; + } else if (df == sd->tb) { + dp->ib_vec = &bi_getvec; + dp->ib_ndx = &bi_getndx; + dp->ob_vec = &fo_getvec; + dp->ob_ndx = &fo_getndx; } else if (df == sd->rf) { - dp->ib_vec = ab_getvec; - dp->ib_ndx = ab_getndx; - dp->ob_vec = ab_getvec; - dp->ob_ndx = ab_getndx; + dp->ib_vec = &fi_getvec; + dp->ib_ndx = &fi_getndx; + dp->ob_vec = &fo_getvec; + dp->ob_ndx = &fo_getndx; } else /* df == sd->rb */ { - dp->ib_vec = ab_getvecR; - dp->ib_ndx = ab_getndxR; - dp->ob_vec = ab_getvecR; - dp->ob_ndx = ab_getndxR; + dp->ib_vec = &bi_getvec; + dp->ib_ndx = &bi_getndx; + dp->ob_vec = &bo_getvec; + dp->ob_ndx = &bo_getndx; } - dp->ib_ohm = ab_getohm; - dp->ob_ohm = ab_getohm; - df->comp[0].cspec[0] = c_dfcolor; /* XXX monochrome for now */ - df->comp[0].dist = dp; - df->comp[0].func = &SDhandleMtx; + dp->ib_ohm = &io_getohm; + dp->ob_ohm = &io_getohm; + df->comp[ct].dist = dp; + df->comp[ct].func = &SDhandleMtx; /* read BSDF data */ - sdata = ezxml_txt(ezxml_child(wdb,"ScatteringData")); + sdata = ezxml_txt(ezxml_child(wdb, "ScatteringData")); if (!sdata || !*sdata) { sprintf(SDerrorDetail, "Missing BSDF ScatteringData in '%s'", sd->name); return RC_FORMERR; } for (i = 0; i < dp->ninc*dp->nout; i++) { - char *sdnext = f_skip(sdata); + char *sdnext = fskip(sdata); + double val; if (sdnext == NULL) { sprintf(SDerrorDetail, "Bad/missing BSDF ScatteringData in '%s'", sd->name); return RC_FORMERR; } - while (*sdnext && isspace(*sdnext)) + while (isspace(*sdnext)) sdnext++; if (*sdnext == ',') sdnext++; + if ((val = atof(sdata)) < 0) + val = 0; /* don't allow negative values */ if (rowinc) { int r = i/dp->nout; - int c = i - c*dp->nout; - mBSDF_value(dp,r,c) = atof(sdata); + int c = i - r*dp->nout; + mBSDF_value(dp,r,c) = val; } else - dp->bsdf[i] = atof(sdata); + dp->bsdf[i] = val; sdata = sdnext; } - return get_extrema(df); + return (ct == mtx_Y) ? get_extrema(df) : RC_GOOD; } -/* Subtract minimum (diffuse) scattering amount from BSDF */ +/* copy our RGB (x,y) primary chromaticities */ +static void +copy_RGB_prims(C_COLOR cspec[]) +{ + if (mtx_RGB_coef[1] < .001) { /* need to initialize */ + int i = 3; + while (i--) { + float rgb[3]; + rgb[0] = rgb[1] = rgb[2] = .0f; + rgb[i] = 1.f; + mtx_RGB_coef[i] = c_fromSharpRGB(rgb, &mtx_RGB_prim[i]); + } + } + memcpy(cspec, mtx_RGB_prim, sizeof(mtx_RGB_prim)); +} + +/* encode chromaticity if XYZ -- reduce to one channel in any case */ +static SDSpectralDF * +encode_chroma(SDSpectralDF *df) +{ + SDMat *mpx, *mpy, *mpz; + int n; + + if (df == NULL || df->ncomp != 3) + return df; + + mpy = (SDMat *)df->comp[mtx_Y].dist; + if (mpy == NULL) { + free(df); + return NULL; + } + mpx = (SDMat *)df->comp[mtx_X].dist; + mpz = (SDMat *)df->comp[mtx_Z].dist; + if (mpx == NULL || (mpx->ninc != mpy->ninc) | (mpx->nout != mpy->nout)) + goto done; + if (mpz == NULL || (mpz->ninc != mpy->ninc) | (mpz->nout != mpy->nout)) + goto done; + mpy->chroma = (C_CHROMA *)malloc(sizeof(C_CHROMA)*mpy->ninc*mpy->nout); + if (mpy->chroma == NULL) + goto done; /* XXX punt */ + /* encode chroma values */ + for (n = mpy->ninc*mpy->nout; n--; ) { + const double sum = mpx->bsdf[n] + mpy->bsdf[n] + mpz->bsdf[n]; + C_COLOR cxy; + if (sum > .0) + c_cset(&cxy, mpx->bsdf[n]/sum, mpy->bsdf[n]/sum); + else + c_cset(&cxy, 1./3., 1./3.); + mpy->chroma[n] = c_encodeChroma(&cxy); + } +done: /* free X & Z channels */ + if (mpx != NULL) SDfreeMatrix(mpx); + if (mpz != NULL) SDfreeMatrix(mpz); + if (mpy->chroma == NULL) /* grayscale after all? */ + df->comp[0].cspec[0] = c_dfcolor; + else /* else copy RGB primaries */ + copy_RGB_prims(df->comp[0].cspec); + df->ncomp = 1; /* return resized struct */ + return (SDSpectralDF *)realloc(df, sizeof(SDSpectralDF)); +} + +/* subtract minimum (diffuse) scattering amount from BSDF */ static double -subtract_min(SDMat *sm) +subtract_min(C_COLOR *cs, SDMat *sm) { - float minv = sm->bsdf[0]; - int n = sm->ninc*sm->nout; - int i; + const int ncomp = 1 + 2*(sm->chroma != NULL); + float min_coef[3], ymin, coef[3]; + int i, o, c; - for (i = n; --i; ) - if (sm->bsdf[i] < minv) - minv = sm->bsdf[i]; - for (i = n; i--; ) - sm->bsdf[i] -= minv; + min_coef[0] = min_coef[1] = min_coef[2] = FHUGE; + for (i = 0; i < sm->ninc; i++) + for (o = 0; o < sm->nout; o++) { + c = mBSDF_color(coef, sm, i, o); + while (c--) + if (coef[c] < min_coef[c]) + min_coef[c] = coef[c]; + } + ymin = 0; + for (c = ncomp; c--; ) + ymin += min_coef[c]; + if (ymin <= .01/M_PI) /* not worth bothering about? */ + return .0; + if (ncomp == 1) { /* subtract grayscale minimum */ + for (i = sm->ninc*sm->nout; i--; ) + sm->bsdf[i] -= ymin; + *cs = c_dfcolor; + return M_PI*ymin; + } + /* else subtract colored minimum */ + for (i = 0; i < sm->ninc; i++) + for (o = 0; o < sm->nout; o++) { + C_COLOR cxy; + c = mBSDF_color(coef, sm, i, o); + while (c--) + coef[c] = (coef[c] - min_coef[c]) / + mtx_RGB_coef[c]; + if (c_fromSharpRGB(coef, &cxy) > 1e-5) + sm->chroma[o*sm->ninc + i] = c_encodeChroma(&cxy); + mBSDF_value(sm,i,o) -= ymin; + } + /* return colored minimum */ + for (i = 3; i--; ) + coef[i] = min_coef[i]/mtx_RGB_coef[i]; + c_fromSharpRGB(coef, cs); - return minv*M_PI; /* be sure to include multiplier */ + return M_PI*ymin; } -/* Extract and separate diffuse portion of BSDF */ -static void +/* Extract and separate diffuse portion of BSDF & convert color */ +static SDSpectralDF * extract_diffuse(SDValue *dv, SDSpectralDF *df) { - int n; + df = encode_chroma(df); /* reduce XYZ to Y + chroma */ if (df == NULL || df->ncomp <= 0) { dv->spec = c_dfcolor; dv->cieY = .0; - return; + return df; } - dv->spec = df->comp[0].cspec[0]; - dv->cieY = subtract_min((SDMat *)df->comp[0].dist); - /* in case of multiple components */ - for (n = df->ncomp; --n; ) { - double ymin = subtract_min((SDMat *)df->comp[n].dist); - c_cmix(&dv->spec, dv->cieY, &dv->spec, ymin, &df->comp[n].cspec[0]); - dv->cieY += ymin; - } - df->maxHemi -= dv->cieY; /* correct minimum hemispherical */ - dv->spec.clock++; /* make sure everything is set */ + /* subtract minimum value */ + dv->cieY = subtract_min(&dv->spec, (SDMat *)df->comp[0].dist); + df->maxHemi -= dv->cieY; /* adjust maximum hemispherical */ + /* make sure everything is set */ c_ccvt(&dv->spec, C_CSXY+C_CSSPEC); + return df; } /* Load a BSDF matrix from an open XML file */ SDError -SDloadMtx(SDData *sd, ezxml_t fl) +SDloadMtx(SDData *sd, ezxml_t wtl) { - ezxml_t wtl, wld, wdb; - int rowIn; - struct BSDF_data *dp; - char *txt; - int rval; - - if (strcmp(ezxml_name(fl), "WindowElement")) { + ezxml_t wld, wdb; + int rowIn; + char *txt; + int rval; + /* basic checks and data ordering */ + txt = ezxml_txt(ezxml_child(ezxml_child(wtl, + "DataDefinition"), "IncidentDataStructure")); + if (txt == NULL || !*txt) { sprintf(SDerrorDetail, - "BSDF \"%s\": top level node not 'WindowElement'", + "BSDF \"%s\": missing IncidentDataStructure", sd->name); return SDEformat; } - wtl = ezxml_child(ezxml_child(fl, "Optical"), "Layer"); - txt = ezxml_txt(ezxml_child(ezxml_child(wtl, - "DataDefinition"), "IncidentDataStructure")); if (!strcasecmp(txt, "Rows")) rowIn = 1; else if (!strcasecmp(txt, "Columns")) @@ -592,51 +679,53 @@ SDloadMtx(SDData *sd, ezxml_t fl) sd->name); return SDEsupport; } - /* get angle basis */ - rval = load_angle_basis(ezxml_child(ezxml_child(wtl, - "DataDefinition"), "AngleBasis")); - if (rval < 0) - goto err_return; - /* load BSDF components */ + /* get angle bases */ + for (wld = ezxml_child(ezxml_child(wtl, "DataDefinition"), "AngleBasis"); + wld != NULL; wld = wld->next) { + rval = load_angle_basis(wld); + if (rval < 0) + return convert_errcode(rval); + } + /* load BSDF components */ for (wld = ezxml_child(wtl, "WavelengthData"); wld != NULL; wld = wld->next) { - if (strcasecmp(ezxml_txt(ezxml_child(wld,"Wavelength")), - "Visible")) - continue; /* just visible for now */ + const char *cnm = ezxml_txt(ezxml_child(wld,"Wavelength")); + int ct = -1; + if (!strcasecmp(cnm, "Visible")) + ct = mtx_Y; + else if (!strcasecmp(cnm, "CIE-X")) + ct = mtx_X; + else if (!strcasecmp(cnm, "CIE-Z")) + ct = mtx_Z; + else + continue; for (wdb = ezxml_child(wld, "WavelengthDataBlock"); wdb != NULL; wdb = wdb->next) - if ((rval = load_bsdf_data(sd, wdb, rowIn)) < 0) - goto err_return; + if ((rval = load_bsdf_data(sd, wdb, ct, rowIn)) < 0) + return convert_errcode(rval); } - /* separate diffuse components */ - extract_diffuse(&sd->rLambFront, sd->rf); - extract_diffuse(&sd->rLambBack, sd->rb); - extract_diffuse(&sd->tLamb, sd->tf); - /* return success */ + /* separate diffuse components */ + sd->rf = extract_diffuse(&sd->rLambFront, sd->rf); + sd->rb = extract_diffuse(&sd->rLambBack, sd->rb); + if (sd->tf != NULL) + sd->tf = extract_diffuse(&sd->tLamb, sd->tf); + if (sd->tb != NULL) + sd->tb = extract_diffuse(&sd->tLamb, sd->tb); + /* return success */ return SDEnone; -err_return: /* jump here on failure */ - if (sd->rf != NULL) { - SDfreeSpectralDF(sd->rf); - sd->rf = NULL; - } - if (sd->rb != NULL) { - SDfreeSpectralDF(sd->rb); - sd->rb = NULL; - } - if (sd->tf != NULL) { - SDfreeSpectralDF(sd->tf); - sd->tf = NULL; - } - return convert_errcode(rval); } /* Get Matrix BSDF value */ static int SDgetMtxBSDF(float coef[SDmaxCh], const FVECT outVec, - const FVECT inVec, const void *dist) + const FVECT inVec, SDComponent *sdc) { - const SDMat *dp = (const SDMat *)dist; + const SDMat *dp; int i_ndx, o_ndx; + /* check arguments */ + if ((coef == NULL) | (outVec == NULL) | (inVec == NULL) | (sdc == NULL) + || (dp = (SDMat *)sdc->dist) == NULL) + return 0; /* get angle indices */ i_ndx = mBSDF_incndx(dp, inVec); o_ndx = mBSDF_outndx(dp, outVec); @@ -647,63 +736,57 @@ SDgetMtxBSDF(float coef[SDmaxCh], const FVECT outVec, } if ((i_ndx < 0) | (o_ndx < 0)) return 0; /* nothing from this component */ - coef[0] = mBSDF_value(dp, i_ndx, o_ndx); - return 1; /* XXX monochrome for now */ + + return mBSDF_color(coef, dp, i_ndx, o_ndx); } -/* Query solid angle for vector */ +/* Query solid angle for vector(s) */ static SDError -SDqueryMtxProjSA(double *psa, const FVECT vec, int qflags, const void *dist) +SDqueryMtxProjSA(double *psa, const FVECT v1, const RREAL *v2, + int qflags, SDComponent *sdc) { - const SDMat *dp = (const SDMat *)dist; - - if (!(qflags & SDqueryInc+SDqueryOut)) + const SDMat *dp; + double inc_psa, out_psa; + /* check arguments */ + if ((psa == NULL) | (v1 == NULL) | (sdc == NULL) || + (dp = (SDMat *)sdc->dist) == NULL) return SDEargument; - if (qflags & SDqueryInc) { - double inc_psa = mBSDF_incohm(dp, mBSDF_incndx(dp, vec)); - if (inc_psa < .0) - return SDEinternal; - switch (qflags & SDqueryMin+SDqueryMax) { - case SDqueryMax: - if (inc_psa > psa[0]) - psa[0] = inc_psa; - break; - case SDqueryMin+SDqueryMax: - if (inc_psa > psa[1]) - psa[1] = inc_psa; - /* fall through */ - case SDqueryMin: - if (inc_psa < psa[0]) - psa[0] = inc_psa; - break; - case 0: + if (v2 == NULL) + v2 = v1; + /* get projected solid angles */ + out_psa = mBSDF_outohm(dp, mBSDF_outndx(dp, v1)); + inc_psa = mBSDF_incohm(dp, mBSDF_incndx(dp, v2)); + if ((v1 != v2) & (out_psa <= 0) & (inc_psa <= 0)) { + inc_psa = mBSDF_outohm(dp, mBSDF_outndx(dp, v2)); + out_psa = mBSDF_incohm(dp, mBSDF_incndx(dp, v1)); + } + + switch (qflags) { /* record based on flag settings */ + case SDqueryMax: + if (inc_psa > psa[0]) psa[0] = inc_psa; - break; - } + if (out_psa > psa[0]) + psa[0] = out_psa; + break; + case SDqueryMin+SDqueryMax: + if (inc_psa > psa[1]) + psa[1] = inc_psa; + if (out_psa > psa[1]) + psa[1] = out_psa; + /* fall through */ + case SDqueryVal: + if (qflags == SDqueryVal) + psa[0] = M_PI; + /* fall through */ + case SDqueryMin: + if ((inc_psa > 0) & (inc_psa < psa[0])) + psa[0] = inc_psa; + if ((out_psa > 0) & (out_psa < psa[0])) + psa[0] = out_psa; + break; } - if (qflags & SDqueryOut) { - double out_psa = mBSDF_outohm(dp, mBSDF_outndx(dp, vec)); - if (out_psa < .0) - return SDEinternal; - switch (qflags & SDqueryMin+SDqueryMax) { - case SDqueryMax: - if (out_psa > psa[0]) - psa[0] = out_psa; - break; - case SDqueryMin+SDqueryMax: - if (out_psa > psa[1]) - psa[1] = out_psa; - /* fall through */ - case SDqueryMin: - if (out_psa < psa[0]) - psa[0] = out_psa; - break; - case 0: - psa[(qflags&SDqueryInc)!=0] = out_psa; - break; - } - } - return SDEnone; + /* make sure it's legal */ + return (psa[0] <= 0) ? SDEinternal : SDEnone; } /* Compute new cumulative distribution from BSDF */ @@ -741,12 +824,13 @@ make_cdist(SDMatCDst *cd, const FVECT inVec, SDMat *dp static const SDCDst * SDgetMtxCDist(const FVECT inVec, SDComponent *sdc) { - SDMat *dp = (SDMat *)sdc->dist; + SDMat *dp; int reverse; SDMatCDst myCD; SDMatCDst *cd, *cdlast; - - if (dp == NULL) + /* check arguments */ + if ((inVec == NULL) | (sdc == NULL) || + (dp = (SDMat *)sdc->dist) == NULL) return NULL; memset(&myCD, 0, sizeof(myCD)); myCD.indx = mBSDF_incndx(dp, inVec); @@ -765,17 +849,15 @@ SDgetMtxCDist(const FVECT inVec, SDComponent *sdc) reverse = 1; } cdlast = NULL; /* check for it in cache list */ - for (cd = (SDMatCDst *)sdc->cdList; - cd != NULL; cd = (SDMatCDst *)cd->next) { + for (cd = (SDMatCDst *)sdc->cdList; cd != NULL; + cdlast = cd, cd = cd->next) if (cd->indx == myCD.indx && (cd->calen == myCD.calen) & (cd->ob_priv == myCD.ob_priv) & (cd->ob_vec == myCD.ob_vec)) break; - cdlast = cd; - } if (cd == NULL) { /* need to allocate new entry */ cd = (SDMatCDst *)malloc(sizeof(SDMatCDst) + - myCD.calen*sizeof(myCD.carr[0])); + sizeof(myCD.carr[0])*myCD.calen); if (cd == NULL) return NULL; *cd = myCD; /* compute cumulative distribution */ @@ -787,7 +869,7 @@ SDgetMtxCDist(const FVECT inVec, SDComponent *sdc) } if (cdlast != NULL) { /* move entry to head of cache list */ cdlast->next = cd->next; - cd->next = sdc->cdList; + cd->next = (SDMatCDst *)sdc->cdList; sdc->cdList = (SDCDst *)cd; } return (SDCDst *)cd; /* ready to go */ @@ -795,16 +877,19 @@ SDgetMtxCDist(const FVECT inVec, SDComponent *sdc) /* Sample cumulative distribution */ static SDError -SDsampMtxCDist(FVECT outVec, double randX, const SDCDst *cdp) +SDsampMtxCDist(FVECT ioVec, double randX, const SDCDst *cdp) { const unsigned maxval = ~0; const SDMatCDst *mcd = (const SDMatCDst *)cdp; const unsigned target = randX*maxval; int i, iupper, ilower; + /* check arguments */ + if ((ioVec == NULL) | (mcd == NULL)) + return SDEargument; /* binary search to find index */ ilower = 0; iupper = mcd->calen; while ((i = (iupper + ilower) >> 1) != ilower) - if ((long)target >= (long)mcd->carr[i]) + if (target >= mcd->carr[i]) ilower = i; else iupper = i; @@ -812,9 +897,9 @@ SDsampMtxCDist(FVECT outVec, double randX, const SDCDs randX = (randX*maxval - mcd->carr[ilower]) / (double)(mcd->carr[iupper] - mcd->carr[ilower]); /* convert index to vector */ - if ((*mcd->ob_vec)(outVec, i, randX, mcd->ob_priv)) + if ((*mcd->ob_vec)(ioVec, i+randX, mcd->ob_priv)) return SDEnone; - strcpy(SDerrorDetail, "BSDF sampling fault"); + strcpy(SDerrorDetail, "Matrix BSDF sampling fault"); return SDEinternal; }