--- ray/src/cv/checkBSDF.c 2021/12/14 02:33:18 2.1 +++ ray/src/cv/checkBSDF.c 2023/09/11 18:43:36 2.10 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: checkBSDF.c,v 2.1 2021/12/14 02:33:18 greg Exp $"; +static const char RCSid[] = "$Id: checkBSDF.c,v 2.10 2023/09/11 18:43:36 greg Exp $"; #endif /* * checkBSDF.c @@ -10,6 +10,7 @@ static const char RCSid[] = "$Id: checkBSDF.c,v 2.1 20 #define _USE_MATH_DEFINES #include #include "rtio.h" +#include "random.h" #include "bsdf.h" #include "bsdf_m.h" #include "bsdf_t.h" @@ -19,6 +20,17 @@ static const char RCSid[] = "$Id: checkBSDF.c,v 2.1 20 #define F_MATRIX 0x4 #define F_TTREE 0x8 +typedef struct { + double vmin, vmax; /* extrema */ + double vsum; /* straight sum */ + long nvals; /* number of values */ +} SimpleStats; + +const SimpleStats SSinit = {FHUGE, -FHUGE, .0, 0}; + + /* relative difference formula */ +#define rdiff(a,b) ((a)>(b) ? ((a)-(b))/((a)+FTINY) : ((b)-(a))/((b)+FTINY)) + /* Figure out BSDF type (and optionally determine if in color) */ const char * getBSDFtype(const SDData *bsdf, int *flags) @@ -68,25 +80,100 @@ getBSDFtype(const SDData *bsdf, int *flags) void detailComponent(const char *nm, const SDValue *lamb, const SDSpectralDF *df) { - printf("%s\t%4.1f %4.1f %4.1f\t\t", nm, 100.*lamb->cieY*lamb->spec.cx/lamb->spec.cy, + double Lamb = 0; + + fputs(nm, stdout); + if (lamb->spec.flags) { + printf("\t%4.1f %4.1f %4.1f\t\t", + 100.*lamb->cieY*lamb->spec.cx/lamb->spec.cy, 100.*lamb->cieY, 100.*lamb->cieY*(1.f - lamb->spec.cx - lamb->spec.cy)/lamb->spec.cy); + Lamb = lamb->cieY; + } else + fputs("\t 0 0 0\t\t", stdout); if (df) - printf("%5.1f%%\t\t%.2f deg\n", 100.*df->maxHemi, + printf("%5.1f%%\t\t%.2f deg\n", 100.*(Lamb+df->maxHemi), sqrt(df->minProjSA/M_PI)*(360./M_PI)); else - puts("0%\t\t180"); + printf("%5.1f%%\t\t180 deg\n", Lamb); } +/* Add a value to stats */ +void +addStat(SimpleStats *ssp, double v) +{ + if (v < ssp->vmin) ssp->vmin = v; + if (v > ssp->vmax) ssp->vmax = v; + ssp->vsum += v; + ssp->nvals++; +} + +/* Sample a BSDF hemisphere with callback (quadtree recursion) */ +int +qtSampBSDF(double xleft, double ytop, double siz, + const SDData *bsdf, const int side, const RREAL *v0, + int (*cf)(const SDData *b, const FVECT v1, const RREAL *v0, void *p), + void *cdata) +{ + if (siz < 0.124) { /* make sure we subdivide, first */ + FVECT vsmp; + double sa; + square2disk(vsmp, xleft + frandom()*siz, ytop + frandom()*siz); + vsmp[2] = 1. - vsmp[0]*vsmp[0] - vsmp[1]*vsmp[1]; + if (vsmp[2] <= 0) return 0; + vsmp[2] = side * sqrt(vsmp[2]); + if (SDreportError( SDsizeBSDF(&sa, vsmp, v0, SDqueryMin, bsdf), stderr)) + return 0; + if (sa >= M_PI*siz*siz - FTINY) /* no further division needed */ + return (*cf)(bsdf, vsmp, v0, cdata); + } + siz *= .5; /* 4-branch recursion */ + return( qtSampBSDF(xleft, ytop, siz, bsdf, side, v0, cf, cdata) && + qtSampBSDF(xleft+siz, ytop, siz, bsdf, side, v0, cf, cdata) && + qtSampBSDF(xleft, ytop+siz, siz, bsdf, side, v0, cf, cdata) && + qtSampBSDF(xleft+siz, ytop+siz, siz, bsdf, side, v0, cf, cdata) ); +} + +#define sampBSDFhemi(b,s,v0,cf,cd) qtSampBSDF(0,0,1,b,s,v0,cf,cd) + +/* Call-back to compute reciprocity difference */ +int +diffRecip(const SDData *bsdf, const FVECT v1, const RREAL *v0, void *p) +{ + SDValue sdv; + double otherY; + + if (SDreportError( SDevalBSDF(&sdv, v0, v1, bsdf), stderr)) + return 0; + otherY = sdv.cieY; + if (SDreportError( SDevalBSDF(&sdv, v1, v0, bsdf), stderr)) + return 0; + + addStat((SimpleStats *)p, rdiff(sdv.cieY, otherY)); + return 1; +} + +/* Call-back to compute reciprocity over reflected hemisphere */ +int +reflHemi(const SDData *bsdf, const FVECT v1, const RREAL *v0, void *p) +{ + return sampBSDFhemi(bsdf, 1 - 2*(v1[2]<0), v1, &diffRecip, p); +} + +/* Call-back to compute reciprocity over transmitted hemisphere */ +int +transHemi(const SDData *bsdf, const FVECT v1, const RREAL *v0, void *p) +{ + return sampBSDFhemi(bsdf, 1 - 2*(v1[2]>0), v1, &diffRecip, p); +} + /* Report reciprocity errors for the given directions */ void checkReciprocity(const char *nm, const int side1, const int side2, const SDData *bsdf, const int fl) { + SimpleStats myStats = SSinit; const SDSpectralDF *df = bsdf->tf; - double emin=FHUGE, emax=0, esum=0; - int ntested=0; - int ec; if (side1 == side2) { df = (side1 > 0) ? bsdf->rf : bsdf->rb; @@ -97,9 +184,15 @@ checkReciprocity(const char *nm, const int side1, cons if (fl & F_MATRIX) { /* special case for matrix BSDF */ const SDMat *m = (const SDMat *)df->comp[0].dist; int i = m->ninc; + double diffuseY; FVECT vin, vout; - double rerr; - SDValue vrev; + double fwdY; + SDValue rev; + if (side1 == side2) + diffuseY = (side1 > 0) ? bsdf->rLambFront.cieY : bsdf->rLambBack.cieY; + else + diffuseY = (side1 > 0) ? bsdf->tLambFront.cieY : bsdf->tLambBack.cieY; + diffuseY /= M_PI; while (i--) { int o = m->nout; if (!mBSDF_incvec(vin, m, i+.5)) @@ -107,26 +200,36 @@ checkReciprocity(const char *nm, const int side1, cons while (o--) { if (!mBSDF_outvec(vout, m, o+.5)) continue; - rerr = mBSDF_value(m, o, i); - if (rerr <= FTINY) + fwdY = mBSDF_value(m, o, i) + diffuseY; + if (fwdY <= 1e-4) continue; - if (SDreportError( SDevalBSDF(&vrev, vout, vin, bsdf), stderr)) + if (SDreportError( SDevalBSDF(&rev, vout, vin, bsdf), stderr)) return; - rerr = 100.*fabs(rerr - vrev.cieY)/rerr; - if (rerr < emin) emin = rerr; - if (rerr > emax) emax = rerr; - esum += rerr; - ++ntested; + if (rev.cieY > 1e-4) + addStat(&myStats, rdiff(fwdY, rev.cieY)); } } - } else { - } - if (ntested) { - printf("%s\t%.1f\t%.1f\t%.1f\n", nm, emin, esum/(double)ntested, emax); + } else if (fl & F_ISOTROPIC) { /* isotropic case */ + const double stepSize = sqrt(df->minProjSA/M_PI); + FVECT vin; + vin[1] = 0; + for (vin[0] = 0.5*stepSize; vin[0] < 1; vin[0] += stepSize) { + vin[2] = side1*sqrt(1. - vin[0]*vin[0]); + if (!sampBSDFhemi(bsdf, side2, vin, &diffRecip, &myStats)) + return; + } + } else if (!sampBSDFhemi(bsdf, side1, NULL, + (side1==side2) ? &reflHemi : &transHemi, &myStats)) return; + if (myStats.nvals) { + printf("%s\t%5.1f\t%5.1f\t%5.1f\n", nm, + 100.*myStats.vmin, + 100.*myStats.vsum/(double)myStats.nvals, + 100.*myStats.vmax); + return; } nothing2do: - printf("%s\t0\t0\t0\n", nm); + printf("%s\t 0\t 0\t 0\n", nm); } /* Report on the given BSDF XML file */ @@ -134,19 +237,19 @@ int checkXML(char *fname) { int flags; - SDError ec; SDData myBSDF; char *pth; + puts("====================================================="); printf("File: '%s'\n", fname); SDclearBSDF(&myBSDF, fname); - pth = getpath(fname, getrlibpath(), R_OK); + pth = getpath(fname, getrlibpath(), 0); if (!pth) { fprintf(stderr, "Cannot find file '%s'\n", fname); return 0; } - ec = SDloadFile(&myBSDF, pth); - if (ec) goto err; + if (SDreportError( SDloadFile(&myBSDF, pth), stderr)) + return 0; printf("Manufacturer: '%s'\n", myBSDF.makr); printf("BSDF Name: '%s'\n", myBSDF.matn); printf("Dimensions (W x H x Thickness): %g x %g x %g cm\n", 100.*myBSDF.dim[0], @@ -154,21 +257,17 @@ checkXML(char *fname) printf("Type: %s\n", getBSDFtype(&myBSDF, &flags)); printf("Color: %d\n", (flags & F_IN_COLOR) != 0); printf("Has Geometry: %d\n", (myBSDF.mgf != NULL)); - puts("Component\tLambertian XYZ %\tMax. Dir\tMin. Angle"); - detailComponent("Internal Refl", &myBSDF.rLambFront, myBSDF.rf); - detailComponent("External Refl", &myBSDF.rLambBack, myBSDF.rb); + puts("Component\tLambertian XYZ (%)\tMax. Tot. Dir\tMin. Angle"); + detailComponent("Interior Refl", &myBSDF.rLambFront, myBSDF.rf); + detailComponent("Exterior Refl", &myBSDF.rLambBack, myBSDF.rb); detailComponent("Int->Ext Trans", &myBSDF.tLambFront, myBSDF.tf); detailComponent("Ext->Int Trans", &myBSDF.tLambBack, myBSDF.tb); - puts("Component\tReciprocity Error (min/avg/max %)"); - checkReciprocity("Front Refl", 1, 1, &myBSDF, flags); - checkReciprocity("Back Refl", -1, -1, &myBSDF, flags); + puts("Component\tReciprocity Error (min avg max %)"); + checkReciprocity("Interior Refl", 1, 1, &myBSDF, flags); + checkReciprocity("Exterior Refl", -1, -1, &myBSDF, flags); checkReciprocity("Transmission", -1, 1, &myBSDF, flags); SDfreeBSDF(&myBSDF); return 1; -err: - SDreportError(ec, stderr); - SDfreeBSDF(&myBSDF); - return 0; } int @@ -180,10 +279,8 @@ main(int argc, char *argv[]) fprintf(stderr, "Usage: %s bsdf.xml ..\n", argv[0]); return 1; } - for (i = 1; i < argc; i++) { - puts("====================================================="); + for (i = 1; i < argc; i++) if (!checkXML(argv[i])) return 1; - } return 0; }