--- ray/src/common/bsdf_m.c 2011/02/18 02:41:55 3.3 +++ ray/src/common/bsdf_m.c 2012/05/07 23:17:59 3.22 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: bsdf_m.c,v 3.3 2011/02/18 02:41:55 greg Exp $"; +static const char RCSid[] = "$Id: bsdf_m.c,v 3.22 2012/05/07 23:17:59 greg Exp $"; #endif /* * bsdf_m.c @@ -11,6 +11,7 @@ static const char RCSid[] = "$Id: bsdf_m.c,v 3.3 2011/ * */ +#define _USE_MATH_DEFINES #include "rtio.h" #include #include @@ -19,10 +20,6 @@ static const char RCSid[] = "$Id: bsdf_m.c,v 3.3 2011/ #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 @@ -49,7 +46,7 @@ typedef struct { static ANGLE_BASIS abase_list[MAXABASES] = { { "LBNL/Klems Full", 145, - { {-5., 1}, + { {0., 1}, {5., 8}, {15., 16}, {25., 20}, @@ -61,7 +58,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}, @@ -71,7 +68,7 @@ static ANGLE_BASIS abase_list[MAXABASES] = { {90., 0} } }, { "LBNL/Klems Quarter", 41, - { {-9., 1}, + { {0., 1}, {9., 8}, {27., 12}, {46., 12}, @@ -85,24 +82,12 @@ static int nabases = 3; /* current number of defined b 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 * @@ -161,16 +146,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; @@ -185,17 +172,17 @@ 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; + 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]); @@ -214,9 +201,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; @@ -231,21 +218,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 */ +/* 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]; @@ -255,9 +238,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; @@ -265,9 +248,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) @@ -299,13 +333,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 += @@ -315,7 +348,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; } } @@ -350,8 +383,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) @@ -364,41 +396,43 @@ get_extrema(SDSpectralDF *df) 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")) { + 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) 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) { @@ -407,11 +441,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")); @@ -421,41 +454,48 @@ 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; /* 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); @@ -469,12 +509,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); @@ -494,6 +534,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; @@ -519,30 +563,28 @@ 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 maximum 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; - 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")) @@ -553,12 +595,14 @@ 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")), @@ -567,37 +611,27 @@ 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 */ + /* separate diffuse components */ extract_diffuse(&sd->rLambFront, sd->rf); extract_diffuse(&sd->rLambBack, sd->rb); extract_diffuse(&sd->tLamb, sd->tf); - /* return success */ + /* 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); @@ -612,59 +646,53 @@ 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; - - 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 */ @@ -702,12 +730,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); @@ -726,17 +755,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 */ @@ -748,7 +775,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 */ @@ -756,16 +783,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; @@ -773,9 +803,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; }