ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/cmatrix.c
Revision: 2.19
Committed: Tue Apr 10 22:11:30 2018 UTC (6 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.18: +2 -2 lines
Log Message:
Fixed inconsistency introduced in 2014 (rev. 2.8)

File Contents

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