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.10 by greg, Fri Oct 18 02:49:30 2013 UTC

# Line 14 | Line 14 | static const char RCSid[] = "$Id$";
14   #include <math.h>
15   #include "bsdfrep.h"
16  
17 < #ifndef RSCA
18 < #define RSCA            2.7             /* radius scaling factor (empirical) */
17 > #ifndef MINRSCA
18 > #define MINRSCA         0.5             /* minimum radius scaling factor */
19   #endif
20 + #ifndef MAXRSCA
21 + #define MAXRSCA         2.7             /* maximum radius scaling factor */
22 + #endif
23 + #ifndef DIFFTHRESH
24 + #define DIFFTHRESH      0.2             /* culling difference threshold */
25 + #endif
26 + #ifndef MAXFRAC
27 + #define MAXFRAC         0.5             /* maximum contribution to neighbor */
28 + #endif
29 + #ifndef NNEIGH
30 + #define NNEIGH          10              /* number of neighbors to consider */
31 + #endif
32                                  /* our loaded grid for this incident angle */
33   GRIDVAL                 dsf_grid[GRIDRES][GRIDRES];
34  
# Line 48 | Line 60 | add_bsdf_data(double theta_out, double phi_out, double
60          ovec[1] = sin((M_PI/180.)*phi_out) * ovec[2];
61          ovec[2] = sqrt(1. - ovec[2]*ovec[2]);
62  
63 <        if (!isDSF)
63 >        if (val <= 0)                   /* truncate to zero */
64 >                val = 0;
65 >        else if (!isDSF)
66                  val *= ovec[2];         /* convert from BSDF to DSF */
67  
68 +                                        /* update BSDF histogram */
69 +        if (val < BSDF2BIG*ovec[2] && val > BSDF2SML*ovec[2])
70 +                ++bsdf_hist[histndx(val/ovec[2])];
71 +
72          pos_from_vec(pos, ovec);
73  
74          dsf_grid[pos[0]][pos[1]].vsum += val;
75          dsf_grid[pos[0]][pos[1]].nval++;
76   }
77  
78 + /* Check if the two DSF values are significantly different */
79 + static int
80 + big_diff(double ref, double tst)
81 + {
82 +        if (ref > 0) {
83 +                tst = tst/ref - 1.;
84 +                if (tst < 0) tst = -tst;
85 +        } else
86 +                tst *= 50.;
87 +        return(tst > DIFFTHRESH);
88 + }
89 +
90   /* Compute radii for non-empty bins */
91 < /* (distance to furthest empty bin for which non-empty bin is the closest) */
91 > /* (distance to furthest empty bin for which non-empty test bin is closest) */
92   static void
93   compute_radii(void)
94   {
95 +        const int       cradmin = ANG2R(.5*M_PI/GRIDRES);
96          unsigned int    fill_grid[GRIDRES][GRIDRES];
97          unsigned short  fill_cnt[GRIDRES][GRIDRES];
98          FVECT           ovec0, ovec1;
# Line 106 | Line 137 | compute_radii(void)
137                                                  /* next search radius */
138                  r = ang2*(2.*GRIDRES/M_PI) + 3;
139              }
140 +        for (i = 0; i < GRIDRES; i++)           /* grow radii where uniform */
141 +            for (j = 0; j < GRIDRES; j++) {
142 +                double  midmean;
143 +                if (!dsf_grid[i][j].nval)
144 +                        continue;
145 +                midmean = dsf_grid[i][j].vsum / (double)dsf_grid[i][j].nval;
146 +                r = R2ANG(dsf_grid[i][j].crad)*(MAXRSCA*GRIDRES/M_PI);
147 +                while (++r < GRIDRES) {         /* attempt to grow perimeter */
148 +                    for (ii = i-r; ii <= i+r; ii++) {
149 +                        int     jstep = 1;
150 +                        if (ii < 0) continue;
151 +                        if (ii >= GRIDRES) break;
152 +                        if ((i-r < ii) & (ii < i+r))
153 +                                jstep = r<<1;
154 +                        for (jj = j-r; jj <= j+r; jj += jstep) {
155 +                                if (jj < 0) continue;
156 +                                if (jj >= GRIDRES) break;
157 +                                if (dsf_grid[ii][jj].nval && big_diff(midmean,
158 +                                                dsf_grid[ii][jj].vsum /
159 +                                                (double)dsf_grid[ii][jj].nval))
160 +                                        goto hit_diff;
161 +                        }
162 +                    }
163 +                }
164 + hit_diff:       --r;
165 +                dsf_grid[i][j].crad = ANG2R(r*(M_PI/MAXRSCA/GRIDRES));
166 +                if (dsf_grid[i][j].crad < cradmin)
167 +                        dsf_grid[i][j].crad = cradmin;
168 +            }
169                                                  /* blur radii over hemisphere */
170          memset(fill_grid, 0, sizeof(fill_grid));
171          memset(fill_cnt, 0, sizeof(fill_cnt));
172          for (i = 0; i < GRIDRES; i++)
173              for (j = 0; j < GRIDRES; j++) {
174 <                if (!dsf_grid[i][j].crad)
175 <                        continue;               /* missing distance */
176 <                r = R2ANG(dsf_grid[i][j].crad)*(2.*RSCA*GRIDRES/M_PI);
174 >                if (!dsf_grid[i][j].nval)
175 >                        continue;               /* not part of this */
176 >                r = R2ANG(dsf_grid[i][j].crad)*(2.*MAXRSCA*GRIDRES/M_PI);
177                  for (ii = i-r; ii <= i+r; ii++) {
178                      if (ii < 0) continue;
179                      if (ii >= GRIDRES) break;
# Line 150 | Line 210 | cull_values(void)
210                          continue;               /* shouldn't happen */
211                  ovec_from_pos(ovec0, i, j);
212                  maxang = 2.*R2ANG(dsf_grid[i][j].crad);
213 <                if (maxang > ovec0[2])          /* clamp near horizon */
214 <                        maxang = ovec0[2];
213 >                                                /* clamp near horizon */
214 >                if (maxang > output_orient*ovec0[2])
215 >                        maxang = output_orient*ovec0[2];
216                  r = maxang*(2.*GRIDRES/M_PI) + 1;
217                  maxang2 = maxang*maxang;
218                  for (ii = i-r; ii <= i+r; ii++) {
219                      if (ii < 0) continue;
220                      if (ii >= GRIDRES) break;
221                      for (jj = j-r; jj <= j+r; jj++) {
222 +                        if ((ii == i) & (jj == j))
223 +                                continue;       /* don't get self-absorbed */
224                          if (jj < 0) continue;
225                          if (jj >= GRIDRES) break;
226                          if (!dsf_grid[ii][jj].nval)
227                                  continue;
165 –                        if ((ii == i) & (jj == j))
166 –                                continue;       /* don't get self-absorbed */
228                          ovec_from_pos(ovec1, ii, jj);
229                          if (2. - 2.*DOT(ovec0,ovec1) >= maxang2)
230                                  continue;
# Line 185 | Line 246 | cull_values(void)
246                  }
247   }
248  
249 + /* Compute minimum BSDF from histogram and clear it */
250 + static void
251 + comp_bsdf_min()
252 + {
253 +        int     cnt;
254 +        int     i, target;
255 +
256 +        cnt = 0;
257 +        for (i = HISTLEN; i--; )
258 +                cnt += bsdf_hist[i];
259 +        if (!cnt) {                             /* shouldn't happen */
260 +                bsdf_min = 0;
261 +                return;
262 +        }
263 +        target = cnt/100;                       /* ignore bottom 1% */
264 +        cnt = 0;
265 +        for (i = 0; cnt <= target; i++)
266 +                cnt += bsdf_hist[i];
267 +        bsdf_min = histval(i-1);
268 +        memset(bsdf_hist, 0, sizeof(bsdf_hist));
269 + }
270 +
271 + /* Find n nearest sub-sampled neighbors to the given grid position */
272 + static int
273 + get_neighbors(int neigh[][2], int n, const int i, const int j)
274 + {
275 +        int     k = 0;
276 +        int     r;
277 +                                                /* search concentric squares */
278 +        for (r = 1; r < GRIDRES; r++) {
279 +                int     ii, jj;
280 +                for (ii = i-r; ii <= i+r; ii++) {
281 +                        int     jstep = 1;
282 +                        if (ii < 0) continue;
283 +                        if (ii >= GRIDRES) break;
284 +                        if ((i-r < ii) & (ii < i+r))
285 +                                jstep = r<<1;
286 +                        for (jj = j-r; jj <= j+r; jj += jstep) {
287 +                                if (jj < 0) continue;
288 +                                if (jj >= GRIDRES) break;
289 +                                if (dsf_grid[ii][jj].nval) {
290 +                                        neigh[k][0] = ii;
291 +                                        neigh[k][1] = jj;
292 +                                        if (++k >= n)
293 +                                                return(n);
294 +                                }
295 +                        }
296 +                }
297 +        }
298 +        return(k);
299 + }
300 +
301 + /* Adjust coded radius for the given grid position based on neighborhood */
302 + static int
303 + adj_coded_radius(const int i, const int j)
304 + {
305 +        const double    rad0 = R2ANG(dsf_grid[i][j].crad);
306 +        const double    minrad = MINRSCA * rad0;
307 +        double          currad = MAXRSCA * rad0;
308 +        int             neigh[NNEIGH][2];
309 +        int             n;
310 +        FVECT           our_dir;
311 +
312 +        ovec_from_pos(our_dir, i, j);
313 +        n = get_neighbors(neigh, NNEIGH, i, j);
314 +        while (n--) {
315 +                FVECT   their_dir;
316 +                double  max_ratio, rad_ok2;
317 +                                                /* check our value at neighbor */
318 +                ovec_from_pos(their_dir, neigh[n][0], neigh[n][1]);
319 +                max_ratio = MAXFRAC * dsf_grid[neigh[n][0]][neigh[n][1]].vsum
320 +                                / dsf_grid[i][j].vsum;
321 +                if (max_ratio >= 1)
322 +                        continue;
323 +                rad_ok2 = (DOT(their_dir,our_dir) - 1.)/log(max_ratio);
324 +                if (rad_ok2 >= currad*currad)
325 +                        continue;               /* value fraction OK */
326 +                currad = sqrt(rad_ok2);         /* else reduce lobe radius */
327 +                if (currad <= minrad)           /* limit how small we'll go */
328 +                        return(ANG2R(minrad));
329 +        }
330 +        return(ANG2R(currad));                  /* encode selected radius */
331 + }
332 +
333   /* Count up filled nodes and build RBF representation from current grid */
334   RBFNODE *
335   make_rbfrep(void)
336   {
337 +        long    cradsum = 0, ocradsum = 0;
338          int     niter = 16;
339          double  lastVar, thisVar = 100.;
340          int     nn;
341          RBFNODE *newnode;
342 +        RBFVAL  *itera;
343          int     i, j;
344 +
345 + #ifdef DEBUG
346 + {
347 +        int     maxcnt = 0, nempty = 0;
348 +        long    cntsum = 0;
349 +        for (i = 0; i < GRIDRES; i++)
350 +            for (j = 0; j < GRIDRES; j++)
351 +                if (!dsf_grid[i][j].nval) {
352 +                        ++nempty;
353 +                } else {
354 +                        if (dsf_grid[i][j].nval > maxcnt)
355 +                                maxcnt = dsf_grid[i][j].nval;
356 +                        cntsum += dsf_grid[i][j].nval;
357 +                }
358 +        fprintf(stderr, "Average, maximum bin count: %d, %d (%.1f%% empty)\n",
359 +                        (int)(cntsum/((GRIDRES*GRIDRES)-nempty)), maxcnt,
360 +                        100./(GRIDRES*GRIDRES)*nempty);
361 + }
362 + #endif
363                                  /* compute RBF radii */
364          compute_radii();
365                                  /* coagulate lobes */
# Line 202 | Line 368 | make_rbfrep(void)
368          for (i = 0; i < GRIDRES; i++)
369              for (j = 0; j < GRIDRES; j++)
370                  nn += dsf_grid[i][j].nval;
371 +                                /* compute minimum BSDF */
372 +        comp_bsdf_min();
373                                  /* allocate RBF array */
374          newnode = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(nn-1));
375 <        if (newnode == NULL) {
376 <                fprintf(stderr, "%s: Out of memory in make_rbfrep()\n", progname);
209 <                exit(1);
210 <        }
375 >        if (newnode == NULL)
376 >                goto memerr;
377          newnode->ord = -1;
378          newnode->next = NULL;
379          newnode->ejl = NULL;
# Line 222 | Line 388 | make_rbfrep(void)
388              for (j = 0; j < GRIDRES; j++)
389                  if (dsf_grid[i][j].nval) {
390                          newnode->rbfa[nn].peak = dsf_grid[i][j].vsum;
391 <                        newnode->rbfa[nn].crad = RSCA*dsf_grid[i][j].crad + .5;
391 >                        ocradsum += dsf_grid[i][j].crad;
392 >                        cradsum +=
393 >                        newnode->rbfa[nn].crad = adj_coded_radius(i, j);
394                          newnode->rbfa[nn].gx = i;
395                          newnode->rbfa[nn].gy = j;
396                          ++nn;
397                  }
398 + #ifdef DEBUG
399 +        fprintf(stderr,
400 +        "Average radius reduced from %.2f to %.2f degrees for %d lobes\n",
401 +                        180./M_PI*MAXRSCA*R2ANG(ocradsum/newnode->nrbf),
402 +                        180./M_PI*R2ANG(cradsum/newnode->nrbf), newnode->nrbf);
403 + #endif
404                                  /* iterate to improve interpolation accuracy */
405 +        itera = (RBFVAL *)malloc(sizeof(RBFVAL)*newnode->nrbf);
406 +        if (itera == NULL)
407 +                goto memerr;
408 +        memcpy(itera, newnode->rbfa, sizeof(RBFVAL)*newnode->nrbf);
409          do {
410                  double  dsum = 0, dsum2 = 0;
411                  nn = 0;
# Line 237 | Line 415 | make_rbfrep(void)
415                                  FVECT   odir;
416                                  double  corr;
417                                  ovec_from_pos(odir, i, j);
418 <                                newnode->rbfa[nn++].peak *= corr =
418 >                                itera[nn++].peak *= corr =
419                                          dsf_grid[i][j].vsum /
420                                                  eval_rbfrep(newnode, odir);
421 <                                dsum += corr - 1.;
422 <                                dsum2 += (corr-1.)*(corr-1.);
421 >                                dsum += 1. - corr;
422 >                                dsum2 += (1.-corr)*(1.-corr);
423                          }
424 +                memcpy(newnode->rbfa, itera, sizeof(RBFVAL)*newnode->nrbf);
425                  lastVar = thisVar;
426                  thisVar = dsum2/(double)nn;
427   #ifdef DEBUG
# Line 252 | Line 431 | make_rbfrep(void)
431   #endif
432          } while (--niter > 0 && lastVar-thisVar > 0.02*lastVar);
433  
434 +        free(itera);
435          nn = 0;                 /* compute sum for normalization */
436          while (nn < newnode->nrbf)
437                  newnode->vtotal += rbf_volume(&newnode->rbfa[nn++]);
438 <
438 > #ifdef DEBUG
439 >        fprintf(stderr, "Integrated DSF at (%.1f,%.1f) deg. is %.2f\n",
440 >                        get_theta180(newnode->invec), get_phi360(newnode->invec),
441 >                        newnode->vtotal);
442 > #endif
443          insert_dsf(newnode);
444  
445          return(newnode);
446 + memerr:
447 +        fprintf(stderr, "%s: Out of memory in make_rbfrep()\n", progname);
448 +        exit(1);
449   }

Diff Legend

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