ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/pvsum.c
Revision: 2.1
Committed: Thu Mar 27 01:26:55 2025 UTC (7 weeks, 3 days ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
feat(pvsum): Added new pvsum tool for spectral data, similar to dctimestep

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2     static const char RCSid[] = "$Id$";
3     #endif
4     /*
5     * pvsum.c - add together spectral and/or float pictures
6     * based on vector or matrix, similar to dctimestep
7     */
8    
9     #include <math.h>
10     #include "rtio.h"
11     #include "resolu.h"
12     #include "platform.h"
13     #include "paths.h"
14     #include "rmatrix.h"
15     #if !defined(_WIN32) && !defined(_WIN64)
16     #include <sys/mman.h>
17     #endif
18    
19     int nprocs = 1; /* # of calculation processes (Unix) */
20     int in_type = DTfromHeader; /* input data type */
21     int out_type = DTfromHeader; /* output data type */
22     char *in_spec = NULL; /* input specification */
23     char *out_spec = NULL; /* output file specification */
24    
25     int iswapped = 0; /* input data is byte-swapped? */
26     int ncomp = 3; /* # input components */
27     int xres=0, yres=0; /* input image dimensions */
28    
29     RMATRIX *cmtx = NULL; /* coefficient matrix */
30    
31     /* does the given spec contain integer format? */
32     int
33     hasFormat(const char *s)
34     {
35     restart:
36     if (s) s = strchr(s, '%');
37     if (!s) return(0);
38     if (s[1] == '%') { /* "%%" ? */
39     s += 2;
40     goto restart;
41     }
42     while (*++s) {
43     if (strchr("diouxX", *s))
44     return(1); /* integer format */
45     if (strchr("%fFeEgGaAcsb", *s))
46     break; /* non-integer format */
47     }
48     return(0);
49     }
50    
51     /* get first input header data we'll need */
52     int
53     iheadline(char *s, void *p)
54     {
55     char fmt[MAXFMTLEN];
56     int i;
57    
58     if (!strncmp(s, "NCOLS=", 6)) {
59     xres = atoi(s+6);
60     return(1);
61     }
62     if (!strncmp(s, "NROWS=", 6)) {
63     yres = atoi(s+6);
64     return(1);
65     }
66     if (isncomp(s)) {
67     ncomp = ncompval(s);
68     return(1);
69     }
70     if (isexpos(s)) {
71     if (fabs(1. - exposval(s)) > 0.04)
72     fputs("Warning - ignoring EXPOSURE setting\n", stderr);
73     return(1);
74     }
75     if (formatval(fmt, s)) {
76     for (in_type = DTend; --in_type > DTfromHeader; )
77     if (!strcmp(fmt, cm_fmt_id[in_type]))
78     return(1);
79     return(-1);
80     }
81     i = isbigendian(s);
82     if (i >= 0) {
83     iswapped = (i != nativebigendian());
84     return(1);
85     }
86     return(0);
87     }
88    
89     /* open initial file and get relevant dimensions and type */
90     int
91     get_iotypes(void)
92     {
93     char fbuf[256];
94     FILE *fp;
95    
96     sprintf(fbuf, in_spec, 0);
97     fp = fopen(fbuf, "rb");
98     if (!fp) {
99     fprintf(stderr, "%s: cannot open for reading\n", fbuf);
100     return(0);
101     }
102     if (getheader(fp, iheadline, NULL) < 0) {
103     fprintf(stderr, "%s: bad header - wrong format?\n", fbuf);
104     fclose(fp);
105     return(0);
106     }
107     if ((xres <= 0) | (yres <= 0) && !fscnresolu(&xres, &yres, fp)) {
108     fprintf(stderr, "%s: missing input resolution\n", fbuf);
109     fclose(fp);
110     return(0);
111     }
112     if (nprocs > 1 && (in_type==DTrgbe) | (in_type==DTxyze)) {
113     long data_start = ftell(fp); /* make sure input flat */
114     off_t dend = lseek(fileno(fp), 0, SEEK_END);
115     if (dend < data_start + 4L*xres*yres) {
116     fputs("Warning - multi-processing requires flat input files\n",
117     stderr);
118     nprocs = 1;
119     }
120     }
121     fclose(fp);
122     if ((cmtx->ncomp == 1) & (ncomp != 1)) {
123     double xfm[MAXCSAMP];
124     RMATRIX *nmtx;
125     int i;
126     for (i = ncomp; i--; )
127     xfm[i] = 1.;
128     nmtx = rmx_transform(cmtx, ncomp, xfm);
129     if (!nmtx)
130     return(0);
131     rmx_free(cmtx);
132     cmtx = nmtx;
133     } else if (cmtx->ncomp != ncomp) {
134     fprintf(stderr, "operation %s needs %d components, has %d\n",
135     cmtx->ncols == 1 ? "vector" : "matrix",
136     ncomp, cmtx->ncomp);
137     return(0);
138     }
139     if ((in_type != DTrgbe) & (in_type != DTxyze) & (in_type != DTspec) &
140     (in_type != DTfloat)) {
141     fprintf(stderr, "%s unsupported input data type: %s\n",
142     fbuf, cm_fmt_id[in_type]);
143     return(0);
144     }
145     if ((out_type == DTrgbe) & (ncomp > 3))
146     out_type = DTspec;
147     else if (out_type == DTfromHeader ||
148     (out_type == DTrgbe) & (in_type != DTfloat))
149     out_type = in_type;
150     return(1);
151     }
152    
153     /* check subsequent headers match initial file */
154     int
155     checkline(char *s, void *p)
156     {
157     int *xyres = (int *)p;
158     char fmt[MAXFMTLEN];
159    
160     if (!strncmp(s, "NCOLS=", 6)) {
161     xyres[0] = atoi(s+6);
162     if (xyres[0] <= 0)
163     return(-1);
164     return(1);
165     }
166     if (!strncmp(s, "NROWS=", 6)) {
167     xyres[1] = atoi(s+6);
168     if (xyres[1] <= 0)
169     return(-1);
170     return(1);
171     }
172     if (isncomp(s)) {
173     if (ncompval(s) != ncomp)
174     return(-1);
175     return(1);
176     }
177     if (isexpos(s)) {
178     if (fabs(1. - exposval(s)) > 0.04)
179     fputs("Warning - ignoring EXPOSURE setting\n", stderr);
180     return(1);
181     }
182     if (formatval(fmt, s)) {
183     if (strcmp(fmt, cm_fmt_id[in_type]))
184     return(-1);
185     return(1);
186     }
187     return(0);
188     }
189    
190     /* open and check input file */
191     FILE *
192     open_input(char *fname)
193     {
194     int xyres[2];
195     FILE *fp = fopen(fname, "rb");
196    
197     if (!fp) {
198     fprintf(stderr, "%s: cannot open for reading\n", fname);
199     return(NULL);
200     }
201     xyres[0] = xyres[1] = 0;
202     if (getheader(fp, checkline, xyres) < 0) {
203     fprintf(stderr, "%s: bad/inconsistent header\n", fname);
204     fclose(fp);
205     return(NULL);
206     }
207     if ((xyres[0] <= 0) | (xyres[1] <= 0) &&
208     !fscnresolu(&xyres[0], &xyres[1], fp)) {
209     fprintf(stderr, "%s: missing resolution\n", fname);
210     fclose(fp);
211     return(NULL);
212     }
213     if ((xyres[0] != xres) | (xyres[1] != yres)) {
214     fprintf(stderr, "%s: mismatched resolution\n", fname);
215     fclose(fp);
216     return(NULL);
217     }
218     return(fp);
219     }
220    
221     /* open output file or command (if !NULL) and write info header */
222     FILE *
223     open_output(char *ospec, int fno)
224     {
225     FILE *fp;
226    
227     if (!ospec) {
228     ospec = "<stdout>";
229     fp = stdout;
230     SET_FILE_BINARY(fp);
231     } else if (ospec[0] == '!') {
232     if (!(fp = popen(ospec+1, "w"))) {
233     fprintf(stderr, "Cannot start: %s\n", ospec);
234     return(NULL);
235     }
236     SET_FILE_BINARY(fp);
237     } else if (!(fp = fopen(ospec, "wb"))) {
238     fprintf(stderr, "%s: cannot open for writing\n", ospec);
239     return(NULL);
240     }
241     newheader("RADIANCE", fp);
242     fputnow(fp);
243     if (fno >= 0)
244     fprintf(fp, "FRAME=%d\n", fno);
245     switch (out_type) {
246     case DTfloat:
247     case DTdouble:
248     fprintf(fp, "NCOLS=%d\nNROWS=%d\n", xres, yres);
249     fputncomp(ncomp, fp);
250     fputendian(fp);
251     fputformat(cm_fmt_id[out_type], fp);
252     fputc('\n', fp);
253     break;
254     case DTrgbe:
255     fputformat(COLRFMT, fp);
256     fputc('\n', fp);
257     fprtresolu(xres, yres, fp);
258     break;
259     case DTxyze:
260     fputformat(CIEFMT, fp);
261     fputc('\n', fp);
262     fprtresolu(xres, yres, fp);
263     break;
264     case DTspec:
265     fputncomp(ncomp, fp);
266     fputwlsplit(cmtx->wlpart, fp);
267     fputformat(SPECFMT, fp);
268     fputc('\n', fp);
269     fprtresolu(xres, yres, fp);
270     break;
271     default:
272     fputs("Unsupported output type!\n", stderr);
273     return(NULL);
274     }
275     if (fflush(fp) < 0) {
276     fprintf(stderr, "%s: write error\n", ospec);
277     fclose(fp);
278     return(NULL);
279     }
280     return(fp);
281     }
282    
283     /* run calculation from a single process */
284     int
285     solo_process(void)
286     {
287     float *osum = (float *)calloc((size_t)xres*yres, sizeof(float)*ncomp);
288     COLORV *iscan = (COLORV *)malloc(sizeof(COLORV)*ncomp*xres);
289     char fbuf[512];
290     int c;
291    
292     if (!osum | !iscan) {
293     fprintf(stderr, "Cannot allocate %dx%d %d-component accumulator\n",
294     xres, yres, ncomp);
295     return(0);
296     }
297     if (sizeof(float) != sizeof(COLORV)) {
298     fputs("Code Error 1 in solo_process()\n", stderr);
299     return(0);
300     }
301     for (c = 0; c < cmtx->ncols; c++) { /* run through each column/output */
302     FILE *fout;
303     int y;
304     int rc = cmtx->nrows;
305     if (c > 0) /* clear accumulator? */
306     memset(osum, 0, sizeof(float)*ncomp*xres*yres);
307     while (rc-- > 0) { /* run through each input file */
308     const int r = c&1 ? rc : cmtx->nrows-1 - rc;
309     const rmx_dtype *cval = rmx_val(cmtx, r, c);
310     FILE *finp;
311     int i, x;
312     for (i = ncomp; i--; )
313     if (cval[i] != 0) break;
314     if (i < 0) /* this coefficient is zero, skip */
315     continue;
316     sprintf(fbuf, in_spec, r);
317     finp = open_input(fbuf);
318     if (!finp)
319     return(0);
320     for (y = 0; y < yres; y++) {
321     float *dst = osum + y*xres*ncomp;
322     if (in_type == DTfloat ? getbinary(iscan, sizeof(float)*ncomp,
323     xres, finp) != xres
324     : freadsscan(iscan, ncomp, xres, finp) < 0)
325     goto readerr;
326     if ((in_type == DTfloat) & iswapped)
327     swap32((char *)iscan, ncomp*xres);
328     for (x = 0; x < xres; x++, dst += ncomp)
329     for (i = ncomp; i--; )
330     dst[i] += cval[i]*iscan[x*ncomp + i];
331     }
332     fclose(finp);
333     } /* write out accumulated column result... */
334     if (out_spec) { /* ...to file or command */
335     if (cmtx->ncols > 1 && !hasFormat(out_spec)) {
336     fputs("Sequential result must go to stdout\n", stderr);
337     return(0);
338     }
339     sprintf(fbuf, out_spec, c);
340     fout = open_output(fbuf, c-(cmtx->ncols==1));
341     } else { /* ...to stdout */
342     if ((out_type == DTfloat) & (cmtx->ncols > 1)) {
343     fputs("Float outputs must have separate destinations\n",
344     stderr);
345     return(0);
346     }
347     strcpy(fbuf, "<stdout>");
348     fout = open_output(NULL, c-(cmtx->ncols==1));
349     }
350     if (!fout)
351     return(0); /* assume error was reported */
352     if (out_type != DTfloat) {
353     for (y = 0; y < yres; y++)
354     if (fwritesscan(osum + (size_t)y*xres*ncomp,
355     ncomp, xres, fout) < 0)
356     goto writerr;
357     } else if (fwrite(osum, sizeof(float)*ncomp, (size_t)xres*yres, fout) !=
358     (size_t)xres*yres)
359     goto writerr;
360    
361     if (fbuf[0] == '!') {
362     if (pclose(fout) != 0) {
363     fprintf(stderr, "Bad status from: %s\n", fbuf);
364     return(0);
365     }
366     } else if (fout != stdout && fclose(fout) == EOF)
367     goto writerr;
368     }
369     free(osum); /* clean up on success */
370     free(iscan);
371     return(1);
372     readerr:
373     fprintf(stderr, "%s: read error\n", fbuf);
374     return(0);
375     writerr:
376     fprintf(stderr, "%s: write error\n", fbuf);
377     return(0);
378     }
379    
380     /* run calculation on multiple processes, using memory maps and fork() */
381     int
382     multi_process(void)
383     {
384     #if defined(_WIN32) || defined(_WIN64)
385     fputs("Bad call to multi_process()\n", stderr);
386     return(0);
387     #else
388     int coff = nprocs;
389     int odd = 0;
390     char fbuf[512];
391     float *osum;
392     int c;
393     /* sanity check */
394     if (sizeof(float) != sizeof(COLORV)) {
395     fputs("Code Error 1 in multi_process()\n", stderr);
396     return(0);
397     }
398     while (--coff > 0) { /* parent births children */
399     int pid = fork();
400     if (pid < 0) {
401     fputs("fork() call failed!\n", stderr);
402     return(0);
403     }
404     if (pid == 0) break; /* child gets to work */
405     }
406     osum = (float *)calloc((size_t)xres*yres, sizeof(float)*ncomp);
407     if (!osum) {
408     fprintf(stderr, "Cannot allocate %dx%d %d-component accumulator\n",
409     xres, yres, ncomp);
410     return(0);
411     }
412     /* run through our unique set of columns */
413     for (c = coff; c < cmtx->ncols; c += nprocs) {
414     FILE *fout;
415     int y;
416     int rc = cmtx->nrows;
417     if (c > coff) /* clear accumulator? */
418     memset(osum, 0, sizeof(float)*ncomp*xres*yres);
419     while (rc-- > 0) { /* map & sum each input file */
420     const int r = odd ? rc : cmtx->nrows-1 - rc;
421     const rmx_dtype *cval = rmx_val(cmtx, r, c);
422     long dstart, n;
423     size_t maplen;
424     void *imap;
425     FILE *finp;
426     float *dst;
427     int i;
428     for (i = ncomp; i--; )
429     if (cval[i] != 0) break;
430     if (i < 0) /* this coefficient is zero, skip */
431     continue;
432     sprintf(fbuf, in_spec, r);
433     finp = open_input(fbuf);
434     if (!finp)
435     return(0);
436     dstart = ftell(finp);
437     if (dstart < 0) {
438     fprintf(stderr, "%s: ftell() failed!\n", fbuf);
439     return(0);
440     }
441     if (in_type == DTfloat && dstart%sizeof(float)) {
442     fprintf(stderr, "%s: float header misalignment\n", fbuf);
443     return(0);
444     }
445     i = in_type==DTfloat ? ncomp*(int)sizeof(float) : ncomp+1;
446     maplen = dstart + yres*xres*i;
447     imap = mmap(NULL, maplen, PROT_READ,
448     MAP_FILE|MAP_SHARED, fileno(finp), 0);
449     fclose(finp); /* will load from map */
450     if (imap == MAP_FAILED) {
451     fprintf(stderr, "%s: unable to map input file\n", fbuf);
452     return(0);
453     }
454     dst = osum; /* -> weighted image sum */
455     if (in_type == DTfloat) {
456     const float *fvp = (float *)((char *)imap + dstart);
457     for (n = (long)yres*xres; n-- > 0;
458     dst += ncomp, fvp += ncomp)
459     for (i = ncomp; i--; )
460     dst[i] += cval[i]*fvp[i];
461     } else {
462     const COLRV *cvp = (COLRV *)((char *)imap + dstart);
463     for (n = (long)yres*xres; n-- > 0;
464     dst += ncomp, cvp += ncomp+1) {
465     const rmx_dtype fe = cxponent[cvp[ncomp]];
466     for (i = ncomp; i--; )
467     dst[i] += cval[i]*(cvp[i]+(rmx_dtype).5)*fe;
468     }
469     }
470     munmap(imap, maplen);
471     } /* write out accumulated column result */
472     sprintf(fbuf, out_spec, c);
473     fout = open_output(fbuf, c);
474     if (!fout)
475     return(0); /* assume error was reported */
476     if (out_type != DTfloat) {
477     for (y = 0; y < yres; y++)
478     if (fwritesscan(osum + (size_t)y*xres*ncomp,
479     ncomp, xres, fout) < 0)
480     goto writerr;
481     } else if (fwrite(osum, sizeof(float)*ncomp, (size_t)xres*yres, fout) !=
482     (size_t)xres*yres)
483     goto writerr;
484    
485     if (fbuf[0] == '!') {
486     if (pclose(fout) != 0) {
487     fprintf(stderr, "Bad status from: %s\n", fbuf);
488     return(0);
489     }
490     } else if (fclose(fout) == EOF)
491     goto writerr;
492     odd = !odd; /* go back & forth to milk page cache */
493     }
494     free(osum);
495     if (coff) /* children return here... */
496     return(1);
497     c = 0; /* ...but parent waits for children */
498     while (++coff < nprocs) {
499     int st;
500     if (wait(&st) < 0)
501     break;
502     if (st) c = st;
503     }
504     return(c == 0);
505     writerr:
506     fprintf(stderr, "%s: write error\n", fbuf);
507     return(0);
508     #endif
509     }
510    
511     int
512     main(int argc, char *argv[])
513     {
514     int a;
515    
516     for (a = 1; a < argc-1 && argv[a][0] == '-'; a++)
517     switch (argv[a][1]) {
518     case 'o': /* output spec/format */
519     switch (argv[a][2]) {
520     case '\0':
521     out_spec = argv[++a];
522     break;
523     case 'f':
524     out_type = DTfloat;
525     break;
526     case 'c':
527     out_type = DTrgbe;
528     break;
529     default:
530     goto badopt;
531     }
532     break;
533     case 'N': /* number of desired processes */
534     case 'n': /* quietly supported alternate */
535     nprocs = atoi(argv[++a]);
536     if (nprocs <= 0)
537     goto userr;
538     break;
539     default:;
540     badopt: fprintf(stderr, "%s: bad option: %s\n", argv[0], argv[a]);
541     goto userr;
542     return(1);
543     }
544     if ((argc-a < 1) | (argc-a > 2) || argv[a][0] == '-')
545     goto userr;
546     in_spec = argv[a];
547     cmtx = rmx_load(argv[a+1], RMPnone); /* may load from stdin */
548     if (cmtx == NULL)
549     return(1); /* error reported */
550     #if defined(_WIN32) || defined(_WIN64)
551     if (nprocs > 1) {
552     fprintf(stderr, "%s: warning - Windows only allows -N 1\n", argv[0]);
553     nprocs = 1;
554     }
555     #else
556     if (nprocs > cmtx->ncols)
557     nprocs = cmtx->ncols;
558     if ((nprocs > 1) & !out_spec) {
559     fprintf(stderr, "%s: multi-processing result cannot go to stdout\n",
560     argv[0]);
561     nprocs = 1;
562     }
563     if ((nprocs > 1) & iswapped && (in_type==DTfloat) | (in_type==DTdouble)) {
564     fprintf(stderr, "%s: multi-processing unsupported on swapped input\n",
565     argv[0]);
566     nprocs = 1;
567     }
568     #endif
569     if (cmtx->nrows > 1 && !hasFormat(in_spec)) {
570     fprintf(stderr, "%s: input specification '%s' needs %%d format\n",
571     argv[0], in_spec);
572     goto userr;
573     }
574     if (!get_iotypes())
575     return(1);
576     if (!(nprocs == 1 ? solo_process() : multi_process()))
577     return(1);
578     return(0);
579     userr:
580     fprintf(stderr, "Usage: %s [-oc | -of][-o ospec][-N nproc] inpspec [mtx]\n",
581     argv[0]);
582     return(1);
583     }