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 (13 years, 5 months ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
Created dctimestep utility for daylight coefficient calculations

File Contents

# User Rev Content
1 greg 2.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     }