--- ray/src/common/bsdf_m.c 2011/02/21 22:50:37 3.6 +++ ray/src/common/bsdf_m.c 2013/04/21 22:58:40 3.25 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: bsdf_m.c,v 3.6 2011/02/21 22:50:37 greg Exp $"; +static const char RCSid[] = "$Id: bsdf_m.c,v 3.25 2013/04/21 22:58:40 greg Exp $"; #endif /* * bsdf_m.c @@ -11,6 +11,7 @@ static const char RCSid[] = "$Id: bsdf_m.c,v 3.6 2011/ * */ +#define _USE_MATH_DEFINES #include "rtio.h" #include #include @@ -28,24 +29,10 @@ static const char RCSid[] = "$Id: bsdf_m.c,v 3.6 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 +44,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 +54,7 @@ static ANGLE_BASIS abase_list[MAXABASES] = { {90., 0} } }, { "LBNL/Klems Quarter", 41, - { {-9., 1}, + { {0., 1}, {9., 8}, {27., 12}, {46., 12}, @@ -76,29 +63,17 @@ static ANGLE_BASIS abase_list[MAXABASES] = { } }; -static int nabases = 3; /* current number of defined bases */ +int nabases = 3; /* current number of defined bases */ 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 */ +/* Returns the given tag's character content or empty string if none */ #ifdef ezxml_txt #undef ezxml_txt static char * @@ -157,16 +132,18 @@ SDnewMatrix(int ni, int no) /* Free a BSDF matrix */ #define SDfreeMatrix free -/* get vector for this angle basis index */ -static int -ab_getvec(FVECT v, int ndx, double randX, void *p) +/* 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; - 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; @@ -181,17 +158,17 @@ ab_getvec(FVECT v, int ndx, double randX, void *p) 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]); azi = 180.0/M_PI*atan2(v[1], v[0]); @@ -210,9 +187,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 */ -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; @@ -227,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]; @@ -251,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; @@ -261,9 +234,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) */ +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); +} + /* load custom BSDF angle basis */ static int load_angle_basis(ezxml_t wab) @@ -295,13 +319,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 += @@ -311,7 +334,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; } } @@ -366,21 +389,31 @@ load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc) int i; /* allocate BSDF component */ sdata = ezxml_txt(ezxml_child(wdb, "WavelengthDataDirection")); - if (!strcasecmp(sdata, "Transmission Back") || (sd->tf == NULL && - !strcasecmp(sdata, "Transmission Front"))) { + if (!sdata) + return RC_FAIL; + /* + * 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) + 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) return RC_MEMERR; df = sd->tf; } else if (!strcasecmp(sdata, "Reflection Front")) { - if (sd->rb != NULL) /* note back-front reversal */ + if (sd->rb != NULL) 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 */ + if (sd->rf != NULL) SDfreeSpectralDF(sd->rf); if ((sd->rf = SDnewSpectralDF(1)) == NULL) return RC_MEMERR; @@ -422,28 +455,33 @@ 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) { - 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; + 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; /* 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); @@ -457,12 +495,12 @@ load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc) sd->name); return RC_FORMERR; } - while (*sdnext && isspace(*sdnext)) + while (isspace(*sdnext)) sdnext++; if (*sdnext == ',') sdnext++; if (rowinc) { int r = i/dp->nout; - int c = i - c*dp->nout; + int c = i - r*dp->nout; mBSDF_value(dp,r,c) = atof(sdata); } else dp->bsdf[i] = atof(sdata); @@ -482,6 +520,10 @@ subtract_min(SDMat *sm) for (i = n; --i; ) if (sm->bsdf[i] < minv) minv = sm->bsdf[i]; + + if (minv <= FTINY) + return .0; + for (i = n; i--; ) sm->bsdf[i] -= minv; @@ -507,7 +549,7 @@ 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; /* adjust minimum hemispherical */ + df->maxHemi -= dv->cieY; /* adjust maximum hemispherical */ /* make sure everything is set */ c_ccvt(&dv->spec, C_CSXY+C_CSSPEC); } @@ -516,12 +558,11 @@ extract_diffuse(SDValue *dv, SDSpectralDF *df) 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) { @@ -540,12 +581,14 @@ 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")), @@ -556,21 +599,25 @@ SDloadMtx(SDData *sd, ezxml_t wtl) if ((rval = load_bsdf_data(sd, wdb, rowIn)) < 0) return convert_errcode(rval); } - /* separate diffuse components */ + /* separate diffuse components */ extract_diffuse(&sd->rLambFront, sd->rf); extract_diffuse(&sd->rLambBack, sd->rb); - extract_diffuse(&sd->tLamb, sd->tf); - /* return success */ + extract_diffuse(&sd->tLamb, (sd->tf != NULL) ? sd->tf : sd->tb); + /* return success */ return SDEnone; } /* 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); @@ -585,23 +632,28 @@ SDgetMtxBSDF(float coef[SDmaxCh], const FVECT outVec, return 1; /* XXX monochrome for now */ } -/* 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; + const SDMat *dp; double inc_psa, out_psa; /* check arguments */ - if ((psa == NULL) | (vec == NULL) | (dp == NULL)) + if ((psa == NULL) | (v1 == NULL) | (sdc == NULL) || + (dp = (SDMat *)sdc->dist) == NULL) return SDEargument; + if (v2 == NULL) + v2 = v1; /* get projected solid angles */ - inc_psa = mBSDF_incohm(dp, mBSDF_incndx(dp, vec)); - out_psa = mBSDF_outohm(dp, mBSDF_outndx(dp, vec)); + 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 SDqueryVal: - psa[0] = .0; - /* fall through */ case SDqueryMax: if (inc_psa > psa[0]) psa[0] = inc_psa; @@ -609,20 +661,24 @@ SDqueryMtxProjSA(double *psa, const FVECT vec, int qfl psa[0] = out_psa; break; case SDqueryMin+SDqueryMax: - if (inc_psa > psa[0]) + if (inc_psa > psa[1]) psa[1] = inc_psa; - if (out_psa > psa[0]) + 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])) + if ((inc_psa > 0) & (inc_psa < psa[0])) psa[0] = inc_psa; - if ((out_psa > .0) & (out_psa < psa[0])) + if ((out_psa > 0) & (out_psa < psa[0])) psa[0] = out_psa; break; } /* make sure it's legal */ - return (psa[0] <= .0) ? SDEinternal : SDEnone; + return (psa[0] <= 0) ? SDEinternal : SDEnone; } /* Compute new cumulative distribution from BSDF */ @@ -660,12 +716,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; /* check arguments */ - if ((inVec == NULL) | (dp == NULL)) + if ((inVec == NULL) | (sdc == NULL) || + (dp = (SDMat *)sdc->dist) == NULL) return NULL; memset(&myCD, 0, sizeof(myCD)); myCD.indx = mBSDF_incndx(dp, inVec); @@ -684,17 +741,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 */ @@ -706,7 +761,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 */ @@ -714,19 +769,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 ((outVec == NULL) | (mcd == NULL)) + 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; @@ -734,9 +789,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; }