| 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 = 99.;          /* target culling percentile */ | 
| 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)] | 
| 110 |  |  | 
| 111 |  | /* Set our trimming threshold */ | 
| 112 |  | static void | 
| 113 | < | set_threshold() | 
| 113 | > | set_threshold(void) | 
| 114 |  | { | 
| 115 |  | int     hsum = 0; | 
| 116 |  | int     i; | 
| 152 |  | for (i = 0; i < 1<<ttrank; i++) { | 
| 153 |  | float   val; | 
| 154 |  | for (j = ttrank; j--; ) | 
| 155 | < | bkmin[j] = bmin[j] + (i>>j & 1); | 
| 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); | 
| 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() | 
| 261 | > | load_data(void) | 
| 262 |  | { | 
| 263 |  | int     (*readf)(float *, int) = NULL; | 
| 264 |  |  | 
| 281 |  | error(SYSTEM, "out of memory in load_data"); | 
| 282 |  | if (ttrank == 3) { | 
| 283 |  | int     ix, ox; | 
| 284 | < | for (ix = 0; ix < 1<<log2g; ix++) | 
| 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-1)); | 
| 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), | 
| 294 | > | (((((ix<<log2g)+iy)<<log2g)+ox)<<log2g), | 
| 295 |  | 1<<log2g); | 
| 296 |  | } | 
| 297 |  | (*readf)(NULL, 0);      /* releases any buffers */ | 
| 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; | 
| 387 |  | if (i < argc-1) | 
| 388 |  | goto userr; | 
| 389 |  | /* load input data */ | 
| 390 | < | if (i == argc-1 && freopen(argv[i], "rb", stdin) == NULL) { | 
| 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); | 
| 420 |  | */ | 
| 421 |  | return(0); | 
| 422 |  | userr: | 
| 423 | < | fprintf(stderr, "Usage: %s [-h][-f{a|f|d}][-r {3|4}][-g log2grid][-t trim%%] [input]\n", | 
| 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 |  | } |