src/util/rmtxop.c

#ifndef lint
static const char RCSid[] = "$Id: rmtxop.c,v 2.21 2023/11/21 01:30:20 greg Exp $";
#endif
/*
 * General component matrix operations.
 */

#include <stdlib.h>
#include <errno.h>
#include <ctype.h>
#include "rtio.h"
#include "resolu.h"
#include "rmatrix.h"
#include "platform.h"

#define MAXCOMP         MAXCSAMP        /* #components we support */

/* Unary matrix operation(s) */
typedef struct {
        double          sca[MAXCOMP];           /* scalar coefficients */
        double          cmat[MAXCOMP*MAXCOMP];  /* component transformation */
        short           nsf;                    /* number of scalars */
        short           clen;                   /* number of coefficients */
        char            csym[11];               /* symbolic coefficients */
        char            transpose;              /* do transpose? */
} RUNARYOP;

/* Matrix input source and requested operation(s) */
typedef struct {
        const char      *inspec;                /* input specification */
        RMPref          rmp;                    /* matrix preference */
        RUNARYOP        preop;                  /* unary operation(s) */
        RMATRIX         *mtx;                   /* original matrix if loaded */
        int             binop;                  /* binary op with next (or 0) */
} ROPMAT;

int     verbose = 0;                    /* verbose reporting? */

/* Load matrix */
static int
loadmatrix(ROPMAT *rop)
{
        if (rop->mtx != NULL)           /* already loaded? */
                return(0);

        rop->mtx = rmx_load(rop->inspec, rop->rmp);
        if (rop->mtx == NULL) {
                fputs(rop->inspec, stderr);
                fputs(": cannot load matrix\n", stderr);
                return(-1);
        }
        return(1);
}

/* Compute conversion row from spectrum to one channel of RGB */
static void
rgbrow(ROPMAT *rop, int r, int p)
{
        const int       nc = rop->mtx->ncomp;
        const float *   wlp = rop->mtx->wlpart;
        int             i;

        for (i = nc; i--; ) {
                int     nmEnd = wlp[0] + (wlp[3] - wlp[0])*i/nc;
                int     nmStart = wlp[0] + (wlp[3] - wlp[0])*(i+1)/nc;
                COLOR   crgb;
                spec_rgb(crgb, nmStart, nmEnd);
                rop->preop.cmat[r*nc+i] = crgb[p];
        }
}

/* Compute conversion row from spectrum to one channel of XYZ */
static void
xyzrow(ROPMAT *rop, int r, int p)
{
        const int       nc = rop->mtx->ncomp;
        const float *   wlp = rop->mtx->wlpart;
        int             i;

        for (i = nc; i--; ) {
                int     nmEnd = wlp[0] + (wlp[3] - wlp[0])*i/nc;
                int     nmStart = wlp[0] + (wlp[3] - wlp[0])*(i+1)/nc;
                COLOR   cxyz;
                spec_cie(cxyz, nmStart, nmEnd);
                rop->preop.cmat[r*nc+i] = cxyz[p];
        }
}

/* Use the spectral sensitivity function to compute matrix coefficients */
static void
sensrow(ROPMAT *rop, int r, double (*sf)(SCOLOR sc, int ncs, const float wlpt[4]))
{
        const int       nc = rop->mtx->ncomp;
        int             i;

        for (i = nc; i--; ) {
                SCOLOR  sclr;
                scolorblack(sclr);
                sclr[i] = 1;
                rop->preop.cmat[r*nc+i] = (*sf)(sclr, nc, rop->mtx->wlpart);
        }
}

/* Check/set symbolic transform */
static int
checksymbolic(ROPMAT *rop)
{
        const int       nc = rop->mtx->ncomp;
        const int       dt = rop->mtx->dtype;
        int             i, j;

        if (nc < 3) {
                fprintf(stderr, "%s: -c '%s' requires at least 3 components\n",
                                rop->inspec, rop->preop.csym);
                return(-1);
        }
        rop->preop.clen = strlen(rop->preop.csym) * nc;
        if (rop->preop.clen > MAXCOMP*MAXCOMP) {
                fprintf(stderr, "%s: -c '%s' results in too many components\n",
                                rop->inspec, rop->preop.csym);
                return(-1);
        }
        for (j = 0; rop->preop.csym[j]; j++) {
                int     comp = 0;
                switch (rop->preop.csym[j]) {
                case 'B':
                        ++comp;
                        /* fall through */
                case 'G':
                        ++comp;
                        /* fall through */
                case 'R':
                        if (dt == DTxyze) {
                                for (i = 3; i--; )
                                        rop->preop.cmat[j*nc+i] = 1./WHTEFFICACY *
                                                        xyz2rgbmat[comp][i];
                        } else if (nc == 3)
                                rop->preop.cmat[j*nc+comp] = 1.;
                        else
                                rgbrow(rop, j, comp);
                        break;
                case 'Z':
                        ++comp;
                        /* fall through */
                case 'Y':
                        ++comp;
                        /* fall through */
                case 'X':
                        if (dt == DTxyze) {
                                rop->preop.cmat[j*nc+comp] = 1.;
                        } else if (nc == 3) {
                                for (i = 3; i--; )
                                        rop->preop.cmat[j*nc+i] =
                                                        rgb2xyzmat[comp][i];
                        } else if (comp == CIEY)
                                sensrow(rop, j, scolor2photopic);
                        else
                                xyzrow(rop, j, comp);

                        for (i = nc*(dt != DTxyze); i--; )
                                rop->preop.cmat[j*nc+i] *= WHTEFFICACY;
                        break;
                case 'S':
                        sensrow(rop, j, scolor2scotopic);
                        for (i = nc; i--; )
                                rop->preop.cmat[j*nc+i] *= WHTSCOTOPIC;
                        break;
                case 'M':
                        sensrow(rop, j, scolor2melanopic);
                        for (i = nc; i--; )
                                rop->preop.cmat[j*nc+i] *= WHTMELANOPIC;
                        break;
                default:
                        fprintf(stderr, "%s: -c '%c' unsupported\n",
                                rop->inspec, rop->preop.csym[j]);
                        return(-1);
                }
        }
                                        /* return recommended output type */
        if (!strcmp(rop->preop.csym, "XYZ")) {
                if (dt <= DTspec)
                        return(DTxyze);
        } else if (!strcmp(rop->preop.csym, "RGB")) {
                if (dt <= DTspec)
                        return(DTrgbe);
        }
        if ((nc > 3) & (dt <= DTspec))
                return(DTfloat);        /* probably not actual spectrum */
        return(0);
}

