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

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.20 static const char RCSid[] = "$Id: rttree_reduce.c,v 2.19 2020/05/11 20:26:12 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 greg 2.3 #include <math.h>
15 greg 2.1
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 greg 2.4 double pctcull = 95.; /* target culling percentile */
21 greg 2.2
22 greg 2.19 #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 greg 2.1
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 greg 2.3 #define HISTLEN 300 /* histogram resolution */
33 greg 2.18 #define HISTMAX 4. /* maximum recorded measure in histogram */
34 greg 2.3
35     int histo[HISTLEN]; /* histogram freq. of variance measure */
36    
37     double tthresh; /* acceptance threshold (TBD) */
38    
39 greg 2.18 #define var_measure(tp) sqrt( ((tp)->vmax - (tp)->vmin) / \
40 greg 2.3 (sqrt((tp)->vavg) + .03) )
41     #define above_threshold(tp) (var_measure(tp) > tthresh)
42    
43 greg 2.1 /* Allocate a new set of children for the given node (no checks) */
44     static void
45     new_kids(TNODE *pn)
46     {
47 greg 2.14 pn->kid = (TNODE *)calloc((size_t)1<<ttrank, sizeof(TNODE));
48 greg 2.1 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 greg 2.2 /* record stats */
91 greg 2.3 i = (HISTLEN/HISTMAX) * var_measure(tp);
92 greg 2.2 if (i >= HISTLEN) i = HISTLEN-1;
93     ++histo[i];
94 greg 2.1 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 greg 2.2 }
110    
111     /* Set our trimming threshold */
112     static void
113 greg 2.10 set_threshold(void)
114 greg 2.2 {
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 greg 2.1 }
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 greg 2.5 bkmin[j] = bmin[j] + (i>>(ttrank-1-j) & 1);
156 greg 2.1 val = (ttrank == 3) ? dval3(bkmin[0],bkmin[1],bkmin[2])
157     : dval4(bkmin[0],bkmin[1],bkmin[2],bkmin[3]);
158 greg 2.12 printf((0.001<=val)&(val<10.) ? " %.7f" : " %.3e", val);
159 greg 2.1 }
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 greg 2.16 nread = getbinary(rowp, sizeof(float), n, stdin);
205 greg 2.1 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 greg 2.17 if (rblen) free(rowbuf);
227     rowbuf = (double *)malloc(sizeof(double)*(rblen=n));
228 greg 2.1 if (rowbuf == NULL)
229     error(SYSTEM, "out of memory in read_double");
230     }
231 greg 2.16 nread = getbinary(rowbuf, sizeof(double), n, stdin);
232 greg 2.1 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 greg 2.10 /* 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 greg 2.1 /* Load data array, filling zeroes for rank 3 demi-tensor */
261     static void
262 greg 2.10 load_data(void)
263 greg 2.1 {
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 schorsch 2.15 datarr = (float *)calloc((size_t)1<<(log2g*ttrank), sizeof(float));
281 greg 2.1 if (datarr == NULL)
282     error(SYSTEM, "out of memory in load_data");
283     if (ttrank == 3) {
284     int ix, ox;
285 greg 2.4 for (ix = 0; ix < 1<<(log2g-1); ix++)
286 greg 2.1 for (ox = 0; ox < 1<<log2g; ox++)
287 greg 2.17 (*readf)(&dval3(ix,ox,0), 1<<log2g);
288 greg 2.1 } 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 greg 2.17 (*readf)(&dval4(ix,iy,ox,0), 1<<log2g);
294 greg 2.1 }
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 greg 2.10
312     noneg(datarr, 1<<(log2g*ttrank)); /* take precautions */
313 greg 2.1 }
314    
315 greg 2.6 /* Enforce reciprocity by averaging data values */
316     static void
317 greg 2.10 do_reciprocity(void)
318 greg 2.6 {
319 greg 2.8 const int siz = 1<<log2g;
320     float *v1p, *v2p;
321 greg 2.6
322     if (ttrank == 3) {
323     int ix, ox, oy;
324 greg 2.8 for (ix = 0; ix < siz>>1; ix++)
325     for (ox = 0; ox < siz; ox++)
326     for (oy = 0; oy < siz>>1; oy++) {
327 greg 2.6 v1p = &dval3(ix,ox,oy);
328 greg 2.8 v2p = &dval3(ix,ox,siz-1-oy);
329 greg 2.6 *v1p = *v2p = .5f*( *v1p + *v2p );
330     }
331     } else /* ttrank == 4 */ {
332     int ix, iy, ox, oy;
333 greg 2.17 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 greg 2.8 v1p = &dval4(siz-1-ix,siz-1-iy,ox,oy);
339     v2p = &dval4(siz-1-ox,siz-1-oy,ix,iy);
340 greg 2.6 *v1p = *v2p = .5f*( *v1p + *v2p );
341 greg 2.17 if (--cnt <= 0)
342     break;
343 greg 2.6 }
344 greg 2.17 }
345 greg 2.6 }
346     }
347    
348 greg 2.1 /* 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 greg 2.6 int recipavg = 0;
354 greg 2.1 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 greg 2.6 case 'a':
361     recipavg = !recipavg;
362     break;
363 greg 2.1 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 greg 2.2 pctcull = atof(argv[++i]);
378     if ((pctcull < 0) | (pctcull >= 100.))
379 greg 2.1 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 greg 2.9 if (i == argc-1 && freopen(argv[i], "r", stdin) == NULL) {
393 greg 2.1 sprintf(errmsg, "cannot open input file '%s'", argv[i]);
394     error(SYSTEM, errmsg);
395     }
396     if (infmt != 'a')
397     SET_FILE_BINARY(stdin);
398 greg 2.9 #ifdef getc_unlocked /* avoid lock/unlock overhead */
399     flockfile(stdin);
400     #endif
401 greg 2.1 load_data();
402 greg 2.6 if (recipavg)
403     do_reciprocity();
404 greg 2.1 if (doheader) {
405 greg 2.20 newheader("RADIANCE", stdout);
406     printargs(argc, argv, stdout);
407     fputnow(stdout);
408     fputformat("ascii", stdout);
409 greg 2.1 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 greg 2.2 /* compute threshold & trim tree */
415     set_threshold();
416     trim_tree(&gtree);
417 greg 2.1 /* 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 greg 2.7 fprintf(stderr, "Usage: %s [-h][-a][-f{a|f|d}][-r {3|4}][-g log2grid][-t trim%%] [input]\n",
426 greg 2.1 argv[0]);
427     return(1);
428     }