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

# User Rev Content
1 greg 2.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     }