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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines