ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/cmatrix.c
Revision: 2.21
Committed: Mon Aug 12 18:28:37 2019 UTC (4 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.20: +27 -39 lines
Log Message:
Made error checking more consistent and eliminated unused function call

File Contents

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