ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/pabopto2bsdf.c
Revision: 2.37
Committed: Mon Feb 15 23:20:22 2021 UTC (3 years, 3 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.36: +49 -32 lines
Log Message:
feat: Added "up" symmetry for duplicating data @ phi += 180 degrees

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.37 static const char RCSid[] = "$Id: pabopto2bsdf.c,v 2.36 2019/07/19 17:37:56 greg Exp $";
3 greg 2.1 #endif
4     /*
5     * Load measured BSDF data in PAB-Opto format.
6 greg 2.29 * Assumes that surface-normal (Z-axis) faces into room unless -t option given.
7 greg 2.1 *
8     * G. Ward
9     */
10    
11 greg 2.20 #define _USE_MATH_DEFINES
12 greg 2.2 #include <stdlib.h>
13 greg 2.1 #include <ctype.h>
14 greg 2.6 #include <math.h>
15 greg 2.36 #include "rtio.h"
16 greg 2.1 #include "platform.h"
17     #include "bsdfrep.h"
18     /* global argv[0] */
19     char *progname;
20    
21 greg 2.2 typedef struct {
22     const char *fname; /* input file path */
23 greg 2.27 double theta, phi; /* input angles (degrees) */
24     double up_phi; /* azimuth for "up" direction */
25 greg 2.20 int igp[2]; /* input grid position */
26 greg 2.2 int isDSF; /* data is DSF (rather than BSDF)? */
27 greg 2.28 int nspec; /* number of spectral samples */
28 greg 2.2 long dstart; /* data start offset in file */
29     } PGINPUT;
30    
31     PGINPUT *inpfile; /* input files sorted by incidence */
32    
33 greg 2.29 int rev_orient = 0; /* shall we reverse surface orientation? */
34    
35 greg 2.2 /* Compare incident angles */
36     static int
37 greg 2.20 cmp_indir(const void *p1, const void *p2)
38 greg 2.2 {
39     const PGINPUT *inp1 = (const PGINPUT *)p1;
40     const PGINPUT *inp2 = (const PGINPUT *)p2;
41 greg 2.20 int ydif = inp1->igp[1] - inp2->igp[1];
42 greg 2.2
43 greg 2.20 if (ydif)
44     return(ydif);
45    
46     return(inp1->igp[0] - inp2->igp[0]);
47 greg 2.2 }
48    
49 greg 2.37 /* Assign grid position from theta and phi */
50     static void
51     set_grid_pos(PGINPUT *pip)
52     {
53     FVECT dv;
54    
55     if (pip->theta <= FTINY) {
56     pip->igp[0] = pip->igp[1] = grid_res/2 - 1;
57     return;
58     }
59     dv[2] = sin(M_PI/180.*pip->theta);
60     dv[0] = cos(M_PI/180.*pip->phi)*dv[2];
61     dv[1] = sin(M_PI/180.*pip->phi)*dv[2];
62     dv[2] = sqrt(1. - dv[2]*dv[2]);
63     pos_from_vec(pip->igp, dv);
64     }
65    
66 greg 2.2 /* Prepare a PAB-Opto input file by reading its header */
67 greg 2.1 static int
68 greg 2.2 init_pabopto_inp(const int i, const char *fname)
69 greg 2.1 {
70     FILE *fp = fopen(fname, "r");
71     char buf[2048];
72 greg 2.2 int c;
73 greg 2.1
74     if (fp == NULL) {
75     fputs(fname, stderr);
76     fputs(": cannot open\n", stderr);
77     return(0);
78     }
79 greg 2.2 inpfile[i].fname = fname;
80     inpfile[i].isDSF = -1;
81 greg 2.28 inpfile[i].nspec = 0;
82 greg 2.27 inpfile[i].up_phi = 0;
83 greg 2.2 inpfile[i].theta = inpfile[i].phi = -10001.;
84 greg 2.1 /* read header information */
85     while ((c = getc(fp)) == '#' || c == EOF) {
86 greg 2.28 char typ[64];
87 greg 2.1 if (fgets(buf, sizeof(buf), fp) == NULL) {
88     fputs(fname, stderr);
89     fputs(": unexpected EOF\n", stderr);
90     fclose(fp);
91     return(0);
92     }
93 greg 2.10 if (sscanf(buf, "sample_name \"%[^\"]\"", bsdf_name) == 1)
94     continue;
95 greg 2.28 if (sscanf(buf, "colorimetry: %s", typ) == 1) {
96     if (!strcasecmp(typ, "CIE-XYZ"))
97     inpfile[i].nspec = 3;
98     else if (!strcasecmp(typ, "CIE-Y"))
99     inpfile[i].nspec = 1;
100     continue;
101     }
102 greg 2.5 if (sscanf(buf, "format: theta phi %s", typ) == 1) {
103 greg 2.28 if (!strcasecmp(typ, "DSF"))
104 greg 2.5 inpfile[i].isDSF = 1;
105 greg 2.28 else if (!strcasecmp(typ, "BSDF") ||
106 greg 2.23 !strcasecmp(typ, "BRDF") ||
107 greg 2.28 !strcasecmp(typ, "BTDF"))
108 greg 2.5 inpfile[i].isDSF = 0;
109 greg 2.28 continue;
110 greg 2.1 }
111 greg 2.27 if (sscanf(buf, "upphi %lf", &inpfile[i].up_phi) == 1)
112     continue;
113 greg 2.2 if (sscanf(buf, "intheta %lf", &inpfile[i].theta) == 1)
114 greg 2.1 continue;
115 greg 2.2 if (sscanf(buf, "inphi %lf", &inpfile[i].phi) == 1)
116 greg 2.1 continue;
117     if (sscanf(buf, "incident_angle %lf %lf",
118 greg 2.2 &inpfile[i].theta, &inpfile[i].phi) == 2)
119 greg 2.1 continue;
120     }
121 greg 2.2 inpfile[i].dstart = ftell(fp) - 1;
122     fclose(fp);
123     if (inpfile[i].isDSF < 0) {
124 greg 2.1 fputs(fname, stderr);
125     fputs(": unknown format\n", stderr);
126     return(0);
127     }
128 greg 2.2 if ((inpfile[i].theta < -10000.) | (inpfile[i].phi < -10000.)) {
129     fputs(fname, stderr);
130     fputs(": unknown incident angle\n", stderr);
131     return(0);
132     }
133 greg 2.29 if (rev_orient) { /* reverse Z-axis to face outside */
134     inpfile[i].theta = 180. - inpfile[i].theta;
135     inpfile[i].phi = 360. - inpfile[i].phi;
136     }
137 greg 2.27 /* convert to Y-up orientation */
138     inpfile[i].phi += 90.-inpfile[i].up_phi;
139 greg 2.20 /* convert angle to grid position */
140 greg 2.37 set_grid_pos(&inpfile[i]);
141 greg 2.2 return(1);
142     }
143    
144     /* Load a set of measurements corresponding to a particular incident angle */
145     static int
146     add_pabopto_inp(const int i)
147     {
148 greg 2.19 FILE *fp = fopen(inpfile[i].fname, "r");
149 greg 2.28 double theta_out, phi_out, val[3];
150 greg 2.19 int n, c;
151 greg 2.2
152     if (fp == NULL || fseek(fp, inpfile[i].dstart, 0) == EOF) {
153     fputs(inpfile[i].fname, stderr);
154     fputs(": cannot open\n", stderr);
155     return(0);
156     }
157 greg 2.1 /* prepare input grid */
158 greg 2.20 if (!i || cmp_indir(&inpfile[i-1], &inpfile[i])) {
159     if (i) /* process previous incidence */
160 greg 2.2 make_rbfrep();
161     #ifdef DEBUG
162 greg 2.22 fprintf(stderr, "New incident (theta,phi)=(%.1f,%.1f)\n",
163 greg 2.2 inpfile[i].theta, inpfile[i].phi);
164     #endif
165 greg 2.28 if (inpfile[i].nspec)
166     set_spectral_samples(inpfile[i].nspec);
167 greg 2.2 new_bsdf_data(inpfile[i].theta, inpfile[i].phi);
168     }
169 greg 2.3 #ifdef DEBUG
170     fprintf(stderr, "Loading measurements from '%s'...\n", inpfile[i].fname);
171     #endif
172 greg 2.2 /* read scattering data */
173 greg 2.28 while (fscanf(fp, "%lf %lf %lf", &theta_out, &phi_out, val) == 3) {
174     for (n = 1; n < inpfile[i].nspec; n++)
175     if (fscanf(fp, "%lf", val+n) != 1) {
176     fprintf(stderr, "%s: warning: unexpected EOF\n",
177     inpfile[i].fname);
178     fclose(fp);
179     return(1);
180     }
181 greg 2.29 if (rev_orient) { /* reverse Z-axis to face outside */
182     theta_out = 180. - theta_out;
183     phi_out = 360. - phi_out;
184     }
185 greg 2.37 phi_out += 90.-inpfile[i].up_phi;
186     add_bsdf_data(theta_out, phi_out, val, inpfile[i].isDSF);
187 greg 2.28 }
188 greg 2.1 n = 0;
189     while ((c = getc(fp)) != EOF)
190     n += !isspace(c);
191     if (n)
192     fprintf(stderr,
193     "%s: warning: %d unexpected characters past EOD\n",
194 greg 2.2 inpfile[i].fname, n);
195 greg 2.1 fclose(fp);
196     return(1);
197     }
198    
199 greg 2.8 #ifndef TEST_MAIN
200 greg 2.35
201     #define SYM_ILL '?' /* illegal symmetry value */
202     #define SYM_ISO 'I' /* isotropic */
203     #define SYM_QUAD 'Q' /* quadrilateral symmetry */
204     #define SYM_BILAT 'B' /* bilateral symmetry */
205     #define SYM_ANISO 'A' /* anisotropic */
206 greg 2.37 #define SYM_UP 'U' /* "up-down" (180°) symmetry */
207 greg 2.35
208     static const char quadrant_rep[16][16] = {
209     "in-plane","0-90","90-180","0-180",
210     "180-270","0-90+180-270","90-270",
211     "0-270","270-360","270-90",
212     "90-180+270-360","270-180","180-360",
213     "180-90","90-360","0-360"
214     };
215     static const char quadrant_sym[16] = {
216     SYM_ISO, SYM_QUAD, SYM_QUAD, SYM_BILAT,
217     SYM_QUAD, SYM_ILL, SYM_BILAT, SYM_ILL,
218     SYM_QUAD, SYM_BILAT, SYM_ILL, SYM_ILL,
219     SYM_BILAT, SYM_ILL, SYM_ILL, SYM_ANISO
220     };
221    
222 greg 2.1 /* Read in PAB-Opto BSDF files and output RBF interpolant */
223     int
224     main(int argc, char *argv[])
225     {
226 greg 2.35 extern int nprocs;
227 greg 2.37 const char *symmetry = "0";
228     int ninpfiles, totinc;
229     int a, i;
230 greg 2.1 /* start header */
231     SET_FILE_BINARY(stdout);
232     newheader("RADIANCE", stdout);
233     printargs(argc, argv, stdout);
234     fputnow(stdout);
235     progname = argv[0]; /* get options */
236 greg 2.37 for (a = 1; a < argc && argv[a][0] == '-'; a++)
237     switch (argv[a][1]) {
238 greg 2.29 case 't':
239     rev_orient = !rev_orient;
240     break;
241 greg 2.1 case 'n':
242 greg 2.37 nprocs = atoi(argv[++a]);
243 greg 2.1 break;
244 greg 2.33 case 's':
245 greg 2.37 symmetry = argv[++a];
246 greg 2.33 break;
247 greg 2.1 default:
248     goto userr;
249     }
250 greg 2.37 totinc = ninpfiles = argc - a; /* initialize & sort inputs */
251 greg 2.2 if (ninpfiles < 2)
252 greg 2.1 goto userr;
253 greg 2.37 if (toupper(symmetry[0]) == SYM_UP) /* special case for "up" symmetry */
254     totinc += ninpfiles;
255     inpfile = (PGINPUT *)malloc(sizeof(PGINPUT)*totinc);
256 greg 2.2 if (inpfile == NULL)
257     return(1);
258     for (i = 0; i < ninpfiles; i++)
259 greg 2.37 if (!init_pabopto_inp(i, argv[a+i]))
260 greg 2.1 return(1);
261 greg 2.37
262     for (i = ninpfiles; i < totinc; i++) { /* copy for "up" symmetry */
263     inpfile[i] = inpfile[i-ninpfiles];
264     inpfile[i].phi += 180.; /* invert duplicate data */
265     inpfile[i].up_phi -= 180.;
266     set_grid_pos(&inpfile[i]); /* grid location for sorting */
267     }
268 greg 2.25 qsort(inpfile, ninpfiles, sizeof(PGINPUT), cmp_indir);
269 greg 2.2 /* compile measurements */
270 greg 2.37 for (i = 0; i < totinc; i++)
271 greg 2.2 if (!add_pabopto_inp(i))
272     return(1);
273     make_rbfrep(); /* process last data set */
274 greg 2.33 /* check input symmetry */
275     switch (toupper(symmetry[0])) {
276 greg 2.37 case '0': /* unspecified symmetry */
277 greg 2.35 if (quadrant_sym[inp_coverage] != SYM_ILL)
278     break; /* anything legal goes */
279     fprintf(stderr, "%s: unsupported phi coverage (%s)\n",
280     progname, quadrant_rep[inp_coverage]);
281     return(1);
282 greg 2.37 case SYM_UP: /* faux "up" symmetry */
283     if (quadrant_sym[inp_coverage] == SYM_ANISO)
284     break;
285     /* fall through */
286     case SYM_ISO: /* usual symmetry types */
287 greg 2.35 case SYM_QUAD:
288     case SYM_BILAT:
289     case SYM_ANISO:
290     if (quadrant_sym[inp_coverage] == toupper(symmetry[0]))
291     break; /* matches spec */
292     fprintf(stderr,
293     "%s: phi coverage (%s) does not match requested '%s' symmetry\n",
294     progname, quadrant_rep[inp_coverage], symmetry);
295     return(1);
296 greg 2.33 default:
297     fprintf(stderr,
298 greg 2.37 "%s: -s option must be Isotropic, Quadrilateral, Bilateral, Up, or Anisotropic\n",
299 greg 2.33 progname);
300     return(1);
301     }
302     #ifdef DEBUG
303 greg 2.35 fprintf(stderr, "Input phi coverage (%s) has '%c' symmetry\n",
304     quadrant_rep[inp_coverage],
305     quadrant_sym[inp_coverage]);
306 greg 2.33 #endif
307 greg 2.1 build_mesh(); /* create interpolation */
308     save_bsdf_rep(stdout); /* write it out */
309     return(0);
310     userr:
311 greg 2.33 fprintf(stderr, "Usage: %s [-t][-n nproc][-s symmetry] meas1.dat meas2.dat .. > bsdf.sir\n",
312 greg 2.1 progname);
313     return(1);
314     }
315 greg 2.35
316     #else /* TEST_MAIN */
317    
318 greg 2.6 /* Test main produces a Radiance model from the given input file */
319     int
320     main(int argc, char *argv[])
321     {
322     PGINPUT pginp;
323     char buf[128];
324     FILE *pfp;
325     double bsdf, min_log;
326     FVECT dir;
327     int i, j, n;
328    
329 greg 2.21 progname = argv[0];
330 greg 2.6 if (argc != 2) {
331 greg 2.21 fprintf(stderr, "Usage: %s input.dat > output.rad\n", progname);
332 greg 2.6 return(1);
333     }
334     inpfile = &pginp;
335     if (!init_pabopto_inp(0, argv[1]) || !add_pabopto_inp(0))
336     return(1);
337     /* reduce data set */
338 greg 2.21 if (make_rbfrep() == NULL) {
339     fprintf(stderr, "%s: nothing to plot!\n", progname);
340     exit(1);
341     }
342 greg 2.9 #ifdef DEBUG
343     fprintf(stderr, "Minimum BSDF = %.4f\n", bsdf_min);
344     #endif
345 greg 2.13 min_log = log(bsdf_min*.5 + 1e-5);
346 greg 2.11 #if 1 /* produce spheres at meas. */
347 greg 2.6 puts("void plastic yellow\n0\n0\n5 .6 .4 .01 .04 .08\n");
348     n = 0;
349 greg 2.28 for (i = 0; i < grid_res; i++)
350     for (j = 0; j < grid_res; j++)
351 greg 2.16 if (dsf_grid[i][j].sum.n > 0) {
352 greg 2.6 ovec_from_pos(dir, i, j);
353 greg 2.16 bsdf = dsf_grid[i][j].sum.v /
354 greg 2.24 ((double)dsf_grid[i][j].sum.n*output_orient*dir[2]);
355 greg 2.6 if (bsdf <= bsdf_min*.6)
356     continue;
357 greg 2.14 bsdf = log(bsdf + 1e-5) - min_log;
358 greg 2.7 ovec_from_pos(dir, i, j);
359     printf("yellow sphere s%04d\n0\n0\n", ++n);
360     printf("4 %.6g %.6g %.6g %.6g\n\n",
361 greg 2.6 dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf,
362 greg 2.7 .007*bsdf);
363 greg 2.6 }
364 greg 2.11 #endif
365     #if 1 /* spheres at RBF peaks */
366     puts("void plastic red\n0\n0\n5 .8 .01 .01 .04 .08\n");
367     for (n = 0; n < dsf_list->nrbf; n++) {
368     RBFVAL *rbf = &dsf_list->rbfa[n];
369     ovec_from_pos(dir, rbf->gx, rbf->gy);
370 greg 2.26 bsdf = eval_rbfrep(dsf_list, dir);
371 greg 2.14 bsdf = log(bsdf + 1e-5) - min_log;
372 greg 2.11 printf("red sphere p%04d\n0\n0\n", ++n);
373     printf("4 %.6g %.6g %.6g %.6g\n\n",
374     dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf,
375     .011*bsdf);
376     }
377     #endif
378     #if 1 /* output continuous surface */
379 greg 2.6 puts("void trans tgreen\n0\n0\n7 .7 1 .7 .04 .04 .9 1\n");
380     fflush(stdout);
381 greg 2.28 sprintf(buf, "gensurf tgreen bsdf - - - %d %d", grid_res-1, grid_res-1);
382 greg 2.6 pfp = popen(buf, "w");
383     if (pfp == NULL) {
384 greg 2.21 fprintf(stderr, "%s: cannot open '| %s'\n", progname, buf);
385 greg 2.6 return(1);
386     }
387 greg 2.28 for (i = 0; i < grid_res; i++)
388     for (j = 0; j < grid_res; j++) {
389 greg 2.6 ovec_from_pos(dir, i, j);
390 greg 2.26 bsdf = eval_rbfrep(dsf_list, dir);
391 greg 2.14 bsdf = log(bsdf + 1e-5) - min_log;
392 greg 2.6 fprintf(pfp, "%.8e %.8e %.8e\n",
393     dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
394     }
395 greg 2.11 if (pclose(pfp) != 0)
396     return(1);
397     #endif
398     return(0);
399 greg 2.6 }
400 greg 2.35
401     #endif /* TEST_MAIN */