ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmatrix.c
Revision: 2.53
Committed: Sat Mar 5 01:45:21 2022 UTC (2 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.52: +65 -52 lines
Log Message:
refactor(rmtxop): Changed matrix array access macros

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.53 static const char RCSid[] = "$Id: rmatrix.c,v 2.52 2022/03/04 17:17:28 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     long pos; /* map memory to file 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     }
218     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.53 double *dp = rmx_lval(rm,i,j);
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     setcolr(scan[j], dp[0], dp[1], dp[2]);
440     }
441 greg 2.27 if (fwritecolrs(scan, rm->ncols, fp) < 0) {
442 greg 2.1 free(scan);
443     return(0);
444     }
445     }
446     free(scan);
447     return(1);
448     }
449    
450 greg 2.48 /* Check if CIE XYZ primaries were specified */
451     static int
452     findCIEprims(const char *info)
453     {
454     RGBPRIMS prims;
455    
456     if (!info)
457     return(0);
458     info = strstr(info, PRIMARYSTR);
459     if (!info || !primsval(prims, info))
460     return(0);
461    
462     return((prims[RED][CIEX] > .99) & (prims[RED][CIEY] < .01) &&
463     (prims[GRN][CIEX] < .01) & (prims[GRN][CIEY] > .99) &&
464     (prims[BLU][CIEX] < .01) & (prims[BLU][CIEY] < .01));
465     }
466    
467 greg 2.6 /* Write matrix to file type indicated by dtype */
468 greg 2.10 int
469 greg 2.1 rmx_write(const RMATRIX *rm, int dtype, FILE *fp)
470     {
471     int ok = 1;
472    
473 greg 2.52 if (!rm | !fp || !rm->mtx)
474 greg 2.1 return(0);
475 greg 2.17 #ifdef getc_unlocked
476     flockfile(fp);
477     #endif
478 greg 2.1 /* complete header */
479 greg 2.5 if (rm->info)
480     fputs(rm->info, fp);
481 greg 2.6 if (dtype == DTfromHeader)
482     dtype = rm->dtype;
483 greg 2.48 else if (dtype == DTrgbe && (rm->dtype == DTxyze ||
484     findCIEprims(rm->info)))
485 greg 2.6 dtype = DTxyze;
486 greg 2.9 else if ((dtype == DTxyze) & (rm->dtype == DTrgbe))
487 greg 2.6 dtype = DTrgbe;
488 greg 2.1 if ((dtype != DTrgbe) & (dtype != DTxyze)) {
489     fprintf(fp, "NROWS=%d\n", rm->nrows);
490     fprintf(fp, "NCOLS=%d\n", rm->ncols);
491     fprintf(fp, "NCOMP=%d\n", rm->ncomp);
492     } else if (rm->ncomp != 3) { /* wrong # components? */
493 greg 2.52 CMATRIX *cm; /* convert & write */
494 greg 2.1 if (rm->ncomp != 1) /* only convert grayscale */
495     return(0);
496 greg 2.52 if (!(cm = cm_from_rmatrix(rm)))
497 greg 2.1 return(0);
498 greg 2.52 fputformat(cm_fmt_id[dtype], fp);
499     fputc('\n', fp);
500     ok = cm_write(cm, dtype, fp);
501     cm_free(cm);
502     return(ok);
503 greg 2.1 }
504 greg 2.35 if ((dtype == DTfloat) | (dtype == DTdouble))
505     fputendian(fp); /* important to record */
506 greg 2.49 fputformat(cm_fmt_id[dtype], fp);
507 greg 2.1 fputc('\n', fp);
508     switch (dtype) { /* write data */
509     case DTascii:
510     ok = rmx_write_ascii(rm, fp);
511     break;
512     case DTfloat:
513     ok = rmx_write_float(rm, fp);
514     break;
515     case DTdouble:
516     ok = rmx_write_double(rm, fp);
517     break;
518     case DTrgbe:
519     case DTxyze:
520     fprtresolu(rm->ncols, rm->nrows, fp);
521     ok = rmx_write_rgbe(rm, fp);
522     break;
523     default:
524     return(0);
525     }
526     ok &= (fflush(fp) == 0);
527 greg 2.17 #ifdef getc_unlocked
528     funlockfile(fp);
529     #endif
530 greg 2.10 return(ok);
531 greg 2.1 }
532    
533     /* Allocate and assign square identity matrix with n components */
534     RMATRIX *
535     rmx_identity(const int dim, const int n)
536     {
537     RMATRIX *rid = rmx_alloc(dim, dim, n);
538 greg 2.12 int i, k;
539 greg 2.1
540 greg 2.33 if (!rid)
541 greg 2.1 return(NULL);
542 greg 2.52 memset(rid->mtx, 0, array_size(rid));
543 greg 2.53 for (i = dim; i--; ) {
544     double *dp = rmx_lval(rid,i,i);
545 greg 2.12 for (k = n; k--; )
546 greg 2.53 dp[k] = 1.;
547     }
548 greg 2.1 return(rid);
549     }
550    
551     /* Duplicate the given matrix */
552     RMATRIX *
553     rmx_copy(const RMATRIX *rm)
554     {
555     RMATRIX *dnew;
556    
557 greg 2.33 if (!rm)
558 greg 2.1 return(NULL);
559     dnew = rmx_alloc(rm->nrows, rm->ncols, rm->ncomp);
560 greg 2.33 if (!dnew)
561 greg 2.1 return(NULL);
562 greg 2.5 rmx_addinfo(dnew, rm->info);
563 greg 2.6 dnew->dtype = rm->dtype;
564 greg 2.52 memcpy(dnew->mtx, rm->mtx, array_size(dnew));
565 greg 2.1 return(dnew);
566     }
567    
568 greg 2.2 /* Allocate and assign transposed matrix */
569     RMATRIX *
570     rmx_transpose(const RMATRIX *rm)
571 greg 2.1 {
572 greg 2.2 RMATRIX *dnew;
573 greg 2.1 int i, j, k;
574    
575 greg 2.33 if (!rm)
576 greg 2.1 return(0);
577 greg 2.31 if ((rm->nrows == 1) | (rm->ncols == 1)) {
578     dnew = rmx_copy(rm);
579 greg 2.33 if (!dnew)
580 greg 2.32 return(NULL);
581 greg 2.31 dnew->nrows = rm->ncols;
582     dnew->ncols = rm->nrows;
583     return(dnew);
584     }
585 greg 2.2 dnew = rmx_alloc(rm->ncols, rm->nrows, rm->ncomp);
586 greg 2.33 if (!dnew)
587 greg 2.2 return(NULL);
588 greg 2.5 if (rm->info) {
589     rmx_addinfo(dnew, rm->info);
590     rmx_addinfo(dnew, "Transposed rows and columns\n");
591     }
592 greg 2.6 dnew->dtype = rm->dtype;
593 greg 2.2 for (i = dnew->nrows; i--; )
594     for (j = dnew->ncols; j--; )
595 greg 2.53 memcpy(rmx_lval(dnew,i,j), rmx_lval(rm,i,j),
596     sizeof(double)*dnew->ncomp);
597 greg 2.2 return(dnew);
598 greg 2.1 }
599    
600     /* Multiply (concatenate) two matrices and allocate the result */
601     RMATRIX *
602     rmx_multiply(const RMATRIX *m1, const RMATRIX *m2)
603     {
604     RMATRIX *mres;
605     int i, j, k, h;
606    
607 greg 2.33 if (!m1 | !m2 || (m1->ncomp != m2->ncomp) | (m1->ncols != m2->nrows))
608 greg 2.1 return(NULL);
609     mres = rmx_alloc(m1->nrows, m2->ncols, m1->ncomp);
610 greg 2.33 if (!mres)
611 greg 2.1 return(NULL);
612 greg 2.6 i = rmx_newtype(m1->dtype, m2->dtype);
613     if (i)
614     mres->dtype = i;
615     else
616     rmx_addinfo(mres, rmx_mismatch_warn);
617 greg 2.1 for (i = mres->nrows; i--; )
618     for (j = mres->ncols; j--; )
619 greg 2.8 for (k = mres->ncomp; k--; ) {
620 greg 2.53 double d = 0;
621 greg 2.8 for (h = m1->ncols; h--; )
622 greg 2.53 d += rmx_lval(m1,i,h)[k] * rmx_lval(m2,h,j)[k];
623     rmx_lval(mres,i,j)[k] = d;
624 greg 2.1 }
625     return(mres);
626     }
627    
628 greg 2.25 /* Element-wise multiplication (or division) of m2 into m1 */
629     int
630     rmx_elemult(RMATRIX *m1, const RMATRIX *m2, int divide)
631     {
632     int zeroDivides = 0;
633     int i, j, k;
634    
635 greg 2.33 if (!m1 | !m2 || (m1->ncols != m2->ncols) | (m1->nrows != m2->nrows))
636 greg 2.25 return(0);
637     if ((m2->ncomp > 1) & (m2->ncomp != m1->ncomp))
638     return(0);
639     i = rmx_newtype(m1->dtype, m2->dtype);
640     if (i)
641     m1->dtype = i;
642     else
643     rmx_addinfo(m1, rmx_mismatch_warn);
644     for (i = m1->nrows; i--; )
645     for (j = m1->ncols; j--; )
646     if (divide) {
647     double d;
648     if (m2->ncomp == 1) {
649 greg 2.53 d = rmx_lval(m2,i,j)[0];
650 greg 2.25 if (d == 0) {
651     ++zeroDivides;
652     for (k = m1->ncomp; k--; )
653 greg 2.53 rmx_lval(m1,i,j)[k] = 0;
654 greg 2.25 } else {
655     d = 1./d;
656     for (k = m1->ncomp; k--; )
657 greg 2.53 rmx_lval(m1,i,j)[k] *= d;
658 greg 2.25 }
659     } else
660     for (k = m1->ncomp; k--; ) {
661 greg 2.53 d = rmx_lval(m2,i,j)[k];
662 greg 2.25 if (d == 0) {
663     ++zeroDivides;
664 greg 2.53 rmx_lval(m1,i,j)[k] = 0;
665 greg 2.25 } else
666 greg 2.53 rmx_lval(m1,i,j)[k] /= d;
667 greg 2.25 }
668     } else {
669     if (m2->ncomp == 1) {
670 greg 2.53 const double d = rmx_lval(m2,i,j)[0];
671 greg 2.25 for (k = m1->ncomp; k--; )
672 greg 2.53 rmx_lval(m1,i,j)[k] *= d;
673 greg 2.25 } else
674     for (k = m1->ncomp; k--; )
675 greg 2.53 rmx_lval(m1,i,j)[k] *= rmx_lval(m2,i,j)[k];
676 greg 2.25 }
677     if (zeroDivides) {
678 greg 2.34 rmx_addinfo(m1, "WARNING: zero divide(s) corrupted results\n");
679 greg 2.25 errno = ERANGE;
680     }
681     return(1);
682     }
683    
684 greg 2.1 /* Sum second matrix into first, applying scale factor beforehand */
685     int
686     rmx_sum(RMATRIX *msum, const RMATRIX *madd, const double sf[])
687     {
688     double *mysf = NULL;
689     int i, j, k;
690    
691 greg 2.33 if (!msum | !madd ||
692 greg 2.1 (msum->nrows != madd->nrows) |
693     (msum->ncols != madd->ncols) |
694     (msum->ncomp != madd->ncomp))
695     return(0);
696 greg 2.33 if (!sf) {
697 greg 2.1 mysf = (double *)malloc(sizeof(double)*msum->ncomp);
698 greg 2.33 if (!mysf)
699 greg 2.1 return(0);
700     for (k = msum->ncomp; k--; )
701     mysf[k] = 1;
702     sf = mysf;
703     }
704 greg 2.6 i = rmx_newtype(msum->dtype, madd->dtype);
705     if (i)
706     msum->dtype = i;
707     else
708     rmx_addinfo(msum, rmx_mismatch_warn);
709 greg 2.1 for (i = msum->nrows; i--; )
710 greg 2.53 for (j = msum->ncols; j--; ) {
711     const double *da = rmx_lval(madd,i,j);
712     double *ds = rmx_lval(msum,i,j);
713 greg 2.1 for (k = msum->ncomp; k--; )
714 greg 2.53 ds[k] += sf[k] * da[k];
715     }
716 greg 2.33 if (mysf)
717     free(mysf);
718 greg 2.1 return(1);
719     }
720    
721     /* Scale the given matrix by the indicated scalar component vector */
722     int
723     rmx_scale(RMATRIX *rm, const double sf[])
724     {
725     int i, j, k;
726    
727 greg 2.33 if (!rm | !sf)
728 greg 2.1 return(0);
729     for (i = rm->nrows; i--; )
730 greg 2.53 for (j = rm->ncols; j--; ) {
731     double *dp = rmx_lval(rm,i,j);
732 greg 2.1 for (k = rm->ncomp; k--; )
733 greg 2.53 dp[k] *= sf[k];
734     }
735 greg 2.28 if (rm->info)
736     rmx_addinfo(rm, "Applied scalar\n");
737 greg 2.1 return(1);
738     }
739    
740     /* Allocate new matrix and apply component transformation */
741     RMATRIX *
742     rmx_transform(const RMATRIX *msrc, int n, const double cmat[])
743     {
744     int i, j, ks, kd;
745     RMATRIX *dnew;
746    
747 greg 2.33 if (!msrc | (n <= 0) | !cmat)
748 greg 2.1 return(NULL);
749     dnew = rmx_alloc(msrc->nrows, msrc->ncols, n);
750 greg 2.33 if (!dnew)
751 greg 2.1 return(NULL);
752 greg 2.28 if (msrc->info) {
753     char buf[128];
754 greg 2.38 sprintf(buf, "Applied %dx%d component transform\n",
755 greg 2.28 dnew->ncomp, msrc->ncomp);
756     rmx_addinfo(dnew, msrc->info);
757     rmx_addinfo(dnew, buf);
758     }
759 greg 2.6 dnew->dtype = msrc->dtype;
760 greg 2.1 for (i = dnew->nrows; i--; )
761 greg 2.53 for (j = dnew->ncols; j--; ) {
762     const double *ds = rmx_lval(msrc,i,j);
763 greg 2.1 for (kd = dnew->ncomp; kd--; ) {
764     double d = 0;
765     for (ks = msrc->ncomp; ks--; )
766 greg 2.53 d += cmat[kd*msrc->ncomp + ks] * ds[ks];
767     rmx_lval(dnew,i,j)[kd] = d;
768 greg 2.1 }
769 greg 2.53 }
770 greg 2.1 return(dnew);
771     }
772    
773     /* Convert a color matrix to newly allocated RMATRIX buffer */
774     RMATRIX *
775     rmx_from_cmatrix(const CMATRIX *cm)
776     {
777     int i, j;
778     RMATRIX *dnew;
779    
780 greg 2.33 if (!cm)
781 greg 2.1 return(NULL);
782     dnew = rmx_alloc(cm->nrows, cm->ncols, 3);
783 greg 2.33 if (!dnew)
784 greg 2.1 return(NULL);
785 greg 2.6 dnew->dtype = DTfloat;
786 greg 2.1 for (i = dnew->nrows; i--; )
787     for (j = dnew->ncols; j--; ) {
788     const COLORV *cv = cm_lval(cm,i,j);
789 greg 2.53 double *dp = rmx_lval(dnew,i,j);
790     dp[0] = cv[0];
791     dp[1] = cv[1];
792     dp[2] = cv[2];
793 greg 2.1 }
794     return(dnew);
795     }
796    
797     /* Convert general matrix to newly allocated CMATRIX buffer */
798     CMATRIX *
799     cm_from_rmatrix(const RMATRIX *rm)
800     {
801     int i, j;
802     CMATRIX *cnew;
803    
804 greg 2.52 if (!rm || !rm->mtx | ((rm->ncomp != 3) & (rm->ncomp != 1)))
805 greg 2.1 return(NULL);
806     cnew = cm_alloc(rm->nrows, rm->ncols);
807 greg 2.33 if (!cnew)
808 greg 2.1 return(NULL);
809     for (i = cnew->nrows; i--; )
810     for (j = cnew->ncols; j--; ) {
811 greg 2.53 const double *dp = rmx_lval(rm,i,j);
812     COLORV *cv = cm_lval(cnew,i,j);
813 greg 2.52 if (rm->ncomp == 1)
814 greg 2.53 setcolor(cv, dp[0], dp[0], dp[0]);
815     else
816     setcolor(cv, dp[0], dp[1], dp[2]);
817 greg 2.1 }
818     return(cnew);
819     }