ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/pabopto2xyz.c
Revision: 2.7
Committed: Wed Dec 15 01:38:50 2021 UTC (2 years, 4 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R4, HEAD
Changes since 2.6: +5 -6 lines
Log Message:
refactor: removed prefix from SDdisk2square() and SDsquare2disk() & made consistent

File Contents

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