ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/src/util/rmtxop.c
(Generate patch)

Comparing src/util/rmtxop.c (file contents):
Revision 2.17 by greg, Sat Dec 28 18:05:14 2019 UTC vs.
Revision 2.27 by greg, Sat Dec 2 00:42:21 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>
9   #include "rtio.h"
10   #include "resolu.h"
11   #include "rmatrix.h"
12   #include "platform.h"
13  
14 < #define MAXCOMP         16              /* #components we support */
14 > #define MAXCOMP         MAXCSAMP        /* #components we support */
15  
16 < static const char       stdin_name[] = "<stdin>";
18 <
19 < /* unary matrix operation(s) */
16 > /* Unary matrix operation(s) */
17   typedef struct {
21 –        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 +        short           nsf;                    /* number of scalars */
23          short           transpose;              /* do transpose? */
24   } RUNARYOP;
25  
26 < /* matrix input source and requested operation(s) */
26 > /* Matrix input source and requested operation(s) */
27   typedef struct {
28          const char      *inspec;                /* input specification */
29 +        RMPref          rmp;                    /* matrix preference */
30          RUNARYOP        preop;                  /* unary operation(s) */
31          RMATRIX         *mtx;                   /* original matrix if loaded */
32          int             binop;                  /* binary op with next (or 0) */
# Line 39 | Line 38 | int    verbose = 0;                    /* verbose reporting? */
38   static int
39   loadmatrix(ROPMAT *rop)
40   {
41 <        if (rop->mtx != NULL)
41 >        if (rop->mtx != NULL)           /* already loaded? */
42                  return(0);
43  
44 <        rop->mtx = rmx_load(rop->inspec == stdin_name ?
45 <                                (const char *)NULL : rop->inspec);
46 <        if (rop->mtx == NULL) {
47 <                fputs(rop->inspec, stderr);
48 <                fputs(": cannot load matrix\n", stderr);
44 >        rop->mtx = rmx_load(rop->inspec, rop->rmp);
45 >
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) & (refm.dtype != DTspec)) {
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 <        return(1);
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 >        return(0);
107   }
108  
109 + /* Compute conversion row from spectrum to one channel of RGB */
110 + static void
111 + rgbrow(ROPMAT *rop, int r, int p)
112 + {
113 +        const int       nc = rop->mtx->ncomp;
114 +        const float *   wlp = rop->mtx->wlpart;
115 +        int             i;
116 +
117 +        for (i = nc; i--; ) {
118 +                int     nmEnd = wlp[0] + (wlp[3] - wlp[0])*i/nc;
119 +                int     nmStart = wlp[0] + (wlp[3] - wlp[0])*(i+1)/nc;
120 +                COLOR   crgb;
121 +                spec_rgb(crgb, nmStart, nmEnd);
122 +                rop->preop.cmat[r*nc+i] = crgb[p];
123 +        }
124 + }
125 +
126 + /* Compute conversion row from spectrum to one channel of XYZ */
127 + static void
128 + xyzrow(ROPMAT *rop, int r, int p)
129 + {
130 +        const int       nc = rop->mtx->ncomp;
131 +        const float *   wlp = rop->mtx->wlpart;
132 +        int             i;
133 +
134 +        for (i = nc; i--; ) {
135 +                int     nmEnd = wlp[0] + (wlp[3] - wlp[0])*i/nc;
136 +                int     nmStart = wlp[0] + (wlp[3] - wlp[0])*(i+1)/nc;
137 +                COLOR   cxyz;
138 +                spec_cie(cxyz, nmStart, nmEnd);
139 +                rop->preop.cmat[r*nc+i] = cxyz[p];
140 +        }
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]))
146 + {
147 +        const int       nc = rop->mtx->ncomp;
148 +        int             i;
149 +
150 +        for (i = nc; i--; ) {
151 +                SCOLOR  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
160 + checksymbolic(ROPMAT *rop)
161 + {
162 +        const int       nc = rop->mtx->ncomp;
163 +        const int       dt = rop->mtx->dtype;
164 +        int             i, j;
165 +                                        /* check suffix => reference file */
166 +        if (strchr(rop->preop.csym, '.') > rop->preop.csym)
167 +                return(checkreffile(rop));
168 +
169 +        if (nc < 3) {
170 +                fprintf(stderr, "%s: -c '%s' requires at least 3 components\n",
171 +                                rop->inspec, rop->preop.csym);
172 +                return(-1);
173 +        }
174 +        rop->preop.clen = strlen(rop->preop.csym) * nc;
175 +        if (rop->preop.clen > MAXCOMP*MAXCOMP) {
176 +                fprintf(stderr, "%s: -c '%s' results in too many components\n",
177 +                                rop->inspec, rop->preop.csym);
178 +                return(-1);
179 +        }
180 +        for (j = 0; rop->preop.csym[j]; j++) {
181 +                int     comp = 0;
182 +                switch (rop->preop.csym[j]) {
183 +                case 'B':
184 +                        ++comp;
185 +                        /* fall through */
186 +                case 'G':
187 +                        ++comp;
188 +                        /* fall through */
189 +                case 'R':
190 +                        if (dt == DTxyze) {
191 +                                for (i = 3; i--; )
192 +                                        rop->preop.cmat[j*nc+i] = 1./WHTEFFICACY *
193 +                                                        xyz2rgbmat[comp][i];
194 +                        } else if (nc == 3)
195 +                                rop->preop.cmat[j*nc+comp] = 1.;
196 +                        else
197 +                                rgbrow(rop, j, comp);
198 +                        break;
199 +                case 'Z':
200 +                        ++comp;
201 +                        /* fall through */
202 +                case 'Y':
203 +                        ++comp;
204 +                        /* fall through */
205 +                case 'X':
206 +                        if (dt == DTxyze) {
207 +                                rop->preop.cmat[j*nc+comp] = 1.;
208 +                        } else if (nc == 3) {
209 +                                for (i = 3; i--; )
210 +                                        rop->preop.cmat[j*nc+i] =
211 +                                                        rgb2xyzmat[comp][i];
212 +                        } else if (comp == CIEY)
213 +                                sensrow(rop, j, scolor2photopic);
214 +                        else
215 +                                xyzrow(rop, j, comp);
216 +
217 +                        for (i = nc*(dt != DTxyze); i--; )
218 +                                rop->preop.cmat[j*nc+i] *= WHTEFFICACY;
219 +                        break;
220 +                case 'S':               /* scotopic (il)luminance */
221 +                        sensrow(rop, j, scolor2scotopic);
222 +                        for (i = nc; i--; )
223 +                                rop->preop.cmat[j*nc+i] *= WHTSCOTOPIC;
224 +                        break;
225 +                case 'M':               /* melanopic (il)luminance */
226 +                        sensrow(rop, j, scolor2melanopic);
227 +                        for (i = nc; i--; )
228 +                                rop->preop.cmat[j*nc+i] *= WHTMELANOPIC;
229 +                        break;
230 +                case 'A':               /* average component */
231 +                        for (i = nc; i--; )
232 +                                rop->preop.cmat[j*nc+i] = 1./(double)nc;
233 +                        break;
234 +                default:
235 +                        fprintf(stderr, "%s: -c '%c' unsupported\n",
236 +                                rop->inspec, rop->preop.csym[j]);
237 +                        return(-1);
238 +                }
239 +        }
240 +                                        /* return recommended output type */
241 +        if (!strcmp(rop->preop.csym, "XYZ")) {
242 +                if (dt <= DTspec)
243 +                        return(DTxyze);
244 +        } else if (!strcmp(rop->preop.csym, "RGB")) {
245 +                if (dt <= DTspec)
246 +                        return(DTrgbe);
247 +        }
248 +        if ((nc > 3) & (dt <= DTspec))
249 +                return(DTfloat);        /* probably not actual spectrum */
250 +        return(0);
251 + }
252 +
253   /* Get matrix and perform unary operations */
254   static RMATRIX *
255   loadop(ROPMAT *rop)
256   {
257 +        int     outtype = 0;
258          RMATRIX *mres;
259 <        int     i;
259 >        int     i, j;
260  
261          if (loadmatrix(rop) < 0)                /* make sure we're loaded */
262                  return(NULL);
263  
264 <        if (rop->preop.nsf > 0) {               /* apply scalar(s) */
265 <                if (rop->preop.clen > 0) {
266 <                        fputs("Options -s and -c are exclusive\n", stderr);
264 >        if (rop->preop.csym &&                  /* symbolic transform? */
265 >                        (outtype = checksymbolic(rop)) < 0)
266 >                goto failure;
267 >        if (rop->preop.clen > 0) {              /* apply component transform? */
268 >                if (rop->preop.clen % rop->mtx->ncomp) {
269 >                        fprintf(stderr, "%s: -c must have N x %d coefficients\n",
270 >                                        rop->inspec, rop->mtx->ncomp);
271                          goto failure;
272                  }
273 +                if (rop->preop.nsf > 0) {       /* scale transform, first */
274 +                        if (rop->preop.nsf == 1) {
275 +                                for (i = rop->preop.clen; i--; )
276 +                                        rop->preop.cmat[i] *= rop->preop.sca[0];
277 +                        } else if (rop->preop.nsf != rop->mtx->ncomp) {
278 +                                fprintf(stderr, "%s: -s must have one or %d factors\n",
279 +                                                rop->inspec, rop->mtx->ncomp);
280 +                                goto failure;
281 +                        } else {
282 +                                for (j = rop->preop.clen/rop->preop.nsf; j--; )
283 +                                        for (i = rop->preop.nsf; i--; )
284 +                                                rop->preop.cmat[j*rop->preop.nsf+i] *=
285 +                                                                rop->preop.sca[i];
286 +                        }
287 +                }
288 +                mres = rmx_transform(rop->mtx, rop->preop.clen/rop->mtx->ncomp,
289 +                                        rop->preop.cmat);
290 +                if (mres == NULL) {
291 +                        fprintf(stderr, "%s: matrix transform failed\n",
292 +                                                rop->inspec);
293 +                        goto failure;
294 +                }
295 +                if (verbose)
296 +                        fprintf(stderr, "%s: applied %d x %d transform%s\n",
297 +                                        rop->inspec, mres->ncomp,
298 +                                        rop->mtx->ncomp,
299 +                                        rop->preop.nsf ? " (* scalar)" : "");
300 +                rop->preop.nsf = 0;             /* now folded in */
301 +                if ((mres->ncomp > 3) & (mres->dtype <= DTspec))
302 +                        outtype = DTfloat;      /* probably not actual spectrum */
303 +                rmx_free(rop->mtx);
304 +                rop->mtx = mres;
305 +        }
306 +        if (rop->preop.nsf > 0) {               /* apply scalar(s)? */
307                  if (rop->preop.nsf == 1) {
308                          for (i = rop->mtx->ncomp; --i; )
309                                  rop->preop.sca[i] = rop->preop.sca[0];
# Line 88 | Line 325 | loadop(ROPMAT *rop)
325                          fputs(" )\n", stderr);
326                  }
327          }
328 <        if (rop->preop.clen > 0) {      /* apply transform */
92 <                if (rop->preop.clen % rop->mtx->ncomp) {
93 <                        fprintf(stderr, "%s: -c must have N x %d coefficients\n",
94 <                                        rop->inspec, rop->mtx->ncomp);
95 <                        goto failure;
96 <                }
97 <                mres = rmx_transform(rop->mtx, rop->preop.clen/rop->mtx->ncomp,
98 <                                        rop->preop.cmat);
99 <                if (mres == NULL) {
100 <                        fprintf(stderr, "%s: matrix transform failed\n",
101 <                                                rop->inspec);
102 <                        goto failure;
103 <                }
104 <                if (verbose)
105 <                        fprintf(stderr, "%s: applied %d x %d transform\n",
106 <                                        rop->inspec, mres->ncomp,
107 <                                        rop->mtx->ncomp);
108 <                rmx_free(rop->mtx);
109 <                rop->mtx = mres;
110 <        }
111 <        if (rop->preop.transpose) {     /* transpose matrix? */
328 >        if (rop->preop.transpose) {             /* transpose matrix? */
329                  mres = rmx_transpose(rop->mtx);
330                  if (mres == NULL) {
331                          fputs(rop->inspec, stderr);
# Line 124 | Line 341 | loadop(ROPMAT *rop)
341          }
342          mres = rop->mtx;
343          rop->mtx = NULL;
344 +        if (outtype)
345 +                mres->dtype = outtype;
346          return(mres);
347   failure:
348          rmx_free(rop->mtx);
# Line 241 | Line 460 | op_right2left(ROPMAT *mop)
460          while (rpos-- > 0) {
461                  if (mright == NULL)
462                          break;
463 <                mright = binaryop(mop[rpos].inspec,
463 >                mright = binaryop(mop[rpos+1].inspec,
464                                  loadop(mop+rpos), mop[rpos].binop, mright);
465          }
466          return(mright);
# Line 291 | Line 510 | get_factors(double da[], int n, char *av[])
510   }
511  
512   static ROPMAT *
513 < grow_moparray(ROPMAT *mop, int n2alloc)
513 > resize_moparr(ROPMAT *mop, int n2alloc)
514   {
515          int     nmats = 0;
516 +        int     i;
517  
518          while (mop[nmats++].binop)
519                  ;
520 +        for (i = nmats; i > n2alloc; i--)
521 +                rmx_free(mop[i].mtx);
522          mop = (ROPMAT *)realloc(mop, n2alloc*sizeof(ROPMAT));
523          if (mop == NULL) {
524 <                fputs("Out of memory in grow_moparray()\n", stderr);
524 >                fputs("Out of memory in resize_moparr()\n", stderr);
525                  exit(1);
526          }
527          if (n2alloc > nmats)
# Line 311 | Line 533 | grow_moparray(ROPMAT *mop, int n2alloc)
533   int
534   main(int argc, char *argv[])
535   {
536 <        int     outfmt = DTfromHeader;
537 <        int     nall = 2;
538 <        ROPMAT  *mop = (ROPMAT *)calloc(nall, sizeof(ROPMAT));
539 <        int     nmats = 0;
540 <        RMATRIX *mres = NULL;
541 <        int     stdin_used = 0;
542 <        int     i;
536 >        int             outfmt = DTfromHeader;
537 >        const char      *defCsym = NULL;
538 >        int             nall = 2;
539 >        ROPMAT          *mop = (ROPMAT *)calloc(nall, sizeof(ROPMAT));
540 >        int             nmats = 0;
541 >        RMATRIX         *mres = NULL;
542 >        int             stdin_used = 0;
543 >        int             i;
544                                          /* get options and arguments */
545          for (i = 1; i < argc; i++) {
546                  if (argv[i][0] && !argv[i][1] &&
# Line 340 | Line 563 | main(int argc, char *argv[])
563                                  mop[nmats].inspec = stdin_name;
564                          } else
565                                  mop[nmats].inspec = argv[i];
566 +                        if (!mop[nmats].preop.csym)
567 +                                mop[nmats].preop.csym = defCsym;
568                          if (nmats > 0 && !mop[nmats-1].binop)
569                                  mop[nmats-1].binop = '.';
570                          nmats++;
# Line 375 | Line 600 | main(int argc, char *argv[])
600                                  i += mop[nmats].preop.nsf =
601                                          get_factors(mop[nmats].preop.sca,
602                                                          n, argv+i+1);
603 +                                if (mop[nmats].preop.nsf <= 0) {
604 +                                        fprintf(stderr, "%s: -s missing arguments\n",
605 +                                                        argv[0]);
606 +                                        goto userr;
607 +                                }
608                                  break;
609 +                        case 'C':
610 +                                if (!n || isflt(argv[i+1]))
611 +                                        goto userr;
612 +                                defCsym = mop[nmats].preop.csym = argv[++i];
613 +                                mop[nmats].preop.clen = 0;
614 +                                break;
615                          case 'c':
616 +                                if (n && !isflt(argv[i+1])) {
617 +                                        mop[nmats].preop.csym = argv[++i];
618 +                                        mop[nmats].preop.clen = 0;
619 +                                        break;
620 +                                }
621                                  if (n > MAXCOMP*MAXCOMP) n = MAXCOMP*MAXCOMP;
622                                  i += mop[nmats].preop.clen =
623                                          get_factors(mop[nmats].preop.cmat,
624                                                          n, argv+i+1);
625 +                                if (mop[nmats].preop.clen <= 0) {
626 +                                        fprintf(stderr, "%s: -c missing arguments\n",
627 +                                                        argv[0]);
628 +                                        goto userr;
629 +                                }
630 +                                mop[nmats].preop.csym = NULL;
631                                  break;
632 +                        case 'r':
633 +                                if (argv[i][2] == 'f')
634 +                                        mop[nmats].rmp = RMPreflF;
635 +                                else if (argv[i][2] == 'b')
636 +                                        mop[nmats].rmp = RMPreflB;
637 +                                else
638 +                                        goto userr;
639 +                                break;
640                          default:
641                                  fprintf(stderr, "%s: unknown operation '%s'\n",
642                                                  argv[0], argv[i]);
# Line 389 | Line 644 | main(int argc, char *argv[])
644                          }
645                  }
646                  if (nmats >= nall)
647 <                        mop = grow_moparray(mop, nall += 2);
647 >                        mop = resize_moparr(mop, nall += 2);
648          }
649          if (mop[0].inspec == NULL)      /* nothing to do? */
650                  goto userr;
# Line 409 | Line 664 | main(int argc, char *argv[])
664          mres = loadop(mop+nmats);
665          if (mres == NULL)
666                  return(1);
667 <                                        /* write result to stdout */
413 <        if (outfmt == DTfromHeader)
667 >        if (outfmt == DTfromHeader)     /* check data type */
668                  outfmt = mres->dtype;
669 <        if (outfmt != DTascii)
670 <                SET_FILE_BINARY(stdout);
671 <        newheader("RADIANCE", stdout);
672 <        printargs(argc, argv, stdout);
673 <        if (!rmx_write(mres, outfmt, stdout)) {
420 <                fprintf(stderr, "%s: error writing result matrix\n", argv[0]);
421 <                return(1);
669 >        if (outfmt == DTrgbe) {
670 >                if (mres->ncomp > 3)
671 >                        outfmt = DTspec;
672 >                else if (mres->dtype == DTxyze)
673 >                        outfmt = DTxyze;
674          }
675 <        /* rmx_free(mres); free(mop); */
676 <        return(0);
675 >        newheader("RADIANCE", stdout);  /* write result to stdout */
676 >        printargs(argc, argv, stdout);
677 >        return(rmx_write(mres, outfmt, stdout) ? 0 : 1);
678   userr:
679          fprintf(stderr,
680 <        "Usage: %s [-v][-f[adfc][-t][-s sf .. | -c ce ..] m1 [.+*/] .. > mres\n",
680 >        "Usage: %s [-v][-f{adfc}][-t][-s sf .. | -c ce ..][-rf|-rb] m1 [.+*/] .. > mres\n",
681                          argv[0]);
682          return(1);
683   }

Diff Legend

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