ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rttree_reduce.c
Revision: 2.2
Committed: Tue May 31 20:50:26 2011 UTC (12 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.1: +50 -9 lines
Log Message:
Fixed bug in missing persist file and changed to percentile-based threshold

File Contents

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