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

# Content
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 }