enh: multiple masks and options to compcor

Author: mgxd

nipype/algorithms/confounds.py

@@ -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
 
@@ -301,23 +301,29 @@ def _list_outputs(self):
 
 class CompCorInputSpec(BaseInterfaceInputSpec):
     realigned_file = File(exists=True, mandatory=True,
-                          desc='already realigned brain image (4D)')
-    mask_file = InputMultiPath(File(exists=True, desc='mask file(s) 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_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'],
+        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, value=0, xor=['merge_method'], usedefault=True,
+        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):
     '''
@@ -328,7 +334,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
@@ -350,17 +356,46 @@ class CompCor(BaseInterface):
                    }]
 
     def _run_interface(self, runtime):
-        imgseries = nb.load(self.inputs.realigned_file, mmap=NUMPY_MMAP).get_data()
+        imgseries = nb.load(self.inputs.realigned_file,
+                            mmap=NUMPY_MMAP).get_data()
         components = None
 
-        for mask_file in self.inputs.mask_file:
-            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.mask_index) and 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]
+
+        if isdefined(self.inputs.mask_files):
+            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.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
@@ -369,7 +404,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
@@ -386,18 +422,18 @@ def _run_interface(self, runtime):
             if components is None:
                 components = u[:, :self.inputs.num_components]
             else:
-                components = np.hstack((components, u[:, :self.inputs.num_components]))
+                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='')
         return runtime
 
     def _list_outputs(self):
         outputs = self._outputs().get()
-        outputs['components_file'] = os.path.abspath(self.inputs.components_file)
+        outputs['components_file'] = self.inputs.components_file
         return outputs
 
     def _compute_tSTD(self, M, x, axis=0):
@@ -408,7 +444,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 = []
@@ -441,7 +478,8 @@ class TCompCorInputSpec(CompCorInputSpec):
 
 class TCompCorOutputSpec(CompCorInputSpec):
     # and all the fields in CompCorInputSpec
-    high_variance_mask = InputMultiPath(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):
     '''
@@ -452,7 +490,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
@@ -461,54 +499,99 @@ class TCompCor(CompCor):
 
     input_spec = TCompCorInputSpec
     output_spec = TCompCorOutputSpec
+    _out_masks = []
 
     def _run_interface(self, runtime):
-        imgseries = nb.load(self.inputs.realigned_file, mmap=NUMPY_MMAP).get_data()
+        imgseries = nb.load(self.inputs.realigned_file,
+                            mmap=NUMPY_MMAP).get_data()
 
         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[0], 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 isdefined(self.inputs.mask_index) and 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.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
+                mask_file = os.path.abspath('mask{}.nii'.format(i))
+                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._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,
+                  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._out_masks.append(mask_file)
+            self._set_header('tCompCor')
 
+        self.inputs.mask_files = self._out_masks
         super(TCompCor, self)._run_interface(runtime)
         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):

nipype/algorithms/tests/test_compcor.py

@@ -28,7 +28,7 @@ def setup_class(self, tmpdir):
         mask = np.ones(self.fake_data.shape[:3])
         mask[0, 0, 0] = 0
         mask[0, 0, 1] = 0
-        self.mask_file = utils.save_toy_nii(mask, self.filenames['masknii'])
+        self.mask_files = utils.save_toy_nii(mask, self.filenames['masknii'])
 
     def test_compcor(self):
         expected_components = [['-0.1989607212', '-0.5753813646'],
@@ -37,10 +37,10 @@ def test_compcor(self):
                                ['0.4206466244', '-0.3361270124'],
                                ['-0.1246655485', '-0.1235705610']]
 
-        self.run_cc(CompCor(realigned_file=self.realigned_file, mask_file=self.mask_file),
+        self.run_cc(CompCor(realigned_file=self.realigned_file, mask_files=self.mask_files),
                     expected_components)
 
-        self.run_cc(ACompCor(realigned_file=self.realigned_file, mask_file=self.mask_file,
+        self.run_cc(ACompCor(realigned_file=self.realigned_file, mask_files=self.mask_files,
                              components_file='acc_components_file'),
                     expected_components, 'aCompCor')
 
@@ -62,7 +62,7 @@ def test_tcompcor_no_percentile(self):
         assert num_nonmasked_voxels == 1
 
     def test_compcor_no_regress_poly(self):
-        self.run_cc(CompCor(realigned_file=self.realigned_file, mask_file=self.mask_file,
+        self.run_cc(CompCor(realigned_file=self.realigned_file, mask_files=self.mask_files,
                             use_regress_poly=False), [['0.4451946442', '-0.7683311482'],
                                                       ['-0.4285129505', '-0.0926034137'],
                                                       ['0.5721540256', '0.5608764842'],
@@ -77,12 +77,12 @@ def test_tcompcor_asymmetric_dim(self):
         assert nb.load('mask.nii').get_data().shape == asymmetric_shape[:3]
 
     def test_compcor_bad_input_shapes(self):
-        shape_less_than = (1, 2, 2, 5) # dim 0 is < dim 0 of self.mask_file (2)
-        shape_more_than = (3, 3, 3, 5) # dim 0 is > dim 0 of self.mask_file (2)
+        shape_less_than = (1, 2, 2, 5) # dim 0 is < dim 0 of self.mask_files (2)
+        shape_more_than = (3, 3, 3, 5) # dim 0 is > dim 0 of self.mask_files (2)
 
         for data_shape in (shape_less_than, shape_more_than):
             data_file = utils.save_toy_nii(np.zeros(data_shape), 'temp.nii')
-            interface = CompCor(realigned_file=data_file, mask_file=self.mask_file)
+            interface = CompCor(realigned_file=data_file, mask_files=self.mask_files)
             with pytest.raises(ValueError, message="Dimension mismatch"): interface.run()
 
     def test_tcompcor_bad_input_dim(self):