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

Diff Legend

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