ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/Development/ray/src/util/pvsum.c
Revision: 2.8
Committed: Fri Oct 24 22:41:10 2025 UTC (7 days, 7 hours ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.7: +8 -1 lines
Log Message:
feat(pvsum): Add command arguments to output header(s)

File Contents

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