ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/cmatrix.c
Revision: 2.10
Committed: Fri May 8 18:25:03 2015 UTC (9 years, 5 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.9: +2 -1 lines
Log Message:
Added rtprocess.h needed for popen() under Windows

File Contents

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