ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmtxop.c
Revision: 2.1
Committed: Sat May 31 05:02:37 2014 UTC (9 years, 10 months ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
Created rmtxop command, need to test and debug

File Contents

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