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.18 by greg, Tue Jan 19 23:32:00 2021 UTC vs.
Revision 2.23 by greg, Tue Nov 28 16:36:50 2023 UTC

# Line 7 | Line 7 | static const char RCSid[] = "$Id$";
7  
8   #include <stdlib.h>
9   #include <errno.h>
10 + #include <ctype.h>
11   #include "rtio.h"
12   #include "resolu.h"
13   #include "rmatrix.h"
14   #include "platform.h"
15  
16 < #define MAXCOMP         16              /* #components we support */
16 > #define MAXCOMP         MAXCSAMP        /* #components we support */
17  
18 < /* unary matrix operation(s) */
18 > /* Unary matrix operation(s) */
19   typedef struct {
20          double          sca[MAXCOMP];           /* scalar coefficients */
21          double          cmat[MAXCOMP*MAXCOMP];  /* component transformation */
22          short           nsf;                    /* number of scalars */
23          short           clen;                   /* number of coefficients */
24 <        short           transpose;              /* do transpose? */
24 >        char            csym[11];               /* symbolic coefficients */
25 >        char            transpose;              /* do transpose? */
26   } RUNARYOP;
27  
28 < /* matrix input source and requested operation(s) */
28 > /* Matrix input source and requested operation(s) */
29   typedef struct {
30          const char      *inspec;                /* input specification */
31          RMPref          rmp;                    /* matrix preference */
# Line 38 | Line 40 | int    verbose = 0;                    /* verbose reporting? */
40   static int
41   loadmatrix(ROPMAT *rop)
42   {
43 <        if (rop->mtx != NULL)
43 >        if (rop->mtx != NULL)           /* already loaded? */
44                  return(0);
45  
46          rop->mtx = rmx_load(rop->inspec, rop->rmp);
47 <        if (rop->mtx == NULL) {
48 <                fputs(rop->inspec, stderr);
49 <                fputs(": cannot load matrix\n", stderr);
47 >
48 >        return(!rop->mtx ? -1 : 1);
49 > }
50 >
51 > /* Compute conversion row from spectrum to one channel of RGB */
52 > static void
53 > rgbrow(ROPMAT *rop, int r, int p)
54 > {
55 >        const int       nc = rop->mtx->ncomp;
56 >        const float *   wlp = rop->mtx->wlpart;
57 >        int             i;
58 >
59 >        for (i = nc; i--; ) {
60 >                int     nmEnd = wlp[0] + (wlp[3] - wlp[0])*i/nc;
61 >                int     nmStart = wlp[0] + (wlp[3] - wlp[0])*(i+1)/nc;
62 >                COLOR   crgb;
63 >                spec_rgb(crgb, nmStart, nmEnd);
64 >                rop->preop.cmat[r*nc+i] = crgb[p];
65 >        }
66 > }
67 >
68 > /* Compute conversion row from spectrum to one channel of XYZ */
69 > static void
70 > xyzrow(ROPMAT *rop, int r, int p)
71 > {
72 >        const int       nc = rop->mtx->ncomp;
73 >        const float *   wlp = rop->mtx->wlpart;
74 >        int             i;
75 >
76 >        for (i = nc; i--; ) {
77 >                int     nmEnd = wlp[0] + (wlp[3] - wlp[0])*i/nc;
78 >                int     nmStart = wlp[0] + (wlp[3] - wlp[0])*(i+1)/nc;
79 >                COLOR   cxyz;
80 >                spec_cie(cxyz, nmStart, nmEnd);
81 >                rop->preop.cmat[r*nc+i] = cxyz[p];
82 >        }
83 > }
84 >
85 > /* Use the spectral sensitivity function to compute matrix coefficients */
86 > static void
87 > sensrow(ROPMAT *rop, int r, double (*sf)(SCOLOR sc, int ncs, const float wlpt[4]))
88 > {
89 >        const int       nc = rop->mtx->ncomp;
90 >        int             i;
91 >
92 >        for (i = nc; i--; ) {
93 >                SCOLOR  sclr;
94 >                scolorblack(sclr);
95 >                sclr[i] = 1;
96 >                rop->preop.cmat[r*nc+i] = (*sf)(sclr, nc, rop->mtx->wlpart);
97 >        }
98 > }
99 >
100 > /* Check/set symbolic transform */
101 > static int
102 > checksymbolic(ROPMAT *rop)
103 > {
104 >        const int       nc = rop->mtx->ncomp;
105 >        const int       dt = rop->mtx->dtype;
106 >        int             i, j;
107 >
108 >        if (nc < 3) {
109 >                fprintf(stderr, "%s: -c '%s' requires at least 3 components\n",
110 >                                rop->inspec, rop->preop.csym);
111                  return(-1);
112          }
113 <        return(1);
113 >        rop->preop.clen = strlen(rop->preop.csym) * nc;
114 >        if (rop->preop.clen > MAXCOMP*MAXCOMP) {
115 >                fprintf(stderr, "%s: -c '%s' results in too many components\n",
116 >                                rop->inspec, rop->preop.csym);
117 >                return(-1);
118 >        }
119 >        for (j = 0; rop->preop.csym[j]; j++) {
120 >                int     comp = 0;
121 >                switch (rop->preop.csym[j]) {
122 >                case 'B':
123 >                        ++comp;
124 >                        /* fall through */
125 >                case 'G':
126 >                        ++comp;
127 >                        /* fall through */
128 >                case 'R':
129 >                        if (dt == DTxyze) {
130 >                                for (i = 3; i--; )
131 >                                        rop->preop.cmat[j*nc+i] = 1./WHTEFFICACY *
132 >                                                        xyz2rgbmat[comp][i];
133 >                        } else if (nc == 3)
134 >                                rop->preop.cmat[j*nc+comp] = 1.;
135 >                        else
136 >                                rgbrow(rop, j, comp);
137 >                        break;
138 >                case 'Z':
139 >                        ++comp;
140 >                        /* fall through */
141 >                case 'Y':
142 >                        ++comp;
143 >                        /* fall through */
144 >                case 'X':
145 >                        if (dt == DTxyze) {
146 >                                rop->preop.cmat[j*nc+comp] = 1.;
147 >                        } else if (nc == 3) {
148 >                                for (i = 3; i--; )
149 >                                        rop->preop.cmat[j*nc+i] =
150 >                                                        rgb2xyzmat[comp][i];
151 >                        } else if (comp == CIEY)
152 >                                sensrow(rop, j, scolor2photopic);
153 >                        else
154 >                                xyzrow(rop, j, comp);
155 >
156 >                        for (i = nc*(dt != DTxyze); i--; )
157 >                                rop->preop.cmat[j*nc+i] *= WHTEFFICACY;
158 >                        break;
159 >                case 'S':
160 >                        sensrow(rop, j, scolor2scotopic);
161 >                        for (i = nc; i--; )
162 >                                rop->preop.cmat[j*nc+i] *= WHTSCOTOPIC;
163 >                        break;
164 >                case 'M':
165 >                        sensrow(rop, j, scolor2melanopic);
166 >                        for (i = nc; i--; )
167 >                                rop->preop.cmat[j*nc+i] *= WHTMELANOPIC;
168 >                        break;
169 >                default:
170 >                        fprintf(stderr, "%s: -c '%c' unsupported\n",
171 >                                rop->inspec, rop->preop.csym[j]);
172 >                        return(-1);
173 >                }
174 >        }
175 >                                        /* return recommended output type */
176 >        if (!strcmp(rop->preop.csym, "XYZ")) {
177 >                if (dt <= DTspec)
178 >                        return(DTxyze);
179 >        } else if (!strcmp(rop->preop.csym, "RGB")) {
180 >                if (dt <= DTspec)
181 >                        return(DTrgbe);
182 >        }
183 >        if ((nc > 3) & (dt <= DTspec))
184 >                return(DTfloat);        /* probably not actual spectrum */
185 >        return(0);
186   }
187  
188   /* Get matrix and perform unary operations */
189   static RMATRIX *
190   loadop(ROPMAT *rop)
191   {
192 +        int     outtype = 0;
193          RMATRIX *mres;
194 <        int     i;
194 >        int     i, j;
195  
196          if (loadmatrix(rop) < 0)                /* make sure we're loaded */
197                  return(NULL);
198  
199 <        if (rop->preop.nsf > 0) {               /* apply scalar(s) */
200 <                if (rop->preop.clen > 0) {
201 <                        fputs("Options -s and -c are exclusive\n", stderr);
199 >        if (rop->preop.csym[0] &&               /* symbolic transform? */
200 >                        (outtype = checksymbolic(rop)) < 0)
201 >                goto failure;
202 >        if (rop->preop.clen > 0) {              /* apply component transform? */
203 >                if (rop->preop.clen % rop->mtx->ncomp) {
204 >                        fprintf(stderr, "%s: -c must have N x %d coefficients\n",
205 >                                        rop->inspec, rop->mtx->ncomp);
206                          goto failure;
207                  }
208 +                if (rop->preop.nsf > 0) {       /* scale transform, first */
209 +                        if (rop->preop.nsf == 1) {
210 +                                for (i = rop->preop.clen; i--; )
211 +                                        rop->preop.cmat[i] *= rop->preop.sca[0];
212 +                        } else if (rop->preop.nsf != rop->mtx->ncomp) {
213 +                                fprintf(stderr, "%s: -s must have one or %d factors\n",
214 +                                                rop->inspec, rop->mtx->ncomp);
215 +                                goto failure;
216 +                        } else {
217 +                                for (j = rop->preop.clen/rop->preop.nsf; j--; )
218 +                                        for (i = rop->preop.nsf; i--; )
219 +                                                rop->preop.cmat[j*rop->preop.nsf+i] *=
220 +                                                                rop->preop.sca[i];
221 +                        }
222 +                }
223 +                mres = rmx_transform(rop->mtx, rop->preop.clen/rop->mtx->ncomp,
224 +                                        rop->preop.cmat);
225 +                if (mres == NULL) {
226 +                        fprintf(stderr, "%s: matrix transform failed\n",
227 +                                                rop->inspec);
228 +                        goto failure;
229 +                }
230 +                if (verbose)
231 +                        fprintf(stderr, "%s: applied %d x %d transform%s\n",
232 +                                        rop->inspec, mres->ncomp,
233 +                                        rop->mtx->ncomp,
234 +                                        rop->preop.nsf ? " (* scalar)" : "");
235 +                rop->preop.nsf = 0;
236 +                if ((mres->ncomp > 3) & (mres->dtype <= DTspec))
237 +                        outtype = DTfloat;      /* probably not actual spectrum */
238 +                rmx_free(rop->mtx);
239 +                rop->mtx = mres;
240 +        }
241 +        if (rop->preop.nsf > 0) {               /* apply scalar(s)? */
242                  if (rop->preop.nsf == 1) {
243                          for (i = rop->mtx->ncomp; --i; )
244                                  rop->preop.sca[i] = rop->preop.sca[0];
# Line 86 | Line 260 | loadop(ROPMAT *rop)
260                          fputs(" )\n", stderr);
261                  }
262          }
263 <        if (rop->preop.clen > 0) {      /* apply transform */
90 <                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? */
263 >        if (rop->preop.transpose) {             /* transpose matrix? */
264                  mres = rmx_transpose(rop->mtx);
265                  if (mres == NULL) {
266                          fputs(rop->inspec, stderr);
# Line 122 | Line 276 | loadop(ROPMAT *rop)
276          }
277          mres = rop->mtx;
278          rop->mtx = NULL;
279 +        if (outtype)
280 +                mres->dtype = outtype;
281          return(mres);
282   failure:
283          rmx_free(rop->mtx);
# Line 239 | Line 395 | op_right2left(ROPMAT *mop)
395          while (rpos-- > 0) {
396                  if (mright == NULL)
397                          break;
398 <                mright = binaryop(mop[rpos].inspec,
398 >                mright = binaryop(mop[rpos+1].inspec,
399                                  loadop(mop+rpos), mop[rpos].binop, mright);
400          }
401          return(mright);
# Line 373 | Line 529 | main(int argc, char *argv[])
529                                  i += mop[nmats].preop.nsf =
530                                          get_factors(mop[nmats].preop.sca,
531                                                          n, argv+i+1);
532 +                                if (mop[nmats].preop.nsf <= 0) {
533 +                                        fprintf(stderr, "%s: -s missing arguments\n",
534 +                                                        argv[0]);
535 +                                        goto userr;
536 +                                }
537                                  break;
538                          case 'c':
539 +                                if (n && isupper(argv[i+1][0])) {
540 +                                        strlcpy(mop[nmats].preop.csym,
541 +                                                argv[++i],
542 +                                                sizeof(mop[0].preop.csym));
543 +                                        mop[nmats].preop.clen = 0;
544 +                                        break;
545 +                                }
546                                  if (n > MAXCOMP*MAXCOMP) n = MAXCOMP*MAXCOMP;
547                                  i += mop[nmats].preop.clen =
548                                          get_factors(mop[nmats].preop.cmat,
549                                                          n, argv+i+1);
550 +                                if (mop[nmats].preop.clen <= 0) {
551 +                                        fprintf(stderr, "%s: -c missing arguments\n",
552 +                                                        argv[0]);
553 +                                        goto userr;
554 +                                }
555 +                                mop[nmats].preop.csym[0] = '\0';
556                                  break;
557                          case 'r':
558                                  if (argv[i][2] == 'f')
# Line 415 | Line 589 | main(int argc, char *argv[])
589          mres = loadop(mop+nmats);
590          if (mres == NULL)
591                  return(1);
592 <                                        /* write result to stdout */
419 <        if (outfmt == DTfromHeader)
592 >        if (outfmt == DTfromHeader)     /* check data type */
593                  outfmt = mres->dtype;
594 <        if (outfmt != DTascii)
594 >        if (outfmt == DTrgbe) {
595 >                if (mres->ncomp > 3)
596 >                        outfmt = DTspec;
597 >                else if (mres->dtype == DTxyze)
598 >                        outfmt = DTxyze;
599 >        }
600 >        if (outfmt != DTascii)          /* write result to stdout */
601                  SET_FILE_BINARY(stdout);
602          newheader("RADIANCE", stdout);
603          printargs(argc, argv, stdout);
604 <        if (!rmx_write(mres, outfmt, stdout)) {
605 <                fprintf(stderr, "%s: error writing result matrix\n", argv[0]);
427 <                return(1);
428 <        }
429 <        /* rmx_free(mres); free(mop); */
430 <        return(0);
604 >
605 >        return(rmx_write(mres, outfmt, stdout) ? 0 : 1);
606   userr:
607          fprintf(stderr,
608 <        "Usage: %s [-v][-f[adfc][-t][-s sf .. | -c ce ..][-r[fb]] m1 [.+*/] .. > mres\n",
608 >        "Usage: %s [-v][-f[adfc][-t][-s sf .. | -c ce ..][-rf|-rb] m1 [.+*/] .. > mres\n",
609                          argv[0]);
610          return(1);
611   }

Diff Legend

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