--- ray/src/common/bsdf_m.c 2014/03/21 17:49:53 3.28 +++ ray/src/common/bsdf_m.c 2021/03/27 17:50:18 3.41 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: bsdf_m.c,v 3.28 2014/03/21 17:49:53 greg Exp $"; +static const char RCSid[] = "$Id: bsdf_m.c,v 3.41 2021/03/27 17:50:18 greg Exp $"; #endif /* * bsdf_m.c @@ -13,7 +13,6 @@ static const char RCSid[] = "$Id: bsdf_m.c,v 3.28 2014 #define _USE_MATH_DEFINES #include "rtio.h" -#include #include #include #include "ezxml.h" @@ -65,6 +64,12 @@ ANGLE_BASIS abase_list[MAXABASES] = { 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) { @@ -73,19 +78,7 @@ fequal(double a, double b) return (a <= 1e-6) & (a >= -1e-6); } -/* 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) { @@ -106,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) { @@ -130,8 +123,18 @@ SDnewMatrix(int ni, int no) } /* Free a BSDF matrix */ -#define SDfreeMatrix free +void +SDfreeMatrix(void *ptr) +{ + SDMat *mp = (SDMat *)ptr; + 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) @@ -141,17 +144,17 @@ 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; @@ -184,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) */ 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; @@ -201,8 +202,9 @@ 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; theta = M_PI/180. * ab->lat[li].tmin; theta1 = M_PI/180. * ab->lat[li+1].tmin; @@ -288,6 +290,24 @@ 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; + + coef[0] = mBSDF_value(dp, o, i); + 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) @@ -363,7 +383,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; } @@ -380,7 +400,7 @@ 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; @@ -395,32 +415,28 @@ load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc) * Remember that front and back are reversed from WINDOW 6 orientations */ if (!strcasecmp(sdata, "Transmission Front")) { - if (sd->tb != NULL) - SDfreeSpectralDF(sd->tb); - if ((sd->tb = 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) - SDfreeSpectralDF(sd->tf); - if ((sd->tf = SDnewSpectralDF(1)) == NULL) + 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) - 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) - 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) { @@ -477,9 +493,8 @@ 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")); if (!sdata || !*sdata) { @@ -504,57 +519,139 @@ load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc) if (rowinc) { int r = i/dp->nout; int c = i - r*dp->nout; - mBSDF_value(dp,r,c) = val; + mBSDF_value(dp,c,r) = val; } else 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]; - - if (minv <= FTINY) + 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); - for (i = n; i--; ) - sm->bsdf[i] -= minv; - - 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; - } + /* 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); + + c_ccvt(&dv->spec, C_CSXY); /* make sure (x,y) is set */ + return df; } /* Load a BSDF matrix from an open XML file */ @@ -594,21 +691,26 @@ SDloadMtx(SDData *sd, ezxml_t wtl) /* 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); - if (sd->tf != NULL) - extract_diffuse(&sd->tLamb, sd->tf); - if (sd->tb != NULL) - extract_diffuse(&sd->tLamb, sd->tb); + sd->rf = extract_diffuse(&sd->rLambFront, sd->rf); + sd->rb = extract_diffuse(&sd->rLambBack, sd->rb); + sd->tf = extract_diffuse(&sd->tLambFront, sd->tf); + sd->tb = extract_diffuse(&sd->tLambBack, sd->tb); /* return success */ return SDEnone; } @@ -634,8 +736,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) */ @@ -701,10 +803,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]; } @@ -747,6 +849,7 @@ SDgetMtxCDist(const FVECT inVec, SDComponent *sdc) reverse = 1; } cdlast = NULL; /* check for it in cache list */ + /* 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) & @@ -770,6 +873,7 @@ SDgetMtxCDist(const FVECT inVec, SDComponent *sdc) cd->next = (SDMatCDst *)sdc->cdList; sdc->cdList = (SDCDst *)cd; } + /* END MUTEX LOCK */ return (SDCDst *)cd; /* ready to go */ } @@ -802,7 +906,7 @@ SDsampMtxCDist(FVECT ioVec, double randX, const SDCDst } /* Fixed resolution BSDF methods */ -SDFunc SDhandleMtx = { +const SDFunc SDhandleMtx = { &SDgetMtxBSDF, &SDqueryMtxProjSA, &SDgetMtxCDist,