ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmatrix.c
Revision: 2.25
Committed: Mon Aug 28 15:59:46 2017 UTC (6 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.24: +59 -1 lines
Log Message:
Added element-wise multiplication and division to rmtxop command

File Contents

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