src/util/pvsum.c

#ifndef lint
static const char RCSid[] = "$Id: pvsum.c,v 2.7 2025/10/24 16:31:18 greg Exp $";
#endif
/*
 *      pvsum.c - add together spectral and/or float pictures
 *                      based on vector or matrix, similar to dctimestep
 */

#include <math.h>
#include "rtio.h"
#include "resolu.h"
#include "platform.h"
#include "random.h"
#include "rmatrix.h"
#if !defined(_WIN32) && !defined(_WIN64)
#include <sys/mman.h>
#include <sys/wait.h>
#endif

#define  VIEWSTR        "VIEW="         /* borrowed from common/view.h */
#define  VIEWSTRL       5

int     nprocs = 1;                     /* # of calculation processes (Unix) */
int     in_type = DTfromHeader;         /* input data type */
int     out_type = DTfromHeader;        /* output data type */
char    *in_spec = NULL;                /* input specification */
char    *out_spec = NULL;               /* output file specification */

int     iswapped = 0;                   /* input data is byte-swapped? */
int     ncomp = 3;                      /* # input components */
int     xres=0, yres=0;                 /* input image dimensions */
char    viewspec[128] = "";             /* VIEW= line from first header */
char    pixasp[48] = "";                /* PIXASPECT= line from header */

int     gargc;                          /* global argc */
char    **gargv;                        /* global argv */

RMATRIX *cmtx = NULL;                   /* coefficient matrix */

/* does the given spec contain integer format? */
int
hasFormat(const char *s)
{
restart:
        if (s) s = strchr(s, '%');
        if (!s) return(0);
        if (s[1] == '%') {              /* "%%" ? */
                s += 2;
                goto restart;
        }
        while (*++s) {
                if (strchr("diouxX", *s))
                        return(1);      /* integer format */
                if (strchr("%fFeEgGaAcsb", *s))
                        break;          /* non-integer format */
        }
        return(0);
}

/* get first input header data we'll need */
int
iheadline(char *s, void *p)
{
        char    fmt[MAXFMTLEN];
        int     i;

        if (!strncmp(s, "NCOLS=", 6)) {
                xres = atoi(s+6);
                return(1);
        }
        if (!strncmp(s, "NROWS=", 6)) {
                yres = atoi(s+6);
                return(1);
        }
        if (isncomp(s)) {
                ncomp = ncompval(s);
                return(1);
        }
        if (formatval(fmt, s)) {
                for (in_type = DTend; --in_type > DTfromHeader; )
                        if (!strcmp(fmt, cm_fmt_id[in_type]))
                                return(1);
                return(-1);
        }
        i = isbigendian(s);
        if (i >= 0) {
                iswapped = (i != nativebigendian());
                return(1);
        }
        if (!strncmp(s, VIEWSTR, VIEWSTRL)) {
                strcpy(viewspec, s);
                return(1);
        }
        if (isaspect(s)) {
                strcpy(pixasp, s);
                return(1);
        }
        return(0);
}

