ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmtxop.c
(Generate patch)

Comparing ray/src/util/rmtxop.c (file contents):
Revision 2.12 by greg, Mon Aug 27 23:03:05 2018 UTC vs.
Revision 2.26 by greg, Fri Dec 1 02:05:00 2023 UTC

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

Diff Legend

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