ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rttree_reduce.c
Revision: 2.1
Committed: Thu May 26 15:32:02 2011 UTC (12 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
Log Message:
Initial rttree_reduce implementation

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2     static const char RCSid[] = "$Id$";
3     #endif
4     /*
5     * A utility called by genBSDF.pl to reduce tensor tree samples and output
6     * in a standard format as required by XML specification for variable
7     * resolution BSDF data. We are not meant to be run by the user directly.
8     */
9    
10     #include "rtio.h"
11     #include "rterror.h"
12     #include "platform.h"
13     #include <stdlib.h>
14    
15     float *datarr; /* our loaded BSDF data array */
16     int ttrank = 4; /* tensor tree rank */
17     int log2g = 4; /* log2 of grid resolution */
18     int infmt = 'a'; /* input format ('a','f','d') */
19     double tthresh = .05; /* relative acceptance threshold */
20    
21     #define dval3(ix,ox,oy) datarr[((((ix)<<log2g)+(ox))<<log2g)+(oy)]
22     #define dval4(ix,iy,ox,oy) datarr[((((((ix)<<log2g)+(iy))<<log2g)+(ox))<<log2g)+(oy)]
23    
24     #define above_threshold(tp) ((tp)->vmax - (tp)->vmin > 2.*tthresh*(tp)->vavg)
25    
26     /* Tensor tree node */
27     typedef struct ttree_s {
28     float vmin, vmax; /* value extrema */
29     float vavg; /* average */
30     struct ttree_s *kid; /* 2^ttrank children */
31     } TNODE;
32    
33     /* Allocate a new set of children for the given node (no checks) */
34     static void
35     new_kids(TNODE *pn)
36     {
37     pn->kid = (TNODE *)calloc(1<<ttrank, sizeof(TNODE));
38     if (pn->kid == NULL)
39     error(SYSTEM, "out of memory in new_kids");
40     }
41    
42     /* Free children for this node */
43     static void
44     free_kids(TNODE *pn)
45     {
46     int i;
47    
48     if (pn->kid == NULL)
49     return;
50     for (i = 1<<ttrank; i--; )
51     free_kids(pn->kid+i);
52     free(pn->kid);
53     pn->kid = NULL;
54     }
55    
56     /* Build a tensor tree starting from the given hypercube */
57     static void
58     build_tree(TNODE *tp, const int bmin[], int l2s)
59     {
60     int bkmin[4];
61     int i, j;
62    
63     tp->vmin = 1e20;
64     tp->vmax = 0;
65     tp->vavg = 0;
66     if (l2s <= 1) { /* reached upper leaves */
67     for (i = 1<<ttrank; i--; ) {
68     float val;
69     for (j = ttrank; j--; )
70     bkmin[j] = bmin[j] + (i>>j & 1);
71     val = (ttrank == 3) ? dval3(bkmin[0],bkmin[1],bkmin[2])
72     : dval4(bkmin[0],bkmin[1],bkmin[2],bkmin[3]);
73     if (val < tp->vmin)
74     tp->vmin = val;
75     if (val > tp->vmax)
76     tp->vmax = val;
77     tp->vavg += val;
78     }
79     tp->vavg /= (float)(1<<ttrank);
80     return;
81     }
82     --l2s; /* else still branching */
83     new_kids(tp); /* grow recursively */
84     for (i = 1<<ttrank; i--; ) {
85     for (j = ttrank; j--; )
86     bkmin[j] = bmin[j] + ((i>>j & 1)<<l2s);
87     build_tree(tp->kid+i, bkmin, l2s);
88     if (tp->kid[i].vmin < tp->vmin)
89     tp->vmin = tp->kid[i].vmin;
90     if (tp->kid[i].vmax > tp->vmax)
91     tp->vmax = tp->kid[i].vmax;
92     tp->vavg += tp->kid[i].vavg;
93     }
94     tp->vavg /= (float)(1<<ttrank);
95     /* is variation above threshold? */
96     if (!above_threshold(tp))
97     free_kids(tp); /* if not, trim branches */
98     }
99    
100     /* Print a tensor tree from the given hypercube */
101     static void
102     print_tree(const TNODE *tp, const int bmin[], int l2s)
103     {
104     int bkmin[4];
105     int i, j;
106    
107     for (j = log2g-l2s; j--; ) /* indent based on branch level */
108     fputs(" ", stdout);
109     fputc('{', stdout); /* special case for upper leaves */
110     if (l2s <= 1 && above_threshold(tp)) {
111     for (i = 0; i < 1<<ttrank; i++) {
112     float val;
113     for (j = ttrank; j--; )
114     bkmin[j] = bmin[j] + (i>>j & 1);
115     val = (ttrank == 3) ? dval3(bkmin[0],bkmin[1],bkmin[2])
116     : dval4(bkmin[0],bkmin[1],bkmin[2],bkmin[3]);
117     printf(" %.4e", val);
118     }
119     fputs(" }\n", stdout);
120     return;
121     }
122     if (tp->kid == NULL) { /* trimmed limb */
123     printf(" %.6e }\n", tp->vavg);
124     return;
125     }
126     --l2s; /* else still branching */
127     fputc('\n', stdout);
128     for (i = 0; i < 1<<ttrank; i++) {
129     for (j = ttrank; j--; )
130     bkmin[j] = bmin[j] + ((i>>j & 1)<<l2s);
131     print_tree(tp->kid+i, bkmin, l2s);
132     }
133     ++l2s;
134     for (j = log2g-l2s; j--; )
135     fputs(" ", stdout);
136     fputs("}\n", stdout);
137     }
138    
139     /* Read a row of data in ASCII */
140     static int
141     read_ascii(float *rowp, int n)
142     {
143     int n2go;
144    
145     if ((rowp == NULL) | (n <= 0))
146     return(0);
147     for (n2go = n; n2go; n2go--)
148     if (scanf("%f", rowp++) != 1)
149     break;
150     if (n2go)
151     error(USER, "unexpected EOD on ascii input");
152     return(n-n2go);
153     }
154    
155     /* Read a row of float data */
156     static int
157     read_float(float *rowp, int n)
158     {
159     int nread;
160    
161     if ((rowp == NULL) | (n <= 0))
162     return(0);
163     nread = fread(rowp, sizeof(float), n, stdin);
164     if (nread != n)
165     error(USER, "unexpected EOF on float input");
166     return(nread);
167     }
168    
169     /* Read a row of double data */
170     static int
171     read_double(float *rowp, int n)
172     {
173     static double *rowbuf = NULL;
174     static int rblen = 0;
175     int nread, i;
176    
177     if ((rowp == NULL) | (n <= 0)) {
178     if (rblen) {
179     free(rowbuf);
180     rowbuf = NULL; rblen = 0;
181     }
182     return(0);
183     }
184     if (rblen < n) {
185     rowbuf = (double *)realloc(rowbuf, sizeof(double)*(rblen=n));
186     if (rowbuf == NULL)
187     error(SYSTEM, "out of memory in read_double");
188     }
189     nread = fread(rowbuf, sizeof(double), n, stdin);
190     if (nread != n)
191     error(USER, "unexpected EOF on double input");
192     for (i = 0; i < nread; i++)
193     *rowp++ = rowbuf[i];
194     return(nread);
195     }
196    
197     /* Load data array, filling zeroes for rank 3 demi-tensor */
198     static void
199     load_data()
200     {
201     int (*readf)(float *, int) = NULL;
202    
203     switch (infmt) {
204     case 'a':
205     readf = &read_ascii;
206     break;
207     case 'f':
208     readf = &read_float;
209     break;
210     case 'd':
211     readf = &read_double;
212     break;
213     default:
214     error(COMMAND, "unsupported input format");
215     break;
216     }
217     datarr = (float *)calloc(1<<(log2g*ttrank), sizeof(float));
218     if (datarr == NULL)
219     error(SYSTEM, "out of memory in load_data");
220     if (ttrank == 3) {
221     int ix, ox;
222     for (ix = 0; ix < 1<<log2g; ix++)
223     for (ox = 0; ox < 1<<log2g; ox++)
224     (*readf)(datarr+((((ix)<<log2g)+(ox))<<log2g),
225     1<<(log2g-1));
226     } else /* ttrank == 4 */ {
227     int ix, iy, ox;
228     for (ix = 0; ix < 1<<log2g; ix++)
229     for (iy = 0; iy < 1<<log2g; iy++)
230     for (ox = 0; ox < 1<<log2g; ox++)
231     (*readf)(datarr +
232     ((((((ix)<<log2g)+(iy))<<log2g)+(ox))<<log2g),
233     1<<log2g);
234     }
235     (*readf)(NULL, 0); /* releases any buffers */
236     if (infmt == 'a') {
237     int c;
238     while ((c = getc(stdin)) != EOF) {
239     switch (c) {
240     case ' ':
241     case '\t':
242     case '\r':
243     case '\n':
244     continue;
245     }
246     error(WARNING, "data past end of expected input");
247     break;
248     }
249     } else if (getc(stdin) != EOF)
250     error(WARNING, "binary data past end of expected input");
251     }
252    
253     /* Load BSDF array, coalesce uniform regions and format as tensor tree */
254     int
255     main(int argc, char *argv[])
256     {
257     int doheader = 1;
258     int bmin[4];
259     TNODE gtree;
260     int i;
261     /* get options and parameters */
262     for (i = 1; i < argc && argv[i][0] == '-'; i++)
263     switch (argv[i][1]) {
264     case 'h':
265     doheader = !doheader;
266     break;
267     case 'r':
268     ttrank = atoi(argv[++i]);
269     if (ttrank != 3 && ttrank != 4)
270     goto userr;
271     break;
272     case 'g':
273     log2g = atoi(argv[++i]);
274     if (log2g <= 1)
275     goto userr;
276     break;
277     case 't':
278     tthresh = atof(argv[++i]);
279     if (tthresh <= 0)
280     goto userr;
281     break;
282     case 'f':
283     infmt = argv[i][2];
284     if (!infmt || strchr("afd", infmt) == NULL)
285     goto userr;
286     break;
287     default:
288     goto userr;
289     }
290     if (i < argc-1)
291     goto userr;
292     /* load input data */
293     if (i == argc-1 && freopen(argv[i], "rb", stdin) == NULL) {
294     sprintf(errmsg, "cannot open input file '%s'", argv[i]);
295     error(SYSTEM, errmsg);
296     }
297     if (infmt != 'a')
298     SET_FILE_BINARY(stdin);
299     load_data();
300     if (doheader) {
301     for (i = 0; i < argc; i++) {
302     fputs(argv[i], stdout);
303     fputc(i < argc-1 ? ' ' : '\n', stdout);
304     }
305     fputc('\n', stdout);
306     }
307     gtree.kid = NULL; /* create our tree */
308     bmin[0] = bmin[1] = bmin[2] = bmin[3] = 0;
309     build_tree(&gtree, bmin, log2g);
310     /* format to stdout */
311     print_tree(&gtree, bmin, log2g);
312     /* Clean up isn't necessary for main()...
313     free_kids(&gtree);
314     free(datarr);
315     */
316     return(0);
317     userr:
318     fprintf(stderr, "Usage: %s [-h][-f{a|f|d}][-r {3|4}][-g log2grid][-t thresh] [input]\n",
319     argv[0]);
320     return(1);
321     }