/* Get matrix and perform unary operations */
static RMATRIX *
loadop(ROPMAT *rop)
{
        int     outtype = 0;
        RMATRIX *mres;
        int     i, j;

        if (loadmatrix(rop) < 0)                /* make sure we're loaded */
                return(NULL);

        if (rop->preop.csym[0] &&               /* symbolic transform? */
                        (outtype = checksymbolic(rop)) < 0)
                goto failure;
        if (rop->preop.clen > 0) {              /* apply component transform? */
                if (rop->preop.clen % rop->mtx->ncomp) {
                        fprintf(stderr, "%s: -c must have N x %d coefficients\n",
                                        rop->inspec, rop->mtx->ncomp);
                        goto failure;
                }
                if (rop->preop.nsf > 0) {       /* scale transform, first */
                        if (rop->preop.nsf == 1) {
                                for (i = rop->preop.clen; i--; )
                                        rop->preop.cmat[i] *= rop->preop.sca[0];
                        } else if (rop->preop.nsf != rop->mtx->ncomp) {
                                fprintf(stderr, "%s: -s must have one or %d factors\n",
                                                rop->inspec, rop->mtx->ncomp);
                                goto failure;
                        } else {
                                for (j = rop->preop.clen/rop->preop.nsf; j--; )
                                        for (i = rop->preop.nsf; i--; )
                                                rop->preop.cmat[j*rop->preop.nsf+i] *=
                                                                rop->preop.sca[i];
                        }
                }
                mres = rmx_transform(rop->mtx, rop->preop.clen/rop->mtx->ncomp,
                                        rop->preop.cmat);
                if (mres == NULL) {
                        fprintf(stderr, "%s: matrix transform failed\n",
                                                rop->inspec);
                        goto failure;
                }
                if (verbose)
                        fprintf(stderr, "%s: applied %d x %d transform%s\n",
                                        rop->inspec, mres->ncomp,
                                        rop->mtx->ncomp,
                                        rop->preop.nsf ? " (* scalar)" : "");
                rop->preop.nsf = 0;
                if ((mres->ncomp > 3) & (mres->dtype <= DTspec))
                        outtype = DTfloat;      /* probably not actual spectrum */
                rmx_free(rop->mtx);
                rop->mtx = mres;
        }
        if (rop->preop.nsf > 0) {               /* apply scalar(s)? */
                if (rop->preop.nsf == 1) {
                        for (i = rop->mtx->ncomp; --i; )
                                rop->preop.sca[i] = rop->preop.sca[0];
                } else if (rop->preop.nsf != rop->mtx->ncomp) {
                        fprintf(stderr, "%s: -s must have one or %d factors\n",
                                        rop->inspec, rop->mtx->ncomp);
                        goto failure;
                }
                if (!rmx_scale(rop->mtx, rop->preop.sca)) {
                        fputs(rop->inspec, stderr);
                        fputs(": scalar operation failed\n", stderr);
                        goto failure;
                }
                if (verbose) {
                        fputs(rop->inspec, stderr);
                        fputs(": applied scalar (", stderr);
                        for (i = 0; i < rop->preop.nsf; i++)
                                fprintf(stderr, " %f", rop->preop.sca[i]);
                        fputs(" )\n", stderr);
                }
        }
        if (rop->preop.transpose) {             /* transpose matrix? */
                mres = rmx_transpose(rop->mtx);
                if (mres == NULL) {
                        fputs(rop->inspec, stderr);
                        fputs(": transpose failed\n", stderr);
                        goto failure;
                }
                if (verbose) {
                        fputs(rop->inspec, stderr);
                        fputs(": transposed rows and columns\n", stderr);
                }
                rmx_free(rop->mtx);
                rop->mtx = mres;
        }
        mres = rop->mtx;
        rop->mtx = NULL;
        if (outtype)
                mres->dtype = outtype;
        return(mres);
failure:
        rmx_free(rop->mtx);
        return(rop->mtx = NULL);
}

