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.23 by greg, Tue Nov 28 16:36:50 2023 UTC

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

Diff Legend

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