ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rmatrix.c
Revision: 2.4
Committed: Thu Jul 24 16:28:17 2014 UTC (9 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R2
Changes since 2.3: +9 -3 lines
Log Message:
Made NROWS/NCOLS more consistent with resolution string

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: rmatrix.c,v 2.3 2014/07/08 16:39:41 greg Exp $";
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.nrows <= 0) | (dinfo.ncols <= 0)) {
181 if (!fscnresolu(&dinfo.ncols, &dinfo.nrows, fp)) {
182 fclose(fp);
183 return(NULL);
184 }
185 if (dinfo.ncomp <= 0)
186 dinfo.ncomp = 3;
187 else if ((dinfo.dtype == DTrgbe) | (dinfo.dtype == DTxyze) &&
188 dinfo.ncomp != 3) {
189 fclose(fp);
190 return(NULL);
191 }
192 }
193 dnew = rmx_alloc(dinfo.nrows, dinfo.ncols, dinfo.ncomp);
194 if (dnew == NULL) {
195 fclose(fp);
196 return(NULL);
197 }
198 switch (dinfo.dtype) {
199 case DTascii:
200 if (!rmx_load_ascii(dnew, fp))
201 goto loaderr;
202 break;
203 case DTfloat:
204 if (!rmx_load_float(dnew, fp))
205 goto loaderr;
206 break;
207 case DTdouble:
208 if (!rmx_load_double(dnew, fp))
209 goto loaderr;
210 break;
211 case DTrgbe:
212 case DTxyze:
213 if (!rmx_load_rgbe(dnew, fp))
214 goto loaderr;
215 break;
216 default:
217 goto loaderr;
218 }
219 if (fp != stdin)
220 fclose(fp);
221 return(dnew);
222 loaderr: /* should report error? */
223 fclose(fp);
224 rmx_free(dnew);
225 return(NULL);
226 }
227
228 static int
229 rmx_write_ascii(const RMATRIX *rm, FILE *fp)
230 {
231 int i, j, k;
232 #ifdef _WIN32
233 _setmode(fileno(fp), _O_TEXT);
234 #endif
235 for (i = 0; i < rm->nrows; i++) {
236 for (j = 0; j < rm->ncols; j++) {
237 for (k = 0; k < rm->ncomp; k++)
238 fprintf(fp, " %.15e", rmx_lval(rm,i,j,k));
239 fputc('\t', fp);
240 }
241 fputc('\n', fp);
242 }
243 return(1);
244 }
245
246 static int
247 rmx_write_float(const RMATRIX *rm, FILE *fp)
248 {
249 int i, j, k;
250 float val[100];
251
252 if (rm->ncomp > 100) {
253 fputs("Unsupported # components in rmx_write_float()\n", stderr);
254 exit(1);
255 }
256 for (i = 0; i < rm->nrows; i++)
257 for (j = 0; j < rm->ncols; j++) {
258 for (k = rm->ncomp; k--; )
259 val[k] = (float)rmx_lval(rm,i,j,k);
260 if (fwrite(val, sizeof(val[0]), rm->ncomp, fp) != rm->ncomp)
261 return(0);
262 }
263 return(1);
264 }
265
266 static int
267 rmx_write_double(const RMATRIX *rm, FILE *fp)
268 {
269 int i, j, k;
270 double val[100];
271
272 if (rm->ncomp > 100) {
273 fputs("Unsupported # components in rmx_write_double()\n", stderr);
274 exit(1);
275 }
276 for (i = 0; i < rm->nrows; i++)
277 for (j = 0; j < rm->ncols; j++) {
278 for (k = rm->ncomp; k--; )
279 val[k] = rmx_lval(rm,i,j,k);
280 if (fwrite(val, sizeof(val[0]), rm->ncomp, fp) != rm->ncomp)
281 return(0);
282 }
283 return(1);
284 }
285
286 static int
287 rmx_write_rgbe(const RMATRIX *rm, FILE *fp)
288 {
289 COLOR *scan = (COLOR *)malloc(sizeof(COLOR)*rm->ncols);
290 int i, j;
291
292 if (scan == NULL)
293 return(0);
294 for (i = 0; i < rm->nrows; i++) {
295 for (j = rm->ncols; j--; )
296 setcolor(scan[j], rmx_lval(rm,i,j,0),
297 rmx_lval(rm,i,j,1),
298 rmx_lval(rm,i,j,2) );
299 if (fwritescan(scan, rm->ncols, fp) < 0) {
300 free(scan);
301 return(0);
302 }
303 }
304 free(scan);
305 return(1);
306 }
307
308 /* Write matrix to file type indicated by dt */
309 long
310 rmx_write(const RMATRIX *rm, int dtype, FILE *fp)
311 {
312 RMATRIX *mydm = NULL;
313 int ok = 1;
314
315 if ((rm == NULL) | (fp == NULL))
316 return(0);
317 /* complete header */
318 if ((dtype != DTrgbe) & (dtype != DTxyze)) {
319 fprintf(fp, "NROWS=%d\n", rm->nrows);
320 fprintf(fp, "NCOLS=%d\n", rm->ncols);
321 fprintf(fp, "NCOMP=%d\n", rm->ncomp);
322 } else if (rm->ncomp != 3) { /* wrong # components? */
323 double cmtx[3];
324 if (rm->ncomp != 1) /* only convert grayscale */
325 return(0);
326 cmtx[0] = cmtx[1] = cmtx[2] = 1;
327 mydm = rmx_transform(rm, 3, cmtx);
328 if (mydm == NULL)
329 return(0);
330 rm = mydm;
331 }
332 fputformat((char *)cm_fmt_id[dtype], fp);
333 fputc('\n', fp);
334 switch (dtype) { /* write data */
335 case DTascii:
336 ok = rmx_write_ascii(rm, fp);
337 break;
338 case DTfloat:
339 ok = rmx_write_float(rm, fp);
340 break;
341 case DTdouble:
342 ok = rmx_write_double(rm, fp);
343 break;
344 case DTrgbe:
345 case DTxyze:
346 fprtresolu(rm->ncols, rm->nrows, fp);
347 ok = rmx_write_rgbe(rm, fp);
348 break;
349 default:
350 return(0);
351 }
352 ok &= (fflush(fp) == 0);
353 rmx_free(mydm);
354 return(ftell(fp) * ok); /* return # bytes written */
355 }
356
357 /* Allocate and assign square identity matrix with n components */
358 RMATRIX *
359 rmx_identity(const int dim, const int n)
360 {
361 RMATRIX *rid = rmx_alloc(dim, dim, n);
362 int i;
363
364 if (rid == NULL)
365 return(NULL);
366 memset(rid->mtx, 0, sizeof(rid->mtx[0])*dim*dim);
367 for (i = dim; i--; )
368 rmx_lval(rid,i,i,0) = 1;
369 for (i = n; --i; )
370 memcpy(rid->mtx+i*(dim*dim), rid->mtx,
371 sizeof(rid->mtx[0])*dim*dim);
372 return(rid);
373 }
374
375 /* Duplicate the given matrix */
376 RMATRIX *
377 rmx_copy(const RMATRIX *rm)
378 {
379 RMATRIX *dnew;
380
381 if (rm == NULL)
382 return(NULL);
383 dnew = rmx_alloc(rm->nrows, rm->ncols, rm->ncomp);
384 if (dnew == NULL)
385 return(NULL);
386 memcpy(dnew->mtx, rm->mtx,
387 sizeof(rm->mtx[0])*rm->ncomp*rm->nrows*rm->ncols);
388 return(dnew);
389 }
390
391 /* Allocate and assign transposed matrix */
392 RMATRIX *
393 rmx_transpose(const RMATRIX *rm)
394 {
395 RMATRIX *dnew;
396 int i, j, k;
397
398 if (rm == NULL)
399 return(0);
400 dnew = rmx_alloc(rm->ncols, rm->nrows, rm->ncomp);
401 if (dnew == NULL)
402 return(NULL);
403 for (i = dnew->nrows; i--; )
404 for (j = dnew->ncols; j--; )
405 for (k = dnew->ncomp; k--; )
406 rmx_lval(dnew,i,j,k) = rmx_lval(rm,j,i,k);
407 return(dnew);
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 }