ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmtxop.c
Revision: 2.7
Committed: Fri Jan 23 01:03:17 2015 UTC (9 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.6: +2 -2 lines
Log Message:
Minor fix

File Contents

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