ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/cmatrix.c
Revision: 2.5
Committed: Fri May 30 00:00:54 2014 UTC (9 years, 10 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.4: +62 -18 lines
Log Message:
Added NROWS, NCOLS and NCOMP to matrix headers

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.5 static const char RCSid[] = "$Id: cmatrix.c,v 2.4 2014/05/29 17:28:09 greg Exp $";
3 greg 2.1 #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 greg 2.2 #include "resolu.h"
15    
16     const char *cm_fmt_id[] = {
17     "unknown", "ascii", "float", "double",
18     COLRFMT, CIEFMT
19     };
20    
21     const int cm_elem_size[] = {
22     0, 0, 3*sizeof(float), 3*sizeof(double), 4, 4
23     };
24 greg 2.1
25     /* Allocate a color coefficient matrix */
26     CMATRIX *
27     cm_alloc(int nrows, int ncols)
28     {
29     CMATRIX *cm;
30    
31     if ((nrows <= 0) | (ncols <= 0))
32     error(USER, "attempt to create empty matrix");
33     cm = (CMATRIX *)malloc(sizeof(CMATRIX) +
34 greg 2.4 sizeof(COLOR)*(nrows*ncols - 1));
35 greg 2.1 if (cm == NULL)
36     error(SYSTEM, "out of memory in cm_alloc()");
37     cm->nrows = nrows;
38     cm->ncols = ncols;
39     return(cm);
40     }
41    
42     /* Resize color coefficient matrix */
43     CMATRIX *
44     cm_resize(CMATRIX *cm, int nrows)
45     {
46     if (nrows == cm->nrows)
47     return(cm);
48     if (nrows <= 0) {
49     cm_free(cm);
50     return(NULL);
51     }
52     cm = (CMATRIX *)realloc(cm, sizeof(CMATRIX) +
53 greg 2.4 sizeof(COLOR)*(nrows*cm->ncols - 1));
54 greg 2.1 if (cm == NULL)
55     error(SYSTEM, "out of memory in cm_resize()");
56     cm->nrows = nrows;
57     return(cm);
58     }
59    
60 greg 2.5 typedef struct {
61     int dtype; /* data type */
62     int nrows, ncols; /* matrix size */
63     char *err; /* error message */
64     } CMINFO; /* header info record */
65    
66 greg 2.1 static int
67 greg 2.5 get_cminfo(char *s, void *p)
68 greg 2.1 {
69 greg 2.5 CMINFO *ip = (CMINFO *)p;
70 greg 2.1 char fmt[32];
71 greg 2.2 int i;
72 greg 2.5
73     if (!strncmp(s, "NCOMP=", 6) && atoi(s+6) != 3) {
74     ip->err = "unexpected # components (must be 3)";
75     return(-1);
76     }
77     if (!strncmp(s, "NROWS=", 6)) {
78     ip->nrows = atoi(s+6);
79     return(0);
80     }
81     if (!strncmp(s, "NCOLS=", 6)) {
82     ip->ncols = atoi(s+6);
83     return(0);
84     }
85 greg 2.2 if (!formatval(fmt, s))
86     return(0);
87     for (i = 1; i < DTend; i++)
88     if (!strcmp(fmt, cm_fmt_id[i]))
89 greg 2.5 ip->dtype = i;
90 greg 2.1 return(0);
91     }
92    
93 greg 2.5 /* Load header to obtain/check data type and number of columns */
94     char *
95     cm_getheader(int *dt, int *nr, int *nc, FILE *fp)
96 greg 2.1 {
97 greg 2.5 CMINFO cmi;
98     /* read header */
99     cmi.dtype = DTfromHeader;
100     cmi.nrows = cmi.ncols = 0;
101     cmi.err = "unexpected EOF in header";
102     if (getheader(fp, &get_cminfo, &cmi) < 0)
103     return(cmi.err);
104     if (dt != NULL) { /* get/check data type? */
105     if (cmi.dtype == DTfromHeader) {
106     if (*dt == DTfromHeader)
107     return("missing/unknown data format in header");
108     } else if (*dt == DTfromHeader)
109     *dt = cmi.dtype;
110     else if (*dt != cmi.dtype)
111     return("unexpected data format in header");
112     }
113     if (nr != NULL) { /* get/check #rows? */
114     if (*nr <= 0)
115     *nr = cmi.nrows;
116     else if ((cmi.nrows > 0) & (*nr != cmi.nrows))
117     return("unexpected row count in header");
118     }
119     if (nc != NULL) { /* get/check #columns? */
120     if (*nc <= 0)
121     *nc = cmi.ncols;
122     else if ((cmi.ncols > 0) & (*nc != cmi.ncols))
123     return("unexpected column count in header");
124     }
125     return(NULL);
126 greg 2.1 }
127    
128     /* Allocate and load a matrix from the given file (or stdin if NULL) */
129     CMATRIX *
130     cm_load(const char *fname, int nrows, int ncols, int dtype)
131     {
132     FILE *fp = stdin;
133     CMATRIX *cm;
134    
135     if (fname == NULL)
136     fname = "<stdin>";
137     else if ((fp = fopen(fname, "r")) == NULL) {
138     sprintf(errmsg, "cannot open file '%s'", fname);
139     error(SYSTEM, errmsg);
140     }
141     #ifdef getc_unlocked
142     flockfile(fp);
143     #endif
144     if (dtype != DTascii)
145     SET_FILE_BINARY(fp); /* doesn't really work */
146 greg 2.5 if (!dtype | !ncols) { /* expecting header? */
147     char *err = cm_getheader(&dtype, &nrows, &ncols, fp);
148     if (err != NULL)
149     error(USER, err);
150     if (ncols <= 0)
151     error(USER, "unspecified number of columns");
152     }
153 greg 2.1 switch (dtype) {
154     case DTascii:
155     case DTfloat:
156     case DTdouble:
157     break;
158     default:
159     error(USER, "unexpected data type in cm_load()");
160     }
161     if (nrows <= 0) { /* don't know length? */
162     int guessrows = 147; /* usually big enough */
163     if ((dtype != DTascii) & (fp != stdin)) {
164     long startpos = ftell(fp);
165     if (fseek(fp, 0L, SEEK_END) == 0) {
166     long endpos = ftell(fp);
167     long elemsiz = 3*(dtype==DTfloat ?
168     sizeof(float) : sizeof(double));
169    
170     if ((endpos - startpos) % (ncols*elemsiz)) {
171     sprintf(errmsg,
172     "improper length for binary file '%s'",
173     fname);
174     error(USER, errmsg);
175     }
176     guessrows = (endpos - startpos)/(ncols*elemsiz);
177     if (fseek(fp, startpos, SEEK_SET) < 0) {
178     sprintf(errmsg,
179     "fseek() error on file '%s'",
180     fname);
181     error(SYSTEM, errmsg);
182     }
183     nrows = guessrows; /* we're confident */
184     }
185     }
186     cm = cm_alloc(guessrows, ncols);
187     } else
188     cm = cm_alloc(nrows, ncols);
189     if (cm == NULL) /* XXX never happens */
190     return(NULL);
191     if (dtype == DTascii) { /* read text file */
192     int maxrow = (nrows > 0 ? nrows : 32000);
193     int r, c;
194     for (r = 0; r < maxrow; r++) {
195     if (r >= cm->nrows) /* need more space? */
196     cm = cm_resize(cm, 2*cm->nrows);
197     for (c = 0; c < ncols; c++) {
198     COLORV *cv = cm_lval(cm,r,c);
199     if (fscanf(fp, COLSPEC, cv, cv+1, cv+2) != 3)
200     if ((nrows <= 0) & (r > 0) & !c) {
201     cm = cm_resize(cm, maxrow=r);
202     break;
203     } else
204     goto EOFerror;
205     }
206     }
207     while ((c = getc(fp)) != EOF)
208     if (!isspace(c)) {
209     sprintf(errmsg,
210     "unexpected data at end of ascii file %s",
211     fname);
212     error(WARNING, errmsg);
213     break;
214     }
215     } else { /* read binary file */
216 greg 2.2 if (sizeof(COLOR) == cm_elem_size[dtype]) {
217 greg 2.1 int nread = 0;
218     do { /* read all we can */
219     nread += fread(cm->cmem + 3*nread,
220 greg 2.4 sizeof(COLOR),
221 greg 2.1 cm->nrows*cm->ncols - nread,
222     fp);
223     if (nrows <= 0) { /* unknown length */
224     if (nread == cm->nrows*cm->ncols)
225     /* need more space? */
226     cm = cm_resize(cm, 2*cm->nrows);
227     else if (nread && !(nread % cm->ncols))
228     /* seem to be done */
229     cm = cm_resize(cm, nread/cm->ncols);
230     else /* ended mid-row */
231     goto EOFerror;
232     } else if (nread < cm->nrows*cm->ncols)
233     goto EOFerror;
234     } while (nread < cm->nrows*cm->ncols);
235    
236     } else if (dtype == DTdouble) {
237     double dc[3]; /* load from double */
238     COLORV *cvp = cm->cmem;
239     int n = nrows*ncols;
240    
241     if (n <= 0)
242     goto not_handled;
243     while (n--) {
244     if (fread(dc, sizeof(double), 3, fp) != 3)
245     goto EOFerror;
246     copycolor(cvp, dc);
247     cvp += 3;
248     }
249     } else /* dtype == DTfloat */ {
250     float fc[3]; /* load from float */
251     COLORV *cvp = cm->cmem;
252     int n = nrows*ncols;
253    
254     if (n <= 0)
255     goto not_handled;
256     while (n--) {
257     if (fread(fc, sizeof(float), 3, fp) != 3)
258     goto EOFerror;
259     copycolor(cvp, fc);
260     cvp += 3;
261     }
262     }
263     if (fgetc(fp) != EOF) {
264     sprintf(errmsg,
265     "unexpected data at end of binary file %s",
266     fname);
267     error(WARNING, errmsg);
268     }
269     }
270     if (fp != stdin)
271     fclose(fp);
272     #ifdef getc_unlocked
273     else
274     funlockfile(fp);
275     #endif
276     return(cm);
277     EOFerror:
278     sprintf(errmsg, "unexpected EOF reading %s", fname);
279     error(USER, errmsg);
280     not_handled:
281     error(INTERNAL, "unhandled data size or length in cm_load()");
282     return(NULL); /* gratis return */
283     }
284    
285     /* Extract a column vector from a matrix */
286     CMATRIX *
287     cm_column(const CMATRIX *cm, int c)
288     {
289     CMATRIX *cvr;
290     int dr;
291    
292     if ((c < 0) | (c >= cm->ncols))
293     error(INTERNAL, "column requested outside matrix");
294     cvr = cm_alloc(cm->nrows, 1);
295     if (cvr == NULL)
296     return(NULL);
297     for (dr = 0; dr < cm->nrows; dr++) {
298     const COLORV *sp = cm_lval(cm,dr,c);
299     COLORV *dp = cv_lval(cvr,dr);
300     dp[0] = sp[0];
301     dp[1] = sp[1];
302     dp[2] = sp[2];
303     }
304     return(cvr);
305     }
306    
307     /* Scale a matrix by a single value */
308     CMATRIX *
309     cm_scale(const CMATRIX *cm1, const COLOR sca)
310     {
311     CMATRIX *cmr;
312     int dr, dc;
313    
314     cmr = cm_alloc(cm1->nrows, cm1->ncols);
315     if (cmr == NULL)
316     return(NULL);
317     for (dr = 0; dr < cmr->nrows; dr++)
318     for (dc = 0; dc < cmr->ncols; dc++) {
319     const COLORV *sp = cm_lval(cm1,dr,dc);
320     COLORV *dp = cm_lval(cmr,dr,dc);
321     dp[0] = sp[0] * sca[0];
322     dp[1] = sp[1] * sca[1];
323     dp[2] = sp[2] * sca[2];
324     }
325     return(cmr);
326     }
327    
328     /* Multiply two matrices (or a matrix and a vector) and allocate the result */
329     CMATRIX *
330     cm_multiply(const CMATRIX *cm1, const CMATRIX *cm2)
331     {
332     char *rowcheck=NULL, *colcheck=NULL;
333     CMATRIX *cmr;
334     int dr, dc, i;
335    
336     if ((cm1->ncols <= 0) | (cm1->ncols != cm2->nrows))
337     error(INTERNAL, "matrix dimension mismatch in cm_multiply()");
338     cmr = cm_alloc(cm1->nrows, cm2->ncols);
339     if (cmr == NULL)
340     return(NULL);
341     /* optimization: check for zero rows & cols */
342     if (((cm1->nrows > 5) | (cm2->ncols > 5)) & (cm1->ncols > 5)) {
343     static const COLOR czero;
344     rowcheck = (char *)calloc(cmr->nrows, 1);
345     for (dr = cm1->nrows*(rowcheck != NULL); dr--; )
346     for (dc = cm1->ncols; dc--; )
347     if (memcmp(cm_lval(cm1,dr,dc), czero, sizeof(COLOR))) {
348     rowcheck[dr] = 1;
349     break;
350     }
351     colcheck = (char *)calloc(cmr->ncols, 1);
352     for (dc = cm2->ncols*(colcheck != NULL); dc--; )
353     for (dr = cm2->nrows; dr--; )
354     if (memcmp(cm_lval(cm2,dr,dc), czero, sizeof(COLOR))) {
355     colcheck[dc] = 1;
356     break;
357     }
358     }
359     for (dr = 0; dr < cmr->nrows; dr++)
360     for (dc = 0; dc < cmr->ncols; dc++) {
361     COLORV *dp = cm_lval(cmr,dr,dc);
362 greg 2.4 double res[3];
363 greg 2.1 dp[0] = dp[1] = dp[2] = 0;
364     if (rowcheck != NULL && !rowcheck[dr])
365     continue;
366     if (colcheck != NULL && !colcheck[dc])
367     continue;
368 greg 2.4 res[0] = res[1] = res[2] = 0;
369 greg 2.1 for (i = 0; i < cm1->ncols; i++) {
370     const COLORV *cp1 = cm_lval(cm1,dr,i);
371     const COLORV *cp2 = cm_lval(cm2,i,dc);
372 greg 2.4 res[0] += cp1[0] * cp2[0];
373     res[1] += cp1[1] * cp2[1];
374     res[2] += cp1[2] * cp2[2];
375 greg 2.1 }
376 greg 2.4 copycolor(dp, res);
377 greg 2.1 }
378     if (rowcheck != NULL) free(rowcheck);
379     if (colcheck != NULL) free(colcheck);
380     return(cmr);
381     }
382    
383 greg 2.2 /* write out matrix to file (precede by resolution string if picture) */
384     int
385     cm_write(const CMATRIX *cm, int dtype, FILE *fp)
386 greg 2.1 {
387 greg 2.3 static const char tabEOL[2] = {'\t','\n'};
388     const COLORV *mp = cm->cmem;
389     int r, c;
390 greg 2.2
391     switch (dtype) {
392     case DTascii:
393 greg 2.3 for (r = 0; r < cm->nrows; r++)
394 greg 2.2 for (c = 0; c < cm->ncols; c++, mp += 3)
395 greg 2.3 fprintf(fp, "%.6e %.6e %.6e%c",
396     mp[0], mp[1], mp[2],
397     tabEOL[c >= cm->ncols-1]);
398 greg 2.2 break;
399     case DTfloat:
400     case DTdouble:
401     if (sizeof(COLOR) == cm_elem_size[dtype]) {
402     r = cm->ncols*cm->nrows;
403     while (r > 0) {
404     c = fwrite(mp, sizeof(COLOR), r, fp);
405     if (c <= 0)
406     return(0);
407     mp += 3*c;
408     r -= c;
409     }
410     } else if (dtype == DTdouble) {
411     double dc[3];
412     r = cm->ncols*cm->nrows;
413     while (r--) {
414     copycolor(dc, mp);
415     if (fwrite(dc, sizeof(double), 3, fp) != 3)
416     return(0);
417     mp += 3;
418     }
419     } else /* dtype == DTfloat */ {
420     float fc[3];
421     r = cm->ncols*cm->nrows;
422     while (r--) {
423     copycolor(fc, mp);
424     if (fwrite(fc, sizeof(float), 3, fp) != 3)
425     return(0);
426     mp += 3;
427     }
428     }
429     break;
430     case DTrgbe:
431     case DTxyze:
432     fprtresolu(cm->ncols, cm->nrows, fp);
433     for (r = 0; r < cm->nrows; r++, mp += 3*cm->ncols)
434     if (fwritescan((COLOR *)mp, cm->ncols, fp) < 0)
435     return(0);
436     break;
437     default:
438     fputs("Unsupported data type in cm_write()!\n", stderr);
439     return(0);
440 greg 2.1 }
441 greg 2.2 return(fflush(fp) == 0);
442 greg 2.1 }