ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmatrix.c
Revision: 2.50
Committed: Fri Mar 4 01:27:12 2022 UTC (2 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.49: +95 -53 lines
Log Message:
perf(rmtxop): Added memory-mapping for double matrices

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.50 static const char RCSid[] = "$Id: rmatrix.c,v 2.49 2022/03/03 03:55: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     double d = exposval(s);
141 greg 2.50 scalecolor(ip->cexp, d);
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     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 greg 2.50 if (!rmx_prepare(rm))
188     return(0);
189 greg 2.1 for (i = 0; i < rm->nrows; i++)
190     for (j = 0; j < rm->ncols; j++) {
191 greg 2.21 if (getbinary(val, sizeof(val[0]), rm->ncomp, fp) != rm->ncomp)
192 greg 2.1 return(0);
193 greg 2.35 if (rm->swapin)
194     swap32((char *)val, rm->ncomp);
195 greg 2.1 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 greg 2.41 int i;
205 greg 2.50 #ifdef MAP_FILE
206     long pos; /* map memory to file if possible */
207     if (!rm->swapin && (pos = ftell(fp)) >= 0 && !(pos % sizeof(double))) {
208     rm->mapped = mmap(NULL, array_size(rm)+pos, PROT_READ|PROT_WRITE,
209     MAP_PRIVATE, fileno(fp), 0);
210     if (rm->mapped != MAP_FAILED) {
211     rm->mtx = (double *)rm->mapped + pos/sizeof(double);
212     return(1);
213     }
214     rm->mapped = NULL;
215 greg 2.41 }
216 greg 2.50 #endif
217     if (!rmx_prepare(rm))
218     return(0);
219 greg 2.41 for (i = 0; i < rm->nrows; i++) {
220     if (getbinary(&rmx_lval(rm,i,0,0), sizeof(double)*rm->ncomp,
221     rm->ncols, fp) != rm->ncols)
222     return(0);
223 greg 2.35 if (rm->swapin)
224 greg 2.41 swap64((char *)&rmx_lval(rm,i,0,0), rm->ncols*rm->ncomp);
225     }
226 greg 2.1 return(1);
227     }
228    
229     static int
230     rmx_load_rgbe(RMATRIX *rm, FILE *fp)
231     {
232     COLOR *scan = (COLOR *)malloc(sizeof(COLOR)*rm->ncols);
233     int i, j;
234    
235 greg 2.33 if (!scan)
236 greg 2.1 return(0);
237 greg 2.50 if (!rmx_prepare(rm))
238     return(0);
239 greg 2.1 for (i = 0; i < rm->nrows; i++) {
240     if (freadscan(scan, rm->ncols, fp) < 0) {
241     free(scan);
242     return(0);
243     }
244     for (j = rm->ncols; j--; ) {
245     rmx_lval(rm,i,j,0) = colval(scan[j],RED);
246     rmx_lval(rm,i,j,1) = colval(scan[j],GRN);
247     rmx_lval(rm,i,j,2) = colval(scan[j],BLU);
248     }
249     }
250     free(scan);
251     return(1);
252     }
253    
254     /* Load matrix from supported file type */
255     RMATRIX *
256 greg 2.46 rmx_load(const char *inspec, RMPref rmp)
257 greg 2.1 {
258 greg 2.46 FILE *fp;
259 greg 2.1 RMATRIX *dnew;
260    
261 greg 2.47 if (!inspec)
262     inspec = stdin_name;
263     else if (!*inspec)
264 greg 2.46 return(NULL);
265     if (inspec == stdin_name) { /* reading from stdin? */
266     fp = stdin;
267 greg 2.13 } else if (inspec[0] == '!') {
268 greg 2.33 if (!(fp = popen(inspec+1, "r")))
269 greg 2.13 return(NULL);
270 greg 2.1 } else {
271 greg 2.13 const char *sp = inspec; /* check suffix */
272 greg 2.1 while (*sp)
273     ++sp;
274 greg 2.13 while (sp > inspec && sp[-1] != '.')
275 greg 2.1 --sp;
276     if (!strcasecmp(sp, "XML")) { /* assume it's a BSDF */
277 greg 2.46 CMATRIX *cm = rmp==RMPtrans ? cm_loadBTDF(inspec) :
278     cm_loadBRDF(inspec, rmp==RMPreflB) ;
279 greg 2.33 if (!cm)
280 greg 2.1 return(NULL);
281     dnew = rmx_from_cmatrix(cm);
282     cm_free(cm);
283 greg 2.18 dnew->dtype = DTascii;
284 greg 2.1 return(dnew);
285     }
286     /* else open it ourselves */
287 greg 2.46 if (!(fp = fopen(inspec, "r")))
288 greg 2.1 return(NULL);
289     }
290 greg 2.46 SET_FILE_BINARY(fp);
291 greg 2.1 #ifdef getc_unlocked
292     flockfile(fp);
293     #endif
294 greg 2.50 if (!(dnew = rmx_new(0,0,3))) {
295 greg 2.1 fclose(fp);
296     return(NULL);
297     }
298 greg 2.50 dnew->dtype = DTascii; /* assumed w/o FORMAT */
299     dnew->cexp[0] = dnew->cexp[1] = dnew->cexp[2] = 1.f;
300     if (getheader(fp, get_dminfo, dnew) < 0) {
301     fclose(fp);
302     return(NULL);
303     }
304     if ((dnew->nrows <= 0) | (dnew->ncols <= 0)) {
305     if (!fscnresolu(&dnew->ncols, &dnew->nrows, fp)) {
306 greg 2.1 fclose(fp);
307     return(NULL);
308     }
309 greg 2.50 if ((dnew->dtype == DTrgbe) | (dnew->dtype == DTxyze) &&
310     dnew->ncomp != 3) {
311 greg 2.4 fclose(fp);
312     return(NULL);
313     }
314 greg 2.1 }
315 greg 2.50 switch (dnew->dtype) {
316 greg 2.1 case DTascii:
317 greg 2.37 SET_FILE_TEXT(fp);
318 greg 2.1 if (!rmx_load_ascii(dnew, fp))
319     goto loaderr;
320 greg 2.11 dnew->dtype = DTascii; /* should leave double? */
321 greg 2.1 break;
322     case DTfloat:
323     if (!rmx_load_float(dnew, fp))
324     goto loaderr;
325 greg 2.6 dnew->dtype = DTfloat;
326 greg 2.1 break;
327     case DTdouble:
328     if (!rmx_load_double(dnew, fp))
329     goto loaderr;
330 greg 2.6 dnew->dtype = DTdouble;
331 greg 2.1 break;
332     case DTrgbe:
333     case DTxyze:
334     if (!rmx_load_rgbe(dnew, fp))
335     goto loaderr;
336 greg 2.50 /* undo exposure? */
337     if ((dnew->cexp[0] != 1.f) | (dnew->cexp[1] != 1.f) |
338     (dnew->cexp[2] != 1.f)) {
339     double cmlt[3];
340     cmlt[0] = 1./dnew->cexp[0];
341     cmlt[1] = 1./dnew->cexp[1];
342     cmlt[2] = 1./dnew->cexp[2];
343     rmx_scale(dnew, cmlt);
344 greg 2.40 }
345 greg 2.50 dnew->swapin = 0;
346 greg 2.1 break;
347     default:
348     goto loaderr;
349     }
350 greg 2.13 if (fp != stdin) {
351     if (inspec[0] == '!')
352     pclose(fp);
353     else
354     fclose(fp);
355     }
356     #ifdef getc_unlocked
357     else
358     funlockfile(fp);
359     #endif
360 greg 2.1 return(dnew);
361     loaderr: /* should report error? */
362 greg 2.13 if (inspec[0] == '!')
363     pclose(fp);
364     else
365     fclose(fp);
366 greg 2.1 rmx_free(dnew);
367     return(NULL);
368     }
369    
370     static int
371     rmx_write_ascii(const RMATRIX *rm, FILE *fp)
372     {
373 greg 2.30 const char *fmt = (rm->dtype == DTfloat) ? " %.7e" :
374     (rm->dtype == DTrgbe) | (rm->dtype == DTxyze) ? " %.3e" :
375     " %.15e" ;
376 greg 2.1 int i, j, k;
377 greg 2.17
378 greg 2.1 for (i = 0; i < rm->nrows; i++) {
379     for (j = 0; j < rm->ncols; j++) {
380     for (k = 0; k < rm->ncomp; k++)
381 greg 2.30 fprintf(fp, fmt, rmx_lval(rm,i,j,k));
382 greg 2.1 fputc('\t', fp);
383     }
384     fputc('\n', fp);
385     }
386     return(1);
387     }
388    
389     static int
390     rmx_write_float(const RMATRIX *rm, FILE *fp)
391     {
392     int i, j, k;
393     float val[100];
394    
395     if (rm->ncomp > 100) {
396     fputs("Unsupported # components in rmx_write_float()\n", stderr);
397     exit(1);
398     }
399     for (i = 0; i < rm->nrows; i++)
400     for (j = 0; j < rm->ncols; j++) {
401     for (k = rm->ncomp; k--; )
402     val[k] = (float)rmx_lval(rm,i,j,k);
403 greg 2.21 if (putbinary(val, sizeof(val[0]), rm->ncomp, fp) != rm->ncomp)
404 greg 2.1 return(0);
405     }
406     return(1);
407     }
408    
409     static int
410     rmx_write_double(const RMATRIX *rm, FILE *fp)
411     {
412 greg 2.26 int i, j;
413 greg 2.1
414     for (i = 0; i < rm->nrows; i++)
415 greg 2.26 for (j = 0; j < rm->ncols; j++)
416     if (putbinary(&rmx_lval(rm,i,j,0), sizeof(double), rm->ncomp, fp) != rm->ncomp)
417 greg 2.1 return(0);
418     return(1);
419     }
420    
421     static int
422     rmx_write_rgbe(const RMATRIX *rm, FILE *fp)
423     {
424 greg 2.27 COLR *scan = (COLR *)malloc(sizeof(COLR)*rm->ncols);
425 greg 2.1 int i, j;
426    
427 greg 2.33 if (!scan)
428 greg 2.1 return(0);
429     for (i = 0; i < rm->nrows; i++) {
430     for (j = rm->ncols; j--; )
431 greg 2.27 setcolr(scan[j], rmx_lval(rm,i,j,0),
432 greg 2.1 rmx_lval(rm,i,j,1),
433     rmx_lval(rm,i,j,2) );
434 greg 2.27 if (fwritecolrs(scan, rm->ncols, fp) < 0) {
435 greg 2.1 free(scan);
436     return(0);
437     }
438     }
439     free(scan);
440     return(1);
441     }
442    
443 greg 2.48 /* Check if CIE XYZ primaries were specified */
444     static int
445     findCIEprims(const char *info)
446     {
447     RGBPRIMS prims;
448    
449     if (!info)
450     return(0);
451     info = strstr(info, PRIMARYSTR);
452     if (!info || !primsval(prims, info))
453     return(0);
454    
455     return((prims[RED][CIEX] > .99) & (prims[RED][CIEY] < .01) &&
456     (prims[GRN][CIEX] < .01) & (prims[GRN][CIEY] > .99) &&
457     (prims[BLU][CIEX] < .01) & (prims[BLU][CIEY] < .01));
458     }
459    
460 greg 2.6 /* Write matrix to file type indicated by dtype */
461 greg 2.10 int
462 greg 2.1 rmx_write(const RMATRIX *rm, int dtype, FILE *fp)
463     {
464     RMATRIX *mydm = NULL;
465     int ok = 1;
466    
467 greg 2.33 if (!rm | !fp)
468 greg 2.1 return(0);
469 greg 2.17 #ifdef getc_unlocked
470     flockfile(fp);
471     #endif
472 greg 2.1 /* complete header */
473 greg 2.5 if (rm->info)
474     fputs(rm->info, fp);
475 greg 2.6 if (dtype == DTfromHeader)
476     dtype = rm->dtype;
477 greg 2.48 else if (dtype == DTrgbe && (rm->dtype == DTxyze ||
478     findCIEprims(rm->info)))
479 greg 2.6 dtype = DTxyze;
480 greg 2.9 else if ((dtype == DTxyze) & (rm->dtype == DTrgbe))
481 greg 2.6 dtype = DTrgbe;
482 greg 2.1 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     double cmtx[3];
488     if (rm->ncomp != 1) /* only convert grayscale */
489     return(0);
490     cmtx[0] = cmtx[1] = cmtx[2] = 1;
491     mydm = rmx_transform(rm, 3, cmtx);
492 greg 2.33 if (!mydm)
493 greg 2.1 return(0);
494     rm = mydm;
495     }
496 greg 2.35 if ((dtype == DTfloat) | (dtype == DTdouble))
497     fputendian(fp); /* important to record */
498 greg 2.49 fputformat(cm_fmt_id[dtype], fp);
499 greg 2.1 fputc('\n', fp);
500     switch (dtype) { /* write data */
501     case DTascii:
502     ok = rmx_write_ascii(rm, fp);
503     break;
504     case DTfloat:
505     ok = rmx_write_float(rm, fp);
506     break;
507     case DTdouble:
508     ok = rmx_write_double(rm, fp);
509     break;
510     case DTrgbe:
511     case DTxyze:
512     fprtresolu(rm->ncols, rm->nrows, fp);
513     ok = rmx_write_rgbe(rm, fp);
514     break;
515     default:
516     return(0);
517     }
518     ok &= (fflush(fp) == 0);
519 greg 2.17 #ifdef getc_unlocked
520     funlockfile(fp);
521     #endif
522 greg 2.33 if (mydm)
523     rmx_free(mydm);
524 greg 2.10 return(ok);
525 greg 2.1 }
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 greg 2.12 int i, k;
533 greg 2.1
534 greg 2.33 if (!rid)
535 greg 2.1 return(NULL);
536 greg 2.12 memset(rid->mtx, 0, sizeof(rid->mtx[0])*n*dim*dim);
537 greg 2.1 for (i = dim; i--; )
538 greg 2.12 for (k = n; k--; )
539     rmx_lval(rid,i,i,k) = 1;
540 greg 2.1 return(rid);
541     }
542    
543     /* Duplicate the given matrix */
544     RMATRIX *
545     rmx_copy(const RMATRIX *rm)
546     {
547     RMATRIX *dnew;
548    
549 greg 2.33 if (!rm)
550 greg 2.1 return(NULL);
551     dnew = rmx_alloc(rm->nrows, rm->ncols, rm->ncomp);
552 greg 2.33 if (!dnew)
553 greg 2.1 return(NULL);
554 greg 2.5 rmx_addinfo(dnew, rm->info);
555 greg 2.6 dnew->dtype = rm->dtype;
556 greg 2.1 memcpy(dnew->mtx, rm->mtx,
557     sizeof(rm->mtx[0])*rm->ncomp*rm->nrows*rm->ncols);
558     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.1 int i, j, k;
567    
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.2 for (i = dnew->nrows; i--; )
587     for (j = dnew->ncols; j--; )
588     for (k = dnew->ncomp; k--; )
589     rmx_lval(dnew,i,j,k) = rmx_lval(rm,j,i,k);
590     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.1 long double d = 0;
614 greg 2.8 for (h = m1->ncols; h--; )
615 greg 2.24 d += rmx_lval(m1,i,h,k) * rmx_lval(m2,h,j,k);
616 greg 2.1 rmx_lval(mres,i,j,k) = (double)d;
617     }
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     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 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     for (j = msum->ncols; j--; )
704     for (k = msum->ncomp; k--; )
705     rmx_lval(msum,i,j,k) += sf[k] * rmx_lval(madd,i,j,k);
706 greg 2.33 if (mysf)
707     free(mysf);
708 greg 2.1 return(1);
709     }
710    
711     /* Scale the given matrix by the indicated scalar component vector */
712     int
713     rmx_scale(RMATRIX *rm, const double sf[])
714     {
715     int i, j, k;
716    
717 greg 2.33 if (!rm | !sf)
718 greg 2.1 return(0);
719     for (i = rm->nrows; i--; )
720     for (j = rm->ncols; j--; )
721     for (k = rm->ncomp; k--; )
722     rmx_lval(rm,i,j,k) *= sf[k];
723    
724 greg 2.28 if (rm->info)
725     rmx_addinfo(rm, "Applied scalar\n");
726 greg 2.1 return(1);
727     }
728    
729     /* Allocate new matrix and apply component transformation */
730     RMATRIX *
731     rmx_transform(const RMATRIX *msrc, int n, const double cmat[])
732     {
733     int i, j, ks, kd;
734     RMATRIX *dnew;
735    
736 greg 2.33 if (!msrc | (n <= 0) | !cmat)
737 greg 2.1 return(NULL);
738     dnew = rmx_alloc(msrc->nrows, msrc->ncols, n);
739 greg 2.33 if (!dnew)
740 greg 2.1 return(NULL);
741 greg 2.28 if (msrc->info) {
742     char buf[128];
743 greg 2.38 sprintf(buf, "Applied %dx%d component transform\n",
744 greg 2.28 dnew->ncomp, msrc->ncomp);
745     rmx_addinfo(dnew, msrc->info);
746     rmx_addinfo(dnew, buf);
747     }
748 greg 2.6 dnew->dtype = msrc->dtype;
749 greg 2.1 for (i = dnew->nrows; i--; )
750     for (j = dnew->ncols; j--; )
751     for (kd = dnew->ncomp; kd--; ) {
752     double d = 0;
753     for (ks = msrc->ncomp; ks--; )
754     d += cmat[kd*msrc->ncomp + ks] * rmx_lval(msrc,i,j,ks);
755     rmx_lval(dnew,i,j,kd) = d;
756     }
757     return(dnew);
758     }
759    
760     /* Convert a color matrix to newly allocated RMATRIX buffer */
761     RMATRIX *
762     rmx_from_cmatrix(const CMATRIX *cm)
763     {
764     int i, j;
765     RMATRIX *dnew;
766    
767 greg 2.33 if (!cm)
768 greg 2.1 return(NULL);
769     dnew = rmx_alloc(cm->nrows, cm->ncols, 3);
770 greg 2.33 if (!dnew)
771 greg 2.1 return(NULL);
772 greg 2.6 dnew->dtype = DTfloat;
773 greg 2.1 for (i = dnew->nrows; i--; )
774     for (j = dnew->ncols; j--; ) {
775     const COLORV *cv = cm_lval(cm,i,j);
776     rmx_lval(dnew,i,j,0) = cv[0];
777     rmx_lval(dnew,i,j,1) = cv[1];
778     rmx_lval(dnew,i,j,2) = cv[2];
779     }
780     return(dnew);
781     }
782    
783     /* Convert general matrix to newly allocated CMATRIX buffer */
784     CMATRIX *
785     cm_from_rmatrix(const RMATRIX *rm)
786     {
787     int i, j;
788     CMATRIX *cnew;
789    
790 greg 2.33 if (!rm || rm->ncomp != 3)
791 greg 2.1 return(NULL);
792     cnew = cm_alloc(rm->nrows, rm->ncols);
793 greg 2.33 if (!cnew)
794 greg 2.1 return(NULL);
795     for (i = cnew->nrows; i--; )
796     for (j = cnew->ncols; j--; ) {
797     COLORV *cv = cm_lval(cnew,i,j);
798     cv[0] = (COLORV)rmx_lval(rm,i,j,0);
799     cv[1] = (COLORV)rmx_lval(rm,i,j,1);
800     cv[2] = (COLORV)rmx_lval(rm,i,j,2);
801     }
802     return(cnew);
803     }