ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/pvsum.c
Revision: 2.5
Committed: Fri Jun 6 18:26:22 2025 UTC (26 hours, 59 minutes ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 2.4: +2 -1 lines
Log Message:
fix(pvsum): Added missing include

File Contents

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