Skip to content

enh: multiple masks for compcor #1968

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.

Sign up for GitHub

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

Merged
merged 12 commits into from
Apr 28, 2017
Merged

enh: multiple masks for compcor #1968

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
adb3623
fix: allow multiple masks for interface
mgxd Mar 8, 2017
5bf6aec
Merge branch 'master' of github.com:nipy/nipype into fix/compcor
mgxd Apr 19, 2017
2dd15e7
Merge branch 'master' of github.com:nipy/nipype into fix/compcor
mgxd Apr 24, 2017
7ff5148
fix: remove temp slurm fix
mgxd Apr 24, 2017
af1fd73
fix: tcompcor output
mgxd Apr 24, 2017
c442383
fix: grab first mask for tcompcor
mgxd Apr 24, 2017
2b24f31
enh: multiple masks and options to compcor
mgxd Apr 27, 2017
2eb0176
fix: remove mask_file check
mgxd Apr 27, 2017
138afc9
Merge branch 'master' of git://github.com/nipy/nipype into fix/compcor
mgxd Apr 27, 2017
a9cc329
fix: fixed failing tests
mgxd Apr 27, 2017
164a78e
fix: deprecation cycle and basic testing
mgxd Apr 28, 2017
668cf8a
fix: sidestep Undefined issue
mgxd Apr 28, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
236 changes: 172 additions & 64 deletions nipype/algorithms/confounds.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@
from ..external.due import BibTeX
from ..interfaces.base import (traits, TraitedSpec, BaseInterface,
BaseInterfaceInputSpec, File, isdefined,
InputMultiPath)
InputMultiPath, OutputMultiPath)
from ..utils import NUMPY_MMAP
from ..utils.misc import normalize_mc_params

Expand Down Expand Up @@ -301,23 +301,33 @@ def _list_outputs(self):

class CompCorInputSpec(BaseInterfaceInputSpec):
realigned_file = File(exists=True, mandatory=True,
desc='already realigned brain image (4D)')
mask_file = File(exists=True, desc='mask file that determines ROI (3D)')
components_file = File('components_file.txt', exists=False,
usedefault=True,
desc='filename to store physiological components')
desc='already realigned brain image (4D)')
mask_file = InputMultiPath(File(exists=True, deprecated='0.13',
new_name='mask_files',
desc='One or more mask files that determines ROI (3D)'))
mask_files = InputMultiPath(File(exists=True,
desc='One or more mask files that determines ROI (3D)'))
merge_method = traits.Enum('union', 'intersect', 'none', xor=['mask_index'],
requires=['mask_files'],
desc='Merge method if multiple masks are present - `union` aggregates '
'all masks, `intersect` computes the truth value of all masks, `none` '
'performs CompCor on each mask individually')
mask_index = traits.Range(0, xor=['merge_method'], requires=['mask_files'],
desc='Position of mask in `mask_files` to use - first is the default')
components_file = File('components_file.txt', exists=False, usedefault=True,
desc='filename to store physiological components')
num_components = traits.Int(6, usedefault=True) # 6 for BOLD, 4 for ASL
use_regress_poly = traits.Bool(True, usedefault=True,
desc='use polynomial regression'
'pre-component extraction')
desc='use polynomial regression pre-component extraction')
regress_poly_degree = traits.Range(low=1, default=1, usedefault=True,
desc='the degree polynomial to use')
header = traits.Str(desc='the desired header for the output tsv file (one column).'
'If undefined, will default to "CompCor"')
desc='the degree polynomial to use')
header = traits.Str(
desc='the desired header for the output tsv file (one column).'
'If undefined, will default to "CompCor"')

class CompCorOutputSpec(TraitedSpec):
components_file = File(exists=True,
desc='text file containing the noise components')
desc='text file containing the noise components')

