ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/checkBSDF.c
Revision: 2.10
Committed: Mon Sep 11 18:43:36 2023 UTC (9 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R4, HEAD
Changes since 2.9: +2 -2 lines
Log Message:
fix(checkBSDF): added missing newline in output

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: checkBSDF.c,v 2.9 2022/08/12 23:55:00 greg Exp $";
3 #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 #include "random.h"
14 #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 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 #define rdiff(a,b) ((a)>(b) ? ((a)-(b))/((a)+FTINY) : ((b)-(a))/((b)+FTINY))
33
34 /* 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 double Lamb = 0;
84
85 fputs(nm, stdout);
86 if (lamb->spec.flags) {
87 printf("\t%4.1f %4.1f %4.1f\t\t",
88 100.*lamb->cieY*lamb->spec.cx/lamb->spec.cy,
89 100.*lamb->cieY,
90 100.*lamb->cieY*(1.f - lamb->spec.cx - lamb->spec.cy)/lamb->spec.cy);
91 Lamb = lamb->cieY;
92 } else
93 fputs("\t 0 0 0\t\t", stdout);
94 if (df)
95 printf("%5.1f%%\t\t%.2f deg\n", 100.*(Lamb+df->maxHemi),
96 sqrt(df->minProjSA/M_PI)*(360./M_PI));
97 else
98 printf("%5.1f%%\t\t180 deg\n", Lamb);
99 }
100
101 /* 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_PI*siz*siz - 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 /* 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 SimpleStats myStats = SSinit;
176 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 double diffuseY;
188 FVECT vin, vout;
189 double fwdY;
190 SDValue rev;
191 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 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 fwdY = mBSDF_value(m, o, i) + diffuseY;
204 if (fwdY <= 1e-4)
205 continue;
206 if (SDreportError( SDevalBSDF(&rev, vout, vin, bsdf), stderr))
207 return;
208 if (rev.cieY > 1e-4)
209 addStat(&myStats, rdiff(fwdY, rev.cieY));
210 }
211 }
212 } else if (fl & F_ISOTROPIC) { /* isotropic case */
213 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] = side1*sqrt(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 return;
230 }
231 nothing2do:
232 printf("%s\t 0\t 0\t 0\n", nm);
233 }
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 puts("=====================================================");
244 printf("File: '%s'\n", fname);
245 SDclearBSDF(&myBSDF, fname);
246 pth = getpath(fname, getrlibpath(), 0);
247 if (!pth) {
248 fprintf(stderr, "Cannot find file '%s'\n", fname);
249 return 0;
250 }
251 if (SDreportError( SDloadFile(&myBSDF, pth), stderr))
252 return 0;
253 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 puts("Component\tLambertian XYZ (%)\tMax. Tot. Dir\tMin. Angle");
261 detailComponent("Interior Refl", &myBSDF.rLambFront, myBSDF.rf);
262 detailComponent("Exterior Refl", &myBSDF.rLambBack, myBSDF.rb);
263 detailComponent("Int->Ext Trans", &myBSDF.tLambFront, myBSDF.tf);
264 detailComponent("Ext->Int Trans", &myBSDF.tLambBack, myBSDF.tb);
265 puts("Component\tReciprocity Error (min avg max %)");
266 checkReciprocity("Interior Refl", 1, 1, &myBSDF, flags);
267 checkReciprocity("Exterior Refl", -1, -1, &myBSDF, flags);
268 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 for (i = 1; i < argc; i++)
283 if (!checkXML(argv[i]))
284 return 1;
285 return 0;
286 }