ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmatrix.c
Revision: 2.5
Committed: Fri Aug 1 23:37:24 2014 UTC (9 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.4: +47 -3 lines
Log Message:
Maintain header information in matrix structure

File Contents

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