ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/cmatrix.c
Revision: 2.1
Committed: Mon Jan 20 21:29:04 2014 UTC (10 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
Split out color coefficient matrix routines from dctimestep

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2     static const char RCSid[] = "$Id$";
3     #endif
4     /*
5     * Color matrix routines.
6     *
7     * G. Ward
8     */
9    
10     #include <ctype.h>
11     #include "standard.h"
12     #include "cmatrix.h"
13     #include "platform.h"
14    
15     /* Allocate a color coefficient matrix */
16     CMATRIX *
17     cm_alloc(int nrows, int ncols)
18     {
19     CMATRIX *cm;
20    
21     if ((nrows <= 0) | (ncols <= 0))
22     error(USER, "attempt to create empty matrix");
23     cm = (CMATRIX *)malloc(sizeof(CMATRIX) +
24     3*sizeof(COLORV)*(nrows*ncols - 1));
25     if (cm == NULL)
26     error(SYSTEM, "out of memory in cm_alloc()");
27     cm->nrows = nrows;
28     cm->ncols = ncols;
29     return(cm);
30     }
31    
32     /* Resize color coefficient matrix */
33     CMATRIX *
34     cm_resize(CMATRIX *cm, int nrows)
35     {
36     if (nrows == cm->nrows)
37     return(cm);
38     if (nrows <= 0) {
39     cm_free(cm);
40     return(NULL);
41     }
42     cm = (CMATRIX *)realloc(cm, sizeof(CMATRIX) +
43     3*sizeof(COLORV)*(nrows*cm->ncols - 1));
44     if (cm == NULL)
45     error(SYSTEM, "out of memory in cm_resize()");
46     cm->nrows = nrows;
47     return(cm);
48     }
49    
50     static int
51     getDT(char *s, void *p)
52     {
53     char fmt[32];
54    
55     if (formatval(fmt, s)) {
56     if (!strcmp(fmt, "ascii"))
57     *((int *)p) = DTascii;
58     else if (!strcmp(fmt, "float"))
59     *((int *)p) = DTfloat;
60     else if (!strcmp(fmt, "double"))
61     *((int *)p) = DTdouble;
62     else if (!strcmp(fmt, COLRFMT))
63     *((int *)p) = DTrgbe;
64     else if (!strcmp(fmt, CIEFMT))
65     *((int *)p) = DTxyze;
66     }
67     return(0);
68     }
69    
70     /* Load header to obtain data type */
71     int
72     getDTfromHeader(FILE *fp)
73     {
74     int dt = DTfromHeader;
75    
76     if (getheader(fp, getDT, &dt) < 0)
77     error(SYSTEM, "header read error");
78     if (dt == DTfromHeader)
79     error(USER, "missing data format in header");
80     return(dt);
81     }
82    
83     /* Allocate and load a matrix from the given file (or stdin if NULL) */
84     CMATRIX *
85     cm_load(const char *fname, int nrows, int ncols, int dtype)
86     {
87     FILE *fp = stdin;
88     CMATRIX *cm;
89    
90     if (ncols <= 0)
91     error(USER, "Non-positive number of columns");
92     if (fname == NULL)
93     fname = "<stdin>";
94     else if ((fp = fopen(fname, "r")) == NULL) {
95     sprintf(errmsg, "cannot open file '%s'", fname);
96     error(SYSTEM, errmsg);
97     }
98     #ifdef getc_unlocked
99     flockfile(fp);
100     #endif
101     if (dtype != DTascii)
102     SET_FILE_BINARY(fp); /* doesn't really work */
103     if (dtype == DTfromHeader)
104     dtype = getDTfromHeader(fp);
105     switch (dtype) {
106     case DTascii:
107     case DTfloat:
108     case DTdouble:
109     break;
110     default:
111     error(USER, "unexpected data type in cm_load()");
112     }
113     if (nrows <= 0) { /* don't know length? */
114     int guessrows = 147; /* usually big enough */
115     if ((dtype != DTascii) & (fp != stdin)) {
116     long startpos = ftell(fp);
117     if (fseek(fp, 0L, SEEK_END) == 0) {
118     long endpos = ftell(fp);
119     long elemsiz = 3*(dtype==DTfloat ?
120     sizeof(float) : sizeof(double));
121    
122     if ((endpos - startpos) % (ncols*elemsiz)) {
123     sprintf(errmsg,
124     "improper length for binary file '%s'",
125     fname);
126     error(USER, errmsg);
127     }
128     guessrows = (endpos - startpos)/(ncols*elemsiz);
129     if (fseek(fp, startpos, SEEK_SET) < 0) {
130     sprintf(errmsg,
131     "fseek() error on file '%s'",
132     fname);
133     error(SYSTEM, errmsg);
134     }
135     nrows = guessrows; /* we're confident */
136     }
137     }
138     cm = cm_alloc(guessrows, ncols);
139     } else
140     cm = cm_alloc(nrows, ncols);
141     if (cm == NULL) /* XXX never happens */
142     return(NULL);
143     if (dtype == DTascii) { /* read text file */
144     int maxrow = (nrows > 0 ? nrows : 32000);
145     int r, c;
146     for (r = 0; r < maxrow; r++) {
147     if (r >= cm->nrows) /* need more space? */
148     cm = cm_resize(cm, 2*cm->nrows);
149     for (c = 0; c < ncols; c++) {
150     COLORV *cv = cm_lval(cm,r,c);
151     if (fscanf(fp, COLSPEC, cv, cv+1, cv+2) != 3)
152     if ((nrows <= 0) & (r > 0) & !c) {
153     cm = cm_resize(cm, maxrow=r);
154     break;
155     } else
156     goto EOFerror;
157     }
158     }
159     while ((c = getc(fp)) != EOF)
160     if (!isspace(c)) {
161     sprintf(errmsg,
162     "unexpected data at end of ascii file %s",
163     fname);
164     error(WARNING, errmsg);
165     break;
166     }
167     } else { /* read binary file */
168     if (sizeof(COLORV) == (dtype==DTfloat ? sizeof(float) :
169     sizeof(double))) {
170     int nread = 0;
171     do { /* read all we can */
172     nread += fread(cm->cmem + 3*nread,
173     3*sizeof(COLORV),
174     cm->nrows*cm->ncols - nread,
175     fp);
176     if (nrows <= 0) { /* unknown length */
177     if (nread == cm->nrows*cm->ncols)
178     /* need more space? */
179     cm = cm_resize(cm, 2*cm->nrows);
180     else if (nread && !(nread % cm->ncols))
181     /* seem to be done */
182     cm = cm_resize(cm, nread/cm->ncols);
183     else /* ended mid-row */
184     goto EOFerror;
185     } else if (nread < cm->nrows*cm->ncols)
186     goto EOFerror;
187     } while (nread < cm->nrows*cm->ncols);
188    
189     } else if (dtype == DTdouble) {
190     double dc[3]; /* load from double */
191     COLORV *cvp = cm->cmem;
192     int n = nrows*ncols;
193    
194     if (n <= 0)
195     goto not_handled;
196     while (n--) {
197     if (fread(dc, sizeof(double), 3, fp) != 3)
198     goto EOFerror;
199     copycolor(cvp, dc);
200     cvp += 3;
201     }
202     } else /* dtype == DTfloat */ {
203     float fc[3]; /* load from float */
204     COLORV *cvp = cm->cmem;
205     int n = nrows*ncols;
206    
207     if (n <= 0)
208     goto not_handled;
209     while (n--) {
210     if (fread(fc, sizeof(float), 3, fp) != 3)
211     goto EOFerror;
212     copycolor(cvp, fc);
213     cvp += 3;
214     }
215     }
216     if (fgetc(fp) != EOF) {
217     sprintf(errmsg,
218     "unexpected data at end of binary file %s",
219     fname);
220     error(WARNING, errmsg);
221     }
222     }
223     if (fp != stdin)
224     fclose(fp);
225     #ifdef getc_unlocked
226     else
227     funlockfile(fp);
228     #endif
229     return(cm);
230     EOFerror:
231     sprintf(errmsg, "unexpected EOF reading %s", fname);
232     error(USER, errmsg);
233     not_handled:
234     error(INTERNAL, "unhandled data size or length in cm_load()");
235     return(NULL); /* gratis return */
236     }
237    
238     /* Extract a column vector from a matrix */
239     CMATRIX *
240     cm_column(const CMATRIX *cm, int c)
241     {
242     CMATRIX *cvr;
243     int dr;
244    
245     if ((c < 0) | (c >= cm->ncols))
246     error(INTERNAL, "column requested outside matrix");
247     cvr = cm_alloc(cm->nrows, 1);
248     if (cvr == NULL)
249     return(NULL);
250     for (dr = 0; dr < cm->nrows; dr++) {
251     const COLORV *sp = cm_lval(cm,dr,c);
252     COLORV *dp = cv_lval(cvr,dr);
253     dp[0] = sp[0];
254     dp[1] = sp[1];
255     dp[2] = sp[2];
256     }
257     return(cvr);
258     }
259    
260     /* Scale a matrix by a single value */
261     CMATRIX *
262     cm_scale(const CMATRIX *cm1, const COLOR sca)
263     {
264     CMATRIX *cmr;
265     int dr, dc;
266    
267     cmr = cm_alloc(cm1->nrows, cm1->ncols);
268     if (cmr == NULL)
269     return(NULL);
270     for (dr = 0; dr < cmr->nrows; dr++)
271     for (dc = 0; dc < cmr->ncols; dc++) {
272     const COLORV *sp = cm_lval(cm1,dr,dc);
273     COLORV *dp = cm_lval(cmr,dr,dc);
274     dp[0] = sp[0] * sca[0];
275     dp[1] = sp[1] * sca[1];
276     dp[2] = sp[2] * sca[2];
277     }
278     return(cmr);
279     }
280    
281     /* Multiply two matrices (or a matrix and a vector) and allocate the result */
282     CMATRIX *
283     cm_multiply(const CMATRIX *cm1, const CMATRIX *cm2)
284     {
285     char *rowcheck=NULL, *colcheck=NULL;
286     CMATRIX *cmr;
287     int dr, dc, i;
288    
289     if ((cm1->ncols <= 0) | (cm1->ncols != cm2->nrows))
290     error(INTERNAL, "matrix dimension mismatch in cm_multiply()");
291     cmr = cm_alloc(cm1->nrows, cm2->ncols);
292     if (cmr == NULL)
293     return(NULL);
294     /* optimization: check for zero rows & cols */
295     if (((cm1->nrows > 5) | (cm2->ncols > 5)) & (cm1->ncols > 5)) {
296     static const COLOR czero;
297     rowcheck = (char *)calloc(cmr->nrows, 1);
298     for (dr = cm1->nrows*(rowcheck != NULL); dr--; )
299     for (dc = cm1->ncols; dc--; )
300     if (memcmp(cm_lval(cm1,dr,dc), czero, sizeof(COLOR))) {
301     rowcheck[dr] = 1;
302     break;
303     }
304     colcheck = (char *)calloc(cmr->ncols, 1);
305     for (dc = cm2->ncols*(colcheck != NULL); dc--; )
306     for (dr = cm2->nrows; dr--; )
307     if (memcmp(cm_lval(cm2,dr,dc), czero, sizeof(COLOR))) {
308     colcheck[dc] = 1;
309     break;
310     }
311     }
312     for (dr = 0; dr < cmr->nrows; dr++)
313     for (dc = 0; dc < cmr->ncols; dc++) {
314     COLORV *dp = cm_lval(cmr,dr,dc);
315     dp[0] = dp[1] = dp[2] = 0;
316     if (rowcheck != NULL && !rowcheck[dr])
317     continue;
318     if (colcheck != NULL && !colcheck[dc])
319     continue;
320     for (i = 0; i < cm1->ncols; i++) {
321     const COLORV *cp1 = cm_lval(cm1,dr,i);
322     const COLORV *cp2 = cm_lval(cm2,i,dc);
323     dp[0] += cp1[0] * cp2[0];
324     dp[1] += cp1[1] * cp2[1];
325     dp[2] += cp1[2] * cp2[2];
326     }
327     }
328     if (rowcheck != NULL) free(rowcheck);
329     if (colcheck != NULL) free(colcheck);
330     return(cmr);
331     }
332    
333     /* print out matrix as ASCII text -- no header */
334     void
335     cm_print(const CMATRIX *cm, FILE *fp)
336     {
337     int r, c;
338     const COLORV *mp = cm->cmem;
339    
340     for (r = 0; r < cm->nrows; r++) {
341     for (c = 0; c < cm->ncols; c++, mp += 3)
342     fprintf(fp, "\t%.6e %.6e %.6e", mp[0], mp[1], mp[2]);
343     fputc('\n', fp);
344     }
345     }