class CompCor(BaseInterface):
'''
Expand All @@ -328,7 +338,7 @@ class CompCor(BaseInterface):

>>> ccinterface = CompCor()
>>> ccinterface.inputs.realigned_file = 'functional.nii'
>>> ccinterface.inputs.mask_file = 'mask.nii'
>>> ccinterface.inputs.mask_files = 'mask.nii'
>>> ccinterface.inputs.num_components = 1
>>> ccinterface.inputs.use_regress_poly = True
>>> ccinterface.inputs.regress_poly_degree = 2
Expand All @@ -349,15 +359,53 @@ class CompCor(BaseInterface):
'tags': ['method', 'implementation']
}]

def _run_interface(self, runtime):
imgseries = nb.load(self.inputs.realigned_file, mmap=NUMPY_MMAP).get_data()
mask = nb.load(self.inputs.mask_file, mmap=NUMPY_MMAP).get_data()

if imgseries.shape[:3] != mask.shape:
raise ValueError('Inputs for CompCor, func {} and mask {}, do not have matching '
'spatial dimensions ({} and {}, respectively)'
.format(self.inputs.realigned_file, self.inputs.mask_file,
imgseries.shape[:3], mask.shape))
def _run_interface(self, runtime, tCompCor_mask=False):

imgseries = nb.load(self.inputs.realigned_file,
mmap=NUMPY_MMAP).get_data()
components = None

if isdefined(self.inputs.mask_files) or isdefined(self.inputs.mask_file):
if (not isdefined(self.inputs.mask_index) and
not isdefined(self.inputs.merge_method)):
self.inputs.mask_index = 0

if isdefined(self.inputs.mask_index):
if self.inputs.mask_index < len(self.inputs.mask_files):
self.inputs.mask_files = [
self.inputs.mask_files[self.inputs.mask_index]]
else:
self.inputs.mask_files = self.inputs.mask_files[0]
if not tCompCor_mask:
RuntimeWarning('Mask index exceeded number of masks, using '
'mask {} instead'.format(self.inputs.mask_files[0]))

for mask_file in self.inputs.mask_files:
mask = nb.load(mask_file, mmap=NUMPY_MMAP).get_data()

if imgseries.shape[:3] != mask.shape:
raise ValueError('Inputs for CompCor, func {} and mask {}, '
'do not have matching spatial dimensions '
'({} and {}, respectively)'.format(
self.inputs.realigned_file, mask_file,
imgseries.shape[:3], mask.shape))

if (isdefined(self.inputs.merge_method) and
self.inputs.merge_method != 'none' and
len(self.inputs.mask_files) > 1):
if mask_file == self.inputs.mask_files[0]:
new_mask = mask
continue
else:
if self.inputs.merge_method == 'union':
new_mask = np.logical_or(new_mask, mask).astype(int)
elif self.inputs.merge_method == 'intersect':
new_mask = np.logical_and(new_mask, mask).astype(int)

if mask_file != self.inputs.mask_files[-1]:
continue
else: # merge complete
mask = new_mask

voxel_timecourses = imgseries[mask > 0]
# Zero-out any bad values
Expand All @@ -366,7 +414,8 @@ def _run_interface(self, runtime):
# from paper:
# "The constant and linear trends of the columns in the matrix M were
# removed [prior to ...]"
degree = self.inputs.regress_poly_degree if self.inputs.use_regress_poly else 0
degree = (self.inputs.regress_poly_degree if
self.inputs.use_regress_poly else 0)
voxel_timecourses = regress_poly(degree, voxel_timecourses)

# "Voxel time series from the noise ROI (either anatomical or tSTD) were
Expand All @@ -380,9 +429,13 @@ def _run_interface(self, runtime):
# "The covariance matrix C = MMT was constructed and decomposed into its
# principal components using a singular value decomposition."
u, _, _ = linalg.svd(M, full_matrices=False)
components = u[:, :self.inputs.num_components]
components_file = os.path.join(os.getcwd(), self.inputs.components_file)
if components is None:
components = u[:, :self.inputs.num_components]
else:
components = np.hstack((components,
u[:, :self.inputs.num_components]))

components_file = os.path.join(os.getcwd(), self.inputs.components_file)
self._set_header()
np.savetxt(components_file, components, fmt=b"%.10f", delimiter='\t',
header=self._make_headers(components.shape[1]), comments='')
Expand All @@ -401,7 +454,8 @@ def _compute_tSTD(self, M, x, axis=0):
return stdM

def _set_header(self, header='CompCor'):
self.inputs.header = self.inputs.header if isdefined(self.inputs.header) else header
self.inputs.header = (self.inputs.header if isdefined(self.inputs.header)
else header)

def _make_headers(self, num_col):
headers = []
Expand Down Expand Up @@ -434,7 +488,8 @@ class TCompCorInputSpec(CompCorInputSpec):

class TCompCorOutputSpec(CompCorInputSpec):
# and all the fields in CompCorInputSpec
high_variance_mask = File(exists=True, desc="voxels excedding the variance threshold")
high_variance_masks = OutputMultiPath(File(exists=True,
desc="voxels excedding the variance threshold"))

class TCompCor(CompCor):
'''
Expand All @@ -445,7 +500,7 @@ class TCompCor(CompCor):

>>> ccinterface = TCompCor()
>>> ccinterface.inputs.realigned_file = 'functional.nii'
>>> ccinterface.inputs.mask_file = 'mask.nii'
>>> ccinterface.inputs.mask_files = 'mask.nii'
>>> ccinterface.inputs.num_components = 1
>>> ccinterface.inputs.use_regress_poly = True
>>> ccinterface.inputs.regress_poly_degree = 2
Expand All @@ -456,52 +511,105 @@ class TCompCor(CompCor):
output_spec = TCompCorOutputSpec

def _run_interface(self, runtime):
imgseries = nb.load(self.inputs.realigned_file, mmap=NUMPY_MMAP).get_data()

_out_masks = []
img = nb.load(self.inputs.realigned_file, mmap=NUMPY_MMAP)
imgseries = img.get_data()
aff = img.affine


if imgseries.ndim != 4:
raise ValueError('tCompCor expected a 4-D nifti file. Input {} has {} dimensions '
'(shape {})'
.format(self.inputs.realigned_file, imgseries.ndim, imgseries.shape))

