--- ray/src/common/bsdf_m.c 2011/02/18 00:41:44 3.2 +++ ray/src/common/bsdf_m.c 2011/04/08 18:13:48 3.9 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: bsdf_m.c,v 3.2 2011/02/18 00:41:44 greg Exp $"; +static const char RCSid[] = "$Id: bsdf_m.c,v 3.9 2011/04/08 18:13:48 greg Exp $"; #endif /* * bsdf_m.c @@ -11,19 +11,14 @@ static const char RCSid[] = "$Id: bsdf_m.c,v 3.2 2011/ * */ -#include +#include "rtio.h" #include #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 @@ -91,7 +86,7 @@ fequal(double a, double b) return (a <= 1e-6) & (a >= -1e-6); } -/* returns the name of the given tag */ +/* Returns the name of the given tag */ #ifdef ezxml_name #undef ezxml_name static char * @@ -103,7 +98,7 @@ ezxml_name(ezxml_t xml) } #endif -/* returns the given tag's character content or empty string if none */ +/* Returns the given tag's character content or empty string if none */ #ifdef ezxml_txt #undef ezxml_txt static char * @@ -162,16 +157,18 @@ SDnewMatrix(int ni, int no) /* Free a BSDF matrix */ #define SDfreeMatrix free -/* get vector for this angle basis index */ +/* get vector for this angle basis index (front exiting) */ static int -ab_getvec(FVECT v, int ndx, double randX, void *p) +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; - 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; @@ -186,11 +183,11 @@ ab_getvec(FVECT v, int ndx, double randX, void *p) return RC_GOOD; } -/* get index corresponding to the given vector */ +/* get index corresponding to the given vector (front exiting) */ static int -ab_getndx(const FVECT v, void *p) +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; @@ -215,9 +212,9 @@ ab_getndx(const FVECT v, void *p) /* compute square of real value */ static double sq(double x) { return x*x; } -/* get projected solid angle for this angle basis index */ +/* get projected solid angle for this angle basis index (universal) */ static double -ab_getohm(int ndx, void *p) +io_getohm(int ndx, void *p) { static int last_li = -1; static double last_ohm; @@ -240,13 +237,11 @@ ab_getohm(int ndx, void *p) (double)ab->lat[li].nphis; } -/* get reverse vector for this angle basis index */ +/* get vector for this angle basis index (back incident) */ static int -ab_getvecR(FVECT v, int ndx, double randX, void *p) +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]; @@ -256,9 +251,9 @@ ab_getvecR(FVECT v, int ndx, double randX, void *p) return RC_GOOD; } -/* get index corresponding to the reverse vector */ +/* get index corresponding to the vector (back incident) */ static int -ab_getndxR(const FVECT v, void *p) +bi_getndx(const FVECT v, void *p) { FVECT v2; @@ -266,9 +261,60 @@ 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) */ +static 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) */ +static 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) */ +static 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) */ +static 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); +} + /* load custom BSDF angle basis */ static int load_angle_basis(ezxml_t wab) @@ -351,8 +397,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) @@ -361,86 +406,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) { - char *cbasis = ezxml_txt(ezxml_child(wdb,"ColumnAngleBasis")); - char *rbasis = ezxml_txt(ezxml_child(wdb,"RowAngleBasis")); SDSpectralDF *df; SDMat *dp; char *sdata; + int tfront; 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 (!strcasecmp(sdata, "Transmission Front")) { + /* + * 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) 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) - return RC_MEMERR; - df = sd->rf; - } else if (!strcasecmp(sdata, "Reflection Back")) { - if (sd->rb != NULL) + if (sd->rb != NULL) /* note back-front reversal */ SDfreeSpectralDF(sd->rb); if ((sd->rb = SDnewSpectralDF(1)) == 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) + return RC_MEMERR; + df = sd->rf; } else return RC_FAIL; + /* XXX should also check "ScatteringDataType" for consistency? */ /* get angle bases */ sdata = ezxml_txt(ezxml_child(wdb,"ColumnAngleBasis")); if (!sdata || !*sdata) { @@ -449,11 +453,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")); @@ -463,36 +466,43 @@ 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; + 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; + } } 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; + 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; @@ -504,7 +514,7 @@ load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc) return RC_FORMERR; } for (i = 0; i < dp->ninc*dp->nout; i++) { - char *sdnext = f_skip(sdata); + char *sdnext = fskip(sdata); if (sdnext == NULL) { sprintf(SDerrorDetail, "Bad/missing BSDF ScatteringData in '%s'", @@ -561,30 +571,29 @@ extract_diffuse(SDValue *dv, SDSpectralDF *df) 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 */ + df->maxHemi -= dv->cieY; /* adjust minimum hemispherical */ + /* make sure everything is set */ c_ccvt(&dv->spec, C_CSXY+C_CSSPEC); } /* 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; + ezxml_t wld, wdb; int rowIn; struct BSDF_data *dp; char *txt; int rval; - if (strcmp(ezxml_name(fl), "WindowElement")) { + 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")) @@ -599,7 +608,7 @@ SDloadMtx(SDData *sd, ezxml_t fl) rval = load_angle_basis(ezxml_child(ezxml_child(wtl, "DataDefinition"), "AngleBasis")); if (rval < 0) - goto err_return; + return convert_errcode(rval); /* load BSDF components */ for (wld = ezxml_child(wtl, "WavelengthData"); wld != NULL; wld = wld->next) { @@ -609,7 +618,7 @@ SDloadMtx(SDData *sd, ezxml_t fl) for (wdb = ezxml_child(wld, "WavelengthDataBlock"); wdb != NULL; wdb = wdb->next) if ((rval = load_bsdf_data(sd, wdb, rowIn)) < 0) - goto err_return; + return convert_errcode(rval); } /* separate diffuse components */ extract_diffuse(&sd->rLambFront, sd->rf); @@ -617,20 +626,6 @@ SDloadMtx(SDData *sd, ezxml_t fl) extract_diffuse(&sd->tLamb, sd->tf); /* 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 */ @@ -659,54 +654,39 @@ static SDError SDqueryMtxProjSA(double *psa, const FVECT vec, int qflags, const void *dist) { const SDMat *dp = (const SDMat *)dist; - - if (!(qflags & SDqueryInc+SDqueryOut)) + double inc_psa, out_psa; + /* check arguments */ + if ((psa == NULL) | (vec == NULL) | (dp == 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: + /* get projected solid angles */ + inc_psa = mBSDF_incohm(dp, mBSDF_incndx(dp, vec)); + out_psa = mBSDF_outohm(dp, mBSDF_outndx(dp, vec)); + + switch (qflags) { /* record based on flag settings */ + case SDqueryVal: + psa[0] = .0; + /* fall through */ + 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[0]) + psa[1] = inc_psa; + if (out_psa > psa[0]) + psa[1] = out_psa; + /* 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 */ @@ -748,8 +728,8 @@ SDgetMtxCDist(const FVECT inVec, SDComponent *sdc) int reverse; SDMatCDst myCD; SDMatCDst *cd, *cdlast; - - if (dp == NULL) + /* check arguments */ + if ((inVec == NULL) | (dp == NULL)) return NULL; memset(&myCD, 0, sizeof(myCD)); myCD.indx = mBSDF_incndx(dp, inVec); @@ -804,6 +784,9 @@ SDsampMtxCDist(FVECT outVec, double randX, const SDCDs const SDMatCDst *mcd = (const SDMatCDst *)cdp; const unsigned target = randX*maxval; int i, iupper, ilower; + /* check arguments */ + if ((outVec == NULL) | (mcd == NULL)) + return SDEargument; /* binary search to find index */ ilower = 0; iupper = mcd->calen; while ((i = (iupper + ilower) >> 1) != ilower) @@ -815,7 +798,7 @@ 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)(outVec, i+randX, mcd->ob_priv)) return SDEnone; strcpy(SDerrorDetail, "BSDF sampling fault"); return SDEinternal;