| 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)] | 
| 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 { | 
| 30 |  | } TNODE; | 
| 31 |  |  | 
| 32 |  | #define HISTLEN         300     /* histogram resolution */ | 
| 33 | < | #define HISTMAX         10.     /* maximum recorded measure in histogram */ | 
| 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)         ( ((tp)->vmax - (tp)->vmin) / \ | 
| 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 |  |  | 
| 44 |  | static void | 
| 45 |  | new_kids(TNODE *pn) | 
| 46 |  | { | 
| 47 | < | pn->kid = (TNODE *)calloc(1<<ttrank, sizeof(TNODE)); | 
| 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 |  | } | 
| 110 |  |  | 
| 111 |  | /* Set our trimming threshold */ | 
| 112 |  | static void | 
| 113 | < | set_threshold() | 
| 113 | > | set_threshold(void) | 
| 114 |  | { | 
| 115 |  | int     hsum = 0; | 
| 116 |  | int     i; | 
| 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); | 
| 158 | > | printf((0.001<=val)&(val<10.) ? " %.7f" : " %.3e", val); | 
| 159 |  | } | 
| 160 |  | fputs(" }\n", stdout); | 
| 161 |  | return; | 
| 201 |  |  | 
| 202 |  | if ((rowp == NULL) | (n <= 0)) | 
| 203 |  | return(0); | 
| 204 | < | nread = fread(rowp, sizeof(float), n, stdin); | 
| 204 | > | nread = getbinary(rowp, sizeof(float), n, stdin); | 
| 205 |  | if (nread != n) | 
| 206 |  | error(USER, "unexpected EOF on float input"); | 
| 207 |  | return(nread); | 
| 223 |  | return(0); | 
| 224 |  | } | 
| 225 |  | if (rblen < n) { | 
| 226 | < | rowbuf = (double *)realloc(rowbuf, sizeof(double)*(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 = fread(rowbuf, sizeof(double), n, stdin); | 
| 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++) | 
| 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 |  |  | 
| 277 |  | error(COMMAND, "unsupported input format"); | 
| 278 |  | break; | 
| 279 |  | } | 
| 280 | < | datarr = (float *)calloc(1<<(log2g*ttrank), sizeof(float)); | 
| 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)(datarr+(((ix<<log2g)+ox)<<log2g), | 
| 266 | < | 1<<log2g); | 
| 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)(datarr + | 
| 273 | < | (((((ix<<log2g)+iy)<<log2g)+ox)<<log2g), | 
| 274 | < | 1<<log2g); | 
| 293 | > | (*readf)(&dval4(ix,iy,ox,0), 1<<log2g); | 
| 294 |  | } | 
| 295 |  | (*readf)(NULL, 0);      /* releases any buffers */ | 
| 296 |  | if (infmt == 'a') { | 
| 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() | 
| 317 | > | do_reciprocity(void) | 
| 318 |  | { | 
| 319 |  | const int       siz = 1<<log2g; | 
| 320 |  | float           *v1p, *v2p; | 
| 330 |  | } | 
| 331 |  | } else /* ttrank == 4 */ { | 
| 332 |  | int     ix, iy, ox, oy; | 
| 333 | < | for (ix = 1; ix < siz; ix++) | 
| 334 | < | for (iy = 1; iy < siz; iy++) | 
| 335 | < | for (ox = 0; ox < ix; ox++) | 
| 336 | < | for (oy = 0; oy < iy; 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 |  |  | 
| 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(); |