[ViewVC] Diff of: radiance/ray/src/util/rmtxop.c

Comparing ray/src/util/rmtxop.c (file contents):
Revision 2.12 by greg, Mon Aug 27 23:03:05 2018 UTC vs.
Revision 2.29 by greg, Sun Dec 3 03:44:42 2023 UTC

+ * General component matrix operations.
+ */
-–
+#include <stdio.h>
-–
+#include <stdlib.h>
+#include <errno.h>
+#include "rtio.h"
+#include "resolu.h"
+#include "rmatrix.h"
+#include "platform.h"
-<
+#define MAXCOMP         50              /* #components we support */
->
+#define MAXCOMP         MAXCSAMP        /* #components we support */
-+
+/* Unary matrix operation(s) */
+typedef struct {
-–
+        double          sca[MAXCOMP];           /* scalar coefficients */
-–
+        int             nsf;                    /* number of scalars */
+        double          cmat[MAXCOMP*MAXCOMP];  /* component transformation */
-<
+        int             clen;                   /* number of coefficients */
-<
+        int             transpose;              /* do transpose? */
-<
+        int             op;                     /* '*' or '+' */
-<
+} ROPERAT;                              /* matrix operation */
->
+        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 */
-+
+static int
-+
+loadmatrix(ROPMAT *rop)
-+
+{
-+
+        if (rop->mtx != NULL)           /* already loaded? */
-+
+                return(0);
-+
-+
+        rop->mtx = rmx_load(rop->inspec, rop->rmp);
-+
-+
+        return(!rop->mtx ? -1 : 1);
-+
+}
-+
-+
+static int      checksymbolic(ROPMAT *rop);
-+
-+
+/* Check/set transform based on a reference input file */
-+
+static 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) & (refm.dtype != DTspec)) {
-+
+                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 */
+static void
-<
+op_default(ROPERAT *op)
->
+rgbrow(ROPMAT *rop, int r, int p)
-<
+        memset(op, 0, sizeof(ROPERAT));
-<
+        op->op = '.';
->
+        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;
-+
+                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 */
-+
+static int
-+
+checksymbolic(ROPMAT *rop)
-+
+{
-+
+        const int       nc = rop->mtx->ncomp;
-+
+        const int       dt = rop->mtx->dtype;
-+
+        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':
-+
+                        ++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':               /* scotopic (il)luminance */
-+
+                        sensrow(rop, j, scolor2scotopic);
-+
+                        for (i = nc; i--; )
-+
+                                rop->preop.cmat[j*nc+i] *= WHTSCOTOPIC;
-+
+                        break;
-+
+                case 'M':               /* melanopic (il)luminance */
-+
+                        sensrow(rop, j, scolor2melanopic);
-+
+                        for (i = nc; i--; )
-+
+                                rop->preop.cmat[j*nc+i] *= WHTMELANOPIC;
-+
+                        break;
-+
+                case 'A':               /* average component */
-+
+                        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 (!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 *
-<
+operate(RMATRIX *mleft, ROPERAT *op, const char *fname)
->
+loadop(ROPMAT *rop)
-<
+        RMATRIX *mright = rmx_load(fname);
-<
+        RMATRIX *mtmp;
-<
+        int     i;
->
+        int     outtype = 0;
->
+        RMATRIX *mres;
->
+        int     i, j;
-<
+        if (fname == NULL)
-<
+                fname = "<stdin>";
-<
+        if (mright == NULL) {
-<
+                fputs(fname, stderr);
-<
+                fputs(": cannot load matrix\n", stderr);
->
+        if (loadmatrix(rop) < 0)                /* make sure we're loaded */
+                return(NULL);
-<
+        }
-<
+        if (op->nsf > 0) {              /* apply scalar(s) */
-<
+                if (op->clen > 0) {
-<
+                        fputs("Options -s and -c are exclusive\n", stderr);
-<
+                        rmx_free(mright);
-<
+                        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 (op->nsf == 1) {
-<
+                        for (i = mright->ncomp; --i; )
-<
+                                op->sca[i] = op->sca[0];
-<
+                } else if (op->nsf != mright->ncomp) {
->
+                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",
-<
+                                        fname, mright->ncomp);
-<
+                        rmx_free(mright);
-<
+                        return(NULL);
->
+                                        rop->inspec, rop->mtx->ncomp);
->
+                        goto failure;
-<
+                if ((mleft == NULL) | (op->op != '+') &&
-<
+                                !rmx_scale(mright, op->sca)) {
-<
+                        fputs(fname, stderr);
->
+                if (!rmx_scale(rop->mtx, rop->preop.sca)) {
->
+                        fputs(rop->inspec, stderr);
+                        fputs(": scalar operation failed\n", stderr);
-<
+                        rmx_free(mright);
-<
+                        return(NULL);
->
+                        goto failure;
+                if (verbose) {
-<
+                        fputs(fname, stderr);
->
+                        fputs(rop->inspec, stderr);
+                        fputs(": applied scalar (", stderr);
-<
+                        for (i = 0; i < op->nsf; i++)
-<
+                                fprintf(stderr, " %f", op->sca[i]);
->
+                        for (i = 0; i < rop->preop.nsf; i++)
->
+                                fprintf(stderr, " %f", rop->preop.sca[i]);
+                        fputs(" )\n", stderr);
-<
+        if (op->clen > 0) {             /* apply transform */
-<
+                if (op->clen % mright->ncomp) {
-<
+                        fprintf(stderr, "%s: -c must have N x %d coefficients\n",
-<
+                                        fname, mright->ncomp);
-<
+                        rmx_free(mright);
-<
+                        return(NULL);
-<
+                }
-<
+                mtmp = rmx_transform(mright, op->clen/mright->ncomp, op->cmat);
-<
+                if (mtmp == NULL) {
-<
+                        fprintf(stderr, "%s: matrix transform failed\n", fname);
-<
+                        rmx_free(mright);
-<
+                        return(NULL);
-<
+                }
-<
+                if (verbose)
-<
+                        fprintf(stderr, "%s: applied %d x %d transform\n",
-<
+                                        fname, mtmp->ncomp, mright->ncomp);
-<
+                rmx_free(mright);
-<
+                mright = mtmp;
-<
+        }
-<
+        if (op->transpose) {            /* transpose matrix? */
-<
+                mtmp = rmx_transpose(mright);
-<
+                if (mtmp == NULL) {
-<
+                        fputs(fname, stderr);
->
+        if (rop->preop.transpose) {             /* transpose matrix? */
->
+                mres = rmx_transpose(rop->mtx);
->
+                if (mres == NULL) {
->
+                        fputs(rop->inspec, stderr);
+                        fputs(": transpose failed\n", stderr);
-<
+                        rmx_free(mright);
-<
+                        return(NULL);
->
+                        goto failure;
+                if (verbose) {
-<
+                        fputs(fname, stderr);
->
+                        fputs(rop->inspec, stderr);
+                        fputs(": transposed rows and columns\n", stderr);
-<
+                rmx_free(mright);
-<
+                mright = mtmp;
->
+                rmx_free(rop->mtx);
->
+                rop->mtx = mres;
-<
+        if (mleft == NULL)              /* just one matrix */
-<
+                return(mright);
-<
+        if (op->op == '.') {            /* concatenate */
-<
+                RMATRIX *mres = rmx_multiply(mleft, mright);
->
+        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(fname, stderr);
-<
+                        if (mleft->ncols != mright->nrows)
-<
+                                fputs(": mismatched dimensions for multiply\n",
-<
+                                                stderr);
-<
+                        else
-<
+                                fputs(": concatenation failed\n", stderr);
-<
+                        rmx_free(mright);
->
+                        fputs(inspec, stderr);
->
+                        fputs(": concatenation failed\n", stderr);
+                        return(NULL);
+                if (verbose) {
-<
+                        fputs(fname, stderr);
->
+                        fputs(inspec, stderr);
+                        fputs(": concatenated matrix\n", stderr);
-<
+                rmx_free(mright);
-<
+                rmx_free(mleft);
-<
+                mleft = mres;
-<
+        } else if (op->op == '+') {
-<
+                if (!rmx_sum(mleft, mright, op->nsf ? op->sca : (double *)NULL)) {
-<
+                        fputs(fname, 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(fname, stderr);
->
+                        fputs(inspec, stderr);
+                        fputs(": added in matrix\n", stderr);
+                rmx_free(mright);
-<
+        } else if ((op->op == '*') | (op->op == '/')) {
-<
+                const char *    tnam = (op->op == '/') ?
->
+                mres = mleft;
->
+                break;
->
+        case '*':
->
+        case '/': {
->
+                const char *    tnam = (op == '/') ?
+                                        "division" : "multiplication";
+                errno = 0;
-<
+                if (!rmx_elemult(mleft, mright, (op->op == '/'))) {
->
+                if (!rmx_elemult(mleft, mright, (op == '/'))) {
+                        fprintf(stderr, "%s: element-wise %s failed\n",
-<
+                                        fname, tnam);
->
+                                        inspec, tnam);
->
+                        rmx_free(mleft);
+                        rmx_free(mright);
+                        return(NULL);
+                if (errno)
+                        fprintf(stderr,
+                                "%s: warning - error during element-wise %s\n",
-<
+                                        fname, tnam);
->
+                                        inspec, tnam);
+                else if (verbose)
-<
+                        fprintf(stderr, "%s: element-wise %s\n", fname, tnam);
->
+                        fprintf(stderr, "%s: element-wise %s\n", inspec, tnam);
+                rmx_free(mright);
-<
+        } else {
-<
+                fprintf(stderr, "%s: unknown operation '%c'\n", fname, op->op);
->
+                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;
+        return(ac);
-+
+static 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;
-<
+        RMATRIX *mres = NULL;
-<
+        ROPERAT op;
-<
+        int     i;
-<
+                                        /* initialize */
-<
+        op_default(&op);
->
+        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++)
->
+        for (i = 1; i < argc; i++) {
+                if (argv[i][0] && !argv[i][1] &&
-<
+                                strchr("+*/", argv[i][0]) != NULL) {
-<
+                        op.op = argv[i][0];
-<
+                } else if (argv[i][0] != '-' || !argv[i][1]) {
-<
+                        char    *fname = NULL;  /* load matrix */
-<
+                        if (argv[i][0] != '-')
-<
+                                fname = argv[i];
-<
+                        mres = operate(mres, &op, fname);
-<
+                        if (mres == NULL) {
-<
+                                fprintf(stderr, "%s: operation failed on '%s'\n",
-<
+                                                argv[0], argv[i]);
-<
+                                return(0);
->
+                                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);
-<
+                        op_default(&op);        /* reset operator */
->
+                        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 = !verbose;
->
+                                verbose++;
+                                break;
+                        case 'f':
+                                switch (argv[i][2]) {
+                                break;
+                        case 't':
-<
+                                op.transpose = 1;
->
+                                mop[nmats].preop.transpose = 1;
+                                break;
+                        case 's':
+                                if (n > MAXCOMP) n = MAXCOMP;
-<
+                                op.nsf = get_factors(op.sca, n, argv+i+1);
-<
+                                i += op.nsf;
->
+                                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;
-<
+                                op.clen = get_factors(op.cmat, n, argv+i+1);
-<
+                                i += op.clen;
->
+                                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 (mres == NULL)               /* check that we got something */
->
+                if (nmats >= nall)
->
+                        mop = resize_moparr(mop, nall += 2);
->
+        }
->
+        if (mop[0].inspec == NULL)      /* nothing to do? */
+                goto userr;
-<
+                                        /* write result to stdout */
-<
+        if (outfmt == DTfromHeader)
-<
+                outfmt = mres->dtype;
-<
+        if (outfmt != DTascii)
-<
+                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]);
->
+        if (mop[nmats-1].binop) {
->
+                fprintf(stderr,
->
+                        "%s: missing matrix argument after '%c' operation\n",
->
+                                argv[0], mop[nmats-1].binop);
+                return(1);
-<
+        /* rmx_free(mres); mres = NULL; */
-<
+        return(0);
->
+                                        /* 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;
->
+        }
->
+        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 ..] m1 [+*/] .. > mres\n",
->
+        "Usage: %s [-v][-f{adfc}][-t][-s sf .. | -c ce ..][-rf|-rb] m1 [.+*/] .. > mres\n",
+                        argv[0]);
+        return(1);

Diff Legend

-–
+Removed lines
-+
+Added lines
-<
+Changed lines
->
+Changed lines

Comparing ray/src/util/rmtxop.c (file contents): Revision 2.12 by greg, Mon Aug 27 23:03:05 2018 UTC vs. Revision 2.29 by greg, Sun Dec 3 03:44:42 2023 UTC

Diff Legend

Comparing ray/src/util/rmtxop.c (file contents):
Revision 2.12 by greg, Mon Aug 27 23:03:05 2018 UTC vs.
Revision 2.29 by greg, Sun Dec 3 03:44:42 2023 UTC