ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmatrix.c
Revision: 2.51
Committed: Fri Mar 4 02:07:34 2022 UTC (2 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.50: +3 -2 lines
Log Message:
perf(rmtxop): Turned off memory-mapped input for files < 1 MByte

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: rmatrix.c,v 2.50 2022/03/04 01:27:12 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 double d = exposval(s);
141 scalecolor(ip->cexp, d);
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, j;
414
415 for (i = 0; i < rm->nrows; i++)
416 for (j = 0; j < rm->ncols; j++)
417 if (putbinary(&rmx_lval(rm,i,j,0), sizeof(double), rm->ncomp, fp) != rm->ncomp)
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 RMATRIX *mydm = NULL;
466 int ok = 1;
467
468 if (!rm | !fp)
469 return(0);
470 #ifdef getc_unlocked
471 flockfile(fp);
472 #endif
473 /* complete header */
474 if (rm->info)
475 fputs(rm->info, fp);
476 if (dtype == DTfromHeader)
477 dtype = rm->dtype;
478 else if (dtype == DTrgbe && (rm->dtype == DTxyze ||
479 findCIEprims(rm->info)))
480 dtype = DTxyze;
481 else if ((dtype == DTxyze) & (rm->dtype == DTrgbe))
482 dtype = DTrgbe;
483 if ((dtype != DTrgbe) & (dtype != DTxyze)) {
484 fprintf(fp, "NROWS=%d\n", rm->nrows);
485 fprintf(fp, "NCOLS=%d\n", rm->ncols);
486 fprintf(fp, "NCOMP=%d\n", rm->ncomp);
487 } else if (rm->ncomp != 3) { /* wrong # components? */
488 double cmtx[3];
489 if (rm->ncomp != 1) /* only convert grayscale */
490 return(0);
491 cmtx[0] = cmtx[1] = cmtx[2] = 1;
492 mydm = rmx_transform(rm, 3, cmtx);
493 if (!mydm)
494 return(0);
495 rm = mydm;
496 }
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 if (mydm)
524 rmx_free(mydm);
525 return(ok);
526 }
527
528 /* Allocate and assign square identity matrix with n components */
529 RMATRIX *
530 rmx_identity(const int dim, const int n)
531 {
532 RMATRIX *rid = rmx_alloc(dim, dim, n);
533 int i, k;
534
535 if (!rid)
536 return(NULL);
537 memset(rid->mtx, 0, sizeof(rid->mtx[0])*n*dim*dim);
538 for (i = dim; i--; )
539 for (k = n; k--; )
540 rmx_lval(rid,i,i,k) = 1;
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,
558 sizeof(rm->mtx[0])*rm->ncomp*rm->nrows*rm->ncols);
559 return(dnew);
560 }
561
562 /* Allocate and assign transposed matrix */
563 RMATRIX *
564 rmx_transpose(const RMATRIX *rm)
565 {
566 RMATRIX *dnew;
567 int i, j, k;
568
569 if (!rm)
570 return(0);
571 if ((rm->nrows == 1) | (rm->ncols == 1)) {
572 dnew = rmx_copy(rm);
573 if (!dnew)
574 return(NULL);
575 dnew->nrows = rm->ncols;
576 dnew->ncols = rm->nrows;
577 return(dnew);
578 }
579 dnew = rmx_alloc(rm->ncols, rm->nrows, rm->ncomp);
580 if (!dnew)
581 return(NULL);
582 if (rm->info) {
583 rmx_addinfo(dnew, rm->info);
584 rmx_addinfo(dnew, "Transposed rows and columns\n");
585 }
586 dnew->dtype = rm->dtype;
587 for (i = dnew->nrows; i--; )
588 for (j = dnew->ncols; j--; )
589 for (k = dnew->ncomp; k--; )
590 rmx_lval(dnew,i,j,k) = rmx_lval(rm,j,i,k);
591 return(dnew);
592 }
593
594 /* Multiply (concatenate) two matrices and allocate the result */
595 RMATRIX *
596 rmx_multiply(const RMATRIX *m1, const RMATRIX *m2)
597 {
598 RMATRIX *mres;
599 int i, j, k, h;
600
601 if (!m1 | !m2 || (m1->ncomp != m2->ncomp) | (m1->ncols != m2->nrows))
602 return(NULL);
603 mres = rmx_alloc(m1->nrows, m2->ncols, m1->ncomp);
604 if (!mres)
605 return(NULL);
606 i = rmx_newtype(m1->dtype, m2->dtype);
607 if (i)
608 mres->dtype = i;
609 else
610 rmx_addinfo(mres, rmx_mismatch_warn);
611 for (i = mres->nrows; i--; )
612 for (j = mres->ncols; j--; )
613 for (k = mres->ncomp; k--; ) {
614 long double d = 0;
615 for (h = m1->ncols; h--; )
616 d += rmx_lval(m1,i,h,k) * rmx_lval(m2,h,j,k);
617 rmx_lval(mres,i,j,k) = (double)d;
618 }
619 return(mres);
620 }
621
622 /* Element-wise multiplication (or division) of m2 into m1 */
623 int
624 rmx_elemult(RMATRIX *m1, const RMATRIX *m2, int divide)
625 {
626 int zeroDivides = 0;
627 int i, j, k;
628
629 if (!m1 | !m2 || (m1->ncols != m2->ncols) | (m1->nrows != m2->nrows))
630 return(0);
631 if ((m2->ncomp > 1) & (m2->ncomp != m1->ncomp))
632 return(0);
633 i = rmx_newtype(m1->dtype, m2->dtype);
634 if (i)
635 m1->dtype = i;
636 else
637 rmx_addinfo(m1, rmx_mismatch_warn);
638 for (i = m1->nrows; i--; )
639 for (j = m1->ncols; j--; )
640 if (divide) {
641 double d;
642 if (m2->ncomp == 1) {
643 d = rmx_lval(m2,i,j,0);
644 if (d == 0) {
645 ++zeroDivides;
646 for (k = m1->ncomp; k--; )
647 rmx_lval(m1,i,j,k) = 0;
648 } else {
649 d = 1./d;
650 for (k = m1->ncomp; k--; )
651 rmx_lval(m1,i,j,k) *= d;
652 }
653 } else
654 for (k = m1->ncomp; k--; ) {
655 d = rmx_lval(m2,i,j,k);
656 if (d == 0) {
657 ++zeroDivides;
658 rmx_lval(m1,i,j,k) = 0;
659 } else
660 rmx_lval(m1,i,j,k) /= d;
661 }
662 } else {
663 if (m2->ncomp == 1) {
664 const double d = rmx_lval(m2,i,j,0);
665 for (k = m1->ncomp; k--; )
666 rmx_lval(m1,i,j,k) *= d;
667 } else
668 for (k = m1->ncomp; k--; )
669 rmx_lval(m1,i,j,k) *= rmx_lval(m2,i,j,k);
670 }
671 if (zeroDivides) {
672 rmx_addinfo(m1, "WARNING: zero divide(s) corrupted results\n");
673 errno = ERANGE;
674 }
675 return(1);
676 }
677
678 /* Sum second matrix into first, applying scale factor beforehand */
679 int
680 rmx_sum(RMATRIX *msum, const RMATRIX *madd, const double sf[])
681 {
682 double *mysf = NULL;
683 int i, j, k;
684
685 if (!msum | !madd ||
686 (msum->nrows != madd->nrows) |
687 (msum->ncols != madd->ncols) |
688 (msum->ncomp != madd->ncomp))
689 return(0);
690 if (!sf) {
691 mysf = (double *)malloc(sizeof(double)*msum->ncomp);
692 if (!mysf)
693 return(0);
694 for (k = msum->ncomp; k--; )
695 mysf[k] = 1;
696 sf = mysf;
697 }
698 i = rmx_newtype(msum->dtype, madd->dtype);
699 if (i)
700 msum->dtype = i;
701 else
702 rmx_addinfo(msum, rmx_mismatch_warn);
703 for (i = msum->nrows; i--; )
704 for (j = msum->ncols; j--; )
705 for (k = msum->ncomp; k--; )
706 rmx_lval(msum,i,j,k) += sf[k] * rmx_lval(madd,i,j,k);
707 if (mysf)
708 free(mysf);
709 return(1);
710 }
711
712 /* Scale the given matrix by the indicated scalar component vector */
713 int
714 rmx_scale(RMATRIX *rm, const double sf[])
715 {
716 int i, j, k;
717
718 if (!rm | !sf)
719 return(0);
720 for (i = rm->nrows; i--; )
721 for (j = rm->ncols; j--; )
722 for (k = rm->ncomp; k--; )
723 rmx_lval(rm,i,j,k) *= sf[k];
724
725 if (rm->info)
726 rmx_addinfo(rm, "Applied scalar\n");
727 return(1);
728 }
729
730 /* Allocate new matrix and apply component transformation */
731 RMATRIX *
732 rmx_transform(const RMATRIX *msrc, int n, const double cmat[])
733 {
734 int i, j, ks, kd;
735 RMATRIX *dnew;
736
737 if (!msrc | (n <= 0) | !cmat)
738 return(NULL);
739 dnew = rmx_alloc(msrc->nrows, msrc->ncols, n);
740 if (!dnew)
741 return(NULL);
742 if (msrc->info) {
743 char buf[128];
744 sprintf(buf, "Applied %dx%d component transform\n",
745 dnew->ncomp, msrc->ncomp);
746 rmx_addinfo(dnew, msrc->info);
747 rmx_addinfo(dnew, buf);
748 }
749 dnew->dtype = msrc->dtype;
750 for (i = dnew->nrows; i--; )
751 for (j = dnew->ncols; j--; )
752 for (kd = dnew->ncomp; kd--; ) {
753 double d = 0;
754 for (ks = msrc->ncomp; ks--; )
755 d += cmat[kd*msrc->ncomp + ks] * rmx_lval(msrc,i,j,ks);
756 rmx_lval(dnew,i,j,kd) = d;
757 }
758 return(dnew);
759 }
760
761 /* Convert a color matrix to newly allocated RMATRIX buffer */
762 RMATRIX *
763 rmx_from_cmatrix(const CMATRIX *cm)
764 {
765 int i, j;
766 RMATRIX *dnew;
767
768 if (!cm)
769 return(NULL);
770 dnew = rmx_alloc(cm->nrows, cm->ncols, 3);
771 if (!dnew)
772 return(NULL);
773 dnew->dtype = DTfloat;
774 for (i = dnew->nrows; i--; )
775 for (j = dnew->ncols; j--; ) {
776 const COLORV *cv = cm_lval(cm,i,j);
777 rmx_lval(dnew,i,j,0) = cv[0];
778 rmx_lval(dnew,i,j,1) = cv[1];
779 rmx_lval(dnew,i,j,2) = cv[2];
780 }
781 return(dnew);
782 }
783
784 /* Convert general matrix to newly allocated CMATRIX buffer */
785 CMATRIX *
786 cm_from_rmatrix(const RMATRIX *rm)
787 {
788 int i, j;
789 CMATRIX *cnew;
790
791 if (!rm || rm->ncomp != 3)
792 return(NULL);
793 cnew = cm_alloc(rm->nrows, rm->ncols);
794 if (!cnew)
795 return(NULL);
796 for (i = cnew->nrows; i--; )
797 for (j = cnew->ncols; j--; ) {
798 COLORV *cv = cm_lval(cnew,i,j);
799 cv[0] = (COLORV)rmx_lval(rm,i,j,0);
800 cv[1] = (COLORV)rmx_lval(rm,i,j,1);
801 cv[2] = (COLORV)rmx_lval(rm,i,j,2);
802 }
803 return(cnew);
804 }