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, 1 month 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

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.57 static const char RCSid[] = "$Id: rmatrix.c,v 2.56 2022/03/11 01:11:13 greg Exp $";
3 greg 2.1 #endif
4     /*
5     * General matrix operations.
6     */
7    
8     #include <stdlib.h>
9 greg 2.25 #include <errno.h>
10 greg 2.21 #include "rtio.h"
11 schorsch 2.20 #include "platform.h"
12 greg 2.1 #include "resolu.h"
13 greg 2.19 #include "paths.h"
14 greg 2.1 #include "rmatrix.h"
15 greg 2.50 #if !defined(_WIN32) && !defined(_WIN64)
16     #include <sys/mman.h>
17     #endif
18 greg 2.1
19 greg 2.6 static char rmx_mismatch_warn[] = "WARNING: data type mismatch\n";
20 greg 2.1
21 greg 2.50 #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 greg 2.1 RMATRIX *
26 greg 2.50 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 greg 2.1 rmx_alloc(int nr, int nc, int n)
53     {
54 greg 2.50 RMATRIX *dnew = rmx_new(nr, nc, n);
55 greg 2.1
56 greg 2.50 if (dnew && !rmx_prepare(dnew)) {
57     rmx_free(dnew);
58     dnew = NULL;
59     }
60 greg 2.1 return(dnew);
61     }
62    
63 greg 2.6 /* 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 greg 2.50 #ifdef MAP_FILE
71     if (rm->mapped)
72     munmap(rm->mapped, mapped_size(rm));
73     else
74     #endif
75     free(rm->mtx);
76 greg 2.6 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 greg 2.14 if ((dtyp1==DTxyze) | (dtyp1==DTrgbe) |
84     (dtyp2==DTxyze) | (dtyp2==DTrgbe)
85     && dtyp1 != dtyp2)
86 greg 2.6 return(0);
87     if (dtyp1 < dtyp2)
88     return(dtyp1);
89     return(dtyp2);
90     }
91    
92 greg 2.5 /* Append header information associated with matrix data */
93     int
94     rmx_addinfo(RMATRIX *rm, const char *info)
95     {
96 greg 2.50 int oldlen = 0;
97    
98 greg 2.33 if (!rm || !info || !*info)
99 greg 2.5 return(0);
100 greg 2.8 if (!rm->info) {
101 greg 2.5 rm->info = (char *)malloc(strlen(info)+1);
102 greg 2.8 if (rm->info) rm->info[0] = '\0';
103 greg 2.50 } else {
104     oldlen = strlen(rm->info);
105 greg 2.5 rm->info = (char *)realloc(rm->info,
106 greg 2.50 oldlen+strlen(info)+1);
107     }
108 greg 2.5 if (!rm->info)
109     return(0);
110 greg 2.50 strcpy(rm->info+oldlen, info);
111 greg 2.5 return(1);
112     }
113    
114 greg 2.1 static int
115     get_dminfo(char *s, void *p)
116     {
117 greg 2.6 RMATRIX *ip = (RMATRIX *)p;
118 greg 2.29 char fmt[MAXFMTLEN];
119 greg 2.1 int i;
120    
121 greg 2.6 if (headidval(fmt, s))
122     return(0);
123 greg 2.1 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 greg 2.35 if ((i = isbigendian(s)) >= 0) {
136     ip->swapin = (nativebigendian() != i);
137     return(0);
138     }
139 greg 2.40 if (isexpos(s)) {
140 greg 2.52 float f = exposval(s);
141     scalecolor(ip->cexp, f);
142 greg 2.40 return(0);
143     }
144     if (iscolcor(s)) {
145     COLOR ctmp;
146     colcorval(ctmp, s);
147 greg 2.50 multcolor(ip->cexp, ctmp);
148 greg 2.40 return(0);
149     }
150 greg 2.5 if (!formatval(fmt, s)) {
151 greg 2.6 rmx_addinfo(ip, s);
152 greg 2.1 return(0);
153 greg 2.50 } /* else check format */
154 greg 2.1 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 greg 2.17
167 greg 2.50 if (!rmx_prepare(rm))
168     return(0);
169 greg 2.1 for (i = 0; i < rm->nrows; i++)
170 greg 2.53 for (j = 0; j < rm->ncols; j++) {
171     double *dp = rmx_lval(rm,i,j);
172 greg 2.1 for (k = 0; k < rm->ncomp; k++)
173 greg 2.53 if (fscanf(fp, "%lf", &dp[k]) != 1)
174 greg 2.1 return(0);
175 greg 2.53 }
176 greg 2.1 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 greg 2.50 if (!rmx_prepare(rm))
190     return(0);
191 greg 2.1 for (i = 0; i < rm->nrows; i++)
192     for (j = 0; j < rm->ncols; j++) {
193 greg 2.53 double *dp = rmx_lval(rm,i,j);
194 greg 2.21 if (getbinary(val, sizeof(val[0]), rm->ncomp, fp) != rm->ncomp)
195 greg 2.1 return(0);
196 greg 2.35 if (rm->swapin)
197     swap32((char *)val, rm->ncomp);
198 greg 2.1 for (k = rm->ncomp; k--; )
199 greg 2.53 dp[k] = val[k];
200 greg 2.1 }
201     return(1);
202     }
203    
204     static int
205     rmx_load_double(RMATRIX *rm, FILE *fp)
206     {
207 greg 2.41 int i;
208 greg 2.50 #ifdef MAP_FILE
209 greg 2.56 long pos; /* map memory for file > 1MB if possible */
210 greg 2.51 if (!rm->swapin && array_size(rm) >= 1L<<20 &&
211     (pos = ftell(fp)) >= 0 && !(pos % sizeof(double))) {
212 greg 2.50 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 greg 2.56 } /* else fall back on reading into memory */
218 greg 2.50 rm->mapped = NULL;
219 greg 2.41 }
220 greg 2.50 #endif
221     if (!rmx_prepare(rm))
222     return(0);
223 greg 2.41 for (i = 0; i < rm->nrows; i++) {
224 greg 2.53 if (getbinary(rmx_lval(rm,i,0), sizeof(double)*rm->ncomp,
225 greg 2.41 rm->ncols, fp) != rm->ncols)
226     return(0);
227 greg 2.35 if (rm->swapin)
228 greg 2.53 swap64((char *)rmx_lval(rm,i,0), rm->ncols*rm->ncomp);
229 greg 2.41 }
230 greg 2.1 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 greg 2.33 if (!scan)
240 greg 2.1 return(0);
241 greg 2.50 if (!rmx_prepare(rm))
242     return(0);
243 greg 2.1 for (i = 0; i < rm->nrows; i++) {
244 greg 2.57 double *dp = rmx_lval(rm,i,0);
245 greg 2.1 if (freadscan(scan, rm->ncols, fp) < 0) {
246     free(scan);
247     return(0);
248     }
249 greg 2.53 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 greg 2.1 }
254     }
255     free(scan);
256     return(1);
257     }
258    
259     /* Load matrix from supported file type */
260     RMATRIX *
261 greg 2.46 rmx_load(const char *inspec, RMPref rmp)
262 greg 2.1 {
263 greg 2.46 FILE *fp;
264 greg 2.1 RMATRIX *dnew;
265    
266 greg 2.47 if (!inspec)
267     inspec = stdin_name;
268     else if (!*inspec)
269 greg 2.46 return(NULL);
270     if (inspec == stdin_name) { /* reading from stdin? */
271     fp = stdin;
272 greg 2.13 } else if (inspec[0] == '!') {
273 greg 2.33 if (!(fp = popen(inspec+1, "r")))
274 greg 2.13 return(NULL);
275 greg 2.1 } else {
276 greg 2.13 const char *sp = inspec; /* check suffix */
277 greg 2.1 while (*sp)
278     ++sp;
279 greg 2.13 while (sp > inspec && sp[-1] != '.')
280 greg 2.1 --sp;
281     if (!strcasecmp(sp, "XML")) { /* assume it's a BSDF */
282 greg 2.46 CMATRIX *cm = rmp==RMPtrans ? cm_loadBTDF(inspec) :
283     cm_loadBRDF(inspec, rmp==RMPreflB) ;
284 greg 2.33 if (!cm)
285 greg 2.1 return(NULL);
286     dnew = rmx_from_cmatrix(cm);
287     cm_free(cm);
288 greg 2.18 dnew->dtype = DTascii;
289 greg 2.1 return(dnew);
290     }
291     /* else open it ourselves */
292 greg 2.46 if (!(fp = fopen(inspec, "r")))
293 greg 2.1 return(NULL);
294     }
295 greg 2.46 SET_FILE_BINARY(fp);
296 greg 2.1 #ifdef getc_unlocked
297     flockfile(fp);
298     #endif
299 greg 2.50 if (!(dnew = rmx_new(0,0,3))) {
300 greg 2.1 fclose(fp);
301     return(NULL);
302     }
303 greg 2.50 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 greg 2.1 fclose(fp);
312     return(NULL);
313     }
314 greg 2.50 if ((dnew->dtype == DTrgbe) | (dnew->dtype == DTxyze) &&
315     dnew->ncomp != 3) {
316 greg 2.4 fclose(fp);
317     return(NULL);
318     }
319 greg 2.1 }
320 greg 2.50 switch (dnew->dtype) {
321 greg 2.1 case DTascii:
322 greg 2.37 SET_FILE_TEXT(fp);
323 greg 2.1 if (!rmx_load_ascii(dnew, fp))
324     goto loaderr;
325 greg 2.11 dnew->dtype = DTascii; /* should leave double? */
326 greg 2.1 break;
327     case DTfloat:
328     if (!rmx_load_float(dnew, fp))
329     goto loaderr;
330 greg 2.6 dnew->dtype = DTfloat;
331 greg 2.1 break;
332     case DTdouble:
333     if (!rmx_load_double(dnew, fp))
334     goto loaderr;
335 greg 2.6 dnew->dtype = DTdouble;
336 greg 2.1 break;
337     case DTrgbe:
338     case DTxyze:
339     if (!rmx_load_rgbe(dnew, fp))
340     goto loaderr;
341 greg 2.50 /* 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 greg 2.40 }
350 greg 2.50 dnew->swapin = 0;
351 greg 2.1 break;
352     default:
353     goto loaderr;
354     }
355 greg 2.13 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 greg 2.1 return(dnew);
366     loaderr: /* should report error? */
367 greg 2.13 if (inspec[0] == '!')
368     pclose(fp);
369     else
370     fclose(fp);
371 greg 2.1 rmx_free(dnew);
372     return(NULL);
373     }
374    
375     static int
376     rmx_write_ascii(const RMATRIX *rm, FILE *fp)
377     {
378 greg 2.30 const char *fmt = (rm->dtype == DTfloat) ? " %.7e" :
379     (rm->dtype == DTrgbe) | (rm->dtype == DTxyze) ? " %.3e" :
380     " %.15e" ;
381 greg 2.1 int i, j, k;
382 greg 2.17
383 greg 2.1 for (i = 0; i < rm->nrows; i++) {
384     for (j = 0; j < rm->ncols; j++) {
385 greg 2.53 const double *dp = rmx_lval(rm,i,j);
386 greg 2.1 for (k = 0; k < rm->ncomp; k++)
387 greg 2.53 fprintf(fp, fmt, dp[k]);
388 greg 2.1 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 greg 2.53 const double *dp = rmx_lval(rm,i,j);
408 greg 2.1 for (k = rm->ncomp; k--; )
409 greg 2.53 val[k] = (float)dp[k];
410     if (putbinary(val, sizeof(float), rm->ncomp, fp) != rm->ncomp)
411 greg 2.1 return(0);
412     }
413     return(1);
414     }
415    
416     static int
417     rmx_write_double(const RMATRIX *rm, FILE *fp)
418     {
419 greg 2.52 int i;
420 greg 2.1
421     for (i = 0; i < rm->nrows; i++)
422 greg 2.53 if (putbinary(rmx_lval(rm,i,0), sizeof(double)*rm->ncomp,
423 greg 2.52 rm->ncols, fp) != rm->ncols)
424 greg 2.1 return(0);
425     return(1);
426     }
427    
428     static int
429     rmx_write_rgbe(const RMATRIX *rm, FILE *fp)
430     {
431 greg 2.27 COLR *scan = (COLR *)malloc(sizeof(COLR)*rm->ncols);
432 greg 2.1 int i, j;
433    
434 greg 2.33 if (!scan)
435 greg 2.1 return(0);
436     for (i = 0; i < rm->nrows; i++) {
437 greg 2.53 for (j = rm->ncols; j--; ) {
438     const double *dp = rmx_lval(rm,i,j);
439 greg 2.54 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 greg 2.53 }
444 greg 2.27 if (fwritecolrs(scan, rm->ncols, fp) < 0) {
445 greg 2.1 free(scan);
446     return(0);
447     }
448     }
449     free(scan);
450     return(1);
451     }
452    
453 greg 2.48 /* 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 greg 2.6 /* Write matrix to file type indicated by dtype */
471 greg 2.10 int
472 greg 2.1 rmx_write(const RMATRIX *rm, int dtype, FILE *fp)
473     {
474     int ok = 1;
475    
476 greg 2.52 if (!rm | !fp || !rm->mtx)
477 greg 2.1 return(0);
478 greg 2.17 #ifdef getc_unlocked
479     flockfile(fp);
480     #endif
481 greg 2.1 /* complete header */
482 greg 2.5 if (rm->info)
483     fputs(rm->info, fp);
484 greg 2.6 if (dtype == DTfromHeader)
485     dtype = rm->dtype;
486 greg 2.48 else if (dtype == DTrgbe && (rm->dtype == DTxyze ||
487     findCIEprims(rm->info)))
488 greg 2.6 dtype = DTxyze;
489 greg 2.9 else if ((dtype == DTxyze) & (rm->dtype == DTrgbe))
490 greg 2.6 dtype = DTrgbe;
491 greg 2.1 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 greg 2.54 } else if ((rm->ncomp != 3) & (rm->ncomp != 1))
496     return(0); /* wrong # components */
497 greg 2.35 if ((dtype == DTfloat) | (dtype == DTdouble))
498     fputendian(fp); /* important to record */
499 greg 2.49 fputformat(cm_fmt_id[dtype], fp);
500 greg 2.1 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 greg 2.17 #ifdef getc_unlocked
521     funlockfile(fp);
522     #endif
523 greg 2.10 return(ok);
524 greg 2.1 }
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 greg 2.12 int i, k;
532 greg 2.1
533 greg 2.33 if (!rid)
534 greg 2.1 return(NULL);
535 greg 2.52 memset(rid->mtx, 0, array_size(rid));
536 greg 2.53 for (i = dim; i--; ) {
537     double *dp = rmx_lval(rid,i,i);
538 greg 2.12 for (k = n; k--; )
539 greg 2.53 dp[k] = 1.;
540     }
541 greg 2.1 return(rid);
542     }
543    
544     /* Duplicate the given matrix */
545     RMATRIX *
546     rmx_copy(const RMATRIX *rm)
547     {
548     RMATRIX *dnew;
549    
550 greg 2.33 if (!rm)
551 greg 2.1 return(NULL);
552     dnew = rmx_alloc(rm->nrows, rm->ncols, rm->ncomp);
553 greg 2.33 if (!dnew)
554 greg 2.1 return(NULL);
555 greg 2.5 rmx_addinfo(dnew, rm->info);
556 greg 2.6 dnew->dtype = rm->dtype;
557 greg 2.52 memcpy(dnew->mtx, rm->mtx, array_size(dnew));
558 greg 2.1 return(dnew);
559     }
560    
561 greg 2.2 /* Allocate and assign transposed matrix */
562     RMATRIX *
563     rmx_transpose(const RMATRIX *rm)
564 greg 2.1 {
565 greg 2.2 RMATRIX *dnew;
566 greg 2.55 int i, j;
567 greg 2.1
568 greg 2.33 if (!rm)
569 greg 2.1 return(0);
570 greg 2.31 if ((rm->nrows == 1) | (rm->ncols == 1)) {
571     dnew = rmx_copy(rm);
572 greg 2.33 if (!dnew)
573 greg 2.32 return(NULL);
574 greg 2.31 dnew->nrows = rm->ncols;
575     dnew->ncols = rm->nrows;
576     return(dnew);
577     }
578 greg 2.2 dnew = rmx_alloc(rm->ncols, rm->nrows, rm->ncomp);
579 greg 2.33 if (!dnew)
580 greg 2.2 return(NULL);
581 greg 2.5 if (rm->info) {
582     rmx_addinfo(dnew, rm->info);
583     rmx_addinfo(dnew, "Transposed rows and columns\n");
584     }
585 greg 2.6 dnew->dtype = rm->dtype;
586 greg 2.56 for (j = dnew->ncols; j--; )
587     for (i = dnew->nrows; i--; )
588 greg 2.55 memcpy(rmx_lval(dnew,i,j), rmx_lval(rm,j,i),
589 greg 2.53 sizeof(double)*dnew->ncomp);
590 greg 2.2 return(dnew);
591 greg 2.1 }
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 greg 2.33 if (!m1 | !m2 || (m1->ncomp != m2->ncomp) | (m1->ncols != m2->nrows))
601 greg 2.1 return(NULL);
602     mres = rmx_alloc(m1->nrows, m2->ncols, m1->ncomp);
603 greg 2.33 if (!mres)
604 greg 2.1 return(NULL);
605 greg 2.6 i = rmx_newtype(m1->dtype, m2->dtype);
606     if (i)
607     mres->dtype = i;
608     else
609     rmx_addinfo(mres, rmx_mismatch_warn);
610 greg 2.1 for (i = mres->nrows; i--; )
611     for (j = mres->ncols; j--; )
612 greg 2.8 for (k = mres->ncomp; k--; ) {
613 greg 2.53 double d = 0;
614 greg 2.8 for (h = m1->ncols; h--; )
615 greg 2.53 d += rmx_lval(m1,i,h)[k] * rmx_lval(m2,h,j)[k];
616     rmx_lval(mres,i,j)[k] = d;
617 greg 2.1 }
618     return(mres);
619     }
620    
621 greg 2.25 /* 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 greg 2.33 if (!m1 | !m2 || (m1->ncols != m2->ncols) | (m1->nrows != m2->nrows))
629 greg 2.25 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 greg 2.53 d = rmx_lval(m2,i,j)[0];
643 greg 2.25 if (d == 0) {
644     ++zeroDivides;
645     for (k = m1->ncomp; k--; )
646 greg 2.53 rmx_lval(m1,i,j)[k] = 0;
647 greg 2.25 } else {
648     d = 1./d;
649     for (k = m1->ncomp; k--; )
650 greg 2.53 rmx_lval(m1,i,j)[k] *= d;
651 greg 2.25 }
652     } else
653     for (k = m1->ncomp; k--; ) {
654 greg 2.53 d = rmx_lval(m2,i,j)[k];
655 greg 2.25 if (d == 0) {
656     ++zeroDivides;
657 greg 2.53 rmx_lval(m1,i,j)[k] = 0;
658 greg 2.25 } else
659 greg 2.53 rmx_lval(m1,i,j)[k] /= d;
660 greg 2.25 }
661     } else {
662     if (m2->ncomp == 1) {
663 greg 2.53 const double d = rmx_lval(m2,i,j)[0];
664 greg 2.25 for (k = m1->ncomp; k--; )
665 greg 2.53 rmx_lval(m1,i,j)[k] *= d;
666 greg 2.25 } else
667     for (k = m1->ncomp; k--; )
668 greg 2.53 rmx_lval(m1,i,j)[k] *= rmx_lval(m2,i,j)[k];
669 greg 2.25 }
670     if (zeroDivides) {
671 greg 2.34 rmx_addinfo(m1, "WARNING: zero divide(s) corrupted results\n");
672 greg 2.25 errno = ERANGE;
673     }
674     return(1);
675     }
676    
677 greg 2.1 /* 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 greg 2.33 if (!msum | !madd ||
685 greg 2.1 (msum->nrows != madd->nrows) |
686     (msum->ncols != madd->ncols) |
687     (msum->ncomp != madd->ncomp))
688     return(0);
689 greg 2.33 if (!sf) {
690 greg 2.1 mysf = (double *)malloc(sizeof(double)*msum->ncomp);
691 greg 2.33 if (!mysf)
692 greg 2.1 return(0);
693     for (k = msum->ncomp; k--; )
694     mysf[k] = 1;
695     sf = mysf;
696     }
697 greg 2.6 i = rmx_newtype(msum->dtype, madd->dtype);
698     if (i)
699     msum->dtype = i;
700     else
701     rmx_addinfo(msum, rmx_mismatch_warn);
702 greg 2.1 for (i = msum->nrows; i--; )
703 greg 2.53 for (j = msum->ncols; j--; ) {
704     const double *da = rmx_lval(madd,i,j);
705     double *ds = rmx_lval(msum,i,j);
706 greg 2.1 for (k = msum->ncomp; k--; )
707 greg 2.53 ds[k] += sf[k] * da[k];
708     }
709 greg 2.33 if (mysf)
710     free(mysf);
711 greg 2.1 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 greg 2.33 if (!rm | !sf)
721 greg 2.1 return(0);
722     for (i = rm->nrows; i--; )
723 greg 2.53 for (j = rm->ncols; j--; ) {
724     double *dp = rmx_lval(rm,i,j);
725 greg 2.1 for (k = rm->ncomp; k--; )
726 greg 2.53 dp[k] *= sf[k];
727     }
728 greg 2.28 if (rm->info)
729     rmx_addinfo(rm, "Applied scalar\n");
730 greg 2.1 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 greg 2.33 if (!msrc | (n <= 0) | !cmat)
741 greg 2.1 return(NULL);
742     dnew = rmx_alloc(msrc->nrows, msrc->ncols, n);
743 greg 2.33 if (!dnew)
744 greg 2.1 return(NULL);
745 greg 2.28 if (msrc->info) {
746     char buf[128];
747 greg 2.38 sprintf(buf, "Applied %dx%d component transform\n",
748 greg 2.28 dnew->ncomp, msrc->ncomp);
749     rmx_addinfo(dnew, msrc->info);
750     rmx_addinfo(dnew, buf);
751     }
752 greg 2.6 dnew->dtype = msrc->dtype;
753 greg 2.1 for (i = dnew->nrows; i--; )
754 greg 2.53 for (j = dnew->ncols; j--; ) {
755     const double *ds = rmx_lval(msrc,i,j);
756 greg 2.1 for (kd = dnew->ncomp; kd--; ) {
757     double d = 0;
758     for (ks = msrc->ncomp; ks--; )
759 greg 2.53 d += cmat[kd*msrc->ncomp + ks] * ds[ks];
760     rmx_lval(dnew,i,j)[kd] = d;
761 greg 2.1 }
762 greg 2.53 }
763 greg 2.1 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 greg 2.33 if (!cm)
774 greg 2.1 return(NULL);
775     dnew = rmx_alloc(cm->nrows, cm->ncols, 3);
776 greg 2.33 if (!dnew)
777 greg 2.1 return(NULL);
778 greg 2.6 dnew->dtype = DTfloat;
779 greg 2.1 for (i = dnew->nrows; i--; )
780     for (j = dnew->ncols; j--; ) {
781     const COLORV *cv = cm_lval(cm,i,j);
782 greg 2.53 double *dp = rmx_lval(dnew,i,j);
783     dp[0] = cv[0];
784     dp[1] = cv[1];
785     dp[2] = cv[2];
786 greg 2.1 }
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 greg 2.52 if (!rm || !rm->mtx | ((rm->ncomp != 3) & (rm->ncomp != 1)))
798 greg 2.1 return(NULL);
799     cnew = cm_alloc(rm->nrows, rm->ncols);
800 greg 2.33 if (!cnew)
801 greg 2.1 return(NULL);
802     for (i = cnew->nrows; i--; )
803     for (j = cnew->ncols; j--; ) {
804 greg 2.53 const double *dp = rmx_lval(rm,i,j);
805     COLORV *cv = cm_lval(cnew,i,j);
806 greg 2.52 if (rm->ncomp == 1)
807 greg 2.53 setcolor(cv, dp[0], dp[0], dp[0]);
808     else
809     setcolor(cv, dp[0], dp[1], dp[2]);
810 greg 2.1 }
811     return(cnew);
812     }