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.2 by greg, Tue Nov 13 04:23:38 2012 UTC vs.
Revision 2.9 by greg, Thu Oct 17 19:09:11 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.15            /* minimum radius scaling factor */
19   #endif
20 + #ifndef MAXRSCA
21 + #define MAXRSCA         2.7             /* maximum radius scaling factor */
22 + #endif
23 + #ifndef MAXFRAC
24 + #define MAXFRAC         0.5             /* maximum contribution to neighbor */
25 + #endif
26 + #ifndef NNEIGH
27 + #define NNEIGH          10              /* number of neighbors to consider */
28 + #endif
29                                  /* our loaded grid for this incident angle */
30   GRIDVAL                 dsf_grid[GRIDRES][GRIDRES];
31  
# Line 48 | Line 57 | add_bsdf_data(double theta_out, double phi_out, double
57          ovec[1] = sin((M_PI/180.)*phi_out) * ovec[2];
58          ovec[2] = sqrt(1. - ovec[2]*ovec[2]);
59  
60 <        if (!isDSF)
60 >        if (val <= 0)                   /* truncate to zero */
61 >                val = 0;
62 >        else if (!isDSF)
63                  val *= ovec[2];         /* convert from BSDF to DSF */
64  
65 +                                        /* update BSDF histogram */
66 +        if (val < BSDF2BIG*ovec[2] && val > BSDF2SML*ovec[2])
67 +                ++bsdf_hist[histndx(val/ovec[2])];
68 +
69          pos_from_vec(pos, ovec);
70  
71          dsf_grid[pos[0]][pos[1]].vsum += val;
# Line 58 | Line 73 | add_bsdf_data(double theta_out, double phi_out, double
73   }
74  
75   /* Compute radii for non-empty bins */
76 < /* (distance to furthest empty bin for which non-empty bin is the closest) */
76 > /* (distance to furthest empty bin for which non-empty test bin is closest) */
77   static void
78   compute_radii(void)
79   {
80 +        const int       cradmin = ANG2R(.5*M_PI/GRIDRES);
81          unsigned int    fill_grid[GRIDRES][GRIDRES];
82          unsigned short  fill_cnt[GRIDRES][GRIDRES];
83          FVECT           ovec0, ovec1;
84          double          ang2, lastang2;
85          int             r, i, j, jn, ii, jj, inear, jnear;
86  
87 +        for (i = 0; i < GRIDRES; i++)           /* initialize minimum radii */
88 +            for (j = 0; j < GRIDRES; j++)
89 +                if (dsf_grid[i][j].nval)
90 +                    dsf_grid[i][j].crad = cradmin;
91 +
92          r = GRIDRES/2;                          /* proceed in zig-zag */
93          for (i = 0; i < GRIDRES; i++)
94              for (jn = 0; jn < GRIDRES; jn++) {
# Line 111 | Line 132 | compute_radii(void)
132          memset(fill_cnt, 0, sizeof(fill_cnt));
133          for (i = 0; i < GRIDRES; i++)
134              for (j = 0; j < GRIDRES; j++) {
135 <                if (!dsf_grid[i][j].crad)
136 <                        continue;               /* missing distance */
137 <                r = R2ANG(dsf_grid[i][j].crad)*(2.*RSCA*GRIDRES/M_PI);
135 >                if (!dsf_grid[i][j].nval)
136 >                        continue;               /* not part of this */
137 >                r = R2ANG(dsf_grid[i][j].crad)*(2.*MAXRSCA*GRIDRES/M_PI);
138                  for (ii = i-r; ii <= i+r; ii++) {
139                      if (ii < 0) continue;
140                      if (ii >= GRIDRES) break;
# Line 158 | Line 179 | cull_values(void)
179                      if (ii < 0) continue;
180                      if (ii >= GRIDRES) break;
181                      for (jj = j-r; jj <= j+r; jj++) {
182 +                        if ((ii == i) & (jj == j))
183 +                                continue;       /* don't get self-absorbed */
184                          if (jj < 0) continue;
185                          if (jj >= GRIDRES) break;
186                          if (!dsf_grid[ii][jj].nval)
187                                  continue;
165 –                        if ((ii == i) & (jj == j))
166 –                                continue;       /* don't get self-absorbed */
188                          ovec_from_pos(ovec1, ii, jj);
189                          if (2. - 2.*DOT(ovec0,ovec1) >= maxang2)
190                                  continue;
# Line 185 | Line 206 | cull_values(void)
206                  }
207   }
208  
209 + /* Compute minimum BSDF from histogram and clear it */
210 + static void
211 + comp_bsdf_min()
212 + {
213 +        int     cnt;
214 +        int     i, target;
215 +
216 +        cnt = 0;
217 +        for (i = HISTLEN; i--; )
218 +                cnt += bsdf_hist[i];
219 +        if (!cnt) {                             /* shouldn't happen */
220 +                bsdf_min = 0;
221 +                return;
222 +        }
223 +        target = cnt/100;                       /* ignore bottom 1% */
224 +        cnt = 0;
225 +        for (i = 0; cnt <= target; i++)
226 +                cnt += bsdf_hist[i];
227 +        bsdf_min = histval(i-1);
228 +        memset(bsdf_hist, 0, sizeof(bsdf_hist));
229 + }
230 +
231 + /* Find n nearest sub-sampled neighbors to the given grid position */
232 + static int
233 + get_neighbors(int neigh[][2], int n, const int i, const int j)
234 + {
235 +        int     k = 0;
236 +        int     r;
237 +                                                /* search concentric squares */
238 +        for (r = 1; r < GRIDRES; r++) {
239 +                int     ii, jj;
240 +                for (ii = i-r; ii <= i+r; ii++) {
241 +                        int     jstep = 1;
242 +                        if (ii < 0) continue;
243 +                        if (ii >= GRIDRES) break;
244 +                        if ((i-r < ii) & (ii < i+r))
245 +                                jstep = r<<1;
246 +                        for (jj = j-r; jj <= j+r; jj += jstep) {
247 +                                if (jj < 0) continue;
248 +                                if (jj >= GRIDRES) break;
249 +                                if (dsf_grid[ii][jj].nval) {
250 +                                        neigh[k][0] = ii;
251 +                                        neigh[k][1] = jj;
252 +                                        if (++k >= n)
253 +                                                return(n);
254 +                                }
255 +                        }
256 +                }
257 +        }
258 +        return(k);
259 + }
260 +
261 + /* Adjust coded radius for the given grid position based on neighborhood */
262 + static int
263 + adj_coded_radius(const int i, const int j)
264 + {
265 +        const double    rad0 = R2ANG(dsf_grid[i][j].crad);
266 +        const double    minrad = MINRSCA * rad0;
267 +        double          currad = MAXRSCA * rad0;
268 +        int             neigh[NNEIGH][2];
269 +        int             n;
270 +        FVECT           our_dir;
271 +
272 +        ovec_from_pos(our_dir, i, j);
273 +        n = get_neighbors(neigh, NNEIGH, i, j);
274 +        while (n--) {
275 +                FVECT   their_dir;
276 +                double  max_ratio, rad_ok2;
277 +                                                /* check our value at neighbor */
278 +                ovec_from_pos(their_dir, neigh[n][0], neigh[n][1]);
279 +                max_ratio = MAXFRAC * dsf_grid[neigh[n][0]][neigh[n][1]].vsum
280 +                                / dsf_grid[i][j].vsum;
281 +                if (max_ratio >= 1)
282 +                        continue;
283 +                rad_ok2 = (DOT(their_dir,our_dir) - 1.)/log(max_ratio);
284 +                if (rad_ok2 >= currad*currad)
285 +                        continue;               /* value fraction OK */
286 +                currad = sqrt(rad_ok2);         /* else reduce lobe radius */
287 +                if (currad <= minrad)           /* limit how small we'll go */
288 +                        return(ANG2R(minrad));
289 +        }
290 +        return(ANG2R(currad));                  /* encode selected radius */
291 + }
292 +
293   /* Count up filled nodes and build RBF representation from current grid */
294   RBFNODE *
295   make_rbfrep(void)
296   {
297 +        long    cradsum = 0, ocradsum = 0;
298          int     niter = 16;
299          double  lastVar, thisVar = 100.;
300          int     nn;
301          RBFNODE *newnode;
302          RBFVAL  *itera;
303          int     i, j;
304 +
305 + #ifdef DEBUG
306 + {
307 +        int     maxcnt = 0, nempty = 0;
308 +        long    cntsum = 0;
309 +        for (i = 0; i < GRIDRES; i++)
310 +            for (j = 0; j < GRIDRES; j++)
311 +                if (!dsf_grid[i][j].nval) {
312 +                        ++nempty;
313 +                } else {
314 +                        if (dsf_grid[i][j].nval > maxcnt)
315 +                                maxcnt = dsf_grid[i][j].nval;
316 +                        cntsum += dsf_grid[i][j].nval;
317 +                }
318 +        fprintf(stderr, "Average, maximum bin count: %d, %d (%.1f%% empty)\n",
319 +                        (int)(cntsum/((GRIDRES*GRIDRES)-nempty)), maxcnt,
320 +                        100./(GRIDRES*GRIDRES)*nempty);
321 + }
322 + #endif
323                                  /* compute RBF radii */
324          compute_radii();
325                                  /* coagulate lobes */
# Line 203 | Line 328 | make_rbfrep(void)
328          for (i = 0; i < GRIDRES; i++)
329              for (j = 0; j < GRIDRES; j++)
330                  nn += dsf_grid[i][j].nval;
331 +                                /* compute minimum BSDF */
332 +        comp_bsdf_min();
333                                  /* allocate RBF array */
334          newnode = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(nn-1));
335          if (newnode == NULL)
# Line 221 | Line 348 | make_rbfrep(void)
348              for (j = 0; j < GRIDRES; j++)
349                  if (dsf_grid[i][j].nval) {
350                          newnode->rbfa[nn].peak = dsf_grid[i][j].vsum;
351 <                        newnode->rbfa[nn].crad = RSCA*dsf_grid[i][j].crad + .5;
351 >                        ocradsum += dsf_grid[i][j].crad;
352 >                        cradsum +=
353 >                        newnode->rbfa[nn].crad = adj_coded_radius(i, j);
354                          newnode->rbfa[nn].gx = i;
355                          newnode->rbfa[nn].gy = j;
356                          ++nn;
357                  }
358 + #ifdef DEBUG
359 +        fprintf(stderr,
360 +        "Average radius reduced from %.2f to %.2f degrees for %d lobes\n",
361 +                        180./M_PI*MAXRSCA*R2ANG(ocradsum/newnode->nrbf),
362 +                        180./M_PI*R2ANG(cradsum/newnode->nrbf), newnode->nrbf);
363 + #endif
364                                  /* iterate to improve interpolation accuracy */
365          itera = (RBFVAL *)malloc(sizeof(RBFVAL)*newnode->nrbf);
366          if (itera == NULL)
# Line 260 | Line 395 | make_rbfrep(void)
395          nn = 0;                 /* compute sum for normalization */
396          while (nn < newnode->nrbf)
397                  newnode->vtotal += rbf_volume(&newnode->rbfa[nn++]);
398 <
398 > #ifdef DEBUG
399 >        fprintf(stderr, "Integrated DSF at (%.1f,%.1f) deg. is %.2f\n",
400 >                        get_theta180(newnode->invec), get_phi360(newnode->invec),
401 >                        newnode->vtotal);
402 > #endif
403          insert_dsf(newnode);
404  
405          return(newnode);

Diff Legend

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