ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/dctimestep.c
Revision: 2.1
Committed: Wed Jun 17 20:41:47 2009 UTC (14 years, 9 months ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
Created dctimestep utility for daylight coefficient calculations

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id$";
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 matrix 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 return(NULL);
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 (fname == NULL)
112 fname = "<stdin>";
113 else if ((fp = fopen(fname, "r")) == NULL) {
114 sprintf(errmsg, "cannot open file '%s'", fname);
115 error(SYSTEM, errmsg);
116 }
117 if (dtype != DTascii)
118 SET_FILE_BINARY(fp);
119 if (dtype == DTfromHeader)
120 dtype = getDTfromHeader(fp);
121 switch (dtype) {
122 case DTascii:
123 case DTfloat:
124 case DTdouble:
125 break;
126 default:
127 error(USER, "unexpected data type in cm_load()");
128 return(NULL);
129 }
130 if (nrows <= 0) { /* don't know length? */
131 int guessrows = 147; /* usually big enough */
132 if ((dtype != DTascii) & (fp != stdin)) {
133 long startpos = ftell(fp);
134 if (fseek(fp, 0L, SEEK_END) == 0) {
135 long endpos = ftell(fp);
136 long elemsiz = 3*(dtype==DTfloat ?
137 sizeof(float) : sizeof(double));
138 guessrows = (endpos - startpos)/elemsiz;
139 if (fseek(fp, startpos, SEEK_SET) < 0) {
140 sprintf(errmsg,
141 "fseek() error on file '%s'",
142 fname);
143 error(SYSTEM, errmsg);
144 }
145 if ((endpos - startpos) % elemsiz == 0)
146 nrows = guessrows; /* confident */
147 }
148 }
149 cm = cm_alloc(guessrows, ncols);
150 } else
151 cm = cm_alloc(nrows, ncols);
152 if (cm == NULL)
153 return(NULL);
154 if (dtype == DTascii) { /* read text file */
155 int maxrow = (nrows > 0 ? nrows : 32000);
156 int r, c;
157 for (r = 0; r < maxrow; r++) {
158 if (r >= cm->nrows) /* need more space? */
159 cm = cm_resize(cm, 2*cm->nrows);
160 for (c = 0; c < ncols; c++) {
161 COLORV *cv = cm_lval(cm,r,c);
162 if (fscanf(fp, COLSPEC, cv, cv+1, cv+2) != 3)
163 if ((nrows <= 0) & (r > 0) & (c == 0)) {
164 cm = cm_resize(cm, maxrow=r);
165 break;
166 } else
167 goto EOFerror;
168 }
169 }
170 while ((c = getc(fp)) != EOF)
171 if (!isspace(c)) {
172 sprintf(errmsg,
173 "unexpected data at end of ascii file %s",
174 fname);
175 error(WARNING, errmsg);
176 break;
177 }
178 } else { /* read binary file */
179 if (sizeof(COLORV) == (dtype==DTfloat ? sizeof(float) :
180 sizeof(double))) {
181 int nread = 0;
182 do { /* read all we can */
183 nread += fread(cm->cmem + 3*nread,
184 3*sizeof(COLORV),
185 cm->nrows*cm->ncols - nread,
186 fp);
187 if (nrows <= 0) { /* unknown length */
188 if (nread == cm->nrows*cm->ncols)
189 /* need more space? */
190 cm = cm_resize(cm, 2*cm->nrows);
191 else if (nread % cm->ncols == 0)
192 /* seem to be done */
193 cm = cm_resize(cm, nread/cm->ncols);
194 else /* ended mid-row */
195 goto EOFerror;
196 } else if (nread < cm->nrows*cm->ncols)
197 goto EOFerror;
198 } while (nread < cm->nrows*cm->ncols);
199
200 } else if (dtype == DTdouble) {
201 double dc[3]; /* load from double */
202 COLORV *cvp = cm->cmem;
203 int n = nrows*ncols;
204
205 if (n <= 0)
206 goto not_handled;
207 while (n--) {
208 if (fread(dc, sizeof(double), 3, fp) != 3)
209 goto EOFerror;
210 copycolor(cvp, dc);
211 cvp += 3;
212 }
213
214 } else /* dtype == DTfloat */ {
215 float fc[3]; /* load from float */
216 COLORV *cvp = cm->cmem;
217 int n = nrows*ncols;
218
219 if (n <= 0)
220 goto not_handled;
221 while (n--) {
222 if (fread(fc, sizeof(float), 3, fp) != 3)
223 goto EOFerror;
224 copycolor(cvp, fc);
225 cvp += 3;
226 }
227 }
228 if (getc(fp) != EOF) {
229 sprintf(errmsg,
230 "unexpected data at end of binary file %s",
231 fname);
232 error(WARNING, errmsg);
233 }
234 }
235 if (fp != stdin)
236 fclose(fp);
237 return(cm);
238 EOFerror:
239 sprintf(errmsg, "unexpected EOF reading %s",
240 fname);
241 error(USER, errmsg);
242 not_handled:
243 error(INTERNAL, "unhandled data size or length in cm_load()");
244 return(NULL); /* gratis return */
245 }
246
247 /* Multiply two matrices (or a matrix and a vector) and allocate the result*/
248 static CMATRIX *
249 cm_multiply(const CMATRIX *cm1, const CMATRIX *cm2)
250 {
251 CMATRIX *cmr;
252 int dr, dc, i;
253
254 if ((cm1->ncols <= 0) | (cm1->ncols != cm2->nrows))
255 error(INTERNAL, "matrix dimension mismatch in cm_multiply()");
256 cmr = cm_alloc(cm1->nrows, cm2->ncols);
257 if (cmr == NULL)
258 return(NULL);
259 for (dr = 0; dr < cmr->nrows; dr++)
260 for (dc = 0; dc < cmr->ncols; dc++) {
261 COLORV *dp = cm_lval(cmr,dr,dc);
262 dp[0] = dp[1] = dp[2] = 0;
263 for (i = 0; i < cm1->ncols; i++) {
264 const COLORV *cp1 = cm_lval(cm1,dr,i);
265 const COLORV *cp2 = cm_lval(cm2,i,dc);
266 dp[0] += cp1[0] * cp2[0];
267 dp[1] += cp1[1] * cp2[1];
268 dp[2] += cp1[2] * cp2[2];
269 }
270 }
271 return(cmr);
272 }
273
274 /* print out matrix as ASCII text -- no header */
275 static void
276 cm_print(const CMATRIX *cm, FILE *fp)
277 {
278 int r, c;
279 const COLORV *mp = cm->cmem;
280
281 for (r = 0; r < cm->nrows; r++) {
282 for (c = 0; c < cm->ncols; c++, mp += 3)
283 fprintf(fp, "\t%.6e %.6e %.6e", mp[0], mp[1], mp[2]);
284 fputc('\n', fp);
285 }
286 }
287
288 /* convert a BSDF to our matrix representation */
289 static CMATRIX *
290 cm_bsdf(const struct BSDF_data *bsdf)
291 {
292 CMATRIX *cm = cm_alloc(bsdf->nout, bsdf->ninc);
293 COLORV *mp = cm->cmem;
294 int r, c;
295
296 for (r = 0; r < cm->nrows; r++)
297 for (c = 0; c < cm->ncols; c++, mp += 3)
298 mp[0] = mp[1] = mp[2] = BSDF_value(bsdf,c,r);
299 return(cm);
300 }
301
302 /* Sum together a set of images and write result to stdout */
303 static int
304 sum_images(const char *fspec, const CMATRIX *cv)
305 {
306 int myDT = DTfromHeader;
307 CMATRIX *pmat;
308 COLOR *scanline;
309 int myXR, myYR;
310 int i, y;
311
312 if (cv->ncols != 1)
313 error(INTERNAL, "expected vector in sum_images()");
314 for (i = 0; i < cv->nrows; i++) {
315 const COLORV *scv = cv_lval(cv,i);
316 char fname[1024];
317 FILE *fp;
318 int dt, xr, yr;
319 COLORV *psp;
320 /* open next picture */
321 sprintf(fname, fspec, i);
322 if ((fp = fopen(fname, "r")) == NULL) {
323 sprintf(errmsg, "cannot open picture '%s'", fname);
324 error(SYSTEM, errmsg);
325 }
326 SET_FILE_BINARY(fp);
327 dt = getDTfromHeader(fp);
328 if ((dt != DTrgbe) & (dt != DTxyze) ||
329 !fscnresolu(&xr, &yr, fp)) {
330 sprintf(errmsg, "file '%s' not a picture", fname);
331 error(USER, errmsg);
332 }
333 if (myDT == DTfromHeader) { /* on first one */
334 myDT = dt;
335 myXR = xr; myYR = yr;
336 scanline = (COLOR *)malloc(sizeof(COLOR)*myXR);
337 if (scanline == NULL)
338 error(SYSTEM, "out of memory in sum_images()");
339 pmat = cm_alloc(myYR, myXR);
340 memset(pmat->cmem, 0, sizeof(COLOR)*myXR*myYR);
341 } else if ((dt != myDT) | (xr != myXR) | (yr != myYR)) {
342 sprintf(errmsg, "picture '%s' format/size mismatch",
343 fname);
344 error(USER, errmsg);
345 }
346 psp = pmat->cmem;
347 for (y = 0; y < yr; y++) { /* read it in */
348 int x;
349 if (freadscan(scanline, xr, fp) < 0) {
350 sprintf(errmsg, "error reading picture '%s'",
351 fname);
352 error(SYSTEM, errmsg);
353 }
354 /* sum in scanline */
355 for (x = 0; x < xr; x++, psp += 3) {
356 multcolor(scanline[x], scv);
357 addcolor(psp, scanline[x]);
358 }
359 }
360 fclose(fp); /* done this picture */
361 }
362 free(scanline);
363 /* finish header */
364 fputformat(myDT==DTrgbe ? COLRFMT : CIEFMT, stdout);
365 fputc('\n', stdout);
366 fprtresolu(myXR, myYR, stdout);
367 /* write scanlines */
368 for (y = 0; y < myYR; y++)
369 if (fwritescan((COLOR *)cm_lval(pmat, y, 0), myXR, stdout) < 0)
370 return(0);
371 cm_free(pmat); /* all done */
372 return(fflush(stdout) == 0);
373 }
374
375 /* check to see if a string contains a %d specification */
376 int
377 hasDecimalSpec(const char *s)
378 {
379 while (*s && *s != '%')
380 s++;
381 if (!*s)
382 return(0);
383 do
384 ++s;
385 while (isdigit(*s));
386
387 return(*s == 'd');
388 }
389
390 int
391 main(int argc, char *argv[])
392 {
393 CMATRIX *tvec, *Dmat, *Tmat, *vtmp, *cvec;
394 struct BSDF_data *btdf;
395
396 progname = argv[0];
397
398 if ((argc < 4) | (argc > 5)) {
399 fprintf(stderr, "Usage: %s Vspec Tbsdf.xml Dmat.dat [tregvec]\n",
400 progname);
401 return(1);
402 }
403 tvec = cm_load(argv[4], 0, 1, DTascii); /* argv[4]==NULL iff argc==4 */
404 Dmat = cm_load(argv[3], 0, tvec->nrows, DTfromHeader);
405 btdf = load_BSDF(argv[2]);
406 if (btdf == NULL)
407 return(1);
408 if (btdf->ninc != Dmat->nrows) {
409 sprintf(errmsg, "Incoming BTDF dir (%d) mismatch to D (%d)",
410 btdf->ninc, Dmat->nrows);
411 error(USER, errmsg);
412 }
413 /* multiply vector through */
414 vtmp = cm_multiply(Dmat, tvec);
415 cm_free(Dmat); cm_free(tvec);
416 Tmat = cm_bsdf(btdf); /* convert BTDF to matrix */
417 free_BSDF(btdf);
418 cvec = cm_multiply(Tmat, vtmp); /* cvec = component vector */
419 cm_free(Tmat); cm_free(vtmp);
420 if (hasDecimalSpec(argv[1])) { /* generating image */
421 SET_FILE_BINARY(stdout);
422 newheader("RADIANCE", stdout);
423 printargs(argc, argv, stdout);
424 fputnow(stdout);
425 if (!sum_images(argv[1], cvec))
426 return(1);
427 } else { /* generating vector */
428 CMATRIX *Vmat = cm_load(argv[1], 0, cvec->nrows, DTfromHeader);
429 CMATRIX *rvec = cm_multiply(Vmat, cvec);
430 cm_free(Vmat);
431 cm_print(rvec, stdout);
432 cm_free(rvec);
433 }
434 cm_free(cvec); /* final clean-up */
435 return(0);
436 }