search

BIC-MNI / minc-tools

Basic minc-tools from former minc repository
Other
29 stars 25 forks source link

mincstats reports negative variance and -nan for stddev as a result #79

Closed gdevenyi closed 6 years ago

gdevenyi commented 6 years ago

Looking at the ratio of two files generated by minccalc like this that are very very similar:

minccalc -quiet -clobber -expression 'A[0]/A[1]' /tmp/tmp.HeE0mf4YTk/0/bias.mnc /tmp/tmp.HeE0mf4YTk/1/bias.mnc /tmp/tmp.HeE0mf4YTk/1/ratio.mnc

Computing statistics on the ratio, mincstats breaks down:

mincstats /tmp/tmp.HeE0mf4YTk/1/ratio.mnc 
File:              /tmp/tmp.HeE0mf4YTk/1/ratio.mnc
Mask file:         (null)
Total voxels:      6759288
# voxels:          6759288
% of total:        100
Volume (mm3):      6759288
Min:               0.999998033
Max:               1.000001669
Sum:               6759287.534
Sum^2:             6759287.069
Mean:              0.9999999311
Variance:          -7.57812793e-15
Stddev:            -nan
CoM_voxel(z,y,x):  91.00000626 113.500004 80.49999674
CoM_real(x,y,z):   0.4428688544 -14.97257385 9.441359867

Histogram:         (null)
Total voxels:      6759288
# voxels:          6759288
% of total:        100
Median:            0.9999999974
Majority:          1.000000059
BiModalT:          0.9999995828
PctT [  0%]:       0.999998034
Entropy :          4.498000219
bbbxyz commented 6 years ago

Seems like this is caused by cancellation when calculating the variance . It only seems likely to become an issue in volumes that are uniform and close to 1, since the Sum^2 is close to Sum * Sum / #voxels, which might explain why the bug went unnoticed for so long. See: https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Na%C3%AFve_algorithm

I'll work on fixing this using value shifting like suggested in the wiki article.

gdevenyi commented 6 years ago

Wow, this is the second naive numerical implementation I've run into this month :)

Thanks for looking into this!

andrewjanke commented 6 years ago

Indeed a good catch! In my defense when I wrote mincstats (circa 2002) that wikipedia page looked a lot different! :)

gdevenyi commented 6 years ago

We'll overlook it this time @andrewjanke, but you're on thin ice as far as bugs in the 100+ minctools you wrote from scratch.... ;)

gdevenyi commented 6 years ago

Solved by #83