ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmatrix.c
Revision: 2.90
Committed: Fri Apr 4 22:47:56 2025 UTC (3 weeks, 6 days ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.89: +7 -3 lines
Log Message:
fix(rmtxop,rcomb): Minor fixes to auto-typing, preferring native output

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: rmatrix.c,v 2.89 2025/04/04 18:06:48 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, RMPref rmp)
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 { /* check suffix */
407 const char *sp = strrchr(inspec, '.');
408 if (sp > inspec && !strcasecmp(sp+1, "XML")) { /* BSDF? */
409 CMATRIX *cm = rmp==RMPnone ? (CMATRIX *)NULL :
410 rmp==RMPtrans ? cm_loadBTDF(inspec) :
411 cm_loadBRDF(inspec, rmp==RMPreflB) ;
412 if (!cm)
413 return(NULL);
414 dnew = rmx_from_cmatrix(cm);
415 cm_free(cm);
416 dnew->dtype = DTascii;
417 return(dnew); /* return here */
418 } /* else open it ourselves */
419 fp = fopen(inspec, "r");
420 }
421 if (!fp) {
422 fprintf(stderr, "Cannot open for reading: %s\n", inspec);
423 return(NULL);
424 }
425 #ifdef getc_unlocked
426 flockfile(fp);
427 #endif
428 SET_FILE_BINARY(fp); /* load header info */
429 if (!rmx_load_header(dnew = rmx_new(0,0,3), fp)) {
430 fprintf(stderr, "Bad header in: %s\n", inspec);
431 if (inspec[0] == '!') pclose(fp);
432 else fclose(fp);
433 rmx_free(dnew);
434 return(NULL);
435 }
436 ok = rmx_load_data(dnew, fp); /* allocate & load data */
437
438 if (fp != stdin) { /* close input stream */
439 if (inspec[0] == '!')
440 pclose(fp);
441 else
442 fclose(fp);
443 }
444 #ifdef getc_unlocked
445 else
446 funlockfile(fp);
447 #endif
448 if (!ok) { /* load failure? */
449 fprintf(stderr, "Error loading data from: %s\n", inspec);
450 rmx_free(dnew);
451 return(NULL);
452 }
453 /* undo exposure? */
454 if ((dnew->cexp[0] != 1.f) |
455 (dnew->cexp[1] != 1.f) | (dnew->cexp[2] != 1.f)) {
456 double cmlt[MAXCOMP];
457 int i;
458 if (dnew->ncomp > MAXCOMP) {
459 fprintf(stderr, "Excess spectral components in: %s\n",
460 inspec);
461 rmx_free(dnew);
462 return(NULL);
463 }
464 cmlt[0] = 1./dnew->cexp[0];
465 cmlt[1] = 1./dnew->cexp[1];
466 cmlt[2] = 1./dnew->cexp[2];
467 for (i = dnew->ncomp; i-- > 3; )
468 cmlt[i] = cmlt[1]; /* XXX hack! */
469 rmx_scale(dnew, cmlt);
470 setcolor(dnew->cexp, 1.f, 1.f, 1.f);
471 }
472 return(dnew);
473 }
474
475 #if DTrmx_native==DTdouble
476 static int
477 rmx_write_float(const rmx_dtype *dp, int len, FILE *fp)
478 {
479 float val;
480
481 while (len--) {
482 val = (float)*dp++;
483 if (putbinary(&val, sizeof(val), 1, fp) != 1)
484 return(0);
485 }
486 return(1);
487 }
488 #else
489 static int
490 rmx_write_double(const rmx_dtype *dp, int len, FILE *fp)
491 {
492 double val;
493
494 while (len--) {
495 val = *dp++;
496 if (putbinary(&val, sizeof(val), 1, fp) != 1)
497 return(0);
498 }
499 return(1);
500 }
501 #endif
502
503 static int
504 rmx_write_ascii(const rmx_dtype *dp, int nc, int len, FILE *fp)
505 {
506 while (len-- > 0) {
507 int k = nc;
508 while (k-- > 0)
509 fprintf(fp, " %.7e", *dp++);
510 fputc('\t', fp);
511 }
512 return(fputc('\n', fp) != EOF);
513 }
514
515 static int
516 rmx_write_rgbe(const rmx_dtype *dp, int nc, int len, FILE *fp)
517 {
518 COLR *scan;
519 int j;
520
521 if ((nc != 1) & (nc != 3)) return(0);
522 scan = (COLR *)tempbuffer(sizeof(COLR)*len);
523 if (!scan) return(0);
524
525 for (j = 0; j < len; j++, dp += nc)
526 if (nc == 1)
527 setcolr(scan[j], dp[0], dp[0], dp[0]);
528 else
529 setcolr(scan[j], dp[0], dp[1], dp[2]);
530
531 return(fwritecolrs(scan, len, fp) >= 0);
532 }
533
534 static int
535 rmx_write_spec(const rmx_dtype *dp, int nc, int len, FILE *fp)
536 {
537 COLRV *scan;
538 COLORV scol[MAXCOMP];
539 int j, k;
540
541 if ((nc < 3) | (nc > MAXCOMP)) return(0);
542 scan = (COLRV *)tempbuffer((nc+1)*len);
543 if (!scan) return(0);
544 for (j = 0; j < len; j++, dp += nc) {
545 for (k = nc; k--; )
546 scol[k] = dp[k];
547 scolor2scolr(scan+j*(nc+1), scol, nc);
548 }
549 return(fwritescolrs(scan, nc, len, fp) >= 0);
550 }
551
552 /* Check if CIE XYZ primaries were specified */
553 static int
554 findCIEprims(const char *info)
555 {
556 RGBPRIMS prims;
557
558 if (!info)
559 return(0);
560 info = strstr(info, PRIMARYSTR);
561 if (!info || !primsval(prims, info))
562 return(0);
563
564 return((prims[RED][CIEX] > .99) & (prims[RED][CIEY] < .01) &&
565 (prims[GRN][CIEX] < .01) & (prims[GRN][CIEY] > .99) &&
566 (prims[BLU][CIEX] < .01) & (prims[BLU][CIEY] < .01));
567 }
568
569 /* Finish writing header data with resolution and format, returning type used */
570 int
571 rmx_write_header(const RMATRIX *rm, int dtype, FILE *fp)
572 {
573 if (!rm | !fp || rm->ncols <= 0)
574 return(0);
575 if (rm->info)
576 fputs(rm->info, fp);
577 if (dtype == DTfromHeader) {
578 dtype = rm->dtype;
579 #if DTrmx_native==DTfloat
580 if (dtype == DTdouble) /* but stored as float? */
581 dtype = DTfloat;
582 #endif
583 } else if (dtype == DTrgbe && (rm->dtype == DTxyze ||
584 findCIEprims(rm->info)))
585 dtype = DTxyze;
586 else if ((dtype == DTxyze) & (rm->dtype == DTrgbe))
587 dtype = DTrgbe;
588 if ((dtype < DTspec) & (rm->ncomp > 3))
589 dtype = DTspec;
590 else if ((dtype == DTspec) & (rm->ncomp <= 3))
591 return(0);
592
593 if (dtype == DTascii) /* set file type (WINDOWS) */
594 SET_FILE_TEXT(fp);
595 else
596 SET_FILE_BINARY(fp);
597 /* write exposure? */
598 if (rm->ncomp == 3 && (rm->cexp[RED] != rm->cexp[GRN]) |
599 (rm->cexp[GRN] != rm->cexp[BLU]))
600 fputcolcor(rm->cexp, fp);
601 else if (rm->cexp[GRN] != 1.f)
602 fputexpos(rm->cexp[GRN], fp);
603 /* matrix size? */
604 if ((dtype > DTspec) | (rm->nrows <= 0)) {
605 if (rm->nrows > 0)
606 fprintf(fp, "NROWS=%d\n", rm->nrows);
607 fprintf(fp, "NCOLS=%d\n", rm->ncols);
608 }
609 if (dtype >= DTspec) { /* # components & split? */
610 fputncomp(rm->ncomp, fp);
611 if (rm->ncomp > 3 &&
612 memcmp(rm->wlpart, WLPART, sizeof(WLPART)))
613 fputwlsplit(rm->wlpart, fp);
614 } else if ((rm->ncomp != 3) & (rm->ncomp != 1))
615 return(0); /* wrong # components */
616 if ((dtype == DTfloat) | (dtype == DTdouble))
617 fputendian(fp); /* important to record */
618 fputformat(cm_fmt_id[dtype], fp);
619 fputc('\n', fp); /* end of header */
620 if ((dtype <= DTspec) & (rm->nrows > 0))
621 fprtresolu(rm->ncols, rm->nrows, fp);
622 return(dtype);
623 }
624
625 /* Write out matrix data (usually by row) */
626 int
627 rmx_write_data(const rmx_dtype *dp, int nc, int len, int dtype, FILE *fp)
628 {
629 switch (dtype) {
630 #if DTrmx_native==DTdouble
631 case DTfloat:
632 return(rmx_write_float(dp, nc*len, fp));
633 #else
634 case DTdouble:
635 return(rmx_write_double(dp, nc*len, fp));
636 #endif
637 case DTrmx_native:
638 return(putbinary(dp, sizeof(*dp)*nc, len, fp) == len);
639 case DTascii:
640 return(rmx_write_ascii(dp, nc, len, fp));
641 case DTrgbe:
642 case DTxyze:
643 return(rmx_write_rgbe(dp, nc, len, fp));
644 case DTspec:
645 return(rmx_write_spec(dp, nc, len, fp));
646 }
647 return(0);
648 }
649
650 /* Write matrix using file format indicated by dtype */
651 int
652 rmx_write(const RMATRIX *rm, int dtype, FILE *fp)
653 {
654 int ok = 0;
655 int i;
656 /* complete header */
657 dtype = rmx_write_header(rm, dtype, fp);
658 if (dtype <= 0)
659 return(0);
660 #ifdef getc_unlocked
661 flockfile(fp);
662 #endif
663 if (dtype == DTrmx_native) /* write all at once? */
664 ok = rmx_write_data(rm->mtx, rm->ncomp,
665 rm->nrows*rm->ncols, dtype, fp);
666 else /* else row by row */
667 for (i = 0; i < rm->nrows; i++) {
668 ok = rmx_write_data(rmx_val(rm,i,0), rm->ncomp,
669 rm->ncols, dtype, fp);
670 if (!ok) break;
671 }
672
673 if (ok) ok = (fflush(fp) == 0);
674 #ifdef getc_unlocked
675 funlockfile(fp);
676 #endif
677 if (!ok) fputs("Error writing matrix\n", stderr);
678 return(ok);
679 }
680
681 /* Allocate and assign square identity matrix with n components */
682 RMATRIX *
683 rmx_identity(const int dim, const int n)
684 {
685 RMATRIX *rid = rmx_alloc(dim, dim, n);
686 int i, k;
687
688 if (!rid)
689 return(NULL);
690 memset(rid->mtx, 0, rmx_array_size(rid));
691 for (i = dim; i--; ) {
692 rmx_dtype *dp = rmx_lval(rid,i,i);
693 for (k = n; k--; )
694 dp[k] = 1.;
695 }
696 return(rid);
697 }
698
699 /* Duplicate the given matrix (may be unallocated) */
700 RMATRIX *
701 rmx_copy(const RMATRIX *rm)
702 {
703 RMATRIX *dnew;
704
705 if (!rm)
706 return(NULL);
707 dnew = rmx_new(rm->nrows, rm->ncols, rm->ncomp);
708 if (!dnew)
709 return(NULL);
710 if (rm->mtx) {
711 if (!rmx_prepare(dnew)) {
712 rmx_free(dnew);
713 return(NULL);
714 }
715 memcpy(dnew->mtx, rm->mtx, rmx_array_size(dnew));
716 }
717 rmx_addinfo(dnew, rm->info);
718 dnew->dtype = rm->dtype;
719 copycolor(dnew->cexp, rm->cexp);
720 memcpy(dnew->wlpart, rm->wlpart, sizeof(dnew->wlpart));
721 return(dnew);
722 }
723
724 /* Replace data in first matrix with data from second */
725 int
726 rmx_transfer_data(RMATRIX *rdst, RMATRIX *rsrc, int dometa)
727 {
728 if (!rdst | !rsrc)
729 return(0);
730 if (dometa) { /* transfer everything? */
731 rmx_reset(rdst);
732 *rdst = *rsrc;
733 rsrc->info = NULL; rsrc->mapped = NULL; rsrc->mtx = NULL;
734 return(1);
735 }
736 /* just matrix data -- leave metadata */
737 if ((rdst->nrows != rsrc->nrows) |
738 (rdst->ncols != rsrc->ncols) |
739 (rdst->ncomp != rsrc->ncomp))
740 return(0);
741 #ifdef MAP_FILE
742 if (rdst->mapped)
743 munmap(rdst->mapped, rmx_mapped_size(rdst));
744 else
745 #endif
746 if (rdst->pflags & RMF_FREEMEM) {
747 free(rdst->mtx);
748 rdst->pflags &= ~RMF_FREEMEM;
749 }
750 rdst->mapped = rsrc->mapped;
751 rdst->mtx = rsrc->mtx;
752 rdst->pflags |= rsrc->pflags & RMF_FREEMEM;
753 rsrc->mapped = NULL; rsrc->mtx = NULL;
754 return(1);
755 }
756
757 /* Transpose the given matrix */
758 int
759 rmx_transpose(RMATRIX *rm)
760 {
761 RMATRIX dnew;
762 int i, j;
763
764 if (!rm || !rm->mtx | (rm->ncomp > MAXCOMP))
765 return(0);
766 if (rm->info)
767 rmx_addinfo(rm, "Transposed rows and columns\n");
768 if ((rm->nrows == 1) | (rm->ncols == 1)) { /* vector? */
769 j = rm->ncols;
770 rm->ncols = rm->nrows;
771 rm->nrows = j;
772 return(1);
773 }
774 if (rm->nrows == rm->ncols) { /* square matrix case */
775 rmx_dtype val[MAXCOMP];
776 for (j = rm->ncols; j--; )
777 for (i = rm->nrows; i--; ) {
778 if (i == j) continue;
779 memcpy(val, rmx_val(rm,i,j),
780 sizeof(rmx_dtype)*rm->ncomp);
781 memcpy(rmx_lval(rm,i,j), rmx_val(rm,j,i),
782 sizeof(rmx_dtype)*rm->ncomp);
783 memcpy(rmx_val(rm,j,i), val,
784 sizeof(rmx_dtype)*rm->ncomp);
785 }
786 return(1);
787 }
788 memset(&dnew, 0, sizeof(dnew));
789 dnew.ncols = rm->nrows; dnew.nrows = rm->ncols;
790 dnew.ncomp = rm->ncomp;
791 if (!rmx_prepare(&dnew))
792 return(0);
793 rmx_addinfo(&dnew, rm->info);
794 dnew.dtype = rm->dtype;
795 copycolor(dnew.cexp, rm->cexp);
796 memcpy(dnew.wlpart, rm->wlpart, sizeof(dnew.wlpart));
797 for (j = dnew.ncols; j--; )
798 for (i = dnew.nrows; i--; )
799 memcpy(rmx_lval(&dnew,i,j), rmx_val(rm,j,i),
800 sizeof(rmx_dtype)*dnew.ncomp);
801 /* and reassign result */
802 return(rmx_transfer_data(rm, &dnew, 1));
803 }
804
805 /* Multiply (concatenate) two matrices and allocate the result */
806 RMATRIX *
807 rmx_multiply(const RMATRIX *m1, const RMATRIX *m2)
808 {
809 RMATRIX *mres;
810 int i, j, k, h;
811
812 if (!m1 | !m2 || !m1->mtx | !m2->mtx |
813 (m1->ncomp != m2->ncomp) | (m1->ncols != m2->nrows))
814 return(NULL);
815 mres = rmx_alloc(m1->nrows, m2->ncols, m1->ncomp);
816 if (!mres)
817 return(NULL);
818 i = rmx_newtype(m1->dtype, m2->dtype);
819 if (i)
820 mres->dtype = i;
821 else
822 rmx_addinfo(mres, rmx_mismatch_warn);
823 for (i = mres->nrows; i--; )
824 for (j = mres->ncols; j--; )
825 for (k = mres->ncomp; k--; ) {
826 double d = 0;
827 for (h = m1->ncols; h--; )
828 d += (double)rmx_val(m1,i,h)[k] *
829 rmx_val(m2,h,j)[k];
830 rmx_lval(mres,i,j)[k] = (rmx_dtype)d;
831 }
832 return(mres);
833 }
834
835 /* Element-wise multiplication (or division) of m2 into m1 */
836 int
837 rmx_elemult(RMATRIX *m1, const RMATRIX *m2, int divide)
838 {
839 int zeroDivides = 0;
840 int i, j, k;
841
842 if (!m1 | !m2 || !m1->mtx | !m2->mtx |
843 (m1->ncols != m2->ncols) | (m1->nrows != m2->nrows))
844 return(0);
845 if ((m2->ncomp > 1) & (m2->ncomp != m1->ncomp))
846 return(0);
847 i = rmx_newtype(m1->dtype, m2->dtype);
848 if (i)
849 m1->dtype = i;
850 else
851 rmx_addinfo(m1, rmx_mismatch_warn);
852 for (i = m1->nrows; i--; )
853 for (j = m1->ncols; j--; )
854 if (divide) {
855 rmx_dtype d;
856 if (m2->ncomp == 1) {
857 d = rmx_val(m2,i,j)[0];
858 if (d == 0) {
859 ++zeroDivides;
860 for (k = m1->ncomp; k--; )
861 rmx_lval(m1,i,j)[k] = 0;
862 } else {
863 d = 1./d;
864 for (k = m1->ncomp; k--; )
865 rmx_lval(m1,i,j)[k] *= d;
866 }
867 } else
868 for (k = m1->ncomp; k--; ) {
869 d = rmx_val(m2,i,j)[k];
870 if (d == 0) {
871 ++zeroDivides;
872 rmx_lval(m1,i,j)[k] = 0;
873 } else
874 rmx_lval(m1,i,j)[k] /= d;
875 }
876 } else {
877 if (m2->ncomp == 1) {
878 const rmx_dtype d = rmx_val(m2,i,j)[0];
879 for (k = m1->ncomp; k--; )
880 rmx_lval(m1,i,j)[k] *= d;
881 } else
882 for (k = m1->ncomp; k--; )
883 rmx_lval(m1,i,j)[k] *= rmx_val(m2,i,j)[k];
884 }
885 if (zeroDivides) {
886 rmx_addinfo(m1, "WARNING: zero divide(s) corrupted results\n");
887 errno = ERANGE;
888 }
889 return(1);
890 }
891
892 /* Sum second matrix into first, applying scale factor beforehand */
893 int
894 rmx_sum(RMATRIX *msum, const RMATRIX *madd, const double sf[])
895 {
896 double *mysf = NULL;
897 int i, j, k;
898
899 if (!msum | !madd || !msum->mtx | !madd->mtx |
900 (msum->nrows != madd->nrows) |
901 (msum->ncols != madd->ncols) |
902 (msum->ncomp != madd->ncomp))
903 return(0);
904 if (!sf) {
905 mysf = (double *)malloc(sizeof(double)*msum->ncomp);
906 if (!mysf)
907 return(0);
908 for (k = msum->ncomp; k--; )
909 mysf[k] = 1;
910 sf = mysf;
911 }
912 i = rmx_newtype(msum->dtype, madd->dtype);
913 if (i)
914 msum->dtype = i;
915 else
916 rmx_addinfo(msum, rmx_mismatch_warn);
917 for (i = msum->nrows; i--; )
918 for (j = msum->ncols; j--; ) {
919 const rmx_dtype *da = rmx_val(madd,i,j);
920 rmx_dtype *ds = rmx_lval(msum,i,j);
921 for (k = msum->ncomp; k--; )
922 ds[k] += (rmx_dtype)sf[k] * da[k];
923 }
924 if (mysf)
925 free(mysf);
926 return(1);
927 }
928
929 /* Scale the given matrix by the indicated scalar component vector */
930 int
931 rmx_scale(RMATRIX *rm, const double sf[])
932 {
933 int i, j, k;
934
935 if (!rm | !sf || !rm->mtx)
936 return(0);
937 for (i = rm->nrows; i--; )
938 for (j = rm->ncols; j--; ) {
939 rmx_dtype *dp = rmx_lval(rm,i,j);
940 for (k = rm->ncomp; k--; )
941 dp[k] *= (rmx_dtype)sf[k];
942 }
943 if (rm->info)
944 rmx_addinfo(rm, "Applied scalar\n");
945 /* XXX: should record as exposure for COLR and SCOLR types? */
946 return(1);
947 }
948
949 /* Allocate new matrix and apply component transformation */
950 RMATRIX *
951 rmx_transform(const RMATRIX *msrc, int n, const double cmat[])
952 {
953 int i, j, ks, kd;
954 RMATRIX *dnew;
955
956 if (!msrc | (n <= 0) | !cmat || !msrc->mtx)
957 return(NULL);
958 dnew = rmx_alloc(msrc->nrows, msrc->ncols, n);
959 if (!dnew)
960 return(NULL);
961 if (msrc->info) {
962 char buf[128];
963 sprintf(buf, "Applied %dx%d component transform\n",
964 dnew->ncomp, msrc->ncomp);
965 rmx_addinfo(dnew, msrc->info);
966 rmx_addinfo(dnew, buf);
967 }
968 dnew->dtype = msrc->dtype;
969 for (i = dnew->nrows; i--; )
970 for (j = dnew->ncols; j--; ) {
971 const rmx_dtype *ds = rmx_val(msrc,i,j);
972 for (kd = dnew->ncomp; kd--; ) {
973 double d = 0;
974 for (ks = msrc->ncomp; ks--; )
975 d += cmat[kd*msrc->ncomp + ks] * ds[ks];
976 rmx_lval(dnew,i,j)[kd] = (rmx_dtype)d;
977 }
978 }
979 return(dnew);
980 }
981
982 /* Convert a color matrix to newly allocated RMATRIX buffer */
983 RMATRIX *
984 rmx_from_cmatrix(const CMATRIX *cm)
985 {
986 RMATRIX *dnew;
987
988 if (!cm)
989 return(NULL);
990 dnew = rmx_alloc(cm->nrows, cm->ncols, 3);
991 if (!dnew)
992 return(NULL);
993
994 dnew->dtype = sizeof(COLORV)==sizeof(float) ?
995 DTfloat : DTdouble;
996
997 if (sizeof(COLORV) == sizeof(rmx_dtype)) {
998 memcpy(dnew->mtx, cm->cmem, rmx_array_size(dnew));
999 } else {
1000 int i, j;
1001 for (i = dnew->nrows; i--; )
1002 for (j = dnew->ncols; j--; ) {
1003 const COLORV *cv = cm_lval(cm,i,j);
1004 rmx_dtype *dp = rmx_lval(dnew,i,j);
1005 dp[0] = cv[0];
1006 dp[1] = cv[1];
1007 dp[2] = cv[2];
1008 }
1009 }
1010 return(dnew);
1011 }
1012
1013 /* Convert general matrix to newly allocated CMATRIX buffer */
1014 CMATRIX *
1015 cm_from_rmatrix(const RMATRIX *rm)
1016 {
1017 CMATRIX *cnew;
1018
1019 if (!rm || !rm->mtx | (rm->ncomp == 2) | (rm->ncomp > MAXCOMP))
1020 return(NULL);
1021 cnew = cm_alloc(rm->nrows, rm->ncols);
1022 if (!cnew)
1023 return(NULL);
1024 if ((sizeof(COLORV) == sizeof(rmx_dtype)) & (rm->ncomp == 3)) {
1025 memcpy(cnew->cmem, rm->mtx, rmx_array_size(rm));
1026 } else {
1027 int i, j;
1028 for (i = cnew->nrows; i--; )
1029 for (j = cnew->ncols; j--; ) {
1030 const rmx_dtype *dp = rmx_val(rm,i,j);
1031 COLORV *cv = cm_lval(cnew,i,j);
1032 switch (rm->ncomp) {
1033 case 3:
1034 setcolor(cv, dp[0], dp[1], dp[2]);
1035 break;
1036 case 1:
1037 setcolor(cv, dp[0], dp[0], dp[0]);
1038 break;
1039 default: {
1040 COLORV scol[MAXCOMP];
1041 int k;
1042 for (k = rm->ncomp; k--; )
1043 scol[k] = dp[k];
1044 scolor2color(cv, scol, rm->ncomp, rm->wlpart);
1045 } break;
1046 }
1047 }
1048 }
1049 return(cnew);
1050 }