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

Comparing ray/src/cv/pabopto2xml.c (file contents):
Revision 2.1 by greg, Fri Aug 24 15:20:18 2012 UTC vs.
Revision 2.2 by greg, Fri Aug 24 20:55:28 2012 UTC

# Line 20 | Line 20 | static const char RCSid[] = "$Id$";
20   #define GRIDRES         200             /* max. grid resolution per side */
21   #endif
22  
23 < #define RSCA            1.9             /* radius scaling factor (empirical) */
24 < #define MSCA            .12             /* magnitude scaling (empirical) */
23 > #define RSCA            3.              /* radius scaling factor (empirical) */
24 > #define MSCA            .2              /* magnitude scaling (empirical) */
25  
26 + #define R2ANG(c)        (((c)+.5)*(M_PI/(1<<16)))
27 + #define ANG2R(r)        (int)((r)*((1<<16)/M_PI))
28 +
29   typedef struct {
30          float           vsum;           /* BSDF sum */
31          unsigned short  nval;           /* number of values in sum */
32 <        unsigned short  hrad2;          /* half radius squared */
32 >        unsigned short  crad;           /* radius (coded angle) */
33   } GRIDVAL;                      /* grid value */
34  
35   typedef struct {
36 <        float           bsdf;           /* BSDF value at peak */
37 <        unsigned short  rad;            /* radius */
36 >        float           bsdf;           /* lobe value at peak */
37 >        unsigned short  crad;           /* radius (coded angle) */
38          unsigned char   gx, gy;         /* grid position */
39   } RBFVAL;                       /* radial basis function value */
40  
# Line 77 | Line 80 | make_rbfrep(void)
80                  if (bsdf_grid[i][j].nval) {
81                          newnode->rbfa[nn].bsdf = MSCA*bsdf_grid[i][j].vsum /
82                                                  (double)bsdf_grid[i][j].nval;
83 <                        newnode->rbfa[nn].rad =
81 <                                (int)(2.*RSCA*sqrt((double)bsdf_grid[i][j].hrad2) + .5);
83 >                        newnode->rbfa[nn].crad = RSCA*bsdf_grid[i][j].crad + .5;
84                          newnode->rbfa[nn].gx = i;
85                          newnode->rbfa[nn].gy = j;
86                          ++nn;
# Line 116 | Line 118 | vec_from_pos(FVECT vec, int xpos, int ypos)
118          vec[2] = 1. - r2;
119   }
120  
119 – /* Evaluate RBF at this grid position */
120 – static double
121 – eval_rbfrep2(const RBFLIST *rp, int xi, int yi)
122 – {
123 –        double          res = .0;
124 –        const RBFVAL    *rbfp;
125 –        double          sig2;
126 –        int             x2, y2;
127 –        int             n;
128 –
129 –        rbfp = rp->rbfa;
130 –        for (n = rp->nrbf; n--; rbfp++) {
131 –                x2 = (signed)rbfp->gx - xi;
132 –                x2 *= x2;
133 –                y2 = (signed)rbfp->gy - yi;
134 –                y2 *= y2;
135 –                sig2 = -.5*(x2 + y2)/(double)(rbfp->rad*rbfp->rad);
136 –                if (sig2 > -19.)
137 –                        res += rbfp->bsdf * exp(sig2);
138 –        }
139 –        return(res);
140 – }
141 –
121   /* Evaluate RBF for BSDF at the given normalized outgoing direction */
122   static double
123   eval_rbfrep(const RBFLIST *rp, const FVECT outvec)
# Line 152 | Line 131 | eval_rbfrep(const RBFLIST *rp, const FVECT outvec)
131          rbfp = rp->rbfa;
132          for (n = rp->nrbf; n--; rbfp++) {
133                  vec_from_pos(odir, rbfp->gx, rbfp->gy);
134 <                sig2 = (DOT(odir, outvec) - 1.) /
135 <                                ((M_PI*M_PI/(double)(GRIDRES*GRIDRES)) *
157 <                                                rbfp->rad*rbfp->rad);
134 >                sig2 = R2ANG(rbfp->crad);
135 >                sig2 = (DOT(odir,outvec) - 1.) / (sig2*sig2);
136                  if (sig2 > -19.)
137                          res += rbfp->bsdf * exp(sig2);
138          }
# Line 241 | Line 219 | load_bsdf_meas(const char *fname)
219   static void
220   compute_radii(void)
221   {
222 <        unsigned char   fill_grid[GRIDRES][GRIDRES];
223 <        int             r, r2, lastr2;
224 <        int             i, j, jn, ii, jj, inear, jnear;
225 <                                                /* proceed in zig-zag */
226 <        lastr2 = GRIDRES*GRIDRES;
222 >        unsigned short  fill_grid[GRIDRES][GRIDRES];
223 >        FVECT           ovec0, ovec1;
224 >        double          ang2, lastang2;
225 >        int             r2, lastr2;
226 >        int             r, i, j, jn, ii, jj, inear, jnear;
227 >
228 >        r = GRIDRES/2;                          /* proceed in zig-zag */
229          for (i = 0; i < GRIDRES; i++)
230              for (jn = 0; jn < GRIDRES; jn++) {
231                  j = (i&1) ? jn : GRIDRES-1-jn;
232                  if (bsdf_grid[i][j].nval)       /* find empty grid pos. */
233                          continue;
234 <                r = (int)sqrt((double)lastr2) + 2;
234 >                vec_from_pos(ovec0, i, j);
235                  inear = jnear = -1;             /* find nearest non-empty */
236 <                lastr2 = 2*GRIDRES*GRIDRES;
236 >                lastang2 = M_PI*M_PI;
237                  for (ii = i-r; ii <= i+r; ii++) {
238                      if (ii < 0) continue;
239                      if (ii >= GRIDRES) break;
# Line 262 | Line 242 | compute_radii(void)
242                          if (jj >= GRIDRES) break;
243                          if (!bsdf_grid[ii][jj].nval)
244                                  continue;
245 <                        r2 = (ii-i)*(ii-i) + (jj-j)*(jj-j);
246 <                        if (r2 >= lastr2)
245 >                        vec_from_pos(ovec1, ii, jj);
246 >                        ang2 = 2. - 2.*DOT(ovec0,ovec1);
247 >                        if (ang2 >= lastang2)
248                                  continue;
249 <                        lastr2 = r2;
249 >                        lastang2 = ang2;
250                          inear = ii; jnear = jj;
251                      }
252                  }
253 <                                                /* record if > previous */
254 <                if (bsdf_grid[inear][jnear].hrad2 < lastr2)
255 <                        bsdf_grid[inear][jnear].hrad2 = lastr2;
253 >                if (inear < 0) {
254 >                        fputs("Could not find non-empty neighbor!\n", stderr);
255 >                        exit(1);
256 >                }
257 >                ang2 = sqrt(lastang2);
258 >                r = ANG2R(ang2);                /* record if > previous */
259 >                if (r > bsdf_grid[inear][jnear].crad)
260 >                        bsdf_grid[inear][jnear].crad = r;
261 >                                                /* next search radius */
262 >                r = ang2*(2.*GRIDRES/M_PI) + 1;
263              }
264 <                                                /* fill in others */
264 >                                                /* fill in neighbors */
265          memset(fill_grid, 0, sizeof(fill_grid));
266          for (i = 0; i < GRIDRES; i++)
267              for (j = 0; j < GRIDRES; j++) {
268                  if (!bsdf_grid[i][j].nval)
269 <                        continue;
270 <                if (bsdf_grid[i][j].hrad2)
271 <                        continue;
269 >                        continue;               /* no value -- skip */
270 >                if (bsdf_grid[i][j].crad)
271 >                        continue;               /* has distance already */
272                  r = GRIDRES/20;
273 <                lastr2 = 2*r*r;
273 >                lastr2 = 2*r*r + 1;
274                  for (ii = i-r; ii <= i+r; ii++) {
275                      if (ii < 0) continue;
276                      if (ii >= GRIDRES) break;
277                      for (jj = j-r; jj <= j+r; jj++) {
278                          if (jj < 0) continue;
279                          if (jj >= GRIDRES) break;
280 <                        if (!bsdf_grid[ii][jj].hrad2)
280 >                        if (!bsdf_grid[ii][jj].crad)
281                                  continue;
282 +                                                /* OK to use approx. closest */
283                          r2 = (ii-i)*(ii-i) + (jj-j)*(jj-j);
284                          if (r2 >= lastr2)
285                                  continue;
286 <                        fill_grid[i][j] = bsdf_grid[ii][jj].hrad2;
286 >                        fill_grid[i][j] = bsdf_grid[ii][jj].crad;
287                          lastr2 = r2;
288                      }
289                  }
290              }
291 +                                                /* copy back filled entries */
292          for (i = 0; i < GRIDRES; i++)
293              for (j = 0; j < GRIDRES; j++)
294                  if (fill_grid[i][j])
295 <                        bsdf_grid[i][j].hrad2 = fill_grid[i][j];
295 >                        bsdf_grid[i][j].crad = fill_grid[i][j];
296   }
297  
298   /* Cull points for more uniform distribution */
299   static void
300   cull_values(void)
301   {
302 <        int     i, j, ii, jj, r, r2;
302 >        FVECT   ovec0, ovec1;
303 >        double  maxang, maxang2;
304 >        int     i, j, ii, jj, r;
305                                                  /* simple greedy algorithm */
306          for (i = 0; i < GRIDRES; i++)
307              for (j = 0; j < GRIDRES; j++) {
308                  if (!bsdf_grid[i][j].nval)
309                          continue;
310 <                if (!bsdf_grid[i][j].hrad2)
311 <                        continue;
312 <                r = (int)(2.*sqrt((double)bsdf_grid[i][j].hrad2) + .9999);
310 >                if (!bsdf_grid[i][j].crad)
311 >                        continue;               /* shouldn't happen */
312 >                vec_from_pos(ovec0, i, j);
313 >                maxang = 2.*R2ANG(bsdf_grid[i][j].crad);
314 >                if (maxang > ovec0[2])          /* clamp near horizon */
315 >                        maxang = ovec0[2];
316 >                r = maxang*(2.*GRIDRES/M_PI) + 1;
317 >                maxang2 = maxang*maxang;
318                  for (ii = i-r; ii <= i+r; ii++) {
319                      if (ii < 0) continue;
320                      if (ii >= GRIDRES) break;
# Line 326 | Line 323 | cull_values(void)
323                          if (jj >= GRIDRES) break;
324                          if (!bsdf_grid[ii][jj].nval)
325                                  continue;
326 <                        r2 = (ii-i)*(ii-i) + (jj-j)*(jj-j);
327 <                        if (!r2 | (r2 > r*r))
326 >                        if ((ii == i) & (jj == j))
327 >                                continue;       /* don't get self-absorbed */
328 >                        vec_from_pos(ovec1, ii, jj);
329 >                        if (2. - 2.*DOT(ovec0,ovec1) >= maxang2)
330                                  continue;
331 <                                                /* absorb victim's value */
331 >                                                /* absorb sum */
332                          bsdf_grid[i][j].vsum += bsdf_grid[ii][jj].vsum;
333                          bsdf_grid[i][j].nval += bsdf_grid[ii][jj].nval;
334 <                        memset(&bsdf_grid[ii][jj], 0, sizeof(GRIDVAL));
334 >                                                /* keep value, though */
335 >                        bsdf_grid[ii][jj].vsum /= (double)bsdf_grid[ii][jj].nval;
336 >                        bsdf_grid[ii][jj].nval = 0;
337                      }
338                  }
339              }

Diff Legend

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