/* open initial file and get relevant dimensions and type */
int
get_iotypes(void)
{
        char    fbuf[256];
        FILE    *fp;

        sprintf(fbuf, in_spec, 0);
        fp = fopen(fbuf, "rb");
        if (!fp) {
                fprintf(stderr, "%s: cannot open for reading\n", fbuf);
                return(0);
        }
        if (getheader(fp, iheadline, NULL) < 0) {
                fprintf(stderr, "%s: bad header - wrong format?\n", fbuf);
                fclose(fp);
                return(0);
        }
        if ((xres <= 0) | (yres <= 0) && !fscnresolu(&xres, &yres, fp)) {
                fprintf(stderr, "%s: missing input resolution\n", fbuf);
                fclose(fp);
                return(0);
        }
        if (nprocs > 1 && (in_type==DTrgbe) | (in_type==DTxyze)) {
                long    data_start = ftell(fp);         /* make sure input flat */
                off_t   dend = lseek(fileno(fp), 0, SEEK_END);
                if (dend < data_start + 4L*xres*yres) {
                        fputs("Warning - multi-processing requires flat input files\n",
                                        stderr);
                        nprocs = 1;
                }
        }
        fclose(fp);
        if ((cmtx->ncomp == 1) & (ncomp != 1)) {
                double  xfm[MAXCSAMP];
                RMATRIX *nmtx;
                int     i;
                for (i = ncomp; i--; )
                        xfm[i] = 1.;
                nmtx = rmx_transform(cmtx, ncomp, xfm);
                if (!nmtx)
                        return(0);
                rmx_free(cmtx);
                cmtx = nmtx;
        } else if (cmtx->ncomp != ncomp) {
                fprintf(stderr, "operation %s needs %d components, has %d\n",
                                cmtx->ncols == 1 ? "vector" : "matrix",
                                ncomp, cmtx->ncomp);
                return(0);
        }
        if ((in_type != DTrgbe) & (in_type != DTxyze) & (in_type != DTspec) &
                        (in_type != DTfloat)) {
                fprintf(stderr, "%s unsupported input data type: %s\n",
                                fbuf, cm_fmt_id[in_type]);
                return(0);
        }
        if ((out_type == DTrgbe) & (ncomp > 3))
                out_type = DTspec;
        else if (out_type == DTfromHeader ||
                        (out_type == DTrgbe) & (in_type != DTfloat))
                out_type = in_type;
        return(1);
}

/* check subsequent headers match initial file */
int
checkline(char *s, void *p)
{
        static int      exposWarned = 0;
        int             *xyres = (int *)p;
        char            fmt[MAXFMTLEN];

        if (!strncmp(s, "NCOLS=", 6)) {
                xyres[0] = atoi(s+6);
                if (xyres[0] <= 0)
                        return(-1);
                return(1);
        }
        if (!strncmp(s, "NROWS=", 6)) {
                xyres[1] = atoi(s+6);
                if (xyres[1] <= 0)
                        return(-1);
                return(1);
        }
        if (isncomp(s)) {
                if (ncompval(s) != ncomp)
                        return(-1);
                return(1);
        }
        if (isexpos(s)) {
                if (!exposWarned && fabs(1. - exposval(s)) > 0.04) {
                        fputs("Warning - ignoring EXPOSURE setting(s)\n",
                                        stderr);
                        exposWarned++;
                }
                return(1);
        }
        if (formatval(fmt, s)) {
                if (strcmp(fmt, cm_fmt_id[in_type]))
                        return(-1);
                return(1);
        }
        return(0);
}

/* open and check input file */
FILE *
open_input(char *fname)
{
        int     xyres[2];
        FILE    *fp = fopen(fname, "rb");

        if (!fp) {
                fprintf(stderr, "%s: cannot open for reading\n", fname);
                return(NULL);
        }
        xyres[0] = xyres[1] = 0;
        if (getheader(fp, checkline, xyres) < 0) {
                fprintf(stderr, "%s: bad/inconsistent header\n", fname);
                fclose(fp);
                return(NULL);
        }
        if ((xyres[0] <= 0) | (xyres[1] <= 0) &&
                        !fscnresolu(&xyres[0], &xyres[1], fp)) {
                fprintf(stderr, "%s: missing resolution\n", fname);
                fclose(fp);
                return(NULL);
        }
        if ((xyres[0] != xres) | (xyres[1] != yres)) {
                fprintf(stderr, "%s: mismatched resolution\n", fname);
                fclose(fp);
                return(NULL);
        }
        return(fp);
}

