ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/pabopto2xml.c
Revision: 2.6
Committed: Sun Aug 26 19:40:02 2012 UTC (11 years, 9 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.5: +26 -19 lines
Log Message:
Changed to adaptive iteration based on convergence

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.6 static const char RCSid[] = "$Id: pabopto2xml.c,v 2.5 2012/08/25 22:39:03 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 greg 2.6 /* convert to/from coded radians */
26     #define ANG2R(r) (int)((r)*((1<<16)/M_PI))
27 greg 2.2 #define R2ANG(c) (((c)+.5)*(M_PI/(1<<16)))
28 greg 2.1
29     typedef struct {
30 greg 2.5 float vsum; /* DSF sum */
31 greg 2.1 unsigned short nval; /* number of values in sum */
32 greg 2.2 unsigned short crad; /* radius (coded angle) */
33 greg 2.1 } GRIDVAL; /* grid value */
34    
35     typedef struct {
36 greg 2.5 float peak; /* lobe value at peak */
37 greg 2.2 unsigned short crad; /* radius (coded angle) */
38 greg 2.1 unsigned char gx, gy; /* grid position */
39     } RBFVAL; /* radial basis function value */
40    
41     typedef struct s_rbflist {
42     struct s_rbflist *next; /* next in our RBF list */
43     FVECT invec; /* incident vector direction */
44     int nrbf; /* number of RBFs */
45     RBFVAL rbfa[1]; /* RBF array (extends struct) */
46 greg 2.5 } RBFLIST; /* RBF representation of DSF @ 1 incidence */
47 greg 2.1
48     /* our loaded grid for this incident angle */
49     static double theta_in_deg, phi_in_deg;
50 greg 2.5 static GRIDVAL dsf_grid[GRIDRES][GRIDRES];
51 greg 2.1
52 greg 2.5 /* processed incident DSF measurements */
53     static RBFLIST *dsf_list = NULL;
54 greg 2.1
55 greg 2.3 /* Compute outgoing vector from grid position */
56     static void
57     vec_from_pos(FVECT vec, int xpos, int ypos)
58 greg 2.1 {
59 greg 2.3 double uv[2];
60     double r2;
61    
62     SDsquare2disk(uv, (1./GRIDRES)*(xpos+.5), (1./GRIDRES)*(ypos+.5));
63     /* uniform hemispherical projection */
64     r2 = uv[0]*uv[0] + uv[1]*uv[1];
65     vec[0] = vec[1] = sqrt(2. - r2);
66     vec[0] *= uv[0];
67     vec[1] *= uv[1];
68     vec[2] = 1. - r2;
69 greg 2.1 }
70    
71     /* Compute grid position from normalized outgoing vector */
72     static void
73     pos_from_vec(int pos[2], const FVECT vec)
74     {
75     double sq[2]; /* uniform hemispherical projection */
76     double norm = 1./sqrt(1. + vec[2]);
77    
78     SDdisk2square(sq, vec[0]*norm, vec[1]*norm);
79    
80     pos[0] = (int)(sq[0]*GRIDRES);
81     pos[1] = (int)(sq[1]*GRIDRES);
82     }
83    
84 greg 2.5 /* Evaluate RBF for DSF at the given normalized outgoing direction */
85 greg 2.1 static double
86     eval_rbfrep(const RBFLIST *rp, const FVECT outvec)
87     {
88     double res = .0;
89     const RBFVAL *rbfp;
90     FVECT odir;
91     double sig2;
92     int n;
93    
94     rbfp = rp->rbfa;
95     for (n = rp->nrbf; n--; rbfp++) {
96     vec_from_pos(odir, rbfp->gx, rbfp->gy);
97 greg 2.2 sig2 = R2ANG(rbfp->crad);
98     sig2 = (DOT(odir,outvec) - 1.) / (sig2*sig2);
99 greg 2.1 if (sig2 > -19.)
100 greg 2.5 res += rbfp->peak * exp(sig2);
101 greg 2.1 }
102     return(res);
103     }
104    
105 greg 2.3 /* Count up filled nodes and build RBF representation from current grid */
106     static RBFLIST *
107     make_rbfrep(void)
108     {
109 greg 2.6 int niter = 16;
110     double lastVar, thisVar = 100.;
111 greg 2.3 int nn;
112     RBFLIST *newnode;
113     int i, j;
114    
115     nn = 0; /* count selected bins */
116     for (i = 0; i < GRIDRES; i++)
117     for (j = 0; j < GRIDRES; j++)
118 greg 2.6 nn += dsf_grid[i][j].nval;
119 greg 2.3 /* allocate RBF array */
120     newnode = (RBFLIST *)malloc(sizeof(RBFLIST) + sizeof(RBFVAL)*(nn-1));
121     if (newnode == NULL) {
122     fputs("Out of memory in make_rbfrep\n", stderr);
123     exit(1);
124     }
125     newnode->next = NULL;
126     newnode->invec[2] = sin(M_PI/180.*theta_in_deg);
127     newnode->invec[0] = cos(M_PI/180.*phi_in_deg)*newnode->invec[2];
128     newnode->invec[1] = sin(M_PI/180.*phi_in_deg)*newnode->invec[2];
129     newnode->invec[2] = sqrt(1. - newnode->invec[2]*newnode->invec[2]);
130     newnode->nrbf = nn;
131     nn = 0; /* fill RBF array */
132     for (i = 0; i < GRIDRES; i++)
133     for (j = 0; j < GRIDRES; j++)
134 greg 2.5 if (dsf_grid[i][j].nval) {
135 greg 2.6 newnode->rbfa[nn].peak = dsf_grid[i][j].vsum;
136 greg 2.5 newnode->rbfa[nn].crad = RSCA*dsf_grid[i][j].crad + .5;
137 greg 2.3 newnode->rbfa[nn].gx = i;
138     newnode->rbfa[nn].gy = j;
139     ++nn;
140     }
141 greg 2.6 /* iterate to improve interpolation accuracy */
142     do {
143 greg 2.4 double dsum = .0, dsum2 = .0;
144 greg 2.3 nn = 0;
145     for (i = 0; i < GRIDRES; i++)
146     for (j = 0; j < GRIDRES; j++)
147 greg 2.5 if (dsf_grid[i][j].nval) {
148 greg 2.3 FVECT odir;
149 greg 2.6 double corr;
150 greg 2.3 vec_from_pos(odir, i, j);
151 greg 2.6 newnode->rbfa[nn++].peak *= corr =
152 greg 2.5 dsf_grid[i][j].vsum /
153 greg 2.3 eval_rbfrep(newnode, odir);
154 greg 2.4 dsum += corr - 1.;
155     dsum2 += (corr-1.)*(corr-1.);
156 greg 2.3 }
157 greg 2.6 lastVar = thisVar;
158     thisVar = dsum2/(double)nn;
159 greg 2.4 /*
160     fprintf(stderr, "Avg., RMS error: %.1f%% %.1f%%\n",
161     100.*dsum/(double)nn,
162 greg 2.6 100.*sqrt(thisVar));
163 greg 2.4 */
164 greg 2.6 } while (--niter > 0 && lastVar-thisVar > 0.02*lastVar);
165    
166 greg 2.5 newnode->next = dsf_list;
167     return(dsf_list = newnode);
168 greg 2.3 }
169    
170 greg 2.1 /* Load a set of measurements corresponding to a particular incident angle */
171     static int
172     load_bsdf_meas(const char *fname)
173     {
174     FILE *fp = fopen(fname, "r");
175     int inp_is_DSF = -1;
176     double theta_out, phi_out, val;
177     char buf[2048];
178     int n, c;
179    
180     if (fp == NULL) {
181     fputs(fname, stderr);
182     fputs(": cannot open\n", stderr);
183     return(0);
184     }
185 greg 2.5 memset(dsf_grid, 0, sizeof(dsf_grid));
186 greg 2.1 /* read header information */
187     while ((c = getc(fp)) == '#' || c == EOF) {
188     if (fgets(buf, sizeof(buf), fp) == NULL) {
189     fputs(fname, stderr);
190     fputs(": unexpected EOF\n", stderr);
191     fclose(fp);
192     return(0);
193     }
194     if (!strcmp(buf, "format: theta phi DSF\n")) {
195     inp_is_DSF = 1;
196     continue;
197     }
198     if (!strcmp(buf, "format: theta phi BSDF\n")) {
199     inp_is_DSF = 0;
200     continue;
201     }
202     if (sscanf(buf, "intheta %lf", &theta_in_deg) == 1)
203     continue;
204     if (sscanf(buf, "inphi %lf", &phi_in_deg) == 1)
205     continue;
206     if (sscanf(buf, "incident_angle %lf %lf",
207     &theta_in_deg, &phi_in_deg) == 2)
208     continue;
209     }
210     if (inp_is_DSF < 0) {
211     fputs(fname, stderr);
212     fputs(": unknown format\n", stderr);
213     fclose(fp);
214     return(0);
215     }
216     ungetc(c, fp); /* read actual data */
217     while (fscanf(fp, "%lf %lf %lf\n", &theta_out, &phi_out, &val) == 3) {
218     FVECT ovec;
219     int pos[2];
220    
221     ovec[2] = sin(M_PI/180.*theta_out);
222     ovec[0] = cos(M_PI/180.*phi_out) * ovec[2];
223     ovec[1] = sin(M_PI/180.*phi_out) * ovec[2];
224     ovec[2] = sqrt(1. - ovec[2]*ovec[2]);
225    
226 greg 2.5 if (!inp_is_DSF)
227     val *= ovec[2]; /* convert from BSDF to DSF */
228 greg 2.1
229     pos_from_vec(pos, ovec);
230    
231 greg 2.5 dsf_grid[pos[0]][pos[1]].vsum += val;
232     dsf_grid[pos[0]][pos[1]].nval++;
233 greg 2.1 }
234     n = 0;
235     while ((c = getc(fp)) != EOF)
236     n += !isspace(c);
237     if (n)
238     fprintf(stderr,
239     "%s: warning: %d unexpected characters past EOD\n",
240     fname, n);
241     fclose(fp);
242     return(1);
243     }
244    
245     /* Compute radii for non-empty bins */
246     /* (distance to furthest empty bin for which non-empty bin is the closest) */
247     static void
248     compute_radii(void)
249     {
250 greg 2.4 unsigned int fill_grid[GRIDRES][GRIDRES];
251     unsigned short fill_cnt[GRIDRES][GRIDRES];
252 greg 2.2 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 greg 2.1 for (i = 0; i < GRIDRES; i++)
258     for (jn = 0; jn < GRIDRES; jn++) {
259     j = (i&1) ? jn : GRIDRES-1-jn;
260 greg 2.5 if (dsf_grid[i][j].nval) /* find empty grid pos. */
261 greg 2.1 continue;
262 greg 2.2 vec_from_pos(ovec0, i, j);
263 greg 2.1 inear = jnear = -1; /* find nearest non-empty */
264 greg 2.2 lastang2 = M_PI*M_PI;
265 greg 2.1 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 greg 2.5 if (!dsf_grid[ii][jj].nval)
272 greg 2.1 continue;
273 greg 2.2 vec_from_pos(ovec1, ii, jj);
274     ang2 = 2. - 2.*DOT(ovec0,ovec1);
275     if (ang2 >= lastang2)
276 greg 2.1 continue;
277 greg 2.2 lastang2 = ang2;
278 greg 2.1 inear = ii; jnear = jj;
279     }
280     }
281 greg 2.2 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 greg 2.5 if (r > dsf_grid[inear][jnear].crad)
288     dsf_grid[inear][jnear].crad = r;
289 greg 2.2 /* next search radius */
290     r = ang2*(2.*GRIDRES/M_PI) + 1;
291 greg 2.1 }
292 greg 2.4 /* blur radii over hemisphere */
293 greg 2.1 memset(fill_grid, 0, sizeof(fill_grid));
294 greg 2.4 memset(fill_cnt, 0, sizeof(fill_cnt));
295 greg 2.1 for (i = 0; i < GRIDRES; i++)
296     for (j = 0; j < GRIDRES; j++) {
297 greg 2.5 if (!dsf_grid[i][j].crad)
298 greg 2.4 continue; /* missing distance */
299 greg 2.5 r = R2ANG(dsf_grid[i][j].crad)*(2.*RSCA*GRIDRES/M_PI);
300 greg 2.1 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 greg 2.4 if ((ii-i)*(ii-i) + (jj-j)*(jj-j) > r*r)
307 greg 2.1 continue;
308 greg 2.5 fill_grid[ii][jj] += dsf_grid[i][j].crad;
309 greg 2.4 fill_cnt[ii][jj]++;
310 greg 2.1 }
311     }
312     }
313 greg 2.6 /* copy back blurred radii */
314 greg 2.1 for (i = 0; i < GRIDRES; i++)
315     for (j = 0; j < GRIDRES; j++)
316 greg 2.4 if (fill_cnt[i][j])
317 greg 2.5 dsf_grid[i][j].crad = fill_grid[i][j]/fill_cnt[i][j];
318 greg 2.1 }
319    
320 greg 2.6 /* Cull points for more uniform distribution, leave all nval 0 or 1 */
321 greg 2.1 static void
322     cull_values(void)
323     {
324 greg 2.2 FVECT ovec0, ovec1;
325     double maxang, maxang2;
326     int i, j, ii, jj, r;
327 greg 2.1 /* simple greedy algorithm */
328     for (i = 0; i < GRIDRES; i++)
329     for (j = 0; j < GRIDRES; j++) {
330 greg 2.5 if (!dsf_grid[i][j].nval)
331 greg 2.1 continue;
332 greg 2.5 if (!dsf_grid[i][j].crad)
333 greg 2.2 continue; /* shouldn't happen */
334     vec_from_pos(ovec0, i, j);
335 greg 2.5 maxang = 2.*R2ANG(dsf_grid[i][j].crad);
336 greg 2.2 if (maxang > ovec0[2]) /* clamp near horizon */
337     maxang = ovec0[2];
338     r = maxang*(2.*GRIDRES/M_PI) + 1;
339     maxang2 = maxang*maxang;
340 greg 2.1 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 greg 2.5 if (!dsf_grid[ii][jj].nval)
347 greg 2.1 continue;
348 greg 2.2 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 greg 2.1 continue;
353 greg 2.2 /* absorb sum */
354 greg 2.5 dsf_grid[i][j].vsum += dsf_grid[ii][jj].vsum;
355     dsf_grid[i][j].nval += dsf_grid[ii][jj].nval;
356 greg 2.2 /* keep value, though */
357 greg 2.6 dsf_grid[ii][jj].vsum /= (float)dsf_grid[ii][jj].nval;
358 greg 2.5 dsf_grid[ii][jj].nval = 0;
359 greg 2.1 }
360     }
361     }
362 greg 2.6 /* final averaging pass */
363     for (i = 0; i < GRIDRES; i++)
364     for (j = 0; j < GRIDRES; j++)
365     if (dsf_grid[i][j].nval > 1) {
366     dsf_grid[i][j].vsum /= (float)dsf_grid[i][j].nval;
367     dsf_grid[i][j].nval = 1;
368     }
369 greg 2.1 }
370    
371    
372     #if 1
373     /* Test main produces a Radiance model from the given input file */
374     int
375     main(int argc, char *argv[])
376     {
377     char buf[128];
378     FILE *pfp;
379     double bsdf;
380     FVECT dir;
381     int i, j, n;
382    
383     if (argc != 2) {
384     fprintf(stderr, "Usage: %s input.dat > output.rad\n", argv[0]);
385     return(1);
386     }
387     if (!load_bsdf_meas(argv[1]))
388     return(1);
389    
390     compute_radii();
391     cull_values();
392 greg 2.3 make_rbfrep();
393     /* produce spheres at meas. */
394     puts("void plastic yellow\n0\n0\n5 .6 .4 .01 .04 .08\n");
395 greg 2.1 puts("void plastic pink\n0\n0\n5 .5 .05 .9 .04 .08\n");
396     n = 0;
397     for (i = 0; i < GRIDRES; i++)
398     for (j = 0; j < GRIDRES; j++)
399 greg 2.5 if (dsf_grid[i][j].vsum > .0f) {
400 greg 2.1 vec_from_pos(dir, i, j);
401 greg 2.5 bsdf = dsf_grid[i][j].vsum / dir[2];
402     if (dsf_grid[i][j].nval) {
403 greg 2.3 printf("pink cone c%04d\n0\n0\n8\n", ++n);
404     printf("\t%.6g %.6g %.6g\n",
405 greg 2.1 dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
406 greg 2.3 printf("\t%.6g %.6g %.6g\n",
407 greg 2.1 dir[0]*(bsdf+.005), dir[1]*(bsdf+.005),
408     dir[2]*(bsdf+.005));
409 greg 2.3 puts("\t.003\t0\n");
410     } else {
411     vec_from_pos(dir, i, j);
412     printf("yellow sphere s%04d\n0\n0\n", ++n);
413     printf("4 %.6g %.6g %.6g .0015\n\n",
414     dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
415     }
416 greg 2.1 }
417     /* output continuous surface */
418     puts("void trans tgreen\n0\n0\n7 .7 1 .7 .04 .04 .9 .9\n");
419     fflush(stdout);
420 greg 2.5 sprintf(buf, "gensurf tgreen bsdf - - - %d %d", GRIDRES-1, GRIDRES-1);
421 greg 2.1 pfp = popen(buf, "w");
422     if (pfp == NULL) {
423     fputs(buf, stderr);
424     fputs(": cannot start command\n", stderr);
425     return(1);
426     }
427     for (i = 0; i < GRIDRES; i++)
428     for (j = 0; j < GRIDRES; j++) {
429     vec_from_pos(dir, i, j);
430 greg 2.5 bsdf = eval_rbfrep(dsf_list, dir) / dir[2];
431 greg 2.1 fprintf(pfp, "%.8e %.8e %.8e\n",
432     dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
433     }
434     return(pclose(pfp)==0 ? 0 : 1);
435     }
436     #endif