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

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.7 static const char RCSid[] = "$Id: pabopto2xyz.c,v 2.6 2021/02/11 03:05:34 greg Exp $";
3 greg 2.1 #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 greg 2.7 #include "fvect.h"
13 greg 2.1 #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 greg 2.2 int incident_side = 0; /* signum(intheta) */
46     int exiting_side = 0; /* signum(outtheta) */
47    
48 greg 2.1 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 greg 2.7 sq_from_ang(RREAL sq[2], double theta, double phi)
133 greg 2.1 {
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 greg 2.7 disk2square(sq, vec[0]*norm, vec[1]*norm);
144 greg 2.1 }
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 greg 2.7 RREAL sq[2];
151 greg 2.1
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 greg 2.2 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 greg 2.1 /* 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 greg 2.2 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 greg 2.1 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 greg 2.3 #define NSMP 36
325 greg 2.1 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 greg 2.2 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 greg 2.1 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 greg 2.4 if (i != argc-3)
493 greg 2.1 goto userr;
494     for (j = 0; j < 3; j++) { /* prep input channels */
495 greg 2.5 char *flist[MAX_INPUTS];
496 greg 2.1 int k, n;
497     n = wordfile(flist, MAX_INPUTS, argv[i+j]);
498 greg 2.6 if ((n <= 0) | (n >= MAX_INPUTS-1)) {
499 greg 2.1 fputs(argv[i+j], stderr);
500 greg 2.6 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 greg 2.1 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     }