ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rttree_reduce.c
Revision: 2.20
Committed: Tue Jun 30 17:18:22 2020 UTC (3 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R4, rad5R3, HEAD
Changes since 2.19: +5 -5 lines
Log Message:
feat(rttree_reduce): made header conform to more standard practice

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: rttree_reduce.c,v 2.19 2020/05/11 20:26:12 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[((size_t)(((ix)<<log2g)+(ox))<<log2g)+(oy)]
23 #define dval4(ix,iy,ox,oy) datarr[((((size_t)(((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 4. /* 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) sqrt( ((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((size_t)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((0.001<=val)&(val<10.) ? " %.7f" : " %.3e", 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 = getbinary(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 if (rblen) free(rowbuf);
227 rowbuf = (double *)malloc(sizeof(double)*(rblen=n));
228 if (rowbuf == NULL)
229 error(SYSTEM, "out of memory in read_double");
230 }
231 nread = getbinary(rowbuf, sizeof(double), n, stdin);
232 if (nread != n)
233 error(USER, "unexpected EOF on double input");
234 for (i = 0; i < nread; i++)
235 *rowp++ = rowbuf[i];
236 return(nread);
237 }
238
239 /* Truncate any negative values to zero */
240 static void
241 noneg(float *varr, int n)
242 {
243 int nnan = 0;
244
245 while (n-- > 0) {
246 #ifdef isnan
247 if (isnan(*varr)) {
248 *varr = 0;
249 ++nnan;
250 } else
251 #endif
252 if (*varr < 0) *varr = 0;
253 ++varr;
254 }
255 if (nnan)
256 fprintf(stderr, "Warning: BSDF data contains %d NaN values\n",
257 nnan);
258 }
259
260 /* Load data array, filling zeroes for rank 3 demi-tensor */
261 static void
262 load_data(void)
263 {
264 int (*readf)(float *, int) = NULL;
265
266 switch (infmt) {
267 case 'a':
268 readf = &read_ascii;
269 break;
270 case 'f':
271 readf = &read_float;
272 break;
273 case 'd':
274 readf = &read_double;
275 break;
276 default:
277 error(COMMAND, "unsupported input format");
278 break;
279 }
280 datarr = (float *)calloc((size_t)1<<(log2g*ttrank), sizeof(float));
281 if (datarr == NULL)
282 error(SYSTEM, "out of memory in load_data");
283 if (ttrank == 3) {
284 int ix, ox;
285 for (ix = 0; ix < 1<<(log2g-1); ix++)
286 for (ox = 0; ox < 1<<log2g; ox++)
287 (*readf)(&dval3(ix,ox,0), 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)(&dval4(ix,iy,ox,0), 1<<log2g);
294 }
295 (*readf)(NULL, 0); /* releases any buffers */
296 if (infmt == 'a') {
297 int c;
298 while ((c = getc(stdin)) != EOF) {
299 switch (c) {
300 case ' ':
301 case '\t':
302 case '\r':
303 case '\n':
304 continue;
305 }
306 error(WARNING, "data past end of expected input");
307 break;
308 }
309 } else if (getc(stdin) != EOF)
310 error(WARNING, "binary data past end of expected input");
311
312 noneg(datarr, 1<<(log2g*ttrank)); /* take precautions */
313 }
314
315 /* Enforce reciprocity by averaging data values */
316 static void
317 do_reciprocity(void)
318 {
319 const int siz = 1<<log2g;
320 float *v1p, *v2p;
321
322 if (ttrank == 3) {
323 int ix, ox, oy;
324 for (ix = 0; ix < siz>>1; ix++)
325 for (ox = 0; ox < siz; ox++)
326 for (oy = 0; oy < siz>>1; oy++) {
327 v1p = &dval3(ix,ox,oy);
328 v2p = &dval3(ix,ox,siz-1-oy);
329 *v1p = *v2p = .5f*( *v1p + *v2p );
330 }
331 } else /* ttrank == 4 */ {
332 int ix, iy, ox, oy;
333 for (ix = 0; ix < siz; ix++)
334 for (iy = 0; iy < siz; iy++) {
335 int cnt = ix*siz + iy;
336 for (ox = 0; cnt > 0; ox++)
337 for (oy = 0; oy < siz; oy++) {
338 v1p = &dval4(siz-1-ix,siz-1-iy,ox,oy);
339 v2p = &dval4(siz-1-ox,siz-1-oy,ix,iy);
340 *v1p = *v2p = .5f*( *v1p + *v2p );
341 if (--cnt <= 0)
342 break;
343 }
344 }
345 }
346 }
347
348 /* Load BSDF array, coalesce uniform regions and format as tensor tree */
349 int
350 main(int argc, char *argv[])
351 {
352 int doheader = 1;
353 int recipavg = 0;
354 int bmin[4];
355 TNODE gtree;
356 int i;
357 /* get options and parameters */
358 for (i = 1; i < argc && argv[i][0] == '-'; i++)
359 switch (argv[i][1]) {
360 case 'a':
361 recipavg = !recipavg;
362 break;
363 case 'h':
364 doheader = !doheader;
365 break;
366 case 'r':
367 ttrank = atoi(argv[++i]);
368 if (ttrank != 3 && ttrank != 4)
369 goto userr;
370 break;
371 case 'g':
372 log2g = atoi(argv[++i]);
373 if (log2g <= 1)
374 goto userr;
375 break;
376 case 't':
377 pctcull = atof(argv[++i]);
378 if ((pctcull < 0) | (pctcull >= 100.))
379 goto userr;
380 break;
381 case 'f':
382 infmt = argv[i][2];
383 if (!infmt || strchr("afd", infmt) == NULL)
384 goto userr;
385 break;
386 default:
387 goto userr;
388 }
389 if (i < argc-1)
390 goto userr;
391 /* load input data */
392 if (i == argc-1 && freopen(argv[i], "r", stdin) == NULL) {
393 sprintf(errmsg, "cannot open input file '%s'", argv[i]);
394 error(SYSTEM, errmsg);
395 }
396 if (infmt != 'a')
397 SET_FILE_BINARY(stdin);
398 #ifdef getc_unlocked /* avoid lock/unlock overhead */
399 flockfile(stdin);
400 #endif
401 load_data();
402 if (recipavg)
403 do_reciprocity();
404 if (doheader) {
405 newheader("RADIANCE", stdout);
406 printargs(argc, argv, stdout);
407 fputnow(stdout);
408 fputformat("ascii", stdout);
409 fputc('\n', stdout);
410 }
411 gtree.kid = NULL; /* create our tree */
412 bmin[0] = bmin[1] = bmin[2] = bmin[3] = 0;
413 build_tree(&gtree, bmin, log2g);
414 /* compute threshold & trim tree */
415 set_threshold();
416 trim_tree(&gtree);
417 /* format to stdout */
418 print_tree(&gtree, bmin, log2g);
419 /* Clean up isn't necessary for main()...
420 free_kids(&gtree);
421 free(datarr);
422 */
423 return(0);
424 userr:
425 fprintf(stderr, "Usage: %s [-h][-a][-f{a|f|d}][-r {3|4}][-g log2grid][-t trim%%] [input]\n",
426 argv[0]);
427 return(1);
428 }