--- ray/src/common/bsdf_m.c 2011/04/24 20:16:52 3.13 +++ ray/src/common/bsdf_m.c 2022/01/25 01:34:20 3.45 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: bsdf_m.c,v 3.13 2011/04/24 20:16:52 greg Exp $"; +static const char RCSid[] = "$Id: bsdf_m.c,v 3.45 2022/01/25 01:34:20 greg Exp $"; #endif /* * bsdf_m.c @@ -11,8 +11,8 @@ static const char RCSid[] = "$Id: bsdf_m.c,v 3.13 2011 * */ +#define _USE_MATH_DEFINES #include "rtio.h" -#include #include #include #include "ezxml.h" @@ -28,24 +28,10 @@ static const char RCSid[] = "$Id: bsdf_m.c,v 3.13 2011 #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}, @@ -57,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}, @@ -67,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}, @@ -76,8 +62,14 @@ 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) { @@ -86,31 +78,7 @@ fequal(double a, double b) 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) { @@ -131,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) { @@ -155,10 +123,20 @@ 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 (front exiting) */ -static int + 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; @@ -166,36 +144,36 @@ fo_getvec(FVECT v, double ndxr, void *p) double randX = ndxr - ndx; double rx[2]; int li; - double pol, azi, d; + double azi, d; 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. - rx[0])*sq(cos(M_PI/180.*ab->lat[li].tmin)) + + rx[0]*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 (front exiting) */ -static int +/* Get index corresponding to the given vector (front exiting) */ +int fo_getndx(const FVECT v, void *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.)) + if ((v[2] < 0) | (v[2] > 1.00001)) 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,13 +187,11 @@ fo_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 (universal) */ -static double +/* Get projected solid angle for this angle basis index (universal) */ +double io_getohm(int ndx, void *p) { + static void *last_p = NULL; static int last_li = -1; static double last_ohm; ANGLE_BASIS *ab = (ANGLE_BASIS *)p; @@ -226,19 +202,18 @@ io_getohm(int ndx, void *p) return -1.; for (li = 0; ndx >= ab->lat[li].nphis; li++) ndx -= ab->lat[li].nphis; - if (li == last_li) /* cached latitude? */ + if ((p == last_p) & (li == last_li)) /* cached latitude? */ return last_ohm; + last_p = p; 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 vector for this angle basis index (back incident) */ -static int +/* Get vector for this angle basis index (back incident) */ +int bi_getvec(FVECT v, double ndxr, void *p) { if (!fo_getvec(v, ndxr, p)) @@ -251,8 +226,8 @@ bi_getvec(FVECT v, double ndxr, void *p) return RC_GOOD; } -/* get index corresponding to the vector (back incident) */ -static int +/* Get index corresponding to the vector (back incident) */ +int bi_getndx(const FVECT v, void *p) { FVECT v2; @@ -264,8 +239,8 @@ bi_getndx(const FVECT v, void *p) return fo_getndx(v2, p); } -/* get vector for this angle basis index (back exiting) */ -static int +/* Get vector for this angle basis index (back exiting) */ +int bo_getvec(FVECT v, double ndxr, void *p) { if (!fo_getvec(v, ndxr, p)) @@ -276,8 +251,8 @@ bo_getvec(FVECT v, double ndxr, void *p) return RC_GOOD; } -/* get index corresponding to the vector (back exiting) */ -static int +/* Get index corresponding to the vector (back exiting) */ +int bo_getndx(const FVECT v, void *p) { FVECT v2; @@ -289,8 +264,8 @@ bo_getndx(const FVECT v, void *p) return fo_getndx(v2, p); } -/* get vector for this angle basis index (front incident) */ -static int +/* Get vector for this angle basis index (front incident) */ +int fi_getvec(FVECT v, double ndxr, void *p) { if (!fo_getvec(v, ndxr, p)) @@ -302,8 +277,8 @@ fi_getvec(FVECT v, double ndxr, void *p) return RC_GOOD; } -/* get index corresponding to the vector (front incident) */ -static int +/* Get index corresponding to the vector (front incident) */ +int fi_getndx(const FVECT v, void *p) { FVECT v2; @@ -315,6 +290,29 @@ fi_getndx(const FVECT v, void *p) 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; + double d; + + coef[0] = mBSDF_value(dp, o, i); + /* position-specific perturbation */ + d = 2*dp->ninc/(i + .22545) + 4*dp->nout/(o + .70281); + d -= (int)d; + coef[0] *= 1. + 6e-4*(d - .5); + if (dp->chroma == NULL) + return 1; /* grayscale */ + + c_decodeChroma(&cxy, mBSDF_chroma(dp,o,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) @@ -346,8 +344,7 @@ 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)) { @@ -391,7 +388,7 @@ get_extrema(SDSpectralDF *df) for (i = dp->ninc; i--; ) { double hemi = .0; for (o = dp->nout; o--; ) - hemi += ohma[o] * mBSDF_value(dp, i, o); + hemi += ohma[o] * mBSDF_value(dp, o, i); if (hemi > df->maxHemi) df->maxHemi = hemi; } @@ -408,43 +405,43 @@ get_extrema(SDSpectralDF *df) /* 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) { SDSpectralDF *df; SDMat *dp; char *sdata; - int tfront; int inbi, outbi; int i; /* 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 - * Favor their "Front" (incoming light) since that's more often valid */ - tfront = !strcasecmp(sdata, "Transmission Back"); - if (!strcasecmp(sdata, "Transmission Front") || - tfront & (sd->tf == NULL)) { - if (sd->tf != NULL) - SDfreeSpectralDF(sd->tf); - if ((sd->tf = SDnewSpectralDF(1)) == NULL) + if (!strcasecmp(sdata, "Transmission Front")) { + 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->rb != NULL) /* note back-front reversal */ - SDfreeSpectralDF(sd->rb); - if ((sd->rb = SDnewSpectralDF(1)) == NULL) + if (sd->rb == NULL && (sd->rb = SDnewSpectralDF(3)) == NULL) return RC_MEMERR; df = sd->rb; } else if (!strcasecmp(sdata, "Reflection Back")) { - if (sd->rf != NULL) /* note front-back reversal */ - SDfreeSpectralDF(sd->rf); - if ((sd->rf = SDnewSpectralDF(1)) == NULL) + if (sd->rf == NULL && (sd->rf = SDnewSpectralDF(3)) == NULL) return RC_MEMERR; df = sd->rf; } else return RC_FAIL; - /* XXX should also check "ScatteringDataType" for consistency? */ + /* 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) { @@ -479,17 +476,15 @@ load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc) dp->ib_priv = &abase_list[inbi]; dp->ob_priv = &abase_list[outbi]; if (df == sd->tf) { - if (tfront) { - dp->ib_vec = &fi_getvec; - dp->ib_ndx = &fi_getndx; - dp->ob_vec = &bo_getvec; - dp->ob_ndx = &bo_getndx; - } else { - dp->ib_vec = &bi_getvec; - dp->ib_ndx = &bi_getndx; - dp->ob_vec = &fo_getvec; - dp->ob_ndx = &fo_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 = &fi_getvec; dp->ib_ndx = &fi_getndx; @@ -503,89 +498,176 @@ load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc) } dp->ib_ohm = &io_getohm; dp->ob_ohm = &io_getohm; - df->comp[0].cspec[0] = c_dfcolor; /* XXX monochrome for now */ - df->comp[0].dist = dp; - df->comp[0].func = &SDhandleMtx; + 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 = fskip(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,c,r) = 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) + mBSDF_chroma(sm,o,i) = c_encodeChroma(&cxy); + mBSDF_value(sm,o,i) -= 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; /* adjust minimum hemispherical */ - /* make sure everything is set */ - c_ccvt(&dv->spec, C_CSXY+C_CSSPEC); + /* subtract minimum value */ + dv->cieY = subtract_min(&dv->spec, (SDMat *)df->comp[0].dist); + df->maxHemi -= dv->cieY; /* adjust maximum hemispherical */ + + c_ccvt(&dv->spec, C_CSXY); /* make sure (x,y) is set */ + return df; } /* Load a BSDF matrix from an open XML file */ SDError SDloadMtx(SDData *sd, ezxml_t wtl) { - ezxml_t wld, wdb; - int rowIn; - struct BSDF_data *dp; - char *txt; - int rval; - + 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) { @@ -604,34 +686,49 @@ SDloadMtx(SDData *sd, ezxml_t wtl) sd->name); return SDEsupport; } - /* get angle basis */ - rval = load_angle_basis(ezxml_child(ezxml_child(wtl, - "DataDefinition"), "AngleBasis")); - if (rval < 0) - return convert_errcode(rval); - /* 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) + 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); + sd->tf = extract_diffuse(&sd->tLambFront, sd->tf); + if (sd->tb != NULL) { + sd->tb = extract_diffuse(&sd->tLambBack, sd->tb); + if (sd->tf == NULL) + sd->tLambFront = sd->tLambBack; + } else if (sd->tf != NULL) + sd->tLambBack = sd->tLambFront; + /* return success */ return SDEnone; } /* Get Matrix BSDF value */ static int -SDgetMtxBSDF(float coef[SDmaxCh], const FVECT outVec, - const FVECT inVec, SDComponent *sdc) +SDgetMtxBSDF(float coef[SDmaxCh], const FVECT inVec, + const FVECT outVec, SDComponent *sdc) { const SDMat *dp; int i_ndx, o_ndx; @@ -649,8 +746,8 @@ 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(s) */ @@ -687,10 +784,11 @@ SDqueryMtxProjSA(double *psa, const FVECT v1, const RR if (out_psa > psa[1]) psa[1] = out_psa; /* fall through */ - case SDqueryMin: 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])) @@ -715,10 +813,10 @@ make_cdist(SDMatCDst *cd, const FVECT inVec, SDMat *dp cmtab[0] = .0; for (o = 0; o < cd->calen; o++) { if (rev) - cmtab[o+1] = mBSDF_value(dp, o, cd->indx) * + cmtab[o+1] = mBSDF_value(dp, cd->indx, o) * (*dp->ib_ohm)(o, dp->ib_priv); else - cmtab[o+1] = mBSDF_value(dp, cd->indx, o) * + cmtab[o+1] = mBSDF_value(dp, o, cd->indx) * (*dp->ob_ohm)(o, dp->ob_priv); cmtab[o+1] += cmtab[o]; } @@ -761,17 +859,16 @@ 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) { + /* PLACE MUTEX LOCK HERE FOR THREAD-SAFE */ + 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 */ @@ -783,9 +880,10 @@ 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; } + /* END MUTEX LOCK */ return (SDCDst *)cd; /* ready to go */ } @@ -803,7 +901,7 @@ SDsampMtxCDist(FVECT ioVec, double randX, const SDCDst /* 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; @@ -818,7 +916,7 @@ SDsampMtxCDist(FVECT ioVec, double randX, const SDCDst } /* Fixed resolution BSDF methods */ -SDFunc SDhandleMtx = { +const SDFunc SDhandleMtx = { &SDgetMtxBSDF, &SDqueryMtxProjSA, &SDgetMtxCDist,