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

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.4 static const char RCSid[] = "$Id: cmatrix.c,v 2.3 2014/04/08 23:45:33 greg Exp $";
3 greg 2.1 #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 greg 2.2 #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 greg 2.1
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 greg 2.4 sizeof(COLOR)*(nrows*ncols - 1));
35 greg 2.1 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 greg 2.4 sizeof(COLOR)*(nrows*cm->ncols - 1));
54 greg 2.1 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 greg 2.2 int i;
65 greg 2.1
66 greg 2.2 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 greg 2.1 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 greg 2.2 if (sizeof(COLOR) == cm_elem_size[dtype]) {
173 greg 2.1 int nread = 0;
174     do { /* read all we can */
175     nread += fread(cm->cmem + 3*nread,
176 greg 2.4 sizeof(COLOR),
177 greg 2.1 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 greg 2.4 double res[3];
319 greg 2.1 dp[0] = dp[1] = dp[2] = 0;
320     if (rowcheck != NULL && !rowcheck[dr])
321     continue;
322     if (colcheck != NULL && !colcheck[dc])
323     continue;
324 greg 2.4 res[0] = res[1] = res[2] = 0;
325 greg 2.1 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 greg 2.4 res[0] += cp1[0] * cp2[0];
329     res[1] += cp1[1] * cp2[1];
330     res[2] += cp1[2] * cp2[2];
331 greg 2.1 }
332 greg 2.4 copycolor(dp, res);
333 greg 2.1 }
334     if (rowcheck != NULL) free(rowcheck);
335     if (colcheck != NULL) free(colcheck);
336     return(cmr);
337     }
338    
339 greg 2.2 /* write out matrix to file (precede by resolution string if picture) */
340     int
341     cm_write(const CMATRIX *cm, int dtype, FILE *fp)
342 greg 2.1 {
343 greg 2.3 static const char tabEOL[2] = {'\t','\n'};
344     const COLORV *mp = cm->cmem;
345     int r, c;
346 greg 2.2
347     switch (dtype) {
348     case DTascii:
349 greg 2.3 for (r = 0; r < cm->nrows; r++)
350 greg 2.2 for (c = 0; c < cm->ncols; c++, mp += 3)
351 greg 2.3 fprintf(fp, "%.6e %.6e %.6e%c",
352     mp[0], mp[1], mp[2],
353     tabEOL[c >= cm->ncols-1]);
354 greg 2.2 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 greg 2.1 }
397 greg 2.2 return(fflush(fp) == 0);
398 greg 2.1 }