ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rttree_reduce.c
(Generate patch)

Comparing ray/src/util/rttree_reduce.c (file contents):
Revision 2.1 by greg, Thu May 26 15:32:02 2011 UTC vs.
Revision 2.13 by schorsch, Sun Mar 6 01:13:18 2016 UTC

# Line 11 | Line 11 | static const char RCSid[] = "$Id$";
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  tthresh = .05;          /* relative acceptance threshold */
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  
24 – #define above_threshold(tp)     ((tp)->vmax - (tp)->vmin > 2.*tthresh*(tp)->vavg)
25 –
25   /* Tensor tree node */
26   typedef struct ttree_s {
27          float           vmin, vmax;     /* value extrema */
# Line 30 | Line 29 | typedef struct ttree_s {
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 +        /* XXX VC warns about 32 bit shift coerced to 64 bit */
48          pn->kid = (TNODE *)calloc(1<<ttrank, sizeof(TNODE));
49          if (pn->kid == NULL)
50                  error(SYSTEM, "out of memory in new_kids");
# Line 77 | Line 88 | build_tree(TNODE *tp, const int bmin[], int l2s)
88                          tp->vavg += val;
89                  }
90                  tp->vavg /= (float)(1<<ttrank);
91 +                                        /* record stats */
92 +                i = (HISTLEN/HISTMAX) * var_measure(tp);
93 +                if (i >= HISTLEN) i = HISTLEN-1;
94 +                ++histo[i];
95                  return;
96          }
97          --l2s;                          /* else still branching */
# Line 92 | Line 107 | build_tree(TNODE *tp, const int bmin[], int l2s)
107                  tp->vavg += tp->kid[i].vavg;
108          }
109          tp->vavg /= (float)(1<<ttrank);
95 –                                        /* is variation above threshold? */
96 –        if (!above_threshold(tp))
97 –                free_kids(tp);          /* if not, trim branches */
110   }
111  
112 + /* Set our trimming threshold */
113 + static void
114 + set_threshold(void)
115 + {
116 +        int     hsum = 0;
117 +        int     i;
118 +
119 +        for (i = HISTLEN; i--; )
120 +                hsum += histo[i];
121 +        hsum = pctcull*.01 * (double)hsum;
122 +        for (i = 0; hsum > 0; i++)
123 +                hsum -= histo[i];
124 +        tthresh = (HISTMAX/HISTLEN) * i;
125 + }
126 +
127 + /* Trim our tree according to the current threshold */
128 + static void
129 + trim_tree(TNODE *tp)
130 + {
131 +        if (tp->kid == NULL)
132 +                return;
133 +        if (above_threshold(tp)) {      /* keeping branches? */
134 +                int     i = 1<<ttrank;
135 +                while (i--)
136 +                        trim_tree(tp->kid+i);
137 +                return;
138 +        }
139 +        free_kids(tp);                  /* else trim at this point */
140 + }
141 +
142   /* Print a tensor tree from the given hypercube */
143   static void
144   print_tree(const TNODE *tp, const int bmin[], int l2s)
# Line 111 | Line 153 | print_tree(const TNODE *tp, const int bmin[], int l2s)
153                  for (i = 0; i < 1<<ttrank; i++) {
154                          float   val;
155                          for (j = ttrank; j--; )
156 <                                bkmin[j] = bmin[j] + (i>>j & 1);
156 >                                bkmin[j] = bmin[j] + (i>>(ttrank-1-j) & 1);
157                          val = (ttrank == 3) ? dval3(bkmin[0],bkmin[1],bkmin[2])
158                                  : dval4(bkmin[0],bkmin[1],bkmin[2],bkmin[3]);
159 <                        printf(" %.4e", val);
159 >                        printf((0.001<=val)&(val<10.) ? " %.7f" : " %.3e", val);
160                  }
161                  fputs(" }\n", stdout);
162                  return;
# Line 194 | Line 236 | read_double(float *rowp, int n)
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()
262 > load_data(void)
263   {
264          int     (*readf)(float *, int) = NULL;
265          
# Line 214 | Line 277 | load_data()
277                  error(COMMAND, "unsupported input format");
278                  break;
279          }
280 +        /* XXX VC warns about 32 bit shift coerced to 64 bit */
281          datarr = (float *)calloc(1<<(log2g*ttrank), sizeof(float));
282          if (datarr == NULL)
283                  error(SYSTEM, "out of memory in load_data");
284          if (ttrank == 3) {
285                  int     ix, ox;
286 <                for (ix = 0; ix < 1<<log2g; ix++)
286 >                for (ix = 0; ix < 1<<(log2g-1); ix++)
287                          for (ox = 0; ox < 1<<log2g; ox++)
288 <                                (*readf)(datarr+((((ix)<<log2g)+(ox))<<log2g),
289 <                                                1<<(log2g-1));
288 >                                (*readf)(datarr+(((ix<<log2g)+ox)<<log2g),
289 >                                                1<<log2g);
290          } else /* ttrank == 4 */ {
291                  int     ix, iy, ox;
292                  for (ix = 0; ix < 1<<log2g; ix++)
293                      for (iy = 0; iy < 1<<log2g; iy++)
294                          for (ox = 0; ox < 1<<log2g; ox++)
295                                  (*readf)(datarr +
296 <                                ((((((ix)<<log2g)+(iy))<<log2g)+(ox))<<log2g),
296 >                                (((((ix<<log2g)+iy)<<log2g)+ox)<<log2g),
297                                                  1<<log2g);
298          }
299          (*readf)(NULL, 0);      /* releases any buffers */
# Line 248 | Line 312 | load_data()
312                  }
313          } else if (getc(stdin) != EOF)
314                  error(WARNING, "binary data past end of expected input");
315 +
316 +        noneg(datarr, 1<<(log2g*ttrank));       /* take precautions */
317   }
318  
319 + /* Enforce reciprocity by averaging data values */
320 + static void
321 + do_reciprocity(void)
322 + {
323 +        const int       siz = 1<<log2g;
324 +        float           *v1p, *v2p;
325 +
326 +        if (ttrank == 3) {
327 +                int     ix, ox, oy;
328 +                for (ix = 0; ix < siz>>1; ix++)
329 +                    for (ox = 0; ox < siz; ox++)
330 +                        for (oy = 0; oy < siz>>1; oy++) {
331 +                                v1p = &dval3(ix,ox,oy);
332 +                                v2p = &dval3(ix,ox,siz-1-oy);
333 +                                *v1p = *v2p = .5f*( *v1p + *v2p );
334 +                        }
335 +        } else /* ttrank == 4 */ {
336 +                int     ix, iy, ox, oy;
337 +                for (ix = 1; ix < siz; ix++)
338 +                    for (iy = 1; iy < siz; iy++)
339 +                        for (ox = 0; ox < ix; ox++)
340 +                            for (oy = 0; oy < iy; oy++) {
341 +                                v1p = &dval4(siz-1-ix,siz-1-iy,ox,oy);
342 +                                v2p = &dval4(siz-1-ox,siz-1-oy,ix,iy);
343 +                                *v1p = *v2p = .5f*( *v1p + *v2p );
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;
# Line 275 | Line 374 | main(int argc, char *argv[])
374                                  goto userr;
375                          break;
376                  case 't':
377 <                        tthresh = atof(argv[++i]);
378 <                        if (tthresh <= 0)
377 >                        pctcull = atof(argv[++i]);
378 >                        if ((pctcull < 0) | (pctcull >= 100.))
379                                  goto userr;
380                          break;
381                  case 'f':
# Line 290 | Line 389 | main(int argc, char *argv[])
389          if (i < argc-1)
390                  goto userr;
391                                          /* load input data */
392 <        if (i == argc-1 && freopen(argv[i], "rb", stdin) == NULL) {
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                  for (i = 0; i < argc; i++) {
406                          fputs(argv[i], stdout);
# Line 307 | Line 411 | main(int argc, char *argv[])
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()...
# Line 315 | Line 422 | main(int argc, char *argv[])
422          */
423          return(0);
424   userr:
425 <        fprintf(stderr, "Usage: %s [-h][-f{a|f|d}][-r {3|4}][-g log2grid][-t thresh] [input]\n",
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   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines