ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/cmatrix.c
Revision: 2.8
Committed: Sat Aug 2 17:10:43 2014 UTC (9 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R2P2, rad4R2P1
Changes since 2.7: +3 -3 lines
Log Message:
Made rmtxop output the same type as lesser of inputs by default

File Contents

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