Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

mean of 3d nifti does not match calc by fslstats and 3dBrickStat #70

Open
WillForan opened this issue Oct 24, 2023 · 2 comments
Open

mean of 3d nifti does not match calc by fslstats and 3dBrickStat #70

WillForan opened this issue Oct 24, 2023 · 2 comments

Comments

@WillForan
Copy link

Thanks for this package! I'm excited to be able to use Julia in my analysis.

I'm sanity checking and measuring performance here, with this nifti.

The mean of that 3d dataset is 838.206399 (according to fsl and AFNI tools), but mean(niread('mynii.nii.gz')) reports 836.1491.

I'm not sure how to narrow down where the discrepancy is coming from.

fslstats wf-mp2rage-7t_2017087.nii.gz -m
 838.206399 

3dBrickStat -slow -mean wf-mp2rage-7t_2017087.nii.gz
 838.206       


scripts/niimean.jl
 836.1491

cat scripts/niimean.jl 
#!/usr/bin/env julia
using Statistics
using NIfTI
ni = niread("wf-mp2rage-7t_2017087.nii.gz");
print(mean(ni))
@WillForan
Copy link
Author

WillForan commented Oct 25, 2023

Looks like a type issue. Not sure if there's anything NIfTI.jl could do about it.
Is it worth mentioning in the readme? Or is this a common footgun/julia knowledge?

# naive attempt, no care about types. loop matches much faster Statistics::mean 
sum=0; for x in ni; sum=sum+x; end;
print(sum, "\n", sum/prod(size(ni)),"\n")
# 1.5604549e10
# 836.1491

# make sure accumulator is a BigInt, get the expected value.
sum=BigInt(0); for x in ni; sum=sum+x; end; print(sum, "\n", sum/prod(size(ni)),"\n")
# 1.5642943107e+10
# 838.2063993377057613168724279835390946502057613168724279835390946502057613168706

# force each element to bigint before meaning also works (and is slow)
mean(BigInt,ni)
# 838.2063


# scaling input before mean also works and is much faster than the other methods above
mean(ni/1000)*1000
# 838.2065f0

@Tokazama
Copy link
Member

A lot of images store in Int16 so this comes up. You also are looking at the mean of the entire array instead of just the brain unless you average only masked regions explicitly

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants