ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmatrix.c
Revision: 2.97
Committed: Fri Apr 18 23:59:03 2025 UTC (13 days, 12 hours ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.96: +3 -16 lines
Log Message:
refactor(rmtxop,rcomb,pvsum): Removed BSDF library dependency where unneeded

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: rmatrix.c,v 2.96 2025/04/17 23:30:16 greg Exp $";
3 #endif
4 /*
5 * General matrix operations.
6 */
7
8 #include <stdlib.h>
9 #include <errno.h>
10 #include "rtio.h"
11 #include "platform.h"
12 #include "resolu.h"
13 #include "paths.h"
14 #include "rmatrix.h"
15 #if !defined(_WIN32) && !defined(_WIN64)
16 #include <sys/mman.h>
17 #endif
18
19 static const char rmx_mismatch_warn[] = "WARNING: data type mismatch\n";
20
21 /* Initialize a RMATRIX struct but don't allocate array space */
22 RMATRIX *
23 rmx_new(int nr, int nc, int n)
24 {
25 RMATRIX *dnew;
26
27 if (n <= 0)
28 return(NULL);
29
30 dnew = (RMATRIX *)calloc(1, sizeof(RMATRIX));
31 if (!dnew)
32 return(NULL);
33
34 dnew->dtype = DTrmx_native;
35 dnew->nrows = nr;
36 dnew->ncols = nc;
37 dnew->ncomp = n;
38 setcolor(dnew->cexp, 1.f, 1.f, 1.f);
39 memcpy(dnew->wlpart, WLPART, sizeof(dnew->wlpart));
40
41 return(dnew);
42 }
43
44 /* Prepare a RMATRIX for writing (allocate array if needed) */
45 int
46 rmx_prepare(RMATRIX *rm)
47 {
48 if (!rm) return(0);
49 if (rm->mtx) /* assume it's right size */
50 return(1);
51 if ((rm->nrows <= 0) | (rm->ncols <= 0) | (rm->ncomp <= 0))
52 return(0);
53 rm->mtx = (rmx_dtype *)malloc(rmx_array_size(rm));
54 rm->pflags |= RMF_FREEMEM;
55 return(rm->mtx != NULL);
56 }
57
58 /* Call rmx_new() and rmx_prepare() */
59 RMATRIX *
60 rmx_alloc(int nr, int nc, int n)
61 {
62 RMATRIX *dnew = rmx_new(nr, nc, n);
63
64 if (!rmx_prepare(dnew)) {
65 rmx_free(dnew);
66 return(NULL);
67 }
68 return(dnew);
69 }
70
71 /* Clear state by freeing info and matrix data */
72 void
73 rmx_reset(RMATRIX *rm)
74 {
75 if (!rm) return;
76 if (rm->info) {
77 free(rm->info);
78 rm->info = NULL;
79 }
80 #ifdef MAP_FILE
81 if (rm->mapped) {
82 munmap(rm->mapped, rmx_mapped_size(rm));
83 rm->mapped = NULL;
84 } else
85 #endif
86 if (rm->pflags & RMF_FREEMEM) {
87 free(rm->mtx);
88 rm->pflags &= ~RMF_FREEMEM;
89 }
90 rm->mtx = NULL;
91 }
92
93 /* Free an RMATRIX struct and data */
94 void
95 rmx_free(RMATRIX *rm)
96 {
97 if (!rm) return;
98 rmx_reset(rm);
99 free(rm);
100 }
101
102 /* Resolve data type based on two input types (returns 0 for mismatch) */
103 int
104 rmx_newtype(int dtyp1, int dtyp2)
105 {
106 if ((dtyp1==DTxyze) | (dtyp1==DTrgbe) | (dtyp1==DTspec) |
107 (dtyp2==DTxyze) | (dtyp2==DTrgbe) | (dtyp2==DTspec)
108 && dtyp1 != dtyp2)
109 return(0);
110 if (dtyp1 < dtyp2)
111 return(dtyp1);
112 return(dtyp2);
113 }
114
115 /* Append header information associated with matrix data */
116 int
117 rmx_addinfo(RMATRIX *rm, const char *info)
118 {
119 size_t oldlen = 0;
120
121 if (!rm || !info || !*info)
122 return(0);
123 if (!rm->info) {
124 rm->info = (char *)malloc(strlen(info)+1);
125 } else {
126 oldlen = strlen(rm->info);
127 rm->info = (char *)realloc(rm->info,
128 oldlen+strlen(info)+1);
129 }
130 if (!rm->info)
131 return(0);
132 strcpy(rm->info+oldlen, info);
133 return(1);
134 }
135
136 static int
137 get_dminfo(char *s, void *p)
138 {
139 RMATRIX *ip = (RMATRIX *)p;
140 char fmt[MAXFMTLEN];
141 int i;
142
143 if (isheadid(s))
144 return(0);
145 if (isncomp(s)) {
146 ip->ncomp = ncompval(s);
147 return(ip->ncomp - 1);
148 }
149 if (!strncmp(s, "NROWS=", 6)) {
150 ip->nrows = atoi(s+6);
151 return(ip->nrows - 1);
152 }
153 if (!strncmp(s, "NCOLS=", 6)) {
154 ip->ncols = atoi(s+6);
155 return(ip->ncols - 1);
156 }
157 if ((i = isbigendian(s)) >= 0) {
158 if (nativebigendian() != i)
159 ip->pflags |= RMF_SWAPIN;
160 else
161 ip->pflags &= ~RMF_SWAPIN;
162 return(0);
163 }
164 if (isexpos(s)) {
165 float f = exposval(s);
166 scalecolor(ip->cexp, f);
167 return(f > .0 ? 0 : -1);
168 }
169 if (iscolcor(s)) {
170 COLOR ctmp;
171 if (!colcorval(ctmp, s)) return(-1);
172 multcolor(ip->cexp, ctmp);
173 return(0);
174 }
175 if (iswlsplit(s))
176 return(wlsplitval(ip->wlpart, s) - 1);
177
178 if (!formatval(fmt, s)) {
179 rmx_addinfo(ip, s);
180 return(0);
181 } /* else check format */
182 for (i = 1; i < DTend; i++)
183 if (!strcmp(fmt, cm_fmt_id[i])) {
184 ip->dtype = i;
185 return(0);
186 }
187 return(-1); /* bad format */
188 }
189
190 static int
191 rmx_load_ascii(rmx_dtype *drp, const RMATRIX *rm, FILE *fp)
192 {
193 int j, k;
194
195 for (j = 0; j < rm->ncols; j++)
196 for (k = rm->ncomp; k-- > 0; )
197 if (fscanf(fp, rmx_scanfmt, drp++) != 1)
198 return(0);
199 return(1);
200 }
201
202 static int
203 rmx_load_float(rmx_dtype *drp, const RMATRIX *rm, FILE *fp)
204 {
205 #if DTrmx_native==DTfloat
206 if (getbinary(drp, sizeof(*drp)*rm->ncomp, rm->ncols, fp) != rm->ncols)
207 return(0);
208 if (rm->pflags & RMF_SWAPIN)
209 swap32((char *)drp, rm->ncols*rm->ncomp);
210 #else
211 int j, k;
212 float val[MAXCOMP];
213
214 if (rm->ncomp > MAXCOMP) {
215 fputs("Unsupported # components in rmx_load_float()\n", stderr);
216 exit(1);
217 }
218 for (j = 0; j < rm->ncols; j++) {
219 if (getbinary(val, sizeof(val[0]), rm->ncomp, fp) != rm->ncomp)
220 return(0);
221 if (rm->pflags & RMF_SWAPIN)
222 swap32((char *)val, rm->ncomp);
223 for (k = 0; k < rm->ncomp; k++)
224 *drp++ = val[k];
225 }
226 #endif
227 return(1);
228 }
229
230 static int
231 rmx_load_double(rmx_dtype *drp, const RMATRIX *rm, FILE *fp)
232 {
233 #if DTrmx_native==DTdouble
234 if (getbinary(drp, sizeof(*drp)*rm->ncomp, rm->ncols, fp) != rm->ncols)
235 return(0);
236 if (rm->pflags & RMF_SWAPIN)
237 swap64((char *)drp, rm->ncols*rm->ncomp);
238 #else
239 int j, k;
240 double val[MAXCOMP];
241
242 if (rm->ncomp > MAXCOMP) {
243 fputs("Unsupported # components in rmx_load_double()\n", stderr);
244 exit(1);
245 }
246 for (j = 0; j < rm->ncols; j++) {
247 if (getbinary(val, sizeof(val[0]), rm->ncomp, fp) != rm->ncomp)
248 return(0);
249 if (rm->pflags & RMF_SWAPIN)
250 swap64((char *)val, rm->ncomp);
251 for (k = 0; k < rm->ncomp; k++)
252 *drp++ = (float)val[k];
253 }
254 #endif
255 return(1);
256 }
257
258 static int
259 rmx_load_rgbe(rmx_dtype *drp, const RMATRIX *rm, FILE *fp)
260 {
261 COLR *scan;
262 COLOR col;
263 int j;
264
265 if (rm->ncomp != 3)
266 return(0);
267 scan = (COLR *)tempbuffer(sizeof(COLR)*rm->ncols);
268 if (!scan)
269 return(0);
270 if (freadcolrs(scan, rm->ncols, fp) < 0)
271 return(0);
272 for (j = 0; j < rm->ncols; j++) {
273 colr_color(col, scan[j]);
274 *drp++ = colval(col,RED);
275 *drp++ = colval(col,GRN);
276 *drp++ = colval(col,BLU);
277 }
278 return(1);
279 }
280
281 static int
282 rmx_load_spec(rmx_dtype *drp, const RMATRIX *rm, FILE *fp)
283 {
284 COLRV *scan;
285 COLORV scol[MAXCOMP];
286 int j, k;
287
288 if ((rm->ncomp < 3) | (rm->ncomp > MAXCOMP))
289 return(0);
290 scan = (COLRV *)tempbuffer((rm->ncomp+1)*rm->ncols);
291 if (!scan)
292 return(0);
293 if (freadscolrs(scan, rm->ncomp, rm->ncols, fp) < 0)
294 return(0);
295 for (j = 0; j < rm->ncols; j++) {
296 scolr2scolor(scol, scan+j*(rm->ncomp+1), rm->ncomp);
297 for (k = 0; k < rm->ncomp; k++)
298 *drp++ = scol[k];
299 }
300 return(1);
301 }
302
303 /* Read matrix header from input stream (cannot be XML) */
304 int
305 rmx_load_header(RMATRIX *rm, FILE *fp)
306 {
307 if (!rm | !fp)
308 return(0);
309 rmx_reset(rm); /* clear state */
310 if (rm->nrows | rm->ncols | !rm->dtype) {
311 rm->nrows = rm->ncols = 0;
312 rm->ncomp = 3;
313 setcolor(rm->cexp, 1.f, 1.f, 1.f);
314 memcpy(rm->wlpart, WLPART, sizeof(rm->wlpart));
315 rm->pflags = 0;
316 }
317 rm->dtype = DTascii; /* assumed w/o FORMAT */
318 if (getheader(fp, get_dminfo, rm) < 0) {
319 fputs("Bad matrix header\n", stderr);
320 return(0);
321 }
322 if ((rm->dtype == DTrgbe) | (rm->dtype == DTxyze) &&
323 rm->ncomp != 3)
324 return(0);
325 if (rm->ncols <= 0 && /* resolution string? */
326 !fscnresolu(&rm->ncols, &rm->nrows, fp))
327 return(0);
328 if (rm->dtype == DTascii) /* set file type (WINDOWS) */
329 SET_FILE_TEXT(fp);
330 else
331 SET_FILE_BINARY(fp);
332 return(1);
333 }
334
335 /* Load next row as rmx_dtype (cannot be XML) */
336 int
337 rmx_load_row(rmx_dtype *drp, const RMATRIX *rm, FILE *fp)
338 {
339 switch (rm->dtype) {
340 case DTascii:
341 return(rmx_load_ascii(drp, rm, fp));
342 case DTfloat:
343 return(rmx_load_float(drp, rm, fp));
344 case DTdouble:
345 return(rmx_load_double(drp, rm, fp));
346 case DTrgbe:
347 case DTxyze:
348 return(rmx_load_rgbe(drp, rm, fp));
349 case DTspec:
350 return(rmx_load_spec(drp, rm, fp));
351 default:
352 fputs("Unsupported data type in rmx_load_row()\n", stderr);
353 }
354 return(0);
355 }
356
357 /* Allocate & load post-header data from stream given type set in rm->dtype */
358 int
359 rmx_load_data(RMATRIX *rm, FILE *fp)
360 {
361 int i;
362 #ifdef MAP_FILE
363 long pos; /* map memory for file > 1MB if possible */
364 if ((rm->dtype == DTrmx_native) & !(rm->pflags & RMF_SWAPIN) &
365 (rmx_array_size(rm) >= 1L<<20) &&
366 (pos = ftell(fp)) >= 0 && !(pos % sizeof(rmx_dtype))) {
367 rm->mapped = mmap(NULL, rmx_array_size(rm)+pos, PROT_READ|PROT_WRITE,
368 MAP_PRIVATE, fileno(fp), 0);
369 if (rm->mapped != MAP_FAILED) {
370 if (rm->pflags & RMF_FREEMEM)
371 free(rm->mtx);
372 rm->mtx = (rmx_dtype *)rm->mapped + pos/sizeof(rmx_dtype);
373 rm->pflags &= ~RMF_FREEMEM;
374 return(1);
375 } /* else fall back on reading into memory */
376 rm->mapped = NULL;
377 }
378 #endif
379 if (!rmx_prepare(rm)) { /* need in-core matrix array */
380 fprintf(stderr, "Cannot allocate %g MByte matrix array\n",
381 (1./(1L<<20))*(double)rmx_array_size(rm));
382 return(0);
383 }
384 for (i = 0; i < rm->nrows; i++)
385 if (!rmx_load_row(rmx_lval(rm,i,0), rm, fp))
386 return(0);
387 return(1);
388 }
389
390 /* Load matrix from supported file type */
391 RMATRIX *
392 rmx_load(const char *inspec)
393 {
394 FILE *fp;
395 RMATRIX *dnew;
396 int ok;
397
398 if (!inspec)
399 inspec = stdin_name;
400 else if (!*inspec)
401 return(NULL);
402 if (inspec == stdin_name) /* reading from stdin? */
403 fp = stdin;
404 else if (inspec[0] == '!')
405 fp = popen(inspec+1, "r");
406 else
407 fp = fopen(inspec, "r");
408 if (!fp) {
409 fprintf(stderr, "Cannot open for reading: %s\n", inspec);
410 return(NULL);
411 }
412 #ifdef getc_unlocked
413 flockfile(fp);
414 #endif
415 SET_FILE_BINARY(fp); /* load header info */
416 if (!rmx_load_header(dnew = rmx_new(0,0,3), fp)) {
417 fprintf(stderr, "Bad header in: %s\n", inspec);
418 if (inspec[0] == '!') pclose(fp);
419 else fclose(fp);
420 rmx_free(dnew);
421 return(NULL);
422 }
423 ok = rmx_load_data(dnew, fp); /* allocate & load data */
424
425 if (fp != stdin) { /* close input stream */
426 if (inspec[0] == '!')
427 pclose(fp);
428 else
429 fclose(fp);
430 }
431 #ifdef getc_unlocked
432 else
433 funlockfile(fp);
434 #endif
435 if (!ok) { /* load failure? */
436 fprintf(stderr, "Error loading data from: %s\n", inspec);
437 rmx_free(dnew);
438 return(NULL);
439 }
440 /* undo exposure? */
441 if ((dnew->cexp[0] != 1.f) |
442 (dnew->cexp[1] != 1.f) | (dnew->cexp[2] != 1.f)) {
443 double cmlt[MAXCOMP];
444 int i;
445 if (dnew->ncomp > MAXCOMP) {
446 fprintf(stderr, "Excess spectral components in: %s\n",
447 inspec);
448 rmx_free(dnew);
449 return(NULL);
450 }
451 cmlt[0] = 1./dnew->cexp[0];
452 cmlt[1] = 1./dnew->cexp[1];
453 cmlt[2] = 1./dnew->cexp[2];
454 for (i = dnew->ncomp; i-- > 3; )
455 cmlt[i] = cmlt[1]; /* XXX hack! */
456 rmx_scale(dnew, cmlt);
457 setcolor(dnew->cexp, 1.f, 1.f, 1.f);
458 }
459 return(dnew);
460 }
461
462 #if DTrmx_native==DTdouble
463 static int
464 rmx_write_float(const rmx_dtype *dp, int len, FILE *fp)
465 {
466 float val;
467
468 while (len--) {
469 val = (float)*dp++;
470 if (putbinary(&val, sizeof(val), 1, fp) != 1)
471 return(0);
472 }
473 return(1);
474 }
475 #else
476 static int
477 rmx_write_double(const rmx_dtype *dp, int len, FILE *fp)
478 {
479 double val;
480
481 while (len--) {
482 val = *dp++;
483 if (putbinary(&val, sizeof(val), 1, fp) != 1)
484 return(0);
485 }
486 return(1);
487 }
488 #endif
489
490 static int
491 rmx_write_ascii(const rmx_dtype *dp, int nc, int len, FILE *fp)
492 {
493 while (len-- > 0) {
494 int k = nc;
495 while (k-- > 0)
496 fprintf(fp, " %.7e", *dp++);
497 fputc('\t', fp);
498 }
499 return(fputc('\n', fp) != EOF);
500 }
501
502 static int
503 rmx_write_rgbe(const rmx_dtype *dp, int nc, int len, FILE *fp)
504 {
505 COLR *scan;
506 int j;
507
508 if ((nc != 1) & (nc != 3)) return(0);
509 scan = (COLR *)tempbuffer(sizeof(COLR)*len);
510 if (!scan) return(0);
511
512 for (j = 0; j < len; j++, dp += nc)
513 if (nc == 1)
514 setcolr(scan[j], dp[0], dp[0], dp[0]);
515 else
516 setcolr(scan[j], dp[0], dp[1], dp[2]);
517
518 return(fwritecolrs(scan, len, fp) >= 0);
519 }
520
521 static int
522 rmx_write_spec(const rmx_dtype *dp, int nc, int len, FILE *fp)
523 {
524 COLRV *scan;
525 COLORV scol[MAXCOMP];
526 int j, k;
527
528 if ((nc < 3) | (nc > MAXCOMP)) return(0);
529 scan = (COLRV *)tempbuffer((nc+1)*len);
530 if (!scan) return(0);
531 for (j = 0; j < len; j++, dp += nc) {
532 for (k = nc; k--; )
533 scol[k] = dp[k];
534 scolor2scolr(scan+j*(nc+1), scol, nc);
535 }
536 return(fwritescolrs(scan, nc, len, fp) >= 0);
537 }
538
539 /* Check if CIE XYZ primaries were specified */
540 static int
541 findCIEprims(const char *info)
542 {
543 RGBPRIMS prims;
544
545 if (!info)
546 return(0);
547 info = strstr(info, PRIMARYSTR);
548 if (!info || !primsval(prims, info))
549 return(0);
550
551 return((prims[RED][CIEX] > .99) & (prims[RED][CIEY] < .01) &&
552 (prims[GRN][CIEX] < .01) & (prims[GRN][CIEY] > .99) &&
553 (prims[BLU][CIEX] < .01) & (prims[BLU][CIEY] < .01));
554 }
555
556 /* Finish writing header data with resolution and format, returning type used */
557 int
558 rmx_write_header(const RMATRIX *rm, int dtype, FILE *fp)
559 {
560 if (!rm | !fp || rm->ncols <= 0)
561 return(0);
562 if (rm->info)
563 fputs(rm->info, fp);
564 if (dtype == DTfromHeader) {
565 dtype = rm->dtype;
566 #if DTrmx_native==DTfloat
567 if (dtype == DTdouble) /* but stored as float? */
568 dtype = DTfloat;
569 #endif
570 } else if (dtype == DTrgbe && (rm->dtype == DTxyze ||
571 findCIEprims(rm->info)))
572 dtype = DTxyze;
573 else if ((dtype == DTxyze) & (rm->dtype == DTrgbe))
574 dtype = DTrgbe;
575 if ((dtype < DTspec) & (rm->ncomp > 3))
576 dtype = DTspec;
577 else if ((dtype == DTspec) & (rm->ncomp <= 3))
578 return(0);
579
580 if (dtype == DTascii) /* set file type (WINDOWS) */
581 SET_FILE_TEXT(fp);
582 else
583 SET_FILE_BINARY(fp);
584 /* write exposure? */
585 if (rm->ncomp == 3 && (rm->cexp[RED] != rm->cexp[GRN]) |
586 (rm->cexp[GRN] != rm->cexp[BLU]))
587 fputcolcor(rm->cexp, fp);
588 else if (rm->cexp[GRN] != 1.f)
589 fputexpos(rm->cexp[GRN], fp);
590 /* matrix size? */
591 if ((dtype > DTspec) | (rm->nrows <= 0)) {
592 if (rm->nrows > 0)
593 fprintf(fp, "NROWS=%d\n", rm->nrows);
594 fprintf(fp, "NCOLS=%d\n", rm->ncols);
595 }
596 if (dtype >= DTspec) { /* # components & split? */
597 fputncomp(rm->ncomp, fp);
598 if (rm->ncomp > 3 &&
599 memcmp(rm->wlpart, WLPART, sizeof(WLPART)))
600 fputwlsplit(rm->wlpart, fp);
601 } else if ((rm->ncomp != 3) & (rm->ncomp != 1))
602 return(0); /* wrong # components */
603 if ((dtype == DTfloat) | (dtype == DTdouble))
604 fputendian(fp); /* important to record */
605 fputformat(cm_fmt_id[dtype], fp);
606 fputc('\n', fp); /* end of header */
607 if ((dtype <= DTspec) & (rm->nrows > 0))
608 fprtresolu(rm->ncols, rm->nrows, fp);
609 return(dtype);
610 }
611
612 /* Write out matrix data (usually by row) */
613 int
614 rmx_write_data(const rmx_dtype *dp, int nc, int len, int dtype, FILE *fp)
615 {
616 switch (dtype) {
617 #if DTrmx_native==DTdouble
618 case DTfloat:
619 return(rmx_write_float(dp, nc*len, fp));
620 #else
621 case DTdouble:
622 return(rmx_write_double(dp, nc*len, fp));
623 #endif
624 case DTrmx_native:
625 return(putbinary(dp, sizeof(*dp)*nc, len, fp) == len);
626 case DTascii:
627 return(rmx_write_ascii(dp, nc, len, fp));
628 case DTrgbe:
629 case DTxyze:
630 return(rmx_write_rgbe(dp, nc, len, fp));
631 case DTspec:
632 return(rmx_write_spec(dp, nc, len, fp));
633 }
634 return(0);
635 }
636
637 /* Write matrix using file format indicated by dtype */
638 int
639 rmx_write(const RMATRIX *rm, int dtype, FILE *fp)
640 {
641 int ok = 0;
642 int i;
643 /* complete header */
644 dtype = rmx_write_header(rm, dtype, fp);
645 if (dtype <= 0)
646 return(0);
647 #ifdef getc_unlocked
648 flockfile(fp);
649 #endif
650 if (dtype == DTrmx_native) /* write all at once? */
651 ok = rmx_write_data(rm->mtx, rm->ncomp,
652 rm->nrows*rm->ncols, dtype, fp);
653 else /* else row by row */
654 for (i = 0; i < rm->nrows; i++) {
655 ok = rmx_write_data(rmx_val(rm,i,0), rm->ncomp,
656 rm->ncols, dtype, fp);
657 if (!ok) break;
658 }
659
660 if (ok) ok = (fflush(fp) == 0);
661 #ifdef getc_unlocked
662 funlockfile(fp);
663 #endif
664 if (!ok) fputs("Error writing matrix\n", stderr);
665 return(ok);
666 }
667
668 /* Allocate and assign square identity matrix with n components */
669 RMATRIX *
670 rmx_identity(const int dim, const int n)
671 {
672 RMATRIX *rid = rmx_alloc(dim, dim, n);
673 int i, k;
674
675 if (!rid)
676 return(NULL);
677 memset(rid->mtx, 0, rmx_array_size(rid));
678 for (i = dim; i--; ) {
679 rmx_dtype *dp = rmx_lval(rid,i,i);
680 for (k = n; k--; )
681 dp[k] = 1.;
682 }
683 return(rid);
684 }
685
686 /* Duplicate the given matrix (may be unallocated) */
687 RMATRIX *
688 rmx_copy(const RMATRIX *rm)
689 {
690 RMATRIX *dnew;
691
692 if (!rm)
693 return(NULL);
694 dnew = rmx_new(rm->nrows, rm->ncols, rm->ncomp);
695 if (!dnew)
696 return(NULL);
697 if (rm->mtx) {
698 if (!rmx_prepare(dnew)) {
699 rmx_free(dnew);
700 return(NULL);
701 }
702 memcpy(dnew->mtx, rm->mtx, rmx_array_size(dnew));
703 }
704 rmx_addinfo(dnew, rm->info);
705 dnew->dtype = rm->dtype;
706 copycolor(dnew->cexp, rm->cexp);
707 memcpy(dnew->wlpart, rm->wlpart, sizeof(dnew->wlpart));
708 return(dnew);
709 }
710
711 /* Replace data in first matrix with data from second */
712 int
713 rmx_transfer_data(RMATRIX *rdst, RMATRIX *rsrc, int dometa)
714 {
715 if (!rdst | !rsrc)
716 return(0);
717 if (dometa) { /* transfer everything? */
718 rmx_reset(rdst);
719 *rdst = *rsrc;
720 rsrc->info = NULL; rsrc->mapped = NULL; rsrc->mtx = NULL;
721 return(1);
722 }
723 /* just matrix data -- leave metadata */
724 if ((rdst->nrows != rsrc->nrows) |
725 (rdst->ncols != rsrc->ncols) |
726 (rdst->ncomp != rsrc->ncomp))
727 return(0);
728 #ifdef MAP_FILE
729 if (rdst->mapped)
730 munmap(rdst->mapped, rmx_mapped_size(rdst));
731 else
732 #endif
733 if (rdst->pflags & RMF_FREEMEM) {
734 free(rdst->mtx);
735 rdst->pflags &= ~RMF_FREEMEM;
736 }
737 rdst->mapped = rsrc->mapped;
738 rdst->mtx = rsrc->mtx;
739 rdst->pflags |= rsrc->pflags & RMF_FREEMEM;
740 rsrc->mapped = NULL; rsrc->mtx = NULL;
741 return(1);
742 }
743
744 /* Transpose the given matrix */
745 int
746 rmx_transpose(RMATRIX *rm)
747 {
748 uby8 *bmap;
749 rmx_dtype val[MAXCOMP];
750 RMATRIX dold;
751 int i, j;
752
753 if (!rm || !rm->mtx | (rm->ncomp > MAXCOMP))
754 return(0);
755 if (rm->info)
756 rmx_addinfo(rm, "Transposed rows and columns\n");
757 if ((rm->nrows == 1) | (rm->ncols == 1)) { /* vector? */
758 j = rm->ncols;
759 rm->ncols = rm->nrows;
760 rm->nrows = j;
761 return(1);
762 }
763 if (rm->nrows == rm->ncols) { /* square matrix case */
764 for (i = rm->nrows; i--; )
765 for (j = rm->ncols; j--; ) {
766 if (i == j) continue;
767 memcpy(val, rmx_val(rm,i,j),
768 sizeof(rmx_dtype)*rm->ncomp);
769 memcpy(rmx_lval(rm,i,j), rmx_val(rm,j,i),
770 sizeof(rmx_dtype)*rm->ncomp);
771 memcpy(rmx_val(rm,j,i), val,
772 sizeof(rmx_dtype)*rm->ncomp);
773 }
774 return(1);
775 }
776 #define bmbyte(r,c) bmap[((r)*rm->ncols+(c))>>3]
777 #define bmbit(r,c) (1 << ((r)*rm->ncols+(c) & 7))
778 #define bmop(r,c, op) (bmbyte(r,c) op bmbit(r,c))
779 #define bmtest(r,c) bmop(r,c,&)
780 #define bmset(r,c) bmop(r,c,|=)
781 /* loop completion bitmap */
782 bmap = (uby8 *)calloc(((size_t)rm->nrows*rm->ncols+7)>>3, 1);
783 if (!bmap)
784 return(0);
785 dold = *rm;
786 rm->ncols = dold.nrows; rm->nrows = dold.ncols;
787 for (i = rm->nrows; i--; ) /* try every starting point */
788 for (j = rm->ncols; j--; ) {
789 int i0, j0;
790 int i1 = i;
791 size_t j1 = j;
792 if (bmtest(i, j))
793 continue; /* traversed loop earlier */
794 memcpy(val, rmx_val(rm,i,j),
795 sizeof(rmx_dtype)*rm->ncomp);
796 for ( ; ; ) { /* new transpose loop */
797 const rmx_dtype *ds;
798 i0 = i1; j0 = j1;
799 ds = rmx_val(&dold, j0, i0);
800 j1 = (ds - dold.mtx)/dold.ncomp;
801 i1 = j1 / rm->ncols;
802 j1 -= (size_t)i1*rm->ncols;
803 bmset(i1, j1); /* mark as done */
804 if ((i1 == i) & (j1 == j))
805 break; /* back at start */
806 memcpy(rmx_lval(rm,i0,j0), ds,
807 sizeof(rmx_dtype)*rm->ncomp);
808 } /* complete the loop */
809 memcpy(rmx_lval(rm,i0,j0), val,
810 sizeof(rmx_dtype)*rm->ncomp);
811 }
812 free(bmap); /* all done! */
813 return(1);
814 #undef bmbyte
815 #undef bmbit
816 #undef bmop
817 #undef bmtest
818 #undef bmset
819 }
820
821 /* Multiply (concatenate) two matrices and allocate the result */
822 RMATRIX *
823 rmx_multiply(const RMATRIX *m1, const RMATRIX *m2)
824 {
825 RMATRIX *mres;
826 int i, j, k, h;
827
828 if (!m1 | !m2 || !m1->mtx | !m2->mtx |
829 (m1->ncomp != m2->ncomp) | (m1->ncols != m2->nrows))
830 return(NULL);
831 mres = rmx_alloc(m1->nrows, m2->ncols, m1->ncomp);
832 if (!mres)
833 return(NULL);
834 i = rmx_newtype(m1->dtype, m2->dtype);
835 if (i)
836 mres->dtype = i;
837 else
838 rmx_addinfo(mres, rmx_mismatch_warn);
839 for (i = mres->nrows; i--; )
840 for (j = mres->ncols; j--; )
841 for (k = mres->ncomp; k--; ) {
842 double d = 0;
843 for (h = m1->ncols; h--; )
844 d += (double)rmx_val(m1,i,h)[k] *
845 rmx_val(m2,h,j)[k];
846 rmx_lval(mres,i,j)[k] = (rmx_dtype)d;
847 }
848 return(mres);
849 }
850
851 /* Element-wise multiplication (or division) of m2 into m1 */
852 int
853 rmx_elemult(RMATRIX *m1, const RMATRIX *m2, int divide)
854 {
855 int zeroDivides = 0;
856 int i, j, k;
857
858 if (!m1 | !m2 || !m1->mtx | !m2->mtx |
859 (m1->ncols != m2->ncols) | (m1->nrows != m2->nrows))
860 return(0);
861 if ((m2->ncomp > 1) & (m2->ncomp != m1->ncomp))
862 return(0);
863 i = rmx_newtype(m1->dtype, m2->dtype);
864 if (i)
865 m1->dtype = i;
866 else
867 rmx_addinfo(m1, rmx_mismatch_warn);
868 for (i = m1->nrows; i--; )
869 for (j = m1->ncols; j--; )
870 if (divide) {
871 rmx_dtype d;
872 if (m2->ncomp == 1) {
873 d = rmx_val(m2,i,j)[0];
874 if (d == 0) {
875 ++zeroDivides;
876 for (k = m1->ncomp; k--; )
877 rmx_lval(m1,i,j)[k] = 0;
878 } else {
879 d = 1./d;
880 for (k = m1->ncomp; k--; )
881 rmx_lval(m1,i,j)[k] *= d;
882 }
883 } else
884 for (k = m1->ncomp; k--; ) {
885 d = rmx_val(m2,i,j)[k];
886 if (d == 0) {
887 ++zeroDivides;
888 rmx_lval(m1,i,j)[k] = 0;
889 } else
890 rmx_lval(m1,i,j)[k] /= d;
891 }
892 } else {
893 if (m2->ncomp == 1) {
894 const rmx_dtype d = rmx_val(m2,i,j)[0];
895 for (k = m1->ncomp; k--; )
896 rmx_lval(m1,i,j)[k] *= d;
897 } else
898 for (k = m1->ncomp; k--; )
899 rmx_lval(m1,i,j)[k] *= rmx_val(m2,i,j)[k];
900 }
901 if (zeroDivides) {
902 rmx_addinfo(m1, "WARNING: zero divide(s) corrupted results\n");
903 errno = ERANGE;
904 }
905 return(1);
906 }
907
908 /* Sum second matrix into first, applying scale factor beforehand */
909 int
910 rmx_sum(RMATRIX *msum, const RMATRIX *madd, const double sf[])
911 {
912 double *mysf = NULL;
913 int i, j, k;
914
915 if (!msum | !madd || !msum->mtx | !madd->mtx |
916 (msum->nrows != madd->nrows) |
917 (msum->ncols != madd->ncols) |
918 (msum->ncomp != madd->ncomp))
919 return(0);
920 if (!sf) {
921 mysf = (double *)malloc(sizeof(double)*msum->ncomp);
922 if (!mysf)
923 return(0);
924 for (k = msum->ncomp; k--; )
925 mysf[k] = 1;
926 sf = mysf;
927 }
928 i = rmx_newtype(msum->dtype, madd->dtype);
929 if (i)
930 msum->dtype = i;
931 else
932 rmx_addinfo(msum, rmx_mismatch_warn);
933 for (i = msum->nrows; i--; )
934 for (j = msum->ncols; j--; ) {
935 const rmx_dtype *da = rmx_val(madd,i,j);
936 rmx_dtype *ds = rmx_lval(msum,i,j);
937 for (k = msum->ncomp; k--; )
938 ds[k] += (rmx_dtype)sf[k] * da[k];
939 }
940 if (mysf)
941 free(mysf);
942 return(1);
943 }
944
945 /* Scale the given matrix by the indicated scalar component vector */
946 int
947 rmx_scale(RMATRIX *rm, const double sf[])
948 {
949 int i, j, k;
950
951 if (!rm | !sf || !rm->mtx)
952 return(0);
953 for (i = rm->nrows; i--; )
954 for (j = rm->ncols; j--; ) {
955 rmx_dtype *dp = rmx_lval(rm,i,j);
956 for (k = rm->ncomp; k--; )
957 dp[k] *= (rmx_dtype)sf[k];
958 }
959 if (rm->info)
960 rmx_addinfo(rm, "Applied scalar\n");
961 /* XXX: should record as exposure for COLR and SCOLR types? */
962 return(1);
963 }
964
965 /* Allocate new matrix and apply component transformation */
966 RMATRIX *
967 rmx_transform(const RMATRIX *msrc, int n, const double cmat[])
968 {
969 int i, j, ks, kd;
970 RMATRIX *dnew;
971
972 if (!msrc | (n <= 0) | !cmat || !msrc->mtx)
973 return(NULL);
974 dnew = rmx_alloc(msrc->nrows, msrc->ncols, n);
975 if (!dnew)
976 return(NULL);
977 if (msrc->info) {
978 char buf[128];
979 sprintf(buf, "Applied %dx%d component transform\n",
980 dnew->ncomp, msrc->ncomp);
981 rmx_addinfo(dnew, msrc->info);
982 rmx_addinfo(dnew, buf);
983 }
984 dnew->dtype = msrc->dtype;
985 for (i = dnew->nrows; i--; )
986 for (j = dnew->ncols; j--; ) {
987 const rmx_dtype *ds = rmx_val(msrc,i,j);
988 for (kd = dnew->ncomp; kd--; ) {
989 double d = 0;
990 for (ks = msrc->ncomp; ks--; )
991 d += cmat[kd*msrc->ncomp + ks] * ds[ks];
992 rmx_lval(dnew,i,j)[kd] = (rmx_dtype)d;
993 }
994 }
995 return(dnew);
996 }
997
998 /* Convert a color matrix to newly allocated RMATRIX buffer */
999 RMATRIX *
1000 rmx_from_cmatrix(const CMATRIX *cm)
1001 {
1002 RMATRIX *dnew;
1003
1004 if (!cm)
1005 return(NULL);
1006 dnew = rmx_alloc(cm->nrows, cm->ncols, 3);
1007 if (!dnew)
1008 return(NULL);
1009
1010 dnew->dtype = sizeof(COLORV)==sizeof(float) ?
1011 DTfloat : DTdouble;
1012
1013 if (sizeof(COLORV) == sizeof(rmx_dtype)) {
1014 memcpy(dnew->mtx, cm->cmem, rmx_array_size(dnew));
1015 } else {
1016 int i, j;
1017 for (i = dnew->nrows; i--; )
1018 for (j = dnew->ncols; j--; ) {
1019 const COLORV *cv = cm_lval(cm,i,j);
1020 rmx_dtype *dp = rmx_lval(dnew,i,j);
1021 dp[0] = cv[0];
1022 dp[1] = cv[1];
1023 dp[2] = cv[2];
1024 }
1025 }
1026 return(dnew);
1027 }
1028
1029 /* Convert general matrix to newly allocated CMATRIX buffer */
1030 CMATRIX *
1031 cm_from_rmatrix(const RMATRIX *rm)
1032 {
1033 CMATRIX *cnew;
1034
1035 if (!rm || !rm->mtx | (rm->ncomp == 2) | (rm->ncomp > MAXCOMP))
1036 return(NULL);
1037 cnew = cm_alloc(rm->nrows, rm->ncols);
1038 if (!cnew)
1039 return(NULL);
1040 if ((sizeof(COLORV) == sizeof(rmx_dtype)) & (rm->ncomp == 3)) {
1041 memcpy(cnew->cmem, rm->mtx, rmx_array_size(rm));
1042 } else {
1043 int i, j;
1044 for (i = cnew->nrows; i--; )
1045 for (j = cnew->ncols; j--; ) {
1046 const rmx_dtype *dp = rmx_val(rm,i,j);
1047 COLORV *cv = cm_lval(cnew,i,j);
1048 switch (rm->ncomp) {
1049 case 3:
1050 setcolor(cv, dp[0], dp[1], dp[2]);
1051 break;
1052 case 1:
1053 setcolor(cv, dp[0], dp[0], dp[0]);
1054 break;
1055 default:
1056 if (sizeof(COLORV) == sizeof(rmx_dtype)) {
1057 scolor2color(cv, (const COLORV *)dp,
1058 rm->ncomp, rm->wlpart);
1059 } else {
1060 COLORV scol[MAXCOMP];
1061 int k = rm->ncomp;
1062 while (k--) scol[k] = dp[k];
1063 scolor2color(cv, scol, rm->ncomp, rm->wlpart);
1064 }
1065 break;
1066 }
1067 }
1068 }
1069 return(cnew);
1070 }