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.25 by greg, Wed Nov 29 18:57:10 2023 UTC vs.
Revision 2.32 by greg, Tue Dec 19 16:09:20 2023 UTC

# Line 5 | Line 5 | static const char RCSid[] = "$Id$";
5   * General component matrix operations.
6   */
7  
8 – #include <stdlib.h>
8   #include <errno.h>
10 – #include <ctype.h>
9   #include "rtio.h"
10   #include "resolu.h"
11   #include "rmatrix.h"
# Line 17 | Line 15 | static const char RCSid[] = "$Id$";
15  
16   /* Unary matrix operation(s) */
17   typedef struct {
20 –        double          sca[MAXCOMP];           /* scalar coefficients */
18          double          cmat[MAXCOMP*MAXCOMP];  /* component transformation */
19 <        short           nsf;                    /* number of scalars */
19 >        double          sca[MAXCOMP];           /* scalar coefficients */
20 >        const char      *csym;                  /* symbolic coefs or file */
21          short           clen;                   /* number of coefficients */
22 <        char            csym[11];               /* symbolic coefficients */
23 <        char            transpose;              /* do transpose? */
22 >        short           nsf;                    /* number of scalars */
23 >        short           transpose;              /* do transpose? */
24   } RUNARYOP;
25  
26   /* Matrix input source and requested operation(s) */
# Line 48 | Line 46 | loadmatrix(ROPMAT *rop)
46          return(!rop->mtx ? -1 : 1);
47   }
48  
49 + static int      checksymbolic(ROPMAT *rop);
50 +
51 + /* Check/set transform based on a reference input file */
52 + static int
53 + checkreffile(ROPMAT *rop)
54 + {
55 +        static const char       *curRF = NULL;
56 +        static RMATRIX          refm;
57 +        const int               nc = rop->mtx->ncomp;
58 +        int                     i;
59 +
60 +        if (!curRF || strcmp(rop->preop.csym, curRF)) {
61 +                FILE    *fp = fopen(rop->preop.csym, "rb");
62 +                if (!rmx_load_header(&refm, fp)) {
63 +                        fprintf(stderr, "%s: cannot read info header\n",
64 +                                        rop->preop.csym);
65 +                        curRF = NULL;
66 +                        if (fp) fclose(fp);
67 +                        return(-1);
68 +                }
69 +                fclose(fp);
70 +                curRF = rop->preop.csym;
71 +        }
72 +        if (refm.ncomp == 3) {
73 +                rop->preop.csym = (refm.dtype == DTxyze) ? "XYZ" : "RGB";
74 +                return(checksymbolic(rop));
75 +        }
76 +        if (refm.ncomp == 2) {
77 +                fprintf(stderr, "%s: cannot convert to 2 components\n",
78 +                                curRF);
79 +                return(-1);
80 +        }
81 +        if (refm.ncomp == 1) {
82 +                rop->preop.csym = "Y";          /* XXX big assumption */
83 +                return(checksymbolic(rop));
84 +        }
85 +        if (refm.ncomp == nc &&
86 +                        !memcmp(refm.wlpart, rop->mtx->wlpart, sizeof(refm.wlpart)))
87 +                return(0);                      /* nothing to do */
88 +
89 +        if ((nc <= 3) | (nc > MAXCSAMP) | (refm.ncomp > MAXCSAMP)) {
90 +                fprintf(stderr, "%s: cannot resample from %d to %d components\n",
91 +                                curRF, nc, refm.ncomp);
92 +                return(-1);
93 +        }
94 +        rop->preop.clen = refm.ncomp * nc;      /* compute spec to ref */
95 +
96 +        for (i = 0; i < nc; i++) {
97 +                SCOLOR  scstim, scresp;
98 +                int     j;
99 +                memset(scstim, 0, sizeof(COLORV)*nc);
100 +                scstim[i] = 1.f;
101 +                convertscolor(scresp, refm.ncomp, refm.wlpart[0], refm.wlpart[3],
102 +                                scstim, nc, rop->mtx->wlpart[0], rop->mtx->wlpart[3]);
103 +                for (j = refm.ncomp; j-- > 0; )
104 +                        rop->preop.cmat[j*nc + i] = scresp[j];
105 +        }
106 +        memcpy(rop->mtx->wlpart, refm.wlpart, sizeof(rop->mtx->wlpart));
107 +        return(0);
108 + }
109 +
110   /* Compute conversion row from spectrum to one channel of RGB */
111   static void
112   rgbrow(ROPMAT *rop, int r, int p)
# Line 91 | Line 150 | sensrow(ROPMAT *rop, int r, double (*sf)(SCOLOR sc, in
150  
151          for (i = nc; i--; ) {
152                  SCOLOR  sclr;
153 <                scolorblack(sclr);
153 >                memset(sclr, 0, sizeof(COLORV)*nc);
154                  sclr[i] = 1.f;
155                  rop->preop.cmat[r*nc+i] = (*sf)(sclr, nc, rop->mtx->wlpart);
156          }
# Line 103 | Line 162 | checksymbolic(ROPMAT *rop)
162   {
163          const int       nc = rop->mtx->ncomp;
164          const int       dt = rop->mtx->dtype;
165 +        double          cf = 1;
166          int             i, j;
167 +                                        /* check suffix => reference file */
168 +        if (strchr(rop->preop.csym, '.') > rop->preop.csym)
169 +                return(checkreffile(rop));
170  
171          if (nc < 3) {
172                  fprintf(stderr, "%s: -c '%s' requires at least 3 components\n",
# Line 120 | Line 183 | checksymbolic(ROPMAT *rop)
183                  int     comp = 0;
184                  switch (rop->preop.csym[j]) {
185                  case 'B':
186 +                case 'b':
187                          ++comp;
188                          /* fall through */
189                  case 'G':
190 +                case 'g':
191                          ++comp;
192                          /* fall through */
193                  case 'R':
194 +                case 'r':
195 +                        if (rop->preop.csym[j] <= 'Z')
196 +                                cf = 1./WHTEFFICACY;
197                          if (dt == DTxyze) {
198                                  for (i = 3; i--; )
199 <                                        rop->preop.cmat[j*nc+i] = 1./WHTEFFICACY *
132 <                                                        xyz2rgbmat[comp][i];
199 >                                        rop->preop.cmat[j*nc+i] = cf*xyz2rgbmat[comp][i];
200                          } else if (nc == 3)
201                                  rop->preop.cmat[j*nc+comp] = 1.;
202                          else
203                                  rgbrow(rop, j, comp);
204                          break;
205                  case 'Z':
206 +                case 'z':
207                          ++comp;
208                          /* fall through */
209                  case 'Y':
210 +                case 'y':
211                          ++comp;
212                          /* fall through */
213                  case 'X':
214 +                case 'x':
215 +                        if ((rop->preop.csym[j] <= 'Z') & (dt != DTxyze))
216 +                                cf = WHTEFFICACY;
217                          if (dt == DTxyze) {
218                                  rop->preop.cmat[j*nc+comp] = 1.;
219                          } else if (nc == 3) {
# Line 153 | Line 225 | checksymbolic(ROPMAT *rop)
225                          else
226                                  xyzrow(rop, j, comp);
227  
228 <                        for (i = nc*(dt != DTxyze); i--; )
229 <                                rop->preop.cmat[j*nc+i] *= WHTEFFICACY;
228 >                        for (i = nc*(cf != 1); i--; )
229 >                                rop->preop.cmat[j*nc+i] *= cf;
230                          break;
231                  case 'S':               /* scotopic (il)luminance */
232 +                        cf = WHTSCOTOPIC;
233 +                        /* fall through */
234 +                case 's':
235                          sensrow(rop, j, scolor2scotopic);
236 <                        for (i = nc; i--; )
237 <                                rop->preop.cmat[j*nc+i] *= WHTSCOTOPIC;
236 >                        for (i = nc*(cf != 1); i--; )
237 >                                rop->preop.cmat[j*nc+i] *= cf;
238                          break;
239                  case 'M':               /* melanopic (il)luminance */
240 +                        cf = WHTMELANOPIC;
241 +                        /* fall through */
242 +                case 'm':
243                          sensrow(rop, j, scolor2melanopic);
244 <                        for (i = nc; i--; )
245 <                                rop->preop.cmat[j*nc+i] *= WHTMELANOPIC;
244 >                        for (i = nc*(cf != 1); i--; )
245 >                                rop->preop.cmat[j*nc+i] *= cf;
246                          break;
247                  case 'A':               /* average component */
248 +                case 'a':
249                          for (i = nc; i--; )
250                                  rop->preop.cmat[j*nc+i] = 1./(double)nc;
251                          break;
# Line 177 | Line 256 | checksymbolic(ROPMAT *rop)
256                  }
257          }
258                                          /* return recommended output type */
259 <        if (!strcmp(rop->preop.csym, "XYZ")) {
259 >        if (!strcasecmp(rop->preop.csym, "XYZ")) {
260                  if (dt <= DTspec)
261                          return(DTxyze);
262 <        } else if (!strcmp(rop->preop.csym, "RGB")) {
262 >        } else if (!strcasecmp(rop->preop.csym, "RGB")) {
263                  if (dt <= DTspec)
264                          return(DTrgbe);
265 <        }
187 <        if ((nc > 3) & (dt <= DTspec))
265 >        } else if (dt == DTspec)
266                  return(DTfloat);        /* probably not actual spectrum */
267          return(0);
268   }
# Line 200 | Line 278 | loadop(ROPMAT *rop)
278          if (loadmatrix(rop) < 0)                /* make sure we're loaded */
279                  return(NULL);
280  
281 <        if (rop->preop.csym[0] &&               /* symbolic transform? */
281 >        if (rop->preop.csym &&                  /* symbolic transform? */
282                          (outtype = checksymbolic(rop)) < 0)
283                  goto failure;
284          if (rop->preop.clen > 0) {              /* apply component transform? */
# Line 213 | Line 291 | loadop(ROPMAT *rop)
291                          if (rop->preop.nsf == 1) {
292                                  for (i = rop->preop.clen; i--; )
293                                          rop->preop.cmat[i] *= rop->preop.sca[0];
294 <                        } else if (rop->preop.nsf != rop->mtx->ncomp) {
294 >                        } else if (rop->preop.nsf*rop->mtx->ncomp != rop->preop.clen) {
295                                  fprintf(stderr, "%s: -s must have one or %d factors\n",
296 <                                                rop->inspec, rop->mtx->ncomp);
296 >                                                rop->inspec,
297 >                                                rop->preop.clen/rop->mtx->ncomp);
298                                  goto failure;
299                          } else {
300 <                                for (j = rop->preop.clen/rop->preop.nsf; j--; )
301 <                                        for (i = rop->preop.nsf; i--; )
302 <                                                rop->preop.cmat[j*rop->preop.nsf+i] *=
303 <                                                                rop->preop.sca[i];
300 >                                for (i = rop->preop.nsf; i--; )
301 >                                        for (j = rop->mtx->ncomp; j--; )
302 >                                                rop->preop.cmat[i*rop->mtx->ncomp+j]
303 >                                                                *= rop->preop.sca[i];
304                          }
305                  }
306                  mres = rmx_transform(rop->mtx, rop->preop.clen/rop->mtx->ncomp,
# Line 236 | Line 315 | loadop(ROPMAT *rop)
315                                          rop->inspec, mres->ncomp,
316                                          rop->mtx->ncomp,
317                                          rop->preop.nsf ? " (* scalar)" : "");
318 <                rop->preop.nsf = 0;
318 >                rop->preop.nsf = 0;             /* now folded in */
319                  if ((mres->ncomp > 3) & (mres->dtype <= DTspec))
320                          outtype = DTfloat;      /* probably not actual spectrum */
321                  rmx_free(rop->mtx);
# Line 449 | Line 528 | get_factors(double da[], int n, char *av[])
528   }
529  
530   static ROPMAT *
531 < grow_moparray(ROPMAT *mop, int n2alloc)
531 > resize_moparr(ROPMAT *mop, int n2alloc)
532   {
533          int     nmats = 0;
534 +        int     i;
535  
536          while (mop[nmats++].binop)
537                  ;
538 +        for (i = nmats; i > n2alloc; i--)
539 +                rmx_free(mop[i].mtx);
540          mop = (ROPMAT *)realloc(mop, n2alloc*sizeof(ROPMAT));
541          if (mop == NULL) {
542 <                fputs("Out of memory in grow_moparray()\n", stderr);
542 >                fputs("Out of memory in resize_moparr()\n", stderr);
543                  exit(1);
544          }
545          if (n2alloc > nmats)
# Line 469 | Line 551 | grow_moparray(ROPMAT *mop, int n2alloc)
551   int
552   main(int argc, char *argv[])
553   {
554 <        int     outfmt = DTfromHeader;
555 <        int     nall = 2;
556 <        ROPMAT  *mop = (ROPMAT *)calloc(nall, sizeof(ROPMAT));
557 <        int     nmats = 0;
558 <        RMATRIX *mres = NULL;
559 <        int     stdin_used = 0;
560 <        int     i;
554 >        int             outfmt = DTfromHeader;
555 >        const char      *defCsym = NULL;
556 >        int             nall = 2;
557 >        ROPMAT          *mop = (ROPMAT *)calloc(nall, sizeof(ROPMAT));
558 >        int             nmats = 0;
559 >        RMATRIX         *mres = NULL;
560 >        int             stdin_used = 0;
561 >        int             i;
562                                          /* get options and arguments */
563          for (i = 1; i < argc; i++) {
564                  if (argv[i][0] && !argv[i][1] &&
# Line 498 | Line 581 | main(int argc, char *argv[])
581                                  mop[nmats].inspec = stdin_name;
582                          } else
583                                  mop[nmats].inspec = argv[i];
584 +                        if (!mop[nmats].preop.csym)
585 +                                mop[nmats].preop.csym = defCsym;
586                          if (nmats > 0 && !mop[nmats-1].binop)
587                                  mop[nmats-1].binop = '.';
588                          nmats++;
# Line 539 | Line 624 | main(int argc, char *argv[])
624                                          goto userr;
625                                  }
626                                  break;
627 +                        case 'C':
628 +                                if (!n || isflt(argv[i+1]))
629 +                                        goto userr;
630 +                                defCsym = mop[nmats].preop.csym = argv[++i];
631 +                                mop[nmats].preop.clen = 0;
632 +                                break;
633                          case 'c':
634 <                                if (n && isupper(argv[i+1][0])) {
635 <                                        strlcpy(mop[nmats].preop.csym,
545 <                                                argv[++i],
546 <                                                sizeof(mop[0].preop.csym));
634 >                                if (n && !isflt(argv[i+1])) {
635 >                                        mop[nmats].preop.csym = argv[++i];
636                                          mop[nmats].preop.clen = 0;
637                                          break;
638                                  }
# Line 556 | Line 645 | main(int argc, char *argv[])
645                                                          argv[0]);
646                                          goto userr;
647                                  }
648 <                                mop[nmats].preop.csym[0] = '\0';
648 >                                mop[nmats].preop.csym = NULL;
649                                  break;
650                          case 'r':
651                                  if (argv[i][2] == 'f')
# Line 573 | Line 662 | main(int argc, char *argv[])
662                          }
663                  }
664                  if (nmats >= nall)
665 <                        mop = grow_moparray(mop, nall += 2);
665 >                        mop = resize_moparr(mop, nall += 2);
666          }
667          if (mop[0].inspec == NULL)      /* nothing to do? */
668                  goto userr;
# Line 601 | Line 690 | main(int argc, char *argv[])
690                  else if (mres->dtype == DTxyze)
691                          outfmt = DTxyze;
692          }
693 <        if (outfmt != DTascii)          /* write result to stdout */
605 <                SET_FILE_BINARY(stdout);
606 <        newheader("RADIANCE", stdout);
693 >        newheader("RADIANCE", stdout);  /* write result to stdout */
694          printargs(argc, argv, stdout);
608 –
695          return(rmx_write(mres, outfmt, stdout) ? 0 : 1);
696   userr:
697          fprintf(stderr,
698 <        "Usage: %s [-v][-f[adfc][-t][-s sf .. | -c ce ..][-rf|-rb] m1 [.+*/] .. > mres\n",
698 >        "Usage: %s [-v][-f{adfc}][-t][-s sf .. | -c ce ..][-rf|-rb] m1 [.+*/] .. > mres\n",
699                          argv[0]);
700          return(1);
701   }

Diff Legend

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