/* Execute binary operation, free matrix arguments and return new result */
static RMATRIX *
binaryop(const char *inspec, RMATRIX *mleft, int op, RMATRIX *mright)
{
        RMATRIX *mres = NULL;
        int     i;

        if ((mleft == NULL) | (mright == NULL))
                return(NULL);
        switch (op) {
        case '.':                       /* concatenate */
                if (mleft->ncomp != mright->ncomp) {
                        fputs(inspec, stderr);
                        fputs(": # components do not match\n", stderr);
                } else if (mleft->ncols != mright->nrows) {
                        fputs(inspec, stderr);
                        fputs(": mismatched dimensions\n",
                                        stderr);
                } else
                        mres = rmx_multiply(mleft, mright);
                rmx_free(mleft);
                rmx_free(mright);
                if (mres == NULL) {
                        fputs(inspec, stderr);
                        fputs(": concatenation failed\n", stderr);
                        return(NULL);
                }
                if (verbose) {
                        fputs(inspec, stderr);
                        fputs(": concatenated matrix\n", stderr);
                }
                break;
        case '+':
                if (!rmx_sum(mleft, mright, NULL)) {
                        fputs(inspec, stderr);
                        fputs(": matrix sum failed\n", stderr);
                        rmx_free(mleft);
                        rmx_free(mright);
                        return(NULL);
                }
                if (verbose) {
                        fputs(inspec, stderr);
                        fputs(": added in matrix\n", stderr);
                }
                rmx_free(mright);
                mres = mleft;
                break;
        case '*':
        case '/': {
                const char *    tnam = (op == '/') ?
                                        "division" : "multiplication";
                errno = 0;
                if (!rmx_elemult(mleft, mright, (op == '/'))) {
                        fprintf(stderr, "%s: element-wise %s failed\n",
                                        inspec, tnam);
                        rmx_free(mleft);
                        rmx_free(mright);
                        return(NULL);
                }
                if (errno)
                        fprintf(stderr,
                                "%s: warning - error during element-wise %s\n",
                                        inspec, tnam);
                else if (verbose)
                        fprintf(stderr, "%s: element-wise %s\n", inspec, tnam);
                rmx_free(mright);
                mres = mleft;
                } break;
        default:
                fprintf(stderr, "%s: unknown operation '%c'\n", inspec, op);
                rmx_free(mleft);
                rmx_free(mright);
                return(NULL);
        }
        return(mres);
}

/* Perform matrix operations from left to right */
static RMATRIX *
op_left2right(ROPMAT *mop)
{
        RMATRIX *mleft = loadop(mop);

        while (mop->binop) {
                if (mleft == NULL)
                        break;
                mleft = binaryop(mop[1].inspec,
                                mleft, mop->binop, loadop(mop+1));
                mop++;
        }
        return(mleft);
}

/* Perform matrix operations from right to left */
static RMATRIX *
op_right2left(ROPMAT *mop)
{
        RMATRIX *mright;
        int     rpos = 0;
                                        /* find end of list */
        while (mop[rpos].binop)
                if (mop[rpos++].binop != '.') {
                        fputs(
                "Right-to-left evaluation only for matrix multiplication!\n",
                                        stderr);
                        return(NULL);
                }
        mright = loadop(mop+rpos);
        while (rpos-- > 0) {
                if (mright == NULL)
                        break;
                mright = binaryop(mop[rpos+1].inspec,
                                loadop(mop+rpos), mop[rpos].binop, mright);
        }
        return(mright);
}

#define t_nrows(mop)    ((mop)->preop.transpose ? (mop)->mtx->ncols \
                                                : (mop)->mtx->nrows)
#define t_ncols(mop)    ((mop)->preop.transpose ? (mop)->mtx->nrows \
                                                : (mop)->mtx->ncols)

/* Should we prefer concatenating from rightmost matrix towards left? */
static int
prefer_right2left(ROPMAT *mop)
{
        int     mri = 0;

        while (mop[mri].binop)          /* find rightmost matrix */
                if (mop[mri++].binop != '.')
                        return(0);      /* pre-empt reversal for other ops */

        if (mri <= 1)
                return(0);              /* won't matter */

        if (loadmatrix(mop+mri) < 0)    /* load rightmost cat */
                return(1);              /* fail will bail in a moment */

        if (t_ncols(mop+mri) == 1)
                return(1);              /* definitely better R->L */

        if (t_ncols(mop+mri) >= t_nrows(mop+mri))
                return(0);              /* ...probably worse */

        if (loadmatrix(mop) < 0)        /* load leftmost */
                return(0);              /* fail will bail in a moment */

        return(t_ncols(mop+mri) < t_nrows(mop));
}

static int
get_factors(double da[], int n, char *av[])
{
        int     ac;

        for (ac = 0; ac < n && isflt(av[ac]); ac++)
                da[ac] = atof(av[ac]);
        return(ac);
}

static ROPMAT *
grow_moparray(ROPMAT *mop, int n2alloc)
{
        int     nmats = 0;

        while (mop[nmats++].binop)
                ;
        mop = (ROPMAT *)realloc(mop, n2alloc*sizeof(ROPMAT));
        if (mop == NULL) {
                fputs("Out of memory in grow_moparray()\n", stderr);
                exit(1);
        }
        if (n2alloc > nmats)
                memset(mop+nmats, 0, (n2alloc-nmats)*sizeof(ROPMAT));
        return(mop);
}

