ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/pabopto2xml.c
Revision: 2.3
Committed: Fri Aug 24 22:08:50 2012 UTC (12 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.2: +86 -76 lines
Log Message:
Added optimization for better fitting results

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.3 static const char RCSid[] = "$Id: pabopto2xml.c,v 2.2 2012/08/24 20:55:28 greg Exp $";
3 greg 2.1 #endif
4     /*
5     * Convert PAB-Opto measurements to XML format using tensor tree representation
6     * Employs Bonneel et al. Earth Mover's Distance interpolant.
7     *
8     * G.Ward
9     */
10    
11     #define _USE_MATH_DEFINES
12     #include <stdio.h>
13     #include <stdlib.h>
14     #include <string.h>
15     #include <ctype.h>
16     #include <math.h>
17     #include "bsdf.h"
18    
19     #ifndef GRIDRES
20     #define GRIDRES 200 /* max. grid resolution per side */
21     #endif
22    
23 greg 2.3 #define RSCA 2.7 /* radius scaling factor (empirical) */
24 greg 2.2
25     #define R2ANG(c) (((c)+.5)*(M_PI/(1<<16)))
26     #define ANG2R(r) (int)((r)*((1<<16)/M_PI))
27 greg 2.1
28     typedef struct {
29     float vsum; /* BSDF sum */
30     unsigned short nval; /* number of values in sum */
31 greg 2.2 unsigned short crad; /* radius (coded angle) */
32 greg 2.1 } GRIDVAL; /* grid value */
33    
34     typedef struct {
35 greg 2.2 float bsdf; /* lobe value at peak */
36     unsigned short crad; /* radius (coded angle) */
37 greg 2.1 unsigned char gx, gy; /* grid position */
38     } RBFVAL; /* radial basis function value */
39    
40     typedef struct s_rbflist {
41     struct s_rbflist *next; /* next in our RBF list */
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 */
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];
50    
51     /* processed incident BSDF measurements */
52     static RBFLIST *bsdf_list = NULL;
53    
54 greg 2.3 /* Compute outgoing vector from grid position */
55     static void
56     vec_from_pos(FVECT vec, int xpos, int ypos)
57 greg 2.1 {
58 greg 2.3 double uv[2];
59     double r2;
60    
61     SDsquare2disk(uv, (1./GRIDRES)*(xpos+.5), (1./GRIDRES)*(ypos+.5));
62     /* uniform hemispherical projection */
63     r2 = uv[0]*uv[0] + uv[1]*uv[1];
64     vec[0] = vec[1] = sqrt(2. - r2);
65     vec[0] *= uv[0];
66     vec[1] *= uv[1];
67     vec[2] = 1. - r2;
68 greg 2.1 }
69    
70     /* Compute grid position from normalized outgoing vector */
71     static void
72     pos_from_vec(int pos[2], const FVECT vec)
73     {
74     double sq[2]; /* uniform hemispherical projection */
75     double norm = 1./sqrt(1. + vec[2]);
76    
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 */
84     static double
85     eval_rbfrep(const RBFLIST *rp, const FVECT outvec)
86     {
87     double res = .0;
88     const RBFVAL *rbfp;
89     FVECT odir;
90     double sig2;
91     int n;
92    
93     rbfp = rp->rbfa;
94     for (n = rp->nrbf; n--; rbfp++) {
95     vec_from_pos(odir, rbfp->gx, rbfp->gy);
96 greg 2.2 sig2 = R2ANG(rbfp->crad);
97     sig2 = (DOT(odir,outvec) - 1.) / (sig2*sig2);
98 greg 2.1 if (sig2 > -19.)
99     res += rbfp->bsdf * exp(sig2);
100     }
101     return(res);
102     }
103    
104 greg 2.3 /* Count up filled nodes and build RBF representation from current grid */
105     static RBFLIST *
106     make_rbfrep(void)
107     {
108     int niter = 4;
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 += (bsdf_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 (bsdf_grid[i][j].nval) {
133     newnode->rbfa[nn].bsdf =
134     bsdf_grid[i][j].vsum /=
135     (double)bsdf_grid[i][j].nval;
136     bsdf_grid[i][j].nval = 1;
137     newnode->rbfa[nn].crad = RSCA*bsdf_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     nn = 0;
145     for (i = 0; i < GRIDRES; i++)
146     for (j = 0; j < GRIDRES; j++)
147     if (bsdf_grid[i][j].nval) {
148     FVECT odir;
149     vec_from_pos(odir, i, j);
150     newnode->rbfa[nn++].bsdf *=
151     bsdf_grid[i][j].vsum /
152     eval_rbfrep(newnode, odir);
153     }
154     }
155     newnode->next = bsdf_list;
156     return(bsdf_list = newnode);
157     }
158    
159 greg 2.1 /* Load a set of measurements corresponding to a particular incident angle */
160     static int
161     load_bsdf_meas(const char *fname)
162     {
163     FILE *fp = fopen(fname, "r");
164     int inp_is_DSF = -1;
165     double theta_out, phi_out, val;
166     char buf[2048];
167     int n, c;
168    
169     if (fp == NULL) {
170     fputs(fname, stderr);
171     fputs(": cannot open\n", stderr);
172     return(0);
173     }
174     memset(bsdf_grid, 0, sizeof(bsdf_grid));
175     /* read header information */
176     while ((c = getc(fp)) == '#' || c == EOF) {
177     if (fgets(buf, sizeof(buf), fp) == NULL) {
178     fputs(fname, stderr);
179     fputs(": unexpected EOF\n", stderr);
180     fclose(fp);
181     return(0);
182     }
183     if (!strcmp(buf, "format: theta phi DSF\n")) {
184     inp_is_DSF = 1;
185     continue;
186     }
187     if (!strcmp(buf, "format: theta phi BSDF\n")) {
188     inp_is_DSF = 0;
189     continue;
190     }
191     if (sscanf(buf, "intheta %lf", &theta_in_deg) == 1)
192     continue;
193     if (sscanf(buf, "inphi %lf", &phi_in_deg) == 1)
194     continue;
195     if (sscanf(buf, "incident_angle %lf %lf",
196     &theta_in_deg, &phi_in_deg) == 2)
197     continue;
198     }
199     if (inp_is_DSF < 0) {
200     fputs(fname, stderr);
201     fputs(": unknown format\n", stderr);
202     fclose(fp);
203     return(0);
204     }
205     ungetc(c, fp); /* read actual data */
206     while (fscanf(fp, "%lf %lf %lf\n", &theta_out, &phi_out, &val) == 3) {
207     FVECT ovec;
208     int pos[2];
209    
210     ovec[2] = sin(M_PI/180.*theta_out);
211     ovec[0] = cos(M_PI/180.*phi_out) * ovec[2];
212     ovec[1] = sin(M_PI/180.*phi_out) * ovec[2];
213     ovec[2] = sqrt(1. - ovec[2]*ovec[2]);
214    
215     if (inp_is_DSF)
216     val /= ovec[2]; /* convert from DSF to BSDF */
217    
218     pos_from_vec(pos, ovec);
219    
220     bsdf_grid[pos[0]][pos[1]].vsum += val;
221     bsdf_grid[pos[0]][pos[1]].nval++;
222     }
223     n = 0;
224     while ((c = getc(fp)) != EOF)
225     n += !isspace(c);
226     if (n)
227     fprintf(stderr,
228     "%s: warning: %d unexpected characters past EOD\n",
229     fname, n);
230     fclose(fp);
231     return(1);
232     }
233    
234     /* Compute radii for non-empty bins */
235     /* (distance to furthest empty bin for which non-empty bin is the closest) */
236     static void
237     compute_radii(void)
238     {
239 greg 2.2 unsigned short fill_grid[GRIDRES][GRIDRES];
240     FVECT ovec0, ovec1;
241     double ang2, lastang2;
242     int r2, lastr2;
243     int r, i, j, jn, ii, jj, inear, jnear;
244    
245     r = GRIDRES/2; /* proceed in zig-zag */
246 greg 2.1 for (i = 0; i < GRIDRES; i++)
247     for (jn = 0; jn < GRIDRES; jn++) {
248     j = (i&1) ? jn : GRIDRES-1-jn;
249     if (bsdf_grid[i][j].nval) /* find empty grid pos. */
250     continue;
251 greg 2.2 vec_from_pos(ovec0, i, j);
252 greg 2.1 inear = jnear = -1; /* find nearest non-empty */
253 greg 2.2 lastang2 = M_PI*M_PI;
254 greg 2.1 for (ii = i-r; ii <= i+r; ii++) {
255     if (ii < 0) continue;
256     if (ii >= GRIDRES) break;
257     for (jj = j-r; jj <= j+r; jj++) {
258     if (jj < 0) continue;
259     if (jj >= GRIDRES) break;
260     if (!bsdf_grid[ii][jj].nval)
261     continue;
262 greg 2.2 vec_from_pos(ovec1, ii, jj);
263     ang2 = 2. - 2.*DOT(ovec0,ovec1);
264     if (ang2 >= lastang2)
265 greg 2.1 continue;
266 greg 2.2 lastang2 = ang2;
267 greg 2.1 inear = ii; jnear = jj;
268     }
269     }
270 greg 2.2 if (inear < 0) {
271     fputs("Could not find non-empty neighbor!\n", stderr);
272     exit(1);
273     }
274     ang2 = sqrt(lastang2);
275     r = ANG2R(ang2); /* record if > previous */
276     if (r > bsdf_grid[inear][jnear].crad)
277     bsdf_grid[inear][jnear].crad = r;
278     /* next search radius */
279     r = ang2*(2.*GRIDRES/M_PI) + 1;
280 greg 2.1 }
281 greg 2.2 /* fill in neighbors */
282 greg 2.1 memset(fill_grid, 0, sizeof(fill_grid));
283     for (i = 0; i < GRIDRES; i++)
284     for (j = 0; j < GRIDRES; j++) {
285     if (!bsdf_grid[i][j].nval)
286 greg 2.2 continue; /* no value -- skip */
287     if (bsdf_grid[i][j].crad)
288     continue; /* has distance already */
289 greg 2.1 r = GRIDRES/20;
290 greg 2.2 lastr2 = 2*r*r + 1;
291 greg 2.1 for (ii = i-r; ii <= i+r; ii++) {
292     if (ii < 0) continue;
293     if (ii >= GRIDRES) break;
294     for (jj = j-r; jj <= j+r; jj++) {
295     if (jj < 0) continue;
296     if (jj >= GRIDRES) break;
297 greg 2.2 if (!bsdf_grid[ii][jj].crad)
298 greg 2.1 continue;
299 greg 2.2 /* OK to use approx. closest */
300 greg 2.1 r2 = (ii-i)*(ii-i) + (jj-j)*(jj-j);
301     if (r2 >= lastr2)
302     continue;
303 greg 2.2 fill_grid[i][j] = bsdf_grid[ii][jj].crad;
304 greg 2.1 lastr2 = r2;
305     }
306     }
307     }
308 greg 2.2 /* copy back filled entries */
309 greg 2.1 for (i = 0; i < GRIDRES; i++)
310     for (j = 0; j < GRIDRES; j++)
311     if (fill_grid[i][j])
312 greg 2.2 bsdf_grid[i][j].crad = fill_grid[i][j];
313 greg 2.1 }
314    
315     /* Cull points for more uniform distribution */
316     static void
317     cull_values(void)
318     {
319 greg 2.2 FVECT ovec0, ovec1;
320     double maxang, maxang2;
321     int i, j, ii, jj, r;
322 greg 2.1 /* simple greedy algorithm */
323     for (i = 0; i < GRIDRES; i++)
324     for (j = 0; j < GRIDRES; j++) {
325     if (!bsdf_grid[i][j].nval)
326     continue;
327 greg 2.2 if (!bsdf_grid[i][j].crad)
328     continue; /* shouldn't happen */
329     vec_from_pos(ovec0, i, j);
330     maxang = 2.*R2ANG(bsdf_grid[i][j].crad);
331     if (maxang > ovec0[2]) /* clamp near horizon */
332     maxang = ovec0[2];
333     r = maxang*(2.*GRIDRES/M_PI) + 1;
334     maxang2 = maxang*maxang;
335 greg 2.1 for (ii = i-r; ii <= i+r; ii++) {
336     if (ii < 0) continue;
337     if (ii >= GRIDRES) break;
338     for (jj = j-r; jj <= j+r; jj++) {
339     if (jj < 0) continue;
340     if (jj >= GRIDRES) break;
341     if (!bsdf_grid[ii][jj].nval)
342     continue;
343 greg 2.2 if ((ii == i) & (jj == j))
344     continue; /* don't get self-absorbed */
345     vec_from_pos(ovec1, ii, jj);
346     if (2. - 2.*DOT(ovec0,ovec1) >= maxang2)
347 greg 2.1 continue;
348 greg 2.2 /* absorb sum */
349 greg 2.1 bsdf_grid[i][j].vsum += bsdf_grid[ii][jj].vsum;
350     bsdf_grid[i][j].nval += bsdf_grid[ii][jj].nval;
351 greg 2.2 /* keep value, though */
352     bsdf_grid[ii][jj].vsum /= (double)bsdf_grid[ii][jj].nval;
353     bsdf_grid[ii][jj].nval = 0;
354 greg 2.1 }
355     }
356     }
357     }
358    
359    
360     #if 1
361     /* Test main produces a Radiance model from the given input file */
362     int
363     main(int argc, char *argv[])
364     {
365     char buf[128];
366     FILE *pfp;
367     double bsdf;
368     FVECT dir;
369     int i, j, n;
370    
371     if (argc != 2) {
372     fprintf(stderr, "Usage: %s input.dat > output.rad\n", argv[0]);
373     return(1);
374     }
375     if (!load_bsdf_meas(argv[1]))
376     return(1);
377    
378     compute_radii();
379     cull_values();
380 greg 2.3 make_rbfrep();
381     /* produce spheres at meas. */
382     puts("void plastic yellow\n0\n0\n5 .6 .4 .01 .04 .08\n");
383 greg 2.1 puts("void plastic pink\n0\n0\n5 .5 .05 .9 .04 .08\n");
384     n = 0;
385     for (i = 0; i < GRIDRES; i++)
386     for (j = 0; j < GRIDRES; j++)
387 greg 2.3 if (bsdf_grid[i][j].vsum > .0f) {
388     bsdf = bsdf_grid[i][j].vsum;
389 greg 2.1 vec_from_pos(dir, i, j);
390 greg 2.3 if (bsdf_grid[i][j].nval) {
391     printf("pink cone c%04d\n0\n0\n8\n", ++n);
392     printf("\t%.6g %.6g %.6g\n",
393 greg 2.1 dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
394 greg 2.3 printf("\t%.6g %.6g %.6g\n",
395 greg 2.1 dir[0]*(bsdf+.005), dir[1]*(bsdf+.005),
396     dir[2]*(bsdf+.005));
397 greg 2.3 puts("\t.003\t0\n");
398     } else {
399     vec_from_pos(dir, i, j);
400     printf("yellow sphere s%04d\n0\n0\n", ++n);
401     printf("4 %.6g %.6g %.6g .0015\n\n",
402     dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
403     }
404 greg 2.1 }
405     /* output continuous surface */
406     puts("void trans tgreen\n0\n0\n7 .7 1 .7 .04 .04 .9 .9\n");
407     fflush(stdout);
408     sprintf(buf, "gensurf tgreen bsdf - - - %d %d", GRIDRES, GRIDRES);
409     pfp = popen(buf, "w");
410     if (pfp == NULL) {
411     fputs(buf, stderr);
412     fputs(": cannot start command\n", stderr);
413     return(1);
414     }
415     for (i = 0; i < GRIDRES; i++)
416     for (j = 0; j < GRIDRES; j++) {
417     vec_from_pos(dir, i, j);
418     bsdf = eval_rbfrep(bsdf_list, dir);
419     fprintf(pfp, "%.8e %.8e %.8e\n",
420     dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
421     }
422     return(pclose(pfp)==0 ? 0 : 1);
423     }
424     #endif