ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/cmatrix.c
Revision: 2.15
Committed: Wed Feb 17 23:26:06 2016 UTC (8 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.14: +6 -5 lines
Log Message:
Made realloc() calls less aggressive to avoid memory errors

File Contents

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