--- ray/src/common/bsdf.c 2011/04/11 03:47:46 2.20 +++ ray/src/common/bsdf.c 2012/06/23 16:47:39 2.41 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: bsdf.c,v 2.20 2011/04/11 03:47:46 greg Exp $"; +static const char RCSid[] = "$Id: bsdf.c,v 2.41 2012/06/23 16:47:39 greg Exp $"; #endif /* * bsdf.c @@ -10,9 +10,12 @@ static const char RCSid[] = "$Id: bsdf.c,v 2.20 2011/0 * */ +#define _USE_MATH_DEFINES #include #include +#include #include +#include #include "ezxml.h" #include "hilbert.h" #include "bsdf.h" @@ -45,10 +48,14 @@ int SDretainSet = SDretainNone; SDError SDreportEnglish(SDError ec, FILE *fp) { - if (fp == NULL) - return ec; if (!ec) return SDEnone; + if ((ec < SDEnone) | (ec > SDEunknown)) { + SDerrorDetail[0] = '\0'; + ec = SDEunknown; + } + if (fp == NULL) + return ec; fputs(SDerrorEnglish[ec], fp); if (SDerrorDetail[0]) { fputs(": ", fp); @@ -85,7 +92,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")); @@ -95,13 +113,19 @@ SDloadGeometry(SDData *sd, ezxml_t wdb) if ((geom = ezxml_child(wdb, "Thickness")) != NULL) 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 ((sd->dim[0] < 0) | (sd->dim[1] < 0) | (sd->dim[2] < 0)) { + 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, @@ -110,6 +134,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"); @@ -153,27 +179,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); @@ -219,6 +241,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) @@ -245,7 +290,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); } @@ -334,7 +380,7 @@ SDgetCache(const char *bname) sdl->next = SDcacheList; SDcacheList = sdl; - sdl->refcnt++; + sdl->refcnt = 1; return &sdl->bsdf; } @@ -378,7 +424,7 @@ SDfreeCache(const SDData *sd) for (sdl = SDcacheList; sdl != NULL; sdl = (sdLast=sdl)->next) if (&sdl->bsdf == sd) break; - if (sdl == NULL || --sdl->refcnt) + if (sdl == NULL || (sdl->refcnt -= (sdl->refcnt > 0))) return; /* missing or still in use */ /* keep unreferenced data? */ if (SDisLoaded(sd) && SDretainSet) { @@ -401,34 +447,35 @@ SDfreeCache(const SDData *sd) /* Sample an individual BSDF component */ SDError -SDsampComponent(SDValue *sv, FVECT outVec, const FVECT inVec, - double randX, SDComponent *sdc) +SDsampComponent(SDValue *sv, FVECT ioVec, double randX, SDComponent *sdc) { float coef[SDmaxCh]; SDError ec; + FVECT inVec; const SDCDst *cd; double d; int n; /* check arguments */ - if ((sv == NULL) | (outVec == NULL) | (inVec == NULL) | (sdc == NULL)) + if ((sv == NULL) | (ioVec == NULL) | (sdc == NULL)) return SDEargument; /* get cumulative distribution */ + VCOPY(inVec, ioVec); 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(outVec, 0, 3*sizeof(double)); + memset(ioVec, 0, 3*sizeof(double)); return SDEnone; } sv->cieY = cd->cTotal; /* compute sample direction */ - ec = (*sdc->func->sampCDist)(outVec, randX, cd); + ec = (*sdc->func->sampCDist)(ioVec, randX, cd); if (ec) return ec; /* get BSDF color */ - n = (*sdc->func->getBSDFs)(coef, outVec, inVec, sdc->dist); + n = (*sdc->func->getBSDFs)(coef, ioVec, inVec, sdc); if (n <= 0) { strcpy(SDerrorDetail, "BSDF sample value error"); return SDEinternal; @@ -453,7 +500,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; @@ -476,7 +531,7 @@ SDdiffuseSamp(FVECT outVec, int outFront, double randX SDmultiSamp(outVec, 2, randX); SDsquare2disk(outVec, outVec[0], outVec[1]); outVec[2] = 1. - outVec[0]*outVec[0] - outVec[1]*outVec[1]; - if (outVec[2] > .0) /* a bit of paranoia */ + if (outVec[2] > 0) /* a bit of paranoia */ outVec[2] = sqrt(outVec[2]); if (!outFront) /* going out back? */ outVec[2] = -outVec[2]; @@ -484,13 +539,14 @@ SDdiffuseSamp(FVECT outVec, int outFront, double randX /* Query projected solid angle coverage for non-diffuse BSDF direction */ SDError -SDsizeBSDF(double *projSA, const FVECT vec, int qflags, const SDData *sd) +SDsizeBSDF(double *projSA, const FVECT v1, const RREAL *v2, + int qflags, const SDData *sd) { - SDSpectralDF *rdf; + SDSpectralDF *rdf, *tdf; SDError ec; int i; /* check arguments */ - if ((projSA == NULL) | (vec == NULL) | (sd == NULL)) + if ((projSA == NULL) | (v1 == NULL) | (sd == NULL)) return SDEargument; /* initialize extrema */ switch (qflags) { @@ -506,20 +562,26 @@ SDsizeBSDF(double *projSA, const FVECT vec, int qflags case 0: return SDEargument; } - if (vec[2] > .0) /* front surface query? */ + if (v1[2] > 0) /* front surface query? */ rdf = sd->rf; else rdf = sd->rb; + tdf = 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, vec, qflags, - rdf->comp[i].dist); + ec = (*rdf->comp[i].func->queryProjSA)(projSA, v1, v2, + qflags, &rdf->comp[i]); if (ec) return ec; } - for (i = (sd->tf==NULL) ? 0 : sd->tf->ncomp; i--; ) { - ec = (*sd->tf->comp[i].func->queryProjSA)(projSA, vec, qflags, - sd->tf->comp[i].dist); + for (i = (tdf==NULL) ? 0 : tdf->ncomp; i--; ) { + ec = (*tdf->comp[i].func->queryProjSA)(projSA, v1, v2, + qflags, &tdf->comp[i]); if (ec) return ec; } @@ -543,8 +605,8 @@ SDevalBSDF(SDValue *sv, const FVECT outVec, const FVEC if ((sv == NULL) | (outVec == NULL) | (inVec == NULL) | (sd == NULL)) return SDEargument; /* whose side are we on? */ - inFront = (inVec[2] > .0); - outFront = (outVec[2] > .0); + inFront = (inVec[2] > 0); + outFront = (outVec[2] > 0); /* start with diffuse portion */ if (inFront & outFront) { *sv = sd->rLambFront; @@ -561,7 +623,7 @@ SDevalBSDF(SDValue *sv, const FVECT outVec, const FVEC i = (sdf != NULL) ? sdf->ncomp : 0; while (i-- > 0) { nch = (*sdf->comp[i].func->getBSDFs)(coef, outVec, inVec, - sdf->comp[i].dist); + &sdf->comp[i]); while (nch-- > 0) { c_cmix(&sv->spec, sv->cieY, &sv->spec, coef[nch], &sdf->comp[i].cspec[nch]); @@ -585,7 +647,7 @@ SDdirectHemi(const FVECT inVec, int sflags, const SDDa if ((inVec == NULL) | (sd == NULL)) return .0; /* gather diffuse components */ - if (inVec[2] > .0) { + if (inVec[2] > 0) { hsum = sd->rLambFront.cieY; rdf = sd->rf; } else /* !inFront */ { @@ -616,10 +678,10 @@ SDdirectHemi(const FVECT inVec, int sflags, const SDDa /* Sample BSDF direction based on the given random variable */ SDError -SDsampBSDF(SDValue *sv, FVECT outVec, const FVECT inVec, - double randX, int sflags, const SDData *sd) +SDsampBSDF(SDValue *sv, FVECT ioVec, double randX, int sflags, const SDData *sd) { SDError ec; + FVECT inVec; int inFront; SDSpectralDF *rdf; double rdiff; @@ -628,11 +690,12 @@ SDsampBSDF(SDValue *sv, FVECT outVec, const FVECT inVe SDComponent *sdc; const SDCDst **cdarr = NULL; /* check arguments */ - if ((sv == NULL) | (outVec == NULL) | (inVec == NULL) | (sd == NULL) | - (randX < .0) | (randX >= 1.)) + if ((sv == NULL) | (ioVec == NULL) | (sd == NULL) | + (randX < 0) | (randX >= 1.)) return SDEargument; /* whose side are we on? */ - inFront = (inVec[2] > .0); + VCOPY(inVec, ioVec); + inFront = (inVec[2] > 0); /* remember diffuse portions */ if (inFront) { *sv = sd->rLambFront; @@ -670,16 +733,16 @@ SDsampBSDF(SDValue *sv, FVECT outVec, const FVECT inVe } 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(outVec, 0, 3*sizeof(double)); + memset(ioVec, 0, 3*sizeof(double)); return SDEnone; } /* scale random variable */ randX *= sv->cieY; /* diffuse reflection? */ if (randX < rdiff) { - SDdiffuseSamp(outVec, inFront, randX/rdiff); + SDdiffuseSamp(ioVec, inFront, randX/rdiff); goto done; } randX -= rdiff; @@ -687,7 +750,7 @@ SDsampBSDF(SDValue *sv, FVECT outVec, const FVECT inVe if ((sflags & SDsampDf+SDsampT) == SDsampDf+SDsampT) { if (randX < sd->tLamb.cieY) { sv->spec = sd->tLamb.spec; - SDdiffuseSamp(outVec, !inFront, randX/sd->tLamb.cieY); + SDdiffuseSamp(ioVec, !inFront, randX/sd->tLamb.cieY); goto done; } randX -= sd->tLamb.cieY; @@ -699,11 +762,11 @@ SDsampBSDF(SDValue *sv, FVECT outVec, const FVECT inVe return SDEinternal; /* compute sample direction */ sdc = (i < nr) ? &rdf->comp[i] : &sd->tf->comp[i-nr]; - ec = (*sdc->func->sampCDist)(outVec, randX/cdarr[i]->cTotal, cdarr[i]); + ec = (*sdc->func->sampCDist)(ioVec, randX/cdarr[i]->cTotal, cdarr[i]); if (ec) return ec; /* compute color */ - j = (*sdc->func->getBSDFs)(coef, outVec, inVec, sdc->dist); + j = (*sdc->func->getBSDFs)(coef, ioVec, inVec, sdc); if (j <= 0) { sprintf(SDerrorDetail, "BSDF \"%s\" sampling value error", sd->name); @@ -730,10 +793,10 @@ SDcompXform(RREAL vMtx[3][3], const FVECT sNrm, const if ((vMtx == NULL) | (sNrm == NULL) | (uVec == NULL)) return SDEargument; VCOPY(vMtx[2], sNrm); - if (normalize(vMtx[2]) == .0) + if (normalize(vMtx[2]) == 0) return SDEargument; fcross(vMtx[0], uVec, vMtx[2]); - if (normalize(vMtx[0]) == .0) + if (normalize(vMtx[0]) == 0) return SDEargument; fcross(vMtx[1], vMtx[2], vMtx[0]); return SDEnone; @@ -753,7 +816,7 @@ SDinvXform(RREAL iMtx[3][3], RREAL vMtx[3][3]) mTmp[0][1] = vMtx[2][1]*vMtx[0][2] - vMtx[2][2]*vMtx[0][1]; mTmp[0][2] = vMtx[1][2]*vMtx[0][1] - vMtx[1][1]*vMtx[0][2]; d = vMtx[0][0]*mTmp[0][0] + vMtx[1][0]*mTmp[0][1] + vMtx[2][0]*mTmp[0][2]; - if (d == .0) { + if (d == 0) { strcpy(SDerrorDetail, "Zero determinant in matrix inversion"); return SDEargument; } @@ -780,12 +843,12 @@ SDmapDir(FVECT resVec, RREAL vMtx[3][3], const FVECT i if (vMtx == NULL) { /* assume they just want to normalize */ if (resVec != inpVec) VCOPY(resVec, inpVec); - return (normalize(resVec) > .0) ? SDEnone : SDEargument; + return (normalize(resVec) > 0) ? SDEnone : SDEargument; } vTmp[0] = DOT(vMtx[0], inpVec); vTmp[1] = DOT(vMtx[1], inpVec); vTmp[2] = DOT(vMtx[2], inpVec); - if (normalize(vTmp) == .0) + if (normalize(vTmp) == 0) return SDEargument; VCOPY(resVec, vTmp); return SDEnone; @@ -800,7 +863,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 */ @@ -857,7 +919,7 @@ 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)); } @@ -920,7 +982,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); @@ -1160,7 +1222,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; @@ -1175,7 +1237,7 @@ check_bsdf_data( /* check that BSDF data is sane */ hemi_total = .0; for (i = dp->ninc; i--; ) { dom = getBSDF_incohm(dp,i); - if (dom <= .0) { + if (dom <= 0) { error(WARNING, "zero/negative incoming solid angle"); continue; } @@ -1198,7 +1260,7 @@ check_bsdf_data( /* check that BSDF data is sane */ hemi_total = .0; for (o = dp->nout; o--; ) { dom = getBSDF_outohm(dp,o); - if (dom <= .0) { + if (dom <= 0) { error(WARNING, "zero/negative outgoing solid angle"); continue; } @@ -1222,7 +1284,7 @@ check_bsdf_data( /* check that BSDF data is sane */ hemi_total = .0; for (o = dp->nout; o--; ) { double f = BSDF_value(dp,i,o); - if (f >= .0) + if (f >= 0) hemi_total += f*omega_oarr[o]; else { nneg += (f < -FTINY); @@ -1310,9 +1372,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"); @@ -1373,7 +1437,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) @@ -1398,7 +1462,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)