if isdefined(self.inputs.mask_file):
in_mask_data = nb.load(self.inputs.mask_file, mmap=NUMPY_MMAP).get_data()
imgseries = imgseries[in_mask_data != 0, :]

# From the paper:
# "For each voxel time series, the temporal standard deviation is
# defined as the standard deviation of the time series after the removal
# of low-frequency nuisance terms (e.g., linear and quadratic drift)."
imgseries = regress_poly(2, imgseries)

# "To construct the tSTD noise ROI, we sorted the voxels by their
# temporal standard deviation ..."
tSTD = self._compute_tSTD(imgseries, 0, axis=-1)

# use percentile_threshold to pick voxels
threshold_std = np.percentile(tSTD, 100. * (1. - self.inputs.percentile_threshold))
mask = tSTD >= threshold_std

if isdefined(self.inputs.mask_file):
mask_data = np.zeros_like(in_mask_data)
mask_data[in_mask_data != 0] = mask
raise ValueError('tCompCor expected a 4-D nifti file. Input {} has '
'{} dimensions (shape {})'.format(
self.inputs.realigned_file, imgseries.ndim,
imgseries.shape))

if isdefined(self.inputs.mask_files):
if (not isdefined(self.inputs.mask_index) and
not isdefined(self.inputs.merge_method)):
self.inputs.mask_index = 0
if isdefined(self.inputs.mask_index):
if self.inputs.mask_index < len(self.inputs.mask_files):
self.inputs.mask_files = [
self.inputs.mask_files[self.inputs.mask_index]]
else:
RuntimeWarning('Mask index exceeded number of masks, using '
'mask {} instead'.format(self.inputs.mask_files[0]))
self.inputs.mask_files = self.inputs.mask_files[0]

for i, mask_file in enumerate(self.inputs.mask_files, 1):
in_mask = nb.load(mask_file, mmap=NUMPY_MMAP).get_data()
if (isdefined(self.inputs.merge_method) and
self.inputs.merge_method != 'none' and
len(self.inputs.mask_files) > 1):
if mask_file == self.inputs.mask_files[0]:
new_mask = in_mask
continue
else:
if self.inputs.merge_method == 'union':
new_mask = np.logical_or(new_mask,
in_mask).astype(int)
elif self.inputs.merge_method == 'intersect':
new_mask = np.logical_and(new_mask,
in_mask).astype(int)
if mask_file != self.inputs.mask_files[-1]:
continue
else: # merge complete
in_mask = new_mask

imgseries = imgseries[in_mask != 0, :]

# From the paper:
# "For each voxel time series, the temporal standard deviation is
# defined as the standard deviation of the time series after the removal
# of low-frequency nuisance terms (e.g., linear and quadratic drift)."
imgseries = regress_poly(2, imgseries)

# "To construct the tSTD noise ROI, we sorted the voxels by their
# temporal standard deviation ..."
tSTD = self._compute_tSTD(imgseries, 0, axis=-1)

# use percentile_threshold to pick voxels
threshold_std = np.percentile(tSTD, 100. *
(1. - self.inputs.percentile_threshold))
mask = tSTD >= threshold_std

mask_data = np.zeros_like(in_mask)
mask_data[in_mask != 0] = mask
# save mask
if self.inputs.merge_method == 'none':
mask_file = os.path.abspath('mask{}.nii'.format(i))
else:
mask_file = os.path.abspath('mask.nii')
nb.Nifti1Image(mask_data, aff).to_filename(mask_file)
IFLOG.debug('tCompcor computed and saved mask of shape {} to '
'mask_file {}'.format(mask.shape, mask_file))
_out_masks.append(mask_file)
self._set_header('tCompCor')

else:
imgseries = regress_poly(2, imgseries)
tSTD = self._compute_tSTD(imgseries, 0, axis=-1)
threshold_std = np.percentile(tSTD, 100. *
(1. - self.inputs.percentile_threshold))
mask = tSTD >= threshold_std
mask_data = mask.astype(int)

# save mask
mask_file = os.path.abspath('mask.nii')
nb.Nifti1Image(mask_data,
nb.load(self.inputs.realigned_file).affine).to_filename(mask_file)
IFLOG.debug('tCompcor computed and saved mask of shape {} to mask_file {}'
.format(mask.shape, mask_file))
self.inputs.mask_file = mask_file
self._set_header('tCompCor')
# save mask
mask_file = os.path.abspath('mask.nii')
nb.Nifti1Image(mask_data, aff).to_filename(mask_file)
IFLOG.debug('tCompcor computed and saved mask of shape {} to '
'mask_file {}'.format(mask.shape, mask_file))
_out_masks.append(mask_file)
self._set_header('tCompCor')

super(TCompCor, self)._run_interface(runtime)
self.inputs.mask_files = _out_masks
super(TCompCor, self)._run_interface(runtime, tCompCor_mask=True)
return runtime

def _list_outputs(self):
outputs = super(TCompCor, self)._list_outputs()
outputs['high_variance_mask'] = self.inputs.mask_file
outputs['high_variance_masks'] = self.inputs.mask_files
return outputs

class TSNRInputSpec(BaseInterfaceInputSpec):
Expand Down
Loading