ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmatrix.c
Revision: 2.93
Committed: Wed Apr 16 22:47:16 2025 UTC (2 weeks, 1 day ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.92: +8 -5 lines
Log Message:
perf(rmtxop): Minor optimization in last change

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: rmatrix.c,v 2.92 2025/04/16 21:54:15 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 #define bmbyte(r,c) bmap[((r)*rm->ncols+(c))>>3]
762 #define bmbit(r,c) (1 << ((r)*rm->ncols+(c) & 7))
763 #define bmop(r,c, op) (bmbyte(r,c) op bmbit(r,c))
764 #define bmtest(r,c) bmop(r,c,&)
765 #define bmset(r,c) bmop(r,c,|=)
766 #define bmtestandset(r,c) (!bmtest(r,c) && bmset(r,c))
767 uby8 *bmap;
768 rmx_dtype val[MAXCOMP];
769 RMATRIX dold;
770 int i, j;
771
772 if (!rm || !rm->mtx | (rm->ncomp > MAXCOMP))
773 return(0);
774 if (rm->info)
775 rmx_addinfo(rm, "Transposed rows and columns\n");
776 if ((rm->nrows == 1) | (rm->ncols == 1)) { /* vector? */
777 j = rm->ncols;
778 rm->ncols = rm->nrows;
779 rm->nrows = j;
780 return(1);
781 }
782 if (rm->nrows == rm->ncols) { /* square matrix case */
783 for (i = rm->nrows; i--; )
784 for (j = rm->ncols; j--; ) {
785 if (i == j) continue;
786 memcpy(val, rmx_val(rm,i,j),
787 sizeof(rmx_dtype)*rm->ncomp);
788 memcpy(rmx_lval(rm,i,j), rmx_val(rm,j,i),
789 sizeof(rmx_dtype)*rm->ncomp);
790 memcpy(rmx_val(rm,j,i), val,
791 sizeof(rmx_dtype)*rm->ncomp);
792 }
793 return(1);
794 }
795 /* clear completion bitmap */
796 bmap = (uby8 *)calloc(((size_t)rm->nrows*rm->ncols+7)>>3, 1);
797 if (!bmap)
798 return(0);
799 dold = *rm;
800 rm->ncols = dold.nrows; rm->nrows = dold.ncols;
801 for (i = rm->nrows; i--; )
802 for (j = rm->ncols; j--; ) {
803 int i0, j0;
804 int i1 = i;
805 size_t j1 = j;
806 if (!bmtestandset(i,j)) continue;
807 memcpy(val, rmx_val(rm,i,j),
808 sizeof(rmx_dtype)*rm->ncomp);
809 for ( ; ; ) { /* value transpose loop */
810 const rmx_dtype *src;
811 i0 = i1; j0 = (int)j1;
812 src = rmx_val(&dold, j0, i0);
813 j1 = (src - dold.mtx)/dold.ncomp;
814 i1 = j1 / rm->ncols;
815 j1 -= i1*rm->ncols;
816 if (!bmtestandset(i1,j1)) break;
817 memcpy(rmx_lval(rm,i0,j0), src,
818 sizeof(rmx_dtype)*rm->ncomp);
819 }
820 /* close the loop */
821 memcpy(rmx_lval(rm,i0,j0), val,
822 sizeof(rmx_dtype)*rm->ncomp);
823 }
824 free(bmap); /* all done! */
825 return(1);
826 #undef bmbyte
827 #undef bmbit
828 #undef bmop
829 #undef bmtest
830 #undef bmset
831 #undef bmtestandset
832 }
833
834 /* Multiply (concatenate) two matrices and allocate the result */
835 RMATRIX *
836 rmx_multiply(const RMATRIX *m1, const RMATRIX *m2)
837 {
838 RMATRIX *mres;
839 int i, j, k, h;
840
841 if (!m1 | !m2 || !m1->mtx | !m2->mtx |
842 (m1->ncomp != m2->ncomp) | (m1->ncols != m2->nrows))
843 return(NULL);
844 mres = rmx_alloc(m1->nrows, m2->ncols, m1->ncomp);
845 if (!mres)
846 return(NULL);
847 i = rmx_newtype(m1->dtype, m2->dtype);
848 if (i)
849 mres->dtype = i;
850 else
851 rmx_addinfo(mres, rmx_mismatch_warn);
852 for (i = mres->nrows; i--; )
853 for (j = mres->ncols; j--; )
854 for (k = mres->ncomp; k--; ) {
855 double d = 0;
856 for (h = m1->ncols; h--; )
857 d += (double)rmx_val(m1,i,h)[k] *
858 rmx_val(m2,h,j)[k];
859 rmx_lval(mres,i,j)[k] = (rmx_dtype)d;
860 }
861 return(mres);
862 }
863
864 /* Element-wise multiplication (or division) of m2 into m1 */
865 int
866 rmx_elemult(RMATRIX *m1, const RMATRIX *m2, int divide)
867 {
868 int zeroDivides = 0;
869 int i, j, k;
870
871 if (!m1 | !m2 || !m1->mtx | !m2->mtx |
872 (m1->ncols != m2->ncols) | (m1->nrows != m2->nrows))
873 return(0);
874 if ((m2->ncomp > 1) & (m2->ncomp != m1->ncomp))
875 return(0);
876 i = rmx_newtype(m1->dtype, m2->dtype);
877 if (i)
878 m1->dtype = i;
879 else
880 rmx_addinfo(m1, rmx_mismatch_warn);
881 for (i = m1->nrows; i--; )
882 for (j = m1->ncols; j--; )
883 if (divide) {
884 rmx_dtype d;
885 if (m2->ncomp == 1) {
886 d = rmx_val(m2,i,j)[0];
887 if (d == 0) {
888 ++zeroDivides;
889 for (k = m1->ncomp; k--; )
890 rmx_lval(m1,i,j)[k] = 0;
891 } else {
892 d = 1./d;
893 for (k = m1->ncomp; k--; )
894 rmx_lval(m1,i,j)[k] *= d;
895 }
896 } else
897 for (k = m1->ncomp; k--; ) {
898 d = rmx_val(m2,i,j)[k];
899 if (d == 0) {
900 ++zeroDivides;
901 rmx_lval(m1,i,j)[k] = 0;
902 } else
903 rmx_lval(m1,i,j)[k] /= d;
904 }
905 } else {
906 if (m2->ncomp == 1) {
907 const rmx_dtype d = rmx_val(m2,i,j)[0];
908 for (k = m1->ncomp; k--; )
909 rmx_lval(m1,i,j)[k] *= d;
910 } else
911 for (k = m1->ncomp; k--; )
912 rmx_lval(m1,i,j)[k] *= rmx_val(m2,i,j)[k];
913 }
914 if (zeroDivides) {
915 rmx_addinfo(m1, "WARNING: zero divide(s) corrupted results\n");
916 errno = ERANGE;
917 }
918 return(1);
919 }
920
921 /* Sum second matrix into first, applying scale factor beforehand */
922 int
923 rmx_sum(RMATRIX *msum, const RMATRIX *madd, const double sf[])
924 {
925 double *mysf = NULL;
926 int i, j, k;
927
928 if (!msum | !madd || !msum->mtx | !madd->mtx |
929 (msum->nrows != madd->nrows) |
930 (msum->ncols != madd->ncols) |
931 (msum->ncomp != madd->ncomp))
932 return(0);
933 if (!sf) {
934 mysf = (double *)malloc(sizeof(double)*msum->ncomp);
935 if (!mysf)
936 return(0);
937 for (k = msum->ncomp; k--; )
938 mysf[k] = 1;
939 sf = mysf;
940 }
941 i = rmx_newtype(msum->dtype, madd->dtype);
942 if (i)
943 msum->dtype = i;
944 else
945 rmx_addinfo(msum, rmx_mismatch_warn);
946 for (i = msum->nrows; i--; )
947 for (j = msum->ncols; j--; ) {
948 const rmx_dtype *da = rmx_val(madd,i,j);
949 rmx_dtype *ds = rmx_lval(msum,i,j);
950 for (k = msum->ncomp; k--; )
951 ds[k] += (rmx_dtype)sf[k] * da[k];
952 }
953 if (mysf)
954 free(mysf);
955 return(1);
956 }
957
958 /* Scale the given matrix by the indicated scalar component vector */
959 int
960 rmx_scale(RMATRIX *rm, const double sf[])
961 {
962 int i, j, k;
963
964 if (!rm | !sf || !rm->mtx)
965 return(0);
966 for (i = rm->nrows; i--; )
967 for (j = rm->ncols; j--; ) {
968 rmx_dtype *dp = rmx_lval(rm,i,j);
969 for (k = rm->ncomp; k--; )
970 dp[k] *= (rmx_dtype)sf[k];
971 }
972 if (rm->info)
973 rmx_addinfo(rm, "Applied scalar\n");
974 /* XXX: should record as exposure for COLR and SCOLR types? */
975 return(1);
976 }
977
978 /* Allocate new matrix and apply component transformation */
979 RMATRIX *
980 rmx_transform(const RMATRIX *msrc, int n, const double cmat[])
981 {
982 int i, j, ks, kd;
983 RMATRIX *dnew;
984
985 if (!msrc | (n <= 0) | !cmat || !msrc->mtx)
986 return(NULL);
987 dnew = rmx_alloc(msrc->nrows, msrc->ncols, n);
988 if (!dnew)
989 return(NULL);
990 if (msrc->info) {
991 char buf[128];
992 sprintf(buf, "Applied %dx%d component transform\n",
993 dnew->ncomp, msrc->ncomp);
994 rmx_addinfo(dnew, msrc->info);
995 rmx_addinfo(dnew, buf);
996 }
997 dnew->dtype = msrc->dtype;
998 for (i = dnew->nrows; i--; )
999 for (j = dnew->ncols; j--; ) {
1000 const rmx_dtype *ds = rmx_val(msrc,i,j);
1001 for (kd = dnew->ncomp; kd--; ) {
1002 double d = 0;
1003 for (ks = msrc->ncomp; ks--; )
1004 d += cmat[kd*msrc->ncomp + ks] * ds[ks];
1005 rmx_lval(dnew,i,j)[kd] = (rmx_dtype)d;
1006 }
1007 }
1008 return(dnew);
1009 }
1010
1011 /* Convert a color matrix to newly allocated RMATRIX buffer */
1012 RMATRIX *
1013 rmx_from_cmatrix(const CMATRIX *cm)
1014 {
1015 RMATRIX *dnew;
1016
1017 if (!cm)
1018 return(NULL);
1019 dnew = rmx_alloc(cm->nrows, cm->ncols, 3);
1020 if (!dnew)
1021 return(NULL);
1022
1023 dnew->dtype = sizeof(COLORV)==sizeof(float) ?
1024 DTfloat : DTdouble;
1025
1026 if (sizeof(COLORV) == sizeof(rmx_dtype)) {
1027 memcpy(dnew->mtx, cm->cmem, rmx_array_size(dnew));
1028 } else {
1029 int i, j;
1030 for (i = dnew->nrows; i--; )
1031 for (j = dnew->ncols; j--; ) {
1032 const COLORV *cv = cm_lval(cm,i,j);
1033 rmx_dtype *dp = rmx_lval(dnew,i,j);
1034 dp[0] = cv[0];
1035 dp[1] = cv[1];
1036 dp[2] = cv[2];
1037 }
1038 }
1039 return(dnew);
1040 }
1041
1042 /* Convert general matrix to newly allocated CMATRIX buffer */
1043 CMATRIX *
1044 cm_from_rmatrix(const RMATRIX *rm)
1045 {
1046 CMATRIX *cnew;
1047
1048 if (!rm || !rm->mtx | (rm->ncomp == 2) | (rm->ncomp > MAXCOMP))
1049 return(NULL);
1050 cnew = cm_alloc(rm->nrows, rm->ncols);
1051 if (!cnew)
1052 return(NULL);
1053 if ((sizeof(COLORV) == sizeof(rmx_dtype)) & (rm->ncomp == 3)) {
1054 memcpy(cnew->cmem, rm->mtx, rmx_array_size(rm));
1055 } else {
1056 int i, j;
1057 for (i = cnew->nrows; i--; )
1058 for (j = cnew->ncols; j--; ) {
1059 const rmx_dtype *dp = rmx_val(rm,i,j);
1060 COLORV *cv = cm_lval(cnew,i,j);
1061 switch (rm->ncomp) {
1062 case 3:
1063 setcolor(cv, dp[0], dp[1], dp[2]);
1064 break;
1065 case 1:
1066 setcolor(cv, dp[0], dp[0], dp[0]);
1067 break;
1068 default:
1069 if (sizeof(COLORV) == sizeof(rmx_dtype)) {
1070 scolor2color(cv, (const COLORV *)dp,
1071 rm->ncomp, rm->wlpart);
1072 } else {
1073 COLORV scol[MAXCOMP];
1074 int k = rm->ncomp;
1075 while (k--) scol[k] = dp[k];
1076 scolor2color(cv, scol, rm->ncomp, rm->wlpart);
1077 }
1078 break;
1079 }
1080 }
1081 }
1082 return(cnew);
1083 }