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.23 by greg, Tue Nov 28 16:36:50 2023 UTC vs.
Revision 2.41 by greg, Fri Apr 18 23:59:03 2025 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"
12 – #include "resolu.h"
10   #include "rmatrix.h"
11   #include "platform.h"
12  
13 < #define MAXCOMP         MAXCSAMP        /* #components we support */
13 > /* Preferred BSDF component:
14 >        none, transmission, reflection front (normal side), reflection back */
15 > typedef enum {RMPnone=-1, RMPtrans=0, RMPreflF, RMPreflB} RMPref;
16  
17   /* Unary matrix operation(s) */
18   typedef struct {
20 –        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 37 | 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? */
42 >        if (rop->mtx)                   /* already loaded? */
43                  return(0);
44 +                                        /* check for BSDF input */
45 +        if ((rop->inspec[0] != '!') & (rop->rmp != RMPnone)) {
46 +                const char      *sp = strrchr(rop->inspec, '.');
47 +                if (sp > rop->inspec && !strcasecmp(sp+1, "XML")) {
48 +                        CMATRIX *cm = rop->rmp==RMPtrans ? cm_loadBTDF(rop->inspec) :
49 +                                        cm_loadBRDF(rop->inspec, rop->rmp==RMPreflB) ;
50 +                        if (!cm)
51 +                                return(-1);
52 +                        rop->mtx = rmx_from_cmatrix(cm);
53 +                        cm_free(cm);
54 +                        rop->mtx->dtype = DTascii;
55 +                        return(1);      /* loaded BSDF XML file */
56 +                }
57 +        }                               /* else load regular matrix */
58 +        rop->mtx = rmx_load(rop->inspec);
59  
60 <        rop->mtx = rmx_load(rop->inspec, rop->rmp);
60 >        return(rop->mtx ? 1 : -1);
61 > }
62  
63 <        return(!rop->mtx ? -1 : 1);
63 > extern int      checksymbolic(ROPMAT *rop);
64 >
65 > /* Check/set transform based on a reference input file */
66 > int
67 > checkreffile(ROPMAT *rop)
68 > {
69 >        static const char       *curRF = NULL;
70 >        static RMATRIX          refm;
71 >        const int               nc = rop->mtx->ncomp;
72 >        int                     i;
73 >
74 >        if (!curRF || strcmp(rop->preop.csym, curRF)) {
75 >                FILE    *fp = fopen(rop->preop.csym, "rb");
76 >                if (!rmx_load_header(&refm, fp)) {
77 >                        fprintf(stderr, "%s: cannot read info header\n",
78 >                                        rop->preop.csym);
79 >                        curRF = NULL;
80 >                        if (fp) fclose(fp);
81 >                        return(-1);
82 >                }
83 >                fclose(fp);
84 >                curRF = rop->preop.csym;
85 >        }
86 >        if (refm.ncomp == 3) {
87 >                rop->preop.csym = (refm.dtype == DTxyze) ? "XYZ" : "RGB";
88 >                return(checksymbolic(rop));
89 >        }
90 >        if (refm.ncomp == 2) {
91 >                fprintf(stderr, "%s: cannot convert to 2 components\n",
92 >                                curRF);
93 >                return(-1);
94 >        }
95 >        if (refm.ncomp == 1) {
96 >                rop->preop.csym = "Y";          /* XXX big assumption */
97 >                return(checksymbolic(rop));
98 >        }
99 >        if (refm.ncomp == nc &&
100 >                        !memcmp(refm.wlpart, rop->mtx->wlpart, sizeof(refm.wlpart)))
101 >                return(0);                      /* nothing to do */
102 >
103 >        if ((nc <= 3) | (nc > MAXCSAMP) | (refm.ncomp > MAXCSAMP)) {
104 >                fprintf(stderr, "%s: cannot resample from %d to %d components\n",
105 >                                curRF, nc, refm.ncomp);
106 >                return(-1);
107 >        }
108 >        rop->preop.clen = refm.ncomp * nc;      /* compute spec to ref */
109 >
110 >        for (i = 0; i < nc; i++) {
111 >                SCOLOR  scstim, scresp;
112 >                int     j;
113 >                memset(scstim, 0, sizeof(COLORV)*nc);
114 >                scstim[i] = 1.f;
115 >                convertscolor(scresp, refm.ncomp, refm.wlpart[0], refm.wlpart[3],
116 >                                scstim, nc, rop->mtx->wlpart[0], rop->mtx->wlpart[3]);
117 >                for (j = refm.ncomp; j-- > 0; )
118 >                        rop->preop.cmat[j*nc + i] = scresp[j];
119 >        }
120 >        memcpy(rop->mtx->wlpart, refm.wlpart, sizeof(rop->mtx->wlpart));
121 >        return(0);
122   }
123  
124   /* Compute conversion row from spectrum to one channel of RGB */
125 < static void
125 > void
126   rgbrow(ROPMAT *rop, int r, int p)
127   {
128          const int       nc = rop->mtx->ncomp;
# Line 66 | Line 139 | rgbrow(ROPMAT *rop, int r, int p)
139   }
140  
141   /* Compute conversion row from spectrum to one channel of XYZ */
142 < static void
142 > void
143   xyzrow(ROPMAT *rop, int r, int p)
144   {
145          const int       nc = rop->mtx->ncomp;
# Line 83 | Line 156 | xyzrow(ROPMAT *rop, int r, int p)
156   }
157  
158   /* Use the spectral sensitivity function to compute matrix coefficients */
159 < static void
160 < sensrow(ROPMAT *rop, int r, double (*sf)(SCOLOR sc, int ncs, const float wlpt[4]))
159 > void
160 > sensrow(ROPMAT *rop, int r, double (*sf)(const SCOLOR sc, int ncs, const float wlpt[4]))
161   {
162          const int       nc = rop->mtx->ncomp;
163          int             i;
164  
165          for (i = nc; i--; ) {
166                  SCOLOR  sclr;
167 <                scolorblack(sclr);
168 <                sclr[i] = 1;
167 >                memset(sclr, 0, sizeof(COLORV)*nc);
168 >                sclr[i] = 1.f;
169                  rop->preop.cmat[r*nc+i] = (*sf)(sclr, nc, rop->mtx->wlpart);
170          }
171   }
172  
173   /* Check/set symbolic transform */
174 < static int
174 > int
175   checksymbolic(ROPMAT *rop)
176   {
177          const int       nc = rop->mtx->ncomp;
178          const int       dt = rop->mtx->dtype;
179 +        double          cf = 1;
180          int             i, j;
181 +                                        /* check suffix => reference file */
182 +        if (strchr(rop->preop.csym, '.') > rop->preop.csym)
183 +                return(checkreffile(rop));
184  
185          if (nc < 3) {
186                  fprintf(stderr, "%s: -c '%s' requires at least 3 components\n",
# Line 120 | Line 197 | checksymbolic(ROPMAT *rop)
197                  int     comp = 0;
198                  switch (rop->preop.csym[j]) {
199                  case 'B':
200 +                case 'b':
201                          ++comp;
202                          /* fall through */
203                  case 'G':
204 +                case 'g':
205                          ++comp;
206                          /* fall through */
207                  case 'R':
208 +                case 'r':
209 +                        if (rop->preop.csym[j] <= 'Z')
210 +                                cf = 1./WHTEFFICACY;
211                          if (dt == DTxyze) {
212                                  for (i = 3; i--; )
213 <                                        rop->preop.cmat[j*nc+i] = 1./WHTEFFICACY *
132 <                                                        xyz2rgbmat[comp][i];
213 >                                        rop->preop.cmat[j*nc+i] = cf*xyz2rgbmat[comp][i];
214                          } else if (nc == 3)
215                                  rop->preop.cmat[j*nc+comp] = 1.;
216                          else
217                                  rgbrow(rop, j, comp);
218                          break;
219                  case 'Z':
220 +                case 'z':
221                          ++comp;
222                          /* fall through */
223                  case 'Y':
224 +                case 'y':
225                          ++comp;
226                          /* fall through */
227                  case 'X':
228 +                case 'x':
229 +                        if ((rop->preop.csym[j] <= 'Z') & (dt != DTxyze))
230 +                                cf = WHTEFFICACY;
231                          if (dt == DTxyze) {
232                                  rop->preop.cmat[j*nc+comp] = 1.;
233                          } else if (nc == 3) {
# Line 153 | Line 239 | checksymbolic(ROPMAT *rop)
239                          else
240                                  xyzrow(rop, j, comp);
241  
242 <                        for (i = nc*(dt != DTxyze); i--; )
243 <                                rop->preop.cmat[j*nc+i] *= WHTEFFICACY;
242 >                        for (i = nc*(cf != 1); i--; )
243 >                                rop->preop.cmat[j*nc+i] *= cf;
244                          break;
245 <                case 'S':
245 >                case 'S':               /* scotopic (il)luminance */
246 >                        cf = WHTSCOTOPIC;
247 >                        /* fall through */
248 >                case 's':
249                          sensrow(rop, j, scolor2scotopic);
250 <                        for (i = nc; i--; )
251 <                                rop->preop.cmat[j*nc+i] *= WHTSCOTOPIC;
250 >                        for (i = nc*(cf != 1); i--; )
251 >                                rop->preop.cmat[j*nc+i] *= cf;
252                          break;
253 <                case 'M':
253 >                case 'M':               /* melanopic (il)luminance */
254 >                        cf = WHTMELANOPIC;
255 >                        /* fall through */
256 >                case 'm':
257                          sensrow(rop, j, scolor2melanopic);
258 +                        for (i = nc*(cf != 1); i--; )
259 +                                rop->preop.cmat[j*nc+i] *= cf;
260 +                        break;
261 +                case 'A':               /* average component */
262 +                case 'a':
263                          for (i = nc; i--; )
264 <                                rop->preop.cmat[j*nc+i] *= WHTMELANOPIC;
264 >                                rop->preop.cmat[j*nc+i] = 1./(double)nc;
265                          break;
266                  default:
267                          fprintf(stderr, "%s: -c '%c' unsupported\n",
# Line 173 | Line 270 | checksymbolic(ROPMAT *rop)
270                  }
271          }
272                                          /* return recommended output type */
273 <        if (!strcmp(rop->preop.csym, "XYZ")) {
273 >        if (!strcasecmp(rop->preop.csym, "XYZ")) {
274                  if (dt <= DTspec)
275                          return(DTxyze);
276 <        } else if (!strcmp(rop->preop.csym, "RGB")) {
276 >        } else if (!strcasecmp(rop->preop.csym, "RGB")) {
277                  if (dt <= DTspec)
278                          return(DTrgbe);
279 <        }
183 <        if ((nc > 3) & (dt <= DTspec))
279 >        } else if (dt == DTspec)
280                  return(DTfloat);        /* probably not actual spectrum */
281          return(0);
282   }
283  
284   /* Get matrix and perform unary operations */
285 < static RMATRIX *
285 > RMATRIX *
286   loadop(ROPMAT *rop)
287   {
288          int     outtype = 0;
# Line 196 | Line 292 | loadop(ROPMAT *rop)
292          if (loadmatrix(rop) < 0)                /* make sure we're loaded */
293                  return(NULL);
294  
295 <        if (rop->preop.csym[0] &&               /* symbolic transform? */
295 >        if (rop->preop.csym &&                  /* symbolic transform? */
296                          (outtype = checksymbolic(rop)) < 0)
297                  goto failure;
298          if (rop->preop.clen > 0) {              /* apply component transform? */
# Line 209 | Line 305 | loadop(ROPMAT *rop)
305                          if (rop->preop.nsf == 1) {
306                                  for (i = rop->preop.clen; i--; )
307                                          rop->preop.cmat[i] *= rop->preop.sca[0];
308 <                        } else if (rop->preop.nsf != rop->mtx->ncomp) {
308 >                        } else if (rop->preop.nsf*rop->mtx->ncomp != rop->preop.clen) {
309                                  fprintf(stderr, "%s: -s must have one or %d factors\n",
310 <                                                rop->inspec, rop->mtx->ncomp);
310 >                                                rop->inspec,
311 >                                                rop->preop.clen/rop->mtx->ncomp);
312                                  goto failure;
313                          } else {
314 <                                for (j = rop->preop.clen/rop->preop.nsf; j--; )
315 <                                        for (i = rop->preop.nsf; i--; )
316 <                                                rop->preop.cmat[j*rop->preop.nsf+i] *=
317 <                                                                rop->preop.sca[i];
314 >                                for (i = rop->preop.nsf; i--; )
315 >                                        for (j = rop->mtx->ncomp; j--; )
316 >                                                rop->preop.cmat[i*rop->mtx->ncomp+j]
317 >                                                                *= rop->preop.sca[i];
318                          }
319                  }
320                  mres = rmx_transform(rop->mtx, rop->preop.clen/rop->mtx->ncomp,
# Line 232 | Line 329 | loadop(ROPMAT *rop)
329                                          rop->inspec, mres->ncomp,
330                                          rop->mtx->ncomp,
331                                          rop->preop.nsf ? " (* scalar)" : "");
332 <                rop->preop.nsf = 0;
332 >                rop->preop.nsf = 0;             /* now folded in */
333                  if ((mres->ncomp > 3) & (mres->dtype <= DTspec))
334                          outtype = DTfloat;      /* probably not actual spectrum */
335                  rmx_free(rop->mtx);
# Line 261 | Line 358 | loadop(ROPMAT *rop)
358                  }
359          }
360          if (rop->preop.transpose) {             /* transpose matrix? */
361 <                mres = rmx_transpose(rop->mtx);
265 <                if (mres == NULL) {
361 >                if (!rmx_transpose(rop->mtx)) {
362                          fputs(rop->inspec, stderr);
363                          fputs(": transpose failed\n", stderr);
364                          goto failure;
# Line 271 | Line 367 | loadop(ROPMAT *rop)
367                          fputs(rop->inspec, stderr);
368                          fputs(": transposed rows and columns\n", stderr);
369                  }
274 –                rmx_free(rop->mtx);
275 –                rop->mtx = mres;
370          }
371          mres = rop->mtx;
372          rop->mtx = NULL;
# Line 285 | Line 379 | failure:
379   }
380  
381   /* Execute binary operation, free matrix arguments and return new result */
382 < static RMATRIX *
382 > RMATRIX *
383   binaryop(const char *inspec, RMATRIX *mleft, int op, RMATRIX *mright)
384   {
385          RMATRIX *mres = NULL;
# Line 362 | Line 456 | binaryop(const char *inspec, RMATRIX *mleft, int op, R
456   }
457  
458   /* Perform matrix operations from left to right */
459 < static RMATRIX *
459 > RMATRIX *
460   op_left2right(ROPMAT *mop)
461   {
462          RMATRIX *mleft = loadop(mop);
# Line 378 | Line 472 | op_left2right(ROPMAT *mop)
472   }
473  
474   /* Perform matrix operations from right to left */
475 < static RMATRIX *
475 > RMATRIX *
476   op_right2left(ROPMAT *mop)
477   {
478          RMATRIX *mright;
# Line 407 | Line 501 | op_right2left(ROPMAT *mop)
501                                                  : (mop)->mtx->ncols)
502  
503   /* Should we prefer concatenating from rightmost matrix towards left? */
504 < static int
504 > int
505   prefer_right2left(ROPMAT *mop)
506   {
507          int     mri = 0;
# Line 434 | Line 528 | prefer_right2left(ROPMAT *mop)
528          return(t_ncols(mop+mri) < t_nrows(mop));
529   }
530  
531 < static int
531 > int
532   get_factors(double da[], int n, char *av[])
533   {
534          int     ac;
# Line 444 | Line 538 | get_factors(double da[], int n, char *av[])
538          return(ac);
539   }
540  
541 < static ROPMAT *
542 < grow_moparray(ROPMAT *mop, int n2alloc)
541 > ROPMAT *
542 > resize_moparr(ROPMAT *mop, int n2alloc)
543   {
544          int     nmats = 0;
545 +        int     i;
546  
547          while (mop[nmats++].binop)
548                  ;
549 +        for (i = nmats; i >= n2alloc; i--)
550 +                rmx_free(mop[i].mtx);
551          mop = (ROPMAT *)realloc(mop, n2alloc*sizeof(ROPMAT));
552          if (mop == NULL) {
553 <                fputs("Out of memory in grow_moparray()\n", stderr);
553 >                fputs("Out of memory in resize_moparr()\n", stderr);
554                  exit(1);
555          }
556          if (n2alloc > nmats)
# Line 465 | Line 562 | grow_moparray(ROPMAT *mop, int n2alloc)
562   int
563   main(int argc, char *argv[])
564   {
565 <        int     outfmt = DTfromHeader;
566 <        int     nall = 2;
567 <        ROPMAT  *mop = (ROPMAT *)calloc(nall, sizeof(ROPMAT));
568 <        int     nmats = 0;
569 <        RMATRIX *mres = NULL;
570 <        int     stdin_used = 0;
571 <        int     i;
565 >        int             outfmt = DTfromHeader;
566 >        const char      *defCsym = NULL;
567 >        int             nall = 2;
568 >        ROPMAT          *mop = (ROPMAT *)calloc(nall, sizeof(ROPMAT));
569 >        int             nmats = 0;
570 >        RMATRIX         *mres = NULL;
571 >        int             stdin_used = 0;
572 >        int             i;
573                                          /* get options and arguments */
574          for (i = 1; i < argc; i++) {
575                  if (argv[i][0] && !argv[i][1] &&
# Line 494 | Line 592 | main(int argc, char *argv[])
592                                  mop[nmats].inspec = stdin_name;
593                          } else
594                                  mop[nmats].inspec = argv[i];
595 +                        if (!mop[nmats].preop.csym)
596 +                                mop[nmats].preop.csym = defCsym;
597                          if (nmats > 0 && !mop[nmats-1].binop)
598                                  mop[nmats-1].binop = '.';
599                          nmats++;
# Line 535 | Line 635 | main(int argc, char *argv[])
635                                          goto userr;
636                                  }
637                                  break;
638 +                        case 'C':
639 +                                if (!n || isflt(argv[i+1]))
640 +                                        goto userr;
641 +                                defCsym = mop[nmats].preop.csym = argv[++i];
642 +                                mop[nmats].preop.clen = 0;
643 +                                break;
644                          case 'c':
645 <                                if (n && isupper(argv[i+1][0])) {
646 <                                        strlcpy(mop[nmats].preop.csym,
541 <                                                argv[++i],
542 <                                                sizeof(mop[0].preop.csym));
645 >                                if (n && !isflt(argv[i+1])) {
646 >                                        mop[nmats].preop.csym = argv[++i];
647                                          mop[nmats].preop.clen = 0;
648                                          break;
649                                  }
# Line 552 | Line 656 | main(int argc, char *argv[])
656                                                          argv[0]);
657                                          goto userr;
658                                  }
659 <                                mop[nmats].preop.csym[0] = '\0';
659 >                                mop[nmats].preop.csym = NULL;
660                                  break;
661                          case 'r':
662                                  if (argv[i][2] == 'f')
# Line 569 | Line 673 | main(int argc, char *argv[])
673                          }
674                  }
675                  if (nmats >= nall)
676 <                        mop = grow_moparray(mop, nall += 2);
676 >                        mop = resize_moparr(mop, nall += 2);
677          }
678          if (mop[0].inspec == NULL)      /* nothing to do? */
679                  goto userr;
# Line 589 | Line 693 | main(int argc, char *argv[])
693          mres = loadop(mop+nmats);
694          if (mres == NULL)
695                  return(1);
696 <        if (outfmt == DTfromHeader)     /* check data type */
696 >        if ((outfmt == DTfromHeader) & (mres->dtype < DTspec))
697                  outfmt = mres->dtype;
698 <        if (outfmt == DTrgbe) {
698 >        if (outfmt == DTrgbe) {         /* check data type */
699                  if (mres->ncomp > 3)
700                          outfmt = DTspec;
701                  else if (mres->dtype == DTxyze)
702                          outfmt = DTxyze;
703          }
704 <        if (outfmt != DTascii)          /* write result to stdout */
705 <                SET_FILE_BINARY(stdout);
706 <        newheader("RADIANCE", stdout);
704 > #if DTrmx_native==DTfloat
705 >        if (outfmt == DTdouble)
706 >                fprintf(stderr,
707 >                        "%s: warning - writing float result as double\n",
708 >                                argv[0]);
709 > #endif
710 >        newheader("RADIANCE", stdout);  /* write result to stdout */
711          printargs(argc, argv, stdout);
604 –
712          return(rmx_write(mres, outfmt, stdout) ? 0 : 1);
713   userr:
714          fprintf(stderr,
715 <        "Usage: %s [-v][-f[adfc][-t][-s sf .. | -c ce ..][-rf|-rb] m1 [.+*/] .. > mres\n",
715 >        "Usage: %s [-v][-f{adfc}][-t][-s sf .. | -c ce ..][-rf|-rb] m1 [.+*/] .. > mres\n",
716                          argv[0]);
717          return(1);
718   }

Diff Legend

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