ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmtxop.c
Revision: 2.18
Committed: Tue Jan 19 23:32:00 2021 UTC (3 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.17: +12 -6 lines
Log Message:
feat: Added -rf and -rb options to rmtxop to load XML reflectance matrices

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.18 static const char RCSid[] = "$Id: rmtxop.c,v 2.17 2019/12/28 18:05:14 greg Exp $";
3 greg 2.1 #endif
4     /*
5     * General component matrix operations.
6     */
7    
8     #include <stdlib.h>
9 greg 2.11 #include <errno.h>
10 greg 2.1 #include "rtio.h"
11     #include "resolu.h"
12     #include "rmatrix.h"
13 greg 2.10 #include "platform.h"
14 greg 2.1
15 greg 2.13 #define MAXCOMP 16 /* #components we support */
16 greg 2.1
17 greg 2.13 /* unary matrix operation(s) */
18 greg 2.1 typedef struct {
19     double sca[MAXCOMP]; /* scalar coefficients */
20     double cmat[MAXCOMP*MAXCOMP]; /* component transformation */
21 greg 2.13 short nsf; /* number of scalars */
22     short clen; /* number of coefficients */
23     short transpose; /* do transpose? */
24     } RUNARYOP;
25    
26     /* matrix input source and requested operation(s) */
27     typedef struct {
28     const char *inspec; /* input specification */
29 greg 2.18 RMPref rmp; /* matrix preference */
30 greg 2.13 RUNARYOP preop; /* unary operation(s) */
31     RMATRIX *mtx; /* original matrix if loaded */
32     int binop; /* binary op with next (or 0) */
33     } ROPMAT;
34 greg 2.1
35     int verbose = 0; /* verbose reporting? */
36    
37 greg 2.13 /* Load matrix */
38     static int
39     loadmatrix(ROPMAT *rop)
40 greg 2.1 {
41 greg 2.13 if (rop->mtx != NULL)
42     return(0);
43    
44 greg 2.18 rop->mtx = rmx_load(rop->inspec, rop->rmp);
45 greg 2.13 if (rop->mtx == NULL) {
46     fputs(rop->inspec, stderr);
47     fputs(": cannot load matrix\n", stderr);
48     return(-1);
49     }
50     return(1);
51 greg 2.1 }
52    
53 greg 2.13 /* Get matrix and perform unary operations */
54 greg 2.1 static RMATRIX *
55 greg 2.13 loadop(ROPMAT *rop)
56 greg 2.1 {
57 greg 2.13 RMATRIX *mres;
58 greg 2.1 int i;
59    
60 greg 2.13 if (loadmatrix(rop) < 0) /* make sure we're loaded */
61 greg 2.1 return(NULL);
62 greg 2.13
63     if (rop->preop.nsf > 0) { /* apply scalar(s) */
64     if (rop->preop.clen > 0) {
65 greg 2.1 fputs("Options -s and -c are exclusive\n", stderr);
66 greg 2.13 goto failure;
67 greg 2.1 }
68 greg 2.13 if (rop->preop.nsf == 1) {
69     for (i = rop->mtx->ncomp; --i; )
70     rop->preop.sca[i] = rop->preop.sca[0];
71     } else if (rop->preop.nsf != rop->mtx->ncomp) {
72 greg 2.1 fprintf(stderr, "%s: -s must have one or %d factors\n",
73 greg 2.13 rop->inspec, rop->mtx->ncomp);
74     goto failure;
75 greg 2.1 }
76 greg 2.13 if (!rmx_scale(rop->mtx, rop->preop.sca)) {
77     fputs(rop->inspec, stderr);
78 greg 2.1 fputs(": scalar operation failed\n", stderr);
79 greg 2.13 goto failure;
80 greg 2.1 }
81     if (verbose) {
82 greg 2.13 fputs(rop->inspec, stderr);
83 greg 2.1 fputs(": applied scalar (", stderr);
84 greg 2.13 for (i = 0; i < rop->preop.nsf; i++)
85     fprintf(stderr, " %f", rop->preop.sca[i]);
86 greg 2.1 fputs(" )\n", stderr);
87     }
88     }
89 greg 2.13 if (rop->preop.clen > 0) { /* apply transform */
90     if (rop->preop.clen % rop->mtx->ncomp) {
91 greg 2.1 fprintf(stderr, "%s: -c must have N x %d coefficients\n",
92 greg 2.13 rop->inspec, rop->mtx->ncomp);
93     goto failure;
94 greg 2.1 }
95 greg 2.13 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 greg 2.1 }
102     if (verbose)
103     fprintf(stderr, "%s: applied %d x %d transform\n",
104 greg 2.13 rop->inspec, mres->ncomp,
105     rop->mtx->ncomp);
106     rmx_free(rop->mtx);
107     rop->mtx = mres;
108 greg 2.1 }
109 greg 2.13 if (rop->preop.transpose) { /* transpose matrix? */
110     mres = rmx_transpose(rop->mtx);
111     if (mres == NULL) {
112     fputs(rop->inspec, stderr);
113 greg 2.12 fputs(": transpose failed\n", stderr);
114 greg 2.13 goto failure;
115 greg 2.12 }
116     if (verbose) {
117 greg 2.13 fputs(rop->inspec, stderr);
118 greg 2.12 fputs(": transposed rows and columns\n", stderr);
119     }
120 greg 2.13 rmx_free(rop->mtx);
121     rop->mtx = mres;
122     }
123     mres = rop->mtx;
124     rop->mtx = NULL;
125     return(mres);
126     failure:
127     rmx_free(rop->mtx);
128     return(rop->mtx = NULL);
129     }
130    
131     /* Execute binary operation, free matrix arguments and return new result */
132     static RMATRIX *
133     binaryop(const char *inspec, RMATRIX *mleft, int op, RMATRIX *mright)
134     {
135     RMATRIX *mres = NULL;
136     int i;
137    
138     if ((mleft == NULL) | (mright == NULL))
139     return(NULL);
140     switch (op) {
141     case '.': /* concatenate */
142 greg 2.16 if (mleft->ncomp != mright->ncomp) {
143     fputs(inspec, stderr);
144     fputs(": # components do not match\n", stderr);
145     } else if (mleft->ncols != mright->nrows) {
146     fputs(inspec, stderr);
147     fputs(": mismatched dimensions\n",
148     stderr);
149     } else
150     mres = rmx_multiply(mleft, mright);
151 greg 2.13 rmx_free(mleft);
152 greg 2.12 rmx_free(mright);
153 greg 2.1 if (mres == NULL) {
154 greg 2.13 fputs(inspec, stderr);
155 greg 2.16 fputs(": concatenation failed\n", stderr);
156 greg 2.1 return(NULL);
157     }
158     if (verbose) {
159 greg 2.13 fputs(inspec, stderr);
160 greg 2.1 fputs(": concatenated matrix\n", stderr);
161     }
162 greg 2.13 break;
163     case '+':
164     if (!rmx_sum(mleft, mright, NULL)) {
165     fputs(inspec, stderr);
166 greg 2.1 fputs(": matrix sum failed\n", stderr);
167 greg 2.13 rmx_free(mleft);
168 greg 2.1 rmx_free(mright);
169     return(NULL);
170     }
171     if (verbose) {
172 greg 2.13 fputs(inspec, stderr);
173 greg 2.1 fputs(": added in matrix\n", stderr);
174     }
175     rmx_free(mright);
176 greg 2.13 mres = mleft;
177     break;
178     case '*':
179     case '/': {
180     const char * tnam = (op == '/') ?
181 greg 2.11 "division" : "multiplication";
182     errno = 0;
183 greg 2.13 if (!rmx_elemult(mleft, mright, (op == '/'))) {
184 greg 2.11 fprintf(stderr, "%s: element-wise %s failed\n",
185 greg 2.13 inspec, tnam);
186     rmx_free(mleft);
187 greg 2.11 rmx_free(mright);
188     return(NULL);
189     }
190     if (errno)
191     fprintf(stderr,
192     "%s: warning - error during element-wise %s\n",
193 greg 2.13 inspec, tnam);
194 greg 2.11 else if (verbose)
195 greg 2.13 fprintf(stderr, "%s: element-wise %s\n", inspec, tnam);
196 greg 2.11 rmx_free(mright);
197 greg 2.13 mres = mleft;
198     } break;
199     default:
200     fprintf(stderr, "%s: unknown operation '%c'\n", inspec, op);
201     rmx_free(mleft);
202 greg 2.1 rmx_free(mright);
203     return(NULL);
204     }
205 greg 2.13 return(mres);
206     }
207    
208     /* Perform matrix operations from left to right */
209     static RMATRIX *
210     op_left2right(ROPMAT *mop)
211     {
212     RMATRIX *mleft = loadop(mop);
213    
214     while (mop->binop) {
215     if (mleft == NULL)
216     break;
217     mleft = binaryop(mop[1].inspec,
218     mleft, mop->binop, loadop(mop+1));
219     mop++;
220     }
221 greg 2.1 return(mleft);
222     }
223    
224 greg 2.13 /* Perform matrix operations from right to left */
225     static RMATRIX *
226     op_right2left(ROPMAT *mop)
227     {
228     RMATRIX *mright;
229     int rpos = 0;
230     /* find end of list */
231     while (mop[rpos].binop)
232 greg 2.14 if (mop[rpos++].binop != '.') {
233     fputs(
234     "Right-to-left evaluation only for matrix multiplication!\n",
235     stderr);
236     return(NULL);
237     }
238 greg 2.13 mright = loadop(mop+rpos);
239     while (rpos-- > 0) {
240     if (mright == NULL)
241     break;
242     mright = binaryop(mop[rpos].inspec,
243     loadop(mop+rpos), mop[rpos].binop, mright);
244     }
245     return(mright);
246     }
247    
248     #define t_nrows(mop) ((mop)->preop.transpose ? (mop)->mtx->ncols \
249     : (mop)->mtx->nrows)
250     #define t_ncols(mop) ((mop)->preop.transpose ? (mop)->mtx->nrows \
251     : (mop)->mtx->ncols)
252    
253     /* Should we prefer concatenating from rightmost matrix towards left? */
254     static int
255     prefer_right2left(ROPMAT *mop)
256     {
257     int mri = 0;
258    
259     while (mop[mri].binop) /* find rightmost matrix */
260     if (mop[mri++].binop != '.')
261     return(0); /* pre-empt reversal for other ops */
262    
263     if (mri <= 1)
264     return(0); /* won't matter */
265    
266     if (loadmatrix(mop+mri) < 0) /* load rightmost cat */
267     return(1); /* fail will bail in a moment */
268    
269     if (t_ncols(mop+mri) == 1)
270     return(1); /* definitely better R->L */
271    
272     if (t_ncols(mop+mri) >= t_nrows(mop+mri))
273     return(0); /* ...probably worse */
274    
275     if (loadmatrix(mop) < 0) /* load leftmost */
276     return(0); /* fail will bail in a moment */
277    
278     return(t_ncols(mop+mri) < t_nrows(mop));
279     }
280    
281 greg 2.1 static int
282     get_factors(double da[], int n, char *av[])
283     {
284     int ac;
285    
286     for (ac = 0; ac < n && isflt(av[ac]); ac++)
287     da[ac] = atof(av[ac]);
288     return(ac);
289     }
290    
291 greg 2.13 static ROPMAT *
292     grow_moparray(ROPMAT *mop, int n2alloc)
293     {
294     int nmats = 0;
295    
296     while (mop[nmats++].binop)
297     ;
298     mop = (ROPMAT *)realloc(mop, n2alloc*sizeof(ROPMAT));
299     if (mop == NULL) {
300     fputs("Out of memory in grow_moparray()\n", stderr);
301     exit(1);
302     }
303     if (n2alloc > nmats)
304     memset(mop+nmats, 0, (n2alloc-nmats)*sizeof(ROPMAT));
305     return(mop);
306     }
307    
308 greg 2.1 /* Load one or more matrices and operate on them, sending results to stdout */
309     int
310     main(int argc, char *argv[])
311     {
312 greg 2.7 int outfmt = DTfromHeader;
313 greg 2.13 int nall = 2;
314     ROPMAT *mop = (ROPMAT *)calloc(nall, sizeof(ROPMAT));
315     int nmats = 0;
316 greg 2.1 RMATRIX *mres = NULL;
317 greg 2.13 int stdin_used = 0;
318 greg 2.1 int i;
319     /* get options and arguments */
320 greg 2.13 for (i = 1; i < argc; i++) {
321 greg 2.11 if (argv[i][0] && !argv[i][1] &&
322 greg 2.13 strchr(".+*/", argv[i][0]) != NULL) {
323 greg 2.15 if (!nmats || mop[nmats-1].binop) {
324 greg 2.14 fprintf(stderr,
325 greg 2.15 "%s: missing matrix argument before '%c' operation\n",
326 greg 2.14 argv[0], argv[i][0]);
327 greg 2.13 return(1);
328     }
329 greg 2.15 mop[nmats-1].binop = argv[i][0];
330 greg 2.1 } else if (argv[i][0] != '-' || !argv[i][1]) {
331 greg 2.13 if (argv[i][0] == '-') {
332     if (stdin_used++) {
333     fprintf(stderr,
334     "%s: standard input used for more than one matrix\n",
335     argv[0]);
336     return(1);
337     }
338     mop[nmats].inspec = stdin_name;
339     } else
340     mop[nmats].inspec = argv[i];
341     if (nmats > 0 && !mop[nmats-1].binop)
342     mop[nmats-1].binop = '.';
343     nmats++;
344 greg 2.1 } else {
345     int n = argc-1 - i;
346     switch (argv[i][1]) { /* get option */
347     case 'v':
348 greg 2.14 verbose++;
349 greg 2.1 break;
350     case 'f':
351     switch (argv[i][2]) {
352     case 'd':
353     outfmt = DTdouble;
354     break;
355     case 'f':
356     outfmt = DTfloat;
357     break;
358     case 'a':
359     outfmt = DTascii;
360     break;
361     case 'c':
362     outfmt = DTrgbe;
363     break;
364     default:
365     goto userr;
366     }
367     break;
368     case 't':
369 greg 2.13 mop[nmats].preop.transpose = 1;
370 greg 2.1 break;
371     case 's':
372     if (n > MAXCOMP) n = MAXCOMP;
373 greg 2.13 i += mop[nmats].preop.nsf =
374     get_factors(mop[nmats].preop.sca,
375     n, argv+i+1);
376 greg 2.1 break;
377     case 'c':
378     if (n > MAXCOMP*MAXCOMP) n = MAXCOMP*MAXCOMP;
379 greg 2.13 i += mop[nmats].preop.clen =
380     get_factors(mop[nmats].preop.cmat,
381     n, argv+i+1);
382 greg 2.1 break;
383 greg 2.18 case 'r':
384     if (argv[i][2] == 'f')
385     mop[nmats].rmp = RMPreflF;
386     else if (argv[i][2] == 'b')
387     mop[nmats].rmp = RMPreflB;
388     else
389     goto userr;
390     break;
391 greg 2.1 default:
392     fprintf(stderr, "%s: unknown operation '%s'\n",
393     argv[0], argv[i]);
394     goto userr;
395     }
396     }
397 greg 2.14 if (nmats >= nall)
398     mop = grow_moparray(mop, nall += 2);
399 greg 2.13 }
400     if (mop[0].inspec == NULL) /* nothing to do? */
401 greg 2.1 goto userr;
402 greg 2.15 if (mop[nmats-1].binop) {
403     fprintf(stderr,
404     "%s: missing matrix argument after '%c' operation\n",
405     argv[0], mop[nmats-1].binop);
406     return(1);
407     }
408 greg 2.13 /* favor quicker concatenation */
409 greg 2.14 mop[nmats].mtx = prefer_right2left(mop) ? op_right2left(mop)
410     : op_left2right(mop);
411     if (mop[nmats].mtx == NULL)
412     return(1);
413     /* apply trailing unary operations */
414     mop[nmats].inspec = "trailing_ops";
415     mres = loadop(mop+nmats);
416     if (mres == NULL)
417 greg 2.13 return(1);
418 greg 2.1 /* write result to stdout */
419 greg 2.6 if (outfmt == DTfromHeader)
420     outfmt = mres->dtype;
421 greg 2.4 if (outfmt != DTascii)
422 greg 2.10 SET_FILE_BINARY(stdout);
423 greg 2.1 newheader("RADIANCE", stdout);
424     printargs(argc, argv, stdout);
425 greg 2.5 if (!rmx_write(mres, outfmt, stdout)) {
426 greg 2.1 fprintf(stderr, "%s: error writing result matrix\n", argv[0]);
427     return(1);
428     }
429 greg 2.13 /* rmx_free(mres); free(mop); */
430 greg 2.1 return(0);
431     userr:
432     fprintf(stderr,
433 greg 2.18 "Usage: %s [-v][-f[adfc][-t][-s sf .. | -c ce ..][-r[fb]] m1 [.+*/] .. > mres\n",
434 greg 2.1 argv[0]);
435     return(1);
436     }