ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/pabopto2xml.c
Revision: 2.1
Committed: Fri Aug 24 15:20:18 2012 UTC (11 years, 9 months ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
Introducing program for interpolating BSDF data

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2     static const char RCSid[] = "$Id$";
3     #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     #define RSCA 1.9 /* radius scaling factor (empirical) */
24     #define MSCA .12 /* magnitude scaling (empirical) */
25    
26     typedef struct {
27     float vsum; /* BSDF sum */
28     unsigned short nval; /* number of values in sum */
29     unsigned short hrad2; /* half radius squared */
30     } GRIDVAL; /* grid value */
31    
32     typedef struct {
33     float bsdf; /* BSDF value at peak */
34     unsigned short rad; /* radius */
35     unsigned char gx, gy; /* grid position */
36     } RBFVAL; /* radial basis function value */
37    
38     typedef struct s_rbflist {
39     struct s_rbflist *next; /* next in our RBF list */
40     FVECT invec; /* incident vector direction */
41     int nrbf; /* number of RBFs */
42     RBFVAL rbfa[1]; /* RBF array (extends struct) */
43     } RBFLIST; /* RBF representation of BSDF @ 1 incidence */
44    
45     /* our loaded grid for this incident angle */
46     static double theta_in_deg, phi_in_deg;
47     static GRIDVAL bsdf_grid[GRIDRES][GRIDRES];
48    
49     /* processed incident BSDF measurements */
50     static RBFLIST *bsdf_list = NULL;
51    
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    
103     /* Compute outgoing vector from grid position */
104     static void
105     vec_from_pos(FVECT vec, int xpos, int ypos)
106     {
107     double uv[2];
108     double r2;
109    
110     SDsquare2disk(uv, (1./GRIDRES)*(xpos+.5), (1./GRIDRES)*(ypos+.5));
111     /* uniform hemispherical projection */
112     r2 = uv[0]*uv[0] + uv[1]*uv[1];
113     vec[0] = vec[1] = sqrt(2. - r2);
114     vec[0] *= uv[0];
115     vec[1] *= uv[1];
116     vec[2] = 1. - r2;
117     }
118    
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    
142     /* Evaluate RBF for BSDF at the given normalized outgoing direction */
143     static double
144     eval_rbfrep(const RBFLIST *rp, const FVECT outvec)
145     {
146     double res = .0;
147     const RBFVAL *rbfp;
148     FVECT odir;
149     double sig2;
150     int n;
151    
152     rbfp = rp->rbfa;
153     for (n = rp->nrbf; n--; rbfp++) {
154     vec_from_pos(odir, rbfp->gx, rbfp->gy);
155     sig2 = (DOT(odir, outvec) - 1.) /
156     ((M_PI*M_PI/(double)(GRIDRES*GRIDRES)) *
157     rbfp->rad*rbfp->rad);
158     if (sig2 > -19.)
159     res += rbfp->bsdf * exp(sig2);
160     }
161     return(res);
162     }
163    
164     /* Load a set of measurements corresponding to a particular incident angle */
165     static int
166     load_bsdf_meas(const char *fname)
167     {
168     FILE *fp = fopen(fname, "r");
169     int inp_is_DSF = -1;
170     double theta_out, phi_out, val;
171     char buf[2048];
172     int n, c;
173    
174     if (fp == NULL) {
175     fputs(fname, stderr);
176     fputs(": cannot open\n", stderr);
177     return(0);
178     }
179     memset(bsdf_grid, 0, sizeof(bsdf_grid));
180     /* read header information */
181     while ((c = getc(fp)) == '#' || c == EOF) {
182     if (fgets(buf, sizeof(buf), fp) == NULL) {
183     fputs(fname, stderr);
184     fputs(": unexpected EOF\n", stderr);
185     fclose(fp);
186     return(0);
187     }
188     if (!strcmp(buf, "format: theta phi DSF\n")) {
189     inp_is_DSF = 1;
190     continue;
191     }
192     if (!strcmp(buf, "format: theta phi BSDF\n")) {
193     inp_is_DSF = 0;
194     continue;
195     }
196     if (sscanf(buf, "intheta %lf", &theta_in_deg) == 1)
197     continue;
198     if (sscanf(buf, "inphi %lf", &phi_in_deg) == 1)
199     continue;
200     if (sscanf(buf, "incident_angle %lf %lf",
201     &theta_in_deg, &phi_in_deg) == 2)
202     continue;
203     }
204     if (inp_is_DSF < 0) {
205     fputs(fname, stderr);
206     fputs(": unknown format\n", stderr);
207     fclose(fp);
208     return(0);
209     }
210     ungetc(c, fp); /* read actual data */
211     while (fscanf(fp, "%lf %lf %lf\n", &theta_out, &phi_out, &val) == 3) {
212     FVECT ovec;
213     int pos[2];
214    
215     ovec[2] = sin(M_PI/180.*theta_out);
216     ovec[0] = cos(M_PI/180.*phi_out) * ovec[2];
217     ovec[1] = sin(M_PI/180.*phi_out) * ovec[2];
218     ovec[2] = sqrt(1. - ovec[2]*ovec[2]);
219    
220     if (inp_is_DSF)
221     val /= ovec[2]; /* convert from DSF to BSDF */
222    
223     pos_from_vec(pos, ovec);
224    
225     bsdf_grid[pos[0]][pos[1]].vsum += val;
226     bsdf_grid[pos[0]][pos[1]].nval++;
227     }
228     n = 0;
229     while ((c = getc(fp)) != EOF)
230     n += !isspace(c);
231     if (n)
232     fprintf(stderr,
233     "%s: warning: %d unexpected characters past EOD\n",
234     fname, n);
235     fclose(fp);
236     return(1);
237     }
238    
239     /* Compute radii for non-empty bins */
240     /* (distance to furthest empty bin for which non-empty bin is the closest) */
241     static void
242     compute_radii(void)
243     {
244     unsigned char fill_grid[GRIDRES][GRIDRES];
245     int r, r2, lastr2;
246     int i, j, jn, ii, jj, inear, jnear;
247     /* proceed in zig-zag */
248     lastr2 = GRIDRES*GRIDRES;
249     for (i = 0; i < GRIDRES; i++)
250     for (jn = 0; jn < GRIDRES; jn++) {
251     j = (i&1) ? jn : GRIDRES-1-jn;
252     if (bsdf_grid[i][j].nval) /* find empty grid pos. */
253     continue;
254     r = (int)sqrt((double)lastr2) + 2;
255     inear = jnear = -1; /* find nearest non-empty */
256     lastr2 = 2*GRIDRES*GRIDRES;
257     for (ii = i-r; ii <= i+r; ii++) {
258     if (ii < 0) continue;
259     if (ii >= GRIDRES) break;
260     for (jj = j-r; jj <= j+r; jj++) {
261     if (jj < 0) continue;
262     if (jj >= GRIDRES) break;
263     if (!bsdf_grid[ii][jj].nval)
264     continue;
265     r2 = (ii-i)*(ii-i) + (jj-j)*(jj-j);
266     if (r2 >= lastr2)
267     continue;
268     lastr2 = r2;
269     inear = ii; jnear = jj;
270     }
271     }
272     /* record if > previous */
273     if (bsdf_grid[inear][jnear].hrad2 < lastr2)
274     bsdf_grid[inear][jnear].hrad2 = lastr2;
275     }
276     /* fill in others */
277     memset(fill_grid, 0, sizeof(fill_grid));
278     for (i = 0; i < GRIDRES; i++)
279     for (j = 0; j < GRIDRES; j++) {
280     if (!bsdf_grid[i][j].nval)
281     continue;
282     if (bsdf_grid[i][j].hrad2)
283     continue;
284     r = GRIDRES/20;
285     lastr2 = 2*r*r;
286     for (ii = i-r; ii <= i+r; ii++) {
287     if (ii < 0) continue;
288     if (ii >= GRIDRES) break;
289     for (jj = j-r; jj <= j+r; jj++) {
290     if (jj < 0) continue;
291     if (jj >= GRIDRES) break;
292     if (!bsdf_grid[ii][jj].hrad2)
293     continue;
294     r2 = (ii-i)*(ii-i) + (jj-j)*(jj-j);
295     if (r2 >= lastr2)
296     continue;
297     fill_grid[i][j] = bsdf_grid[ii][jj].hrad2;
298     lastr2 = r2;
299     }
300     }
301     }
302     for (i = 0; i < GRIDRES; i++)
303     for (j = 0; j < GRIDRES; j++)
304     if (fill_grid[i][j])
305     bsdf_grid[i][j].hrad2 = fill_grid[i][j];
306     }
307    
308     /* Cull points for more uniform distribution */
309     static void
310     cull_values(void)
311     {
312     int i, j, ii, jj, r, r2;
313     /* simple greedy algorithm */
314     for (i = 0; i < GRIDRES; i++)
315     for (j = 0; j < GRIDRES; j++) {
316     if (!bsdf_grid[i][j].nval)
317     continue;
318     if (!bsdf_grid[i][j].hrad2)
319     continue;
320     r = (int)(2.*sqrt((double)bsdf_grid[i][j].hrad2) + .9999);
321     for (ii = i-r; ii <= i+r; ii++) {
322     if (ii < 0) continue;
323     if (ii >= GRIDRES) break;
324     for (jj = j-r; jj <= j+r; jj++) {
325     if (jj < 0) continue;
326     if (jj >= GRIDRES) break;
327     if (!bsdf_grid[ii][jj].nval)
328     continue;
329     r2 = (ii-i)*(ii-i) + (jj-j)*(jj-j);
330     if (!r2 | (r2 > r*r))
331     continue;
332     /* absorb victim's value */
333     bsdf_grid[i][j].vsum += bsdf_grid[ii][jj].vsum;
334     bsdf_grid[i][j].nval += bsdf_grid[ii][jj].nval;
335     memset(&bsdf_grid[ii][jj], 0, sizeof(GRIDVAL));
336     }
337     }
338     }
339     }
340    
341    
342     #if 1
343     /* Test main produces a Radiance model from the given input file */
344     int
345     main(int argc, char *argv[])
346     {
347     char buf[128];
348     FILE *pfp;
349     double bsdf;
350     FVECT dir;
351     int i, j, n;
352    
353     if (argc != 2) {
354     fprintf(stderr, "Usage: %s input.dat > output.rad\n", argv[0]);
355     return(1);
356     }
357     if (!load_bsdf_meas(argv[1]))
358     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;
368    
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     }
374     compute_radii();
375     cull_values();
376     /* highlight chosen values */
377     puts("void plastic pink\n0\n0\n5 .5 .05 .9 .04 .08\n");
378     n = 0;
379     for (i = 0; i < GRIDRES; i++)
380     for (j = 0; j < GRIDRES; j++)
381     if (bsdf_grid[i][j].nval) {
382     bsdf = bsdf_grid[i][j].vsum /
383     (double)bsdf_grid[i][j].nval;
384     vec_from_pos(dir, i, j);
385     printf("pink cone c%04d\n0\n0\n8\n", ++n);
386     printf("\t%.6g %.6g %.6g\n",
387     dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
388     printf("\t%.6g %.6g %.6g\n",
389     dir[0]*(bsdf+.005), dir[1]*(bsdf+.005),
390     dir[2]*(bsdf+.005));
391     puts("\t.003\t0\n");
392     }
393     /* output continuous surface */
394     make_rbfrep();
395     puts("void trans tgreen\n0\n0\n7 .7 1 .7 .04 .04 .9 .9\n");
396     fflush(stdout);
397     sprintf(buf, "gensurf tgreen bsdf - - - %d %d", GRIDRES, GRIDRES);
398     pfp = popen(buf, "w");
399     if (pfp == NULL) {
400     fputs(buf, stderr);
401     fputs(": cannot start command\n", stderr);
402     return(1);
403     }
404     for (i = 0; i < GRIDRES; i++)
405     for (j = 0; j < GRIDRES; j++) {
406     vec_from_pos(dir, i, j);
407     bsdf = eval_rbfrep(bsdf_list, dir);
408     fprintf(pfp, "%.8e %.8e %.8e\n",
409     dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
410     }
411     return(pclose(pfp)==0 ? 0 : 1);
412     }
413     #endif