src/cv/checkBSDF.c

#ifndef lint
static const char RCSid[] = "$Id: checkBSDF.c,v 2.5 2022/01/26 16:19:30 greg Exp $";
#endif
/*
 *  checkBSDF.c
 *
 *  Load BSDF XML file and check Helmholtz reciprocity
 */

#define _USE_MATH_DEFINES
#include <math.h>
#include "rtio.h"
#include "random.h"
#include "bsdf.h"
#include "bsdf_m.h"
#include "bsdf_t.h"

#define F_IN_COLOR      0x1
#define F_ISOTROPIC     0x2
#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)
{
        const SDSpectralDF      *df = bsdf->tb;

        if (flags) *flags = 0;
        if (!df) df = bsdf->tf;
        if (!df) df = bsdf->rf;
        if (!df) df = bsdf->rb;
        if (!df) return "Pure_Lambertian";
        if (df->comp[0].func == &SDhandleMtx) {
                const SDMat     *m = (const SDMat *)df->comp[0].dist;
                if (flags) {
                        *flags |= F_MATRIX;
                        *flags |= F_IN_COLOR*(m->chroma != NULL);
                }
                switch (m->ninc) {
                case 145:
                        return "Klems_Full";
                case  73:
                        return "Klems_Half";
                case  41:
                        return "Klems_Quarter";
                }
                return "Unknown_Matrix";
        }
        if (df->comp[0].func == &SDhandleTre) {
                const SDTre     *t = (const SDTre *)df->comp[0].dist;
                if (flags) {
                        *flags |= F_TTREE;
                        *flags |= F_IN_COLOR*(t->stc[1] != NULL);
                }
                switch (t->stc[0]->ndim) {
                case 4:
                        return "Anisotropic_Tensor_Tree";
                case 3:
                        if (flags) *flags |= F_ISOTROPIC;
                        return "Isotropic_Tensor_Tree";
                }
                return "Unknown_Tensor_Tree";
        }
        return "Unknown";
}

/* Report details related to one hemisphere distribution */
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,
                        100.*lamb->cieY,
                        100.*lamb->cieY*(1.f - lamb->spec.cx - lamb->spec.cy)/lamb->spec.cy);
        if (df)
                printf("%5.1f%%\t\t%.2f deg\n", 100.*df->maxHemi,
                                sqrt(df->minProjSA/M_PI)*(360./M_PI));
        else
                puts("0%\t\t180");
}

/* 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;

        if (side1 == side2) {
                df = (side1 > 0) ? bsdf->rf : bsdf->rb;
                if (!df) goto nothing2do;
        } else if (!bsdf->tf | !bsdf->tb)
                goto nothing2do;

        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          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))
                        continue;
                    while (o--) {
                        if (!mBSDF_outvec(vout, m, o+.5))
                                continue;
                        fwdY = mBSDF_value(m, o, i) + diffuseY;
                        if (fwdY <= 1e-4)
                                continue;
                        if (SDreportError( SDevalBSDF(&rev, vout, vin, bsdf), stderr))
                                return;
                        if (rev.cieY > 1e-4)
                                addStat(&myStats, rdiff(fwdY, rev.cieY));
                    }
                }
        } 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);
}

/* Report on the given BSDF XML file */
int
checkXML(char *fname)
{
        int             flags;
        SDData          myBSDF;
        char            *pth;

        puts("=====================================================");
        printf("File: '%s'\n", fname);
        SDclearBSDF(&myBSDF, fname);
        pth = getpath(fname, getrlibpath(), 0);
        if (!pth) {
                fprintf(stderr, "Cannot find file '%s'\n", fname);
                return 0;
        }
        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],
                                100.*myBSDF.dim[1], 100.*myBSDF.dim[2]);
        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("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("Interior Refl", 1, 1, &myBSDF, flags);
        checkReciprocity("Exterior Refl", -1, -1, &myBSDF, flags);
        checkReciprocity("Transmission", -1, 1, &myBSDF, flags);
        SDfreeBSDF(&myBSDF);
        return 1;
}

