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 (13 years, 8 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

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.13 static const char RCSid[] = "$Id: dctimestep.c,v 2.12 2010/07/01 21:54:56 greg Exp $";
3 greg 2.1 #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 greg 2.2 /* Data types for file loading */
20 greg 2.1 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 greg 2.10 error(USER, "attempt to create empty matrix");
42 greg 2.1 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 greg 2.12 if (ncols <= 0)
112     error(USER, "Non-positive number of columns");
113 greg 2.1 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 greg 2.2
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 greg 2.1 if (fseek(fp, startpos, SEEK_SET) < 0) {
148     sprintf(errmsg,
149     "fseek() error on file '%s'",
150     fname);
151     error(SYSTEM, errmsg);
152     }
153 greg 2.2 nrows = guessrows; /* we're confident */
154 greg 2.1 }
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 greg 2.6 if (r >= cm->nrows) /* need more space? */
166 greg 2.1 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 greg 2.6 if ((nrows <= 0) & (r > 0) & !c) {
171 greg 2.1 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 greg 2.6 else if (nread && !(nread % cm->ncols))
199 greg 2.1 /* 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 greg 2.3 int nbadohm = 0;
300     int nneg = 0;
301 greg 2.1 int r, c;
302    
303 greg 2.8 for (c = 0; c < cm->ncols; c++) {
304 greg 2.7 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 greg 2.3 float f = BSDF_value(bsdf,c,r);
317 greg 2.9 COLORV *mp = cm_lval(cm,r,c);
318 greg 2.7
319 greg 2.3 if (f <= .0) {
320     nneg += (f < -FTINY);
321 greg 2.8 f = .0f;
322 greg 2.3 }
323 greg 2.7 mp[0] = mp[1] = mp[2] = f * dom;
324 greg 2.3 }
325 greg 2.7 }
326 greg 2.4 if (nneg || nbadohm) {
327 greg 2.3 sprintf(errmsg,
328 greg 2.4 "BTDF has %d negatives and %d bad incoming solid angles",
329     nneg, nbadohm);
330 greg 2.3 error(WARNING, errmsg);
331     }
332 greg 2.1 return(cm);
333     }
334    
335     /* Sum together a set of images and write result to stdout */
336     static int
337 greg 2.13 sum_images(const char *fspec, const CMATRIX *cv, FILE *fout)
338 greg 2.1 {
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 greg 2.2 /* finish header */
375 greg 2.13 fputformat(myDT==DTrgbe ? COLRFMT : CIEFMT, fout);
376     fputc('\n', fout);
377     fprtresolu(myXR, myYR, fout);
378     fflush(fout);
379 greg 2.1 } 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 greg 2.13 if (fwritescan((COLOR *)cm_lval(pmat, y, 0), myXR, fout) < 0)
404 greg 2.1 return(0);
405     cm_free(pmat); /* all done */
406 greg 2.13 return(fflush(fout) == 0);
407 greg 2.1 }
408    
409 greg 2.11 /* check to see if a string contains a %d or %o specification */
410     static int
411     hasNumberFormat(const char *s)
412 greg 2.1 {
413     while (*s && *s != '%')
414     s++;
415     if (!*s)
416     return(0);
417     do
418     ++s;
419     while (isdigit(*s));
420    
421 greg 2.11 return(*s == 'd' | *s == 'i' | *s == 'o' | *s == 'x' | *s == 'X');
422 greg 2.1 }
423    
424     int
425     main(int argc, char *argv[])
426     {
427 greg 2.12 CMATRIX *cvec;
428 greg 2.1
429     progname = argv[0];
430    
431 greg 2.12 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 greg 2.1 return(1);
436     }
437 greg 2.12
438 greg 2.13 if (argc > 3) { /* VTDs expression */
439 greg 2.12 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 greg 2.5 /* load Daylight matrix */
448 greg 2.12 Dmat = cm_load(argv[3], btdf->ninc, svec->nrows, DTfromHeader);
449 greg 2.1 /* multiply vector through */
450 greg 2.12 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 greg 2.11 if (hasNumberFormat(argv[1])) { /* generating image */
461 greg 2.1 SET_FILE_BINARY(stdout);
462     newheader("RADIANCE", stdout);
463     printargs(argc, argv, stdout);
464     fputnow(stdout);
465 greg 2.12 if (!sum_images(argv[1], cvec, stdout))
466 greg 2.1 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     }