ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/cmatrix.c
Revision: 2.4
Committed: Thu May 29 17:28:09 2014 UTC (9 years, 10 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.3: +10 -7 lines
Log Message:
Improved accuracy of matrix multiplication by using double accumulator

File Contents

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