ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmatrix.c
Revision: 2.1
Committed: Sat May 31 05:02:37 2014 UTC (9 years, 10 months ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
Created rmtxop command, need to test and debug

File Contents

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