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.5 by greg, Sat Aug 25 22:39:03 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            2.7             /* radius scaling factor (empirical) */
24  
25 + #define R2ANG(c)        (((c)+.5)*(M_PI/(1<<16)))
26 + #define ANG2R(r)        (int)((r)*((1<<16)/M_PI))
27 +
28   typedef struct {
29 <        float           vsum;           /* BSDF sum */
29 >        float           vsum;           /* DSF sum */
30          unsigned short  nval;           /* number of values in sum */
31 <        unsigned short  hrad2;          /* half radius squared */
31 >        unsigned short  crad;           /* radius (coded angle) */
32   } GRIDVAL;                      /* grid value */
33  
34   typedef struct {
35 <        float           bsdf;           /* BSDF value at peak */
36 <        unsigned short  rad;            /* radius */
35 >        float           peak;           /* lobe value at peak */
36 >        unsigned short  crad;           /* radius (coded angle) */
37          unsigned char   gx, gy;         /* grid position */
38   } RBFVAL;                       /* radial basis function value */
39  
# Line 40 | Line 42 | typedef struct s_rbflist {
42          FVECT                   invec;          /* incident vector direction */
43          int                     nrbf;           /* number of RBFs */
44          RBFVAL                  rbfa[1];        /* RBF array (extends struct) */
45 < } RBFLIST;                      /* RBF representation of BSDF @ 1 incidence */
45 > } RBFLIST;                      /* RBF representation of DSF @ 1 incidence */
46  
47                                  /* our loaded grid for this incident angle */
48   static double   theta_in_deg, phi_in_deg;
49 < static GRIDVAL  bsdf_grid[GRIDRES][GRIDRES];
49 > static GRIDVAL  dsf_grid[GRIDRES][GRIDRES];
50  
51 <                                /* processed incident BSDF measurements */
52 < static RBFLIST  *bsdf_list = NULL;
51 >                                /* processed incident DSF measurements */
52 > static RBFLIST  *dsf_list = NULL;
53  
52 – /* Count up non-empty nodes and build RBF representation from current grid */
53 – static RBFLIST *
54 – make_rbfrep(void)
55 – {
56 –        int     nn = 0;
57 –        RBFLIST *newnode;
58 –        int     i, j;
59 –                                /* count non-empty bins */
60 –        for (i = 0; i < GRIDRES; i++)
61 –            for (j = 0; j < GRIDRES; j++)
62 –                nn += (bsdf_grid[i][j].nval > 0);
63 –                                /* allocate RBF array */
64 –        newnode = (RBFLIST *)malloc(sizeof(RBFLIST) + sizeof(RBFVAL)*(nn-1));
65 –        if (newnode == NULL) {
66 –                fputs("Out of memory in make_rbfrep\n", stderr);
67 –                exit(1);
68 –        }
69 –        newnode->invec[2] = sin(M_PI/180.*theta_in_deg);
70 –        newnode->invec[0] = cos(M_PI/180.*phi_in_deg)*newnode->invec[2];
71 –        newnode->invec[1] = sin(M_PI/180.*phi_in_deg)*newnode->invec[2];
72 –        newnode->invec[2] = sqrt(1. - newnode->invec[2]*newnode->invec[2]);
73 –        newnode->nrbf = nn;
74 –        nn = 0;                 /* fill RBF array */
75 –        for (i = 0; i < GRIDRES; i++)
76 –            for (j = 0; j < GRIDRES; j++)
77 –                if (bsdf_grid[i][j].nval) {
78 –                        newnode->rbfa[nn].bsdf = MSCA*bsdf_grid[i][j].vsum /
79 –                                                (double)bsdf_grid[i][j].nval;
80 –                        newnode->rbfa[nn].rad =
81 –                                (int)(2.*RSCA*sqrt((double)bsdf_grid[i][j].hrad2) + .5);
82 –                        newnode->rbfa[nn].gx = i;
83 –                        newnode->rbfa[nn].gy = j;
84 –                        ++nn;
85 –                }
86 –        newnode->next = bsdf_list;
87 –        return(bsdf_list = newnode);
88 – }
89 –
90 – /* Compute grid position from normalized outgoing vector */
91 – static void
92 – pos_from_vec(int pos[2], const FVECT vec)
93 – {
94 –        double  sq[2];          /* uniform hemispherical projection */
95 –        double  norm = 1./sqrt(1. + vec[2]);
96 –
97 –        SDdisk2square(sq, vec[0]*norm, vec[1]*norm);
98 –
99 –        pos[0] = (int)(sq[0]*GRIDRES);
100 –        pos[1] = (int)(sq[1]*GRIDRES);
101 – }
102 –
54   /* Compute outgoing vector from grid position */
55   static void
56   vec_from_pos(FVECT vec, int xpos, int ypos)
# Line 116 | Line 67 | vec_from_pos(FVECT vec, int xpos, int ypos)
67          vec[2] = 1. - r2;
68   }
69  
70 < /* Evaluate RBF at this grid position */
71 < static double
72 < eval_rbfrep2(const RBFLIST *rp, int xi, int yi)
70 > /* Compute grid position from normalized outgoing vector */
71 > static void
72 > pos_from_vec(int pos[2], const FVECT vec)
73   {
74 <        double          res = .0;
75 <        const RBFVAL    *rbfp;
125 <        double          sig2;
126 <        int             x2, y2;
127 <        int             n;
74 >        double  sq[2];          /* uniform hemispherical projection */
75 >        double  norm = 1./sqrt(1. + vec[2]);
76  
77 <        rbfp = rp->rbfa;
78 <        for (n = rp->nrbf; n--; rbfp++) {
79 <                x2 = (signed)rbfp->gx - xi;
80 <                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);
77 >        SDdisk2square(sq, vec[0]*norm, vec[1]*norm);
78 >
79 >        pos[0] = (int)(sq[0]*GRIDRES);
80 >        pos[1] = (int)(sq[1]*GRIDRES);
81   }
82  
83 < /* Evaluate RBF for BSDF at the given normalized outgoing direction */
83 > /* Evaluate RBF for DSF at the given normalized outgoing direction */
84   static double
85   eval_rbfrep(const RBFLIST *rp, const FVECT outvec)
86   {
# Line 152 | Line 93 | eval_rbfrep(const RBFLIST *rp, const FVECT outvec)
93          rbfp = rp->rbfa;
94          for (n = rp->nrbf; n--; rbfp++) {
95                  vec_from_pos(odir, rbfp->gx, rbfp->gy);
96 <                sig2 = (DOT(odir, outvec) - 1.) /
97 <                                ((M_PI*M_PI/(double)(GRIDRES*GRIDRES)) *
157 <                                                rbfp->rad*rbfp->rad);
96 >                sig2 = R2ANG(rbfp->crad);
97 >                sig2 = (DOT(odir,outvec) - 1.) / (sig2*sig2);
98                  if (sig2 > -19.)
99 <                        res += rbfp->bsdf * exp(sig2);
99 >                        res += rbfp->peak * exp(sig2);
100          }
101          return(res);
102   }
103  
104 + /* Count up filled nodes and build RBF representation from current grid */
105 + static RBFLIST *
106 + make_rbfrep(void)
107 + {
108 +        int     niter = 6;
109 +        int     nn;
110 +        RBFLIST *newnode;
111 +        int     i, j;
112 +        
113 +        nn = 0;                 /* count selected bins */
114 +        for (i = 0; i < GRIDRES; i++)
115 +            for (j = 0; j < GRIDRES; j++)
116 +                nn += (dsf_grid[i][j].nval > 0);
117 +                                /* allocate RBF array */
118 +        newnode = (RBFLIST *)malloc(sizeof(RBFLIST) + sizeof(RBFVAL)*(nn-1));
119 +        if (newnode == NULL) {
120 +                fputs("Out of memory in make_rbfrep\n", stderr);
121 +                exit(1);
122 +        }
123 +        newnode->next = NULL;
124 +        newnode->invec[2] = sin(M_PI/180.*theta_in_deg);
125 +        newnode->invec[0] = cos(M_PI/180.*phi_in_deg)*newnode->invec[2];
126 +        newnode->invec[1] = sin(M_PI/180.*phi_in_deg)*newnode->invec[2];
127 +        newnode->invec[2] = sqrt(1. - newnode->invec[2]*newnode->invec[2]);
128 +        newnode->nrbf = nn;
129 +        nn = 0;                 /* fill RBF array */
130 +        for (i = 0; i < GRIDRES; i++)
131 +            for (j = 0; j < GRIDRES; j++)
132 +                if (dsf_grid[i][j].nval) {
133 +                        newnode->rbfa[nn].peak =
134 +                                        dsf_grid[i][j].vsum /=
135 +                                                (double)dsf_grid[i][j].nval;
136 +                        dsf_grid[i][j].nval = 1;
137 +                        newnode->rbfa[nn].crad = RSCA*dsf_grid[i][j].crad + .5;
138 +                        newnode->rbfa[nn].gx = i;
139 +                        newnode->rbfa[nn].gy = j;
140 +                        ++nn;
141 +                }
142 +                                /* iterate for better convergence */
143 +        while (niter--) {
144 +                double  dsum = .0, dsum2 = .0;
145 +                nn = 0;
146 +                for (i = 0; i < GRIDRES; i++)
147 +                    for (j = 0; j < GRIDRES; j++)
148 +                        if (dsf_grid[i][j].nval) {
149 +                                FVECT   odir;
150 +                                /* double       corr; */
151 +                                vec_from_pos(odir, i, j);
152 +                                newnode->rbfa[nn++].peak *= /* corr = */
153 +                                        dsf_grid[i][j].vsum /
154 +                                                eval_rbfrep(newnode, odir);
155 +                                /*
156 +                                dsum += corr - 1.;
157 +                                dsum2 += (corr-1.)*(corr-1.);
158 +                                */
159 +                        }
160 +                /*
161 +                fprintf(stderr, "Avg., RMS error: %.1f%%  %.1f%%\n",
162 +                                        100.*dsum/(double)nn,
163 +                                        100.*sqrt(dsum2/(double)nn));
164 +                */
165 +        }
166 +        newnode->next = dsf_list;
167 +        return(dsf_list = newnode);
168 + }
169 +
170   /* Load a set of measurements corresponding to a particular incident angle */
171   static int
172   load_bsdf_meas(const char *fname)
# Line 176 | Line 182 | load_bsdf_meas(const char *fname)
182                  fputs(": cannot open\n", stderr);
183                  return(0);
184          }
185 <        memset(bsdf_grid, 0, sizeof(bsdf_grid));
185 >        memset(dsf_grid, 0, sizeof(dsf_grid));
186                                  /* read header information */
187          while ((c = getc(fp)) == '#' || c == EOF) {
188                  if (fgets(buf, sizeof(buf), fp) == NULL) {
# Line 217 | Line 223 | load_bsdf_meas(const char *fname)
223                  ovec[1] = sin(M_PI/180.*phi_out) * ovec[2];
224                  ovec[2] = sqrt(1. - ovec[2]*ovec[2]);
225  
226 <                if (inp_is_DSF)
227 <                        val /= ovec[2]; /* convert from DSF to BSDF */
226 >                if (!inp_is_DSF)
227 >                        val *= ovec[2]; /* convert from BSDF to DSF */
228  
229                  pos_from_vec(pos, ovec);
230  
231 <                bsdf_grid[pos[0]][pos[1]].vsum += val;
232 <                bsdf_grid[pos[0]][pos[1]].nval++;
231 >                dsf_grid[pos[0]][pos[1]].vsum += val;
232 >                dsf_grid[pos[0]][pos[1]].nval++;
233          }
234          n = 0;
235          while ((c = getc(fp)) != EOF)
# Line 241 | Line 247 | load_bsdf_meas(const char *fname)
247   static void
248   compute_radii(void)
249   {
250 <        unsigned char   fill_grid[GRIDRES][GRIDRES];
251 <        int             r, r2, lastr2;
252 <        int             i, j, jn, ii, jj, inear, jnear;
253 <                                                /* proceed in zig-zag */
254 <        lastr2 = GRIDRES*GRIDRES;
250 >        unsigned int    fill_grid[GRIDRES][GRIDRES];
251 >        unsigned short  fill_cnt[GRIDRES][GRIDRES];
252 >        FVECT           ovec0, ovec1;
253 >        double          ang2, lastang2;
254 >        int             r, i, j, jn, ii, jj, inear, jnear;
255 >
256 >        r = GRIDRES/2;                          /* proceed in zig-zag */
257          for (i = 0; i < GRIDRES; i++)
258              for (jn = 0; jn < GRIDRES; jn++) {
259                  j = (i&1) ? jn : GRIDRES-1-jn;
260 <                if (bsdf_grid[i][j].nval)       /* find empty grid pos. */
260 >                if (dsf_grid[i][j].nval)        /* find empty grid pos. */
261                          continue;
262 <                r = (int)sqrt((double)lastr2) + 2;
262 >                vec_from_pos(ovec0, i, j);
263                  inear = jnear = -1;             /* find nearest non-empty */
264 <                lastr2 = 2*GRIDRES*GRIDRES;
264 >                lastang2 = M_PI*M_PI;
265                  for (ii = i-r; ii <= i+r; ii++) {
266                      if (ii < 0) continue;
267                      if (ii >= GRIDRES) break;
268                      for (jj = j-r; jj <= j+r; jj++) {
269                          if (jj < 0) continue;
270                          if (jj >= GRIDRES) break;
271 <                        if (!bsdf_grid[ii][jj].nval)
271 >                        if (!dsf_grid[ii][jj].nval)
272                                  continue;
273 <                        r2 = (ii-i)*(ii-i) + (jj-j)*(jj-j);
274 <                        if (r2 >= lastr2)
273 >                        vec_from_pos(ovec1, ii, jj);
274 >                        ang2 = 2. - 2.*DOT(ovec0,ovec1);
275 >                        if (ang2 >= lastang2)
276                                  continue;
277 <                        lastr2 = r2;
277 >                        lastang2 = ang2;
278                          inear = ii; jnear = jj;
279                      }
280                  }
281 <                                                /* record if > previous */
282 <                if (bsdf_grid[inear][jnear].hrad2 < lastr2)
283 <                        bsdf_grid[inear][jnear].hrad2 = lastr2;
281 >                if (inear < 0) {
282 >                        fputs("Could not find non-empty neighbor!\n", stderr);
283 >                        exit(1);
284 >                }
285 >                ang2 = sqrt(lastang2);
286 >                r = ANG2R(ang2);                /* record if > previous */
287 >                if (r > dsf_grid[inear][jnear].crad)
288 >                        dsf_grid[inear][jnear].crad = r;
289 >                                                /* next search radius */
290 >                r = ang2*(2.*GRIDRES/M_PI) + 1;
291              }
292 <                                                /* fill in others */
292 >                                                /* blur radii over hemisphere */
293          memset(fill_grid, 0, sizeof(fill_grid));
294 +        memset(fill_cnt, 0, sizeof(fill_cnt));
295          for (i = 0; i < GRIDRES; i++)
296              for (j = 0; j < GRIDRES; j++) {
297 <                if (!bsdf_grid[i][j].nval)
298 <                        continue;
299 <                if (bsdf_grid[i][j].hrad2)
283 <                        continue;
284 <                r = GRIDRES/20;
285 <                lastr2 = 2*r*r;
297 >                if (!dsf_grid[i][j].crad)
298 >                        continue;               /* missing distance */
299 >                r = R2ANG(dsf_grid[i][j].crad)*(2.*RSCA*GRIDRES/M_PI);
300                  for (ii = i-r; ii <= i+r; ii++) {
301                      if (ii < 0) continue;
302                      if (ii >= GRIDRES) break;
303                      for (jj = j-r; jj <= j+r; jj++) {
304                          if (jj < 0) continue;
305                          if (jj >= GRIDRES) break;
306 <                        if (!bsdf_grid[ii][jj].hrad2)
306 >                        if ((ii-i)*(ii-i) + (jj-j)*(jj-j) > r*r)
307                                  continue;
308 <                        r2 = (ii-i)*(ii-i) + (jj-j)*(jj-j);
309 <                        if (r2 >= lastr2)
296 <                                continue;
297 <                        fill_grid[i][j] = bsdf_grid[ii][jj].hrad2;
298 <                        lastr2 = r2;
308 >                        fill_grid[ii][jj] += dsf_grid[i][j].crad;
309 >                        fill_cnt[ii][jj]++;
310                      }
311                  }
312              }
313 +                                                /* copy back averaged radii */
314          for (i = 0; i < GRIDRES; i++)
315              for (j = 0; j < GRIDRES; j++)
316 <                if (fill_grid[i][j])
317 <                        bsdf_grid[i][j].hrad2 = fill_grid[i][j];
316 >                if (fill_cnt[i][j])
317 >                        dsf_grid[i][j].crad = fill_grid[i][j]/fill_cnt[i][j];
318   }
319  
320   /* Cull points for more uniform distribution */
321   static void
322   cull_values(void)
323   {
324 <        int     i, j, ii, jj, r, r2;
324 >        FVECT   ovec0, ovec1;
325 >        double  maxang, maxang2;
326 >        int     i, j, ii, jj, r;
327                                                  /* simple greedy algorithm */
328          for (i = 0; i < GRIDRES; i++)
329              for (j = 0; j < GRIDRES; j++) {
330 <                if (!bsdf_grid[i][j].nval)
330 >                if (!dsf_grid[i][j].nval)
331                          continue;
332 <                if (!bsdf_grid[i][j].hrad2)
333 <                        continue;
334 <                r = (int)(2.*sqrt((double)bsdf_grid[i][j].hrad2) + .9999);
332 >                if (!dsf_grid[i][j].crad)
333 >                        continue;               /* shouldn't happen */
334 >                vec_from_pos(ovec0, i, j);
335 >                maxang = 2.*R2ANG(dsf_grid[i][j].crad);
336 >                if (maxang > ovec0[2])          /* clamp near horizon */
337 >                        maxang = ovec0[2];
338 >                r = maxang*(2.*GRIDRES/M_PI) + 1;
339 >                maxang2 = maxang*maxang;
340                  for (ii = i-r; ii <= i+r; ii++) {
341                      if (ii < 0) continue;
342                      if (ii >= GRIDRES) break;
343                      for (jj = j-r; jj <= j+r; jj++) {
344                          if (jj < 0) continue;
345                          if (jj >= GRIDRES) break;
346 <                        if (!bsdf_grid[ii][jj].nval)
346 >                        if (!dsf_grid[ii][jj].nval)
347                                  continue;
348 <                        r2 = (ii-i)*(ii-i) + (jj-j)*(jj-j);
349 <                        if (!r2 | (r2 > r*r))
348 >                        if ((ii == i) & (jj == j))
349 >                                continue;       /* don't get self-absorbed */
350 >                        vec_from_pos(ovec1, ii, jj);
351 >                        if (2. - 2.*DOT(ovec0,ovec1) >= maxang2)
352                                  continue;
353 <                                                /* absorb victim's value */
354 <                        bsdf_grid[i][j].vsum += bsdf_grid[ii][jj].vsum;
355 <                        bsdf_grid[i][j].nval += bsdf_grid[ii][jj].nval;
356 <                        memset(&bsdf_grid[ii][jj], 0, sizeof(GRIDVAL));
353 >                                                /* absorb sum */
354 >                        dsf_grid[i][j].vsum += dsf_grid[ii][jj].vsum;
355 >                        dsf_grid[i][j].nval += dsf_grid[ii][jj].nval;
356 >                                                /* keep value, though */
357 >                        dsf_grid[ii][jj].vsum /= (double)dsf_grid[ii][jj].nval;
358 >                        dsf_grid[ii][jj].nval = 0;
359                      }
360                  }
361              }
# Line 356 | Line 379 | main(int argc, char *argv[])
379          }
380          if (!load_bsdf_meas(argv[1]))
381                  return(1);
359 –                                                /* produce spheres at meas. */
360 –        puts("void plastic orange\n0\n0\n5 .6 .4 .01 .04 .08\n");
361 –        n = 0;
362 –        for (i = 0; i < GRIDRES; i++)
363 –            for (j = 0; j < GRIDRES; j++)
364 –                if (bsdf_grid[i][j].nval) {
365 –                        double  bsdf = bsdf_grid[i][j].vsum /
366 –                                        (double)bsdf_grid[i][j].nval;
367 –                        FVECT   dir;
382  
369 –                        vec_from_pos(dir, i, j);
370 –                        printf("orange sphere s%04d\n0\n0\n", ++n);
371 –                        printf("4 %.6g %.6g %.6g .0015\n\n",
372 –                                        dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
373 –                }
383          compute_radii();
384          cull_values();
385 <                                                /* highlight chosen values */
385 >        make_rbfrep();
386 >                                                /* produce spheres at meas. */
387 >        puts("void plastic yellow\n0\n0\n5 .6 .4 .01 .04 .08\n");
388          puts("void plastic pink\n0\n0\n5 .5 .05 .9 .04 .08\n");
389          n = 0;
390          for (i = 0; i < GRIDRES; i++)
391              for (j = 0; j < GRIDRES; j++)
392 <                if (bsdf_grid[i][j].nval) {
382 <                        bsdf = bsdf_grid[i][j].vsum /
383 <                                        (double)bsdf_grid[i][j].nval;
392 >                if (dsf_grid[i][j].vsum > .0f) {
393                          vec_from_pos(dir, i, j);
394 <                        printf("pink cone c%04d\n0\n0\n8\n", ++n);
395 <                        printf("\t%.6g %.6g %.6g\n",
394 >                        bsdf = dsf_grid[i][j].vsum / dir[2];
395 >                        if (dsf_grid[i][j].nval) {
396 >                                printf("pink cone c%04d\n0\n0\n8\n", ++n);
397 >                                printf("\t%.6g %.6g %.6g\n",
398                                          dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
399 <                        printf("\t%.6g %.6g %.6g\n",
399 >                                printf("\t%.6g %.6g %.6g\n",
400                                          dir[0]*(bsdf+.005), dir[1]*(bsdf+.005),
401                                          dir[2]*(bsdf+.005));
402 <                        puts("\t.003\t0\n");
402 >                                puts("\t.003\t0\n");
403 >                        } else {
404 >                                vec_from_pos(dir, i, j);
405 >                                printf("yellow sphere s%04d\n0\n0\n", ++n);
406 >                                printf("4 %.6g %.6g %.6g .0015\n\n",
407 >                                        dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
408 >                        }
409                  }
410                                                  /* output continuous surface */
394 –        make_rbfrep();
411          puts("void trans tgreen\n0\n0\n7 .7 1 .7 .04 .04 .9 .9\n");
412          fflush(stdout);
413 <        sprintf(buf, "gensurf tgreen bsdf - - - %d %d", GRIDRES, GRIDRES);
413 >        sprintf(buf, "gensurf tgreen bsdf - - - %d %d", GRIDRES-1, GRIDRES-1);
414          pfp = popen(buf, "w");
415          if (pfp == NULL) {
416                  fputs(buf, stderr);
# Line 404 | Line 420 | main(int argc, char *argv[])
420          for (i = 0; i < GRIDRES; i++)
421              for (j = 0; j < GRIDRES; j++) {
422                  vec_from_pos(dir, i, j);
423 <                bsdf = eval_rbfrep(bsdf_list, dir);
423 >                bsdf = eval_rbfrep(dsf_list, dir) / dir[2];
424                  fprintf(pfp, "%.8e %.8e %.8e\n",
425                                  dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
426              }

Diff Legend

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