ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/dctimestep.c
Revision: 2.13
Committed: Thu Jul 1 22:44:25 2010 UTC (12 years, 5 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.12: +9 -9 lines
Log Message:
Bug fix in last change

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: dctimestep.c,v 2.12 2010/07/01 21:54:56 greg Exp $";
3 #endif
4 /*
5 * Compute time-step result using Daylight Coefficient method.
6 *
7 * G. Ward
8 */
9
10 #include <ctype.h>
11 #include "standard.h"
12 #include "platform.h"
13 #include "color.h"
14 #include "resolu.h"
15 #include "bsdf.h"
16
17 char *progname; /* global argv[0] */
18
19 /* Data types for file loading */
20 enum {DTfromHeader, DTascii, DTfloat, DTdouble, DTrgbe, DTxyze};
21
22 /* A color coefficient matrix -- vectors have ncols==1 */
23 typedef struct {
24 int nrows, ncols;
25 COLORV cmem[3]; /* extends struct */
26 } CMATRIX;
27
28 #define COLSPEC (sizeof(COLORV)==sizeof(float) ? "%f %f %f" : "%lf %lf %lf")
29
30 #define cm_lval(cm,r,c) ((cm)->cmem + 3*((r)*(cm)->ncols + (c)))
31
32 #define cv_lval(cm,i) ((cm)->cmem + 3*(i))
33
34 /* Allocate a color coefficient matrix */
35 static CMATRIX *
36 cm_alloc(int nrows, int ncols)
37 {
38 CMATRIX *cm;
39
40 if ((nrows <= 0) | (ncols <= 0))
41 error(USER, "attempt to create empty matrix");
42 cm = (CMATRIX *)malloc(sizeof(CMATRIX) +
43 3*sizeof(COLORV)*(nrows*ncols - 1));
44 if (cm == NULL)
45 error(SYSTEM, "out of memory in cm_alloc()");
46 cm->nrows = nrows;
47 cm->ncols = ncols;
48 return(cm);
49 }
50
51 #define cm_free(cm) free(cm)
52
53 /* Resize color coefficient matrix */
54 static CMATRIX *
55 cm_resize(CMATRIX *cm, int nrows)
56 {
57 if (nrows == cm->nrows)
58 return(cm);
59 if (nrows <= 0) {
60 cm_free(cm);
61 return(NULL);
62 }
63 cm = (CMATRIX *)realloc(cm, sizeof(CMATRIX) +
64 3*sizeof(COLORV)*(nrows*cm->ncols - 1));
65 if (cm == NULL)
66 error(SYSTEM, "out of memory in cm_resize()");
67 cm->nrows = nrows;
68 return(cm);
69 }
70
71 /* Load header to obtain data type */
72 static int
73 getDT(char *s, void *p)
74 {
75 char fmt[32];
76
77 if (formatval(fmt, s)) {
78 if (!strcmp(fmt, "ascii"))
79 *((int *)p) = DTascii;
80 else if (!strcmp(fmt, "float"))
81 *((int *)p) = DTfloat;
82 else if (!strcmp(fmt, "double"))
83 *((int *)p) = DTdouble;
84 else if (!strcmp(fmt, COLRFMT))
85 *((int *)p) = DTrgbe;
86 else if (!strcmp(fmt, CIEFMT))
87 *((int *)p) = DTxyze;
88 }
89 return(0);
90 }
91
92 static int
93 getDTfromHeader(FILE *fp)
94 {
95 int dt = DTfromHeader;
96
97 if (getheader(fp, getDT, &dt) < 0)
98 error(SYSTEM, "header read error");
99 if (dt == DTfromHeader)
100 error(USER, "missing data format in header");
101 return(dt);
102 }
103
104 /* Allocate and load a matrix from the given file (or stdin if NULL) */
105 static CMATRIX *
106 cm_load(const char *fname, int nrows, int ncols, int dtype)
107 {
108 CMATRIX *cm;
109 FILE *fp = stdin;
110
111 if (ncols <= 0)
112 error(USER, "Non-positive number of columns");
113 if (fname == NULL)
114 fname = "<stdin>";
115 else if ((fp = fopen(fname, "r")) == NULL) {
116 sprintf(errmsg, "cannot open file '%s'", fname);
117 error(SYSTEM, errmsg);
118 }
119 if (dtype != DTascii)
120 SET_FILE_BINARY(fp);
121 if (dtype == DTfromHeader)
122 dtype = getDTfromHeader(fp);
123 switch (dtype) {
124 case DTascii:
125 case DTfloat:
126 case DTdouble:
127 break;
128 default:
129 error(USER, "unexpected data type in cm_load()");
130 }
131 if (nrows <= 0) { /* don't know length? */
132 int guessrows = 147; /* usually big enough */
133 if ((dtype != DTascii) & (fp != stdin)) {
134 long startpos = ftell(fp);
135 if (fseek(fp, 0L, SEEK_END) == 0) {
136 long endpos = ftell(fp);
137 long elemsiz = 3*(dtype==DTfloat ?
138 sizeof(float) : sizeof(double));
139
140 if ((endpos - startpos) % (ncols*elemsiz)) {
141 sprintf(errmsg,
142 "improper length for binary file '%s'",
143 fname);
144 error(USER, errmsg);
145 }
146 guessrows = (endpos - startpos)/(ncols*elemsiz);
147 if (fseek(fp, startpos, SEEK_SET) < 0) {
148 sprintf(errmsg,
149 "fseek() error on file '%s'",
150 fname);
151 error(SYSTEM, errmsg);
152 }
153 nrows = guessrows; /* we're confident */
154 }
155 }
156 cm = cm_alloc(guessrows, ncols);
157 } else
158 cm = cm_alloc(nrows, ncols);
159 if (cm == NULL)
160 return(NULL);
161 if (dtype == DTascii) { /* read text file */
162 int maxrow = (nrows > 0 ? nrows : 32000);
163 int r, c;
164 for (r = 0; r < maxrow; r++) {
165 if (r >= cm->nrows) /* need more space? */
166 cm = cm_resize(cm, 2*cm->nrows);
167 for (c = 0; c < ncols; c++) {
168 COLORV *cv = cm_lval(cm,r,c);
169 if (fscanf(fp, COLSPEC, cv, cv+1, cv+2) != 3)
170 if ((nrows <= 0) & (r > 0) & !c) {
171 cm = cm_resize(cm, maxrow=r);
172 break;
173 } else
174 goto EOFerror;
175 }
176 }
177 while ((c = getc(fp)) != EOF)
178 if (!isspace(c)) {
179 sprintf(errmsg,
180 "unexpected data at end of ascii file %s",
181 fname);
182 error(WARNING, errmsg);
183 break;
184 }
185 } else { /* read binary file */
186 if (sizeof(COLORV) == (dtype==DTfloat ? sizeof(float) :
187 sizeof(double))) {
188 int nread = 0;
189 do { /* read all we can */
190 nread += fread(cm->cmem + 3*nread,
191 3*sizeof(COLORV),
192 cm->nrows*cm->ncols - nread,
193 fp);
194 if (nrows <= 0) { /* unknown length */
195 if (nread == cm->nrows*cm->ncols)
196 /* need more space? */
197 cm = cm_resize(cm, 2*cm->nrows);
198 else if (nread && !(nread % cm->ncols))
199 /* seem to be done */
200 cm = cm_resize(cm, nread/cm->ncols);
201 else /* ended mid-row */
202 goto EOFerror;
203 } else if (nread < cm->nrows*cm->ncols)
204 goto EOFerror;
205 } while (nread < cm->nrows*cm->ncols);
206
207 } else if (dtype == DTdouble) {
208 double dc[3]; /* load from double */
209 COLORV *cvp = cm->cmem;
210 int n = nrows*ncols;
211
212 if (n <= 0)
213 goto not_handled;
214 while (n--) {
215 if (fread(dc, sizeof(double), 3, fp) != 3)
216 goto EOFerror;
217 copycolor(cvp, dc);
218 cvp += 3;
219 }
220 } else /* dtype == DTfloat */ {
221 float fc[3]; /* load from float */
222 COLORV *cvp = cm->cmem;
223 int n = nrows*ncols;
224
225 if (n <= 0)
226 goto not_handled;
227 while (n--) {
228 if (fread(fc, sizeof(float), 3, fp) != 3)
229 goto EOFerror;
230 copycolor(cvp, fc);
231 cvp += 3;
232 }
233 }
234 if (getc(fp) != EOF) {
235 sprintf(errmsg,
236 "unexpected data at end of binary file %s",
237 fname);
238 error(WARNING, errmsg);
239 }
240 }
241 if (fp != stdin)
242 fclose(fp);
243 return(cm);
244 EOFerror:
245 sprintf(errmsg, "unexpected EOF reading %s",
246 fname);
247 error(USER, errmsg);
248 not_handled:
249 error(INTERNAL, "unhandled data size or length in cm_load()");
250 return(NULL); /* gratis return */
251 }
252
253 /* Multiply two matrices (or a matrix and a vector) and allocate the result*/
254 static CMATRIX *
255 cm_multiply(const CMATRIX *cm1, const CMATRIX *cm2)
256 {
257 CMATRIX *cmr;
258 int dr, dc, i;
259
260 if ((cm1->ncols <= 0) | (cm1->ncols != cm2->nrows))
261 error(INTERNAL, "matrix dimension mismatch in cm_multiply()");
262 cmr = cm_alloc(cm1->nrows, cm2->ncols);
263 if (cmr == NULL)
264 return(NULL);
265 for (dr = 0; dr < cmr->nrows; dr++)
266 for (dc = 0; dc < cmr->ncols; dc++) {
267 COLORV *dp = cm_lval(cmr,dr,dc);
268 dp[0] = dp[1] = dp[2] = 0;
269 for (i = 0; i < cm1->ncols; i++) {
270 const COLORV *cp1 = cm_lval(cm1,dr,i);
271 const COLORV *cp2 = cm_lval(cm2,i,dc);
272 dp[0] += cp1[0] * cp2[0];
273 dp[1] += cp1[1] * cp2[1];
274 dp[2] += cp1[2] * cp2[2];
275 }
276 }
277 return(cmr);
278 }
279
280 /* print out matrix as ASCII text -- no header */
281 static void
282 cm_print(const CMATRIX *cm, FILE *fp)
283 {
284 int r, c;
285 const COLORV *mp = cm->cmem;
286
287 for (r = 0; r < cm->nrows; r++) {
288 for (c = 0; c < cm->ncols; c++, mp += 3)
289 fprintf(fp, "\t%.6e %.6e %.6e", mp[0], mp[1], mp[2]);
290 fputc('\n', fp);
291 }
292 }
293
294 /* convert a BSDF to our matrix representation */
295 static CMATRIX *
296 cm_bsdf(const struct BSDF_data *bsdf)
297 {
298 CMATRIX *cm = cm_alloc(bsdf->nout, bsdf->ninc);
299 int nbadohm = 0;
300 int nneg = 0;
301 int r, c;
302
303 for (c = 0; c < cm->ncols; c++) {
304 float dom = getBSDF_incohm(bsdf,c);
305 FVECT v;
306
307 if (dom <= .0) {
308 nbadohm++;
309 continue;
310 }
311 if (!getBSDF_incvec(v,bsdf,c) || v[2] > FTINY)
312 error(USER, "illegal incoming BTDF direction");
313 dom *= -v[2];
314
315 for (r = 0; r < cm->nrows; r++) {
316 float f = BSDF_value(bsdf,c,r);
317 COLORV *mp = cm_lval(cm,r,c);
318
319 if (f <= .0) {
320 nneg += (f < -FTINY);
321 f = .0f;
322 }
323 mp[0] = mp[1] = mp[2] = f * dom;
324 }
325 }
326 if (nneg || nbadohm) {
327 sprintf(errmsg,
328 "BTDF has %d negatives and %d bad incoming solid angles",
329 nneg, nbadohm);
330 error(WARNING, errmsg);
331 }
332 return(cm);
333 }
334
335 /* Sum together a set of images and write result to stdout */
336 static int
337 sum_images(const char *fspec, const CMATRIX *cv, FILE *fout)
338 {
339 int myDT = DTfromHeader;
340 CMATRIX *pmat;
341 COLOR *scanline;
342 int myXR, myYR;
343 int i, y;
344
345 if (cv->ncols != 1)
346 error(INTERNAL, "expected vector in sum_images()");
347 for (i = 0; i < cv->nrows; i++) {
348 const COLORV *scv = cv_lval(cv,i);
349 char fname[1024];
350 FILE *fp;
351 int dt, xr, yr;
352 COLORV *psp;
353 /* open next picture */
354 sprintf(fname, fspec, i);
355 if ((fp = fopen(fname, "r")) == NULL) {
356 sprintf(errmsg, "cannot open picture '%s'", fname);
357 error(SYSTEM, errmsg);
358 }
359 SET_FILE_BINARY(fp);
360 dt = getDTfromHeader(fp);
361 if ((dt != DTrgbe) & (dt != DTxyze) ||
362 !fscnresolu(&xr, &yr, fp)) {
363 sprintf(errmsg, "file '%s' not a picture", fname);
364 error(USER, errmsg);
365 }
366 if (myDT == DTfromHeader) { /* on first one */
367 myDT = dt;
368 myXR = xr; myYR = yr;
369 scanline = (COLOR *)malloc(sizeof(COLOR)*myXR);
370 if (scanline == NULL)
371 error(SYSTEM, "out of memory in sum_images()");
372 pmat = cm_alloc(myYR, myXR);
373 memset(pmat->cmem, 0, sizeof(COLOR)*myXR*myYR);
374 /* finish header */
375 fputformat(myDT==DTrgbe ? COLRFMT : CIEFMT, fout);
376 fputc('\n', fout);
377 fprtresolu(myXR, myYR, fout);
378 fflush(fout);
379 } else if ((dt != myDT) | (xr != myXR) | (yr != myYR)) {
380 sprintf(errmsg, "picture '%s' format/size mismatch",
381 fname);
382 error(USER, errmsg);
383 }
384 psp = pmat->cmem;
385 for (y = 0; y < yr; y++) { /* read it in */
386 int x;
387 if (freadscan(scanline, xr, fp) < 0) {
388 sprintf(errmsg, "error reading picture '%s'",
389 fname);
390 error(SYSTEM, errmsg);
391 }
392 /* sum in scanline */
393 for (x = 0; x < xr; x++, psp += 3) {
394 multcolor(scanline[x], scv);
395 addcolor(psp, scanline[x]);
396 }
397 }
398 fclose(fp); /* done this picture */
399 }
400 free(scanline);
401 /* write scanlines */
402 for (y = 0; y < myYR; y++)
403 if (fwritescan((COLOR *)cm_lval(pmat, y, 0), myXR, fout) < 0)
404 return(0);
405 cm_free(pmat); /* all done */
406 return(fflush(fout) == 0);
407 }
408
409 /* check to see if a string contains a %d or %o specification */
410 static int
411 hasNumberFormat(const char *s)
412 {
413 while (*s && *s != '%')
414 s++;
415 if (!*s)
416 return(0);
417 do
418 ++s;
419 while (isdigit(*s));
420
421 return(*s == 'd' | *s == 'i' | *s == 'o' | *s == 'x' | *s == 'X');
422 }
423
424 int
425 main(int argc, char *argv[])
426 {
427 CMATRIX *cvec;
428
429 progname = argv[0];
430
431 if ((argc < 2) | (argc > 5)) {
432 fprintf(stderr, "Usage: %s DCspec [tregvec]\n", progname);
433 fprintf(stderr, " or: %s Vspec Tbsdf.xml Dmat.dat [tregvec]\n",
434 progname);
435 return(1);
436 }
437
438 if (argc > 3) { /* VTDs expression */
439 CMATRIX *svec, *Dmat, *Tmat, *ivec;
440 struct BSDF_data *btdf;
441 /* get sky vector */
442 svec = cm_load(argv[4], 0, 1, DTascii);
443 /* load BSDF */
444 btdf = load_BSDF(argv[2]);
445 if (btdf == NULL)
446 return(1);
447 /* load Daylight matrix */
448 Dmat = cm_load(argv[3], btdf->ninc, svec->nrows, DTfromHeader);
449 /* multiply vector through */
450 ivec = cm_multiply(Dmat, svec);
451 cm_free(Dmat); cm_free(svec);
452 Tmat = cm_bsdf(btdf); /* convert BTDF to matrix */
453 free_BSDF(btdf);
454 cvec = cm_multiply(Tmat, ivec); /* cvec = component vector */
455 cm_free(Tmat); cm_free(ivec);
456 } else { /* else just use sky vector */
457 cvec = cm_load(argv[2], 0, 1, DTascii);
458 }
459
460 if (hasNumberFormat(argv[1])) { /* generating image */
461 SET_FILE_BINARY(stdout);
462 newheader("RADIANCE", stdout);
463 printargs(argc, argv, stdout);
464 fputnow(stdout);
465 if (!sum_images(argv[1], cvec, stdout))
466 return(1);
467 } else { /* generating vector */
468 CMATRIX *Vmat = cm_load(argv[1], 0, cvec->nrows, DTfromHeader);
469 CMATRIX *rvec = cm_multiply(Vmat, cvec);
470 cm_free(Vmat);
471 cm_print(rvec, stdout);
472 cm_free(rvec);
473 }
474 cm_free(cvec); /* final clean-up */
475 return(0);
476 }