ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/dctimestep.c
Revision: 2.7
Committed: Sun Jun 21 21:42:12 2009 UTC (14 years, 9 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.6: +17 -14 lines
Log Message:
Improved BSDF conversion efficiency

File Contents

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