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 (8 weeks 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

# Content
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 }