/* Load one or more matrices and operate on them, sending results to stdout */
int
main(int argc, char *argv[])
{
        int     outfmt = DTfromHeader;
        int     nall = 2;
        ROPMAT  *mop = (ROPMAT *)calloc(nall, sizeof(ROPMAT));
        int     nmats = 0;
        RMATRIX *mres = NULL;
        int     stdin_used = 0;
        int     i;
                                        /* get options and arguments */
        for (i = 1; i < argc; i++) {
                if (argv[i][0] && !argv[i][1] &&
                                strchr(".+*/", argv[i][0]) != NULL) {
                        if (!nmats || mop[nmats-1].binop) {
                                fprintf(stderr,
                        "%s: missing matrix argument before '%c' operation\n",
                                                argv[0], argv[i][0]);
                                return(1);
                        }
                        mop[nmats-1].binop = argv[i][0];
                } else if (argv[i][0] != '-' || !argv[i][1]) {
                        if (argv[i][0] == '-') {
                                if (stdin_used++) {
                                        fprintf(stderr,
                        "%s: standard input used for more than one matrix\n",
                                                argv[0]);
                                        return(1);
                                }
                                mop[nmats].inspec = stdin_name;
                        } else
                                mop[nmats].inspec = argv[i];
                        if (nmats > 0 && !mop[nmats-1].binop)
                                mop[nmats-1].binop = '.';
                        nmats++;
                } else {
                        int     n = argc-1 - i;
                        switch (argv[i][1]) {   /* get option */
                        case 'v':
                                verbose++;
                                break;
                        case 'f':
                                switch (argv[i][2]) {
                                case 'd':
                                        outfmt = DTdouble;
                                        break;
                                case 'f':
                                        outfmt = DTfloat;
                                        break;
                                case 'a':
                                        outfmt = DTascii;
                                        break;
                                case 'c':
                                        outfmt = DTrgbe;
                                        break;
                                default:
                                        goto userr;
                                }
                                break;
                        case 't':
                                mop[nmats].preop.transpose = 1;
                                break;
                        case 's':
                                if (n > MAXCOMP) n = MAXCOMP;
                                i += mop[nmats].preop.nsf =
                                        get_factors(mop[nmats].preop.sca,
                                                        n, argv+i+1);
                                if (mop[nmats].preop.nsf <= 0) {
                                        fprintf(stderr, "%s: -s missing arguments\n",
                                                        argv[0]);
                                        goto userr;
                                }
                                break;
                        case 'c':
                                if (n && isupper(argv[i+1][0])) {
                                        strlcpy(mop[nmats].preop.csym,
                                                argv[++i],
                                                sizeof(mop[0].preop.csym));
                                        mop[nmats].preop.clen = 0;
                                        break;
                                }
                                if (n > MAXCOMP*MAXCOMP) n = MAXCOMP*MAXCOMP;
                                i += mop[nmats].preop.clen =
                                        get_factors(mop[nmats].preop.cmat,
                                                        n, argv+i+1);
                                if (mop[nmats].preop.clen <= 0) {
                                        fprintf(stderr, "%s: -c missing arguments\n",
                                                        argv[0]);
                                        goto userr;
                                }
                                mop[nmats].preop.csym[0] = '\0';
                                break;
                        case 'r':
                                if (argv[i][2] == 'f')
                                        mop[nmats].rmp = RMPreflF;
                                else if (argv[i][2] == 'b')
                                        mop[nmats].rmp = RMPreflB;
                                else
                                        goto userr;
                                break;
                        default:
                                fprintf(stderr, "%s: unknown operation '%s'\n",
                                                argv[0], argv[i]);
                                goto userr;
                        }
                }
                if (nmats >= nall)
                        mop = grow_moparray(mop, nall += 2);
        }
        if (mop[0].inspec == NULL)      /* nothing to do? */
                goto userr;
        if (mop[nmats-1].binop) {
                fprintf(stderr,
                        "%s: missing matrix argument after '%c' operation\n",
                                argv[0], mop[nmats-1].binop);
                return(1);
        }
                                        /* favor quicker concatenation */
        mop[nmats].mtx = prefer_right2left(mop) ? op_right2left(mop)
                                                : op_left2right(mop);
        if (mop[nmats].mtx == NULL)
                return(1);
                                        /* apply trailing unary operations */
        mop[nmats].inspec = "trailing_ops";
        mres = loadop(mop+nmats);
        if (mres == NULL)
                return(1);
        if (outfmt == DTfromHeader)     /* check data type */
                outfmt = mres->dtype;
        if (outfmt == DTrgbe) {
                if (mres->ncomp > 3)
                        outfmt = DTspec;
                else if (mres->dtype == DTxyze)
                        outfmt = DTxyze;
        }
        if (outfmt != DTascii)          /* write result to stdout */
                SET_FILE_BINARY(stdout);
        newheader("RADIANCE", stdout);
        printargs(argc, argv, stdout);
        if (!rmx_write(mres, outfmt, stdout)) {
                fprintf(stderr, "%s: error writing result matrix\n", argv[0]);
                return(1);
        }
        /* rmx_free(mres); free(mop); */
        return(0);
userr:
        fprintf(stderr,
        "Usage: %s [-v][-f[adfc][-t][-s sf .. | -c ce ..][-rf|-rb] m1 [.+*/] .. > mres\n",
                        argv[0]);
        return(1);
}
Revision:	2.22
Committed:	Mon Nov 27 22:04:45 2023 UTC (23 months, 1 week ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.21:	+216 -35 lines
Log Message:	feat(rmtxop): Added symbolic transforms to -c option and made -s apply after any such transform, so both may be used simultaneously
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.22	static const char RCSid[] = "$Id: rmtxop.c,v 2.21 2023/11/21 01:30:20 greg Exp $";
3	greg	2.1	#endif
4			/*
5			* General component matrix operations.
6			*/
7
8			#include <stdlib.h>
9	greg	2.11	#include <errno.h>
10	greg	2.22	#include <ctype.h>
11	greg	2.1	#include "rtio.h"
12			#include "resolu.h"
13			#include "rmatrix.h"
14	greg	2.10	#include "platform.h"
15	greg	2.1
16	greg	2.22	#define MAXCOMP MAXCSAMP /* #components we support */
17	greg	2.1
18	greg	2.22	/* Unary matrix operation(s) */
19	greg	2.1	typedef struct {
20			double sca[MAXCOMP]; /* scalar coefficients */
21			double cmat[MAXCOMPMAXCOMP]; / component transformation */
22	greg	2.13	short nsf; /* number of scalars */
23			short clen; /* number of coefficients */
24	greg	2.22	char csym[11]; /* symbolic coefficients */
25			char transpose; /* do transpose? */
26	greg	2.13	} RUNARYOP;
27
28	greg	2.22	/* Matrix input source and requested operation(s) */
29	greg	2.13	typedef struct {
30			const char inspec; / input specification */
31	greg	2.18	RMPref rmp; /* matrix preference */
32	greg	2.13	RUNARYOP preop; /* unary operation(s) */
33			RMATRIX mtx; / original matrix if loaded */
34			int binop; /* binary op with next (or 0) */
35			} ROPMAT;
36	greg	2.1
37			int verbose = 0; /* verbose reporting? */
38
39	greg	2.13	/* Load matrix */
40			static int
41			loadmatrix(ROPMAT *rop)
42	greg	2.1	{
43	greg	2.20	if (rop->mtx != NULL) /* already loaded? */
44	greg	2.13	return(0);
45
46	greg	2.18	rop->mtx = rmx_load(rop->inspec, rop->rmp);
47	greg	2.13	if (rop->mtx == NULL) {
48			fputs(rop->inspec, stderr);
49			fputs(": cannot load matrix\n", stderr);
50			return(-1);
51			}
52			return(1);
53	greg	2.1	}
54
55	greg	2.22	/* Compute conversion row from spectrum to one channel of RGB */
56			static void
57			rgbrow(ROPMAT *rop, int r, int p)
58			{
59			const int nc = rop->mtx->ncomp;
60			const float * wlp = rop->mtx->wlpart;
61			int i;
62
63			for (i = nc; i--; ) {
64			int nmEnd = wlp[0] + (wlp[3] - wlp[0])*i/nc;
65			int nmStart = wlp[0] + (wlp[3] - wlp[0])*(i+1)/nc;
66			COLOR crgb;
67			spec_rgb(crgb, nmStart, nmEnd);
68			rop->preop.cmat[r*nc+i] = crgb[p];
69			}
70			}
71
72			/* Compute conversion row from spectrum to one channel of XYZ */
73			static void
74			xyzrow(ROPMAT *rop, int r, int p)
75			{
76			const int nc = rop->mtx->ncomp;
77			const float * wlp = rop->mtx->wlpart;
78			int i;
79
80			for (i = nc; i--; ) {
81			int nmEnd = wlp[0] + (wlp[3] - wlp[0])*i/nc;
82			int nmStart = wlp[0] + (wlp[3] - wlp[0])*(i+1)/nc;
83			COLOR cxyz;
84			spec_cie(cxyz, nmStart, nmEnd);
85			rop->preop.cmat[r*nc+i] = cxyz[p];
86			}
87			}
88
89			/* Use the spectral sensitivity function to compute matrix coefficients */
90			static void
91			sensrow(ROPMAT rop, int r, double (sf)(SCOLOR sc, int ncs, const float wlpt[4]))
92			{
93			const int nc = rop->mtx->ncomp;
94			int i;
95
96			for (i = nc; i--; ) {
97			SCOLOR sclr;
98			scolorblack(sclr);
99			sclr[i] = 1;
100			rop->preop.cmat[rnc+i] = (sf)(sclr, nc, rop->mtx->wlpart);
101			}
102			}
103
104			/* Check/set symbolic transform */
105			static int
106			checksymbolic(ROPMAT *rop)
107			{
108			const int nc = rop->mtx->ncomp;
109			const int dt = rop->mtx->dtype;
110			int i, j;
111
112			if (nc < 3) {
113			fprintf(stderr, "%s: -c '%s' requires at least 3 components\n",
114			rop->inspec, rop->preop.csym);
115			return(-1);
116			}
117			rop->preop.clen = strlen(rop->preop.csym) * nc;
118			if (rop->preop.clen > MAXCOMP*MAXCOMP) {
119			fprintf(stderr, "%s: -c '%s' results in too many components\n",
120			rop->inspec, rop->preop.csym);
121			return(-1);
122			}
123			for (j = 0; rop->preop.csym[j]; j++) {
124			int comp = 0;
125			switch (rop->preop.csym[j]) {
126			case 'B':
127			++comp;
128			/* fall through */
129			case 'G':
130			++comp;
131			/* fall through */
132			case 'R':
133			if (dt == DTxyze) {
134			for (i = 3; i--; )
135			rop->preop.cmat[jnc+i] = 1./WHTEFFICACY
136			xyz2rgbmat[comp][i];
137			} else if (nc == 3)
138			rop->preop.cmat[j*nc+comp] = 1.;
139			else
140			rgbrow(rop, j, comp);
141			break;
142			case 'Z':
143			++comp;
144			/* fall through */
145			case 'Y':
146			++comp;
147			/* fall through */
148			case 'X':
149			if (dt == DTxyze) {
150			rop->preop.cmat[j*nc+comp] = 1.;
151			} else if (nc == 3) {
152			for (i = 3; i--; )
153			rop->preop.cmat[j*nc+i] =
154			rgb2xyzmat[comp][i];
155			} else if (comp == CIEY)
156			sensrow(rop, j, scolor2photopic);
157			else
158			xyzrow(rop, j, comp);
159
160			for (i = nc*(dt != DTxyze); i--; )
161			rop->preop.cmat[jnc+i] = WHTEFFICACY;
162			break;
163			case 'S':
164			sensrow(rop, j, scolor2scotopic);
165			for (i = nc; i--; )
166			rop->preop.cmat[jnc+i] = WHTSCOTOPIC;
167			break;
168			case 'M':
169			sensrow(rop, j, scolor2melanopic);
170			for (i = nc; i--; )
171			rop->preop.cmat[jnc+i] = WHTMELANOPIC;
172			break;
173			default:
174			fprintf(stderr, "%s: -c '%c' unsupported\n",
175			rop->inspec, rop->preop.csym[j]);
176			return(-1);
177			}
178			}
179			/* return recommended output type */
180			if (!strcmp(rop->preop.csym, "XYZ")) {
181			if (dt <= DTspec)
182			return(DTxyze);
183			} else if (!strcmp(rop->preop.csym, "RGB")) {
184			if (dt <= DTspec)
185			return(DTrgbe);
186			}
187			if ((nc > 3) & (dt <= DTspec))
188			return(DTfloat); /* probably not actual spectrum */
189			return(0);
190			}
191
192	greg	2.13	/* Get matrix and perform unary operations */
193	greg	2.1	static RMATRIX *
194	greg	2.13	loadop(ROPMAT *rop)
195	greg	2.1	{
196	greg	2.22	int outtype = 0;
197	greg	2.13	RMATRIX *mres;
198	greg	2.22	int i, j;
199	greg	2.1
200	greg	2.13	if (loadmatrix(rop) < 0) /* make sure we're loaded */
201	greg	2.1	return(NULL);
202	greg	2.13
203	greg	2.22	if (rop->preop.csym[0] && /* symbolic transform? */
204			(outtype = checksymbolic(rop)) < 0)
205			goto failure;
206			if (rop->preop.clen > 0) { /* apply component transform? */
207			if (rop->preop.clen % rop->mtx->ncomp) {
208			fprintf(stderr, "%s: -c must have N x %d coefficients\n",
209			rop->inspec, rop->mtx->ncomp);
210			goto failure;
211			}
212			if (rop->preop.nsf > 0) { /* scale transform, first */
213			if (rop->preop.nsf == 1) {
214			for (i = rop->preop.clen; i--; )
215			rop->preop.cmat[i] *= rop->preop.sca[0];
216			} else if (rop->preop.nsf != rop->mtx->ncomp) {
217			fprintf(stderr, "%s: -s must have one or %d factors\n",
218			rop->inspec, rop->mtx->ncomp);
219			goto failure;
220			} else {
221			for (j = rop->preop.clen/rop->preop.nsf; j--; )
222			for (i = rop->preop.nsf; i--; )
223			rop->preop.cmat[jrop->preop.nsf+i] =
224			rop->preop.sca[i];
225			}
226			}
227			mres = rmx_transform(rop->mtx, rop->preop.clen/rop->mtx->ncomp,
228			rop->preop.cmat);
229			if (mres == NULL) {
230			fprintf(stderr, "%s: matrix transform failed\n",
231			rop->inspec);
232	greg	2.13	goto failure;
233	greg	2.1	}
234	greg	2.22	if (verbose)
235			fprintf(stderr, "%s: applied %d x %d transform%s\n",
236			rop->inspec, mres->ncomp,
237			rop->mtx->ncomp,
238			rop->preop.nsf ? " (* scalar)" : "");
239			rop->preop.nsf = 0;
240			if ((mres->ncomp > 3) & (mres->dtype <= DTspec))
241			outtype = DTfloat; /* probably not actual spectrum */
242			rmx_free(rop->mtx);
243			rop->mtx = mres;
244			}
245			if (rop->preop.nsf > 0) { /* apply scalar(s)? */
246	greg	2.13	if (rop->preop.nsf == 1) {
247			for (i = rop->mtx->ncomp; --i; )
248			rop->preop.sca[i] = rop->preop.sca[0];
249			} else if (rop->preop.nsf != rop->mtx->ncomp) {
250	greg	2.1	fprintf(stderr, "%s: -s must have one or %d factors\n",
251	greg	2.13	rop->inspec, rop->mtx->ncomp);
252			goto failure;
253	greg	2.1	}
254	greg	2.13	if (!rmx_scale(rop->mtx, rop->preop.sca)) {
255			fputs(rop->inspec, stderr);
256	greg	2.1	fputs(": scalar operation failed\n", stderr);
257	greg	2.13	goto failure;
258	greg	2.1	}
259			if (verbose) {
260	greg	2.13	fputs(rop->inspec, stderr);
261	greg	2.1	fputs(": applied scalar (", stderr);
262	greg	2.13	for (i = 0; i < rop->preop.nsf; i++)
263			fprintf(stderr, " %f", rop->preop.sca[i]);
264	greg	2.1	fputs(" )\n", stderr);
265			}
266			}
267	greg	2.22	if (rop->preop.transpose) { /* transpose matrix? */
268	greg	2.13	mres = rmx_transpose(rop->mtx);
269			if (mres == NULL) {
270			fputs(rop->inspec, stderr);
271	greg	2.12	fputs(": transpose failed\n", stderr);
272	greg	2.13	goto failure;
273	greg	2.12	}
274			if (verbose) {
275	greg	2.13	fputs(rop->inspec, stderr);
276	greg	2.12	fputs(": transposed rows and columns\n", stderr);
277			}
278	greg	2.13	rmx_free(rop->mtx);
279			rop->mtx = mres;
280			}
281			mres = rop->mtx;
282			rop->mtx = NULL;
283	greg	2.22	if (outtype)
284			mres->dtype = outtype;
285	greg	2.13	return(mres);
286			failure:
287			rmx_free(rop->mtx);
288			return(rop->mtx = NULL);
289			}
290
291			/* Execute binary operation, free matrix arguments and return new result */
292			static RMATRIX *
293			binaryop(const char inspec, RMATRIX mleft, int op, RMATRIX *mright)
294			{
295			RMATRIX *mres = NULL;
296			int i;
297
298			if ((mleft == NULL) \| (mright == NULL))
299			return(NULL);
300			switch (op) {
301			case '.': /* concatenate */
302	greg	2.16	if (mleft->ncomp != mright->ncomp) {
303			fputs(inspec, stderr);
304			fputs(": # components do not match\n", stderr);
305			} else if (mleft->ncols != mright->nrows) {
306			fputs(inspec, stderr);
307			fputs(": mismatched dimensions\n",
308			stderr);
309			} else
310			mres = rmx_multiply(mleft, mright);
311	greg	2.13	rmx_free(mleft);
312	greg	2.12	rmx_free(mright);
313	greg	2.1	if (mres == NULL) {
314	greg	2.13	fputs(inspec, stderr);
315	greg	2.16	fputs(": concatenation failed\n", stderr);
316	greg	2.1	return(NULL);
317			}
318			if (verbose) {
319	greg	2.13	fputs(inspec, stderr);
320	greg	2.1	fputs(": concatenated matrix\n", stderr);
321			}
322	greg	2.13	break;
323			case '+':
324			if (!rmx_sum(mleft, mright, NULL)) {
325			fputs(inspec, stderr);
326	greg	2.1	fputs(": matrix sum failed\n", stderr);
327	greg	2.13	rmx_free(mleft);
328	greg	2.1	rmx_free(mright);
329			return(NULL);
330			}
331			if (verbose) {
332	greg	2.13	fputs(inspec, stderr);
333	greg	2.1	fputs(": added in matrix\n", stderr);
334			}
335			rmx_free(mright);
336	greg	2.13	mres = mleft;
337			break;
338			case '*':
339			case '/': {
340			const char * tnam = (op == '/') ?
341	greg	2.11	"division" : "multiplication";
342			errno = 0;
343	greg	2.13	if (!rmx_elemult(mleft, mright, (op == '/'))) {
344	greg	2.11	fprintf(stderr, "%s: element-wise %s failed\n",
345	greg	2.13	inspec, tnam);
346			rmx_free(mleft);
347	greg	2.11	rmx_free(mright);
348			return(NULL);
349			}
350			if (errno)
351			fprintf(stderr,
352			"%s: warning - error during element-wise %s\n",
353	greg	2.13	inspec, tnam);
354	greg	2.11	else if (verbose)
355	greg	2.13	fprintf(stderr, "%s: element-wise %s\n", inspec, tnam);
356	greg	2.11	rmx_free(mright);
357	greg	2.13	mres = mleft;
358			} break;
359			default:
360			fprintf(stderr, "%s: unknown operation '%c'\n", inspec, op);
361			rmx_free(mleft);
362	greg	2.1	rmx_free(mright);
363			return(NULL);
364			}
365	greg	2.13	return(mres);
366			}
367
368			/* Perform matrix operations from left to right */
369			static RMATRIX *
370			op_left2right(ROPMAT *mop)
371			{
372			RMATRIX *mleft = loadop(mop);
373
374			while (mop->binop) {
375			if (mleft == NULL)
376			break;
377			mleft = binaryop(mop[1].inspec,
378			mleft, mop->binop, loadop(mop+1));
379			mop++;
380			}
381	greg	2.1	return(mleft);
382			}
383
384	greg	2.13	/* Perform matrix operations from right to left */
385			static RMATRIX *
386			op_right2left(ROPMAT *mop)
387			{
388			RMATRIX *mright;
389			int rpos = 0;
390			/* find end of list */
391			while (mop[rpos].binop)
392	greg	2.14	if (mop[rpos++].binop != '.') {
393			fputs(
394			"Right-to-left evaluation only for matrix multiplication!\n",
395			stderr);
396			return(NULL);
397			}
398	greg	2.13	mright = loadop(mop+rpos);
399			while (rpos-- > 0) {
400			if (mright == NULL)
401			break;
402	greg	2.20	mright = binaryop(mop[rpos+1].inspec,
403	greg	2.13	loadop(mop+rpos), mop[rpos].binop, mright);
404			}
405			return(mright);
406			}
407
408			#define t_nrows(mop) ((mop)->preop.transpose ? (mop)->mtx->ncols \
409			: (mop)->mtx->nrows)
410			#define t_ncols(mop) ((mop)->preop.transpose ? (mop)->mtx->nrows \
411			: (mop)->mtx->ncols)
412
413			/* Should we prefer concatenating from rightmost matrix towards left? */
414			static int
415			prefer_right2left(ROPMAT *mop)
416			{
417			int mri = 0;
418
419			while (mop[mri].binop) /* find rightmost matrix */
420			if (mop[mri++].binop != '.')
421			return(0); /* pre-empt reversal for other ops */
422
423			if (mri <= 1)
424			return(0); /* won't matter */
425
426			if (loadmatrix(mop+mri) < 0) /* load rightmost cat */
427			return(1); /* fail will bail in a moment */
428
429			if (t_ncols(mop+mri) == 1)
430			return(1); /* definitely better R->L */
431
432			if (t_ncols(mop+mri) >= t_nrows(mop+mri))
433			return(0); /* ...probably worse */
434
435			if (loadmatrix(mop) < 0) /* load leftmost */
436			return(0); /* fail will bail in a moment */
437
438			return(t_ncols(mop+mri) < t_nrows(mop));
439			}
440
441	greg	2.1	static int
442			get_factors(double da[], int n, char *av[])
443			{
444			int ac;
445
446			for (ac = 0; ac < n && isflt(av[ac]); ac++)
447			da[ac] = atof(av[ac]);
448			return(ac);
449			}
450
451	greg	2.13	static ROPMAT *
452			grow_moparray(ROPMAT *mop, int n2alloc)
453			{
454			int nmats = 0;
455
456			while (mop[nmats++].binop)
457			;
458			mop = (ROPMAT )realloc(mop, n2allocsizeof(ROPMAT));
459			if (mop == NULL) {
460			fputs("Out of memory in grow_moparray()\n", stderr);
461			exit(1);
462			}
463			if (n2alloc > nmats)
464			memset(mop+nmats, 0, (n2alloc-nmats)*sizeof(ROPMAT));
465			return(mop);
466			}
467
468	greg	2.1	/* Load one or more matrices and operate on them, sending results to stdout */
469			int
470			main(int argc, char *argv[])
471			{
472	greg	2.7	int outfmt = DTfromHeader;
473	greg	2.13	int nall = 2;
474			ROPMAT mop = (ROPMAT )calloc(nall, sizeof(ROPMAT));
475			int nmats = 0;
476	greg	2.1	RMATRIX *mres = NULL;
477	greg	2.13	int stdin_used = 0;
478	greg	2.1	int i;
479			/* get options and arguments */
480	greg	2.13	for (i = 1; i < argc; i++) {
481	greg	2.11	if (argv[i][0] && !argv[i][1] &&
482	greg	2.13	strchr(".+*/", argv[i][0]) != NULL) {
483	greg	2.15	if (!nmats \|\| mop[nmats-1].binop) {
484	greg	2.14	fprintf(stderr,
485	greg	2.15	"%s: missing matrix argument before '%c' operation\n",
486	greg	2.14	argv[0], argv[i][0]);
487	greg	2.13	return(1);
488			}
489	greg	2.15	mop[nmats-1].binop = argv[i][0];
490	greg	2.1	} else if (argv[i][0] != '-' \|\| !argv[i][1]) {
491	greg	2.13	if (argv[i][0] == '-') {
492			if (stdin_used++) {
493			fprintf(stderr,
494			"%s: standard input used for more than one matrix\n",
495			argv[0]);
496			return(1);
497			}
498			mop[nmats].inspec = stdin_name;
499			} else
500			mop[nmats].inspec = argv[i];
501			if (nmats > 0 && !mop[nmats-1].binop)
502			mop[nmats-1].binop = '.';
503			nmats++;
504	greg	2.1	} else {
505			int n = argc-1 - i;
506			switch (argv[i][1]) { /* get option */
507			case 'v':
508	greg	2.14	verbose++;
509	greg	2.1	break;
510			case 'f':
511			switch (argv[i][2]) {
512			case 'd':
513			outfmt = DTdouble;
514			break;
515			case 'f':
516			outfmt = DTfloat;
517			break;
518			case 'a':
519			outfmt = DTascii;
520			break;
521			case 'c':
522			outfmt = DTrgbe;
523			break;
524			default:
525			goto userr;
526			}
527			break;
528			case 't':
529	greg	2.13	mop[nmats].preop.transpose = 1;
530	greg	2.1	break;
531			case 's':
532			if (n > MAXCOMP) n = MAXCOMP;
533	greg	2.13	i += mop[nmats].preop.nsf =
534			get_factors(mop[nmats].preop.sca,
535			n, argv+i+1);
536	greg	2.22	if (mop[nmats].preop.nsf <= 0) {
537			fprintf(stderr, "%s: -s missing arguments\n",
538			argv[0]);
539			goto userr;
540			}
541	greg	2.1	break;
542			case 'c':
543	greg	2.22	if (n && isupper(argv[i+1][0])) {
544			strlcpy(mop[nmats].preop.csym,
545			argv[++i],
546			sizeof(mop[0].preop.csym));
547			mop[nmats].preop.clen = 0;
548			break;
549			}
550	greg	2.1	if (n > MAXCOMPMAXCOMP) n = MAXCOMPMAXCOMP;
551	greg	2.13	i += mop[nmats].preop.clen =
552			get_factors(mop[nmats].preop.cmat,
553			n, argv+i+1);
554	greg	2.22	if (mop[nmats].preop.clen <= 0) {
555			fprintf(stderr, "%s: -c missing arguments\n",
556			argv[0]);
557			goto userr;
558			}
559			mop[nmats].preop.csym[0] = '\0';
560	greg	2.1	break;
561	greg	2.18	case 'r':
562			if (argv[i][2] == 'f')
563			mop[nmats].rmp = RMPreflF;
564			else if (argv[i][2] == 'b')
565			mop[nmats].rmp = RMPreflB;
566			else
567			goto userr;
568			break;
569	greg	2.1	default:
570			fprintf(stderr, "%s: unknown operation '%s'\n",
571			argv[0], argv[i]);
572			goto userr;
573			}
574			}
575	greg	2.14	if (nmats >= nall)
576			mop = grow_moparray(mop, nall += 2);
577	greg	2.13	}
578			if (mop[0].inspec == NULL) /* nothing to do? */
579	greg	2.1	goto userr;
580	greg	2.15	if (mop[nmats-1].binop) {
581			fprintf(stderr,
582			"%s: missing matrix argument after '%c' operation\n",
583			argv[0], mop[nmats-1].binop);
584			return(1);
585			}
586	greg	2.13	/* favor quicker concatenation */
587	greg	2.14	mop[nmats].mtx = prefer_right2left(mop) ? op_right2left(mop)
588			: op_left2right(mop);
589			if (mop[nmats].mtx == NULL)
590			return(1);
591			/* apply trailing unary operations */
592			mop[nmats].inspec = "trailing_ops";
593			mres = loadop(mop+nmats);
594			if (mres == NULL)
595	greg	2.13	return(1);
596	greg	2.22	if (outfmt == DTfromHeader) /* check data type */
597	greg	2.6	outfmt = mres->dtype;
598	greg	2.22	if (outfmt == DTrgbe) {
599			if (mres->ncomp > 3)
600			outfmt = DTspec;
601			else if (mres->dtype == DTxyze)
602			outfmt = DTxyze;
603			}
604			if (outfmt != DTascii) /* write result to stdout */
605	greg	2.10	SET_FILE_BINARY(stdout);
606	greg	2.1	newheader("RADIANCE", stdout);
607			printargs(argc, argv, stdout);
608	greg	2.5	if (!rmx_write(mres, outfmt, stdout)) {
609	greg	2.1	fprintf(stderr, "%s: error writing result matrix\n", argv[0]);
610			return(1);
611			}
612	greg	2.13	/* rmx_free(mres); free(mop); */
613	greg	2.1	return(0);
614			userr:
615			fprintf(stderr,
616	greg	2.19	"Usage: %s [-v][-f[adfc][-t][-s sf .. \| -c ce ..][-rf\|-rb] m1 [.+*/] .. > mres\n",
617	greg	2.1	argv[0]);
618			return(1);
619			}