--- ray/src/common/bsdf.c 2011/04/24 20:16:52 2.25 +++ ray/src/common/bsdf.c 2012/09/02 15:33:15 2.42 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: bsdf.c,v 2.25 2011/04/24 20:16:52 greg Exp $"; +static const char RCSid[] = "$Id: bsdf.c,v 2.42 2012/09/02 15:33:15 greg Exp $"; #endif /* * bsdf.c @@ -10,9 +10,12 @@ static const char RCSid[] = "$Id: bsdf.c,v 2.25 2011/0 * */ +#define _USE_MATH_DEFINES #include #include +#include #include +#include #include "ezxml.h" #include "hilbert.h" #include "bsdf.h" @@ -35,6 +38,9 @@ const char *SDerrorEnglish[] = { /* Additional information on last error (ASCII English) */ char SDerrorDetail[256]; +/* Empty distribution for getCDist() calls that fail for some reason */ +const SDCDst SDemptyCD; + /* Cache of loaded BSDFs */ struct SDCache_s *SDcacheList = NULL; @@ -89,7 +95,18 @@ SDloadGeometry(SDData *sd, ezxml_t wdb) if (wdb == NULL) /* no geometry section? */ return SDEnone; + if ((geom = ezxml_child(wdb, "Name")) != NULL) { + strncpy(sd->matn, ezxml_txt(geom), SDnameLn); + if (sd->matn[SDnameLn-1]) + strcpy(sd->matn+(SDnameLn-4), "..."); + } + if ((geom = ezxml_child(wdb, "Manufacturer")) != NULL) { + strncpy(sd->makr, ezxml_txt(geom), SDnameLn); + if (sd->makr[SDnameLn-1]) + strcpy(sd->makr+(SDnameLn-4), "..."); + } sd->dim[0] = sd->dim[1] = sd->dim[2] = .0; + SDerrorDetail[0] = '\0'; if ((geom = ezxml_child(wdb, "Width")) != NULL) sd->dim[0] = atof(ezxml_txt(geom)) * to_meters(ezxml_attr(geom, "unit")); @@ -100,12 +117,18 @@ SDloadGeometry(SDData *sd, ezxml_t wdb) sd->dim[2] = atof(ezxml_txt(geom)) * to_meters(ezxml_attr(geom, "unit")); if ((sd->dim[0] < 0) | (sd->dim[1] < 0) | (sd->dim[2] < 0)) { - sprintf(SDerrorDetail, "Negative size in \"%s\"", sd->name); + if (!SDerrorDetail[0]) + sprintf(SDerrorDetail, "Negative dimension in \"%s\"", + sd->name); return SDEdata; } if ((geom = ezxml_child(wdb, "Geometry")) == NULL || (mgfstr = ezxml_txt(geom)) == NULL) return SDEnone; + while (isspace(*mgfstr)) + ++mgfstr; + if (!*mgfstr) + return SDEnone; if ((fmt = ezxml_attr(geom, "format")) != NULL && strcasecmp(fmt, "MGF")) { sprintf(SDerrorDetail, @@ -114,6 +137,8 @@ SDloadGeometry(SDData *sd, ezxml_t wdb) return SDEsupport; } cfact = to_meters(ezxml_attr(geom, "unit")); + if (cfact <= 0) + return SDEformat; sd->mgf = (char *)malloc(strlen(mgfstr)+32); if (sd->mgf == NULL) { strcpy(SDerrorDetail, "Out of memory in SDloadGeometry"); @@ -157,27 +182,23 @@ SDloadFile(SDData *sd, const char *fname) } wtl = ezxml_child(ezxml_child(fl, "Optical"), "Layer"); if (wtl == NULL) { - sprintf(SDerrorDetail, "BSDF \"%s\": no optical layer'", + sprintf(SDerrorDetail, "BSDF \"%s\": no optical layers'", sd->name); ezxml_free(fl); return SDEformat; } /* load geometry if present */ lastErr = SDloadGeometry(sd, ezxml_child(wtl, "Material")); - if (lastErr) + if (lastErr) { + ezxml_free(fl); return lastErr; + } /* try loading variable resolution data */ lastErr = SDloadTre(sd, wtl); /* check our result */ - switch (lastErr) { - case SDEformat: - case SDEdata: - case SDEsupport: /* possibly we just tried the wrong format */ + if (lastErr == SDEsupport) /* try matrix BSDF if not tree data */ lastErr = SDloadMtx(sd, wtl); - break; - default: /* variable res. OK else serious error */ - break; - } + /* done with XML file */ ezxml_free(fl); @@ -195,6 +216,9 @@ SDloadFile(SDData *sd, const char *fname) if (sd->tf != NULL && sd->tf->maxHemi <= .001) { SDfreeSpectralDF(sd->tf); sd->tf = NULL; } + if (sd->tb != NULL && sd->tb->maxHemi <= .001) { + SDfreeSpectralDF(sd->tb); sd->tb = NULL; + } /* return success */ return SDEnone; } @@ -223,6 +247,29 @@ SDnewSpectralDF(int nc) return df; } +/* Add component(s) to spectral distribution function */ +SDSpectralDF * +SDaddComponent(SDSpectralDF *odf, int nadd) +{ + SDSpectralDF *df; + + if (odf == NULL) + return SDnewSpectralDF(nadd); + if (nadd <= 0) + return odf; + df = (SDSpectralDF *)realloc(odf, sizeof(SDSpectralDF) + + (odf->ncomp+nadd-1)*sizeof(SDComponent)); + if (df == NULL) { + sprintf(SDerrorDetail, + "Cannot add %d component(s) to spectral DF", nadd); + SDfreeSpectralDF(odf); + return NULL; + } + memset(df->comp+df->ncomp, 0, nadd*sizeof(SDComponent)); + df->ncomp += nadd; + return df; +} + /* Free cached cumulative distributions for BSDF component */ void SDfreeCumulativeCache(SDSpectralDF *df) @@ -249,7 +296,8 @@ SDfreeSpectralDF(SDSpectralDF *df) return; SDfreeCumulativeCache(df); for (n = df->ncomp; n-- > 0; ) - (*df->comp[n].func->freeSC)(df->comp[n].dist); + if (df->comp[n].dist != NULL) + (*df->comp[n].func->freeSC)(df->comp[n].dist); free(df); } @@ -305,6 +353,10 @@ SDfreeBSDF(SDData *sd) SDfreeSpectralDF(sd->tf); sd->tf = NULL; } + if (sd->tb != NULL) { + SDfreeSpectralDF(sd->tb); + sd->tb = NULL; + } sd->rLambFront.cieY = .0; sd->rLambFront.spec.flags = 0; sd->rLambBack.cieY = .0; @@ -392,6 +444,7 @@ SDfreeCache(const SDData *sd) SDfreeCumulativeCache(sd->rf); SDfreeCumulativeCache(sd->rb); SDfreeCumulativeCache(sd->tf); + SDfreeCumulativeCache(sd->tb); return; } /* remove from list and free */ @@ -421,7 +474,7 @@ SDsampComponent(SDValue *sv, FVECT ioVec, double randX cd = (*sdc->func->getCDist)(inVec, sdc); if (cd == NULL) return SDEmemory; - if (cd->cTotal <= 1e-7) { /* anything to sample? */ + if (cd->cTotal <= 1e-6) { /* anything to sample? */ sv->spec = c_dfcolor; sv->cieY = .0; memset(ioVec, 0, 3*sizeof(double)); @@ -458,7 +511,15 @@ SDmultiSamp(double t[], int n, double randX) unsigned nBits; double scale; bitmask_t ndx, coord[MS_MAXDIM]; - + + if (n <= 0) /* check corner cases */ + return; + if (randX < 0) randX = 0; + else if (randX >= 1.) randX = 0.999999999999999; + if (n == 1) { + t[0] = randX; + return; + } while (n > MS_MAXDIM) /* punt for higher dimensions */ t[--n] = rand()*(1./(RAND_MAX+.5)); nBits = (8*sizeof(bitmask_t) - 1) / n; @@ -496,7 +557,7 @@ SDsizeBSDF(double *projSA, const FVECT v1, const RREAL SDError ec; int i; /* check arguments */ - if ((projSA == NULL) | (v1 == NULL)) + if ((projSA == NULL) | (v1 == NULL) | (sd == NULL)) return SDEargument; /* initialize extrema */ switch (qflags) { @@ -512,15 +573,18 @@ SDsizeBSDF(double *projSA, const FVECT v1, const RREAL case 0: return SDEargument; } - if (v1[2] > 0) /* front surface query? */ + if (v1[2] > 0) { /* front surface query? */ rdf = sd->rf; - else + tdf = (sd->tf != NULL) ? sd->tf : sd->tb; + } else { rdf = sd->rb; - tdf = NULL; /* transmitted component? */ - if (v2 != NULL && v1[2] > 0 ^ v2[2] > 0) { - rdf = NULL; - tdf = sd->tf; + tdf = (sd->tb != NULL) ? sd->tb : sd->tf; } + if (v2 != NULL) /* bidirectional? */ + if (v1[2] > 0 ^ v2[2] > 0) + rdf = NULL; + else + tdf = NULL; ec = SDEdata; /* run through components */ for (i = (rdf==NULL) ? 0 : rdf->ncomp; i--; ) { ec = (*rdf->comp[i].func->queryProjSA)(projSA, v1, v2, @@ -563,9 +627,12 @@ SDevalBSDF(SDValue *sv, const FVECT outVec, const FVEC } else if (!(inFront | outFront)) { *sv = sd->rLambBack; sdf = sd->rb; - } else /* inFront ^ outFront */ { + } else if (inFront) { *sv = sd->tLamb; - sdf = sd->tf; + sdf = (sd->tf != NULL) ? sd->tf : sd->tb; + } else /* inBack */ { + *sv = sd->tLamb; + sdf = (sd->tb != NULL) ? sd->tb : sd->tf; } sv->cieY *= 1./M_PI; /* add non-diffuse components */ @@ -589,7 +656,7 @@ double SDdirectHemi(const FVECT inVec, int sflags, const SDData *sd) { double hsum; - SDSpectralDF *rdf; + SDSpectralDF *rdf, *tdf; const SDCDst *cd; int i; /* check arguments */ @@ -599,26 +666,28 @@ SDdirectHemi(const FVECT inVec, int sflags, const SDDa if (inVec[2] > 0) { hsum = sd->rLambFront.cieY; rdf = sd->rf; + tdf = (sd->tf != NULL) ? sd->tf : sd->tb; } else /* !inFront */ { hsum = sd->rLambBack.cieY; rdf = sd->rb; + tdf = (sd->tb != NULL) ? sd->tb : sd->tf; } if ((sflags & SDsampDf+SDsampR) != SDsampDf+SDsampR) hsum = .0; if ((sflags & SDsampDf+SDsampT) == SDsampDf+SDsampT) hsum += sd->tLamb.cieY; /* gather non-diffuse components */ - i = ((sflags & SDsampSp+SDsampR) == SDsampSp+SDsampR && - rdf != NULL) ? rdf->ncomp : 0; + i = (((sflags & SDsampSp+SDsampR) == SDsampSp+SDsampR) & + (rdf != NULL)) ? rdf->ncomp : 0; while (i-- > 0) { /* non-diffuse reflection */ cd = (*rdf->comp[i].func->getCDist)(inVec, &rdf->comp[i]); if (cd != NULL) hsum += cd->cTotal; } - i = ((sflags & SDsampSp+SDsampT) == SDsampSp+SDsampT && - sd->tf != NULL) ? sd->tf->ncomp : 0; + i = (((sflags & SDsampSp+SDsampT) == SDsampSp+SDsampT) & + (tdf != NULL)) ? tdf->ncomp : 0; while (i-- > 0) { /* non-diffuse transmission */ - cd = (*sd->tf->comp[i].func->getCDist)(inVec, &sd->tf->comp[i]); + cd = (*tdf->comp[i].func->getCDist)(inVec, &tdf->comp[i]); if (cd != NULL) hsum += cd->cTotal; } @@ -632,7 +701,7 @@ SDsampBSDF(SDValue *sv, FVECT ioVec, double randX, int SDError ec; FVECT inVec; int inFront; - SDSpectralDF *rdf; + SDSpectralDF *rdf, *tdf; double rdiff; float coef[SDmaxCh]; int i, j, n, nr; @@ -649,9 +718,11 @@ SDsampBSDF(SDValue *sv, FVECT ioVec, double randX, int if (inFront) { *sv = sd->rLambFront; rdf = sd->rf; + tdf = (sd->tf != NULL) ? sd->tf : sd->tb; } else /* !inFront */ { *sv = sd->rLambBack; rdf = sd->rb; + tdf = (sd->tb != NULL) ? sd->tb : sd->tf; } if ((sflags & SDsampDf+SDsampR) != SDsampDf+SDsampR) sv->cieY = .0; @@ -659,15 +730,15 @@ SDsampBSDF(SDValue *sv, FVECT ioVec, double randX, int if ((sflags & SDsampDf+SDsampT) == SDsampDf+SDsampT) sv->cieY += sd->tLamb.cieY; /* gather non-diffuse components */ - i = nr = ((sflags & SDsampSp+SDsampR) == SDsampSp+SDsampR && - rdf != NULL) ? rdf->ncomp : 0; - j = ((sflags & SDsampSp+SDsampT) == SDsampSp+SDsampT && - sd->tf != NULL) ? sd->tf->ncomp : 0; + i = nr = (((sflags & SDsampSp+SDsampR) == SDsampSp+SDsampR) & + (rdf != NULL)) ? rdf->ncomp : 0; + j = (((sflags & SDsampSp+SDsampT) == SDsampSp+SDsampT) & + (tdf != NULL)) ? tdf->ncomp : 0; n = i + j; if (n > 0 && (cdarr = (const SDCDst **)malloc(n*sizeof(SDCDst *))) == NULL) return SDEmemory; while (j-- > 0) { /* non-diffuse transmission */ - cdarr[i+j] = (*sd->tf->comp[j].func->getCDist)(inVec, &sd->tf->comp[j]); + cdarr[i+j] = (*tdf->comp[j].func->getCDist)(inVec, &tdf->comp[j]); if (cdarr[i+j] == NULL) { free(cdarr); return SDEmemory; @@ -682,7 +753,7 @@ SDsampBSDF(SDValue *sv, FVECT ioVec, double randX, int } sv->cieY += cdarr[i]->cTotal; } - if (sv->cieY <= 1e-7) { /* anything to sample? */ + if (sv->cieY <= 1e-6) { /* anything to sample? */ sv->cieY = .0; memset(ioVec, 0, 3*sizeof(double)); return SDEnone; @@ -710,7 +781,7 @@ SDsampBSDF(SDValue *sv, FVECT ioVec, double randX, int if (i >= n) return SDEinternal; /* compute sample direction */ - sdc = (i < nr) ? &rdf->comp[i] : &sd->tf->comp[i-nr]; + sdc = (i < nr) ? &rdf->comp[i] : &tdf->comp[i-nr]; ec = (*sdc->func->sampCDist)(ioVec, randX/cdarr[i]->cTotal, cdarr[i]); if (ec) return ec; @@ -812,7 +883,6 @@ SDmapDir(FVECT resVec, RREAL vMtx[3][3], const FVECT i #include "standard.h" #include "paths.h" -#include #define MAXLATS 46 /* maximum number of latitudes */ @@ -932,7 +1002,7 @@ ab_getndx( /* get index corresponding to the given ve { ANGLE_BASIS *ab = (ANGLE_BASIS *)p; int li, ndx; - double pol, azi, d; + double pol, azi; if ((v[2] < -1.0) | (v[2] > 1.0)) return(-1); @@ -1172,7 +1242,7 @@ check_bsdf_data( /* check that BSDF data is sane */ ) { double *omega_iarr, *omega_oarr; - double dom, contrib, hemi_total, full_total; + double dom, hemi_total, full_total; int nneg; FVECT v; int i, o; @@ -1322,9 +1392,11 @@ load_BSDF( /* load BSDF data from file */ error(WARNING, errmsg); ezxml_free(fl); return(NULL); - } - load_angle_basis(ezxml_child(ezxml_child(wtl, - "DataDefinition"), "AngleBasis")); + } + for (wld = ezxml_child(ezxml_child(wtl, + "DataDefinition"), "AngleBasis"); + wld != NULL; wld = wld->next) + load_angle_basis(wld); dp = (struct BSDF_data *)calloc(1, sizeof(struct BSDF_data)); load_geometry(dp, ezxml_child(wtl, "Material")); for (wld = ezxml_child(wtl, "WavelengthData"); @@ -1385,7 +1457,7 @@ r_BSDF_incvec( /* compute random input vector at give if (!getBSDF_incvec(v, b, i)) return(0); rad = sqrt(getBSDF_incohm(b, i) / PI); - multisamp(pert, 3, rv); + SDmultiSamp(pert, 3, rv); for (j = 0; j < 3; j++) v[j] += rad*(2.*pert[j] - 1.); if (xm != NULL) @@ -1410,7 +1482,7 @@ r_BSDF_outvec( /* compute random output vector at giv if (!getBSDF_outvec(v, b, o)) return(0); rad = sqrt(getBSDF_outohm(b, o) / PI); - multisamp(pert, 3, rv); + SDmultiSamp(pert, 3, rv); for (j = 0; j < 3; j++) v[j] += rad*(2.*pert[j] - 1.); if (xm != NULL)