ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/pabopto2xyz.c
Revision: 2.3
Committed: Mon Aug 22 21:03:00 2016 UTC (7 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R2, rad5R1
Changes since 2.2: +2 -2 lines
Log Message:
Increased number of samples during interpolation

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.3 static const char RCSid[] = "$Id: pabopto2xyz.c,v 2.2 2016/07/28 18:21:23 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     #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 greg 2.2 int incident_side = 0; /* signum(intheta) */
47     int exiting_side = 0; /* signum(outtheta) */
48    
49 greg 2.1 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 greg 2.2 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 greg 2.1 /* 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 greg 2.2 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 greg 2.1 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 greg 2.3 #define NSMP 36
326 greg 2.1 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 greg 2.2 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 greg 2.1 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     char *flist[MAX_INPUTS];
469     PGINPUT *slist[3];
470     int i, j;
471     int ndx[3];
472     /* get options */
473     progname = argv[0];
474     for (i = 1; i < argc-3 && argv[i][0] == '-'; i++)
475     switch (argv[i][1]) {
476     case 'n':
477     nprocs = atoi(argv[++i]);
478     break;
479     case 'm':
480     if (i >= argc-(3+9))
481     goto userr;
482     for (j = 0; j < 3; j++) {
483     sens2xyz[j][0] = atof(argv[++i]);
484     sens2xyz[j][1] = atof(argv[++i]);
485     sens2xyz[j][2] = atof(argv[++i]);
486     }
487     break;
488     case 'o':
489     basename = argv[++i];
490     break;
491     default:
492     goto userr;
493     }
494     if (i > argc-3)
495     goto userr;
496     for (j = 0; j < 3; j++) { /* prep input channels */
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     }