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.26 by greg, Fri Dec 1 02:05:00 2023 UTC vs.
Revision 2.35 by greg, Fri Nov 8 17:52:26 2024 UTC

# Line 6 | Line 6 | static const char RCSid[] = "$Id$";
6   */
7  
8   #include <errno.h>
9 – #include <ctype.h>
9   #include "rtio.h"
11 – #include "resolu.h"
10   #include "rmatrix.h"
11   #include "platform.h"
12  
# Line 16 | Line 14 | static const char RCSid[] = "$Id$";
14  
15   /* Unary matrix operation(s) */
16   typedef struct {
19 –        double          sca[MAXCOMP];           /* scalar coefficients */
17          double          cmat[MAXCOMP*MAXCOMP];  /* component transformation */
18 <        short           nsf;                    /* number of scalars */
18 >        double          sca[MAXCOMP];           /* scalar coefficients */
19 >        const char      *csym;                  /* symbolic coefs or file */
20          short           clen;                   /* number of coefficients */
21 <        char            csym[11];               /* symbolic coefficients */
22 <        char            transpose;              /* do transpose? */
21 >        short           nsf;                    /* number of scalars */
22 >        short           transpose;              /* do transpose? */
23   } RUNARYOP;
24  
25   /* Matrix input source and requested operation(s) */
# Line 36 | Line 34 | typedef struct {
34   int     verbose = 0;                    /* verbose reporting? */
35  
36   /* Load matrix */
37 < static int
37 > int
38   loadmatrix(ROPMAT *rop)
39   {
40          if (rop->mtx != NULL)           /* already loaded? */
# Line 47 | Line 45 | loadmatrix(ROPMAT *rop)
45          return(!rop->mtx ? -1 : 1);
46   }
47  
48 + extern int      checksymbolic(ROPMAT *rop);
49 +
50 + /* Check/set transform based on a reference input file */
51 + int
52 + checkreffile(ROPMAT *rop)
53 + {
54 +        static const char       *curRF = NULL;
55 +        static RMATRIX          refm;
56 +        const int               nc = rop->mtx->ncomp;
57 +        int                     i;
58 +
59 +        if (!curRF || strcmp(rop->preop.csym, curRF)) {
60 +                FILE    *fp = fopen(rop->preop.csym, "rb");
61 +                if (!rmx_load_header(&refm, fp)) {
62 +                        fprintf(stderr, "%s: cannot read info header\n",
63 +                                        rop->preop.csym);
64 +                        curRF = NULL;
65 +                        if (fp) fclose(fp);
66 +                        return(-1);
67 +                }
68 +                fclose(fp);
69 +                curRF = rop->preop.csym;
70 +        }
71 +        if (refm.ncomp == 3) {
72 +                rop->preop.csym = (refm.dtype == DTxyze) ? "XYZ" : "RGB";
73 +                return(checksymbolic(rop));
74 +        }
75 +        if (refm.ncomp == 2) {
76 +                fprintf(stderr, "%s: cannot convert to 2 components\n",
77 +                                curRF);
78 +                return(-1);
79 +        }
80 +        if (refm.ncomp == 1) {
81 +                rop->preop.csym = "Y";          /* XXX big assumption */
82 +                return(checksymbolic(rop));
83 +        }
84 +        if (refm.ncomp == nc &&
85 +                        !memcmp(refm.wlpart, rop->mtx->wlpart, sizeof(refm.wlpart)))
86 +                return(0);                      /* nothing to do */
87 +
88 +        if ((nc <= 3) | (nc > MAXCSAMP) | (refm.ncomp > MAXCSAMP)) {
89 +                fprintf(stderr, "%s: cannot resample from %d to %d components\n",
90 +                                curRF, nc, refm.ncomp);
91 +                return(-1);
92 +        }
93 +        rop->preop.clen = refm.ncomp * nc;      /* compute spec to ref */
94 +
95 +        for (i = 0; i < nc; i++) {
96 +                SCOLOR  scstim, scresp;
97 +                int     j;
98 +                memset(scstim, 0, sizeof(COLORV)*nc);
99 +                scstim[i] = 1.f;
100 +                convertscolor(scresp, refm.ncomp, refm.wlpart[0], refm.wlpart[3],
101 +                                scstim, nc, rop->mtx->wlpart[0], rop->mtx->wlpart[3]);
102 +                for (j = refm.ncomp; j-- > 0; )
103 +                        rop->preop.cmat[j*nc + i] = scresp[j];
104 +        }
105 +        memcpy(rop->mtx->wlpart, refm.wlpart, sizeof(rop->mtx->wlpart));
106 +        return(0);
107 + }
108 +
109   /* Compute conversion row from spectrum to one channel of RGB */
110 < static void
110 > void
111   rgbrow(ROPMAT *rop, int r, int p)
112   {
113          const int       nc = rop->mtx->ncomp;
# Line 65 | Line 124 | rgbrow(ROPMAT *rop, int r, int p)
124   }
125  
126   /* Compute conversion row from spectrum to one channel of XYZ */
127 < static void
127 > void
128   xyzrow(ROPMAT *rop, int r, int p)
129   {
130          const int       nc = rop->mtx->ncomp;
# Line 82 | Line 141 | xyzrow(ROPMAT *rop, int r, int p)
141   }
142  
143   /* Use the spectral sensitivity function to compute matrix coefficients */
144 < static void
145 < sensrow(ROPMAT *rop, int r, double (*sf)(SCOLOR sc, int ncs, const float wlpt[4]))
144 > void
145 > sensrow(ROPMAT *rop, int r, double (*sf)(const SCOLOR sc, int ncs, const float wlpt[4]))
146   {
147          const int       nc = rop->mtx->ncomp;
148          int             i;
149  
150          for (i = nc; i--; ) {
151                  SCOLOR  sclr;
152 <                scolorblack(sclr);
152 >                memset(sclr, 0, sizeof(COLORV)*nc);
153                  sclr[i] = 1.f;
154                  rop->preop.cmat[r*nc+i] = (*sf)(sclr, nc, rop->mtx->wlpart);
155          }
156   }
157  
158   /* Check/set symbolic transform */
159 < static int
159 > int
160   checksymbolic(ROPMAT *rop)
161   {
162          const int       nc = rop->mtx->ncomp;
163          const int       dt = rop->mtx->dtype;
164 +        double          cf = 1;
165          int             i, j;
166 +                                        /* check suffix => reference file */
167 +        if (strchr(rop->preop.csym, '.') > rop->preop.csym)
168 +                return(checkreffile(rop));
169  
170          if (nc < 3) {
171                  fprintf(stderr, "%s: -c '%s' requires at least 3 components\n",
# Line 119 | Line 182 | checksymbolic(ROPMAT *rop)
182                  int     comp = 0;
183                  switch (rop->preop.csym[j]) {
184                  case 'B':
185 +                case 'b':
186                          ++comp;
187                          /* fall through */
188                  case 'G':
189 +                case 'g':
190                          ++comp;
191                          /* fall through */
192                  case 'R':
193 +                case 'r':
194 +                        if (rop->preop.csym[j] <= 'Z')
195 +                                cf = 1./WHTEFFICACY;
196                          if (dt == DTxyze) {
197                                  for (i = 3; i--; )
198 <                                        rop->preop.cmat[j*nc+i] = 1./WHTEFFICACY *
131 <                                                        xyz2rgbmat[comp][i];
198 >                                        rop->preop.cmat[j*nc+i] = cf*xyz2rgbmat[comp][i];
199                          } else if (nc == 3)
200                                  rop->preop.cmat[j*nc+comp] = 1.;
201                          else
202                                  rgbrow(rop, j, comp);
203                          break;
204                  case 'Z':
205 +                case 'z':
206                          ++comp;
207                          /* fall through */
208                  case 'Y':
209 +                case 'y':
210                          ++comp;
211                          /* fall through */
212                  case 'X':
213 +                case 'x':
214 +                        if ((rop->preop.csym[j] <= 'Z') & (dt != DTxyze))
215 +                                cf = WHTEFFICACY;
216                          if (dt == DTxyze) {
217                                  rop->preop.cmat[j*nc+comp] = 1.;
218                          } else if (nc == 3) {
# Line 152 | Line 224 | checksymbolic(ROPMAT *rop)
224                          else
225                                  xyzrow(rop, j, comp);
226  
227 <                        for (i = nc*(dt != DTxyze); i--; )
228 <                                rop->preop.cmat[j*nc+i] *= WHTEFFICACY;
227 >                        for (i = nc*(cf != 1); i--; )
228 >                                rop->preop.cmat[j*nc+i] *= cf;
229                          break;
230                  case 'S':               /* scotopic (il)luminance */
231 +                        cf = WHTSCOTOPIC;
232 +                        /* fall through */
233 +                case 's':
234                          sensrow(rop, j, scolor2scotopic);
235 <                        for (i = nc; i--; )
236 <                                rop->preop.cmat[j*nc+i] *= WHTSCOTOPIC;
235 >                        for (i = nc*(cf != 1); i--; )
236 >                                rop->preop.cmat[j*nc+i] *= cf;
237                          break;
238                  case 'M':               /* melanopic (il)luminance */
239 +                        cf = WHTMELANOPIC;
240 +                        /* fall through */
241 +                case 'm':
242                          sensrow(rop, j, scolor2melanopic);
243 <                        for (i = nc; i--; )
244 <                                rop->preop.cmat[j*nc+i] *= WHTMELANOPIC;
243 >                        for (i = nc*(cf != 1); i--; )
244 >                                rop->preop.cmat[j*nc+i] *= cf;
245                          break;
246                  case 'A':               /* average component */
247 +                case 'a':
248                          for (i = nc; i--; )
249                                  rop->preop.cmat[j*nc+i] = 1./(double)nc;
250                          break;
# Line 176 | Line 255 | checksymbolic(ROPMAT *rop)
255                  }
256          }
257                                          /* return recommended output type */
258 <        if (!strcmp(rop->preop.csym, "XYZ")) {
258 >        if (!strcasecmp(rop->preop.csym, "XYZ")) {
259                  if (dt <= DTspec)
260                          return(DTxyze);
261 <        } else if (!strcmp(rop->preop.csym, "RGB")) {
261 >        } else if (!strcasecmp(rop->preop.csym, "RGB")) {
262                  if (dt <= DTspec)
263                          return(DTrgbe);
264 <        }
186 <        if ((nc > 3) & (dt <= DTspec))
264 >        } else if (dt == DTspec)
265                  return(DTfloat);        /* probably not actual spectrum */
266          return(0);
267   }
268  
269   /* Get matrix and perform unary operations */
270 < static RMATRIX *
270 > RMATRIX *
271   loadop(ROPMAT *rop)
272   {
273          int     outtype = 0;
# Line 199 | Line 277 | loadop(ROPMAT *rop)
277          if (loadmatrix(rop) < 0)                /* make sure we're loaded */
278                  return(NULL);
279  
280 <        if (rop->preop.csym[0] &&               /* symbolic transform? */
280 >        if (rop->preop.csym &&                  /* symbolic transform? */
281                          (outtype = checksymbolic(rop)) < 0)
282                  goto failure;
283          if (rop->preop.clen > 0) {              /* apply component transform? */
# Line 212 | Line 290 | loadop(ROPMAT *rop)
290                          if (rop->preop.nsf == 1) {
291                                  for (i = rop->preop.clen; i--; )
292                                          rop->preop.cmat[i] *= rop->preop.sca[0];
293 <                        } else if (rop->preop.nsf != rop->mtx->ncomp) {
293 >                        } else if (rop->preop.nsf*rop->mtx->ncomp != rop->preop.clen) {
294                                  fprintf(stderr, "%s: -s must have one or %d factors\n",
295 <                                                rop->inspec, rop->mtx->ncomp);
295 >                                                rop->inspec,
296 >                                                rop->preop.clen/rop->mtx->ncomp);
297                                  goto failure;
298                          } else {
299 <                                for (j = rop->preop.clen/rop->preop.nsf; j--; )
300 <                                        for (i = rop->preop.nsf; i--; )
301 <                                                rop->preop.cmat[j*rop->preop.nsf+i] *=
302 <                                                                rop->preop.sca[i];
299 >                                for (i = rop->preop.nsf; i--; )
300 >                                        for (j = rop->mtx->ncomp; j--; )
301 >                                                rop->preop.cmat[i*rop->mtx->ncomp+j]
302 >                                                                *= rop->preop.sca[i];
303                          }
304                  }
305                  mres = rmx_transform(rop->mtx, rop->preop.clen/rop->mtx->ncomp,
# Line 288 | Line 367 | failure:
367   }
368  
369   /* Execute binary operation, free matrix arguments and return new result */
370 < static RMATRIX *
370 > RMATRIX *
371   binaryop(const char *inspec, RMATRIX *mleft, int op, RMATRIX *mright)
372   {
373          RMATRIX *mres = NULL;
# Line 365 | Line 444 | binaryop(const char *inspec, RMATRIX *mleft, int op, R
444   }
445  
446   /* Perform matrix operations from left to right */
447 < static RMATRIX *
447 > RMATRIX *
448   op_left2right(ROPMAT *mop)
449   {
450          RMATRIX *mleft = loadop(mop);
# Line 381 | Line 460 | op_left2right(ROPMAT *mop)
460   }
461  
462   /* Perform matrix operations from right to left */
463 < static RMATRIX *
463 > RMATRIX *
464   op_right2left(ROPMAT *mop)
465   {
466          RMATRIX *mright;
# Line 410 | Line 489 | op_right2left(ROPMAT *mop)
489                                                  : (mop)->mtx->ncols)
490  
491   /* Should we prefer concatenating from rightmost matrix towards left? */
492 < static int
492 > int
493   prefer_right2left(ROPMAT *mop)
494   {
495          int     mri = 0;
# Line 437 | Line 516 | prefer_right2left(ROPMAT *mop)
516          return(t_ncols(mop+mri) < t_nrows(mop));
517   }
518  
519 < static int
519 > int
520   get_factors(double da[], int n, char *av[])
521   {
522          int     ac;
# Line 447 | Line 526 | get_factors(double da[], int n, char *av[])
526          return(ac);
527   }
528  
529 < static ROPMAT *
529 > ROPMAT *
530   resize_moparr(ROPMAT *mop, int n2alloc)
531   {
532          int     nmats = 0;
# Line 455 | Line 534 | resize_moparr(ROPMAT *mop, int n2alloc)
534  
535          while (mop[nmats++].binop)
536                  ;
537 <        for (i = nmats; i > n2alloc; i--)
537 >        for (i = nmats; i >= n2alloc; i--)
538                  rmx_free(mop[i].mtx);
539          mop = (ROPMAT *)realloc(mop, n2alloc*sizeof(ROPMAT));
540          if (mop == NULL) {
# Line 471 | Line 550 | resize_moparr(ROPMAT *mop, int n2alloc)
550   int
551   main(int argc, char *argv[])
552   {
553 <        int     outfmt = DTfromHeader;
554 <        int     nall = 2;
555 <        ROPMAT  *mop = (ROPMAT *)calloc(nall, sizeof(ROPMAT));
556 <        int     nmats = 0;
557 <        RMATRIX *mres = NULL;
558 <        int     stdin_used = 0;
559 <        int     i;
553 >        int             outfmt = DTfromHeader;
554 >        const char      *defCsym = NULL;
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;
561                                          /* get options and arguments */
562          for (i = 1; i < argc; i++) {
563                  if (argv[i][0] && !argv[i][1] &&
# Line 500 | Line 580 | main(int argc, char *argv[])
580                                  mop[nmats].inspec = stdin_name;
581                          } else
582                                  mop[nmats].inspec = argv[i];
583 +                        if (!mop[nmats].preop.csym)
584 +                                mop[nmats].preop.csym = defCsym;
585                          if (nmats > 0 && !mop[nmats-1].binop)
586                                  mop[nmats-1].binop = '.';
587                          nmats++;
# Line 541 | Line 623 | main(int argc, char *argv[])
623                                          goto userr;
624                                  }
625                                  break;
626 +                        case 'C':
627 +                                if (!n || isflt(argv[i+1]))
628 +                                        goto userr;
629 +                                defCsym = mop[nmats].preop.csym = argv[++i];
630 +                                mop[nmats].preop.clen = 0;
631 +                                break;
632                          case 'c':
633 <                                if (n && isupper(argv[i+1][0])) {
634 <                                        strlcpy(mop[nmats].preop.csym,
547 <                                                argv[++i],
548 <                                                sizeof(mop[0].preop.csym));
633 >                                if (n && !isflt(argv[i+1])) {
634 >                                        mop[nmats].preop.csym = argv[++i];
635                                          mop[nmats].preop.clen = 0;
636                                          break;
637                                  }
# Line 558 | Line 644 | main(int argc, char *argv[])
644                                                          argv[0]);
645                                          goto userr;
646                                  }
647 <                                mop[nmats].preop.csym[0] = '\0';
647 >                                mop[nmats].preop.csym = NULL;
648                                  break;
649                          case 'r':
650                                  if (argv[i][2] == 'f')

Diff Legend

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