ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rttree_reduce.c
Revision: 2.10
Committed: Tue Mar 18 01:33:07 2014 UTC (10 years ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R2P2, rad4R2, rad4R2P1
Changes since 2.9: +27 -4 lines
Log Message:
Added protection against negative and NaN values on input

File Contents

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