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

Comparing ray/src/cv/bsdfrbf.c (file contents):
Revision 2.1 by greg, Fri Oct 19 04:14:29 2012 UTC vs.
Revision 2.8 by greg, Wed Oct 2 20:38:26 2013 UTC

# Line 17 | Line 17 | static const char RCSid[] = "$Id$";
17   #ifndef RSCA
18   #define RSCA            2.7             /* radius scaling factor (empirical) */
19   #endif
20 + #ifndef MAXFRAC
21 + #define MAXFRAC         0.5             /* maximum contribution to neighbor */
22 + #endif
23 + #ifndef NNEIGH
24 + #define NNEIGH          10              /* number of neighbors to consider */
25 + #endif
26                                  /* our loaded grid for this incident angle */
27   GRIDVAL                 dsf_grid[GRIDRES][GRIDRES];
28  
# Line 48 | Line 54 | add_bsdf_data(double theta_out, double phi_out, double
54          ovec[1] = sin((M_PI/180.)*phi_out) * ovec[2];
55          ovec[2] = sqrt(1. - ovec[2]*ovec[2]);
56  
57 <        if (!isDSF)
57 >        if (val <= 0)                   /* truncate to zero */
58 >                val = 0;
59 >        else if (!isDSF)
60                  val *= ovec[2];         /* convert from BSDF to DSF */
61  
62 +                                        /* update BSDF histogram */
63 +        if (val < BSDF2BIG*ovec[2] && val > BSDF2SML*ovec[2])
64 +                ++bsdf_hist[histndx(val/ovec[2])];
65 +
66          pos_from_vec(pos, ovec);
67  
68          dsf_grid[pos[0]][pos[1]].vsum += val;
# Line 185 | Line 197 | cull_values(void)
197                  }
198   }
199  
200 + /* Compute minimum BSDF from histogram and clear it */
201 + static void
202 + comp_bsdf_min()
203 + {
204 +        int     cnt;
205 +        int     i, target;
206 +
207 +        cnt = 0;
208 +        for (i = HISTLEN; i--; )
209 +                cnt += bsdf_hist[i];
210 +        if (!cnt) {                             /* shouldn't happen */
211 +                bsdf_min = 0;
212 +                return;
213 +        }
214 +        target = cnt/100;                       /* ignore bottom 1% */
215 +        cnt = 0;
216 +        for (i = 0; cnt <= target; i++)
217 +                cnt += bsdf_hist[i];
218 +        bsdf_min = histval(i-1);
219 +        memset(bsdf_hist, 0, sizeof(bsdf_hist));
220 + }
221 +
222 + /* Find n nearest sub-sampled neighbors to the given grid position */
223 + static int
224 + get_neighbors(int neigh[][2], int n, const int i, const int j)
225 + {
226 +        int     k = 0;
227 +        int     r;
228 +                                                /* search concentric squares */
229 +        for (r = 1; r < GRIDRES; r++) {
230 +                int     ii, jj;
231 +                for (ii = i-r; ii <= i+r; ii++) {
232 +                        int     jstep = 1;
233 +                        if (ii < 0) continue;
234 +                        if (ii >= GRIDRES) break;
235 +                        if ((i-r < ii) & (ii < i+r))
236 +                                jstep = r<<1;
237 +                        for (jj = j-r; jj <= j+r; jj += jstep) {
238 +                                if (jj < 0) continue;
239 +                                if (jj >= GRIDRES) break;
240 +                                if (dsf_grid[ii][jj].nval) {
241 +                                        neigh[k][0] = ii;
242 +                                        neigh[k][1] = jj;
243 +                                        if (++k >= n)
244 +                                                return(n);
245 +                                }
246 +                        }
247 +                }
248 +        }
249 +        return(k);
250 + }
251 +
252 + /* Adjust coded radius for the given grid position based on neighborhood */
253 + static int
254 + adj_coded_radius(const int i, const int j)
255 + {
256 +        const double    rad0 = R2ANG(dsf_grid[i][j].crad);
257 +        double          currad = RSCA * rad0;
258 +        int             neigh[NNEIGH][2];
259 +        int             n;
260 +        FVECT           our_dir;
261 +
262 +        ovec_from_pos(our_dir, i, j);
263 +        n = get_neighbors(neigh, NNEIGH, i, j);
264 +        while (n--) {
265 +                FVECT   their_dir;
266 +                double  max_ratio, rad_ok2;
267 +                                                /* check our value at neighbor */
268 +                ovec_from_pos(their_dir, neigh[n][0], neigh[n][1]);
269 +                max_ratio = MAXFRAC * dsf_grid[neigh[n][0]][neigh[n][1]].vsum
270 +                                / dsf_grid[i][j].vsum;
271 +                if (max_ratio >= 1)
272 +                        continue;
273 +                rad_ok2 = (DOT(their_dir,our_dir) - 1.)/log(max_ratio);
274 +                if (rad_ok2 >= currad*currad)
275 +                        continue;               /* value fraction OK */
276 +                currad = sqrt(rad_ok2);         /* else reduce lobe radius */
277 +                if (currad <= rad0)             /* limit how small we'll go */
278 +                        return(dsf_grid[i][j].crad);
279 +        }
280 +        return(ANG2R(currad));                  /* encode selected radius */
281 + }
282 +
283   /* Count up filled nodes and build RBF representation from current grid */
284   RBFNODE *
285   make_rbfrep(void)
# Line 193 | Line 288 | make_rbfrep(void)
288          double  lastVar, thisVar = 100.;
289          int     nn;
290          RBFNODE *newnode;
291 +        RBFVAL  *itera;
292          int     i, j;
293                                  /* compute RBF radii */
294          compute_radii();
# Line 202 | Line 298 | make_rbfrep(void)
298          for (i = 0; i < GRIDRES; i++)
299              for (j = 0; j < GRIDRES; j++)
300                  nn += dsf_grid[i][j].nval;
301 +                                /* compute minimum BSDF */
302 +        comp_bsdf_min();
303                                  /* allocate RBF array */
304          newnode = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(nn-1));
305 <        if (newnode == NULL) {
306 <                fprintf(stderr, "%s: Out of memory in make_rbfrep()\n", progname);
209 <                exit(1);
210 <        }
305 >        if (newnode == NULL)
306 >                goto memerr;
307          newnode->ord = -1;
308          newnode->next = NULL;
309          newnode->ejl = NULL;
# Line 222 | Line 318 | make_rbfrep(void)
318              for (j = 0; j < GRIDRES; j++)
319                  if (dsf_grid[i][j].nval) {
320                          newnode->rbfa[nn].peak = dsf_grid[i][j].vsum;
321 <                        newnode->rbfa[nn].crad = RSCA*dsf_grid[i][j].crad + .5;
321 >                        newnode->rbfa[nn].crad = adj_coded_radius(i, j);
322                          newnode->rbfa[nn].gx = i;
323                          newnode->rbfa[nn].gy = j;
324                          ++nn;
325                  }
326                                  /* iterate to improve interpolation accuracy */
327 +        itera = (RBFVAL *)malloc(sizeof(RBFVAL)*newnode->nrbf);
328 +        if (itera == NULL)
329 +                goto memerr;
330 +        memcpy(itera, newnode->rbfa, sizeof(RBFVAL)*newnode->nrbf);
331          do {
332                  double  dsum = 0, dsum2 = 0;
333                  nn = 0;
# Line 237 | Line 337 | make_rbfrep(void)
337                                  FVECT   odir;
338                                  double  corr;
339                                  ovec_from_pos(odir, i, j);
340 <                                newnode->rbfa[nn++].peak *= corr =
340 >                                itera[nn++].peak *= corr =
341                                          dsf_grid[i][j].vsum /
342                                                  eval_rbfrep(newnode, odir);
343 <                                dsum += corr - 1.;
344 <                                dsum2 += (corr-1.)*(corr-1.);
343 >                                dsum += 1. - corr;
344 >                                dsum2 += (1.-corr)*(1.-corr);
345                          }
346 +                memcpy(newnode->rbfa, itera, sizeof(RBFVAL)*newnode->nrbf);
347                  lastVar = thisVar;
348                  thisVar = dsum2/(double)nn;
349   #ifdef DEBUG
# Line 252 | Line 353 | make_rbfrep(void)
353   #endif
354          } while (--niter > 0 && lastVar-thisVar > 0.02*lastVar);
355  
356 +        free(itera);
357          nn = 0;                 /* compute sum for normalization */
358          while (nn < newnode->nrbf)
359                  newnode->vtotal += rbf_volume(&newnode->rbfa[nn++]);
360 <
360 > #ifdef DEBUG
361 >        fprintf(stderr, "Integrated DSF at (%.1f,%.1f) deg. is %.2f\n",
362 >                        get_theta180(newnode->invec), get_phi360(newnode->invec),
363 >                        newnode->vtotal);
364 > #endif
365          insert_dsf(newnode);
366  
367          return(newnode);
368 + memerr:
369 +        fprintf(stderr, "%s: Out of memory in make_rbfrep()\n", progname);
370 +        exit(1);
371   }

Diff Legend

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