/* open output file or command (if !NULL) and write info header */
FILE *
open_output(char *ospec, int fno)
{
        FILE    *fp;

        if (!ospec) {
                ospec = "<stdout>";
                fp = stdout;
        } else if (ospec[0] == '!') {
                if (!(fp = popen(ospec+1, "w"))) {
                        fprintf(stderr, "Cannot start: %s\n", ospec);
                        return(NULL);
                }
        } else if (!(fp = fopen(ospec, "w"))) {
                fprintf(stderr, "%s: cannot open for writing\n", ospec);
                return(NULL);
        }
        SET_FILE_BINARY(fp);
        newheader("RADIANCE", fp);
        if (cmtx->info)                 /* prepend matrix metadata */
                fputs(cmtx->info, fp);
        else
                fputnow(fp);
        printargs(gargc, gargv, fp);    /* this command */
        if (fno >= 0)
                fprintf(fp, "FRAME=%d\n", fno);
        if (viewspec[0])
                fputs(viewspec, fp);
        if (pixasp[0])
                fputs(pixasp, fp);
        switch (out_type) {
        case DTfloat:
        case DTdouble:
                fprintf(fp, "NCOLS=%d\nNROWS=%d\n", xres, yres);
                fputncomp(ncomp, fp);
                fputendian(fp);
                fputformat(cm_fmt_id[out_type], fp);
                fputc('\n', fp);
                break;
        case DTrgbe:
                fputformat(COLRFMT, fp);
                fputc('\n', fp);
                fprtresolu(xres, yres, fp);
                break;
        case DTxyze:
                fputformat(CIEFMT, fp);
                fputc('\n', fp);
                fprtresolu(xres, yres, fp);
                break;
        case DTspec:
                fputncomp(ncomp, fp);
                fputwlsplit(cmtx->wlpart, fp);
                fputformat(SPECFMT, fp);
                fputc('\n', fp);
                fprtresolu(xres, yres, fp);
                break;
        default:
                fputs("Unsupported output type!\n", stderr);
                return(NULL);
        }
        if (fflush(fp) < 0) {
                fprintf(stderr, "%s: write error\n", ospec);
                fclose(fp);
                return(NULL);
        }
        return(fp);
}

/* run calculation from a single process */
int
solo_process(void)
{
        float   *osum = (float *)calloc((size_t)xres*yres, sizeof(float)*ncomp);
        COLORV  *iscan = (COLORV *)malloc(sizeof(COLORV)*ncomp*xres);
        char    fbuf[512];
        int     c;

        if (!osum | !iscan) {
                fprintf(stderr, "Cannot allocate %dx%d %d-component accumulator\n",
                                xres, yres, ncomp);
                return(0);
        }
        if (sizeof(float) != sizeof(COLORV)) {
                fputs("Code Error 1 in solo_process()\n", stderr);
                return(0);
        }
        for (c = 0; c < cmtx->ncols; c++) {     /* run through each column/output */
                FILE    *fout;
                int     y;
                int     rc = cmtx->nrows;
                if (c > 0)              /* clear accumulator? */
                        memset(osum, 0, sizeof(float)*ncomp*xres*yres);
                while (rc-- > 0) {      /* run through each input file */
                        const int       r = c&1 ? rc : cmtx->nrows-1 - rc;
                        const rmx_dtype *cval = rmx_val(cmtx, r, c);
                        FILE            *finp;
                        int             i, x;
                        for (i = ncomp; i--; )
                                if (cval[i] != 0) break;
                        if (i < 0)      /* this coefficient is zero, skip */
                                continue;
                        sprintf(fbuf, in_spec, r);
                        finp = open_input(fbuf);
                        if (!finp)
                                return(0);
                        for (y = 0; y < yres; y++) {
                                float   *dst = osum + y*xres*ncomp;
                                if (in_type == DTfloat  ? getbinary(iscan, sizeof(float)*ncomp,
                                                                        xres, finp) != xres
                                                        : freadsscan(iscan, ncomp, xres, finp) < 0)
                                        goto readerr;
                                if ((in_type == DTfloat) & iswapped)
                                        swap32((char *)iscan, ncomp*xres);
                                for (x = 0; x < xres; x++, dst += ncomp)
                                        for (i = ncomp; i--; )
                                                dst[i] += cval[i]*iscan[x*ncomp + i];
                        }
                        fclose(finp);
                }                       /* write out accumulated column result... */
                if (out_spec) {         /* ...to file or command */
                        if (cmtx->ncols > 1 && !hasFormat(out_spec)) {
                                fputs("Sequential result must go to stdout\n", stderr);
                                return(0);
                        }
                        sprintf(fbuf, out_spec, c);
                        fout = open_output(fbuf, c-(cmtx->ncols==1));
                } else {                /* ...to stdout */
                        if ((out_type == DTfloat) & (cmtx->ncols > 1)) {
                                fputs("Float outputs must have separate destinations\n",
                                                stderr);
                                return(0);
                        }
                        strcpy(fbuf, "<stdout>");
                        fout = open_output(NULL, c-(cmtx->ncols==1));
                }
                if (!fout)
                        return(0);      /* assume error was reported */
                if (out_type != DTfloat) {
                        for (y = 0; y < yres; y++)
                                if (fwritesscan(osum + (size_t)y*xres*ncomp,
                                                        ncomp, xres, fout) < 0)
                                        goto writerr;
                } else if (fwrite(osum, sizeof(float)*ncomp, (size_t)xres*yres, fout) !=
                                        (size_t)xres*yres)
                        goto writerr;

                if (fbuf[0] == '!') {
                        if (pclose(fout) != 0) {
                                fprintf(stderr, "Bad status from: %s\n", fbuf);
                                return(0);
                        }
                } else if (fout != stdout && fclose(fout) == EOF)
                        goto writerr;
        }
        free(osum);                     /* clean up on success */
        free(iscan);
        return(1);
readerr:
        fprintf(stderr, "%s: read error\n", fbuf);
        return(0);
writerr:
        fprintf(stderr, "%s: write error\n", fbuf);
        return(0);
}

