ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/pvsum.c
Revision: 2.4
Committed: Fri Apr 18 23:59:03 2025 UTC (13 days, 22 hours ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 2.3: +2 -2 lines
Log Message:
refactor(rmtxop,rcomb,pvsum): Removed BSDF library dependency where unneeded

File Contents

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