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.2 by greg, Wed Dec 15 02:13:27 2021 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 98 | Line 178 | checkReciprocity(const char *nm, const int side1, cons
178                  const SDMat     *m = (const SDMat *)df->comp[0].dist;
179                  int             i = m->ninc;
180                  FVECT           vin, vout;
181 <                double          rerr;
182 <                SDValue         vrev;
181 >                double          fwdY;
182 >                SDValue         rev;
183                  while (i--) {
184                      int o = m->nout;
185                      if (!mBSDF_incvec(vin, m, i+.5))
# Line 107 | Line 187 | checkReciprocity(const char *nm, const int side1, cons
187                      while (o--) {
188                          if (!mBSDF_outvec(vout, m, o+.5))
189                                  continue;
190 <                        rerr = mBSDF_value(m, o, i);
191 <                        if (rerr <= FTINY)
190 >                        fwdY = mBSDF_value(m, o, i);
191 >                        if (fwdY <= 1e-4)
192                                  continue;
193 <                        if (SDreportError( SDevalBSDF(&vrev, vout, vin, bsdf), stderr))
193 >                        if (SDreportError( SDevalBSDF(&rev, vout, vin, bsdf), stderr))
194                                  return;
195 <                        rerr = 100.*fabs(rerr - vrev.cieY)/rerr;
196 <                        if (rerr < emin) emin = rerr;
117 <                        if (rerr > emax) emax = rerr;
118 <                        esum += rerr;
119 <                        ++ntested;
195 >                        if (rev.cieY > 1e-4)
196 >                                addStat(&myStats, rdiff(fwdY, rev.cieY));
197                      }
198                  }
199 <        } else {
200 <        }
201 <        if (ntested) {
202 <                printf("%s\t%.1f\t%.1f\t%.1f\n", nm, emin, esum/(double)ntested, emax);
199 >        } if (fl & F_ISOTROPIC) {               /* isotropic case */
200 >                const double    stepSize = sqrt(df->minProjSA/M_PI);
201 >                FVECT           vin;
202 >                vin[1] = 0;
203 >                for (vin[0] = 0.5*stepSize; vin[0] < 1; vin[0] += stepSize) {
204 >                        vin[2] = side1*sqrt(1. - vin[0]*vin[0]);
205 >                        if (!sampBSDFhemi(bsdf, side2, vin, &diffRecip, &myStats))
206 >                                return;
207 >                }
208 >        } else if (!sampBSDFhemi(bsdf, side1, NULL,
209 >                                (side1==side2) ? &reflHemi : &transHemi, &myStats))
210                  return;
211 +        if (myStats.nvals) {
212 +                printf("%s\t%5.1f\t%5.1f\t%5.1f\n", nm,
213 +                                100.*myStats.vmin,
214 +                                100.*myStats.vsum/(double)myStats.nvals,
215 +                                100.*myStats.vmax);
216 +                return;
217          }
218   nothing2do:
219          printf("%s\t0\t0\t0\n", nm);
# Line 134 | Line 224 | int
224   checkXML(char *fname)
225   {
226          int             flags;
137 –        SDError         ec;
227          SDData          myBSDF;
228          char            *pth;
229  
230 +        puts("=====================================================");
231          printf("File: '%s'\n", fname);
232          SDclearBSDF(&myBSDF, fname);
233          pth = getpath(fname, getrlibpath(), R_OK);
# Line 145 | Line 235 | checkXML(char *fname)
235                  fprintf(stderr, "Cannot find file '%s'\n", fname);
236                  return 0;
237          }
238 <        ec = SDloadFile(&myBSDF, pth);
239 <        if (ec) goto err;
238 >        if (SDreportError( SDloadFile(&myBSDF, pth), stderr))
239 >                return 0;
240          printf("Manufacturer: '%s'\n", myBSDF.makr);
241          printf("BSDF Name: '%s'\n", myBSDF.matn);
242          printf("Dimensions (W x H x Thickness): %g x %g x %g cm\n", 100.*myBSDF.dim[0],
# Line 159 | Line 249 | checkXML(char *fname)
249          detailComponent("External Refl", &myBSDF.rLambBack, myBSDF.rb);
250          detailComponent("Int->Ext Trans", &myBSDF.tLambFront, myBSDF.tf);
251          detailComponent("Ext->Int Trans", &myBSDF.tLambBack, myBSDF.tb);
252 <        puts("Component\tReciprocity Error (min/avg/max %)");
252 >        puts("Component\tReciprocity Error (min avg max %)");
253          checkReciprocity("Front Refl", 1, 1, &myBSDF, flags);
254          checkReciprocity("Back Refl", -1, -1, &myBSDF, flags);
255          checkReciprocity("Transmission", -1, 1, &myBSDF, flags);
256          SDfreeBSDF(&myBSDF);
257          return 1;
168 – err:
169 –        SDreportError(ec, stderr);
170 –        SDfreeBSDF(&myBSDF);
171 –        return 0;
258   }
259  
260   int
# Line 180 | Line 266 | main(int argc, char *argv[])
266                  fprintf(stderr, "Usage: %s bsdf.xml ..\n", argv[0]);
267                  return 1;
268          }
269 <        for (i = 1; i < argc; i++) {
184 <                puts("=====================================================");
269 >        for (i = 1; i < argc; i++)
270                  if (!checkXML(argv[i]))
271                          return 1;
187 –        }
272          return 0;
273   }

Diff Legend

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