ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/cmatrix.c
Revision: 2.1
Committed: Mon Jan 20 21:29:04 2014 UTC (10 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
Split out color coefficient matrix routines from dctimestep

File Contents

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