| 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 |
< |
float *v1p, *v2p; |
| 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 < 1<<(log2g-1); ix++) |
| 325 |
< |
for (ox = 0; ox < 1<<log2g; ox++) |
| 326 |
< |
for (oy = 0; oy < 1<<(log2g-1); 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,(1<<log2g)-1-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 = 1; ix < 1<<log2g; ix++) |
| 334 |
< |
for (iy = 1; iy < 1<<log2g; iy++) |
| 335 |
< |
for (ox = 0; ox < ix; ox++) |
| 336 |
< |
for (oy = 0; oy < iy; oy++) { |
| 337 |
< |
v1p = &dval4(ix,iy,ox,oy); |
| 338 |
< |
v2p = &dval4(ox,oy,ix,iy); |
| 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(); |