ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/cmatrix.c
Revision: 2.23
Committed: Wed Oct 23 17:00:14 2019 UTC (4 years, 5 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.22: +33 -1 lines
Log Message:
Added support for RGBE picture i/o to dctimestep

File Contents

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