ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmatrix.c
Revision: 2.52
Committed: Fri Mar 4 17:17:28 2022 UTC (2 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.51: +24 -22 lines
Log Message:
perf(rmtxop): Reduced memory use when writing HDR output from single channel

File Contents

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