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.21 by greg, Tue Nov 21 01:30:20 2023 UTC vs.
Revision 2.39 by greg, Fri Apr 4 18:18:06 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>
9   #include "rtio.h"
11 – #include "resolu.h"
10   #include "rmatrix.h"
11   #include "platform.h"
12  
13 < #define MAXCOMP         16              /* #components we support */
16 <
17 < /* unary matrix operation(s) */
13 > /* Unary matrix operation(s) */
14   typedef struct {
19 –        double          sca[MAXCOMP];           /* scalar coefficients */
15          double          cmat[MAXCOMP*MAXCOMP];  /* component transformation */
16 <        short           nsf;                    /* number of scalars */
16 >        double          sca[MAXCOMP];           /* scalar coefficients */
17 >        const char      *csym;                  /* symbolic coefs or file */
18          short           clen;                   /* number of coefficients */
19 +        short           nsf;                    /* number of scalars */
20          short           transpose;              /* do transpose? */
21   } RUNARYOP;
22  
23 < /* matrix input source and requested operation(s) */
23 > /* Matrix input source and requested operation(s) */
24   typedef struct {
25          const char      *inspec;                /* input specification */
26          RMPref          rmp;                    /* matrix preference */
# Line 35 | Line 32 | typedef struct {
32   int     verbose = 0;                    /* verbose reporting? */
33  
34   /* Load matrix */
35 < static int
35 > int
36   loadmatrix(ROPMAT *rop)
37   {
38          if (rop->mtx != NULL)           /* already loaded? */
39                  return(0);
40  
41          rop->mtx = rmx_load(rop->inspec, rop->rmp);
42 <        if (rop->mtx == NULL) {
43 <                fputs(rop->inspec, stderr);
44 <                fputs(": cannot load matrix\n", stderr);
42 >
43 >        return(!rop->mtx ? -1 : 1);
44 > }
45 >
46 > extern int      checksymbolic(ROPMAT *rop);
47 >
48 > /* Check/set transform based on a reference input file */
49 > int
50 > checkreffile(ROPMAT *rop)
51 > {
52 >        static const char       *curRF = NULL;
53 >        static RMATRIX          refm;
54 >        const int               nc = rop->mtx->ncomp;
55 >        int                     i;
56 >
57 >        if (!curRF || strcmp(rop->preop.csym, curRF)) {
58 >                FILE    *fp = fopen(rop->preop.csym, "rb");
59 >                if (!rmx_load_header(&refm, fp)) {
60 >                        fprintf(stderr, "%s: cannot read info header\n",
61 >                                        rop->preop.csym);
62 >                        curRF = NULL;
63 >                        if (fp) fclose(fp);
64 >                        return(-1);
65 >                }
66 >                fclose(fp);
67 >                curRF = rop->preop.csym;
68 >        }
69 >        if (refm.ncomp == 3) {
70 >                rop->preop.csym = (refm.dtype == DTxyze) ? "XYZ" : "RGB";
71 >                return(checksymbolic(rop));
72 >        }
73 >        if (refm.ncomp == 2) {
74 >                fprintf(stderr, "%s: cannot convert to 2 components\n",
75 >                                curRF);
76                  return(-1);
77          }
78 <        return(1);
78 >        if (refm.ncomp == 1) {
79 >                rop->preop.csym = "Y";          /* XXX big assumption */
80 >                return(checksymbolic(rop));
81 >        }
82 >        if (refm.ncomp == nc &&
83 >                        !memcmp(refm.wlpart, rop->mtx->wlpart, sizeof(refm.wlpart)))
84 >                return(0);                      /* nothing to do */
85 >
86 >        if ((nc <= 3) | (nc > MAXCSAMP) | (refm.ncomp > MAXCSAMP)) {
87 >                fprintf(stderr, "%s: cannot resample from %d to %d components\n",
88 >                                curRF, nc, refm.ncomp);
89 >                return(-1);
90 >        }
91 >        rop->preop.clen = refm.ncomp * nc;      /* compute spec to ref */
92 >
93 >        for (i = 0; i < nc; i++) {
94 >                SCOLOR  scstim, scresp;
95 >                int     j;
96 >                memset(scstim, 0, sizeof(COLORV)*nc);
97 >                scstim[i] = 1.f;
98 >                convertscolor(scresp, refm.ncomp, refm.wlpart[0], refm.wlpart[3],
99 >                                scstim, nc, rop->mtx->wlpart[0], rop->mtx->wlpart[3]);
100 >                for (j = refm.ncomp; j-- > 0; )
101 >                        rop->preop.cmat[j*nc + i] = scresp[j];
102 >        }
103 >        memcpy(rop->mtx->wlpart, refm.wlpart, sizeof(rop->mtx->wlpart));
104 >        return(0);
105   }
106  
107 + /* Compute conversion row from spectrum to one channel of RGB */
108 + void
109 + rgbrow(ROPMAT *rop, int r, int p)
110 + {
111 +        const int       nc = rop->mtx->ncomp;
112 +        const float *   wlp = rop->mtx->wlpart;
113 +        int             i;
114 +
115 +        for (i = nc; i--; ) {
116 +                int     nmEnd = wlp[0] + (wlp[3] - wlp[0])*i/nc;
117 +                int     nmStart = wlp[0] + (wlp[3] - wlp[0])*(i+1)/nc;
118 +                COLOR   crgb;
119 +                spec_rgb(crgb, nmStart, nmEnd);
120 +                rop->preop.cmat[r*nc+i] = crgb[p];
121 +        }
122 + }
123 +
124 + /* Compute conversion row from spectrum to one channel of XYZ */
125 + void
126 + xyzrow(ROPMAT *rop, int r, int p)
127 + {
128 +        const int       nc = rop->mtx->ncomp;
129 +        const float *   wlp = rop->mtx->wlpart;
130 +        int             i;
131 +
132 +        for (i = nc; i--; ) {
133 +                int     nmEnd = wlp[0] + (wlp[3] - wlp[0])*i/nc;
134 +                int     nmStart = wlp[0] + (wlp[3] - wlp[0])*(i+1)/nc;
135 +                COLOR   cxyz;
136 +                spec_cie(cxyz, nmStart, nmEnd);
137 +                rop->preop.cmat[r*nc+i] = cxyz[p];
138 +        }
139 + }
140 +
141 + /* Use the spectral sensitivity function to compute matrix coefficients */
142 + void
143 + sensrow(ROPMAT *rop, int r, double (*sf)(const SCOLOR sc, int ncs, const float wlpt[4]))
144 + {
145 +        const int       nc = rop->mtx->ncomp;
146 +        int             i;
147 +
148 +        for (i = nc; i--; ) {
149 +                SCOLOR  sclr;
150 +                memset(sclr, 0, sizeof(COLORV)*nc);
151 +                sclr[i] = 1.f;
152 +                rop->preop.cmat[r*nc+i] = (*sf)(sclr, nc, rop->mtx->wlpart);
153 +        }
154 + }
155 +
156 + /* Check/set symbolic transform */
157 + int
158 + checksymbolic(ROPMAT *rop)
159 + {
160 +        const int       nc = rop->mtx->ncomp;
161 +        const int       dt = rop->mtx->dtype;
162 +        double          cf = 1;
163 +        int             i, j;
164 +                                        /* check suffix => reference file */
165 +        if (strchr(rop->preop.csym, '.') > rop->preop.csym)
166 +                return(checkreffile(rop));
167 +
168 +        if (nc < 3) {
169 +                fprintf(stderr, "%s: -c '%s' requires at least 3 components\n",
170 +                                rop->inspec, rop->preop.csym);
171 +                return(-1);
172 +        }
173 +        rop->preop.clen = strlen(rop->preop.csym) * nc;
174 +        if (rop->preop.clen > MAXCOMP*MAXCOMP) {
175 +                fprintf(stderr, "%s: -c '%s' results in too many components\n",
176 +                                rop->inspec, rop->preop.csym);
177 +                return(-1);
178 +        }
179 +        for (j = 0; rop->preop.csym[j]; j++) {
180 +                int     comp = 0;
181 +                switch (rop->preop.csym[j]) {
182 +                case 'B':
183 +                case 'b':
184 +                        ++comp;
185 +                        /* fall through */
186 +                case 'G':
187 +                case 'g':
188 +                        ++comp;
189 +                        /* fall through */
190 +                case 'R':
191 +                case 'r':
192 +                        if (rop->preop.csym[j] <= 'Z')
193 +                                cf = 1./WHTEFFICACY;
194 +                        if (dt == DTxyze) {
195 +                                for (i = 3; i--; )
196 +                                        rop->preop.cmat[j*nc+i] = cf*xyz2rgbmat[comp][i];
197 +                        } else if (nc == 3)
198 +                                rop->preop.cmat[j*nc+comp] = 1.;
199 +                        else
200 +                                rgbrow(rop, j, comp);
201 +                        break;
202 +                case 'Z':
203 +                case 'z':
204 +                        ++comp;
205 +                        /* fall through */
206 +                case 'Y':
207 +                case 'y':
208 +                        ++comp;
209 +                        /* fall through */
210 +                case 'X':
211 +                case 'x':
212 +                        if ((rop->preop.csym[j] <= 'Z') & (dt != DTxyze))
213 +                                cf = WHTEFFICACY;
214 +                        if (dt == DTxyze) {
215 +                                rop->preop.cmat[j*nc+comp] = 1.;
216 +                        } else if (nc == 3) {
217 +                                for (i = 3; i--; )
218 +                                        rop->preop.cmat[j*nc+i] =
219 +                                                        rgb2xyzmat[comp][i];
220 +                        } else if (comp == CIEY)
221 +                                sensrow(rop, j, scolor2photopic);
222 +                        else
223 +                                xyzrow(rop, j, comp);
224 +
225 +                        for (i = nc*(cf != 1); i--; )
226 +                                rop->preop.cmat[j*nc+i] *= cf;
227 +                        break;
228 +                case 'S':               /* scotopic (il)luminance */
229 +                        cf = WHTSCOTOPIC;
230 +                        /* fall through */
231 +                case 's':
232 +                        sensrow(rop, j, scolor2scotopic);
233 +                        for (i = nc*(cf != 1); i--; )
234 +                                rop->preop.cmat[j*nc+i] *= cf;
235 +                        break;
236 +                case 'M':               /* melanopic (il)luminance */
237 +                        cf = WHTMELANOPIC;
238 +                        /* fall through */
239 +                case 'm':
240 +                        sensrow(rop, j, scolor2melanopic);
241 +                        for (i = nc*(cf != 1); i--; )
242 +                                rop->preop.cmat[j*nc+i] *= cf;
243 +                        break;
244 +                case 'A':               /* average component */
245 +                case 'a':
246 +                        for (i = nc; i--; )
247 +                                rop->preop.cmat[j*nc+i] = 1./(double)nc;
248 +                        break;
249 +                default:
250 +                        fprintf(stderr, "%s: -c '%c' unsupported\n",
251 +                                rop->inspec, rop->preop.csym[j]);
252 +                        return(-1);
253 +                }
254 +        }
255 +                                        /* return recommended output type */
256 +        if (!strcasecmp(rop->preop.csym, "XYZ")) {
257 +                if (dt <= DTspec)
258 +                        return(DTxyze);
259 +        } else if (!strcasecmp(rop->preop.csym, "RGB")) {
260 +                if (dt <= DTspec)
261 +                        return(DTrgbe);
262 +        } else if (dt == DTspec)
263 +                return(DTfloat);        /* probably not actual spectrum */
264 +        return(0);
265 + }
266 +
267   /* Get matrix and perform unary operations */
268 < static RMATRIX *
268 > RMATRIX *
269   loadop(ROPMAT *rop)
270   {
271 +        int     outtype = 0;
272          RMATRIX *mres;
273 <        int     i;
273 >        int     i, j;
274  
275          if (loadmatrix(rop) < 0)                /* make sure we're loaded */
276                  return(NULL);
277  
278 <        if (rop->preop.nsf > 0) {               /* apply scalar(s) */
279 <                if (rop->preop.clen > 0) {
280 <                        fputs("Options -s and -c are exclusive\n", stderr);
278 >        if (rop->preop.csym &&                  /* symbolic transform? */
279 >                        (outtype = checksymbolic(rop)) < 0)
280 >                goto failure;
281 >        if (rop->preop.clen > 0) {              /* apply component transform? */
282 >                if (rop->preop.clen % rop->mtx->ncomp) {
283 >                        fprintf(stderr, "%s: -c must have N x %d coefficients\n",
284 >                                        rop->inspec, rop->mtx->ncomp);
285                          goto failure;
286                  }
287 +                if (rop->preop.nsf > 0) {       /* scale transform, first */
288 +                        if (rop->preop.nsf == 1) {
289 +                                for (i = rop->preop.clen; i--; )
290 +                                        rop->preop.cmat[i] *= rop->preop.sca[0];
291 +                        } else if (rop->preop.nsf*rop->mtx->ncomp != rop->preop.clen) {
292 +                                fprintf(stderr, "%s: -s must have one or %d factors\n",
293 +                                                rop->inspec,
294 +                                                rop->preop.clen/rop->mtx->ncomp);
295 +                                goto failure;
296 +                        } else {
297 +                                for (i = rop->preop.nsf; i--; )
298 +                                        for (j = rop->mtx->ncomp; j--; )
299 +                                                rop->preop.cmat[i*rop->mtx->ncomp+j]
300 +                                                                *= rop->preop.sca[i];
301 +                        }
302 +                }
303 +                mres = rmx_transform(rop->mtx, rop->preop.clen/rop->mtx->ncomp,
304 +                                        rop->preop.cmat);
305 +                if (mres == NULL) {
306 +                        fprintf(stderr, "%s: matrix transform failed\n",
307 +                                                rop->inspec);
308 +                        goto failure;
309 +                }
310 +                if (verbose)
311 +                        fprintf(stderr, "%s: applied %d x %d transform%s\n",
312 +                                        rop->inspec, mres->ncomp,
313 +                                        rop->mtx->ncomp,
314 +                                        rop->preop.nsf ? " (* scalar)" : "");
315 +                rop->preop.nsf = 0;             /* now folded in */
316 +                if ((mres->ncomp > 3) & (mres->dtype <= DTspec))
317 +                        outtype = DTfloat;      /* probably not actual spectrum */
318 +                rmx_free(rop->mtx);
319 +                rop->mtx = mres;
320 +        }
321 +        if (rop->preop.nsf > 0) {               /* apply scalar(s)? */
322                  if (rop->preop.nsf == 1) {
323                          for (i = rop->mtx->ncomp; --i; )
324                                  rop->preop.sca[i] = rop->preop.sca[0];
# Line 86 | Line 340 | loadop(ROPMAT *rop)
340                          fputs(" )\n", stderr);
341                  }
342          }
343 <        if (rop->preop.clen > 0) {      /* apply transform */
344 <                if (rop->preop.clen % rop->mtx->ncomp) {
91 <                        fprintf(stderr, "%s: -c must have N x %d coefficients\n",
92 <                                        rop->inspec, rop->mtx->ncomp);
93 <                        goto failure;
94 <                }
95 <                mres = rmx_transform(rop->mtx, rop->preop.clen/rop->mtx->ncomp,
96 <                                        rop->preop.cmat);
97 <                if (mres == NULL) {
98 <                        fprintf(stderr, "%s: matrix transform failed\n",
99 <                                                rop->inspec);
100 <                        goto failure;
101 <                }
102 <                if (verbose)
103 <                        fprintf(stderr, "%s: applied %d x %d transform\n",
104 <                                        rop->inspec, mres->ncomp,
105 <                                        rop->mtx->ncomp);
106 <                rmx_free(rop->mtx);
107 <                rop->mtx = mres;
108 <        }
109 <        if (rop->preop.transpose) {     /* transpose matrix? */
110 <                mres = rmx_transpose(rop->mtx);
111 <                if (mres == NULL) {
343 >        if (rop->preop.transpose) {             /* transpose matrix? */
344 >                if (!rmx_transpose(rop->mtx)) {
345                          fputs(rop->inspec, stderr);
346                          fputs(": transpose failed\n", stderr);
347                          goto failure;
# Line 117 | Line 350 | loadop(ROPMAT *rop)
350                          fputs(rop->inspec, stderr);
351                          fputs(": transposed rows and columns\n", stderr);
352                  }
120 –                rmx_free(rop->mtx);
121 –                rop->mtx = mres;
353          }
354          mres = rop->mtx;
355          rop->mtx = NULL;
356 +        if (outtype)
357 +                mres->dtype = outtype;
358          return(mres);
359   failure:
360          rmx_free(rop->mtx);
# Line 129 | Line 362 | failure:
362   }
363  
364   /* Execute binary operation, free matrix arguments and return new result */
365 < static RMATRIX *
365 > RMATRIX *
366   binaryop(const char *inspec, RMATRIX *mleft, int op, RMATRIX *mright)
367   {
368          RMATRIX *mres = NULL;
# Line 206 | Line 439 | binaryop(const char *inspec, RMATRIX *mleft, int op, R
439   }
440  
441   /* Perform matrix operations from left to right */
442 < static RMATRIX *
442 > RMATRIX *
443   op_left2right(ROPMAT *mop)
444   {
445          RMATRIX *mleft = loadop(mop);
# Line 222 | Line 455 | op_left2right(ROPMAT *mop)
455   }
456  
457   /* Perform matrix operations from right to left */
458 < static RMATRIX *
458 > RMATRIX *
459   op_right2left(ROPMAT *mop)
460   {
461          RMATRIX *mright;
# Line 251 | Line 484 | op_right2left(ROPMAT *mop)
484                                                  : (mop)->mtx->ncols)
485  
486   /* Should we prefer concatenating from rightmost matrix towards left? */
487 < static int
487 > int
488   prefer_right2left(ROPMAT *mop)
489   {
490          int     mri = 0;
# Line 278 | Line 511 | prefer_right2left(ROPMAT *mop)
511          return(t_ncols(mop+mri) < t_nrows(mop));
512   }
513  
514 < static int
514 > int
515   get_factors(double da[], int n, char *av[])
516   {
517          int     ac;
# Line 288 | Line 521 | get_factors(double da[], int n, char *av[])
521          return(ac);
522   }
523  
524 < static ROPMAT *
525 < grow_moparray(ROPMAT *mop, int n2alloc)
524 > ROPMAT *
525 > resize_moparr(ROPMAT *mop, int n2alloc)
526   {
527          int     nmats = 0;
528 +        int     i;
529  
530          while (mop[nmats++].binop)
531                  ;
532 +        for (i = nmats; i >= n2alloc; i--)
533 +                rmx_free(mop[i].mtx);
534          mop = (ROPMAT *)realloc(mop, n2alloc*sizeof(ROPMAT));
535          if (mop == NULL) {
536 <                fputs("Out of memory in grow_moparray()\n", stderr);
536 >                fputs("Out of memory in resize_moparr()\n", stderr);
537                  exit(1);
538          }
539          if (n2alloc > nmats)
# Line 309 | Line 545 | grow_moparray(ROPMAT *mop, int n2alloc)
545   int
546   main(int argc, char *argv[])
547   {
548 <        int     outfmt = DTfromHeader;
549 <        int     nall = 2;
550 <        ROPMAT  *mop = (ROPMAT *)calloc(nall, sizeof(ROPMAT));
551 <        int     nmats = 0;
552 <        RMATRIX *mres = NULL;
553 <        int     stdin_used = 0;
554 <        int     i;
548 >        int             outfmt = DTfromHeader;
549 >        const char      *defCsym = NULL;
550 >        int             nall = 2;
551 >        ROPMAT          *mop = (ROPMAT *)calloc(nall, sizeof(ROPMAT));
552 >        int             nmats = 0;
553 >        RMATRIX         *mres = NULL;
554 >        int             stdin_used = 0;
555 >        int             i;
556                                          /* get options and arguments */
557          for (i = 1; i < argc; i++) {
558                  if (argv[i][0] && !argv[i][1] &&
# Line 338 | Line 575 | main(int argc, char *argv[])
575                                  mop[nmats].inspec = stdin_name;
576                          } else
577                                  mop[nmats].inspec = argv[i];
578 +                        if (!mop[nmats].preop.csym)
579 +                                mop[nmats].preop.csym = defCsym;
580                          if (nmats > 0 && !mop[nmats-1].binop)
581                                  mop[nmats-1].binop = '.';
582                          nmats++;
# Line 373 | Line 612 | main(int argc, char *argv[])
612                                  i += mop[nmats].preop.nsf =
613                                          get_factors(mop[nmats].preop.sca,
614                                                          n, argv+i+1);
615 +                                if (mop[nmats].preop.nsf <= 0) {
616 +                                        fprintf(stderr, "%s: -s missing arguments\n",
617 +                                                        argv[0]);
618 +                                        goto userr;
619 +                                }
620                                  break;
621 +                        case 'C':
622 +                                if (!n || isflt(argv[i+1]))
623 +                                        goto userr;
624 +                                defCsym = mop[nmats].preop.csym = argv[++i];
625 +                                mop[nmats].preop.clen = 0;
626 +                                break;
627                          case 'c':
628 +                                if (n && !isflt(argv[i+1])) {
629 +                                        mop[nmats].preop.csym = argv[++i];
630 +                                        mop[nmats].preop.clen = 0;
631 +                                        break;
632 +                                }
633                                  if (n > MAXCOMP*MAXCOMP) n = MAXCOMP*MAXCOMP;
634                                  i += mop[nmats].preop.clen =
635                                          get_factors(mop[nmats].preop.cmat,
636                                                          n, argv+i+1);
637 +                                if (mop[nmats].preop.clen <= 0) {
638 +                                        fprintf(stderr, "%s: -c missing arguments\n",
639 +                                                        argv[0]);
640 +                                        goto userr;
641 +                                }
642 +                                mop[nmats].preop.csym = NULL;
643                                  break;
644                          case 'r':
645                                  if (argv[i][2] == 'f')
# Line 395 | Line 656 | main(int argc, char *argv[])
656                          }
657                  }
658                  if (nmats >= nall)
659 <                        mop = grow_moparray(mop, nall += 2);
659 >                        mop = resize_moparr(mop, nall += 2);
660          }
661          if (mop[0].inspec == NULL)      /* nothing to do? */
662                  goto userr;
# Line 415 | Line 676 | main(int argc, char *argv[])
676          mres = loadop(mop+nmats);
677          if (mres == NULL)
678                  return(1);
679 <                                        /* write result to stdout */
419 <        if (outfmt == DTfromHeader)
679 >        if (outfmt == DTfromHeader)     /* check data type */
680                  outfmt = mres->dtype;
681 <        if ((outfmt == DTrgbe) & (mres->ncomp > 3))
682 <                outfmt = DTspec;
683 <        if (outfmt != DTascii)
684 <                SET_FILE_BINARY(stdout);
685 <        newheader("RADIANCE", stdout);
426 <        printargs(argc, argv, stdout);
427 <        if (!rmx_write(mres, outfmt, stdout)) {
428 <                fprintf(stderr, "%s: error writing result matrix\n", argv[0]);
429 <                return(1);
681 >        if (outfmt == DTrgbe) {
682 >                if (mres->ncomp > 3)
683 >                        outfmt = DTspec;
684 >                else if (mres->dtype == DTxyze)
685 >                        outfmt = DTxyze;
686          }
687 <        /* rmx_free(mres); free(mop); */
688 <        return(0);
687 > #if DTrmx_native==DTfloat
688 >        if (outfmt == DTdouble)
689 >                fprintf(stderr,
690 >                        "%s: warning - writing float result as double\n",
691 >                                argv[0]);
692 > #endif
693 >        newheader("RADIANCE", stdout);  /* write result to stdout */
694 >        printargs(argc, argv, stdout);
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