ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmtxop.c
Revision: 2.24
Committed: Wed Nov 29 18:56:28 2023 UTC (4 months, 4 weeks ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.23: +8 -4 lines
Log Message:
feat(rmtxop): Added 'A' symbolic component for average of all channels

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: rmtxop.c,v 2.23 2023/11/28 16:36:50 greg Exp $";
3 #endif
4 /*
5 * General component matrix operations.
6 */
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 MAXCSAMP /* #components we support */
17
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 char csym[11]; /* symbolic coefficients */
25 char transpose; /* do transpose? */
26 } RUNARYOP;
27
28 /* Matrix input source and requested operation(s) */
29 typedef struct {
30 const char *inspec; /* input specification */
31 RMPref rmp; /* matrix preference */
32 RUNARYOP preop; /* unary operation(s) */
33 RMATRIX *mtx; /* original matrix if loaded */
34 int binop; /* binary op with next (or 0) */
35 } ROPMAT;
36
37 int verbose = 0; /* verbose reporting? */
38
39 /* Load matrix */
40 static int
41 loadmatrix(ROPMAT *rop)
42 {
43 if (rop->mtx != NULL) /* already loaded? */
44 return(0);
45
46 rop->mtx = rmx_load(rop->inspec, rop->rmp);
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.f;
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 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': /* scotopic (il)luminance */
160 sensrow(rop, j, scolor2scotopic);
161 for (i = nc; i--; )
162 rop->preop.cmat[j*nc+i] *= WHTSCOTOPIC;
163 break;
164 case 'M': /* melanopic (il)luminance */
165 sensrow(rop, j, scolor2melanopic);
166 for (i = nc; i--; )
167 rop->preop.cmat[j*nc+i] *= WHTMELANOPIC;
168 break;
169 case 'A': /* average component */
170 for (i = nc; i--; )
171 rop->preop->cmat[j*nc+i] = 1./(double)nc;
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, j;
199
200 if (loadmatrix(rop) < 0) /* make sure we're loaded */
201 return(NULL);
202
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];
249 } else if (rop->preop.nsf != rop->mtx->ncomp) {
250 fprintf(stderr, "%s: -s must have one or %d factors\n",
251 rop->inspec, rop->mtx->ncomp);
252 goto failure;
253 }
254 if (!rmx_scale(rop->mtx, rop->preop.sca)) {
255 fputs(rop->inspec, stderr);
256 fputs(": scalar operation failed\n", stderr);
257 goto failure;
258 }
259 if (verbose) {
260 fputs(rop->inspec, stderr);
261 fputs(": applied scalar (", stderr);
262 for (i = 0; i < rop->preop.nsf; i++)
263 fprintf(stderr, " %f", rop->preop.sca[i]);
264 fputs(" )\n", stderr);
265 }
266 }
267 if (rop->preop.transpose) { /* transpose matrix? */
268 mres = rmx_transpose(rop->mtx);
269 if (mres == NULL) {
270 fputs(rop->inspec, stderr);
271 fputs(": transpose failed\n", stderr);
272 goto failure;
273 }
274 if (verbose) {
275 fputs(rop->inspec, stderr);
276 fputs(": transposed rows and columns\n", stderr);
277 }
278 rmx_free(rop->mtx);
279 rop->mtx = mres;
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);
288 return(rop->mtx = NULL);
289 }
290
291 /* Execute binary operation, free matrix arguments and return new result */
292 static RMATRIX *
293 binaryop(const char *inspec, RMATRIX *mleft, int op, RMATRIX *mright)
294 {
295 RMATRIX *mres = NULL;
296 int i;
297
298 if ((mleft == NULL) | (mright == NULL))
299 return(NULL);
300 switch (op) {
301 case '.': /* concatenate */
302 if (mleft->ncomp != mright->ncomp) {
303 fputs(inspec, stderr);
304 fputs(": # components do not match\n", stderr);
305 } else if (mleft->ncols != mright->nrows) {
306 fputs(inspec, stderr);
307 fputs(": mismatched dimensions\n",
308 stderr);
309 } else
310 mres = rmx_multiply(mleft, mright);
311 rmx_free(mleft);
312 rmx_free(mright);
313 if (mres == NULL) {
314 fputs(inspec, stderr);
315 fputs(": concatenation failed\n", stderr);
316 return(NULL);
317 }
318 if (verbose) {
319 fputs(inspec, stderr);
320 fputs(": concatenated matrix\n", stderr);
321 }
322 break;
323 case '+':
324 if (!rmx_sum(mleft, mright, NULL)) {
325 fputs(inspec, stderr);
326 fputs(": matrix sum failed\n", stderr);
327 rmx_free(mleft);
328 rmx_free(mright);
329 return(NULL);
330 }
331 if (verbose) {
332 fputs(inspec, stderr);
333 fputs(": added in matrix\n", stderr);
334 }
335 rmx_free(mright);
336 mres = mleft;
337 break;
338 case '*':
339 case '/': {
340 const char * tnam = (op == '/') ?
341 "division" : "multiplication";
342 errno = 0;
343 if (!rmx_elemult(mleft, mright, (op == '/'))) {
344 fprintf(stderr, "%s: element-wise %s failed\n",
345 inspec, tnam);
346 rmx_free(mleft);
347 rmx_free(mright);
348 return(NULL);
349 }
350 if (errno)
351 fprintf(stderr,
352 "%s: warning - error during element-wise %s\n",
353 inspec, tnam);
354 else if (verbose)
355 fprintf(stderr, "%s: element-wise %s\n", inspec, tnam);
356 rmx_free(mright);
357 mres = mleft;
358 } break;
359 default:
360 fprintf(stderr, "%s: unknown operation '%c'\n", inspec, op);
361 rmx_free(mleft);
362 rmx_free(mright);
363 return(NULL);
364 }
365 return(mres);
366 }
367
368 /* Perform matrix operations from left to right */
369 static RMATRIX *
370 op_left2right(ROPMAT *mop)
371 {
372 RMATRIX *mleft = loadop(mop);
373
374 while (mop->binop) {
375 if (mleft == NULL)
376 break;
377 mleft = binaryop(mop[1].inspec,
378 mleft, mop->binop, loadop(mop+1));
379 mop++;
380 }
381 return(mleft);
382 }
383
384 /* Perform matrix operations from right to left */
385 static RMATRIX *
386 op_right2left(ROPMAT *mop)
387 {
388 RMATRIX *mright;
389 int rpos = 0;
390 /* find end of list */
391 while (mop[rpos].binop)
392 if (mop[rpos++].binop != '.') {
393 fputs(
394 "Right-to-left evaluation only for matrix multiplication!\n",
395 stderr);
396 return(NULL);
397 }
398 mright = loadop(mop+rpos);
399 while (rpos-- > 0) {
400 if (mright == NULL)
401 break;
402 mright = binaryop(mop[rpos+1].inspec,
403 loadop(mop+rpos), mop[rpos].binop, mright);
404 }
405 return(mright);
406 }
407
408 #define t_nrows(mop) ((mop)->preop.transpose ? (mop)->mtx->ncols \
409 : (mop)->mtx->nrows)
410 #define t_ncols(mop) ((mop)->preop.transpose ? (mop)->mtx->nrows \
411 : (mop)->mtx->ncols)
412
413 /* Should we prefer concatenating from rightmost matrix towards left? */
414 static int
415 prefer_right2left(ROPMAT *mop)
416 {
417 int mri = 0;
418
419 while (mop[mri].binop) /* find rightmost matrix */
420 if (mop[mri++].binop != '.')
421 return(0); /* pre-empt reversal for other ops */
422
423 if (mri <= 1)
424 return(0); /* won't matter */
425
426 if (loadmatrix(mop+mri) < 0) /* load rightmost cat */
427 return(1); /* fail will bail in a moment */
428
429 if (t_ncols(mop+mri) == 1)
430 return(1); /* definitely better R->L */
431
432 if (t_ncols(mop+mri) >= t_nrows(mop+mri))
433 return(0); /* ...probably worse */
434
435 if (loadmatrix(mop) < 0) /* load leftmost */
436 return(0); /* fail will bail in a moment */
437
438 return(t_ncols(mop+mri) < t_nrows(mop));
439 }
440
441 static int
442 get_factors(double da[], int n, char *av[])
443 {
444 int ac;
445
446 for (ac = 0; ac < n && isflt(av[ac]); ac++)
447 da[ac] = atof(av[ac]);
448 return(ac);
449 }
450
451 static ROPMAT *
452 grow_moparray(ROPMAT *mop, int n2alloc)
453 {
454 int nmats = 0;
455
456 while (mop[nmats++].binop)
457 ;
458 mop = (ROPMAT *)realloc(mop, n2alloc*sizeof(ROPMAT));
459 if (mop == NULL) {
460 fputs("Out of memory in grow_moparray()\n", stderr);
461 exit(1);
462 }
463 if (n2alloc > nmats)
464 memset(mop+nmats, 0, (n2alloc-nmats)*sizeof(ROPMAT));
465 return(mop);
466 }
467
468 /* Load one or more matrices and operate on them, sending results to stdout */
469 int
470 main(int argc, char *argv[])
471 {
472 int outfmt = DTfromHeader;
473 int nall = 2;
474 ROPMAT *mop = (ROPMAT *)calloc(nall, sizeof(ROPMAT));
475 int nmats = 0;
476 RMATRIX *mres = NULL;
477 int stdin_used = 0;
478 int i;
479 /* get options and arguments */
480 for (i = 1; i < argc; i++) {
481 if (argv[i][0] && !argv[i][1] &&
482 strchr(".+*/", argv[i][0]) != NULL) {
483 if (!nmats || mop[nmats-1].binop) {
484 fprintf(stderr,
485 "%s: missing matrix argument before '%c' operation\n",
486 argv[0], argv[i][0]);
487 return(1);
488 }
489 mop[nmats-1].binop = argv[i][0];
490 } else if (argv[i][0] != '-' || !argv[i][1]) {
491 if (argv[i][0] == '-') {
492 if (stdin_used++) {
493 fprintf(stderr,
494 "%s: standard input used for more than one matrix\n",
495 argv[0]);
496 return(1);
497 }
498 mop[nmats].inspec = stdin_name;
499 } else
500 mop[nmats].inspec = argv[i];
501 if (nmats > 0 && !mop[nmats-1].binop)
502 mop[nmats-1].binop = '.';
503 nmats++;
504 } else {
505 int n = argc-1 - i;
506 switch (argv[i][1]) { /* get option */
507 case 'v':
508 verbose++;
509 break;
510 case 'f':
511 switch (argv[i][2]) {
512 case 'd':
513 outfmt = DTdouble;
514 break;
515 case 'f':
516 outfmt = DTfloat;
517 break;
518 case 'a':
519 outfmt = DTascii;
520 break;
521 case 'c':
522 outfmt = DTrgbe;
523 break;
524 default:
525 goto userr;
526 }
527 break;
528 case 't':
529 mop[nmats].preop.transpose = 1;
530 break;
531 case 's':
532 if (n > MAXCOMP) n = MAXCOMP;
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')
563 mop[nmats].rmp = RMPreflF;
564 else if (argv[i][2] == 'b')
565 mop[nmats].rmp = RMPreflB;
566 else
567 goto userr;
568 break;
569 default:
570 fprintf(stderr, "%s: unknown operation '%s'\n",
571 argv[0], argv[i]);
572 goto userr;
573 }
574 }
575 if (nmats >= nall)
576 mop = grow_moparray(mop, nall += 2);
577 }
578 if (mop[0].inspec == NULL) /* nothing to do? */
579 goto userr;
580 if (mop[nmats-1].binop) {
581 fprintf(stderr,
582 "%s: missing matrix argument after '%c' operation\n",
583 argv[0], mop[nmats-1].binop);
584 return(1);
585 }
586 /* favor quicker concatenation */
587 mop[nmats].mtx = prefer_right2left(mop) ? op_right2left(mop)
588 : op_left2right(mop);
589 if (mop[nmats].mtx == NULL)
590 return(1);
591 /* apply trailing unary operations */
592 mop[nmats].inspec = "trailing_ops";
593 mres = loadop(mop+nmats);
594 if (mres == NULL)
595 return(1);
596 if (outfmt == DTfromHeader) /* check data type */
597 outfmt = mres->dtype;
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);
608
609 return(rmx_write(mres, outfmt, stdout) ? 0 : 1);
610 userr:
611 fprintf(stderr,
612 "Usage: %s [-v][-f[adfc][-t][-s sf .. | -c ce ..][-rf|-rb] m1 [.+*/] .. > mres\n",
613 argv[0]);
614 return(1);
615 }