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.13 by greg, Mon Aug 12 01:20:26 2019 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         16              /* #components we support */
15 > #define MAXCOMP         MAXCSAMP        /* #components we support */
16  
17 < static const char       stdin_name[] = "<stdin>";
19 <
20 < /* unary matrix operation(s) */
17 > /* Unary matrix operation(s) */
18   typedef struct {
19          double          sca[MAXCOMP];           /* scalar coefficients */
20          double          cmat[MAXCOMP*MAXCOMP];  /* component transformation */
21          short           nsf;                    /* number of scalars */
22          short           clen;                   /* number of coefficients */
23 <        short           transpose;              /* do transpose? */
23 >        char            csym[11];               /* symbolic coefficients */
24 >        char            transpose;              /* do transpose? */
25   } RUNARYOP;
26  
27 < /* matrix input source and requested operation(s) */
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) */
# Line 40 | Line 39 | int    verbose = 0;                    /* verbose reporting? */
39   static int
40   loadmatrix(ROPMAT *rop)
41   {
42 <        if (rop->mtx != NULL)
42 >        if (rop->mtx != NULL)           /* already loaded? */
43                  return(0);
44  
45 <        rop->mtx = rmx_load(rop->inspec == stdin_name ?
46 <                                (const char *)NULL : rop->inspec);
47 <        if (rop->mtx == NULL) {
48 <                fputs(rop->inspec, stderr);
49 <                fputs(": cannot load matrix\n", stderr);
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 > rgbrow(ROPMAT *rop, int r, int p)
53 > {
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 <        return(1);
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   loadop(ROPMAT *rop)
194   {
195 +        int     outtype = 0;
196          RMATRIX *mres;
197 <        int     i;
197 >        int     i, j;
198  
199          if (loadmatrix(rop) < 0)                /* make sure we're loaded */
200                  return(NULL);
201  
202 <        if (rop->preop.nsf > 0) {               /* apply scalar(s) */
203 <                if (rop->preop.clen > 0) {
204 <                        fputs("Options -s and -c are exclusive\n", stderr);
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 (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];
# Line 89 | Line 263 | loadop(ROPMAT *rop)
263                          fputs(" )\n", stderr);
264                  }
265          }
266 <        if (rop->preop.clen > 0) {      /* apply transform */
93 <                if (rop->preop.clen % rop->mtx->ncomp) {
94 <                        fprintf(stderr, "%s: -c must have N x %d coefficients\n",
95 <                                        rop->inspec, rop->mtx->ncomp);
96 <                        goto failure;
97 <                }
98 <                mres = rmx_transform(rop->mtx, rop->preop.clen/rop->mtx->ncomp,
99 <                                        rop->preop.cmat);
100 <                if (mres == NULL) {
101 <                        fprintf(stderr, "%s: matrix transform failed\n",
102 <                                                rop->inspec);
103 <                        goto failure;
104 <                }
105 <                if (verbose)
106 <                        fprintf(stderr, "%s: applied %d x %d transform\n",
107 <                                        rop->inspec, mres->ncomp,
108 <                                        rop->mtx->ncomp);
109 <                rmx_free(rop->mtx);
110 <                rop->mtx = mres;
111 <        }
112 <        if (rop->preop.transpose) {     /* transpose matrix? */
266 >        if (rop->preop.transpose) {             /* transpose matrix? */
267                  mres = rmx_transpose(rop->mtx);
268                  if (mres == NULL) {
269                          fputs(rop->inspec, stderr);
# Line 125 | Line 279 | loadop(ROPMAT *rop)
279          }
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);
# Line 140 | Line 296 | binaryop(const char *inspec, RMATRIX *mleft, int op, R
296  
297          if ((mleft == NULL) | (mright == NULL))
298                  return(NULL);
143 –
299          switch (op) {
300          case '.':                       /* concatenate */
301 <                mres = rmx_multiply(mleft, mright);
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(inspec, stderr);
314 <                        if (mleft->ncols != mright->nrows)
152 <                                fputs(": mismatched dimensions for multiply\n",
153 <                                                stderr);
154 <                        else
155 <                                fputs(": concatenation failed\n", stderr);
314 >                        fputs(": concatenation failed\n", stderr);
315                          return(NULL);
316                  }
317                  if (verbose) {
# Line 229 | Line 388 | op_right2left(ROPMAT *mop)
388          int     rpos = 0;
389                                          /* find end of list */
390          while (mop[rpos].binop)
391 <                rpos++;
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].inspec,
401 >                mright = binaryop(mop[rpos+1].inspec,
402                                  loadop(mop+rpos), mop[rpos].binop, mright);
403          }
404          return(mright);
# Line 284 | Line 448 | get_factors(double da[], int n, char *av[])
448   }
449  
450   static ROPMAT *
451 < grow_moparray(ROPMAT *mop, int n2alloc)
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 grow_moparray()\n", stderr);
462 >                fputs("Out of memory in resize_moparr()\n", stderr);
463                  exit(1);
464          }
465          if (n2alloc > nmats)
# Line 313 | Line 480 | main(int argc, char *argv[])
480          int     i;
481                                          /* get options and arguments */
482          for (i = 1; i < argc; i++) {
316 –                if (nmats >= nall)
317 –                        mop = grow_moparray(mop, nall += 2);
483                  if (argv[i][0] && !argv[i][1] &&
484                                  strchr(".+*/", argv[i][0]) != NULL) {
485 <                        if (mop[nmats].inspec == NULL || mop[nmats].binop) {
486 <                                fprintf(stderr, "%s: missing matrix argument\n",
487 <                                                argv[0]);
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 <                        mop[nmats++].binop = argv[i][0];
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++) {
# Line 341 | Line 507 | main(int argc, char *argv[])
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 369 | Line 535 | main(int argc, char *argv[])
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                                  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 (nmats >= nall)
578 +                        mop = resize_moparr(mop, nall += 2);
579          }
580          if (mop[0].inspec == NULL)      /* nothing to do? */
581                  goto userr;
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                                          /* favor quicker concatenation */
589 <        mres = prefer_right2left(mop) ? op_right2left(mop) : op_left2right(mop);
590 <        if (!mres)
589 >        mop[nmats].mtx = prefer_right2left(mop) ? op_right2left(mop)
590 >                                                : op_left2right(mop);
591 >        if (mop[nmats].mtx == NULL)
592                  return(1);
593 <                                        /* write result to stdout */
594 <        if (outfmt == DTfromHeader)
595 <                outfmt = mres->dtype;
596 <        if (outfmt != DTascii)
396 <                SET_FILE_BINARY(stdout);
397 <        newheader("RADIANCE", stdout);
398 <        printargs(argc, argv, stdout);
399 <        if (!rmx_write(mres, outfmt, stdout)) {
400 <                fprintf(stderr, "%s: error writing result matrix\n", argv[0]);
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 <        /* rmx_free(mres); free(mop); */
607 <        return(0);
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