ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmatrix.c
Revision: 2.16
Committed: Wed Jul 22 04:29:56 2015 UTC (8 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.15: +2 -3 lines
Log Message:
Reduced precision of multiplies in favor of time when little differ found

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: rmatrix.c,v 2.15 2015/05/08 18:25:03 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 "rtprocess.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 #ifdef _WIN32
115 _setmode(fileno(fp), _O_TEXT);
116 #endif
117 for (i = 0; i < rm->nrows; i++)
118 for (j = 0; j < rm->ncols; j++)
119 for (k = 0; k < rm->ncomp; k++)
120 if (fscanf(fp, "%lf", &rmx_lval(rm,i,j,k)) != 1)
121 return(0);
122 return(1);
123 }
124
125 static int
126 rmx_load_float(RMATRIX *rm, FILE *fp)
127 {
128 int i, j, k;
129 float val[100];
130
131 if (rm->ncomp > 100) {
132 fputs("Unsupported # components in rmx_load_float()\n", stderr);
133 exit(1);
134 }
135 for (i = 0; i < rm->nrows; i++)
136 for (j = 0; j < rm->ncols; j++) {
137 if (fread(val, sizeof(val[0]), rm->ncomp, fp) != rm->ncomp)
138 return(0);
139 for (k = rm->ncomp; k--; )
140 rmx_lval(rm,i,j,k) = val[k];
141 }
142 return(1);
143 }
144
145 static int
146 rmx_load_double(RMATRIX *rm, FILE *fp)
147 {
148 int i, j, k;
149 double val[100];
150
151 if (rm->ncomp > 100) {
152 fputs("Unsupported # components in rmx_load_double()\n", stderr);
153 exit(1);
154 }
155 for (i = 0; i < rm->nrows; i++)
156 for (j = 0; j < rm->ncols; j++) {
157 if (fread(val, sizeof(val[0]), rm->ncomp, fp) != rm->ncomp)
158 return(0);
159 for (k = rm->ncomp; k--; )
160 rmx_lval(rm,i,j,k) = val[k];
161 }
162 return(1);
163 }
164
165 static int
166 rmx_load_rgbe(RMATRIX *rm, FILE *fp)
167 {
168 COLOR *scan = (COLOR *)malloc(sizeof(COLOR)*rm->ncols);
169 int i, j;
170
171 if (scan == NULL)
172 return(0);
173 for (i = 0; i < rm->nrows; i++) {
174 if (freadscan(scan, rm->ncols, fp) < 0) {
175 free(scan);
176 return(0);
177 }
178 for (j = rm->ncols; j--; ) {
179 rmx_lval(rm,i,j,0) = colval(scan[j],RED);
180 rmx_lval(rm,i,j,1) = colval(scan[j],GRN);
181 rmx_lval(rm,i,j,2) = colval(scan[j],BLU);
182 }
183 }
184 free(scan);
185 return(1);
186 }
187
188 /* Load matrix from supported file type */
189 RMATRIX *
190 rmx_load(const char *inspec)
191 {
192 FILE *fp = stdin;
193 RMATRIX dinfo;
194 RMATRIX *dnew;
195
196 if (inspec == NULL) { /* reading from stdin? */
197 inspec = "<stdin>";
198 #ifdef _WIN32
199 _setmode(fileno(stdin), _O_BINARY);
200 #endif
201 } else if (inspec[0] == '!') {
202 if ((fp = popen(inspec+1, "r")) == NULL)
203 return(NULL);
204 #ifdef _WIN32
205 _setmode(fileno(fp), _O_BINARY);
206 #endif
207 } else {
208 const char *sp = inspec; /* check suffix */
209 while (*sp)
210 ++sp;
211 while (sp > inspec && sp[-1] != '.')
212 --sp;
213 if (!strcasecmp(sp, "XML")) { /* assume it's a BSDF */
214 CMATRIX *cm = cm_loadBTDF((char *)inspec);
215 if (cm == NULL)
216 return(NULL);
217 dnew = rmx_from_cmatrix(cm);
218 cm_free(cm);
219 return(dnew);
220 }
221 /* else open it ourselves */
222 if ((fp = fopen(inspec, "rb")) == NULL)
223 return(NULL);
224 }
225 #ifdef getc_unlocked
226 flockfile(fp);
227 #endif
228 dinfo.nrows = dinfo.ncols = dinfo.ncomp = 0;
229 dinfo.dtype = DTascii; /* assumed w/o FORMAT */
230 dinfo.info = NULL;
231 if (getheader(fp, get_dminfo, &dinfo) < 0) {
232 fclose(fp);
233 return(NULL);
234 }
235 if ((dinfo.nrows <= 0) | (dinfo.ncols <= 0)) {
236 if (!fscnresolu(&dinfo.ncols, &dinfo.nrows, fp)) {
237 fclose(fp);
238 return(NULL);
239 }
240 if (dinfo.ncomp <= 0)
241 dinfo.ncomp = 3;
242 else if ((dinfo.dtype == DTrgbe) | (dinfo.dtype == DTxyze) &&
243 dinfo.ncomp != 3) {
244 fclose(fp);
245 return(NULL);
246 }
247 }
248 dnew = rmx_alloc(dinfo.nrows, dinfo.ncols, dinfo.ncomp);
249 if (dnew == NULL) {
250 fclose(fp);
251 return(NULL);
252 }
253 dnew->info = dinfo.info;
254 switch (dinfo.dtype) {
255 case DTascii:
256 if (!rmx_load_ascii(dnew, fp))
257 goto loaderr;
258 dnew->dtype = DTascii; /* should leave double? */
259 break;
260 case DTfloat:
261 if (!rmx_load_float(dnew, fp))
262 goto loaderr;
263 dnew->dtype = DTfloat;
264 break;
265 case DTdouble:
266 if (!rmx_load_double(dnew, fp))
267 goto loaderr;
268 dnew->dtype = DTdouble;
269 break;
270 case DTrgbe:
271 case DTxyze:
272 if (!rmx_load_rgbe(dnew, fp))
273 goto loaderr;
274 dnew->dtype = dinfo.dtype;
275 break;
276 default:
277 goto loaderr;
278 }
279 if (fp != stdin) {
280 if (inspec[0] == '!')
281 pclose(fp);
282 else
283 fclose(fp);
284 }
285 #ifdef getc_unlocked
286 else
287 funlockfile(fp);
288 #endif
289 return(dnew);
290 loaderr: /* should report error? */
291 if (inspec[0] == '!')
292 pclose(fp);
293 else
294 fclose(fp);
295 rmx_free(dnew);
296 return(NULL);
297 }
298
299 static int
300 rmx_write_ascii(const RMATRIX *rm, FILE *fp)
301 {
302 int i, j, k;
303 #ifdef _WIN32
304 _setmode(fileno(fp), _O_TEXT);
305 #endif
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 /* complete header */
389 if (rm->info)
390 fputs(rm->info, fp);
391 if (dtype == DTfromHeader)
392 dtype = rm->dtype;
393 else if ((dtype == DTrgbe) & (rm->dtype == DTxyze))
394 dtype = DTxyze;
395 else if ((dtype == DTxyze) & (rm->dtype == DTrgbe))
396 dtype = DTrgbe;
397 if ((dtype != DTrgbe) & (dtype != DTxyze)) {
398 fprintf(fp, "NROWS=%d\n", rm->nrows);
399 fprintf(fp, "NCOLS=%d\n", rm->ncols);
400 fprintf(fp, "NCOMP=%d\n", rm->ncomp);
401 } else if (rm->ncomp != 3) { /* wrong # components? */
402 double cmtx[3];
403 if (rm->ncomp != 1) /* only convert grayscale */
404 return(0);
405 cmtx[0] = cmtx[1] = cmtx[2] = 1;
406 mydm = rmx_transform(rm, 3, cmtx);
407 if (mydm == NULL)
408 return(0);
409 rm = mydm;
410 }
411 fputformat((char *)cm_fmt_id[dtype], fp);
412 fputc('\n', fp);
413 switch (dtype) { /* write data */
414 case DTascii:
415 ok = rmx_write_ascii(rm, fp);
416 break;
417 case DTfloat:
418 ok = rmx_write_float(rm, fp);
419 break;
420 case DTdouble:
421 ok = rmx_write_double(rm, fp);
422 break;
423 case DTrgbe:
424 case DTxyze:
425 fprtresolu(rm->ncols, rm->nrows, fp);
426 ok = rmx_write_rgbe(rm, fp);
427 break;
428 default:
429 return(0);
430 }
431 ok &= (fflush(fp) == 0);
432 rmx_free(mydm);
433 return(ok);
434 }
435
436 /* Allocate and assign square identity matrix with n components */
437 RMATRIX *
438 rmx_identity(const int dim, const int n)
439 {
440 RMATRIX *rid = rmx_alloc(dim, dim, n);
441 int i, k;
442
443 if (rid == NULL)
444 return(NULL);
445 memset(rid->mtx, 0, sizeof(rid->mtx[0])*n*dim*dim);
446 for (i = dim; i--; )
447 for (k = n; k--; )
448 rmx_lval(rid,i,i,k) = 1;
449 return(rid);
450 }
451
452 /* Duplicate the given matrix */
453 RMATRIX *
454 rmx_copy(const RMATRIX *rm)
455 {
456 RMATRIX *dnew;
457
458 if (rm == NULL)
459 return(NULL);
460 dnew = rmx_alloc(rm->nrows, rm->ncols, rm->ncomp);
461 if (dnew == NULL)
462 return(NULL);
463 rmx_addinfo(dnew, rm->info);
464 dnew->dtype = rm->dtype;
465 memcpy(dnew->mtx, rm->mtx,
466 sizeof(rm->mtx[0])*rm->ncomp*rm->nrows*rm->ncols);
467 return(dnew);
468 }
469
470 /* Allocate and assign transposed matrix */
471 RMATRIX *
472 rmx_transpose(const RMATRIX *rm)
473 {
474 RMATRIX *dnew;
475 int i, j, k;
476
477 if (rm == NULL)
478 return(0);
479 dnew = rmx_alloc(rm->ncols, rm->nrows, rm->ncomp);
480 if (dnew == NULL)
481 return(NULL);
482 if (rm->info) {
483 rmx_addinfo(dnew, rm->info);
484 rmx_addinfo(dnew, "Transposed rows and columns\n");
485 }
486 dnew->dtype = rm->dtype;
487 for (i = dnew->nrows; i--; )
488 for (j = dnew->ncols; j--; )
489 for (k = dnew->ncomp; k--; )
490 rmx_lval(dnew,i,j,k) = rmx_lval(rm,j,i,k);
491 return(dnew);
492 }
493
494 /* Multiply (concatenate) two matrices and allocate the result */
495 RMATRIX *
496 rmx_multiply(const RMATRIX *m1, const RMATRIX *m2)
497 {
498 RMATRIX *mres;
499 int i, j, k, h;
500
501 if ((m1 == NULL) | (m2 == NULL) ||
502 (m1->ncomp != m2->ncomp) | (m1->ncols != m2->nrows))
503 return(NULL);
504 mres = rmx_alloc(m1->nrows, m2->ncols, m1->ncomp);
505 if (mres == NULL)
506 return(NULL);
507 i = rmx_newtype(m1->dtype, m2->dtype);
508 if (i)
509 mres->dtype = i;
510 else
511 rmx_addinfo(mres, rmx_mismatch_warn);
512 for (i = mres->nrows; i--; )
513 for (j = mres->ncols; j--; )
514 for (k = mres->ncomp; k--; ) {
515 long double d = 0;
516 for (h = m1->ncols; h--; )
517 d += rmx_lval(m1,i,h,k) * rmx_lval(m2,h,j,k);
518 rmx_lval(mres,i,j,k) = (double)d;
519 }
520 return(mres);
521 }
522
523 /* Sum second matrix into first, applying scale factor beforehand */
524 int
525 rmx_sum(RMATRIX *msum, const RMATRIX *madd, const double sf[])
526 {
527 double *mysf = NULL;
528 int i, j, k;
529
530 if ((msum == NULL) | (madd == NULL) ||
531 (msum->nrows != madd->nrows) |
532 (msum->ncols != madd->ncols) |
533 (msum->ncomp != madd->ncomp))
534 return(0);
535 if (sf == NULL) {
536 mysf = (double *)malloc(sizeof(double)*msum->ncomp);
537 if (mysf == NULL)
538 return(0);
539 for (k = msum->ncomp; k--; )
540 mysf[k] = 1;
541 sf = mysf;
542 }
543 i = rmx_newtype(msum->dtype, madd->dtype);
544 if (i)
545 msum->dtype = i;
546 else
547 rmx_addinfo(msum, rmx_mismatch_warn);
548 for (i = msum->nrows; i--; )
549 for (j = msum->ncols; j--; )
550 for (k = msum->ncomp; k--; )
551 rmx_lval(msum,i,j,k) += sf[k] * rmx_lval(madd,i,j,k);
552
553 free(mysf);
554 return(1);
555 }
556
557 /* Scale the given matrix by the indicated scalar component vector */
558 int
559 rmx_scale(RMATRIX *rm, const double sf[])
560 {
561 int i, j, k;
562
563 if ((rm == NULL) | (sf == NULL))
564 return(0);
565 for (i = rm->nrows; i--; )
566 for (j = rm->ncols; j--; )
567 for (k = rm->ncomp; k--; )
568 rmx_lval(rm,i,j,k) *= sf[k];
569
570 return(1);
571 }
572
573 /* Allocate new matrix and apply component transformation */
574 RMATRIX *
575 rmx_transform(const RMATRIX *msrc, int n, const double cmat[])
576 {
577 int i, j, ks, kd;
578 RMATRIX *dnew;
579
580 if ((msrc == NULL) | (n <= 0) | (cmat == NULL))
581 return(NULL);
582 dnew = rmx_alloc(msrc->nrows, msrc->ncols, n);
583 if (dnew == NULL)
584 return(NULL);
585 dnew->dtype = msrc->dtype;
586 for (i = dnew->nrows; i--; )
587 for (j = dnew->ncols; j--; )
588 for (kd = dnew->ncomp; kd--; ) {
589 double d = 0;
590 for (ks = msrc->ncomp; ks--; )
591 d += cmat[kd*msrc->ncomp + ks] * rmx_lval(msrc,i,j,ks);
592 rmx_lval(dnew,i,j,kd) = d;
593 }
594 return(dnew);
595 }
596
597 /* Convert a color matrix to newly allocated RMATRIX buffer */
598 RMATRIX *
599 rmx_from_cmatrix(const CMATRIX *cm)
600 {
601 int i, j;
602 RMATRIX *dnew;
603
604 if (cm == NULL)
605 return(NULL);
606 dnew = rmx_alloc(cm->nrows, cm->ncols, 3);
607 if (dnew == NULL)
608 return(NULL);
609 dnew->dtype = DTfloat;
610 for (i = dnew->nrows; i--; )
611 for (j = dnew->ncols; j--; ) {
612 const COLORV *cv = cm_lval(cm,i,j);
613 rmx_lval(dnew,i,j,0) = cv[0];
614 rmx_lval(dnew,i,j,1) = cv[1];
615 rmx_lval(dnew,i,j,2) = cv[2];
616 }
617 return(dnew);
618 }
619
620 /* Convert general matrix to newly allocated CMATRIX buffer */
621 CMATRIX *
622 cm_from_rmatrix(const RMATRIX *rm)
623 {
624 int i, j;
625 CMATRIX *cnew;
626
627 if (rm == NULL || rm->ncomp != 3)
628 return(NULL);
629 cnew = cm_alloc(rm->nrows, rm->ncols);
630 if (cnew == NULL)
631 return(NULL);
632 for (i = cnew->nrows; i--; )
633 for (j = cnew->ncols; j--; ) {
634 COLORV *cv = cm_lval(cnew,i,j);
635 cv[0] = (COLORV)rmx_lval(rm,i,j,0);
636 cv[1] = (COLORV)rmx_lval(rm,i,j,1);
637 cv[2] = (COLORV)rmx_lval(rm,i,j,2);
638 }
639 return(cnew);
640 }