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.22 by greg, Mon Nov 27 22:04:45 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 50 | Line 52 | loadmatrix(ROPMAT *rop)
52          return(1);
53   }
54  
55 + /* Compute conversion row from spectrum to one channel of RGB */
56 + static void
57 + rgbrow(ROPMAT *rop, int r, int p)
58 + {
59 +        const int       nc = rop->mtx->ncomp;
60 +        const float *   wlp = rop->mtx->wlpart;
61 +        int             i;
62 +
63 +        for (i = nc; i--; ) {
64 +                int     nmEnd = wlp[0] + (wlp[3] - wlp[0])*i/nc;
65 +                int     nmStart = wlp[0] + (wlp[3] - wlp[0])*(i+1)/nc;
66 +                COLOR   crgb;
67 +                spec_rgb(crgb, nmStart, nmEnd);
68 +                rop->preop.cmat[r*nc+i] = crgb[p];
69 +        }
70 + }
71 +
72 + /* Compute conversion row from spectrum to one channel of XYZ */
73 + static void
74 + xyzrow(ROPMAT *rop, int r, int p)
75 + {
76 +        const int       nc = rop->mtx->ncomp;
77 +        const float *   wlp = rop->mtx->wlpart;
78 +        int             i;
79 +
80 +        for (i = nc; i--; ) {
81 +                int     nmEnd = wlp[0] + (wlp[3] - wlp[0])*i/nc;
82 +                int     nmStart = wlp[0] + (wlp[3] - wlp[0])*(i+1)/nc;
83 +                COLOR   cxyz;
84 +                spec_cie(cxyz, nmStart, nmEnd);
85 +                rop->preop.cmat[r*nc+i] = cxyz[p];
86 +        }
87 + }
88 +
89 + /* Use the spectral sensitivity function to compute matrix coefficients */
90 + static void
91 + sensrow(ROPMAT *rop, int r, double (*sf)(SCOLOR sc, int ncs, const float wlpt[4]))
92 + {
93 +        const int       nc = rop->mtx->ncomp;
94 +        int             i;
95 +
96 +        for (i = nc; i--; ) {
97 +                SCOLOR  sclr;
98 +                scolorblack(sclr);
99 +                sclr[i] = 1;
100 +                rop->preop.cmat[r*nc+i] = (*sf)(sclr, nc, rop->mtx->wlpart);
101 +        }
102 + }
103 +
104 + /* Check/set symbolic transform */
105 + static int
106 + checksymbolic(ROPMAT *rop)
107 + {
108 +        const int       nc = rop->mtx->ncomp;
109 +        const int       dt = rop->mtx->dtype;
110 +        int             i, j;
111 +
112 +        if (nc < 3) {
113 +                fprintf(stderr, "%s: -c '%s' requires at least 3 components\n",
114 +                                rop->inspec, rop->preop.csym);
115 +                return(-1);
116 +        }
117 +        rop->preop.clen = strlen(rop->preop.csym) * nc;
118 +        if (rop->preop.clen > MAXCOMP*MAXCOMP) {
119 +                fprintf(stderr, "%s: -c '%s' results in too many components\n",
120 +                                rop->inspec, rop->preop.csym);
121 +                return(-1);
122 +        }
123 +        for (j = 0; rop->preop.csym[j]; j++) {
124 +                int     comp = 0;
125 +                switch (rop->preop.csym[j]) {
126 +                case 'B':
127 +                        ++comp;
128 +                        /* fall through */
129 +                case 'G':
130 +                        ++comp;
131 +                        /* fall through */
132 +                case 'R':
133 +                        if (dt == DTxyze) {
134 +                                for (i = 3; i--; )
135 +                                        rop->preop.cmat[j*nc+i] = 1./WHTEFFICACY *
136 +                                                        xyz2rgbmat[comp][i];
137 +                        } else if (nc == 3)
138 +                                rop->preop.cmat[j*nc+comp] = 1.;
139 +                        else
140 +                                rgbrow(rop, j, comp);
141 +                        break;
142 +                case 'Z':
143 +                        ++comp;
144 +                        /* fall through */
145 +                case 'Y':
146 +                        ++comp;
147 +                        /* fall through */
148 +                case 'X':
149 +                        if (dt == DTxyze) {
150 +                                rop->preop.cmat[j*nc+comp] = 1.;
151 +                        } else if (nc == 3) {
152 +                                for (i = 3; i--; )
153 +                                        rop->preop.cmat[j*nc+i] =
154 +                                                        rgb2xyzmat[comp][i];
155 +                        } else if (comp == CIEY)
156 +                                sensrow(rop, j, scolor2photopic);
157 +                        else
158 +                                xyzrow(rop, j, comp);
159 +
160 +                        for (i = nc*(dt != DTxyze); i--; )
161 +                                rop->preop.cmat[j*nc+i] *= WHTEFFICACY;
162 +                        break;
163 +                case 'S':
164 +                        sensrow(rop, j, scolor2scotopic);
165 +                        for (i = nc; i--; )
166 +                                rop->preop.cmat[j*nc+i] *= WHTSCOTOPIC;
167 +                        break;
168 +                case 'M':
169 +                        sensrow(rop, j, scolor2melanopic);
170 +                        for (i = nc; i--; )
171 +                                rop->preop.cmat[j*nc+i] *= WHTMELANOPIC;
172 +                        break;
173 +                default:
174 +                        fprintf(stderr, "%s: -c '%c' unsupported\n",
175 +                                rop->inspec, rop->preop.csym[j]);
176 +                        return(-1);
177 +                }
178 +        }
179 +                                        /* return recommended output type */
180 +        if (!strcmp(rop->preop.csym, "XYZ")) {
181 +                if (dt <= DTspec)
182 +                        return(DTxyze);
183 +        } else if (!strcmp(rop->preop.csym, "RGB")) {
184 +                if (dt <= DTspec)
185 +                        return(DTrgbe);
186 +        }
187 +        if ((nc > 3) & (dt <= DTspec))
188 +                return(DTfloat);        /* probably not actual spectrum */
189 +        return(0);
190 + }
191 +
192   /* Get matrix and perform unary operations */
193   static RMATRIX *
194   loadop(ROPMAT *rop)
195   {
196 +        int     outtype = 0;
197          RMATRIX *mres;
198 <        int     i;
198 >        int     i, j;
199  
200          if (loadmatrix(rop) < 0)                /* make sure we're loaded */
201                  return(NULL);
202  
203 <        if (rop->preop.nsf > 0) {               /* apply scalar(s) */
204 <                if (rop->preop.clen > 0) {
205 <                        fputs("Options -s and -c are exclusive\n", stderr);
203 >        if (rop->preop.csym[0] &&               /* symbolic transform? */
204 >                        (outtype = checksymbolic(rop)) < 0)
205 >                goto failure;
206 >        if (rop->preop.clen > 0) {              /* apply component transform? */
207 >                if (rop->preop.clen % rop->mtx->ncomp) {
208 >                        fprintf(stderr, "%s: -c must have N x %d coefficients\n",
209 >                                        rop->inspec, rop->mtx->ncomp);
210                          goto failure;
211                  }
212 +                if (rop->preop.nsf > 0) {       /* scale transform, first */
213 +                        if (rop->preop.nsf == 1) {
214 +                                for (i = rop->preop.clen; i--; )
215 +                                        rop->preop.cmat[i] *= rop->preop.sca[0];
216 +                        } else if (rop->preop.nsf != rop->mtx->ncomp) {
217 +                                fprintf(stderr, "%s: -s must have one or %d factors\n",
218 +                                                rop->inspec, rop->mtx->ncomp);
219 +                                goto failure;
220 +                        } else {
221 +                                for (j = rop->preop.clen/rop->preop.nsf; j--; )
222 +                                        for (i = rop->preop.nsf; i--; )
223 +                                                rop->preop.cmat[j*rop->preop.nsf+i] *=
224 +                                                                rop->preop.sca[i];
225 +                        }
226 +                }
227 +                mres = rmx_transform(rop->mtx, rop->preop.clen/rop->mtx->ncomp,
228 +                                        rop->preop.cmat);
229 +                if (mres == NULL) {
230 +                        fprintf(stderr, "%s: matrix transform failed\n",
231 +                                                rop->inspec);
232 +                        goto failure;
233 +                }
234 +                if (verbose)
235 +                        fprintf(stderr, "%s: applied %d x %d transform%s\n",
236 +                                        rop->inspec, mres->ncomp,
237 +                                        rop->mtx->ncomp,
238 +                                        rop->preop.nsf ? " (* scalar)" : "");
239 +                rop->preop.nsf = 0;
240 +                if ((mres->ncomp > 3) & (mres->dtype <= DTspec))
241 +                        outtype = DTfloat;      /* probably not actual spectrum */
242 +                rmx_free(rop->mtx);
243 +                rop->mtx = mres;
244 +        }
245 +        if (rop->preop.nsf > 0) {               /* apply scalar(s)? */
246                  if (rop->preop.nsf == 1) {
247                          for (i = rop->mtx->ncomp; --i; )
248                                  rop->preop.sca[i] = rop->preop.sca[0];
# Line 86 | Line 264 | loadop(ROPMAT *rop)
264                          fputs(" )\n", stderr);
265                  }
266          }
267 <        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? */
267 >        if (rop->preop.transpose) {             /* transpose matrix? */
268                  mres = rmx_transpose(rop->mtx);
269                  if (mres == NULL) {
270                          fputs(rop->inspec, stderr);
# Line 122 | Line 280 | loadop(ROPMAT *rop)
280          }
281          mres = rop->mtx;
282          rop->mtx = NULL;
283 +        if (outtype)
284 +                mres->dtype = outtype;
285          return(mres);
286   failure:
287          rmx_free(rop->mtx);
# Line 373 | Line 533 | main(int argc, char *argv[])
533                                  i += mop[nmats].preop.nsf =
534                                          get_factors(mop[nmats].preop.sca,
535                                                          n, argv+i+1);
536 +                                if (mop[nmats].preop.nsf <= 0) {
537 +                                        fprintf(stderr, "%s: -s missing arguments\n",
538 +                                                        argv[0]);
539 +                                        goto userr;
540 +                                }
541                                  break;
542                          case 'c':
543 +                                if (n && isupper(argv[i+1][0])) {
544 +                                        strlcpy(mop[nmats].preop.csym,
545 +                                                argv[++i],
546 +                                                sizeof(mop[0].preop.csym));
547 +                                        mop[nmats].preop.clen = 0;
548 +                                        break;
549 +                                }
550                                  if (n > MAXCOMP*MAXCOMP) n = MAXCOMP*MAXCOMP;
551                                  i += mop[nmats].preop.clen =
552                                          get_factors(mop[nmats].preop.cmat,
553                                                          n, argv+i+1);
554 +                                if (mop[nmats].preop.clen <= 0) {
555 +                                        fprintf(stderr, "%s: -c missing arguments\n",
556 +                                                        argv[0]);
557 +                                        goto userr;
558 +                                }
559 +                                mop[nmats].preop.csym[0] = '\0';
560                                  break;
561                          case 'r':
562                                  if (argv[i][2] == 'f')
# Line 415 | Line 593 | main(int argc, char *argv[])
593          mres = loadop(mop+nmats);
594          if (mres == NULL)
595                  return(1);
596 <                                        /* write result to stdout */
419 <        if (outfmt == DTfromHeader)
596 >        if (outfmt == DTfromHeader)     /* check data type */
597                  outfmt = mres->dtype;
598 <        if ((outfmt == DTrgbe) & (mres->ncomp > 3))
599 <                outfmt = DTspec;
600 <        if (outfmt != DTascii)
598 >        if (outfmt == DTrgbe) {
599 >                if (mres->ncomp > 3)
600 >                        outfmt = DTspec;
601 >                else if (mres->dtype == DTxyze)
602 >                        outfmt = DTxyze;
603 >        }
604 >        if (outfmt != DTascii)          /* write result to stdout */
605                  SET_FILE_BINARY(stdout);
606          newheader("RADIANCE", stdout);
607          printargs(argc, argv, stdout);

Diff Legend

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