ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmatrix.c
Revision: 2.5
Committed: Fri Aug 1 23:37:24 2014 UTC (9 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.4: +47 -3 lines
Log Message:
Maintain header information in matrix structure

File Contents

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