ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmatrix.c
Revision: 2.19
Committed: Tue Feb 2 18:02:32 2016 UTC (8 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.18: +2 -2 lines
Log Message:
Moved declaration of popen to paths.h and put convert_command() into module

File Contents

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