ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmatrix.c
Revision: 2.78
Committed: Thu Feb 29 03:11:13 2024 UTC (4 weeks ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 2.77: +2 -2 lines
Log Message:
fix(rcomb): Spectral scanlines were being output in reverse order

File Contents

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