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.27 by greg, Sat Dec 2 00:42:21 2023 UTC vs.
Revision 2.42 by greg, Sat Apr 19 03:58:00 2025 UTC

# Line 7 | Line 7 | static const char RCSid[] = "$Id$";
7  
8   #include <errno.h>
9   #include "rtio.h"
10 – #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 {
# Line 35 | 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 +                        if (!rop->mtx)
55 +                                return(-1);
56 +                        rop->mtx->dtype = DTascii;
57 +                        return(1);      /* loaded BSDF XML file */
58 +                }
59 +        }                               /* else load regular matrix */
60 +        rop->mtx = rmx_load(rop->inspec);
61  
62 <        rop->mtx = rmx_load(rop->inspec, rop->rmp);
45 <
46 <        return(!rop->mtx ? -1 : 1);
62 >        return(rop->mtx ? 1 : -1);
63   }
64  
65 < static int      checksymbolic(ROPMAT *rop);
65 > extern int      checksymbolic(ROPMAT *rop);
66  
67   /* Check/set transform based on a reference input file */
68 < static int
68 > int
69   checkreffile(ROPMAT *rop)
70   {
71          static const char       *curRF = NULL;
# Line 69 | Line 85 | checkreffile(ROPMAT *rop)
85                  fclose(fp);
86                  curRF = rop->preop.csym;
87          }
88 <        if ((refm.ncomp == 3) & (refm.dtype != DTspec)) {
88 >        if (refm.ncomp == 3) {
89                  rop->preop.csym = (refm.dtype == DTxyze) ? "XYZ" : "RGB";
90                  return(checksymbolic(rop));
91          }
# Line 103 | Line 119 | checkreffile(ROPMAT *rop)
119                  for (j = refm.ncomp; j-- > 0; )
120                          rop->preop.cmat[j*nc + i] = scresp[j];
121          }
122 +        memcpy(rop->mtx->wlpart, refm.wlpart, sizeof(rop->mtx->wlpart));
123          return(0);
124   }
125  
126   /* Compute conversion row from spectrum to one channel of RGB */
127 < static void
127 > void
128   rgbrow(ROPMAT *rop, int r, int p)
129   {
130          const int       nc = rop->mtx->ncomp;
# Line 124 | Line 141 | rgbrow(ROPMAT *rop, int r, int p)
141   }
142  
143   /* Compute conversion row from spectrum to one channel of XYZ */
144 < static void
144 > void
145   xyzrow(ROPMAT *rop, int r, int p)
146   {
147          const int       nc = rop->mtx->ncomp;
# Line 141 | Line 158 | xyzrow(ROPMAT *rop, int r, int p)
158   }
159  
160   /* Use the spectral sensitivity function to compute matrix coefficients */
161 < static void
162 < sensrow(ROPMAT *rop, int r, double (*sf)(SCOLOR sc, int ncs, const float wlpt[4]))
161 > void
162 > sensrow(ROPMAT *rop, int r, double (*sf)(const SCOLOR sc, int ncs, const float wlpt[4]))
163   {
164          const int       nc = rop->mtx->ncomp;
165          int             i;
# Line 156 | Line 173 | sensrow(ROPMAT *rop, int r, double (*sf)(SCOLOR sc, in
173   }
174  
175   /* Check/set symbolic transform */
176 < static int
176 > int
177   checksymbolic(ROPMAT *rop)
178   {
179          const int       nc = rop->mtx->ncomp;
180          const int       dt = rop->mtx->dtype;
181 +        double          cf = 1;
182          int             i, j;
183                                          /* check suffix => reference file */
184          if (strchr(rop->preop.csym, '.') > rop->preop.csym)
# Line 181 | Line 199 | checksymbolic(ROPMAT *rop)
199                  int     comp = 0;
200                  switch (rop->preop.csym[j]) {
201                  case 'B':
202 +                case 'b':
203                          ++comp;
204                          /* fall through */
205                  case 'G':
206 +                case 'g':
207                          ++comp;
208                          /* fall through */
209                  case 'R':
210 +                case 'r':
211 +                        if (rop->preop.csym[j] <= 'Z')
212 +                                cf = 1./WHTEFFICACY;
213                          if (dt == DTxyze) {
214                                  for (i = 3; i--; )
215 <                                        rop->preop.cmat[j*nc+i] = 1./WHTEFFICACY *
193 <                                                        xyz2rgbmat[comp][i];
215 >                                        rop->preop.cmat[j*nc+i] = cf*xyz2rgbmat[comp][i];
216                          } else if (nc == 3)
217                                  rop->preop.cmat[j*nc+comp] = 1.;
218                          else
219                                  rgbrow(rop, j, comp);
220                          break;
221                  case 'Z':
222 +                case 'z':
223                          ++comp;
224                          /* fall through */
225                  case 'Y':
226 +                case 'y':
227                          ++comp;
228                          /* fall through */
229                  case 'X':
230 +                case 'x':
231 +                        if ((rop->preop.csym[j] <= 'Z') & (dt != DTxyze))
232 +                                cf = WHTEFFICACY;
233                          if (dt == DTxyze) {
234                                  rop->preop.cmat[j*nc+comp] = 1.;
235                          } else if (nc == 3) {
# Line 214 | Line 241 | checksymbolic(ROPMAT *rop)
241                          else
242                                  xyzrow(rop, j, comp);
243  
244 <                        for (i = nc*(dt != DTxyze); i--; )
245 <                                rop->preop.cmat[j*nc+i] *= WHTEFFICACY;
244 >                        for (i = nc*(cf != 1); i--; )
245 >                                rop->preop.cmat[j*nc+i] *= cf;
246                          break;
247                  case 'S':               /* scotopic (il)luminance */
248 +                        cf = WHTSCOTOPIC;
249 +                        /* fall through */
250 +                case 's':
251                          sensrow(rop, j, scolor2scotopic);
252 <                        for (i = nc; i--; )
253 <                                rop->preop.cmat[j*nc+i] *= WHTSCOTOPIC;
252 >                        for (i = nc*(cf != 1); i--; )
253 >                                rop->preop.cmat[j*nc+i] *= cf;
254                          break;
255                  case 'M':               /* melanopic (il)luminance */
256 +                        cf = WHTMELANOPIC;
257 +                        /* fall through */
258 +                case 'm':
259                          sensrow(rop, j, scolor2melanopic);
260 <                        for (i = nc; i--; )
261 <                                rop->preop.cmat[j*nc+i] *= WHTMELANOPIC;
260 >                        for (i = nc*(cf != 1); i--; )
261 >                                rop->preop.cmat[j*nc+i] *= cf;
262                          break;
263                  case 'A':               /* average component */
264 +                case 'a':
265                          for (i = nc; i--; )
266                                  rop->preop.cmat[j*nc+i] = 1./(double)nc;
267                          break;
# Line 238 | Line 272 | checksymbolic(ROPMAT *rop)
272                  }
273          }
274                                          /* return recommended output type */
275 <        if (!strcmp(rop->preop.csym, "XYZ")) {
275 >        if (!strcasecmp(rop->preop.csym, "XYZ")) {
276                  if (dt <= DTspec)
277                          return(DTxyze);
278 <        } else if (!strcmp(rop->preop.csym, "RGB")) {
278 >        } else if (!strcasecmp(rop->preop.csym, "RGB")) {
279                  if (dt <= DTspec)
280                          return(DTrgbe);
281 <        }
248 <        if ((nc > 3) & (dt <= DTspec))
281 >        } else if (dt == DTspec)
282                  return(DTfloat);        /* probably not actual spectrum */
283          return(0);
284   }
285  
286   /* Get matrix and perform unary operations */
287 < static RMATRIX *
287 > RMATRIX *
288   loadop(ROPMAT *rop)
289   {
290          int     outtype = 0;
# Line 274 | Line 307 | loadop(ROPMAT *rop)
307                          if (rop->preop.nsf == 1) {
308                                  for (i = rop->preop.clen; i--; )
309                                          rop->preop.cmat[i] *= rop->preop.sca[0];
310 <                        } else if (rop->preop.nsf != rop->mtx->ncomp) {
310 >                        } else if (rop->preop.nsf*rop->mtx->ncomp != rop->preop.clen) {
311                                  fprintf(stderr, "%s: -s must have one or %d factors\n",
312 <                                                rop->inspec, rop->mtx->ncomp);
312 >                                                rop->inspec,
313 >                                                rop->preop.clen/rop->mtx->ncomp);
314                                  goto failure;
315                          } else {
316 <                                for (j = rop->preop.clen/rop->preop.nsf; j--; )
317 <                                        for (i = rop->preop.nsf; i--; )
318 <                                                rop->preop.cmat[j*rop->preop.nsf+i] *=
319 <                                                                rop->preop.sca[i];
316 >                                for (i = rop->preop.nsf; i--; )
317 >                                        for (j = rop->mtx->ncomp; j--; )
318 >                                                rop->preop.cmat[i*rop->mtx->ncomp+j]
319 >                                                                *= rop->preop.sca[i];
320                          }
321                  }
322                  mres = rmx_transform(rop->mtx, rop->preop.clen/rop->mtx->ncomp,
# Line 326 | Line 360 | loadop(ROPMAT *rop)
360                  }
361          }
362          if (rop->preop.transpose) {             /* transpose matrix? */
363 <                mres = rmx_transpose(rop->mtx);
330 <                if (mres == NULL) {
363 >                if (!rmx_transpose(rop->mtx)) {
364                          fputs(rop->inspec, stderr);
365                          fputs(": transpose failed\n", stderr);
366                          goto failure;
# Line 336 | Line 369 | loadop(ROPMAT *rop)
369                          fputs(rop->inspec, stderr);
370                          fputs(": transposed rows and columns\n", stderr);
371                  }
339 –                rmx_free(rop->mtx);
340 –                rop->mtx = mres;
372          }
373          mres = rop->mtx;
374          rop->mtx = NULL;
# Line 350 | Line 381 | failure:
381   }
382  
383   /* Execute binary operation, free matrix arguments and return new result */
384 < static RMATRIX *
384 > RMATRIX *
385   binaryop(const char *inspec, RMATRIX *mleft, int op, RMATRIX *mright)
386   {
387          RMATRIX *mres = NULL;
# Line 427 | Line 458 | binaryop(const char *inspec, RMATRIX *mleft, int op, R
458   }
459  
460   /* Perform matrix operations from left to right */
461 < static RMATRIX *
461 > RMATRIX *
462   op_left2right(ROPMAT *mop)
463   {
464          RMATRIX *mleft = loadop(mop);
# Line 443 | Line 474 | op_left2right(ROPMAT *mop)
474   }
475  
476   /* Perform matrix operations from right to left */
477 < static RMATRIX *
477 > RMATRIX *
478   op_right2left(ROPMAT *mop)
479   {
480          RMATRIX *mright;
# Line 472 | Line 503 | op_right2left(ROPMAT *mop)
503                                                  : (mop)->mtx->ncols)
504  
505   /* Should we prefer concatenating from rightmost matrix towards left? */
506 < static int
506 > int
507   prefer_right2left(ROPMAT *mop)
508   {
509          int     mri = 0;
# Line 499 | Line 530 | prefer_right2left(ROPMAT *mop)
530          return(t_ncols(mop+mri) < t_nrows(mop));
531   }
532  
533 < static int
533 > int
534   get_factors(double da[], int n, char *av[])
535   {
536          int     ac;
# Line 509 | Line 540 | get_factors(double da[], int n, char *av[])
540          return(ac);
541   }
542  
543 < static ROPMAT *
543 > ROPMAT *
544   resize_moparr(ROPMAT *mop, int n2alloc)
545   {
546          int     nmats = 0;
# Line 517 | Line 548 | resize_moparr(ROPMAT *mop, int n2alloc)
548  
549          while (mop[nmats++].binop)
550                  ;
551 <        for (i = nmats; i > n2alloc; i--)
551 >        for (i = nmats; i >= n2alloc; i--)
552                  rmx_free(mop[i].mtx);
553          mop = (ROPMAT *)realloc(mop, n2alloc*sizeof(ROPMAT));
554          if (mop == NULL) {
# Line 664 | Line 695 | main(int argc, char *argv[])
695          mres = loadop(mop+nmats);
696          if (mres == NULL)
697                  return(1);
698 <        if (outfmt == DTfromHeader)     /* check data type */
698 >        if ((outfmt == DTfromHeader) & (mres->dtype < DTspec))
699                  outfmt = mres->dtype;
700 <        if (outfmt == DTrgbe) {
700 >        if (outfmt == DTrgbe) {         /* check data type */
701                  if (mres->ncomp > 3)
702                          outfmt = DTspec;
703                  else if (mres->dtype == DTxyze)
704                          outfmt = DTxyze;
705          }
706 + #if DTrmx_native==DTfloat
707 +        if (outfmt == DTdouble)
708 +                fprintf(stderr,
709 +                        "%s: warning - writing float result as double\n",
710 +                                argv[0]);
711 + #endif
712          newheader("RADIANCE", stdout);  /* write result to stdout */
713          printargs(argc, argv, stdout);
714          return(rmx_write(mres, outfmt, stdout) ? 0 : 1);

Diff Legend

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