src/util/rmtxop.c

#ifndef lint
static const char RCSid[] = "$Id: rmtxop.c,v 2.41 2025/04/18 23:59:03 greg Exp $";
#endif
/*
 * General component matrix operations.
 */

#include <errno.h>
#include "rtio.h"
#include "rmatrix.h"
#include "platform.h"

/* Preferred BSDF component:
        none, transmission, reflection front (normal side), reflection back */
typedef enum {RMPnone=-1, RMPtrans=0, RMPreflF, RMPreflB} RMPref;

/* Unary matrix operation(s) */
typedef struct {
        double          cmat[MAXCOMP*MAXCOMP];  /* component transformation */
        double          sca[MAXCOMP];           /* scalar coefficients */
        const char      *csym;                  /* symbolic coefs or file */
        short           clen;                   /* number of coefficients */
        short           nsf;                    /* number of scalars */
        short           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 */
int
loadmatrix(ROPMAT *rop)
{
        if (rop->mtx)                   /* already loaded? */
                return(0);
                                        /* check for BSDF input */
        if ((rop->inspec[0] != '!') & (rop->rmp != RMPnone)) {
                const char      *sp = strrchr(rop->inspec, '.');
                if (sp > rop->inspec && !strcasecmp(sp+1, "XML")) {
                        CMATRIX *cm = rop->rmp==RMPtrans ? cm_loadBTDF(rop->inspec) :
                                        cm_loadBRDF(rop->inspec, rop->rmp==RMPreflB) ;
                        if (!cm)
                                return(-1);
                        rop->mtx = rmx_from_cmatrix(cm);
                        cm_free(cm);
                        if (!rop->mtx)
                                return(-1);
                        rop->mtx->dtype = DTascii;
                        return(1);      /* loaded BSDF XML file */
                }
        }                               /* else load regular matrix */
        rop->mtx = rmx_load(rop->inspec);

        return(rop->mtx ? 1 : -1);
}

extern int      checksymbolic(ROPMAT *rop);

/* Check/set transform based on a reference input file */
int
checkreffile(ROPMAT *rop)
{
        static const char       *curRF = NULL;
        static RMATRIX          refm;
        const int               nc = rop->mtx->ncomp;
        int                     i;

        if (!curRF || strcmp(rop->preop.csym, curRF)) {
                FILE    *fp = fopen(rop->preop.csym, "rb");
                if (!rmx_load_header(&refm, fp)) {
                        fprintf(stderr, "%s: cannot read info header\n",
                                        rop->preop.csym);
                        curRF = NULL;
                        if (fp) fclose(fp);
                        return(-1);
                }
                fclose(fp);
                curRF = rop->preop.csym;
        }
        if (refm.ncomp == 3) {
                rop->preop.csym = (refm.dtype == DTxyze) ? "XYZ" : "RGB";
                return(checksymbolic(rop));
        }
        if (refm.ncomp == 2) {
                fprintf(stderr, "%s: cannot convert to 2 components\n",
                                curRF);
                return(-1);
        }
        if (refm.ncomp == 1) {
                rop->preop.csym = "Y";          /* XXX big assumption */
                return(checksymbolic(rop));
        }
        if (refm.ncomp == nc &&
                        !memcmp(refm.wlpart, rop->mtx->wlpart, sizeof(refm.wlpart)))
                return(0);                      /* nothing to do */

        if ((nc <= 3) | (nc > MAXCSAMP) | (refm.ncomp > MAXCSAMP)) {
                fprintf(stderr, "%s: cannot resample from %d to %d components\n",
                                curRF, nc, refm.ncomp);
                return(-1);
        }
        rop->preop.clen = refm.ncomp * nc;      /* compute spec to ref */

        for (i = 0; i < nc; i++) {
                SCOLOR  scstim, scresp;
                int     j;
                memset(scstim, 0, sizeof(COLORV)*nc);
                scstim[i] = 1.f;
                convertscolor(scresp, refm.ncomp, refm.wlpart[0], refm.wlpart[3],
                                scstim, nc, rop->mtx->wlpart[0], rop->mtx->wlpart[3]);
                for (j = refm.ncomp; j-- > 0; )
                        rop->preop.cmat[j*nc + i] = scresp[j];
        }
        memcpy(rop->mtx->wlpart, refm.wlpart, sizeof(rop->mtx->wlpart));
        return(0);
}

/* Compute conversion row from spectrum to one channel of RGB */
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 */
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 */
void
sensrow(ROPMAT *rop, int r, double (*sf)(const SCOLOR sc, int ncs, const float wlpt[4]))
{
        const int       nc = rop->mtx->ncomp;
        int             i;

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

/* Check/set symbolic transform */
int
checksymbolic(ROPMAT *rop)
{
        const int       nc = rop->mtx->ncomp;
        const int       dt = rop->mtx->dtype;
        double          cf = 1;
        int             i, j;
                                        /* check suffix => reference file */
        if (strchr(rop->preop.csym, '.') > rop->preop.csym)
                return(checkreffile(rop));

        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':
                case 'b':
                        ++comp;
                        /* fall through */
                case 'G':
                case 'g':
                        ++comp;
                        /* fall through */
                case 'R':
                case 'r':
                        if (rop->preop.csym[j] <= 'Z')
                                cf = 1./WHTEFFICACY;
                        if (dt == DTxyze) {
                                for (i = 3; i--; )
                                        rop->preop.cmat[j*nc+i] = cf*xyz2rgbmat[comp][i];
                        } else if (nc == 3)
                                rop->preop.cmat[j*nc+comp] = 1.;
                        else
                                rgbrow(rop, j, comp);
                        break;
                case 'Z':
                case 'z':
                        ++comp;
                        /* fall through */
                case 'Y':
                case 'y':
                        ++comp;
                        /* fall through */
                case 'X':
                case 'x':
                        if ((rop->preop.csym[j] <= 'Z') & (dt != DTxyze))
                                cf = WHTEFFICACY;
                        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*(cf != 1); i--; )
                                rop->preop.cmat[j*nc+i] *= cf;
                        break;
                case 'S':               /* scotopic (il)luminance */
                        cf = WHTSCOTOPIC;
                        /* fall through */
                case 's':
                        sensrow(rop, j, scolor2scotopic);
                        for (i = nc*(cf != 1); i--; )
                                rop->preop.cmat[j*nc+i] *= cf;
                        break;
                case 'M':               /* melanopic (il)luminance */
                        cf = WHTMELANOPIC;
                        /* fall through */
                case 'm':
                        sensrow(rop, j, scolor2melanopic);
                        for (i = nc*(cf != 1); i--; )
                                rop->preop.cmat[j*nc+i] *= cf;
                        break;
                case 'A':               /* average component */
                case 'a':
                        for (i = nc; i--; )
                                rop->preop.cmat[j*nc+i] = 1./(double)nc;
                        break;
                default:
                        fprintf(stderr, "%s: -c '%c' unsupported\n",
                                rop->inspec, rop->preop.csym[j]);
                        return(-1);
                }
        }
                                        /* return recommended output type */
        if (!strcasecmp(rop->preop.csym, "XYZ")) {
                if (dt <= DTspec)
                        return(DTxyze);
        } else if (!strcasecmp(rop->preop.csym, "RGB")) {
                if (dt <= DTspec)
                        return(DTrgbe);
        } else if (dt == DTspec)
                return(DTfloat);        /* probably not actual spectrum */
        return(0);
}

/* Get matrix and perform unary operations */
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 &&                  /* 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 != rop->preop.clen) {
                                fprintf(stderr, "%s: -s must have one or %d factors\n",
                                                rop->inspec,
                                                rop->preop.clen/rop->mtx->ncomp);
                                goto failure;
                        } else {
                                for (i = rop->preop.nsf; i--; )
                                        for (j = rop->mtx->ncomp; j--; )
                                                rop->preop.cmat[i*rop->mtx->ncomp+j]
                                                                *= 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;             /* now folded in */
                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? */
                if (!rmx_transpose(rop->mtx)) {
                        fputs(rop->inspec, stderr);
                        fputs(": transpose failed\n", stderr);
                        goto failure;
                }
                if (verbose) {
                        fputs(rop->inspec, stderr);
                        fputs(": transposed rows and columns\n", stderr);
                }
        }
        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 */
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 */
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 */
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? */
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));
}

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);
}

ROPMAT *
resize_moparr(ROPMAT *mop, int n2alloc)
{
        int     nmats = 0;
        int     i;

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