ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/dctimestep.c
Revision: 2.6
Committed: Sat Jun 20 04:37:43 2009 UTC (14 years, 9 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.5: +4 -4 lines
Log Message:
Added check against empty input

File Contents

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