ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/checkBSDF.c
(Generate patch)

Comparing ray/src/cv/checkBSDF.c (file contents):
Revision 2.1 by greg, Tue Dec 14 02:33:18 2021 UTC vs.
Revision 2.10 by greg, Mon Sep 11 18:43:36 2023 UTC

# Line 10 | Line 10 | static const char RCSid[] = "$Id$";
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"
# Line 19 | Line 20 | static const char RCSid[] = "$Id$";
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)
# Line 68 | Line 80 | getBSDFtype(const SDData *bsdf, int *flags)
80   void
81   detailComponent(const char *nm, const SDValue *lamb, const SDSpectralDF *df)
82   {
83 <        printf("%s\t%4.1f %4.1f %4.1f\t\t", nm, 100.*lamb->cieY*lamb->spec.cx/lamb->spec.cy,
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.*df->maxHemi,
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 <                puts("0%\t\t180");
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;
87 –        double                  emin=FHUGE, emax=0, esum=0;
88 –        int                     ntested=0;
89 –        int                     ec;
177  
178          if (side1 == side2) {
179                  df = (side1 > 0) ? bsdf->rf : bsdf->rb;
# Line 97 | Line 184 | checkReciprocity(const char *nm, const int side1, cons
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          rerr;
190 <                SDValue         vrev;
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))
# Line 107 | Line 200 | checkReciprocity(const char *nm, const int side1, cons
200                      while (o--) {
201                          if (!mBSDF_outvec(vout, m, o+.5))
202                                  continue;
203 <                        rerr = mBSDF_value(m, o, i);
204 <                        if (rerr <= FTINY)
203 >                        fwdY = mBSDF_value(m, o, i) + diffuseY;
204 >                        if (fwdY <= 1e-4)
205                                  continue;
206 <                        if (SDreportError( SDevalBSDF(&vrev, vout, vin, bsdf), stderr))
206 >                        if (SDreportError( SDevalBSDF(&rev, vout, vin, bsdf), stderr))
207                                  return;
208 <                        rerr = 100.*fabs(rerr - vrev.cieY)/rerr;
209 <                        if (rerr < emin) emin = rerr;
117 <                        if (rerr > emax) emax = rerr;
118 <                        esum += rerr;
119 <                        ++ntested;
208 >                        if (rev.cieY > 1e-4)
209 >                                addStat(&myStats, rdiff(fwdY, rev.cieY));
210                      }
211                  }
212 <        } else {
213 <        }
214 <        if (ntested) {
215 <                printf("%s\t%.1f\t%.1f\t%.1f\n", nm, emin, esum/(double)ntested, emax);
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\t0\t0\t0\n", nm);
232 >        printf("%s\t  0\t  0\t  0\n", nm);
233   }
234  
235   /* Report on the given BSDF XML file */
# Line 134 | Line 237 | int
237   checkXML(char *fname)
238   {
239          int             flags;
137 –        SDError         ec;
240          SDData          myBSDF;
241          char            *pth;
242  
243 +        puts("=====================================================");
244          printf("File: '%s'\n", fname);
245          SDclearBSDF(&myBSDF, fname);
246 <        pth = getpath(fname, getrlibpath(), R_OK);
246 >        pth = getpath(fname, getrlibpath(), 0);
247          if (!pth) {
248                  fprintf(stderr, "Cannot find file '%s'\n", fname);
249                  return 0;
250          }
251 <        ec = SDloadFile(&myBSDF, pth);
252 <        if (ec) goto err;
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],
# Line 154 | Line 257 | checkXML(char *fname)
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. Dir\tMin. Angle");
261 <        detailComponent("Internal Refl", &myBSDF.rLambFront, myBSDF.rf);
262 <        detailComponent("External Refl", &myBSDF.rLambBack, myBSDF.rb);
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("Front Refl", 1, 1, &myBSDF, flags);
267 <        checkReciprocity("Back Refl", -1, -1, &myBSDF, flags);
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;
168 – err:
169 –        SDreportError(ec, stderr);
170 –        SDfreeBSDF(&myBSDF);
171 –        return 0;
271   }
272  
273   int
# Line 180 | Line 279 | main(int argc, char *argv[])
279                  fprintf(stderr, "Usage: %s bsdf.xml ..\n", argv[0]);
280                  return 1;
281          }
282 <        for (i = 1; i < argc; i++) {
184 <                puts("=====================================================");
282 >        for (i = 1; i < argc; i++)
283                  if (!checkXML(argv[i]))
284                          return 1;
187 –        }
285          return 0;
286   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines