ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmatrix.c
Revision: 2.57
Committed: Mon Mar 14 23:28:14 2022 UTC (2 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.56: +2 -2 lines
Log Message:
fix(rmtxop): Fixed recently introduced bug in RGBE/XYZE load function

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: rmatrix.c,v 2.56 2022/03/11 01:11:13 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 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 = (RMATRIX *)calloc(1, sizeof(RMATRIX));
29
30 if (dnew) {
31 dnew->dtype = DTdouble;
32 dnew->nrows = nr;
33 dnew->ncols = nc;
34 dnew->ncomp = n;
35 }
36 return(dnew);
37 }
38
39 /* Prepare a RMATRIX for writing (allocate array if needed) */
40 int
41 rmx_prepare(RMATRIX *rm)
42 {
43 if (!rm) return(0);
44 if (rm->mtx)
45 return(1);
46 rm->mtx = (double *)malloc(array_size(rm));
47 return(rm->mtx != NULL);
48 }
49
50 /* Call rmx_new() and rmx_prepare() */
51 RMATRIX *
52 rmx_alloc(int nr, int nc, int n)
53 {
54 RMATRIX *dnew = rmx_new(nr, nc, n);
55
56 if (dnew && !rmx_prepare(dnew)) {
57 rmx_free(dnew);
58 dnew = NULL;
59 }
60 return(dnew);
61 }
62
63 /* Free a RMATRIX array */
64 void
65 rmx_free(RMATRIX *rm)
66 {
67 if (!rm) return;
68 if (rm->info)
69 free(rm->info);
70 #ifdef MAP_FILE
71 if (rm->mapped)
72 munmap(rm->mapped, mapped_size(rm));
73 else
74 #endif
75 free(rm->mtx);
76 free(rm);
77 }
78
79 /* Resolve data type based on two input types (returns 0 for mismatch) */
80 int
81 rmx_newtype(int dtyp1, int dtyp2)
82 {
83 if ((dtyp1==DTxyze) | (dtyp1==DTrgbe) |
84 (dtyp2==DTxyze) | (dtyp2==DTrgbe)
85 && dtyp1 != dtyp2)
86 return(0);
87 if (dtyp1 < dtyp2)
88 return(dtyp1);
89 return(dtyp2);
90 }
91
92 /* Append header information associated with matrix data */
93 int
94 rmx_addinfo(RMATRIX *rm, const char *info)
95 {
96 int oldlen = 0;
97
98 if (!rm || !info || !*info)
99 return(0);
100 if (!rm->info) {
101 rm->info = (char *)malloc(strlen(info)+1);
102 if (rm->info) rm->info[0] = '\0';
103 } else {
104 oldlen = strlen(rm->info);
105 rm->info = (char *)realloc(rm->info,
106 oldlen+strlen(info)+1);
107 }
108 if (!rm->info)
109 return(0);
110 strcpy(rm->info+oldlen, info);
111 return(1);
112 }
113
114 static int
115 get_dminfo(char *s, void *p)
116 {
117 RMATRIX *ip = (RMATRIX *)p;
118 char fmt[MAXFMTLEN];
119 int i;
120
121 if (headidval(fmt, s))
122 return(0);
123 if (!strncmp(s, "NCOMP=", 6)) {
124 ip->ncomp = atoi(s+6);
125 return(0);
126 }
127 if (!strncmp(s, "NROWS=", 6)) {
128 ip->nrows = atoi(s+6);
129 return(0);
130 }
131 if (!strncmp(s, "NCOLS=", 6)) {
132 ip->ncols = atoi(s+6);
133 return(0);
134 }
135 if ((i = isbigendian(s)) >= 0) {
136 ip->swapin = (nativebigendian() != i);
137 return(0);
138 }
139 if (isexpos(s)) {
140 float f = exposval(s);
141 scalecolor(ip->cexp, f);
142 return(0);
143 }
144 if (iscolcor(s)) {
145 COLOR ctmp;
146 colcorval(ctmp, s);
147 multcolor(ip->cexp, ctmp);
148 return(0);
149 }
150 if (!formatval(fmt, s)) {
151 rmx_addinfo(ip, s);
152 return(0);
153 } /* else check format */
154 for (i = 1; i < DTend; i++)
155 if (!strcmp(fmt, cm_fmt_id[i])) {
156 ip->dtype = i;
157 return(0);
158 }
159 return(-1);
160 }
161
162 static int
163 rmx_load_ascii(RMATRIX *rm, FILE *fp)
164 {
165 int i, j, k;
166
167 if (!rmx_prepare(rm))
168 return(0);
169 for (i = 0; i < rm->nrows; i++)
170 for (j = 0; j < rm->ncols; j++) {
171 double *dp = rmx_lval(rm,i,j);
172 for (k = 0; k < rm->ncomp; k++)
173 if (fscanf(fp, "%lf", &dp[k]) != 1)
174 return(0);
175 }
176 return(1);
177 }
178
179 static int
180 rmx_load_float(RMATRIX *rm, FILE *fp)
181 {
182 int i, j, k;
183 float val[100];
184
185 if (rm->ncomp > 100) {
186 fputs("Unsupported # components in rmx_load_float()\n", stderr);
187 exit(1);
188 }
189 if (!rmx_prepare(rm))
190 return(0);
191 for (i = 0; i < rm->nrows; i++)
192 for (j = 0; j < rm->ncols; j++) {
193 double *dp = rmx_lval(rm,i,j);
194 if (getbinary(val, sizeof(val[0]), rm->ncomp, fp) != rm->ncomp)
195 return(0);
196 if (rm->swapin)
197 swap32((char *)val, rm->ncomp);
198 for (k = rm->ncomp; k--; )
199 dp[k] = val[k];
200 }
201 return(1);
202 }
203
204 static int
205 rmx_load_double(RMATRIX *rm, FILE *fp)
206 {
207 int i;
208 #ifdef MAP_FILE
209 long pos; /* map memory for file > 1MB if possible */
210 if (!rm->swapin && array_size(rm) >= 1L<<20 &&
211 (pos = ftell(fp)) >= 0 && !(pos % sizeof(double))) {
212 rm->mapped = mmap(NULL, array_size(rm)+pos, PROT_READ|PROT_WRITE,
213 MAP_PRIVATE, fileno(fp), 0);
214 if (rm->mapped != MAP_FAILED) {
215 rm->mtx = (double *)rm->mapped + pos/sizeof(double);
216 return(1);
217 } /* else fall back on reading into memory */
218 rm->mapped = NULL;
219 }
220 #endif
221 if (!rmx_prepare(rm))
222 return(0);
223 for (i = 0; i < rm->nrows; i++) {
224 if (getbinary(rmx_lval(rm,i,0), sizeof(double)*rm->ncomp,
225 rm->ncols, fp) != rm->ncols)
226 return(0);
227 if (rm->swapin)
228 swap64((char *)rmx_lval(rm,i,0), rm->ncols*rm->ncomp);
229 }
230 return(1);
231 }
232
233 static int
234 rmx_load_rgbe(RMATRIX *rm, FILE *fp)
235 {
236 COLOR *scan = (COLOR *)malloc(sizeof(COLOR)*rm->ncols);
237 int i, j;
238
239 if (!scan)
240 return(0);
241 if (!rmx_prepare(rm))
242 return(0);
243 for (i = 0; i < rm->nrows; i++) {
244 double *dp = rmx_lval(rm,i,0);
245 if (freadscan(scan, rm->ncols, fp) < 0) {
246 free(scan);
247 return(0);
248 }
249 for (j = 0; j < rm->ncols; j++, dp += 3) {
250 dp[0] = colval(scan[j],RED);
251 dp[1] = colval(scan[j],GRN);
252 dp[2] = colval(scan[j],BLU);
253 }
254 }
255 free(scan);
256 return(1);
257 }
258
259 /* Load matrix from supported file type */
260 RMATRIX *
261 rmx_load(const char *inspec, RMPref rmp)
262 {
263 FILE *fp;
264 RMATRIX *dnew;
265
266 if (!inspec)
267 inspec = stdin_name;
268 else if (!*inspec)
269 return(NULL);
270 if (inspec == stdin_name) { /* reading from stdin? */
271 fp = stdin;
272 } else if (inspec[0] == '!') {
273 if (!(fp = popen(inspec+1, "r")))
274 return(NULL);
275 } else {
276 const char *sp = inspec; /* check suffix */
277 while (*sp)
278 ++sp;
279 while (sp > inspec && sp[-1] != '.')
280 --sp;
281 if (!strcasecmp(sp, "XML")) { /* assume it's a BSDF */
282 CMATRIX *cm = rmp==RMPtrans ? cm_loadBTDF(inspec) :
283 cm_loadBRDF(inspec, rmp==RMPreflB) ;
284 if (!cm)
285 return(NULL);
286 dnew = rmx_from_cmatrix(cm);
287 cm_free(cm);
288 dnew->dtype = DTascii;
289 return(dnew);
290 }
291 /* else open it ourselves */
292 if (!(fp = fopen(inspec, "r")))
293 return(NULL);
294 }
295 SET_FILE_BINARY(fp);
296 #ifdef getc_unlocked
297 flockfile(fp);
298 #endif
299 if (!(dnew = rmx_new(0,0,3))) {
300 fclose(fp);
301 return(NULL);
302 }
303 dnew->dtype = DTascii; /* assumed w/o FORMAT */
304 dnew->cexp[0] = dnew->cexp[1] = dnew->cexp[2] = 1.f;
305 if (getheader(fp, get_dminfo, dnew) < 0) {
306 fclose(fp);
307 return(NULL);
308 }
309 if ((dnew->nrows <= 0) | (dnew->ncols <= 0)) {
310 if (!fscnresolu(&dnew->ncols, &dnew->nrows, fp)) {
311 fclose(fp);
312 return(NULL);
313 }
314 if ((dnew->dtype == DTrgbe) | (dnew->dtype == DTxyze) &&
315 dnew->ncomp != 3) {
316 fclose(fp);
317 return(NULL);
318 }
319 }
320 switch (dnew->dtype) {
321 case DTascii:
322 SET_FILE_TEXT(fp);
323 if (!rmx_load_ascii(dnew, fp))
324 goto loaderr;
325 dnew->dtype = DTascii; /* should leave double? */
326 break;
327 case DTfloat:
328 if (!rmx_load_float(dnew, fp))
329 goto loaderr;
330 dnew->dtype = DTfloat;
331 break;
332 case DTdouble:
333 if (!rmx_load_double(dnew, fp))
334 goto loaderr;
335 dnew->dtype = DTdouble;
336 break;
337 case DTrgbe:
338 case DTxyze:
339 if (!rmx_load_rgbe(dnew, fp))
340 goto loaderr;
341 /* undo exposure? */
342 if ((dnew->cexp[0] != 1.f) | (dnew->cexp[1] != 1.f) |
343 (dnew->cexp[2] != 1.f)) {
344 double cmlt[3];
345 cmlt[0] = 1./dnew->cexp[0];
346 cmlt[1] = 1./dnew->cexp[1];
347 cmlt[2] = 1./dnew->cexp[2];
348 rmx_scale(dnew, cmlt);
349 }
350 dnew->swapin = 0;
351 break;
352 default:
353 goto loaderr;
354 }
355 if (fp != stdin) {
356 if (inspec[0] == '!')
357 pclose(fp);
358 else
359 fclose(fp);
360 }
361 #ifdef getc_unlocked
362 else
363 funlockfile(fp);
364 #endif
365 return(dnew);
366 loaderr: /* should report error? */
367 if (inspec[0] == '!')
368 pclose(fp);
369 else
370 fclose(fp);
371 rmx_free(dnew);
372 return(NULL);
373 }
374
375 static int
376 rmx_write_ascii(const RMATRIX *rm, FILE *fp)
377 {
378 const char *fmt = (rm->dtype == DTfloat) ? " %.7e" :
379 (rm->dtype == DTrgbe) | (rm->dtype == DTxyze) ? " %.3e" :
380 " %.15e" ;
381 int i, j, k;
382
383 for (i = 0; i < rm->nrows; i++) {
384 for (j = 0; j < rm->ncols; j++) {
385 const double *dp = rmx_lval(rm,i,j);
386 for (k = 0; k < rm->ncomp; k++)
387 fprintf(fp, fmt, dp[k]);
388 fputc('\t', fp);
389 }
390 fputc('\n', fp);
391 }
392 return(1);
393 }
394
395 static int
396 rmx_write_float(const RMATRIX *rm, FILE *fp)
397 {
398 int i, j, k;
399 float val[100];
400
401 if (rm->ncomp > 100) {
402 fputs("Unsupported # components in rmx_write_float()\n", stderr);
403 exit(1);
404 }
405 for (i = 0; i < rm->nrows; i++)
406 for (j = 0; j < rm->ncols; j++) {
407 const double *dp = rmx_lval(rm,i,j);
408 for (k = rm->ncomp; k--; )
409 val[k] = (float)dp[k];
410 if (putbinary(val, sizeof(float), rm->ncomp, fp) != rm->ncomp)
411 return(0);
412 }
413 return(1);
414 }
415
416 static int
417 rmx_write_double(const RMATRIX *rm, FILE *fp)
418 {
419 int i;
420
421 for (i = 0; i < rm->nrows; i++)
422 if (putbinary(rmx_lval(rm,i,0), sizeof(double)*rm->ncomp,
423 rm->ncols, fp) != rm->ncols)
424 return(0);
425 return(1);
426 }
427
428 static int
429 rmx_write_rgbe(const RMATRIX *rm, FILE *fp)
430 {
431 COLR *scan = (COLR *)malloc(sizeof(COLR)*rm->ncols);
432 int i, j;
433
434 if (!scan)
435 return(0);
436 for (i = 0; i < rm->nrows; i++) {
437 for (j = rm->ncols; j--; ) {
438 const double *dp = rmx_lval(rm,i,j);
439 if (rm->ncomp == 1)
440 setcolr(scan[j], dp[0], dp[0], dp[0]);
441 else
442 setcolr(scan[j], dp[0], dp[1], dp[2]);
443 }
444 if (fwritecolrs(scan, rm->ncols, fp) < 0) {
445 free(scan);
446 return(0);
447 }
448 }
449 free(scan);
450 return(1);
451 }
452
453 /* Check if CIE XYZ primaries were specified */
454 static int
455 findCIEprims(const char *info)
456 {
457 RGBPRIMS prims;
458
459 if (!info)
460 return(0);
461 info = strstr(info, PRIMARYSTR);
462 if (!info || !primsval(prims, info))
463 return(0);
464
465 return((prims[RED][CIEX] > .99) & (prims[RED][CIEY] < .01) &&
466 (prims[GRN][CIEX] < .01) & (prims[GRN][CIEY] > .99) &&
467 (prims[BLU][CIEX] < .01) & (prims[BLU][CIEY] < .01));
468 }
469
470 /* Write matrix to file type indicated by dtype */
471 int
472 rmx_write(const RMATRIX *rm, int dtype, FILE *fp)
473 {
474 int ok = 1;
475
476 if (!rm | !fp || !rm->mtx)
477 return(0);
478 #ifdef getc_unlocked
479 flockfile(fp);
480 #endif
481 /* complete header */
482 if (rm->info)
483 fputs(rm->info, fp);
484 if (dtype == DTfromHeader)
485 dtype = rm->dtype;
486 else if (dtype == DTrgbe && (rm->dtype == DTxyze ||
487 findCIEprims(rm->info)))
488 dtype = DTxyze;
489 else if ((dtype == DTxyze) & (rm->dtype == DTrgbe))
490 dtype = DTrgbe;
491 if ((dtype != DTrgbe) & (dtype != DTxyze)) {
492 fprintf(fp, "NROWS=%d\n", rm->nrows);
493 fprintf(fp, "NCOLS=%d\n", rm->ncols);
494 fprintf(fp, "NCOMP=%d\n", rm->ncomp);
495 } else if ((rm->ncomp != 3) & (rm->ncomp != 1))
496 return(0); /* wrong # components */
497 if ((dtype == DTfloat) | (dtype == DTdouble))
498 fputendian(fp); /* important to record */
499 fputformat(cm_fmt_id[dtype], fp);
500 fputc('\n', fp);
501 switch (dtype) { /* write data */
502 case DTascii:
503 ok = rmx_write_ascii(rm, fp);
504 break;
505 case DTfloat:
506 ok = rmx_write_float(rm, fp);
507 break;
508 case DTdouble:
509 ok = rmx_write_double(rm, fp);
510 break;
511 case DTrgbe:
512 case DTxyze:
513 fprtresolu(rm->ncols, rm->nrows, fp);
514 ok = rmx_write_rgbe(rm, fp);
515 break;
516 default:
517 return(0);
518 }
519 ok &= (fflush(fp) == 0);
520 #ifdef getc_unlocked
521 funlockfile(fp);
522 #endif
523 return(ok);
524 }
525
526 /* Allocate and assign square identity matrix with n components */
527 RMATRIX *
528 rmx_identity(const int dim, const int n)
529 {
530 RMATRIX *rid = rmx_alloc(dim, dim, n);
531 int i, k;
532
533 if (!rid)
534 return(NULL);
535 memset(rid->mtx, 0, array_size(rid));
536 for (i = dim; i--; ) {
537 double *dp = rmx_lval(rid,i,i);
538 for (k = n; k--; )
539 dp[k] = 1.;
540 }
541 return(rid);
542 }
543
544 /* Duplicate the given matrix */
545 RMATRIX *
546 rmx_copy(const RMATRIX *rm)
547 {
548 RMATRIX *dnew;
549
550 if (!rm)
551 return(NULL);
552 dnew = rmx_alloc(rm->nrows, rm->ncols, rm->ncomp);
553 if (!dnew)
554 return(NULL);
555 rmx_addinfo(dnew, rm->info);
556 dnew->dtype = rm->dtype;
557 memcpy(dnew->mtx, rm->mtx, array_size(dnew));
558 return(dnew);
559 }
560
561 /* Allocate and assign transposed matrix */
562 RMATRIX *
563 rmx_transpose(const RMATRIX *rm)
564 {
565 RMATRIX *dnew;
566 int i, j;
567
568 if (!rm)
569 return(0);
570 if ((rm->nrows == 1) | (rm->ncols == 1)) {
571 dnew = rmx_copy(rm);
572 if (!dnew)
573 return(NULL);
574 dnew->nrows = rm->ncols;
575 dnew->ncols = rm->nrows;
576 return(dnew);
577 }
578 dnew = rmx_alloc(rm->ncols, rm->nrows, rm->ncomp);
579 if (!dnew)
580 return(NULL);
581 if (rm->info) {
582 rmx_addinfo(dnew, rm->info);
583 rmx_addinfo(dnew, "Transposed rows and columns\n");
584 }
585 dnew->dtype = rm->dtype;
586 for (j = dnew->ncols; j--; )
587 for (i = dnew->nrows; i--; )
588 memcpy(rmx_lval(dnew,i,j), rmx_lval(rm,j,i),
589 sizeof(double)*dnew->ncomp);
590 return(dnew);
591 }
592
593 /* Multiply (concatenate) two matrices and allocate the result */
594 RMATRIX *
595 rmx_multiply(const RMATRIX *m1, const RMATRIX *m2)
596 {
597 RMATRIX *mres;
598 int i, j, k, h;
599
600 if (!m1 | !m2 || (m1->ncomp != m2->ncomp) | (m1->ncols != m2->nrows))
601 return(NULL);
602 mres = rmx_alloc(m1->nrows, m2->ncols, m1->ncomp);
603 if (!mres)
604 return(NULL);
605 i = rmx_newtype(m1->dtype, m2->dtype);
606 if (i)
607 mres->dtype = i;
608 else
609 rmx_addinfo(mres, rmx_mismatch_warn);
610 for (i = mres->nrows; i--; )
611 for (j = mres->ncols; j--; )
612 for (k = mres->ncomp; k--; ) {
613 double d = 0;
614 for (h = m1->ncols; h--; )
615 d += rmx_lval(m1,i,h)[k] * rmx_lval(m2,h,j)[k];
616 rmx_lval(mres,i,j)[k] = d;
617 }
618 return(mres);
619 }
620
621 /* Element-wise multiplication (or division) of m2 into m1 */
622 int
623 rmx_elemult(RMATRIX *m1, const RMATRIX *m2, int divide)
624 {
625 int zeroDivides = 0;
626 int i, j, k;
627
628 if (!m1 | !m2 || (m1->ncols != m2->ncols) | (m1->nrows != m2->nrows))
629 return(0);
630 if ((m2->ncomp > 1) & (m2->ncomp != m1->ncomp))
631 return(0);
632 i = rmx_newtype(m1->dtype, m2->dtype);
633 if (i)
634 m1->dtype = i;
635 else
636 rmx_addinfo(m1, rmx_mismatch_warn);
637 for (i = m1->nrows; i--; )
638 for (j = m1->ncols; j--; )
639 if (divide) {
640 double d;
641 if (m2->ncomp == 1) {
642 d = rmx_lval(m2,i,j)[0];
643 if (d == 0) {
644 ++zeroDivides;
645 for (k = m1->ncomp; k--; )
646 rmx_lval(m1,i,j)[k] = 0;
647 } else {
648 d = 1./d;
649 for (k = m1->ncomp; k--; )
650 rmx_lval(m1,i,j)[k] *= d;
651 }
652 } else
653 for (k = m1->ncomp; k--; ) {
654 d = rmx_lval(m2,i,j)[k];
655 if (d == 0) {
656 ++zeroDivides;
657 rmx_lval(m1,i,j)[k] = 0;
658 } else
659 rmx_lval(m1,i,j)[k] /= d;
660 }
661 } else {
662 if (m2->ncomp == 1) {
663 const double d = rmx_lval(m2,i,j)[0];
664 for (k = m1->ncomp; k--; )
665 rmx_lval(m1,i,j)[k] *= d;
666 } else
667 for (k = m1->ncomp; k--; )
668 rmx_lval(m1,i,j)[k] *= rmx_lval(m2,i,j)[k];
669 }
670 if (zeroDivides) {
671 rmx_addinfo(m1, "WARNING: zero divide(s) corrupted results\n");
672 errno = ERANGE;
673 }
674 return(1);
675 }
676
677 /* Sum second matrix into first, applying scale factor beforehand */
678 int
679 rmx_sum(RMATRIX *msum, const RMATRIX *madd, const double sf[])
680 {
681 double *mysf = NULL;
682 int i, j, k;
683
684 if (!msum | !madd ||
685 (msum->nrows != madd->nrows) |
686 (msum->ncols != madd->ncols) |
687 (msum->ncomp != madd->ncomp))
688 return(0);
689 if (!sf) {
690 mysf = (double *)malloc(sizeof(double)*msum->ncomp);
691 if (!mysf)
692 return(0);
693 for (k = msum->ncomp; k--; )
694 mysf[k] = 1;
695 sf = mysf;
696 }
697 i = rmx_newtype(msum->dtype, madd->dtype);
698 if (i)
699 msum->dtype = i;
700 else
701 rmx_addinfo(msum, rmx_mismatch_warn);
702 for (i = msum->nrows; i--; )
703 for (j = msum->ncols; j--; ) {
704 const double *da = rmx_lval(madd,i,j);
705 double *ds = rmx_lval(msum,i,j);
706 for (k = msum->ncomp; k--; )
707 ds[k] += sf[k] * da[k];
708 }
709 if (mysf)
710 free(mysf);
711 return(1);
712 }
713
714 /* Scale the given matrix by the indicated scalar component vector */
715 int
716 rmx_scale(RMATRIX *rm, const double sf[])
717 {
718 int i, j, k;
719
720 if (!rm | !sf)
721 return(0);
722 for (i = rm->nrows; i--; )
723 for (j = rm->ncols; j--; ) {
724 double *dp = rmx_lval(rm,i,j);
725 for (k = rm->ncomp; k--; )
726 dp[k] *= sf[k];
727 }
728 if (rm->info)
729 rmx_addinfo(rm, "Applied scalar\n");
730 return(1);
731 }
732
733 /* Allocate new matrix and apply component transformation */
734 RMATRIX *
735 rmx_transform(const RMATRIX *msrc, int n, const double cmat[])
736 {
737 int i, j, ks, kd;
738 RMATRIX *dnew;
739
740 if (!msrc | (n <= 0) | !cmat)
741 return(NULL);
742 dnew = rmx_alloc(msrc->nrows, msrc->ncols, n);
743 if (!dnew)
744 return(NULL);
745 if (msrc->info) {
746 char buf[128];
747 sprintf(buf, "Applied %dx%d component transform\n",
748 dnew->ncomp, msrc->ncomp);
749 rmx_addinfo(dnew, msrc->info);
750 rmx_addinfo(dnew, buf);
751 }
752 dnew->dtype = msrc->dtype;
753 for (i = dnew->nrows; i--; )
754 for (j = dnew->ncols; j--; ) {
755 const double *ds = rmx_lval(msrc,i,j);
756 for (kd = dnew->ncomp; kd--; ) {
757 double d = 0;
758 for (ks = msrc->ncomp; ks--; )
759 d += cmat[kd*msrc->ncomp + ks] * ds[ks];
760 rmx_lval(dnew,i,j)[kd] = d;
761 }
762 }
763 return(dnew);
764 }
765
766 /* Convert a color matrix to newly allocated RMATRIX buffer */
767 RMATRIX *
768 rmx_from_cmatrix(const CMATRIX *cm)
769 {
770 int i, j;
771 RMATRIX *dnew;
772
773 if (!cm)
774 return(NULL);
775 dnew = rmx_alloc(cm->nrows, cm->ncols, 3);
776 if (!dnew)
777 return(NULL);
778 dnew->dtype = DTfloat;
779 for (i = dnew->nrows; i--; )
780 for (j = dnew->ncols; j--; ) {
781 const COLORV *cv = cm_lval(cm,i,j);
782 double *dp = rmx_lval(dnew,i,j);
783 dp[0] = cv[0];
784 dp[1] = cv[1];
785 dp[2] = cv[2];
786 }
787 return(dnew);
788 }
789
790 /* Convert general matrix to newly allocated CMATRIX buffer */
791 CMATRIX *
792 cm_from_rmatrix(const RMATRIX *rm)
793 {
794 int i, j;
795 CMATRIX *cnew;
796
797 if (!rm || !rm->mtx | ((rm->ncomp != 3) & (rm->ncomp != 1)))
798 return(NULL);
799 cnew = cm_alloc(rm->nrows, rm->ncols);
800 if (!cnew)
801 return(NULL);
802 for (i = cnew->nrows; i--; )
803 for (j = cnew->ncols; j--; ) {
804 const double *dp = rmx_lval(rm,i,j);
805 COLORV *cv = cm_lval(cnew,i,j);
806 if (rm->ncomp == 1)
807 setcolor(cv, dp[0], dp[0], dp[0]);
808 else
809 setcolor(cv, dp[0], dp[1], dp[2]);
810 }
811 return(cnew);
812 }