39 lines
1 KiB
Python
Executable file
39 lines
1 KiB
Python
Executable file
#!/usr/bin/python
|
|
"""
|
|
Compute a tSNR volume given an input time series image.
|
|
"""
|
|
|
|
import os
|
|
import sys
|
|
|
|
if not len(sys.argv) == 2:
|
|
print __doc__
|
|
sys.exit(1)
|
|
|
|
import nibabel as nb
|
|
from mvpa2.mappers.detrend import poly_detrend
|
|
from mvpa2.datasets.mri import fmri_dataset
|
|
from mvpa2.datasets.mri import map2nifti
|
|
|
|
infilename = sys.argv[1]
|
|
outfilename = os.path.splitext(infilename)[0]
|
|
if outfilename.endswith('.nii'):
|
|
outfilename = os.path.splitext(outfilename)[0]
|
|
outfilename += '_tsnr.nii.gz'
|
|
|
|
ds = fmri_dataset(infilename)
|
|
# clip the value range, can only be negative due to interpolation
|
|
# regime
|
|
ds.samples[ds.samples < 0] = 0
|
|
# first capture mean -- detrending will bring it close to zero
|
|
m = ds.samples.mean(axis=0)
|
|
# now detrend
|
|
poly_detrend(ds, polyord=1)
|
|
# and get the std after linear trend removal
|
|
std = ds.samples.std(axis=0)
|
|
# compute tSNR
|
|
m[std > 0] /= std[std > 0]
|
|
# protect against ultra-low noise voxels (images are zero'ed at the edges)
|
|
m[m > 1000] = 1000
|
|
|
|
map2nifti(ds, m).to_filename(outfilename)
|