ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/pabopto2xyz.c
Revision: 2.1
Committed: Wed Jul 27 18:05:53 2016 UTC (7 years, 9 months ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
Added pabopto2xyz utility to interpolate CIE-XYZ color values from measurements

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: pabopto2bsdf.c,v 2.29 2016/02/02 00:42:11 greg Exp $";
3 #endif
4 /*
5 * Combine PAB-Opto data files for color (CIE-XYZ) interpolation.
6 *
7 * G. Ward
8 */
9
10 #define _USE_MATH_DEFINES
11 #include "rtio.h"
12 #include <stdlib.h>
13 #include <ctype.h>
14 #include <math.h>
15 #include "interp2d.h"
16
17 /* #define DEBUG 1 */
18
19 #ifndef MAX_INPUTS
20 #define MAX_INPUTS 512 /* can handle this many files per chan */
21 #endif
22
23 #define ANGRES 256 /* resolution per 180 degrees */
24
25 /* Input file header information */
26 typedef struct {
27 const char *fname; /* input file path */
28 double theta, phi; /* input angles (degrees) */
29 double up_phi; /* azimuth for "up" direction */
30 int igp[2]; /* input grid position */
31 int isDSF; /* data is DSF (rather than BSDF)? */
32 int ndata; /* number of data points (or 0) */
33 long dstart; /* data start offset in file */
34 } PGINPUT;
35
36 /* 2-D interpolant for a particular incident direction */
37 typedef struct {
38 INTERP2 *ip2; /* 2-D interpolatiing function */
39 float *va; /* corresponding value array */
40 } PGINTERP;
41
42 extern void SDdisk2square(double sq[2], double diskx, double disky);
43
44 char *progname; /* global argv[0] */
45
46 int nprocs = 1; /* number of processes to run */
47
48 int nchild = 0; /* number of children (-1 in child) */
49
50 char *basename = "pabopto_xyz";
51
52 /* color conversion matrix */
53 double sens2xyz[3][3] = {
54 8.7539e-01, 3.5402e-02, 9.0381e-02,
55 1.3463e+00, -3.9082e-01, 7.2721e-02,
56 2.7306e-01, -2.1111e-01, 1.1014e+00,
57 };
58
59 #if defined(_WIN32) || defined(_WIN64)
60
61 #define await_children(n) (void)(n)
62 #define run_subprocess() 0
63 #define end_subprocess() (void)0
64
65 #else /* ! _WIN32 */
66
67 #include <unistd.h>
68 #include <sys/wait.h>
69
70 /* Wait for the specified number of child processes to complete */
71 static void
72 await_children(int n)
73 {
74 int exit_status = 0;
75
76 if (n > nchild)
77 n = nchild;
78 while (n-- > 0) {
79 int status;
80 if (wait(&status) < 0) {
81 fprintf(stderr, "%s: missing child(ren)!\n", progname);
82 nchild = 0;
83 break;
84 }
85 --nchild;
86 if (status) { /* something wrong */
87 if ((status = WEXITSTATUS(status)))
88 exit_status = status;
89 else
90 exit_status += !exit_status;
91 fprintf(stderr, "%s: subprocess died\n", progname);
92 n = nchild; /* wait for the rest */
93 }
94 }
95 if (exit_status)
96 exit(exit_status);
97 }
98
99 /* Start child process if multiprocessing selected */
100 static pid_t
101 run_subprocess(void)
102 {
103 int status;
104 pid_t pid;
105
106 if (nprocs <= 1) /* any children requested? */
107 return(0);
108 await_children(nchild + 1 - nprocs); /* free up child process */
109 if ((pid = fork())) {
110 if (pid < 0) {
111 fprintf(stderr, "%s: cannot fork subprocess\n",
112 progname);
113 await_children(nchild);
114 exit(1);
115 }
116 ++nchild; /* subprocess started */
117 return(pid);
118 }
119 nchild = -1;
120 return(0); /* put child to work */
121 }
122
123 /* If we are in subprocess, call exit */
124 #define end_subprocess() if (nchild < 0) _exit(0); else
125
126 #endif /* ! _WIN32 */
127
128 /* Compute square location from normalized input/output vector */
129 static void
130 sq_from_ang(double sq[2], double theta, double phi)
131 {
132 double vec[3];
133 double norm;
134 /* uniform hemispherical projection */
135 vec[2] = sin(M_PI/180.*theta);
136 vec[0] = cos(M_PI/180.*phi)*vec[2];
137 vec[1] = sin(M_PI/180.*phi)*vec[2];
138 vec[2] = sqrt(1. - vec[2]*vec[2]);
139 norm = 1./sqrt(1. + vec[2]);
140
141 SDdisk2square(sq, vec[0]*norm, vec[1]*norm);
142 }
143
144 /* Compute quantized grid position from normalized input/output vector */
145 static void
146 pos_from_ang(int gp[2], double theta, double phi)
147 {
148 double sq[2];
149
150 sq_from_ang(sq, theta, phi);
151 gp[0] = (int)(sq[0]*ANGRES);
152 gp[1] = (int)(sq[1]*ANGRES);
153 }
154
155 /* Compare incident angles */
156 static int
157 cmp_indir(const void *p1, const void *p2)
158 {
159 const PGINPUT *inp1 = (const PGINPUT *)p1;
160 const PGINPUT *inp2 = (const PGINPUT *)p2;
161 int ydif;
162
163 if (!inp1->dstart) /* put terminator last */
164 return(inp2->dstart > 0);
165 if (!inp2->dstart)
166 return(-1);
167
168 ydif = inp1->igp[1] - inp2->igp[1];
169 if (ydif)
170 return(ydif);
171
172 return(inp1->igp[0] - inp2->igp[0]);
173 }
174
175 /* Prepare a PAB-Opto input file by reading its header */
176 static int
177 init_pabopto_inp(PGINPUT *pgi, const char *fname)
178 {
179 FILE *fp = fopen(fname, "r");
180 char buf[2048];
181 int c;
182
183 if (fp == NULL) {
184 fputs(fname, stderr);
185 fputs(": cannot open\n", stderr);
186 return(0);
187 }
188 pgi->fname = fname;
189 pgi->isDSF = -1;
190 pgi->ndata = 0;
191 pgi->up_phi = 0;
192 pgi->theta = pgi->phi = -10001.;
193 /* read header information */
194 while ((c = getc(fp)) == '#' || c == EOF) {
195 char typ[64];
196 if (fgets(buf, sizeof(buf), fp) == NULL) {
197 fputs(fname, stderr);
198 fputs(": unexpected EOF\n", stderr);
199 fclose(fp);
200 return(0);
201 }
202 if (sscanf(buf, "colorimetry: %s", typ) == 1) {
203 if (!strcasecmp(typ, "CIE-XYZ")) {
204 fputs(fname, stderr);
205 fputs(": already in XYZ color space!\n", stderr);
206 return(0);
207 }
208 continue;
209 }
210 if (sscanf(buf, "datapoints_in_file:%d", &pgi->ndata) == 1)
211 continue;
212 if (sscanf(buf, "format: theta phi %s", typ) == 1) {
213 if (!strcasecmp(typ, "DSF"))
214 pgi->isDSF = 1;
215 else if (!strcasecmp(typ, "BSDF") ||
216 !strcasecmp(typ, "BRDF") ||
217 !strcasecmp(typ, "BTDF"))
218 pgi->isDSF = 0;
219 continue;
220 }
221 if (sscanf(buf, "upphi %lf", &pgi->up_phi) == 1)
222 continue;
223 if (sscanf(buf, "intheta %lf", &pgi->theta) == 1)
224 continue;
225 if (sscanf(buf, "inphi %lf", &pgi->phi) == 1)
226 continue;
227 if (sscanf(buf, "incident_angle %lf %lf",
228 &pgi->theta, &pgi->phi) == 2)
229 continue;
230 }
231 pgi->dstart = ftell(fp) - 1;
232 fclose(fp);
233 if (pgi->isDSF < 0) {
234 fputs(fname, stderr);
235 fputs(": unknown format\n", stderr);
236 return(0);
237 }
238 if ((pgi->theta < -10000.) | (pgi->phi < -10000.)) {
239 fputs(fname, stderr);
240 fputs(": unknown incident angle\n", stderr);
241 return(0);
242 }
243 /* convert angle to grid position */
244 pos_from_ang(pgi->igp, pgi->theta, pgi->phi);
245 return(1);
246 }
247
248 /* Load interpolating function for a particular incident direction */
249 static int
250 load_interp(PGINTERP *pgint, const int igp[2], const PGINPUT *pginp)
251 {
252 int nv = 0;
253 int nread = 0;
254 pgint->ip2 = NULL;
255 pgint->va = NULL;
256 /* load input file(s) */
257 while (pginp->dstart && (pginp->igp[0] == igp[0]) &
258 (pginp->igp[1] == igp[1])) {
259 FILE *fp = fopen(pginp->fname, "r");
260 double theta, phi;
261 float val;
262
263 if (fp == NULL || fseek(fp, pginp->dstart, SEEK_SET) < 0) {
264 fputs(pginp->fname, stderr);
265 fputs(": cannot re-open input file!\n", stderr);
266 return(0);
267 }
268 #ifdef DEBUG
269 fprintf(stderr, "Loading channel from '%s'\n", pginp->fname);
270 #endif
271 while (fscanf(fp, "%lf %lf %f", &theta, &phi, &val) == 3) {
272 double sq[2];
273 if (nread >= nv) {
274 if (pginp->ndata > 0)
275 nv += pginp->ndata;
276 else
277 nv += (nv>>1) + 1024;
278 pgint->ip2 = interp2_realloc(pgint->ip2, nv);
279 pgint->va = (float *)realloc(pgint->va, sizeof(float)*nv);
280 if ((pgint->ip2 == NULL) | (pgint->va == NULL))
281 goto memerr;
282 }
283 sq_from_ang(sq, theta, phi);
284 pgint->ip2->spt[nread][0] = sq[0]*ANGRES;
285 pgint->ip2->spt[nread][1] = sq[1]*ANGRES;
286 pgint->va[nread++] = val;
287 }
288 fclose(fp);
289 ++pginp; /* advance to next input */
290 }
291 if (nv > nread) { /* fix array sizes */
292 pgint->ip2 = interp2_realloc(pgint->ip2, nread);
293 pgint->va = (float *)realloc(pgint->va, nread);
294 if ((pgint->ip2 == NULL) | (pgint->va == NULL))
295 goto memerr;
296 }
297 return(1);
298 memerr:
299 fputs(progname, stderr);
300 fputs(": Out of memory in load_interp()!\n", stderr);
301 return(0);
302 }
303
304 /* Interpolate value at the given position */
305 static double
306 interp2val(const PGINTERP *pgint, double px, double py)
307 {
308 #define NSMP 12
309 float wt[NSMP];
310 int si[NSMP];
311 int n = interp2_topsamp(wt, si, NSMP, pgint->ip2, px, py);
312 double v = .0;
313
314 while (n-- > 0)
315 v += wt[n]*pgint->va[si[n]];
316
317 return(v);
318 #undef NSMP
319 }
320
321 static void
322 free_interp(PGINTERP *pgint)
323 {
324 interp2_free(pgint->ip2);
325 pgint->ip2 = NULL;
326 free(pgint->va);
327 pgint->va = NULL;
328 }
329
330 /* Interpolate the given incident direction to an output file */
331 static int
332 interp_xyz(const PGINPUT *inp0, int nf, const PGINTERP *int1, const PGINTERP *int2)
333 {
334 int dtotal = 0;
335 char outfname[256];
336 FILE *ofp;
337 int i;
338
339 if (nf <= 0)
340 return(0);
341 /* create output file & write header */
342 sprintf(outfname, "%s_t%03.0fp%03.0f.txt", basename, inp0->theta, inp0->phi);
343 if ((ofp = fopen(outfname, "w")) == NULL) {
344 fputs(outfname, stderr);
345 fputs(": cannot open for writing\n", stderr);
346 return(0);
347 }
348 fprintf(ofp, "#data written using %s\n", progname);
349 fprintf(ofp, "#incident_angle %g %g\n", inp0->theta, inp0->phi);
350 fprintf(ofp, "#intheta %g\n", inp0->theta);
351 fprintf(ofp, "#inphi %g\n", inp0->phi);
352 fprintf(ofp, "#upphi %g\n", inp0->up_phi);
353 fprintf(ofp, "#colorimetry: CIE-XYZ\n");
354 for (i = nf; i--; )
355 dtotal += inp0[i].ndata;
356 if (dtotal > 0)
357 fprintf(ofp, "#datapoints_in_file: %d\n", dtotal);
358 fprintf(ofp, "#format: theta phi %s\n", inp0->isDSF ? "DSF" : "BSDF");
359 dtotal = 0;
360 while (nf-- > 0) { /* load channel 1 file(s) */
361 FILE *ifp = fopen(inp0->fname, "r");
362 double theta, phi, val[3], sq[2];
363
364 if (ifp == NULL || fseek(ifp, inp0->dstart, SEEK_SET) < 0) {
365 fputs(inp0->fname, stderr);
366 fputs(": cannot re-open input file!\n", stderr);
367 return(0);
368 }
369 #ifdef DEBUG
370 fprintf(stderr, "Adding points from '%s' to '%s'\n",
371 inp0->fname, outfname);
372 #endif
373 while (fscanf(ifp, "%lf %lf %lf", &theta, &phi, &val[0]) == 3) {
374 sq_from_ang(sq, theta, phi);
375 val[1] = interp2val(int1, sq[0]*ANGRES, sq[1]*ANGRES);
376 val[2] = interp2val(int2, sq[0]*ANGRES, sq[1]*ANGRES);
377 fprintf(ofp, "%.3f %.3f", theta, phi);
378 for (i = 0; i < 3; i++)
379 fprintf(ofp, " %.4f",
380 sens2xyz[i][0]*val[0] +
381 sens2xyz[i][1]*val[1] +
382 sens2xyz[i][2]*val[2]);
383 fputc('\n', ofp);
384 ++dtotal;
385 }
386 fclose(ifp);
387 ++inp0; /* advance to next input file */
388 }
389 #ifdef DEBUG
390 fprintf(stderr, "Wrote %d values to '%s'\n", dtotal, outfname);
391 #endif
392 return(fclose(ofp) == 0);
393 }
394
395 /* Do next incident direction and advance counters */
396 static int
397 advance_incidence(PGINPUT *slist[3], int ndx[3])
398 {
399 const int i0 = ndx[0];
400 PGINTERP pgi1, pgi2;
401 int c;
402 /* match up directions */
403 while ((c = cmp_indir(slist[1]+ndx[1], slist[0]+ndx[0])) < 0)
404 ndx[1]++;
405 if (c) {
406 fputs(slist[0][ndx[0]].fname, stderr);
407 fputs(": warning - missing channel 2 at this incidence\n", stderr);
408 do
409 ++ndx[0];
410 while (!cmp_indir(slist[0]+ndx[0], slist[0]+i0));
411 return(0);
412 }
413 while ((c = cmp_indir(slist[2]+ndx[2], slist[0]+ndx[0])) < 0)
414 ndx[2]++;
415 if (c) {
416 fputs(slist[0][ndx[0]].fname, stderr);
417 fputs(": warning - missing channel 3 at this incidence\n", stderr);
418 do
419 ++ndx[0];
420 while (!cmp_indir(slist[0]+ndx[0], slist[0]+i0));
421 return(0);
422 }
423 do /* advance first channel index */
424 ++ndx[0];
425 while (!cmp_indir(slist[0]+ndx[0], slist[0]+i0));
426
427 if (run_subprocess()) /* fork here if multiprocessing */
428 return(1);
429 /* interpolate channels 2 & 3 */
430 load_interp(&pgi1, slist[0][i0].igp, slist[1]+ndx[1]);
431 load_interp(&pgi2, slist[0][i0].igp, slist[2]+ndx[2]);
432 if (!interp_xyz(slist[0]+i0, ndx[0]-i0, &pgi1, &pgi2))
433 return(-1);
434 end_subprocess();
435 free_interp(&pgi1);
436 free_interp(&pgi2);
437 return(1);
438 }
439
440 /* Read in single-channel PAB-Opto BSDF files and output XYZ versions */
441 int
442 main(int argc, char *argv[])
443 {
444 char *flist[MAX_INPUTS];
445 PGINPUT *slist[3];
446 int i, j;
447 int ndx[3];
448 /* get options */
449 progname = argv[0];
450 for (i = 1; i < argc-3 && argv[i][0] == '-'; i++)
451 switch (argv[i][1]) {
452 case 'n':
453 nprocs = atoi(argv[++i]);
454 break;
455 case 'm':
456 if (i >= argc-(3+9))
457 goto userr;
458 for (j = 0; j < 3; j++) {
459 sens2xyz[j][0] = atof(argv[++i]);
460 sens2xyz[j][1] = atof(argv[++i]);
461 sens2xyz[j][2] = atof(argv[++i]);
462 }
463 break;
464 case 'o':
465 basename = argv[++i];
466 break;
467 default:
468 goto userr;
469 }
470 if (i > argc-3)
471 goto userr;
472 for (j = 0; j < 3; j++) { /* prep input channels */
473 int k, n;
474 n = wordfile(flist, MAX_INPUTS, argv[i+j]);
475 if (n <= 0) {
476 fputs(argv[i+j], stderr);
477 fputs(": cannot load input file names\n", stderr);
478 return(1);
479 }
480 slist[j] = (PGINPUT *)malloc(sizeof(PGINPUT)*(n+1));
481 if (slist[j] == NULL) {
482 fputs(argv[0], stderr);
483 fputs(": out of memory!\n", stderr);
484 return(1);
485 }
486 #ifdef DEBUG
487 fprintf(stderr, "Checking %d files from '%s'\n", n, argv[i+j]);
488 #endif
489 for (k = 0; k < n; k++)
490 if (!init_pabopto_inp(slist[j]+k, flist[k]))
491 return(1);
492 memset(slist[j]+n, 0, sizeof(PGINPUT));
493 /* sort by incident direction */
494 qsort(slist[j], n, sizeof(PGINPUT), cmp_indir);
495 }
496 ndx[0]=ndx[1]=ndx[2]=0; /* compile measurements */
497 while (slist[0][ndx[0]].dstart)
498 if (advance_incidence(slist, ndx) < 0)
499 return(1);
500 await_children(nchild); /* wait for all to finish */
501 return(0);
502 userr:
503 fputs("Usage: ", stderr);
504 fputs(argv[0], stderr);
505 fputs(" [-m X1 X2 X3 Y1 Y2 Y3 Z1 Z2 Z3][-o basename][-n nprocs]", stderr);
506 fputs(" s1files.txt s2files.txt s3files.txt\n", stderr);
507 return(1);
508 }