#if defined(_WIN32) || defined(_WIN64)
#define multi_process   solo_process
#else

/* allocate a scrambled index array of the specified length */
int *
scramble(int n)
{
        int     *scarr = (int *)malloc(sizeof(int)*n);
        int     i;

        if (!scarr) {
                fprintf(stderr, "Out of memory in scramble(%d)\n", n);
                exit(1);
        }
        for (i = n; i--; )
                scarr[i] = i;
                                        /* perform Fisher-Yates shuffle */
        for (i = 0; i < n-1; i++) {
                int     ix = irandom(n-i) + i;
                int     ndx = scarr[i];
                scarr[i] = scarr[ix];
                scarr[ix] = ndx;
        }
        return(scarr);
}

/* run calculation on multiple processes, using memory maps and fork() */
int
multi_process(void)
{
        int     coff = nprocs;
        int     odd = 0;
        char    fbuf[512];
        float   *osum;
        int     *syarr;
        int     c;
                                        /* sanity check */
        if (sizeof(float) != sizeof(COLORV)) {
                fputs("Code Error 1 in multi_process()\n", stderr);
                return(0);
        }
        while (--coff > 0) {            /* parent births children */
                int     pid = fork();
                if (pid < 0) {
                        fputs("fork() call failed!\n", stderr);
                        return(0);
                }
                if (pid == 0) break;    /* child gets to work */
        }
        osum = (float *)calloc((size_t)xres*yres, sizeof(float)*ncomp);
        if (!osum) {
                fprintf(stderr, "Cannot allocate %dx%d %d-component accumulator\n",
                                xres, yres, ncomp);
                return(0);
        }
        srandom(113*coff + 5669);       /* randomize row access for this process */
        syarr = scramble(yres);
                                        /* run through our unique set of columns */
        for (c = coff; c < cmtx->ncols; c += nprocs) {
                FILE    *fout;
                int     y;
                int     rc = cmtx->nrows;
                if (c > coff)           /* clear accumulator? */
                        memset(osum, 0, sizeof(float)*ncomp*xres*yres);
                while (rc-- > 0) {      /* map & sum each input file */
                        const int       r = odd ? rc : cmtx->nrows-1 - rc;
                        const rmx_dtype *cval = rmx_val(cmtx, r, c);
                        long            dstart;
                        size_t          maplen;
                        void            *imap;
                        FILE            *finp;
                        float           *dst;
                        int             i, x;
                        for (i = ncomp; i--; )
                                if (cval[i] != 0) break;
                        if (i < 0)      /* this coefficient is zero, skip */
                                continue;
                        sprintf(fbuf, in_spec, r);
                        finp = open_input(fbuf);
                        if (!finp)
                                return(0);
                        dstart = ftell(finp);
                        if (dstart < 0) {
                                fprintf(stderr, "%s: ftell() failed!\n", fbuf);
                                return(0);
                        }
                        if (in_type == DTfloat && dstart%sizeof(float)) {
                                fprintf(stderr, "%s: float header misalignment\n", fbuf);
                                return(0);
                        }
                        i = in_type==DTfloat ? ncomp*(int)sizeof(float) : ncomp+1;
                        maplen = dstart + (size_t)yres*xres*i;
                        imap = mmap(NULL, maplen, PROT_READ,
                                        MAP_FILE|MAP_SHARED, fileno(finp), 0);
                        fclose(finp);           /* will read from map (randomly) */
                        if (imap == MAP_FAILED) {
                                fprintf(stderr, "%s: unable to map input file\n", fbuf);
                                return(0);
                        }
                        if (in_type == DTfloat)
                            for (y = yres; y-- > 0; ) {
                                const float     *fvp = (float *)((char *)imap + dstart) +
                                                        (size_t)ncomp*xres*syarr[y];
                                dst = osum + (size_t)ncomp*xres*syarr[y];
                                for (x = xres; x-- > 0; dst += ncomp, fvp += ncomp)
                                    for (i = ncomp; i--; )
                                        dst[i] += cval[i]*fvp[i];
                            }
                        else
                            for (y = yres; y-- > 0; ) {
                                const COLRV     *cvp = (COLRV *)((char *)imap + dstart) +
                                                        (ncomp+1L)*xres*syarr[y];
                                dst = osum + (size_t)ncomp*xres*syarr[y];
                                for (x = xres; x-- > 0; dst += ncomp, cvp += ncomp+1) {
                                    const rmx_dtype     fe = cxponent[cvp[ncomp]];
                                    for (i = ncomp; i--; )
                                        dst[i] += cval[i]*(cvp[i]+(rmx_dtype).5)*fe;
                                }
                            }
                        munmap(imap, maplen);
                }                       /* write accumulated column picture/matrix */
                sprintf(fbuf, out_spec, c);
                fout = open_output(fbuf, c);
                if (!fout)
                        return(0);      /* assume error was reported */
                if (out_type != DTfloat) {
                        for (y = 0; y < yres; y++)
                                if (fwritesscan(osum + (size_t)y*xres*ncomp,
                                                        ncomp, xres, fout) < 0)
                                        goto writerr;
                } else if (fwrite(osum, sizeof(float)*ncomp, (size_t)xres*yres, fout) !=
                                        (size_t)xres*yres)
                        goto writerr;

                if (fbuf[0] == '!') {
                        if (pclose(fout) != 0) {
                                fprintf(stderr, "Bad status from: %s\n", fbuf);
                                return(0);
                        }
                } else if (fclose(fout) == EOF)
                        goto writerr;
                odd = !odd;             /* go back & forth to milk page cache */
        }
        free(osum);
        free(syarr);
        if (coff)                       /* child processes return here... */
                return(1);
        c = 0;                          /* ...but parent waits for children */
        while (++coff < nprocs) {
                int     st;
                if (wait(&st) < 0)
                        break;
                if (st) c = st;
        }
        return(c == 0);
writerr:
        fprintf(stderr, "%s: write error\n", fbuf);
        return(0);
}

