ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/dctimestep.c
Revision: 2.3
Committed: Fri Jun 19 14:21:42 2009 UTC (14 years, 9 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.2: +26 -3 lines
Log Message:
Fixed conversion of BTDF to matrix

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.3 static const char RCSid[] = "$Id: dctimestep.c,v 2.2 2009/06/19 06:49:42 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     if (r >= cm->nrows) /* need more space? */
164     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     if ((nrows <= 0) & (r > 0) & (c == 0)) {
169     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     else if (nread % cm->ncols == 0)
197     /* 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     if (nbadohm | nneg) {
322     sprintf(errmsg,
323     "BTDF has %d negatives and %d bad incoming solid angles");
324     error(WARNING, errmsg);
325     }
326 greg 2.1 return(cm);
327     }
328    
329     /* Sum together a set of images and write result to stdout */
330     static int
331     sum_images(const char *fspec, const CMATRIX *cv)
332     {
333     int myDT = DTfromHeader;
334     CMATRIX *pmat;
335     COLOR *scanline;
336     int myXR, myYR;
337     int i, y;
338    
339     if (cv->ncols != 1)
340     error(INTERNAL, "expected vector in sum_images()");
341     for (i = 0; i < cv->nrows; i++) {
342     const COLORV *scv = cv_lval(cv,i);
343     char fname[1024];
344     FILE *fp;
345     int dt, xr, yr;
346     COLORV *psp;
347     /* open next picture */
348     sprintf(fname, fspec, i);
349     if ((fp = fopen(fname, "r")) == NULL) {
350     sprintf(errmsg, "cannot open picture '%s'", fname);
351     error(SYSTEM, errmsg);
352     }
353     SET_FILE_BINARY(fp);
354     dt = getDTfromHeader(fp);
355     if ((dt != DTrgbe) & (dt != DTxyze) ||
356     !fscnresolu(&xr, &yr, fp)) {
357     sprintf(errmsg, "file '%s' not a picture", fname);
358     error(USER, errmsg);
359     }
360     if (myDT == DTfromHeader) { /* on first one */
361     myDT = dt;
362     myXR = xr; myYR = yr;
363     scanline = (COLOR *)malloc(sizeof(COLOR)*myXR);
364     if (scanline == NULL)
365     error(SYSTEM, "out of memory in sum_images()");
366     pmat = cm_alloc(myYR, myXR);
367     memset(pmat->cmem, 0, sizeof(COLOR)*myXR*myYR);
368 greg 2.2 /* finish header */
369     fputformat(myDT==DTrgbe ? COLRFMT : CIEFMT, stdout);
370     fputc('\n', stdout);
371     fprtresolu(myXR, myYR, stdout);
372     fflush(stdout);
373 greg 2.1 } else if ((dt != myDT) | (xr != myXR) | (yr != myYR)) {
374     sprintf(errmsg, "picture '%s' format/size mismatch",
375     fname);
376     error(USER, errmsg);
377     }
378     psp = pmat->cmem;
379     for (y = 0; y < yr; y++) { /* read it in */
380     int x;
381     if (freadscan(scanline, xr, fp) < 0) {
382     sprintf(errmsg, "error reading picture '%s'",
383     fname);
384     error(SYSTEM, errmsg);
385     }
386     /* sum in scanline */
387     for (x = 0; x < xr; x++, psp += 3) {
388     multcolor(scanline[x], scv);
389     addcolor(psp, scanline[x]);
390     }
391     }
392     fclose(fp); /* done this picture */
393     }
394     free(scanline);
395     /* write scanlines */
396     for (y = 0; y < myYR; y++)
397     if (fwritescan((COLOR *)cm_lval(pmat, y, 0), myXR, stdout) < 0)
398     return(0);
399     cm_free(pmat); /* all done */
400     return(fflush(stdout) == 0);
401     }
402    
403     /* check to see if a string contains a %d specification */
404     int
405     hasDecimalSpec(const char *s)
406     {
407     while (*s && *s != '%')
408     s++;
409     if (!*s)
410     return(0);
411     do
412     ++s;
413     while (isdigit(*s));
414    
415     return(*s == 'd');
416     }
417    
418     int
419     main(int argc, char *argv[])
420     {
421 greg 2.2 CMATRIX *tvec, *Dmat, *Tmat, *ivec, *cvec;
422 greg 2.1 struct BSDF_data *btdf;
423    
424     progname = argv[0];
425    
426     if ((argc < 4) | (argc > 5)) {
427     fprintf(stderr, "Usage: %s Vspec Tbsdf.xml Dmat.dat [tregvec]\n",
428     progname);
429     return(1);
430     }
431     tvec = cm_load(argv[4], 0, 1, DTascii); /* argv[4]==NULL iff argc==4 */
432     Dmat = cm_load(argv[3], 0, tvec->nrows, DTfromHeader);
433     btdf = load_BSDF(argv[2]);
434     if (btdf == NULL)
435     return(1);
436     if (btdf->ninc != Dmat->nrows) {
437     sprintf(errmsg, "Incoming BTDF dir (%d) mismatch to D (%d)",
438     btdf->ninc, Dmat->nrows);
439     error(USER, errmsg);
440     }
441     /* multiply vector through */
442 greg 2.2 ivec = cm_multiply(Dmat, tvec);
443 greg 2.1 cm_free(Dmat); cm_free(tvec);
444     Tmat = cm_bsdf(btdf); /* convert BTDF to matrix */
445     free_BSDF(btdf);
446 greg 2.2 cvec = cm_multiply(Tmat, ivec); /* cvec = component vector */
447     cm_free(Tmat); cm_free(ivec);
448 greg 2.1 if (hasDecimalSpec(argv[1])) { /* generating image */
449     SET_FILE_BINARY(stdout);
450     newheader("RADIANCE", stdout);
451     printargs(argc, argv, stdout);
452     fputnow(stdout);
453     if (!sum_images(argv[1], cvec))
454     return(1);
455     } else { /* generating vector */
456     CMATRIX *Vmat = cm_load(argv[1], 0, cvec->nrows, DTfromHeader);
457     CMATRIX *rvec = cm_multiply(Vmat, cvec);
458     cm_free(Vmat);
459     cm_print(rvec, stdout);
460     cm_free(rvec);
461     }
462     cm_free(cvec); /* final clean-up */
463     return(0);
464     }