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.5 by greg, Fri Jun 28 23:18:51 2013 UTC vs.
Revision 2.11 by greg, Sat Oct 19 00:11:50 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         1.0             /* minimum radius scaling factor */
19   #endif
20 + #ifndef MAXRSCA
21 + #define MAXRSCA         2.7             /* maximum radius scaling factor */
22 + #endif
23 + #ifndef VARTHRESH
24 + #define VARTHRESH       0.0015          /* culling variance threshold */
25 + #endif
26 + #ifndef DIFFMAX2
27 + #define DIFFMAX2        (16.*VARTHRESH) /* maximum ignored sample variance */
28 + #endif
29 + #ifndef MAXFRAC
30 + #define MAXFRAC         0.5             /* maximum contribution to neighbor */
31 + #endif
32 + #ifndef NNEIGH
33 + #define NNEIGH          10              /* number of neighbors to consider */
34 + #endif
35                                  /* our loaded grid for this incident angle */
36   GRIDVAL                 dsf_grid[GRIDRES][GRIDRES];
37  
# Line 48 | Line 63 | add_bsdf_data(double theta_out, double phi_out, double
63          ovec[1] = sin((M_PI/180.)*phi_out) * ovec[2];
64          ovec[2] = sqrt(1. - ovec[2]*ovec[2]);
65  
66 <        if (!isDSF)
66 >        if (val <= 0)                   /* truncate to zero */
67 >                val = 0;
68 >        else if (!isDSF)
69                  val *= ovec[2];         /* convert from BSDF to DSF */
70  
71                                          /* update BSDF histogram */
# Line 62 | Line 79 | add_bsdf_data(double theta_out, double phi_out, double
79   }
80  
81   /* Compute radii for non-empty bins */
82 < /* (distance to furthest empty bin for which non-empty bin is the closest) */
82 > /* (distance to furthest empty bin for which non-empty test bin is closest) */
83   static void
84   compute_radii(void)
85   {
86 +        const int       cradmin = ANG2R(.5*M_PI/GRIDRES);
87          unsigned int    fill_grid[GRIDRES][GRIDRES];
88          unsigned short  fill_cnt[GRIDRES][GRIDRES];
89          FVECT           ovec0, ovec1;
# Line 110 | Line 128 | compute_radii(void)
128                                                  /* next search radius */
129                  r = ang2*(2.*GRIDRES/M_PI) + 3;
130              }
131 +        for (i = 0; i < GRIDRES; i++)           /* grow radii where uniform */
132 +            for (j = 0; j < GRIDRES; j++) {
133 +                double  midmean = 0.0;
134 +                int     nsum = 0;
135 +                if (!dsf_grid[i][j].nval)
136 +                        continue;
137 +                r = 1;                          /* avg. immediate neighbors */
138 +                for (ii = i-r; ii <= i+r; ii++) {
139 +                    if (ii < 0) continue;
140 +                    if (ii >= GRIDRES) break;
141 +                    for (jj = j-r; jj <= j+r; jj++) {
142 +                        if (jj < 0) continue;
143 +                        if (jj >= GRIDRES) break;
144 +                        midmean += dsf_grid[ii][jj].vsum;
145 +                        nsum += dsf_grid[ii][jj].nval;
146 +                    }
147 +                }
148 +                midmean /= (double)nsum;
149 +                while (++r < GRIDRES) {         /* attempt to grow perimeter */
150 +                    double      diff2sum = 0.0;
151 +                    nsum = 0;
152 +                    for (ii = i-r; ii <= i+r; ii++) {
153 +                        int     jstep = 1;
154 +                        if (ii < 0) continue;
155 +                        if (ii >= GRIDRES) break;
156 +                        if ((i-r < ii) & (ii < i+r))
157 +                                jstep = r<<1;
158 +                        for (jj = j-r; jj <= j+r; jj += jstep) {
159 +                                double  d2;
160 +                                if (jj < 0) continue;
161 +                                if (jj >= GRIDRES) break;
162 +                                if (!dsf_grid[ii][jj].nval)
163 +                                        continue;
164 +                                d2 = midmean - dsf_grid[ii][jj].vsum /
165 +                                                (double)dsf_grid[ii][jj].nval;
166 +                                d2 *= d2;
167 +                                if (d2 > DIFFMAX2*midmean*midmean)
168 +                                        goto escape;
169 +                                diff2sum += d2;
170 +                                ++nsum;
171 +                        }
172 +                    }
173 +                    if (diff2sum > VARTHRESH*midmean*midmean*(double)nsum)
174 +                        break;
175 +                }
176 + escape:         --r;
177 +                r = ANG2R(r*(M_PI/MAXRSCA/GRIDRES));
178 +                if (r < cradmin)
179 +                        r = cradmin;
180 +                if (dsf_grid[i][j].crad < r)
181 +                        dsf_grid[i][j].crad = r;
182 +            }
183                                                  /* blur radii over hemisphere */
184          memset(fill_grid, 0, sizeof(fill_grid));
185          memset(fill_cnt, 0, sizeof(fill_cnt));
186          for (i = 0; i < GRIDRES; i++)
187              for (j = 0; j < GRIDRES; j++) {
188 <                if (!dsf_grid[i][j].crad)
189 <                        continue;               /* missing distance */
190 <                r = R2ANG(dsf_grid[i][j].crad)*(2.*RSCA*GRIDRES/M_PI);
188 >                if (!dsf_grid[i][j].nval)
189 >                        continue;               /* not part of this */
190 >                r = R2ANG(dsf_grid[i][j].crad)*(2.*MAXRSCA*GRIDRES/M_PI);
191                  for (ii = i-r; ii <= i+r; ii++) {
192                      if (ii < 0) continue;
193                      if (ii >= GRIDRES) break;
# Line 138 | Line 208 | compute_radii(void)
208                          dsf_grid[i][j].crad = fill_grid[i][j]/fill_cnt[i][j];
209   }
210  
211 + /* Radius comparison for qsort() */
212 + static int
213 + radius_cmp(const void *p1, const void *p2)
214 + {
215 +        return( (int)dsf_grid[0][*(const int *)p1].crad -
216 +                (int)dsf_grid[0][*(const int *)p2].crad );
217 + }
218 +
219   /* Cull points for more uniform distribution, leave all nval 0 or 1 */
220   static void
221   cull_values(void)
222   {
223 +        int     indx[GRIDRES*GRIDRES];
224          FVECT   ovec0, ovec1;
225          double  maxang, maxang2;
226 <        int     i, j, ii, jj, r;
226 >        int     i, j, k, ii, jj, r;
227 >                                                /* sort by radius first */
228 >        for (k = GRIDRES*GRIDRES; k--; )
229 >                indx[k] = k;
230 >        qsort(indx, GRIDRES*GRIDRES, sizeof(int), &radius_cmp);
231                                                  /* simple greedy algorithm */
232 <        for (i = 0; i < GRIDRES; i++)
233 <            for (j = 0; j < GRIDRES; j++) {
232 >        for (k = GRIDRES*GRIDRES; k--; ) {
233 >                i = indx[k]/GRIDRES;            /* from biggest radius down */
234 >                j = indx[k] - i*GRIDRES;
235                  if (!dsf_grid[i][j].nval)
236                          continue;
237                  if (!dsf_grid[i][j].crad)
238 <                        continue;               /* shouldn't happen */
238 >                        break;          /* shouldn't happen */
239                  ovec_from_pos(ovec0, i, j);
240                  maxang = 2.*R2ANG(dsf_grid[i][j].crad);
241 <                if (maxang > ovec0[2])          /* clamp near horizon */
242 <                        maxang = ovec0[2];
241 >                                                /* clamp near horizon */
242 >                if (maxang > output_orient*ovec0[2])
243 >                        maxang = output_orient*ovec0[2];
244                  r = maxang*(2.*GRIDRES/M_PI) + 1;
245                  maxang2 = maxang*maxang;
246                  for (ii = i-r; ii <= i+r; ii++) {
247                      if (ii < 0) continue;
248                      if (ii >= GRIDRES) break;
249                      for (jj = j-r; jj <= j+r; jj++) {
250 +                        if ((ii == i) & (jj == j))
251 +                                continue;       /* don't get self-absorbed */
252                          if (jj < 0) continue;
253                          if (jj >= GRIDRES) break;
254                          if (!dsf_grid[ii][jj].nval)
255                                  continue;
169 –                        if ((ii == i) & (jj == j))
170 –                                continue;       /* don't get self-absorbed */
256                          ovec_from_pos(ovec1, ii, jj);
257                          if (2. - 2.*DOT(ovec0,ovec1) >= maxang2)
258                                  continue;
# Line 179 | Line 264 | cull_values(void)
264                          dsf_grid[ii][jj].nval = 0;
265                      }
266                  }
267 <            }
267 >        }
268                                                  /* final averaging pass */
269          for (i = 0; i < GRIDRES; i++)
270              for (j = 0; j < GRIDRES; j++)
# Line 189 | Line 274 | cull_values(void)
274                  }
275   }
276  
277 < /* Compute minimum BSDF from histogram and clear it */
277 > /* Compute minimum BSDF from histogram (does not clear) */
278   static void
279   comp_bsdf_min()
280   {
# Line 208 | Line 293 | comp_bsdf_min()
293          for (i = 0; cnt <= target; i++)
294                  cnt += bsdf_hist[i];
295          bsdf_min = histval(i-1);
211 –        memset(bsdf_hist, 0, sizeof(bsdf_hist));
296   }
297  
298 + /* Find n nearest sub-sampled neighbors to the given grid position */
299 + static int
300 + get_neighbors(int neigh[][2], int n, const int i, const int j)
301 + {
302 +        int     k = 0;
303 +        int     r;
304 +                                                /* search concentric squares */
305 +        for (r = 1; r < GRIDRES; r++) {
306 +                int     ii, jj;
307 +                for (ii = i-r; ii <= i+r; ii++) {
308 +                        int     jstep = 1;
309 +                        if (ii < 0) continue;
310 +                        if (ii >= GRIDRES) break;
311 +                        if ((i-r < ii) & (ii < i+r))
312 +                                jstep = r<<1;
313 +                        for (jj = j-r; jj <= j+r; jj += jstep) {
314 +                                if (jj < 0) continue;
315 +                                if (jj >= GRIDRES) break;
316 +                                if (dsf_grid[ii][jj].nval) {
317 +                                        neigh[k][0] = ii;
318 +                                        neigh[k][1] = jj;
319 +                                        if (++k >= n)
320 +                                                return(n);
321 +                                }
322 +                        }
323 +                }
324 +        }
325 +        return(k);
326 + }
327 +
328 + /* Adjust coded radius for the given grid position based on neighborhood */
329 + static int
330 + adj_coded_radius(const int i, const int j)
331 + {
332 +        const double    rad0 = R2ANG(dsf_grid[i][j].crad);
333 +        const double    minrad = MINRSCA * rad0;
334 +        double          currad = MAXRSCA * rad0;
335 +        int             neigh[NNEIGH][2];
336 +        int             n;
337 +        FVECT           our_dir;
338 +
339 +        ovec_from_pos(our_dir, i, j);
340 +        n = get_neighbors(neigh, NNEIGH, i, j);
341 +        while (n--) {
342 +                FVECT   their_dir;
343 +                double  max_ratio, rad_ok2;
344 +                                                /* check our value at neighbor */
345 +                ovec_from_pos(their_dir, neigh[n][0], neigh[n][1]);
346 +                max_ratio = MAXFRAC * dsf_grid[neigh[n][0]][neigh[n][1]].vsum
347 +                                / dsf_grid[i][j].vsum;
348 +                if (max_ratio >= 1)
349 +                        continue;
350 +                rad_ok2 = (DOT(their_dir,our_dir) - 1.)/log(max_ratio);
351 +                if (rad_ok2 >= currad*currad)
352 +                        continue;               /* value fraction OK */
353 +                currad = sqrt(rad_ok2);         /* else reduce lobe radius */
354 +                if (currad <= minrad)           /* limit how small we'll go */
355 +                        return(ANG2R(minrad));
356 +        }
357 +        return(ANG2R(currad));                  /* encode selected radius */
358 + }
359 +
360   /* Count up filled nodes and build RBF representation from current grid */
361   RBFNODE *
362   make_rbfrep(void)
363   {
364 +        long    cradsum = 0, ocradsum = 0;
365          int     niter = 16;
366          double  lastVar, thisVar = 100.;
367          int     nn;
368          RBFNODE *newnode;
369          RBFVAL  *itera;
370          int     i, j;
371 +
372 + #ifdef DEBUG
373 + {
374 +        int     maxcnt = 0, nempty = 0;
375 +        long    cntsum = 0;
376 +        for (i = 0; i < GRIDRES; i++)
377 +            for (j = 0; j < GRIDRES; j++)
378 +                if (!dsf_grid[i][j].nval) {
379 +                        ++nempty;
380 +                } else {
381 +                        if (dsf_grid[i][j].nval > maxcnt)
382 +                                maxcnt = dsf_grid[i][j].nval;
383 +                        cntsum += dsf_grid[i][j].nval;
384 +                }
385 +        fprintf(stderr, "Average, maximum bin count: %d, %d (%.1f%% empty)\n",
386 +                        (int)(cntsum/((GRIDRES*GRIDRES)-nempty)), maxcnt,
387 +                        100./(GRIDRES*GRIDRES)*nempty);
388 + }
389 + #endif
390                                  /* compute RBF radii */
391          compute_radii();
392                                  /* coagulate lobes */
# Line 249 | Line 415 | make_rbfrep(void)
415              for (j = 0; j < GRIDRES; j++)
416                  if (dsf_grid[i][j].nval) {
417                          newnode->rbfa[nn].peak = dsf_grid[i][j].vsum;
418 <                        newnode->rbfa[nn].crad = RSCA*dsf_grid[i][j].crad + .5;
418 >                        ocradsum += dsf_grid[i][j].crad;
419 >                        cradsum +=
420 >                        newnode->rbfa[nn].crad = adj_coded_radius(i, j);
421                          newnode->rbfa[nn].gx = i;
422                          newnode->rbfa[nn].gy = j;
423                          ++nn;
424                  }
425 + #ifdef DEBUG
426 +        fprintf(stderr,
427 +        "Average radius reduced from %.2f to %.2f degrees for %d lobes\n",
428 +                        180./M_PI*MAXRSCA*R2ANG(ocradsum/newnode->nrbf),
429 +                        180./M_PI*R2ANG(cradsum/newnode->nrbf), newnode->nrbf);
430 + #endif
431                                  /* iterate to improve interpolation accuracy */
432          itera = (RBFVAL *)malloc(sizeof(RBFVAL)*newnode->nrbf);
433          if (itera == NULL)

Diff Legend

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