ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/cmatrix.c
Revision: 2.35
Committed: Wed Apr 27 01:31:56 2022 UTC (7 months, 2 weeks ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.34: +5 -3 lines
Log Message:
fix(dctimestep): Repaired automatic detection of number of input rows

File Contents

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