ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmtxop.c
Revision: 2.3
Committed: Sat Aug 2 17:10:43 2014 UTC (9 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.2: +2 -2 lines
Log Message:
Made rmtxop output the same type as lesser of inputs by default

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: rmtxop.c,v 2.2 2014/05/31 19:21:21 greg Exp $";
3 #endif
4 /*
5 * General component matrix operations.
6 */
7
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include "rtio.h"
11 #include "resolu.h"
12 #include "rmatrix.h"
13
14 #define MAXCOMP 50 /* #components we support */
15
16 typedef struct {
17 double sca[MAXCOMP]; /* scalar coefficients */
18 int nsf; /* number of scalars */
19 double cmat[MAXCOMP*MAXCOMP]; /* component transformation */
20 int clen; /* number of coefficients */
21 int transpose; /* do transpose? */
22 int op; /* '*' or '+' */
23 } ROPERAT; /* matrix operation */
24
25 int outfmt = DTfromHeader; /* output format */
26 int verbose = 0; /* verbose reporting? */
27
28 static void
29 op_default(ROPERAT *op)
30 {
31 memset(op, 0, sizeof(ROPERAT));
32 op->op = '*';
33 }
34
35 static RMATRIX *
36 operate(RMATRIX *mleft, ROPERAT *op, const char *fname)
37 {
38 RMATRIX *mright = rmx_load(fname);
39 RMATRIX *mtmp;
40 int i;
41
42 if (fname == NULL)
43 fname = "<stdin>";
44 if (mright == NULL) {
45 fputs(fname, stderr);
46 fputs(": cannot load matrix\n", stderr);
47 return(NULL);
48 }
49 if (op->transpose) { /* transpose matrix? */
50 mtmp = rmx_transpose(mright);
51 if (mtmp == NULL) {
52 fputs(fname, stderr);
53 fputs(": transpose failed\n", stderr);
54 rmx_free(mright);
55 return(NULL);
56 }
57 if (verbose) {
58 fputs(fname, stderr);
59 fputs(": transposed rows and columns\n", stderr);
60 }
61 rmx_free(mright);
62 mright = mtmp;
63 }
64 if (op->nsf > 0) { /* apply scalar(s) */
65 if (op->clen > 0) {
66 fputs("Options -s and -c are exclusive\n", stderr);
67 rmx_free(mright);
68 return(NULL);
69 }
70 if (op->nsf == 1) {
71 for (i = mright->ncomp; --i; )
72 op->sca[i] = op->sca[0];
73 } else if (op->nsf != mright->ncomp) {
74 fprintf(stderr, "%s: -s must have one or %d factors\n",
75 fname, mright->ncomp);
76 rmx_free(mright);
77 return(NULL);
78 }
79 if ((mleft == NULL) | (op->op != '+') &&
80 !rmx_scale(mright, op->sca)) {
81 fputs(fname, stderr);
82 fputs(": scalar operation failed\n", stderr);
83 rmx_free(mright);
84 return(NULL);
85 }
86 if (verbose) {
87 fputs(fname, stderr);
88 fputs(": applied scalar (", stderr);
89 for (i = 0; i < op->nsf; i++)
90 fprintf(stderr, " %f", op->sca[i]);
91 fputs(" )\n", stderr);
92 }
93 }
94 if (op->clen > 0) { /* apply transform */
95 if (op->clen % mright->ncomp) {
96 fprintf(stderr, "%s: -c must have N x %d coefficients\n",
97 fname, mright->ncomp);
98 rmx_free(mright);
99 return(NULL);
100 }
101 mtmp = rmx_transform(mright, op->clen/mright->ncomp, op->cmat);
102 if (mtmp == NULL) {
103 fprintf(stderr, "%s: matrix transform failed\n", fname);
104 rmx_free(mright);
105 return(NULL);
106 }
107 if (verbose)
108 fprintf(stderr, "%s: applied %d x %d transform\n",
109 fname, mtmp->ncomp, mright->ncomp);
110 rmx_free(mright);
111 mright = mtmp;
112 }
113 if (mleft == NULL) /* just one matrix */
114 return(mright);
115 if (op->op == '*') { /* concatenate */
116 RMATRIX *mres = rmx_multiply(mleft, mright);
117 if (mres == NULL) {
118 fputs(fname, stderr);
119 if (mleft->ncols != mright->nrows)
120 fputs(": mismatched dimensions for multiply\n",
121 stderr);
122 else
123 fputs(": concatenation failed\n", stderr);
124 rmx_free(mright);
125 return(NULL);
126 }
127 if (verbose) {
128 fputs(fname, stderr);
129 fputs(": concatenated matrix\n", stderr);
130 }
131 rmx_free(mright);
132 rmx_free(mleft);
133 mleft = mres;
134 } else if (op->op == '+') {
135 if (!rmx_sum(mleft, mright, op->nsf ? op->sca : (double *)NULL)) {
136 fputs(fname, stderr);
137 fputs(": matrix sum failed\n", stderr);
138 rmx_free(mright);
139 return(NULL);
140 }
141 if (verbose) {
142 fputs(fname, stderr);
143 fputs(": added in matrix\n", stderr);
144 }
145 rmx_free(mright);
146 } else {
147 fprintf(stderr, "%s: unknown operation '%c'\n", fname, op->op);
148 rmx_free(mright);
149 return(NULL);
150 }
151 return(mleft);
152 }
153
154 static int
155 get_factors(double da[], int n, char *av[])
156 {
157 int ac;
158
159 for (ac = 0; ac < n && isflt(av[ac]); ac++)
160 da[ac] = atof(av[ac]);
161 return(ac);
162 }
163
164 /* Load one or more matrices and operate on them, sending results to stdout */
165 int
166 main(int argc, char *argv[])
167 {
168 RMATRIX *mres = NULL;
169 ROPERAT op;
170 long nbw;
171 int i;
172 /* initialize */
173 op_default(&op);
174 /* get options and arguments */
175 for (i = 1; i < argc; i++)
176 if (argv[i][0] == '+' && !argv[i][1]) {
177 op.op = '+';
178 } else if (argv[i][0] != '-' || !argv[i][1]) {
179 char *fname = NULL; /* load matrix */
180 if (argv[i][0] != '-')
181 fname = argv[i];
182 mres = operate(mres, &op, fname);
183 if (mres == NULL) {
184 fprintf(stderr, "%s: operation failed on '%s'\n",
185 argv[0], argv[i]);
186 return(0);
187 }
188 op_default(&op); /* reset operator */
189 } else {
190 int n = argc-1 - i;
191 switch (argv[i][1]) { /* get option */
192 case 'v':
193 verbose = !verbose;
194 break;
195 case 'f':
196 switch (argv[i][2]) {
197 case 'd':
198 outfmt = DTdouble;
199 break;
200 case 'f':
201 outfmt = DTfloat;
202 break;
203 case 'a':
204 outfmt = DTascii;
205 break;
206 case 'c':
207 outfmt = DTrgbe;
208 break;
209 default:
210 goto userr;
211 }
212 break;
213 case 't':
214 op.transpose = 1;
215 break;
216 case 's':
217 if (n > MAXCOMP) n = MAXCOMP;
218 op.nsf = get_factors(op.sca, n, argv+i+1);
219 i += op.nsf;
220 break;
221 case 'c':
222 if (n > MAXCOMP*MAXCOMP) n = MAXCOMP*MAXCOMP;
223 op.clen = get_factors(op.cmat, n, argv+i+1);
224 i += op.clen;
225 break;
226 default:
227 fprintf(stderr, "%s: unknown operation '%s'\n",
228 argv[0], argv[i]);
229 goto userr;
230 }
231 }
232 if (mres == NULL) /* check that we got something */
233 goto userr;
234 /* write result to stdout */
235 newheader("RADIANCE", stdout);
236 printargs(argc, argv, stdout);
237 nbw = rmx_write(mres, outfmt, stdout);
238 /* rmx_free(mres); mres = NULL; */
239 if (nbw <= 0) {
240 fprintf(stderr, "%s: error writing result matrix\n", argv[0]);
241 return(1);
242 }
243 if (verbose)
244 fprintf(stderr, "%s: %ld bytes written\n", argv[0], nbw);
245 return(0);
246 userr:
247 fprintf(stderr,
248 "Usage: %s [-v][-f[adfc][-t][-s sf .. | -c ce ..] m1 [+] .. > mres\n",
249 argv[0]);
250 return(1);
251 }