ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/pabopto2xyz.c
Revision: 2.5
Committed: Wed Jun 12 00:09:18 2019 UTC (4 years, 10 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R3
Changes since 2.4: +2 -2 lines
Log Message:
Minor optimization

File Contents

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