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

Comparing ray/src/cv/bsdfrbf.c (file contents):
Revision 2.8 by greg, Wed Oct 2 20:38:26 2013 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
# Line 69 | Line 75 | add_bsdf_data(double theta_out, double phi_out, double
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 118 | 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 162 | 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;
177 –                        if ((ii == i) & (jj == j))
178 –                                continue;       /* don't get self-absorbed */
228                          ovec_from_pos(ovec1, ii, jj);
229                          if (2. - 2.*DOT(ovec0,ovec1) >= maxang2)
230                                  continue;
# Line 254 | Line 303 | static int
303   adj_coded_radius(const int i, const int j)
304   {
305          const double    rad0 = R2ANG(dsf_grid[i][j].crad);
306 <        double          currad = RSCA * rad0;
306 >        const double    minrad = MINRSCA * rad0;
307 >        double          currad = MAXRSCA * rad0;
308          int             neigh[NNEIGH][2];
309          int             n;
310          FVECT           our_dir;
# Line 274 | Line 324 | adj_coded_radius(const int i, const int j)
324                  if (rad_ok2 >= currad*currad)
325                          continue;               /* value fraction OK */
326                  currad = sqrt(rad_ok2);         /* else reduce lobe radius */
327 <                if (currad <= rad0)             /* limit how small we'll go */
328 <                        return(dsf_grid[i][j].crad);
327 >                if (currad <= minrad)           /* limit how small we'll go */
328 >                        return(ANG2R(minrad));
329          }
330          return(ANG2R(currad));                  /* encode selected radius */
331   }
# Line 284 | Line 334 | adj_coded_radius(const int i, const int j)
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 318 | 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 +                        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)

Diff Legend

– Removed lines
+ Added lines
< Changed lines (old)
> Changed lines (new)