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.5 by greg, Wed Jan 26 16:19:30 2022 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) : ((b)-(a))/(b))
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 >        printf("%s\t%4.1f %4.1f %4.1f\t\t", nm,
84 >                        100.*lamb->cieY*lamb->spec.cx/lamb->spec.cy,
85                          100.*lamb->cieY,
86                          100.*lamb->cieY*(1.f - lamb->spec.cx - lamb->spec.cy)/lamb->spec.cy);
87          if (df)
# Line 78 | Line 91 | detailComponent(const char *nm, const SDValue *lamb, c
91                  puts("0%\t\t180");
92   }
93  
94 + /* 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_PI*siz*siz - 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   /* 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 +        SimpleStats             myStats = SSinit;
169          const SDSpectralDF      *df = bsdf->tf;
87 –        double                  emin=FHUGE, emax=0, esum=0;
88 –        int                     ntested=0;
89 –        int                     ec;
170  
171          if (side1 == side2) {
172                  df = (side1 > 0) ? bsdf->rf : bsdf->rb;
# Line 97 | Line 177 | checkReciprocity(const char *nm, const int side1, cons
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 +                double          diffuseY;
181                  FVECT           vin, vout;
182 <                double          rerr;
183 <                SDValue         vrev;
182 >                double          fwdY;
183 >                SDValue         rev;
184 >                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                  while (i--) {
190                      int o = m->nout;
191                      if (!mBSDF_incvec(vin, m, i+.5))
# Line 107 | Line 193 | checkReciprocity(const char *nm, const int side1, cons
193                      while (o--) {
194                          if (!mBSDF_outvec(vout, m, o+.5))
195                                  continue;
196 <                        rerr = mBSDF_value(m, o, i);
197 <                        if (rerr <= FTINY)
196 >                        fwdY = mBSDF_value(m, o, i) + diffuseY;
197 >                        if (fwdY <= 1e-4)
198                                  continue;
199 <                        if (SDreportError( SDevalBSDF(&vrev, vout, vin, bsdf), stderr))
199 >                        if (SDreportError( SDevalBSDF(&rev, vout, vin, bsdf), stderr))
200                                  return;
201 <                        rerr = 100.*fabs(rerr - vrev.cieY)/rerr;
202 <                        if (rerr < emin) emin = rerr;
117 <                        if (rerr > emax) emax = rerr;
118 <                        esum += rerr;
119 <                        ++ntested;
201 >                        if (rev.cieY > 1e-4)
202 >                                addStat(&myStats, rdiff(fwdY, rev.cieY));
203                      }
204                  }
205 <        } else {
206 <        }
207 <        if (ntested) {
208 <                printf("%s\t%.1f\t%.1f\t%.1f\n", nm, emin, esum/(double)ntested, emax);
205 >        } else if (fl & F_ISOTROPIC) {          /* isotropic case */
206 >                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] = side1*sqrt(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 +                return;
223          }
224   nothing2do:
225          printf("%s\t0\t0\t0\n", nm);
# Line 134 | Line 230 | int
230   checkXML(char *fname)
231   {
232          int             flags;
137 –        SDError         ec;
233          SDData          myBSDF;
234          char            *pth;
235  
236 +        puts("=====================================================");
237          printf("File: '%s'\n", fname);
238          SDclearBSDF(&myBSDF, fname);
239 <        pth = getpath(fname, getrlibpath(), R_OK);
239 >        pth = getpath(fname, getrlibpath(), 0);
240          if (!pth) {
241                  fprintf(stderr, "Cannot find file '%s'\n", fname);
242                  return 0;
243          }
244 <        ec = SDloadFile(&myBSDF, pth);
245 <        if (ec) goto err;
244 >        if (SDreportError( SDloadFile(&myBSDF, pth), stderr))
245 >                return 0;
246          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],
# Line 155 | Line 251 | checkXML(char *fname)
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 <        detailComponent("Internal Refl", &myBSDF.rLambFront, myBSDF.rf);
255 <        detailComponent("External Refl", &myBSDF.rLambBack, myBSDF.rb);
254 >        detailComponent("Interior Refl", &myBSDF.rLambFront, myBSDF.rf);
255 >        detailComponent("Exterior Refl", &myBSDF.rLambBack, myBSDF.rb);
256          detailComponent("Int->Ext Trans", &myBSDF.tLambFront, myBSDF.tf);
257          detailComponent("Ext->Int Trans", &myBSDF.tLambBack, myBSDF.tb);
258 <        puts("Component\tReciprocity Error (min/avg/max %)");
259 <        checkReciprocity("Front Refl", 1, 1, &myBSDF, flags);
260 <        checkReciprocity("Back Refl", -1, -1, &myBSDF, flags);
258 >        puts("Component\tReciprocity Error (min avg max %)");
259 >        checkReciprocity("Interior Refl", 1, 1, &myBSDF, flags);
260 >        checkReciprocity("Exterior Refl", -1, -1, &myBSDF, flags);
261          checkReciprocity("Transmission", -1, 1, &myBSDF, flags);
262          SDfreeBSDF(&myBSDF);
263          return 1;
168 – err:
169 –        SDreportError(ec, stderr);
170 –        SDfreeBSDF(&myBSDF);
171 –        return 0;
264   }
265  
266   int
# Line 180 | Line 272 | main(int argc, char *argv[])
272                  fprintf(stderr, "Usage: %s bsdf.xml ..\n", argv[0]);
273                  return 1;
274          }
275 <        for (i = 1; i < argc; i++) {
184 <                puts("=====================================================");
275 >        for (i = 1; i < argc; i++)
276                  if (!checkXML(argv[i]))
277                          return 1;
187 –        }
278          return 0;
279   }

Diff Legend

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