--- ray/src/common/bsdf.c 2011/04/24 20:16:52 2.25 +++ 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.25 2011/04/24 20:16:52 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.25 2011/0 * */ +#define _USE_MATH_DEFINES #include #include +#include #include +#include #include "ezxml.h" #include "hilbert.h" #include "bsdf.h" @@ -89,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")); @@ -100,12 +114,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 +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"); @@ -157,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); @@ -223,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) @@ -249,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); } @@ -421,7 +463,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 +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; @@ -496,7 +546,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) { @@ -516,11 +566,12 @@ SDsizeBSDF(double *projSA, const FVECT v1, const RREAL rdf = sd->rf; else rdf = sd->rb; - tdf = NULL; /* transmitted component? */ - if (v2 != NULL && v1[2] > 0 ^ v2[2] > 0) { - rdf = NULL; - tdf = sd->tf; - } + 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, v1, v2, @@ -682,7 +733,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; @@ -812,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 */ @@ -932,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); @@ -1172,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; @@ -1322,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"); @@ -1385,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) @@ -1410,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)