int
main(int argc, char *argv[])
{
        int     i;

        if (argc < 2) {
                fprintf(stderr, "Usage: %s bsdf.xml ..\n", argv[0]);
                return 1;
        }
        for (i = 1; i < argc; i++)
                if (!checkXML(argv[i]))
                        return 1;
        return 0;
}
Revision:	2.6
Committed:	Wed Jan 26 17:30:07 2022 UTC (3 years, 9 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.5:	+2 -2 lines
Log Message:	fix(checkBSDF): Fixed NaN error
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.6	static const char RCSid[] = "$Id: checkBSDF.c,v 2.5 2022/01/26 16:19:30 greg Exp $";
3	greg	2.1	#endif
4			/*
5			* checkBSDF.c
6			*
7			* Load BSDF XML file and check Helmholtz reciprocity
8			*/
9
10			#define _USE_MATH_DEFINES
11			#include <math.h>
12			#include "rtio.h"
13	greg	2.2	#include "random.h"
14	greg	2.1	#include "bsdf.h"
15			#include "bsdf_m.h"
16			#include "bsdf_t.h"
17
18			#define F_IN_COLOR 0x1
19			#define F_ISOTROPIC 0x2
20			#define F_MATRIX 0x4
21			#define F_TTREE 0x8
22
23	greg	2.2	typedef struct {
24			double vmin, vmax; /* extrema */
25			double vsum; /* straight sum */
26			long nvals; /* number of values */
27			} SimpleStats;
28
29			const SimpleStats SSinit = {FHUGE, -FHUGE, .0, 0};
30
31			/* relative difference formula */
32	greg	2.6	#define rdiff(a,b) ((a)>(b) ? ((a)-(b))/((a)+FTINY) : ((b)-(a))/((b)+FTINY))
33	greg	2.2
34	greg	2.1	/* Figure out BSDF type (and optionally determine if in color) */
35			const char *
36			getBSDFtype(const SDData bsdf, int flags)
37			{
38			const SDSpectralDF *df = bsdf->tb;
39
40			if (flags) *flags = 0;
41			if (!df) df = bsdf->tf;
42			if (!df) df = bsdf->rf;
43			if (!df) df = bsdf->rb;
44			if (!df) return "Pure_Lambertian";
45			if (df->comp[0].func == &SDhandleMtx) {
46			const SDMat m = (const SDMat )df->comp[0].dist;
47			if (flags) {
48			*flags \|= F_MATRIX;
49			flags \|= F_IN_COLOR(m->chroma != NULL);
50			}
51			switch (m->ninc) {
52			case 145:
53			return "Klems_Full";
54			case 73:
55			return "Klems_Half";
56			case 41:
57			return "Klems_Quarter";
58			}
59			return "Unknown_Matrix";
60			}
61			if (df->comp[0].func == &SDhandleTre) {
62			const SDTre t = (const SDTre )df->comp[0].dist;
63			if (flags) {
64			*flags \|= F_TTREE;
65			flags \|= F_IN_COLOR(t->stc[1] != NULL);
66			}
67			switch (t->stc[0]->ndim) {
68			case 4:
69			return "Anisotropic_Tensor_Tree";
70			case 3:
71			if (flags) *flags \|= F_ISOTROPIC;
72			return "Isotropic_Tensor_Tree";
73			}
74			return "Unknown_Tensor_Tree";
75			}
76			return "Unknown";
77			}
78
79			/* Report details related to one hemisphere distribution */
80			void
81			detailComponent(const char nm, const SDValue lamb, const SDSpectralDF *df)
82			{
83	greg	2.2	printf("%s\t%4.1f %4.1f %4.1f\t\t", nm,
84			100.lamb->cieYlamb->spec.cx/lamb->spec.cy,
85	greg	2.1	100.*lamb->cieY,
86			100.lamb->cieY(1.f - lamb->spec.cx - lamb->spec.cy)/lamb->spec.cy);
87			if (df)
88			printf("%5.1f%%\t\t%.2f deg\n", 100.*df->maxHemi,
89			sqrt(df->minProjSA/M_PI)*(360./M_PI));
90			else
91			puts("0%\t\t180");
92			}
93
94	greg	2.2	/* Add a value to stats */
95			void
96			addStat(SimpleStats *ssp, double v)
97			{
98			if (v < ssp->vmin) ssp->vmin = v;
99			if (v > ssp->vmax) ssp->vmax = v;
100			ssp->vsum += v;
101			ssp->nvals++;
102			}
103
104			/* Sample a BSDF hemisphere with callback (quadtree recursion) */
105			int
106			qtSampBSDF(double xleft, double ytop, double siz,
107			const SDData bsdf, const int side, const RREAL v0,
108			int (cf)(const SDData b, const FVECT v1, const RREAL v0, void p),
109			void *cdata)
110			{
111			if (siz < 0.124) { /* make sure we subdivide, first */
112			FVECT vsmp;
113			double sa;
114			square2disk(vsmp, xleft + frandom()siz, ytop + frandom()siz);
115			vsmp[2] = 1. - vsmp[0]vsmp[0] - vsmp[1]vsmp[1];
116			if (vsmp[2] <= 0) return 0;
117			vsmp[2] = side * sqrt(vsmp[2]);
118			if (SDreportError( SDsizeBSDF(&sa, vsmp, v0, SDqueryMin, bsdf), stderr))
119			return 0;
120			if (sa >= M_PIsizsiz - FTINY) /* no further division needed */
121			return (*cf)(bsdf, vsmp, v0, cdata);
122			}
123			siz = .5; / 4-branch recursion */
124			return( qtSampBSDF(xleft, ytop, siz, bsdf, side, v0, cf, cdata) &&
125			qtSampBSDF(xleft+siz, ytop, siz, bsdf, side, v0, cf, cdata) &&
126			qtSampBSDF(xleft, ytop+siz, siz, bsdf, side, v0, cf, cdata) &&
127			qtSampBSDF(xleft+siz, ytop+siz, siz, bsdf, side, v0, cf, cdata) );
128			}
129
130			#define sampBSDFhemi(b,s,v0,cf,cd) qtSampBSDF(0,0,1,b,s,v0,cf,cd)
131
132			/* Call-back to compute reciprocity difference */
133			int
134			diffRecip(const SDData bsdf, const FVECT v1, const RREAL v0, void *p)
135			{
136			SDValue sdv;
137			double otherY;
138
139			if (SDreportError( SDevalBSDF(&sdv, v0, v1, bsdf), stderr))
140			return 0;
141			otherY = sdv.cieY;
142			if (SDreportError( SDevalBSDF(&sdv, v1, v0, bsdf), stderr))
143			return 0;
144
145			addStat((SimpleStats *)p, rdiff(sdv.cieY, otherY));
146			return 1;
147			}
148
149			/* Call-back to compute reciprocity over reflected hemisphere */
150			int
151			reflHemi(const SDData bsdf, const FVECT v1, const RREAL v0, void *p)
152			{
153			return sampBSDFhemi(bsdf, 1 - 2*(v1[2]<0), v1, &diffRecip, p);
154			}
155
156			/* Call-back to compute reciprocity over transmitted hemisphere */
157			int
158			transHemi(const SDData bsdf, const FVECT v1, const RREAL v0, void *p)
159			{
160			return sampBSDFhemi(bsdf, 1 - 2*(v1[2]>0), v1, &diffRecip, p);
161			}
162
163	greg	2.1	/* Report reciprocity errors for the given directions */
164			void
165			checkReciprocity(const char *nm, const int side1, const int side2,
166			const SDData *bsdf, const int fl)
167			{
168	greg	2.2	SimpleStats myStats = SSinit;
169	greg	2.1	const SDSpectralDF *df = bsdf->tf;
170
171			if (side1 == side2) {
172			df = (side1 > 0) ? bsdf->rf : bsdf->rb;
173			if (!df) goto nothing2do;
174			} else if (!bsdf->tf \| !bsdf->tb)
175			goto nothing2do;
176
177			if (fl & F_MATRIX) { /* special case for matrix BSDF */
178			const SDMat m = (const SDMat )df->comp[0].dist;
179			int i = m->ninc;
180	greg	2.3	double diffuseY;
181	greg	2.1	FVECT vin, vout;
182	greg	2.2	double fwdY;
183			SDValue rev;
184	greg	2.3	if (side1 == side2)
185			diffuseY = (side1 > 0) ? bsdf->rLambFront.cieY : bsdf->rLambBack.cieY;
186			else
187			diffuseY = (side1 > 0) ? bsdf->tLambFront.cieY : bsdf->tLambBack.cieY;
188			diffuseY /= M_PI;
189	greg	2.1	while (i--) {
190			int o = m->nout;
191			if (!mBSDF_incvec(vin, m, i+.5))
192			continue;
193			while (o--) {
194			if (!mBSDF_outvec(vout, m, o+.5))
195			continue;
196	greg	2.3	fwdY = mBSDF_value(m, o, i) + diffuseY;
197	greg	2.2	if (fwdY <= 1e-4)
198	greg	2.1	continue;
199	greg	2.2	if (SDreportError( SDevalBSDF(&rev, vout, vin, bsdf), stderr))
200	greg	2.1	return;
201	greg	2.2	if (rev.cieY > 1e-4)
202			addStat(&myStats, rdiff(fwdY, rev.cieY));
203	greg	2.1	}
204			}
205	greg	2.3	} else if (fl & F_ISOTROPIC) { /* isotropic case */
206	greg	2.2	const double stepSize = sqrt(df->minProjSA/M_PI);
207			FVECT vin;
208			vin[1] = 0;
209			for (vin[0] = 0.5*stepSize; vin[0] < 1; vin[0] += stepSize) {
210			vin[2] = side1sqrt(1. - vin[0]vin[0]);
211			if (!sampBSDFhemi(bsdf, side2, vin, &diffRecip, &myStats))
212			return;
213			}
214			} else if (!sampBSDFhemi(bsdf, side1, NULL,
215			(side1==side2) ? &reflHemi : &transHemi, &myStats))
216			return;
217			if (myStats.nvals) {
218			printf("%s\t%5.1f\t%5.1f\t%5.1f\n", nm,
219			100.*myStats.vmin,
220			100.*myStats.vsum/(double)myStats.nvals,
221			100.*myStats.vmax);
222	greg	2.1	return;
223			}
224			nothing2do:
225			printf("%s\t0\t0\t0\n", nm);
226			}
227
228			/* Report on the given BSDF XML file */
229			int
230			checkXML(char *fname)
231			{
232			int flags;
233			SDData myBSDF;
234			char *pth;
235
236	greg	2.2	puts("=====================================================");
237	greg	2.1	printf("File: '%s'\n", fname);
238			SDclearBSDF(&myBSDF, fname);
239	greg	2.5	pth = getpath(fname, getrlibpath(), 0);
240	greg	2.1	if (!pth) {
241			fprintf(stderr, "Cannot find file '%s'\n", fname);
242			return 0;
243			}
244	greg	2.2	if (SDreportError( SDloadFile(&myBSDF, pth), stderr))
245			return 0;
246	greg	2.1	printf("Manufacturer: '%s'\n", myBSDF.makr);
247			printf("BSDF Name: '%s'\n", myBSDF.matn);
248			printf("Dimensions (W x H x Thickness): %g x %g x %g cm\n", 100.*myBSDF.dim[0],
249			100.myBSDF.dim[1], 100.myBSDF.dim[2]);
250			printf("Type: %s\n", getBSDFtype(&myBSDF, &flags));
251			printf("Color: %d\n", (flags & F_IN_COLOR) != 0);
252			printf("Has Geometry: %d\n", (myBSDF.mgf != NULL));
253			puts("Component\tLambertian XYZ %\tMax. Dir\tMin. Angle");
254	greg	2.4	detailComponent("Interior Refl", &myBSDF.rLambFront, myBSDF.rf);
255			detailComponent("Exterior Refl", &myBSDF.rLambBack, myBSDF.rb);
256	greg	2.1	detailComponent("Int->Ext Trans", &myBSDF.tLambFront, myBSDF.tf);
257			detailComponent("Ext->Int Trans", &myBSDF.tLambBack, myBSDF.tb);
258	greg	2.2	puts("Component\tReciprocity Error (min avg max %)");
259	greg	2.4	checkReciprocity("Interior Refl", 1, 1, &myBSDF, flags);
260			checkReciprocity("Exterior Refl", -1, -1, &myBSDF, flags);
261	greg	2.1	checkReciprocity("Transmission", -1, 1, &myBSDF, flags);
262			SDfreeBSDF(&myBSDF);
263			return 1;
264			}
265
266			int
267			main(int argc, char *argv[])
268			{
269			int i;
270
271			if (argc < 2) {
272			fprintf(stderr, "Usage: %s bsdf.xml ..\n", argv[0]);
273			return 1;
274			}
275	greg	2.2	for (i = 1; i < argc; i++)
276	greg	2.1	if (!checkXML(argv[i]))
277			return 1;
278			return 0;
279			}