ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmatrix.c
Revision: 2.43
Committed: Thu Mar 26 18:05:44 2020 UTC (4 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.42: +3 -2 lines
Log Message:
Better fix to bug in last change

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.43 static const char RCSid[] = "$Id: rmatrix.c,v 2.42 2020/03/26 18:04:24 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    
16 greg 2.6 static char rmx_mismatch_warn[] = "WARNING: data type mismatch\n";
17 greg 2.1
18     /* Allocate a nr x nc matrix with n components */
19     RMATRIX *
20     rmx_alloc(int nr, int nc, int n)
21     {
22     RMATRIX *dnew;
23    
24     if ((nr <= 0) | (nc <= 0) | (n <= 0))
25     return(NULL);
26     dnew = (RMATRIX *)malloc(sizeof(RMATRIX)-sizeof(dnew->mtx) +
27     sizeof(dnew->mtx[0])*(n*nr*nc));
28 greg 2.33 if (!dnew)
29 greg 2.1 return(NULL);
30     dnew->nrows = nr; dnew->ncols = nc; dnew->ncomp = n;
31 greg 2.6 dnew->dtype = DTdouble;
32 greg 2.36 dnew->swapin = 0;
33 greg 2.5 dnew->info = NULL;
34 greg 2.1 return(dnew);
35     }
36    
37 greg 2.6 /* Free a RMATRIX array */
38     void
39     rmx_free(RMATRIX *rm)
40     {
41     if (!rm) return;
42     if (rm->info)
43     free(rm->info);
44     free(rm);
45     }
46    
47     /* Resolve data type based on two input types (returns 0 for mismatch) */
48     int
49     rmx_newtype(int dtyp1, int dtyp2)
50     {
51 greg 2.14 if ((dtyp1==DTxyze) | (dtyp1==DTrgbe) |
52     (dtyp2==DTxyze) | (dtyp2==DTrgbe)
53     && dtyp1 != dtyp2)
54 greg 2.6 return(0);
55     if (dtyp1 < dtyp2)
56     return(dtyp1);
57     return(dtyp2);
58     }
59    
60 greg 2.5 /* Append header information associated with matrix data */
61     int
62     rmx_addinfo(RMATRIX *rm, const char *info)
63     {
64 greg 2.33 if (!rm || !info || !*info)
65 greg 2.5 return(0);
66 greg 2.8 if (!rm->info) {
67 greg 2.5 rm->info = (char *)malloc(strlen(info)+1);
68 greg 2.8 if (rm->info) rm->info[0] = '\0';
69     } else
70 greg 2.5 rm->info = (char *)realloc(rm->info,
71     strlen(rm->info)+strlen(info)+1);
72     if (!rm->info)
73     return(0);
74     strcat(rm->info, info);
75     return(1);
76     }
77    
78 greg 2.1 static int
79     get_dminfo(char *s, void *p)
80     {
81 greg 2.6 RMATRIX *ip = (RMATRIX *)p;
82 greg 2.29 char fmt[MAXFMTLEN];
83 greg 2.1 int i;
84    
85 greg 2.6 if (headidval(fmt, s))
86     return(0);
87 greg 2.1 if (!strncmp(s, "NCOMP=", 6)) {
88     ip->ncomp = atoi(s+6);
89     return(0);
90     }
91     if (!strncmp(s, "NROWS=", 6)) {
92     ip->nrows = atoi(s+6);
93     return(0);
94     }
95     if (!strncmp(s, "NCOLS=", 6)) {
96     ip->ncols = atoi(s+6);
97     return(0);
98     }
99 greg 2.35 if ((i = isbigendian(s)) >= 0) {
100     ip->swapin = (nativebigendian() != i);
101     return(0);
102     }
103 greg 2.40 if (isexpos(s)) {
104     double d = exposval(s);
105     scalecolor(ip->mtx, d);
106     return(0);
107     }
108     if (iscolcor(s)) {
109     COLOR ctmp;
110     colcorval(ctmp, s);
111     multcolor(ip->mtx, ctmp);
112     return(0);
113     }
114 greg 2.5 if (!formatval(fmt, s)) {
115 greg 2.6 rmx_addinfo(ip, s);
116 greg 2.1 return(0);
117 greg 2.5 }
118 greg 2.1 for (i = 1; i < DTend; i++)
119     if (!strcmp(fmt, cm_fmt_id[i])) {
120     ip->dtype = i;
121     return(0);
122     }
123     return(-1);
124     }
125    
126     static int
127     rmx_load_ascii(RMATRIX *rm, FILE *fp)
128     {
129     int i, j, k;
130 greg 2.17
131 greg 2.1 for (i = 0; i < rm->nrows; i++)
132     for (j = 0; j < rm->ncols; j++)
133     for (k = 0; k < rm->ncomp; k++)
134     if (fscanf(fp, "%lf", &rmx_lval(rm,i,j,k)) != 1)
135     return(0);
136     return(1);
137     }
138    
139     static int
140     rmx_load_float(RMATRIX *rm, FILE *fp)
141     {
142     int i, j, k;
143     float val[100];
144    
145     if (rm->ncomp > 100) {
146     fputs("Unsupported # components in rmx_load_float()\n", stderr);
147     exit(1);
148     }
149     for (i = 0; i < rm->nrows; i++)
150     for (j = 0; j < rm->ncols; j++) {
151 greg 2.21 if (getbinary(val, sizeof(val[0]), rm->ncomp, fp) != rm->ncomp)
152 greg 2.1 return(0);
153 greg 2.35 if (rm->swapin)
154     swap32((char *)val, rm->ncomp);
155 greg 2.1 for (k = rm->ncomp; k--; )
156     rmx_lval(rm,i,j,k) = val[k];
157     }
158     return(1);
159     }
160    
161     static int
162     rmx_load_double(RMATRIX *rm, FILE *fp)
163     {
164 greg 2.41 int i;
165 greg 2.1
166 greg 2.43 if ((char *)&rmx_lval(rm,0,0,0) - (char *)&rmx_lval(rm,1,0,0) !=
167     sizeof(double)*rm->ncols*rm->ncomp) {
168 greg 2.41 fputs("Code error in rmx_load_double()\n", stderr);
169     exit(1);
170     }
171     for (i = 0; i < rm->nrows; i++) {
172     if (getbinary(&rmx_lval(rm,i,0,0), sizeof(double)*rm->ncomp,
173     rm->ncols, fp) != rm->ncols)
174     return(0);
175 greg 2.35 if (rm->swapin)
176 greg 2.41 swap64((char *)&rmx_lval(rm,i,0,0), rm->ncols*rm->ncomp);
177     }
178 greg 2.1 return(1);
179     }
180    
181     static int
182     rmx_load_rgbe(RMATRIX *rm, FILE *fp)
183     {
184     COLOR *scan = (COLOR *)malloc(sizeof(COLOR)*rm->ncols);
185     int i, j;
186    
187 greg 2.33 if (!scan)
188 greg 2.1 return(0);
189     for (i = 0; i < rm->nrows; i++) {
190     if (freadscan(scan, rm->ncols, fp) < 0) {
191     free(scan);
192     return(0);
193     }
194     for (j = rm->ncols; j--; ) {
195     rmx_lval(rm,i,j,0) = colval(scan[j],RED);
196     rmx_lval(rm,i,j,1) = colval(scan[j],GRN);
197     rmx_lval(rm,i,j,2) = colval(scan[j],BLU);
198     }
199     }
200     free(scan);
201     return(1);
202     }
203    
204     /* Load matrix from supported file type */
205     RMATRIX *
206 greg 2.13 rmx_load(const char *inspec)
207 greg 2.1 {
208     FILE *fp = stdin;
209 greg 2.6 RMATRIX dinfo;
210 greg 2.1 RMATRIX *dnew;
211    
212 greg 2.33 if (!inspec) { /* reading from stdin? */
213 greg 2.13 inspec = "<stdin>";
214 greg 2.22 SET_FILE_BINARY(stdin);
215 greg 2.13 } else if (inspec[0] == '!') {
216 greg 2.33 if (!(fp = popen(inspec+1, "r")))
217 greg 2.13 return(NULL);
218 greg 2.37 SET_FILE_BINARY(fp);
219 greg 2.1 } else {
220 greg 2.13 const char *sp = inspec; /* check suffix */
221 greg 2.1 while (*sp)
222     ++sp;
223 greg 2.13 while (sp > inspec && sp[-1] != '.')
224 greg 2.1 --sp;
225     if (!strcasecmp(sp, "XML")) { /* assume it's a BSDF */
226 greg 2.13 CMATRIX *cm = cm_loadBTDF((char *)inspec);
227 greg 2.33 if (!cm)
228 greg 2.1 return(NULL);
229     dnew = rmx_from_cmatrix(cm);
230     cm_free(cm);
231 greg 2.18 dnew->dtype = DTascii;
232 greg 2.1 return(dnew);
233     }
234     /* else open it ourselves */
235 greg 2.33 if (!(fp = fopen(inspec, "rb")))
236 greg 2.1 return(NULL);
237     }
238     #ifdef getc_unlocked
239     flockfile(fp);
240     #endif
241     dinfo.nrows = dinfo.ncols = dinfo.ncomp = 0;
242 greg 2.6 dinfo.dtype = DTascii; /* assumed w/o FORMAT */
243 greg 2.35 dinfo.swapin = 0;
244 greg 2.6 dinfo.info = NULL;
245 greg 2.40 dinfo.mtx[0] = dinfo.mtx[1] = dinfo.mtx[2] = 1.;
246 greg 2.3 if (getheader(fp, get_dminfo, &dinfo) < 0) {
247 greg 2.1 fclose(fp);
248     return(NULL);
249     }
250 greg 2.4 if ((dinfo.nrows <= 0) | (dinfo.ncols <= 0)) {
251 greg 2.1 if (!fscnresolu(&dinfo.ncols, &dinfo.nrows, fp)) {
252     fclose(fp);
253     return(NULL);
254     }
255 greg 2.4 if (dinfo.ncomp <= 0)
256     dinfo.ncomp = 3;
257     else if ((dinfo.dtype == DTrgbe) | (dinfo.dtype == DTxyze) &&
258     dinfo.ncomp != 3) {
259     fclose(fp);
260     return(NULL);
261     }
262 greg 2.1 }
263     dnew = rmx_alloc(dinfo.nrows, dinfo.ncols, dinfo.ncomp);
264 greg 2.33 if (!dnew) {
265 greg 2.1 fclose(fp);
266     return(NULL);
267     }
268 greg 2.6 dnew->info = dinfo.info;
269 greg 2.1 switch (dinfo.dtype) {
270     case DTascii:
271 greg 2.37 SET_FILE_TEXT(fp);
272 greg 2.1 if (!rmx_load_ascii(dnew, fp))
273     goto loaderr;
274 greg 2.11 dnew->dtype = DTascii; /* should leave double? */
275 greg 2.1 break;
276     case DTfloat:
277 greg 2.35 dnew->swapin = dinfo.swapin;
278 greg 2.1 if (!rmx_load_float(dnew, fp))
279     goto loaderr;
280 greg 2.6 dnew->dtype = DTfloat;
281 greg 2.1 break;
282     case DTdouble:
283 greg 2.35 dnew->swapin = dinfo.swapin;
284 greg 2.1 if (!rmx_load_double(dnew, fp))
285     goto loaderr;
286 greg 2.6 dnew->dtype = DTdouble;
287 greg 2.1 break;
288     case DTrgbe:
289     case DTxyze:
290     if (!rmx_load_rgbe(dnew, fp))
291     goto loaderr;
292 greg 2.40 dnew->dtype = dinfo.dtype; /* undo exposure? */
293     if ((dinfo.mtx[0] != 1.) | (dinfo.mtx[1] != 1.) |
294     (dinfo.mtx[2] != 1.)) {
295     dinfo.mtx[0] = 1./dinfo.mtx[0];
296     dinfo.mtx[1] = 1./dinfo.mtx[1];
297     dinfo.mtx[2] = 1./dinfo.mtx[2];
298     rmx_scale(dnew, dinfo.mtx);
299     }
300 greg 2.1 break;
301     default:
302     goto loaderr;
303     }
304 greg 2.13 if (fp != stdin) {
305     if (inspec[0] == '!')
306     pclose(fp);
307     else
308     fclose(fp);
309     }
310     #ifdef getc_unlocked
311     else
312     funlockfile(fp);
313     #endif
314 greg 2.1 return(dnew);
315     loaderr: /* should report error? */
316 greg 2.13 if (inspec[0] == '!')
317     pclose(fp);
318     else
319     fclose(fp);
320 greg 2.1 rmx_free(dnew);
321     return(NULL);
322     }
323    
324     static int
325     rmx_write_ascii(const RMATRIX *rm, FILE *fp)
326     {
327 greg 2.30 const char *fmt = (rm->dtype == DTfloat) ? " %.7e" :
328     (rm->dtype == DTrgbe) | (rm->dtype == DTxyze) ? " %.3e" :
329     " %.15e" ;
330 greg 2.1 int i, j, k;
331 greg 2.17
332 greg 2.1 for (i = 0; i < rm->nrows; i++) {
333     for (j = 0; j < rm->ncols; j++) {
334     for (k = 0; k < rm->ncomp; k++)
335 greg 2.30 fprintf(fp, fmt, rmx_lval(rm,i,j,k));
336 greg 2.1 fputc('\t', fp);
337     }
338     fputc('\n', fp);
339     }
340     return(1);
341     }
342    
343     static int
344     rmx_write_float(const RMATRIX *rm, FILE *fp)
345     {
346     int i, j, k;
347     float val[100];
348    
349     if (rm->ncomp > 100) {
350     fputs("Unsupported # components in rmx_write_float()\n", stderr);
351     exit(1);
352     }
353     for (i = 0; i < rm->nrows; i++)
354     for (j = 0; j < rm->ncols; j++) {
355     for (k = rm->ncomp; k--; )
356     val[k] = (float)rmx_lval(rm,i,j,k);
357 greg 2.21 if (putbinary(val, sizeof(val[0]), rm->ncomp, fp) != rm->ncomp)
358 greg 2.1 return(0);
359     }
360     return(1);
361     }
362    
363     static int
364     rmx_write_double(const RMATRIX *rm, FILE *fp)
365     {
366 greg 2.26 int i, j;
367 greg 2.1
368     for (i = 0; i < rm->nrows; i++)
369 greg 2.26 for (j = 0; j < rm->ncols; j++)
370     if (putbinary(&rmx_lval(rm,i,j,0), sizeof(double), rm->ncomp, fp) != rm->ncomp)
371 greg 2.1 return(0);
372     return(1);
373     }
374    
375     static int
376     rmx_write_rgbe(const RMATRIX *rm, FILE *fp)
377     {
378 greg 2.27 COLR *scan = (COLR *)malloc(sizeof(COLR)*rm->ncols);
379 greg 2.1 int i, j;
380    
381 greg 2.33 if (!scan)
382 greg 2.1 return(0);
383     for (i = 0; i < rm->nrows; i++) {
384     for (j = rm->ncols; j--; )
385 greg 2.27 setcolr(scan[j], rmx_lval(rm,i,j,0),
386 greg 2.1 rmx_lval(rm,i,j,1),
387     rmx_lval(rm,i,j,2) );
388 greg 2.27 if (fwritecolrs(scan, rm->ncols, fp) < 0) {
389 greg 2.1 free(scan);
390     return(0);
391     }
392     }
393     free(scan);
394     return(1);
395     }
396    
397 greg 2.6 /* Write matrix to file type indicated by dtype */
398 greg 2.10 int
399 greg 2.1 rmx_write(const RMATRIX *rm, int dtype, FILE *fp)
400     {
401     RMATRIX *mydm = NULL;
402     int ok = 1;
403    
404 greg 2.33 if (!rm | !fp)
405 greg 2.1 return(0);
406 greg 2.17 #ifdef getc_unlocked
407     flockfile(fp);
408     #endif
409 greg 2.1 /* complete header */
410 greg 2.5 if (rm->info)
411     fputs(rm->info, fp);
412 greg 2.6 if (dtype == DTfromHeader)
413     dtype = rm->dtype;
414     else if ((dtype == DTrgbe) & (rm->dtype == DTxyze))
415     dtype = DTxyze;
416 greg 2.9 else if ((dtype == DTxyze) & (rm->dtype == DTrgbe))
417 greg 2.6 dtype = DTrgbe;
418 greg 2.1 if ((dtype != DTrgbe) & (dtype != DTxyze)) {
419     fprintf(fp, "NROWS=%d\n", rm->nrows);
420     fprintf(fp, "NCOLS=%d\n", rm->ncols);
421     fprintf(fp, "NCOMP=%d\n", rm->ncomp);
422     } else if (rm->ncomp != 3) { /* wrong # components? */
423     double cmtx[3];
424     if (rm->ncomp != 1) /* only convert grayscale */
425     return(0);
426     cmtx[0] = cmtx[1] = cmtx[2] = 1;
427     mydm = rmx_transform(rm, 3, cmtx);
428 greg 2.33 if (!mydm)
429 greg 2.1 return(0);
430     rm = mydm;
431     }
432 greg 2.35 if ((dtype == DTfloat) | (dtype == DTdouble))
433     fputendian(fp); /* important to record */
434 greg 2.1 fputformat((char *)cm_fmt_id[dtype], fp);
435     fputc('\n', fp);
436     switch (dtype) { /* write data */
437     case DTascii:
438     ok = rmx_write_ascii(rm, fp);
439     break;
440     case DTfloat:
441     ok = rmx_write_float(rm, fp);
442     break;
443     case DTdouble:
444     ok = rmx_write_double(rm, fp);
445     break;
446     case DTrgbe:
447     case DTxyze:
448     fprtresolu(rm->ncols, rm->nrows, fp);
449     ok = rmx_write_rgbe(rm, fp);
450     break;
451     default:
452     return(0);
453     }
454     ok &= (fflush(fp) == 0);
455 greg 2.17 #ifdef getc_unlocked
456     funlockfile(fp);
457     #endif
458 greg 2.33 if (mydm)
459     rmx_free(mydm);
460 greg 2.10 return(ok);
461 greg 2.1 }
462    
463     /* Allocate and assign square identity matrix with n components */
464     RMATRIX *
465     rmx_identity(const int dim, const int n)
466     {
467     RMATRIX *rid = rmx_alloc(dim, dim, n);
468 greg 2.12 int i, k;
469 greg 2.1
470 greg 2.33 if (!rid)
471 greg 2.1 return(NULL);
472 greg 2.12 memset(rid->mtx, 0, sizeof(rid->mtx[0])*n*dim*dim);
473 greg 2.1 for (i = dim; i--; )
474 greg 2.12 for (k = n; k--; )
475     rmx_lval(rid,i,i,k) = 1;
476 greg 2.1 return(rid);
477     }
478    
479     /* Duplicate the given matrix */
480     RMATRIX *
481     rmx_copy(const RMATRIX *rm)
482     {
483     RMATRIX *dnew;
484    
485 greg 2.33 if (!rm)
486 greg 2.1 return(NULL);
487     dnew = rmx_alloc(rm->nrows, rm->ncols, rm->ncomp);
488 greg 2.33 if (!dnew)
489 greg 2.1 return(NULL);
490 greg 2.5 rmx_addinfo(dnew, rm->info);
491 greg 2.6 dnew->dtype = rm->dtype;
492 greg 2.1 memcpy(dnew->mtx, rm->mtx,
493     sizeof(rm->mtx[0])*rm->ncomp*rm->nrows*rm->ncols);
494     return(dnew);
495     }
496    
497 greg 2.2 /* Allocate and assign transposed matrix */
498     RMATRIX *
499     rmx_transpose(const RMATRIX *rm)
500 greg 2.1 {
501 greg 2.2 RMATRIX *dnew;
502 greg 2.1 int i, j, k;
503    
504 greg 2.33 if (!rm)
505 greg 2.1 return(0);
506 greg 2.31 if ((rm->nrows == 1) | (rm->ncols == 1)) {
507     dnew = rmx_copy(rm);
508 greg 2.33 if (!dnew)
509 greg 2.32 return(NULL);
510 greg 2.31 dnew->nrows = rm->ncols;
511     dnew->ncols = rm->nrows;
512     return(dnew);
513     }
514 greg 2.2 dnew = rmx_alloc(rm->ncols, rm->nrows, rm->ncomp);
515 greg 2.33 if (!dnew)
516 greg 2.2 return(NULL);
517 greg 2.5 if (rm->info) {
518     rmx_addinfo(dnew, rm->info);
519     rmx_addinfo(dnew, "Transposed rows and columns\n");
520     }
521 greg 2.6 dnew->dtype = rm->dtype;
522 greg 2.2 for (i = dnew->nrows; i--; )
523     for (j = dnew->ncols; j--; )
524     for (k = dnew->ncomp; k--; )
525     rmx_lval(dnew,i,j,k) = rmx_lval(rm,j,i,k);
526     return(dnew);
527 greg 2.1 }
528    
529     /* Multiply (concatenate) two matrices and allocate the result */
530     RMATRIX *
531     rmx_multiply(const RMATRIX *m1, const RMATRIX *m2)
532     {
533     RMATRIX *mres;
534     int i, j, k, h;
535    
536 greg 2.33 if (!m1 | !m2 || (m1->ncomp != m2->ncomp) | (m1->ncols != m2->nrows))
537 greg 2.1 return(NULL);
538     mres = rmx_alloc(m1->nrows, m2->ncols, m1->ncomp);
539 greg 2.33 if (!mres)
540 greg 2.1 return(NULL);
541 greg 2.6 i = rmx_newtype(m1->dtype, m2->dtype);
542     if (i)
543     mres->dtype = i;
544     else
545     rmx_addinfo(mres, rmx_mismatch_warn);
546 greg 2.1 for (i = mres->nrows; i--; )
547     for (j = mres->ncols; j--; )
548 greg 2.8 for (k = mres->ncomp; k--; ) {
549 greg 2.1 long double d = 0;
550 greg 2.8 for (h = m1->ncols; h--; )
551 greg 2.24 d += rmx_lval(m1,i,h,k) * rmx_lval(m2,h,j,k);
552 greg 2.1 rmx_lval(mres,i,j,k) = (double)d;
553     }
554     return(mres);
555     }
556    
557 greg 2.25 /* Element-wise multiplication (or division) of m2 into m1 */
558     int
559     rmx_elemult(RMATRIX *m1, const RMATRIX *m2, int divide)
560     {
561     int zeroDivides = 0;
562     int i, j, k;
563    
564 greg 2.33 if (!m1 | !m2 || (m1->ncols != m2->ncols) | (m1->nrows != m2->nrows))
565 greg 2.25 return(0);
566     if ((m2->ncomp > 1) & (m2->ncomp != m1->ncomp))
567     return(0);
568     i = rmx_newtype(m1->dtype, m2->dtype);
569     if (i)
570     m1->dtype = i;
571     else
572     rmx_addinfo(m1, rmx_mismatch_warn);
573     for (i = m1->nrows; i--; )
574     for (j = m1->ncols; j--; )
575     if (divide) {
576     double d;
577     if (m2->ncomp == 1) {
578     d = rmx_lval(m2,i,j,0);
579     if (d == 0) {
580     ++zeroDivides;
581     for (k = m1->ncomp; k--; )
582     rmx_lval(m1,i,j,k) = 0;
583     } else {
584     d = 1./d;
585     for (k = m1->ncomp; k--; )
586     rmx_lval(m1,i,j,k) *= d;
587     }
588     } else
589     for (k = m1->ncomp; k--; ) {
590     d = rmx_lval(m2,i,j,k);
591     if (d == 0) {
592     ++zeroDivides;
593     rmx_lval(m1,i,j,k) = 0;
594     } else
595     rmx_lval(m1,i,j,k) /= d;
596     }
597     } else {
598     if (m2->ncomp == 1) {
599     const double d = rmx_lval(m2,i,j,0);
600     for (k = m1->ncomp; k--; )
601     rmx_lval(m1,i,j,k) *= d;
602     } else
603     for (k = m1->ncomp; k--; )
604     rmx_lval(m1,i,j,k) *= rmx_lval(m2,i,j,k);
605     }
606     if (zeroDivides) {
607 greg 2.34 rmx_addinfo(m1, "WARNING: zero divide(s) corrupted results\n");
608 greg 2.25 errno = ERANGE;
609     }
610     return(1);
611     }
612    
613 greg 2.1 /* Sum second matrix into first, applying scale factor beforehand */
614     int
615     rmx_sum(RMATRIX *msum, const RMATRIX *madd, const double sf[])
616     {
617     double *mysf = NULL;
618     int i, j, k;
619    
620 greg 2.33 if (!msum | !madd ||
621 greg 2.1 (msum->nrows != madd->nrows) |
622     (msum->ncols != madd->ncols) |
623     (msum->ncomp != madd->ncomp))
624     return(0);
625 greg 2.33 if (!sf) {
626 greg 2.1 mysf = (double *)malloc(sizeof(double)*msum->ncomp);
627 greg 2.33 if (!mysf)
628 greg 2.1 return(0);
629     for (k = msum->ncomp; k--; )
630     mysf[k] = 1;
631     sf = mysf;
632     }
633 greg 2.6 i = rmx_newtype(msum->dtype, madd->dtype);
634     if (i)
635     msum->dtype = i;
636     else
637     rmx_addinfo(msum, rmx_mismatch_warn);
638 greg 2.1 for (i = msum->nrows; i--; )
639     for (j = msum->ncols; j--; )
640     for (k = msum->ncomp; k--; )
641     rmx_lval(msum,i,j,k) += sf[k] * rmx_lval(madd,i,j,k);
642 greg 2.33 if (mysf)
643     free(mysf);
644 greg 2.1 return(1);
645     }
646    
647     /* Scale the given matrix by the indicated scalar component vector */
648     int
649     rmx_scale(RMATRIX *rm, const double sf[])
650     {
651     int i, j, k;
652    
653 greg 2.33 if (!rm | !sf)
654 greg 2.1 return(0);
655     for (i = rm->nrows; i--; )
656     for (j = rm->ncols; j--; )
657     for (k = rm->ncomp; k--; )
658     rmx_lval(rm,i,j,k) *= sf[k];
659    
660 greg 2.28 if (rm->info)
661     rmx_addinfo(rm, "Applied scalar\n");
662 greg 2.1 return(1);
663     }
664    
665     /* Allocate new matrix and apply component transformation */
666     RMATRIX *
667     rmx_transform(const RMATRIX *msrc, int n, const double cmat[])
668     {
669     int i, j, ks, kd;
670     RMATRIX *dnew;
671    
672 greg 2.33 if (!msrc | (n <= 0) | !cmat)
673 greg 2.1 return(NULL);
674     dnew = rmx_alloc(msrc->nrows, msrc->ncols, n);
675 greg 2.33 if (!dnew)
676 greg 2.1 return(NULL);
677 greg 2.28 if (msrc->info) {
678     char buf[128];
679 greg 2.38 sprintf(buf, "Applied %dx%d component transform\n",
680 greg 2.28 dnew->ncomp, msrc->ncomp);
681     rmx_addinfo(dnew, msrc->info);
682     rmx_addinfo(dnew, buf);
683     }
684 greg 2.6 dnew->dtype = msrc->dtype;
685 greg 2.1 for (i = dnew->nrows; i--; )
686     for (j = dnew->ncols; j--; )
687     for (kd = dnew->ncomp; kd--; ) {
688     double d = 0;
689     for (ks = msrc->ncomp; ks--; )
690     d += cmat[kd*msrc->ncomp + ks] * rmx_lval(msrc,i,j,ks);
691     rmx_lval(dnew,i,j,kd) = d;
692     }
693     return(dnew);
694     }
695    
696     /* Convert a color matrix to newly allocated RMATRIX buffer */
697     RMATRIX *
698     rmx_from_cmatrix(const CMATRIX *cm)
699     {
700     int i, j;
701     RMATRIX *dnew;
702    
703 greg 2.33 if (!cm)
704 greg 2.1 return(NULL);
705     dnew = rmx_alloc(cm->nrows, cm->ncols, 3);
706 greg 2.33 if (!dnew)
707 greg 2.1 return(NULL);
708 greg 2.6 dnew->dtype = DTfloat;
709 greg 2.1 for (i = dnew->nrows; i--; )
710     for (j = dnew->ncols; j--; ) {
711     const COLORV *cv = cm_lval(cm,i,j);
712     rmx_lval(dnew,i,j,0) = cv[0];
713     rmx_lval(dnew,i,j,1) = cv[1];
714     rmx_lval(dnew,i,j,2) = cv[2];
715     }
716     return(dnew);
717     }
718    
719     /* Convert general matrix to newly allocated CMATRIX buffer */
720     CMATRIX *
721     cm_from_rmatrix(const RMATRIX *rm)
722     {
723     int i, j;
724     CMATRIX *cnew;
725    
726 greg 2.33 if (!rm || rm->ncomp != 3)
727 greg 2.1 return(NULL);
728     cnew = cm_alloc(rm->nrows, rm->ncols);
729 greg 2.33 if (!cnew)
730 greg 2.1 return(NULL);
731     for (i = cnew->nrows; i--; )
732     for (j = cnew->ncols; j--; ) {
733     COLORV *cv = cm_lval(cnew,i,j);
734     cv[0] = (COLORV)rmx_lval(rm,i,j,0);
735     cv[1] = (COLORV)rmx_lval(rm,i,j,1);
736     cv[2] = (COLORV)rmx_lval(rm,i,j,2);
737     }
738     return(cnew);
739     }