ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/cmatrix.c
Revision: 2.9
Committed: Mon May 4 20:53:21 2015 UTC (8 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.8: +29 -17 lines
Log Message:
Added the ability to read matrix inputs from commands as well as files

File Contents

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