#endif          /* ! Windows */

int
main(int argc, char *argv[])
{
        int     a;

        gargc = argc;                   /* for header output */
        gargv = argv;

        for (a = 1; a < argc-1 && argv[a][0] == '-'; a++)
                switch (argv[a][1]) {
                case 'o':               /* output spec/format */
                        switch (argv[a][2]) {
                        case '\0':
                                out_spec = argv[++a];
                                break;
                        case 'f':
                                out_type = DTfloat;
                                break;
                        case 'c':
                                out_type = DTrgbe;
                                break;
                        default:
                                goto badopt;
                        }
                        break;
                case 'N':               /* number of desired processes */
                case 'n':               /* quietly supported alternate */
                        nprocs = atoi(argv[++a]);
                        if (nprocs <= 0)
                                goto userr;
                        break;
                default:;
badopt:                 fprintf(stderr, "%s: bad option: %s\n", argv[0], argv[a]);
                        goto userr;
                        return(1);
                }
        if ((argc-a < 1) | (argc-a > 2) || argv[a][0] == '-')
                goto userr;
        in_spec = argv[a];
        cmtx = rmx_load(argv[a+1]);     /* loads from stdin if a+1==argc */
        if (cmtx == NULL)
                return(1);              /* error reported */
        if (nprocs > cmtx->ncols)
                nprocs = cmtx->ncols;
#if defined(_WIN32) || defined(_WIN64)
        if (nprocs > 1) {
                fprintf(stderr, "%s: warning - Windows only allows -N 1\n", argv[0]);
                nprocs = 1;
        }
#else
        if ((nprocs > 1) & !out_spec) {
                fprintf(stderr, "%s: multi-processing result cannot go to stdout\n",
                                argv[0]);
                nprocs = 1;
        }
        if ((nprocs > 1)  & iswapped && (in_type==DTfloat) | (in_type==DTdouble)) {
                fprintf(stderr, "%s: multi-processing unsupported on swapped input\n",
                                argv[0]);
                nprocs = 1;
        }
#endif
        if (cmtx->nrows > 1 && !hasFormat(in_spec)) {
                fprintf(stderr, "%s: input specification '%s' needs %%d format\n",
                                argv[0], in_spec);
                goto userr;
        }
        if (!get_iotypes())
                return(1);
        if (!(nprocs == 1 ? solo_process() : multi_process()))
                return(1);
        return(0);
userr:
        fprintf(stderr, "Usage: %s [-oc | -of][-o ospec][-N nproc] inpspec [mtx]\n",
                                argv[0]);
        return(1);
}
Revision:	2.8
Committed:	Fri Oct 24 22:41:10 2025 UTC (25 hours, 32 minutes ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 2.7:	+8 -1 lines
Log Message:	feat(pvsum): Add command arguments to output header(s)
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.8	static const char RCSid[] = "$Id: pvsum.c,v 2.7 2025/10/24 16:31:18 greg Exp $";
3	greg	2.1	#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	greg	2.2	#include "random.h"
14	greg	2.1	#include "rmatrix.h"
15			#if !defined(_WIN32) && !defined(_WIN64)
16			#include <sys/mman.h>
17	greg	2.5	#include <sys/wait.h>
18	greg	2.1	#endif
19
20	greg	2.7	#define VIEWSTR "VIEW=" /* borrowed from common/view.h */
21			#define VIEWSTRL 5
22
23	greg	2.1	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	greg	2.7	char viewspec[128] = ""; /* VIEW= line from first header */
33			char pixasp[48] = ""; /* PIXASPECT= line from header */
34	greg	2.1
35	greg	2.8	int gargc; /* global argc */
36			char *gargv; / global argv */
37
38	greg	2.1	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	greg	2.7	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	greg	2.1	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 + 4Lxresyres) {
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	greg	2.3	static int exposWarned = 0;
170			int xyres = (int )p;
171			char fmt[MAXFMTLEN];
172	greg	2.1
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	greg	2.3	if (!exposWarned && fabs(1. - exposval(s)) > 0.04) {
192			fputs("Warning - ignoring EXPOSURE setting(s)\n",
193			stderr);
194			exposWarned++;
195			}
196	greg	2.1	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	greg	2.3	} else if (!(fp = fopen(ospec, "w"))) {
252	greg	2.1	fprintf(stderr, "%s: cannot open for writing\n", ospec);
253			return(NULL);
254			}
255	greg	2.3	SET_FILE_BINARY(fp);
256	greg	2.1	newheader("RADIANCE", fp);
257	greg	2.2	if (cmtx->info) /* prepend matrix metadata */
258			fputs(cmtx->info, fp);
259			else
260			fputnow(fp);
261	greg	2.8	printargs(gargc, gargv, fp); /* this command */
262	greg	2.1	if (fno >= 0)
263			fprintf(fp, "FRAME=%d\n", fno);
264	greg	2.7	if (viewspec[0])
265			fputs(viewspec, fp);
266			if (pixasp[0])
267			fputs(pixasp, fp);
268	greg	2.1	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)xresyres, sizeof(float)ncomp);
311			COLORV iscan = (COLORV )malloc(sizeof(COLORV)ncompxres);
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)ncompxres*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 + yxres*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, ncompxres);
351			for (x = 0; x < xres; x++, dst += ncomp)
352			for (i = ncomp; i--; )
353			dst[i] += cval[i]iscan[xncomp + 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)yxresncomp,
378			ncomp, xres, fout) < 0)
379			goto writerr;
380			} else if (fwrite(osum, sizeof(float)ncomp, (size_t)xresyres, 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	greg	2.7	#if defined(_WIN32) \|\| defined(_WIN64)
404			#define multi_process solo_process
405			#else
406
407	greg	2.2	/* 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	greg	2.1	/* 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	greg	2.2	int *syarr;
439	greg	2.1	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)xresyres, 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	greg	2.2	srandom(113coff + 5669); / randomize row access for this process */
460			syarr = scramble(yres);
461	greg	2.1	/* 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)ncompxres*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	greg	2.2	long dstart;
472	greg	2.1	size_t maplen;
473			void *imap;
474			FILE *finp;
475			float *dst;
476	greg	2.2	int i, x;
477	greg	2.1	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	greg	2.6	maplen = dstart + (size_t)yresxresi;
496	greg	2.1	imap = mmap(NULL, maplen, PROT_READ,
497			MAP_FILE\|MAP_SHARED, fileno(finp), 0);
498	greg	2.3	fclose(finp); /* will read from map (randomly) */
499	greg	2.1	if (imap == MAP_FAILED) {
500			fprintf(stderr, "%s: unable to map input file\n", fbuf);
501			return(0);
502			}
503	greg	2.2	if (in_type == DTfloat)
504			for (y = yres; y-- > 0; ) {
505			const float fvp = (float )((char *)imap + dstart) +
506			(size_t)ncompxressyarr[y];
507			dst = osum + (size_t)ncompxressyarr[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)xressyarr[y];
516			dst = osum + (size_t)ncompxressyarr[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	greg	2.1	}
523			munmap(imap, maplen);
524	greg	2.3	} /* write accumulated column picture/matrix */
525	greg	2.1	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)yxresncomp,
532			ncomp, xres, fout) < 0)
533			goto writerr;
534			} else if (fwrite(osum, sizeof(float)ncomp, (size_t)xresyres, 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	greg	2.2	free(syarr);
549	greg	2.3	if (coff) /* child processes return here... */
550	greg	2.1	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	greg	2.7	#endif /* ! Windows */
565
566	greg	2.1	int
567			main(int argc, char *argv[])
568			{
569			int a;
570
571	greg	2.8	gargc = argc; /* for header output */
572			gargv = argv;
573
574	greg	2.1	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	greg	2.4	cmtx = rmx_load(argv[a+1]); /* loads from stdin if a+1==argc */
606	greg	2.1	if (cmtx == NULL)
607			return(1); /* error reported */
608	greg	2.3	if (nprocs > cmtx->ncols)
609			nprocs = cmtx->ncols;
610	greg	2.1	#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			}