--- ray/src/common/bsdf.c 2009/06/17 20:41:47 2.1 +++ ray/src/common/bsdf.c 2011/01/03 19:27:00 2.9 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: bsdf.c,v 2.1 2009/06/17 20:41:47 greg Exp $"; +static const char RCSid[] = "$Id: bsdf.c,v 2.9 2011/01/03 19:27:00 greg Exp $"; #endif /* * Routines for handling BSDF data @@ -23,7 +23,7 @@ typedef struct { } lat[MAXLATS+1]; /* latitudes */ } ANGLE_BASIS; -#define MAXABASES 3 /* limit on defined bases */ +#define MAXABASES 7 /* limit on defined bases */ static ANGLE_BASIS abase_list[MAXABASES] = { { @@ -61,8 +61,42 @@ static ANGLE_BASIS abase_list[MAXABASES] = { static int nabases = 3; /* current number of defined bases */ +#define FEQ(a,b) ((a)-(b) <= 1e-6 && (b)-(a) <= 1e-6) static int +fequal(double a, double b) +{ + 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 +#ifdef ezxml_txt +#undef ezxml_txt +static char * +ezxml_txt(ezxml_t xml) +{ + if (xml == NULL) + return(""); + return(xml->txt); +} +#endif + + +static int ab_getvec( /* get vector for this angle basis index */ FVECT v, int ndx, @@ -71,16 +105,16 @@ ab_getvec( /* get vector for this angle basis index * { ANGLE_BASIS *ab = (ANGLE_BASIS *)p; int li; - double alt, azi, d; + double pol, azi, d; if ((ndx < 0) | (ndx >= ab->nangles)) return(0); for (li = 0; ndx >= ab->lat[li].nphis; li++) ndx -= ab->lat[li].nphis; - alt = PI/180.*0.5*(ab->lat[li].tmin + ab->lat[li+1].tmin); + pol = PI/180.*0.5*(ab->lat[li].tmin + ab->lat[li+1].tmin); azi = 2.*PI*ndx/ab->lat[li].nphis; - v[2] = d = cos(alt); - d = sqrt(1. - d*d); /* sin(alt) */ + v[2] = d = cos(pol); + d = sqrt(1. - d*d); /* sin(pol) */ v[0] = cos(azi)*d; v[1] = sin(azi)*d; return(1); @@ -95,14 +129,14 @@ ab_getndx( /* get index corresponding to the given ve { ANGLE_BASIS *ab = (ANGLE_BASIS *)p; int li, ndx; - double alt, azi, d; + double pol, azi, d; if ((v[2] < -1.0) | (v[2] > 1.0)) return(-1); - alt = 180.0/PI*acos(v[2]); + pol = 180.0/PI*acos(v[2]); azi = 180.0/PI*atan2(v[1], v[0]); if (azi < 0.0) azi += 360.0; - for (li = 1; ab->lat[li].tmin <= alt; li++) + for (li = 1; ab->lat[li].tmin <= pol; li++) if (!ab->lat[li].nphis) return(-1); --li; @@ -174,6 +208,103 @@ ab_getndxR( /* get index corresponding to the reverse static void +load_angle_basis( /* load custom BSDF angle basis */ + ezxml_t wab +) +{ + char *abname = ezxml_txt(ezxml_child(wab, "AngleBasisName")); + ezxml_t wbb; + int i; + + if (!abname || !*abname) + return; + for (i = nabases; i--; ) + if (!strcmp(abname, abase_list[i].name)) + return; /* assume it's the same */ + if (nabases >= MAXABASES) + error(INTERNAL, "too many angle bases"); + strcpy(abase_list[nabases].name, abname); + abase_list[nabases].nangles = 0; + for (i = 0, wbb = ezxml_child(wab, "AngleBasisBlock"); + wbb != NULL; i++, wbb = wbb->next) { + if (i >= MAXLATS) + error(INTERNAL, "too many latitudes in custom basis"); + abase_list[nabases].lat[i+1].tmin = atof(ezxml_txt( + ezxml_child(ezxml_child(wbb, + "ThetaBounds"), "UpperTheta"))); + if (!i) + abase_list[nabases].lat[i].tmin = + -abase_list[nabases].lat[i+1].tmin; + else if (!fequal(atof(ezxml_txt(ezxml_child(ezxml_child(wbb, + "ThetaBounds"), "LowerTheta"))), + abase_list[nabases].lat[i].tmin)) + error(WARNING, "theta values disagree in custom basis"); + abase_list[nabases].nangles += + abase_list[nabases].lat[i].nphis = + atoi(ezxml_txt(ezxml_child(wbb, "nPhis"))); + } + abase_list[nabases++].lat[i].nphis = 0; +} + + +static double +to_meters( /* return factor to convert given unit to meters */ + const char *unit +) +{ + if (unit == NULL) return(1.); /* safe assumption? */ + if (!strcasecmp(unit, "Meter")) return(1.); + if (!strcasecmp(unit, "Foot")) return(.3048); + if (!strcasecmp(unit, "Inch")) return(.0254); + if (!strcasecmp(unit, "Centimeter")) return(.01); + if (!strcasecmp(unit, "Millimeter")) return(.001); + sprintf(errmsg, "unknown dimensional unit '%s'", unit); + error(USER, errmsg); +} + + +static void +load_geometry( /* load geometric dimensions and description (if any) */ + struct BSDF_data *dp, + ezxml_t wdb +) +{ + ezxml_t geom; + double cfact; + const char *fmt, *mgfstr; + + dp->dim[0] = dp->dim[1] = dp->dim[2] = 0; + dp->mgf = NULL; + if ((geom = ezxml_child(wdb, "Width")) != NULL) + dp->dim[0] = atof(ezxml_txt(geom)) * + to_meters(ezxml_attr(geom, "unit")); + if ((geom = ezxml_child(wdb, "Height")) != NULL) + dp->dim[1] = atof(ezxml_txt(geom)) * + to_meters(ezxml_attr(geom, "unit")); + if ((geom = ezxml_child(wdb, "Thickness")) != NULL) + dp->dim[2] = atof(ezxml_txt(geom)) * + to_meters(ezxml_attr(geom, "unit")); + if ((geom = ezxml_child(wdb, "Geometry")) == NULL || + (mgfstr = ezxml_txt(geom)) == NULL) + return; + if ((fmt = ezxml_attr(geom, "format")) != NULL && + strcasecmp(fmt, "MGF")) { + sprintf(errmsg, "unrecognized geometry format '%s'", fmt); + error(WARNING, errmsg); + return; + } + cfact = to_meters(ezxml_attr(geom, "unit")); + dp->mgf = (char *)malloc(strlen(mgfstr)+32); + if (dp->mgf == NULL) + error(SYSTEM, "out of memory in load_geometry"); + if (cfact < 0.99 || cfact > 1.01) + sprintf(dp->mgf, "xf -s %.5f\n%s\nxf\n", cfact, mgfstr); + else + strcpy(dp->mgf, mgfstr); +} + + +static void load_bsdf_data( /* load BSDF distribution for this wavelength */ struct BSDF_data *dp, ezxml_t wdb @@ -184,11 +315,10 @@ load_bsdf_data( /* load BSDF distribution for this wa char *sdata; int i; - if ((cbasis == NULL) | (rbasis == NULL)) { + if ((!cbasis || !*cbasis) | (!rbasis || !*rbasis)) { error(WARNING, "missing column/row basis for BSDF"); return; } - /* XXX need to add routines for loading in foreign bases */ for (i = nabases; i--; ) if (!strcmp(cbasis, abase_list[i].name)) { dp->ninc = abase_list[i].nangles; @@ -199,7 +329,7 @@ load_bsdf_data( /* load BSDF distribution for this wa break; } if (i < 0) { - sprintf(errmsg, "unsupported ColumnAngleBasis '%s'", cbasis); + sprintf(errmsg, "undefined ColumnAngleBasis '%s'", cbasis); error(WARNING, errmsg); return; } @@ -213,13 +343,13 @@ load_bsdf_data( /* load BSDF distribution for this wa break; } if (i < 0) { - sprintf(errmsg, "unsupported RowAngleBasis '%s'", cbasis); + sprintf(errmsg, "undefined RowAngleBasis '%s'", cbasis); error(WARNING, errmsg); return; } /* read BSDF data */ sdata = ezxml_txt(ezxml_child(wdb,"ScatteringData")); - if (sdata == NULL) { + if (!sdata || !*sdata) { error(WARNING, "missing BSDF ScatteringData"); return; } @@ -243,7 +373,7 @@ load_bsdf_data( /* load BSDF distribution for this wa sdata++; if (*sdata) { sprintf(errmsg, "%d extra characters after BSDF ScatteringData", - strlen(sdata)); + (int)strlen(sdata)); error(WARNING, errmsg); } } @@ -254,64 +384,112 @@ check_bsdf_data( /* check that BSDF data is sane */ struct BSDF_data *dp ) { - double * omega_arr; - double dom, hemi_total; + double *omega_iarr, *omega_oarr; + double dom, contrib, hemi_total, full_total; int nneg; + FVECT v; int i, o; if (dp == NULL || dp->bsdf == NULL) return(0); - omega_arr = (double *)calloc(dp->nout, sizeof(double)); - if (omega_arr == NULL) + omega_iarr = (double *)calloc(dp->ninc, sizeof(double)); + omega_oarr = (double *)calloc(dp->nout, sizeof(double)); + if ((omega_iarr == NULL) | (omega_oarr == NULL)) error(SYSTEM, "out of memory in check_bsdf_data"); + /* incoming projected solid angles */ hemi_total = .0; + for (i = dp->ninc; i--; ) { + dom = getBSDF_incohm(dp,i); + if (dom <= .0) { + error(WARNING, "zero/negative incoming solid angle"); + continue; + } + if (!getBSDF_incvec(v,dp,i) || v[2] > FTINY) { + error(WARNING, "illegal incoming BSDF direction"); + free(omega_iarr); free(omega_oarr); + return(0); + } + hemi_total += omega_iarr[i] = dom * -v[2]; + } + if ((hemi_total > 1.02*PI) | (hemi_total < 0.98*PI)) { + sprintf(errmsg, "incoming BSDF hemisphere off by %.1f%%", + 100.*(hemi_total/PI - 1.)); + error(WARNING, errmsg); + } + dom = PI / hemi_total; /* fix normalization */ + for (i = dp->ninc; i--; ) + omega_iarr[i] *= dom; + /* outgoing projected solid angles */ + hemi_total = .0; for (o = dp->nout; o--; ) { - FVECT v; dom = getBSDF_outohm(dp,o); if (dom <= .0) { - error(WARNING, "zero/negative solid angle"); + error(WARNING, "zero/negative outgoing solid angle"); continue; } if (!getBSDF_outvec(v,dp,o) || v[2] < -FTINY) { error(WARNING, "illegal outgoing BSDF direction"); - free(omega_arr); + free(omega_iarr); free(omega_oarr); return(0); } - hemi_total += omega_arr[o] = dom*v[2]; + hemi_total += omega_oarr[o] = dom * v[2]; } if ((hemi_total > 1.02*PI) | (hemi_total < 0.98*PI)) { sprintf(errmsg, "outgoing BSDF hemisphere off by %.1f%%", 100.*(hemi_total/PI - 1.)); error(WARNING, errmsg); } - dom = PI / hemi_total; /* normalize solid angles */ + dom = PI / hemi_total; /* fix normalization */ for (o = dp->nout; o--; ) - omega_arr[o] *= dom; - nneg = 0; - for (i = dp->ninc; i--; ) { + omega_oarr[o] *= dom; + nneg = 0; /* check outgoing totals */ + for (i = 0; i < dp->ninc; i++) { hemi_total = .0; for (o = dp->nout; o--; ) { double f = BSDF_value(dp,i,o); - if (f > .0) - hemi_total += f*omega_arr[o]; - else if (f < -FTINY) - ++nneg; + if (f >= .0) + hemi_total += f*omega_oarr[o]; + else { + nneg += (f < -FTINY); + BSDF_value(dp,i,o) = .0f; + } } - if (hemi_total > 1.02) { - sprintf(errmsg, "BSDF direction passes %.1f%% of light", - 100.*hemi_total); + if (hemi_total > 1.01) { + sprintf(errmsg, + "incoming BSDF direction %d passes %.1f%% of light", + i, 100.*hemi_total); error(WARNING, errmsg); } } - free(omega_arr); - if (nneg > 0) { - sprintf(errmsg, "%d negative BSDF values", nneg); + if (nneg) { + sprintf(errmsg, "%d negative BSDF values (ignored)", nneg); error(WARNING, errmsg); - return(0); } + full_total = .0; /* reverse roles and check again */ + for (o = 0; o < dp->nout; o++) { + hemi_total = .0; + for (i = dp->ninc; i--; ) + hemi_total += BSDF_value(dp,i,o) * omega_iarr[i]; + + if (hemi_total > 1.01) { + sprintf(errmsg, + "outgoing BSDF direction %d collects %.1f%% of light", + o, 100.*hemi_total); + error(WARNING, errmsg); + } + full_total += hemi_total*omega_oarr[o]; + } + full_total /= PI; + if (full_total > 1.00001) { + sprintf(errmsg, "BSDF transfers %.4f%% of light", + 100.*full_total); + error(WARNING, errmsg); + } + free(omega_iarr); free(omega_oarr); return(1); } + struct BSDF_data * load_BSDF( /* load BSDF data from file */ char *fname @@ -348,7 +526,10 @@ load_BSDF( /* load BSDF data from file */ return(NULL); } wtl = ezxml_child(ezxml_child(fl, "Optical"), "Layer"); + load_angle_basis(ezxml_child(ezxml_child(wtl, + "DataDefinition"), "AngleBasis")); dp = (struct BSDF_data *)calloc(1, sizeof(struct BSDF_data)); + load_geometry(dp, ezxml_child(wtl, "Material")); for (wld = ezxml_child(wtl, "WavelengthData"); wld != NULL; wld = wld->next) { if (strcmp(ezxml_txt(ezxml_child(wld,"Wavelength")), "Visible")) @@ -379,6 +560,8 @@ free_BSDF( /* free BSDF data structure */ { if (b == NULL) return; + if (b->mgf != NULL) + free(b->mgf); if (b->bsdf != NULL) free(b->bsdf); free(b); @@ -435,8 +618,6 @@ r_BSDF_outvec( /* compute random output vector at giv } -#define FEQ(a,b) ((a)-(b) <= 1e-7 && (b)-(a) <= 1e-7) - static int addrot( /* compute rotation (x,y,z) => (xp,yp,zp) */ char *xfarg[], @@ -487,12 +668,14 @@ int getBSDF_xfm( /* compute BSDF orient. -> world orient. transform */ MAT4 xm, FVECT nrm, - UpDir ud + UpDir ud, + char *xfbuf ) { char *xfargs[7]; XF myxf; FVECT updir, xdest, ydest; + int i; updir[0] = updir[1] = updir[2] = 0.; switch (ud) { @@ -523,5 +706,13 @@ getBSDF_xfm( /* compute BSDF orient. -> world orient. fcross(ydest, nrm, xdest); xf(&myxf, addrot(xfargs, xdest, ydest, nrm), xfargs); copymat4(xm, myxf.xfm); + if (xfbuf == NULL) + return(1); + /* return xf arguments as well */ + for (i = 0; xfargs[i] != NULL; i++) { + *xfbuf++ = ' '; + strcpy(xfbuf, xfargs[i]); + while (*xfbuf) ++xfbuf; + } return(1); }