ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmatrix.c
Revision: 2.3
Committed: Tue Jul 8 16:39:41 2014 UTC (9 years, 9 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.2: +2 -2 lines
Log Message:
Removed unnecessary address operators from function ids.

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: rmatrix.c,v 2.2 2014/05/31 19:21:21 greg Exp $";
3 #endif
4 /*
5 * General matrix operations.
6 */
7
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include <string.h>
11 #include <fcntl.h>
12 #include "resolu.h"
13 #include "rmatrix.h"
14
15 typedef struct {
16 int nrows, ncols, ncomp;
17 int dtype;
18 } DMINFO;
19
20 /* Allocate a nr x nc matrix with n components */
21 RMATRIX *
22 rmx_alloc(int nr, int nc, int n)
23 {
24 RMATRIX *dnew;
25
26 if ((nr <= 0) | (nc <= 0) | (n <= 0))
27 return(NULL);
28 dnew = (RMATRIX *)malloc(sizeof(RMATRIX)-sizeof(dnew->mtx) +
29 sizeof(dnew->mtx[0])*(n*nr*nc));
30 if (dnew == NULL)
31 return(NULL);
32 dnew->nrows = nr; dnew->ncols = nc; dnew->ncomp = n;
33 return(dnew);
34 }
35
36 static int
37 get_dminfo(char *s, void *p)
38 {
39 DMINFO *ip = (DMINFO *)p;
40 char fmt[32];
41 int i;
42
43 if (!strncmp(s, "NCOMP=", 6)) {
44 ip->ncomp = atoi(s+6);
45 return(0);
46 }
47 if (!strncmp(s, "NROWS=", 6)) {
48 ip->nrows = atoi(s+6);
49 return(0);
50 }
51 if (!strncmp(s, "NCOLS=", 6)) {
52 ip->ncols = atoi(s+6);
53 return(0);
54 }
55 if (!formatval(fmt, s))
56 return(0);
57 for (i = 1; i < DTend; i++)
58 if (!strcmp(fmt, cm_fmt_id[i])) {
59 ip->dtype = i;
60 return(0);
61 }
62 return(-1);
63 }
64
65 static int
66 rmx_load_ascii(RMATRIX *rm, FILE *fp)
67 {
68 int i, j, k;
69 #ifdef _WIN32
70 _setmode(fileno(fp), _O_TEXT);
71 #endif
72 for (i = 0; i < rm->nrows; i++)
73 for (j = 0; j < rm->ncols; j++)
74 for (k = 0; k < rm->ncomp; k++)
75 if (fscanf(fp, "%lf", &rmx_lval(rm,i,j,k)) != 1)
76 return(0);
77 return(1);
78 }
79
80 static int
81 rmx_load_float(RMATRIX *rm, FILE *fp)
82 {
83 int i, j, k;
84 float val[100];
85
86 if (rm->ncomp > 100) {
87 fputs("Unsupported # components in rmx_load_float()\n", stderr);
88 exit(1);
89 }
90 for (i = 0; i < rm->nrows; i++)
91 for (j = 0; j < rm->ncols; j++) {
92 if (fread(val, sizeof(val[0]), rm->ncomp, fp) != rm->ncomp)
93 return(0);
94 for (k = rm->ncomp; k--; )
95 rmx_lval(rm,i,j,k) = val[k];
96 }
97 return(1);
98 }
99
100 static int
101 rmx_load_double(RMATRIX *rm, FILE *fp)
102 {
103 int i, j, k;
104 double val[100];
105
106 if (rm->ncomp > 100) {
107 fputs("Unsupported # components in rmx_load_double()\n", stderr);
108 exit(1);
109 }
110 for (i = 0; i < rm->nrows; i++)
111 for (j = 0; j < rm->ncols; j++) {
112 if (fread(val, sizeof(val[0]), rm->ncomp, fp) != rm->ncomp)
113 return(0);
114 for (k = rm->ncomp; k--; )
115 rmx_lval(rm,i,j,k) = val[k];
116 }
117 return(1);
118 }
119
120 static int
121 rmx_load_rgbe(RMATRIX *rm, FILE *fp)
122 {
123 COLOR *scan = (COLOR *)malloc(sizeof(COLOR)*rm->ncols);
124 int i, j;
125
126 if (scan == NULL)
127 return(0);
128 for (i = 0; i < rm->nrows; i++) {
129 if (freadscan(scan, rm->ncols, fp) < 0) {
130 free(scan);
131 return(0);
132 }
133 for (j = rm->ncols; j--; ) {
134 rmx_lval(rm,i,j,0) = colval(scan[j],RED);
135 rmx_lval(rm,i,j,1) = colval(scan[j],GRN);
136 rmx_lval(rm,i,j,2) = colval(scan[j],BLU);
137 }
138 }
139 free(scan);
140 return(1);
141 }
142
143 /* Load matrix from supported file type */
144 RMATRIX *
145 rmx_load(const char *fname)
146 {
147 FILE *fp = stdin;
148 DMINFO dinfo;
149 RMATRIX *dnew;
150
151 if (fname == NULL) { /* reading from stdin? */
152 fname = "<stdin>";
153 } else {
154 const char *sp = fname; /* check suffix */
155 while (*sp)
156 ++sp;
157 while (sp > fname && sp[-1] != '.')
158 --sp;
159 if (!strcasecmp(sp, "XML")) { /* assume it's a BSDF */
160 CMATRIX *cm = cm_loadBTDF((char *)fname);
161 if (cm == NULL)
162 return(NULL);
163 dnew = rmx_from_cmatrix(cm);
164 cm_free(cm);
165 return(dnew);
166 }
167 /* else open it ourselves */
168 if ((fp = fopen(fname, "rb")) == NULL)
169 return(NULL);
170 }
171 #ifdef getc_unlocked
172 flockfile(fp);
173 #endif
174 dinfo.nrows = dinfo.ncols = dinfo.ncomp = 0;
175 dinfo.dtype = DTascii;
176 if (getheader(fp, get_dminfo, &dinfo) < 0) {
177 fclose(fp);
178 return(NULL);
179 }
180 if ((dinfo.dtype == DTrgbe) | (dinfo.dtype == DTxyze)) {
181 if (!fscnresolu(&dinfo.ncols, &dinfo.nrows, fp)) {
182 fclose(fp);
183 return(NULL);
184 }
185 dinfo.ncomp = 3;
186 }
187 dnew = rmx_alloc(dinfo.nrows, dinfo.ncols, dinfo.ncomp);
188 if (dnew == NULL) {
189 fclose(fp);
190 return(NULL);
191 }
192 switch (dinfo.dtype) {
193 case DTascii:
194 if (!rmx_load_ascii(dnew, fp))
195 goto loaderr;
196 break;
197 case DTfloat:
198 if (!rmx_load_float(dnew, fp))
199 goto loaderr;
200 break;
201 case DTdouble:
202 if (!rmx_load_double(dnew, fp))
203 goto loaderr;
204 break;
205 case DTrgbe:
206 case DTxyze:
207 if (!rmx_load_rgbe(dnew, fp))
208 goto loaderr;
209 break;
210 default:
211 goto loaderr;
212 }
213 if (fp != stdin)
214 fclose(fp);
215 return(dnew);
216 loaderr: /* should report error? */
217 fclose(fp);
218 rmx_free(dnew);
219 return(NULL);
220 }
221
222 static int
223 rmx_write_ascii(const RMATRIX *rm, FILE *fp)
224 {
225 int i, j, k;
226 #ifdef _WIN32
227 _setmode(fileno(fp), _O_TEXT);
228 #endif
229 for (i = 0; i < rm->nrows; i++) {
230 for (j = 0; j < rm->ncols; j++) {
231 for (k = 0; k < rm->ncomp; k++)
232 fprintf(fp, " %.15e", rmx_lval(rm,i,j,k));
233 fputc('\t', fp);
234 }
235 fputc('\n', fp);
236 }
237 return(1);
238 }
239
240 static int
241 rmx_write_float(const RMATRIX *rm, FILE *fp)
242 {
243 int i, j, k;
244 float val[100];
245
246 if (rm->ncomp > 100) {
247 fputs("Unsupported # components in rmx_write_float()\n", stderr);
248 exit(1);
249 }
250 for (i = 0; i < rm->nrows; i++)
251 for (j = 0; j < rm->ncols; j++) {
252 for (k = rm->ncomp; k--; )
253 val[k] = (float)rmx_lval(rm,i,j,k);
254 if (fwrite(val, sizeof(val[0]), rm->ncomp, fp) != rm->ncomp)
255 return(0);
256 }
257 return(1);
258 }
259
260 static int
261 rmx_write_double(const RMATRIX *rm, FILE *fp)
262 {
263 int i, j, k;
264 double val[100];
265
266 if (rm->ncomp > 100) {
267 fputs("Unsupported # components in rmx_write_double()\n", stderr);
268 exit(1);
269 }
270 for (i = 0; i < rm->nrows; i++)
271 for (j = 0; j < rm->ncols; j++) {
272 for (k = rm->ncomp; k--; )
273 val[k] = rmx_lval(rm,i,j,k);
274 if (fwrite(val, sizeof(val[0]), rm->ncomp, fp) != rm->ncomp)
275 return(0);
276 }
277 return(1);
278 }
279
280 static int
281 rmx_write_rgbe(const RMATRIX *rm, FILE *fp)
282 {
283 COLOR *scan = (COLOR *)malloc(sizeof(COLOR)*rm->ncols);
284 int i, j;
285
286 if (scan == NULL)
287 return(0);
288 for (i = 0; i < rm->nrows; i++) {
289 for (j = rm->ncols; j--; )
290 setcolor(scan[j], rmx_lval(rm,i,j,0),
291 rmx_lval(rm,i,j,1),
292 rmx_lval(rm,i,j,2) );
293 if (fwritescan(scan, rm->ncols, fp) < 0) {
294 free(scan);
295 return(0);
296 }
297 }
298 free(scan);
299 return(1);
300 }
301
302 /* Write matrix to file type indicated by dt */
303 long
304 rmx_write(const RMATRIX *rm, int dtype, FILE *fp)
305 {
306 RMATRIX *mydm = NULL;
307 int ok = 1;
308
309 if ((rm == NULL) | (fp == NULL))
310 return(0);
311 /* complete header */
312 if ((dtype != DTrgbe) & (dtype != DTxyze)) {
313 fprintf(fp, "NROWS=%d\n", rm->nrows);
314 fprintf(fp, "NCOLS=%d\n", rm->ncols);
315 fprintf(fp, "NCOMP=%d\n", rm->ncomp);
316 } else if (rm->ncomp != 3) { /* wrong # components? */
317 double cmtx[3];
318 if (rm->ncomp != 1) /* only convert grayscale */
319 return(0);
320 cmtx[0] = cmtx[1] = cmtx[2] = 1;
321 mydm = rmx_transform(rm, 3, cmtx);
322 if (mydm == NULL)
323 return(0);
324 rm = mydm;
325 }
326 fputformat((char *)cm_fmt_id[dtype], fp);
327 fputc('\n', fp);
328 switch (dtype) { /* write data */
329 case DTascii:
330 ok = rmx_write_ascii(rm, fp);
331 break;
332 case DTfloat:
333 ok = rmx_write_float(rm, fp);
334 break;
335 case DTdouble:
336 ok = rmx_write_double(rm, fp);
337 break;
338 case DTrgbe:
339 case DTxyze:
340 fprtresolu(rm->ncols, rm->nrows, fp);
341 ok = rmx_write_rgbe(rm, fp);
342 break;
343 default:
344 return(0);
345 }
346 ok &= (fflush(fp) == 0);
347 rmx_free(mydm);
348 return(ftell(fp) * ok); /* return # bytes written */
349 }
350
351 /* Allocate and assign square identity matrix with n components */
352 RMATRIX *
353 rmx_identity(const int dim, const int n)
354 {
355 RMATRIX *rid = rmx_alloc(dim, dim, n);
356 int i;
357
358 if (rid == NULL)
359 return(NULL);
360 memset(rid->mtx, 0, sizeof(rid->mtx[0])*dim*dim);
361 for (i = dim; i--; )
362 rmx_lval(rid,i,i,0) = 1;
363 for (i = n; --i; )
364 memcpy(rid->mtx+i*(dim*dim), rid->mtx,
365 sizeof(rid->mtx[0])*dim*dim);
366 return(rid);
367 }
368
369 /* Duplicate the given matrix */
370 RMATRIX *
371 rmx_copy(const RMATRIX *rm)
372 {
373 RMATRIX *dnew;
374
375 if (rm == NULL)
376 return(NULL);
377 dnew = rmx_alloc(rm->nrows, rm->ncols, rm->ncomp);
378 if (dnew == NULL)
379 return(NULL);
380 memcpy(dnew->mtx, rm->mtx,
381 sizeof(rm->mtx[0])*rm->ncomp*rm->nrows*rm->ncols);
382 return(dnew);
383 }
384
385 /* Allocate and assign transposed matrix */
386 RMATRIX *
387 rmx_transpose(const RMATRIX *rm)
388 {
389 RMATRIX *dnew;
390 int i, j, k;
391
392 if (rm == NULL)
393 return(0);
394 dnew = rmx_alloc(rm->ncols, rm->nrows, rm->ncomp);
395 if (dnew == NULL)
396 return(NULL);
397 for (i = dnew->nrows; i--; )
398 for (j = dnew->ncols; j--; )
399 for (k = dnew->ncomp; k--; )
400 rmx_lval(dnew,i,j,k) = rmx_lval(rm,j,i,k);
401 return(dnew);
402 }
403
404 /* Multiply (concatenate) two matrices and allocate the result */
405 RMATRIX *
406 rmx_multiply(const RMATRIX *m1, const RMATRIX *m2)
407 {
408 RMATRIX *mres;
409 int i, j, k, h;
410
411 if ((m1 == NULL) | (m2 == NULL) ||
412 (m1->ncomp != m2->ncomp) | (m1->ncols != m2->nrows))
413 return(NULL);
414 mres = rmx_alloc(m1->nrows, m2->ncols, m1->ncomp);
415 if (mres == NULL)
416 return(NULL);
417 for (i = mres->nrows; i--; )
418 for (j = mres->ncols; j--; )
419 for (h = m1->ncols; h--; ) {
420 long double d = 0;
421 for (k = mres->ncomp; k--; )
422 d += (long double)rmx_lval(m1,i,h,k) *
423 (long double)rmx_lval(m2,h,j,k);
424 rmx_lval(mres,i,j,k) = (double)d;
425 }
426 return(mres);
427 }
428
429 /* Sum second matrix into first, applying scale factor beforehand */
430 int
431 rmx_sum(RMATRIX *msum, const RMATRIX *madd, const double sf[])
432 {
433 double *mysf = NULL;
434 int i, j, k;
435
436 if ((msum == NULL) | (madd == NULL) ||
437 (msum->nrows != madd->nrows) |
438 (msum->ncols != madd->ncols) |
439 (msum->ncomp != madd->ncomp))
440 return(0);
441 if (sf == NULL) {
442 mysf = (double *)malloc(sizeof(double)*msum->ncomp);
443 if (mysf == NULL)
444 return(0);
445 for (k = msum->ncomp; k--; )
446 mysf[k] = 1;
447 sf = mysf;
448 }
449 for (i = msum->nrows; i--; )
450 for (j = msum->ncols; j--; )
451 for (k = msum->ncomp; k--; )
452 rmx_lval(msum,i,j,k) += sf[k] * rmx_lval(madd,i,j,k);
453
454 free(mysf);
455 return(1);
456 }
457
458 /* Scale the given matrix by the indicated scalar component vector */
459 int
460 rmx_scale(RMATRIX *rm, const double sf[])
461 {
462 int i, j, k;
463
464 if ((rm == NULL) | (sf == NULL))
465 return(0);
466 for (i = rm->nrows; i--; )
467 for (j = rm->ncols; j--; )
468 for (k = rm->ncomp; k--; )
469 rmx_lval(rm,i,j,k) *= sf[k];
470
471 return(1);
472 }
473
474 /* Allocate new matrix and apply component transformation */
475 RMATRIX *
476 rmx_transform(const RMATRIX *msrc, int n, const double cmat[])
477 {
478 int i, j, ks, kd;
479 RMATRIX *dnew;
480
481 if ((msrc == NULL) | (n <= 0) | (cmat == NULL))
482 return(NULL);
483 dnew = rmx_alloc(msrc->nrows, msrc->ncols, n);
484 if (dnew == NULL)
485 return(NULL);
486 for (i = dnew->nrows; i--; )
487 for (j = dnew->ncols; j--; )
488 for (kd = dnew->ncomp; kd--; ) {
489 double d = 0;
490 for (ks = msrc->ncomp; ks--; )
491 d += cmat[kd*msrc->ncomp + ks] * rmx_lval(msrc,i,j,ks);
492 rmx_lval(dnew,i,j,kd) = d;
493 }
494 return(dnew);
495 }
496
497 /* Convert a color matrix to newly allocated RMATRIX buffer */
498 RMATRIX *
499 rmx_from_cmatrix(const CMATRIX *cm)
500 {
501 int i, j;
502 RMATRIX *dnew;
503
504 if (cm == NULL)
505 return(NULL);
506 dnew = rmx_alloc(cm->nrows, cm->ncols, 3);
507 if (dnew == NULL)
508 return(NULL);
509 for (i = dnew->nrows; i--; )
510 for (j = dnew->ncols; j--; ) {
511 const COLORV *cv = cm_lval(cm,i,j);
512 rmx_lval(dnew,i,j,0) = cv[0];
513 rmx_lval(dnew,i,j,1) = cv[1];
514 rmx_lval(dnew,i,j,2) = cv[2];
515 }
516 return(dnew);
517 }
518
519 /* Convert general matrix to newly allocated CMATRIX buffer */
520 CMATRIX *
521 cm_from_rmatrix(const RMATRIX *rm)
522 {
523 int i, j;
524 CMATRIX *cnew;
525
526 if (rm == NULL || rm->ncomp != 3)
527 return(NULL);
528 cnew = cm_alloc(rm->nrows, rm->ncols);
529 if (cnew == NULL)
530 return(NULL);
531 for (i = cnew->nrows; i--; )
532 for (j = cnew->ncols; j--; ) {
533 COLORV *cv = cm_lval(cnew,i,j);
534 cv[0] = (COLORV)rmx_lval(rm,i,j,0);
535 cv[1] = (COLORV)rmx_lval(rm,i,j,1);
536 cv[2] = (COLORV)rmx_lval(rm,i,j,2);
537 }
538 return(cnew);
539 }