src/cv/checkBSDF.c

#ifndef lint
static const char RCSid[] = "$Id: checkBSDF.c,v 2.8 2022/02/01 18:31:28 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)
{
        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.*(Lamb+df->maxHemi),
                                sqrt(df->minProjSA/M_PI)*(360./M_PI));
        else
                printf("%5.1f%%\t\t180 deg", 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;

        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\t  0\t  0\t  0\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. 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("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.9
Committed:	Fri Aug 12 23:55:00 2022 UTC (2 years, 8 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.8:	+9 -6 lines
Log Message:	fix(iso2klems): used solid angle, when projected SA needed, and integral factor
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.9	static const char RCSid[] = "$Id: checkBSDF.c,v 2.8 2022/02/01 18:31:28 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.9	double Lamb = 0;
84
85	greg	2.7	fputs(nm, stdout);
86	greg	2.9	if (lamb->spec.flags) {
87	greg	2.7	printf("\t%4.1f %4.1f %4.1f\t\t",
88	greg	2.2	100.lamb->cieYlamb->spec.cx/lamb->spec.cy,
89	greg	2.1	100.*lamb->cieY,
90			100.lamb->cieY(1.f - lamb->spec.cx - lamb->spec.cy)/lamb->spec.cy);
91	greg	2.9	Lamb = lamb->cieY;
92			} else
93	greg	2.7	fputs("\t 0 0 0\t\t", stdout);
94	greg	2.1	if (df)
95	greg	2.9	printf("%5.1f%%\t\t%.2f deg\n", 100.*(Lamb+df->maxHemi),
96	greg	2.1	sqrt(df->minProjSA/M_PI)*(360./M_PI));
97			else
98	greg	2.9	printf("%5.1f%%\t\t180 deg", Lamb);
99	greg	2.1	}
100
101	greg	2.2	/* Add a value to stats */
102			void
103			addStat(SimpleStats *ssp, double v)
104			{
105			if (v < ssp->vmin) ssp->vmin = v;
106			if (v > ssp->vmax) ssp->vmax = v;
107			ssp->vsum += v;
108			ssp->nvals++;
109			}
110
111			/* Sample a BSDF hemisphere with callback (quadtree recursion) */
112			int
113			qtSampBSDF(double xleft, double ytop, double siz,
114			const SDData bsdf, const int side, const RREAL v0,
115			int (cf)(const SDData b, const FVECT v1, const RREAL v0, void p),
116			void *cdata)
117			{
118			if (siz < 0.124) { /* make sure we subdivide, first */
119			FVECT vsmp;
120			double sa;
121			square2disk(vsmp, xleft + frandom()siz, ytop + frandom()siz);
122			vsmp[2] = 1. - vsmp[0]vsmp[0] - vsmp[1]vsmp[1];
123			if (vsmp[2] <= 0) return 0;
124			vsmp[2] = side * sqrt(vsmp[2]);
125			if (SDreportError( SDsizeBSDF(&sa, vsmp, v0, SDqueryMin, bsdf), stderr))
126			return 0;
127			if (sa >= M_PIsizsiz - FTINY) /* no further division needed */
128			return (*cf)(bsdf, vsmp, v0, cdata);
129			}
130			siz = .5; / 4-branch recursion */
131			return( qtSampBSDF(xleft, ytop, siz, bsdf, side, v0, cf, cdata) &&
132			qtSampBSDF(xleft+siz, ytop, siz, bsdf, side, v0, cf, cdata) &&
133			qtSampBSDF(xleft, ytop+siz, siz, bsdf, side, v0, cf, cdata) &&
134			qtSampBSDF(xleft+siz, ytop+siz, siz, bsdf, side, v0, cf, cdata) );
135			}
136
137			#define sampBSDFhemi(b,s,v0,cf,cd) qtSampBSDF(0,0,1,b,s,v0,cf,cd)
138
139			/* Call-back to compute reciprocity difference */
140			int
141			diffRecip(const SDData bsdf, const FVECT v1, const RREAL v0, void *p)
142			{
143			SDValue sdv;
144			double otherY;
145
146			if (SDreportError( SDevalBSDF(&sdv, v0, v1, bsdf), stderr))
147			return 0;
148			otherY = sdv.cieY;
149			if (SDreportError( SDevalBSDF(&sdv, v1, v0, bsdf), stderr))
150			return 0;
151
152			addStat((SimpleStats *)p, rdiff(sdv.cieY, otherY));
153			return 1;
154			}
155
156			/* Call-back to compute reciprocity over reflected hemisphere */
157			int
158			reflHemi(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			/* Call-back to compute reciprocity over transmitted hemisphere */
164			int
165			transHemi(const SDData bsdf, const FVECT v1, const RREAL v0, void *p)
166			{
167			return sampBSDFhemi(bsdf, 1 - 2*(v1[2]>0), v1, &diffRecip, p);
168			}
169
170	greg	2.1	/* Report reciprocity errors for the given directions */
171			void
172			checkReciprocity(const char *nm, const int side1, const int side2,
173			const SDData *bsdf, const int fl)
174			{
175	greg	2.2	SimpleStats myStats = SSinit;
176	greg	2.1	const SDSpectralDF *df = bsdf->tf;
177
178			if (side1 == side2) {
179			df = (side1 > 0) ? bsdf->rf : bsdf->rb;
180			if (!df) goto nothing2do;
181			} else if (!bsdf->tf \| !bsdf->tb)
182			goto nothing2do;
183
184			if (fl & F_MATRIX) { /* special case for matrix BSDF */
185			const SDMat m = (const SDMat )df->comp[0].dist;
186			int i = m->ninc;
187	greg	2.3	double diffuseY;
188	greg	2.1	FVECT vin, vout;
189	greg	2.2	double fwdY;
190			SDValue rev;
191	greg	2.3	if (side1 == side2)
192			diffuseY = (side1 > 0) ? bsdf->rLambFront.cieY : bsdf->rLambBack.cieY;
193			else
194			diffuseY = (side1 > 0) ? bsdf->tLambFront.cieY : bsdf->tLambBack.cieY;
195			diffuseY /= M_PI;
196	greg	2.1	while (i--) {
197			int o = m->nout;
198			if (!mBSDF_incvec(vin, m, i+.5))
199			continue;
200			while (o--) {
201			if (!mBSDF_outvec(vout, m, o+.5))
202			continue;
203	greg	2.3	fwdY = mBSDF_value(m, o, i) + diffuseY;
204	greg	2.2	if (fwdY <= 1e-4)
205	greg	2.1	continue;
206	greg	2.2	if (SDreportError( SDevalBSDF(&rev, vout, vin, bsdf), stderr))
207	greg	2.1	return;
208	greg	2.2	if (rev.cieY > 1e-4)
209			addStat(&myStats, rdiff(fwdY, rev.cieY));
210	greg	2.1	}
211			}
212	greg	2.3	} else if (fl & F_ISOTROPIC) { /* isotropic case */
213	greg	2.2	const double stepSize = sqrt(df->minProjSA/M_PI);
214			FVECT vin;
215			vin[1] = 0;
216			for (vin[0] = 0.5*stepSize; vin[0] < 1; vin[0] += stepSize) {
217			vin[2] = side1sqrt(1. - vin[0]vin[0]);
218			if (!sampBSDFhemi(bsdf, side2, vin, &diffRecip, &myStats))
219			return;
220			}
221			} else if (!sampBSDFhemi(bsdf, side1, NULL,
222			(side1==side2) ? &reflHemi : &transHemi, &myStats))
223			return;
224			if (myStats.nvals) {
225			printf("%s\t%5.1f\t%5.1f\t%5.1f\n", nm,
226			100.*myStats.vmin,
227			100.*myStats.vsum/(double)myStats.nvals,
228			100.*myStats.vmax);
229	greg	2.1	return;
230			}
231			nothing2do:
232	greg	2.7	printf("%s\t 0\t 0\t 0\n", nm);
233	greg	2.1	}
234
235			/* Report on the given BSDF XML file */
236			int
237			checkXML(char *fname)
238			{
239			int flags;
240			SDData myBSDF;
241			char *pth;
242
243	greg	2.2	puts("=====================================================");
244	greg	2.1	printf("File: '%s'\n", fname);
245			SDclearBSDF(&myBSDF, fname);
246	greg	2.5	pth = getpath(fname, getrlibpath(), 0);
247	greg	2.1	if (!pth) {
248			fprintf(stderr, "Cannot find file '%s'\n", fname);
249			return 0;
250			}
251	greg	2.2	if (SDreportError( SDloadFile(&myBSDF, pth), stderr))
252			return 0;
253	greg	2.1	printf("Manufacturer: '%s'\n", myBSDF.makr);
254			printf("BSDF Name: '%s'\n", myBSDF.matn);
255			printf("Dimensions (W x H x Thickness): %g x %g x %g cm\n", 100.*myBSDF.dim[0],
256			100.myBSDF.dim[1], 100.myBSDF.dim[2]);
257			printf("Type: %s\n", getBSDFtype(&myBSDF, &flags));
258			printf("Color: %d\n", (flags & F_IN_COLOR) != 0);
259			printf("Has Geometry: %d\n", (myBSDF.mgf != NULL));
260	greg	2.9	puts("Component\tLambertian XYZ (%)\tMax. Tot. Dir\tMin. Angle");
261	greg	2.4	detailComponent("Interior Refl", &myBSDF.rLambFront, myBSDF.rf);
262			detailComponent("Exterior Refl", &myBSDF.rLambBack, myBSDF.rb);
263	greg	2.1	detailComponent("Int->Ext Trans", &myBSDF.tLambFront, myBSDF.tf);
264			detailComponent("Ext->Int Trans", &myBSDF.tLambBack, myBSDF.tb);
265	greg	2.2	puts("Component\tReciprocity Error (min avg max %)");
266	greg	2.4	checkReciprocity("Interior Refl", 1, 1, &myBSDF, flags);
267			checkReciprocity("Exterior Refl", -1, -1, &myBSDF, flags);
268	greg	2.1	checkReciprocity("Transmission", -1, 1, &myBSDF, flags);
269			SDfreeBSDF(&myBSDF);
270			return 1;
271			}
272
273			int
274			main(int argc, char *argv[])
275			{
276			int i;
277
278			if (argc < 2) {
279			fprintf(stderr, "Usage: %s bsdf.xml ..\n", argv[0]);
280			return 1;
281			}
282	greg	2.2	for (i = 1; i < argc; i++)
283	greg	2.1	if (!checkXML(argv[i]))
284			return 1;
285			return 0;
286			}