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

Comparing src/cv/bsdfrbf.c (file contents):
Revision 2.14 by greg, Mon Oct 21 22:03:03 2013 UTC vs.
Revision 2.32 by greg, Tue Apr 9 22:39:33 2019 UTC

# Line 29 | Line 29 | static const char RCSid[] = "$Id$";
29   #include "bsdfrep.h"
30  
31   #ifndef RSCA
32 < #define RSCA            2.2             /* radius scaling factor (empirical) */
32 > #define RSCA            2.0             /* radius scaling factor (empirical) */
33   #endif
34 + #ifndef MAXSLOPE
35 + #define MAXSLOPE        200.0           /* maximum slope for smooth region */
36 + #endif
37   #ifndef SMOOTH_MSE
38   #define SMOOTH_MSE      5e-5            /* acceptable mean squared error */
39   #endif
40   #ifndef SMOOTH_MSER
41 < #define SMOOTH_MSER     0.07            /* acceptable relative MSE */
41 > #define SMOOTH_MSER     0.0025          /* acceptable relative MSE */
42   #endif
43   #define MAX_RAD         (GRIDRES/8)     /* maximum lobe radius */
44  
45   #define RBFALLOCB       10              /* RBF allocation block size */
46  
47 <                                /* our loaded grid for this incident angle */
47 >                                        /* loaded grid or comparison DSFs */
48   GRIDVAL                 dsf_grid[GRIDRES][GRIDRES];
49 +                                        /* allocated chrominance sums if any */
50 + float                   (*spec_grid)[GRIDRES][GRIDRES];
51 + int                     nspec_grid = 0;
52  
53 + /* Set up visible spectrum sampling */
54 + void
55 + set_spectral_samples(int nspec)
56 + {
57 +        if (rbf_colorimetry == RBCunknown) {
58 +                if (nspec_grid > 0) {
59 +                        free(spec_grid);
60 +                        spec_grid = NULL;
61 +                        nspec_grid = 0;
62 +                }
63 +                if (nspec == 1) {
64 +                        rbf_colorimetry = RBCphotopic;
65 +                        return;
66 +                }
67 +                if (nspec == 3) {
68 +                        rbf_colorimetry = RBCtristimulus;
69 +                        spec_grid = (float (*)[GRIDRES][GRIDRES])calloc(
70 +                                        2*GRIDRES*GRIDRES, sizeof(float) );
71 +                        if (spec_grid == NULL)
72 +                                goto mem_error;
73 +                        nspec_grid = 2;
74 +                        return;
75 +                }
76 +                fprintf(stderr,
77 +                        "%s: only 1 or 3 spectral samples currently supported\n",
78 +                                progname);
79 +                exit(1);
80 +        }
81 +        if (nspec != nspec_grid+1) {
82 +                fprintf(stderr,
83 +                        "%s: number of spectral samples cannot be changed\n",
84 +                                progname);
85 +                exit(1);
86 +        }
87 +        return;
88 + mem_error:
89 +        fprintf(stderr, "%s: out of memory in set_spectral_samples()\n",
90 +                        progname);
91 +        exit(1);
92 + }
93 +
94   /* Start new DSF input grid */
95   void
96   new_bsdf_data(double new_theta, double new_phi)
# Line 51 | Line 98 | new_bsdf_data(double new_theta, double new_phi)
98          if (!new_input_direction(new_theta, new_phi))
99                  exit(1);
100          memset(dsf_grid, 0, sizeof(dsf_grid));
101 +        if (nspec_grid > 0)
102 +                memset(spec_grid, 0, sizeof(float)*GRIDRES*GRIDRES*nspec_grid);
103   }
104  
105   /* Add BSDF data point */
106   void
107 < add_bsdf_data(double theta_out, double phi_out, double val, int isDSF)
107 > add_bsdf_data(double theta_out, double phi_out, const double val[], int isDSF)
108   {
109          FVECT   ovec;
110 +        double  cfact, Yval;
111          int     pos[2];
112  
113 +        if (nspec_grid > 2) {
114 +                fprintf(stderr, "%s: unsupported color space\n", progname);
115 +                exit(1);
116 +        }
117          if (!output_orient)             /* check output orientation */
118                  output_orient = 1 - 2*(theta_out > 90.);
119          else if (output_orient > 0 ^ theta_out < 90.) {
120 <                fputs("Cannot handle output angles on both sides of surface\n",
121 <                                stderr);
120 >                fprintf(stderr,
121 >                "%s: cannot handle output angles on both sides of surface\n",
122 >                                progname);
123                  exit(1);
124          }
125          ovec[2] = sin((M_PI/180.)*theta_out);
126          ovec[0] = cos((M_PI/180.)*phi_out) * ovec[2];
127          ovec[1] = sin((M_PI/180.)*phi_out) * ovec[2];
128          ovec[2] = sqrt(1. - ovec[2]*ovec[2]);
129 +                                        /* BSDF to DSF correction */
130 +        cfact = isDSF ? 1. : ovec[2];
131  
132 <        if (val <= 0)                   /* truncate to zero */
76 <                val = 0;
77 <        else if (!isDSF)
78 <                val *= ovec[2];         /* convert from BSDF to DSF */
79 <
132 >        Yval = cfact * val[rbf_colorimetry==RBCtristimulus];
133                                          /* update BSDF histogram */
134 <        if (val < BSDF2BIG*ovec[2] && val > BSDF2SML*ovec[2])
135 <                ++bsdf_hist[histndx(val/ovec[2])];
134 >        if (BSDF2SML*ovec[2] < Yval && Yval < BSDF2BIG*ovec[2])
135 >                ++bsdf_hist[histndx(Yval/ovec[2])];
136  
137          pos_from_vec(pos, ovec);
138  
139 <        dsf_grid[pos[0]][pos[1]].vsum += val;
140 <        dsf_grid[pos[0]][pos[1]].nval++;
139 >        dsf_grid[pos[0]][pos[1]].sum.v += Yval;
140 >        dsf_grid[pos[0]][pos[1]].sum.n++;
141 >                                        /* add in X and Z values */
142 >        if (rbf_colorimetry == RBCtristimulus) {
143 >                spec_grid[0][pos[0]][pos[1]] += cfact * val[0];
144 >                spec_grid[1][pos[0]][pos[1]] += cfact * val[2];
145 >        }
146   }
147  
148   /* Compute minimum BSDF from histogram (does not clear) */
149   static void
150   comp_bsdf_min()
151   {
152 <        int     cnt;
153 <        int     i, target;
152 >        unsigned long   cnt, target;
153 >        int             i;
154  
155          cnt = 0;
156          for (i = HISTLEN; i--; )
# Line 116 | Line 174 | empty_region(int x0, int x1, int y0, int y1)
174  
175          for (x = x0; x < x1; x++)
176              for (y = y0; y < y1; y++)
177 <                if (dsf_grid[x][y].nval)
177 >                if (dsf_grid[x][y].sum.n)
178                          return(0);
179          return(1);
180   }
# Line 134 | Line 192 | smooth_region(int x0, int x1, int y0, int y1)
192          memset(xvec, 0, sizeof(xvec));
193          for (x = x0; x < x1; x++)
194              for (y = y0; y < y1; y++)
195 <                if ((n = dsf_grid[x][y].nval) > 0) {
196 <                        double  z = dsf_grid[x][y].vsum;
195 >                if ((n = dsf_grid[x][y].sum.n) > 0) {
196 >                        double  z = dsf_grid[x][y].sum.v;
197                          rMtx[0][0] += x*x*(double)n;
198                          rMtx[0][1] += x*y*(double)n;
199                          rMtx[0][2] += x*(double)n;
# Line 147 | Line 205 | smooth_region(int x0, int x1, int y0, int y1)
205                          xvec[2] += z;
206                  }
207          rMtx[1][0] = rMtx[0][1];
208 +        rMtx[2][0] = rMtx[0][2];
209          rMtx[2][1] = rMtx[1][2];
210          nvs = rMtx[2][2];
211          if (SDinvXform(rMtx, rMtx) != SDEnone)
212 <                return(0);
212 >                return(1);              /* colinear values */
213          A = DOT(rMtx[0], xvec);
214          B = DOT(rMtx[1], xvec);
215 +        if (A*A + B*B > MAXSLOPE*MAXSLOPE)      /* too steep? */
216 +                return(0);
217          C = DOT(rMtx[2], xvec);
218          sqerr = 0.0;                    /* compute mean squared error */
219          for (x = x0; x < x1; x++)
220              for (y = y0; y < y1; y++)
221 <                if ((n = dsf_grid[x][y].nval) > 0) {
222 <                        double  d = A*x + B*y + C - dsf_grid[x][y].vsum/n;
221 >                if ((n = dsf_grid[x][y].sum.n) > 0) {
222 >                        double  d = A*x + B*y + C - dsf_grid[x][y].sum.v/n;
223                          sqerr += n*d*d;
224                  }
225          if (sqerr <= nvs*SMOOTH_MSE)    /* below absolute MSE threshold? */
# Line 168 | Line 229 | smooth_region(int x0, int x1, int y0, int y1)
229   }
230  
231   /* Create new lobe based on integrated samples in region */
232 < static void
232 > static int
233   create_lobe(RBFVAL *rvp, int x0, int x1, int y0, int y1)
234   {
235          double  vtot = 0.0;
236 +        double  CIEXtot = 0.0, CIEZtot = 0.0;
237          int     nv = 0;
238 +        double  wxsum = 0.0, wysum = 0.0, wtsum = 0.0;
239          double  rad;
240          int     x, y;
241                                          /* compute average for region */
242          for (x = x0; x < x1; x++)
243 <            for (y = y0; y < y1; y++) {
244 <                vtot += dsf_grid[x][y].vsum;
245 <                nv += dsf_grid[x][y].nval;
246 <            }
243 >            for (y = y0; y < y1; y++)
244 >                if (dsf_grid[x][y].sum.n) {
245 >                    const double        v = dsf_grid[x][y].sum.v;
246 >                    const int           n = dsf_grid[x][y].sum.n;
247 >
248 >                    if (v > 0) {
249 >                        const double    wt = v / (double)n;
250 >                        wxsum += wt * x;
251 >                        wysum += wt * y;
252 >                        wtsum += wt;
253 >                    }
254 >                    vtot += v;
255 >                    nv += n;
256 >                    if (rbf_colorimetry == RBCtristimulus) {
257 >                        CIEXtot += spec_grid[0][x][y];
258 >                        CIEZtot += spec_grid[1][x][y];
259 >                    }
260 >                }
261          if (!nv) {
262                  fprintf(stderr, "%s: internal - missing samples in create_lobe\n",
263                                  progname);
264                  exit(1);
265          }
266 +        if (vtot <= 0)                  /* only create positive lobes */
267 +                return(0);
268 +                                        /* assign color */
269 +        if (rbf_colorimetry == RBCtristimulus) {
270 +                const double    df = 1.0 / (CIEXtot + vtot + CIEZtot);
271 +                C_COLOR         cclr;
272 +                c_cset(&cclr, CIEXtot*df, vtot*df);
273 +                rvp->chroma = c_encodeChroma(&cclr);
274 +        } else
275 +                rvp->chroma = c_dfchroma;
276                                          /* peak value based on integral */
277          vtot *= (x1-x0)*(y1-y0)*(2.*M_PI/GRIDRES/GRIDRES)/(double)nv;
278          rad = (RSCA/(double)GRIDRES)*(x1-x0);
279          rvp->peak =  vtot / ((2.*M_PI) * rad*rad);
280 <        rvp->crad = ANG2R(rad);
281 <        rvp->gx = (x0+x1)>>1;
282 <        rvp->gy = (y0+y1)>>1;
280 >        rvp->crad = ANG2R(rad);         /* put peak at centroid */
281 >        rvp->gx = (int)(wxsum/wtsum + .5);
282 >        rvp->gy = (int)(wysum/wtsum + .5);
283 >        return(1);
284   }
285  
286   /* Recursive function to build radial basis function representation */
# Line 234 | Line 322 | build_rbfrep(RBFVAL **arp, int *np, int x0, int x1, in
322                          return(-1);
323          }
324                                          /* create lobes for leaves */
325 <        if (!branched[0])
326 <                create_lobe(*arp+(*np)++, x0, xmid, y0, ymid);
327 <        if (!branched[1])
328 <                create_lobe(*arp+(*np)++, xmid, x1, y0, ymid);
329 <        if (!branched[2])
330 <                create_lobe(*arp+(*np)++, x0, xmid, ymid, y1);
331 <        if (!branched[3])
332 <                create_lobe(*arp+(*np)++, xmid, x1, ymid, y1);
333 <        nadded += nleaves;
325 >        if (!branched[0] && create_lobe(*arp+*np, x0, xmid, y0, ymid)) {
326 >                ++(*np); ++nadded;
327 >        }
328 >        if (!branched[1] && create_lobe(*arp+*np, xmid, x1, y0, ymid)) {
329 >                ++(*np); ++nadded;
330 >        }
331 >        if (!branched[2] && create_lobe(*arp+*np, x0, xmid, ymid, y1)) {
332 >                ++(*np); ++nadded;
333 >        }
334 >        if (!branched[3] && create_lobe(*arp+*np, xmid, x1, ymid, y1)) {
335 >                ++(*np); ++nadded;
336 >        }
337          return(nadded);
338   }
339  
# Line 257 | Line 348 | make_rbfrep()
348          comp_bsdf_min();
349                                  /* create RBF node list */
350          rbfarr = NULL; nn = 0;
351 <        if (build_rbfrep(&rbfarr, &nn, 0, GRIDRES, 0, GRIDRES) <= 0)
352 <                goto memerr;
351 >        if (build_rbfrep(&rbfarr, &nn, 0, GRIDRES, 0, GRIDRES) <= 0) {
352 >                if (nn)
353 >                        goto memerr;
354 >                fprintf(stderr,
355 >                        "%s: warning - skipping bad incidence (%.1f,%.1f)\n",
356 >                                progname, theta_in_deg, phi_in_deg);
357 >                return(NULL);
358 >        }
359                                  /* (re)allocate RBF array */
360          newnode = (RBFNODE *)realloc(rbfarr,
361                          sizeof(RBFNODE) + sizeof(RBFVAL)*(nn-1));
# Line 285 | Line 382 | make_rbfrep()
382                          newnode->vtotal);
383   #endif
384          insert_dsf(newnode);
288 –
385          return(newnode);
386   memerr:
387          fprintf(stderr, "%s: Out of memory in make_rbfrep()\n", progname);

Diff Legend

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