Merge pull request #1676 from shoshber/debugging

WIP logging and validations in CompCor, SignalExtraction, ApplyTopUp

Author: Shoshana Berleant

circle.yml

@@ -38,7 +38,7 @@ test:
     - docker run -v /etc/localtime:/etc/localtime:ro -e FSL_COURSE_DATA="/root/examples/nipype-fsl_course_data" -v ~/examples:/root/examples:ro -v ~/scratch:/scratch -w /root/src/nipype nipype/nipype_test:py35 /usr/bin/run_nosetests.sh py35 :
         timeout: 2600
     - docker run -v /etc/localtime:/etc/localtime:ro -e FSL_COURSE_DATA="/root/examples/nipype-fsl_course_data" -v ~/examples:/root/examples:ro -v ~/scratch:/scratch -w /root/src/nipype nipype/nipype_test:py27 /usr/bin/run_nosetests.sh py27 :
-        timeout: 2600
+        timeout: 5200
     - docker run -v /etc/localtime:/etc/localtime:ro -v ~/examples:/root/examples:ro -v ~/scratch:/scratch -w /scratch nipype/nipype_test:py35 /usr/bin/run_examples.sh test_spm Linear /root/examples/ workflow3d :
         timeout: 1600
     - docker run -v /etc/localtime:/etc/localtime:ro -v ~/examples:/root/examples:ro -v ~/scratch:/scratch -w /scratch nipype/nipype_test:py35 /usr/bin/run_examples.sh test_spm Linear /root/examples/ workflow4d :

nipype/algorithms/confounds.py

@@ -32,7 +32,7 @@
 class ComputeDVARSInputSpec(BaseInterfaceInputSpec):
     in_file = File(exists=True, mandatory=True, desc='functional data, after HMC')
     in_mask = File(exists=True, mandatory=True, desc='a brain mask')
-    remove_zerovariance = traits.Bool(False, usedefault=True,
+    remove_zerovariance = traits.Bool(True, usedefault=True,
                                       desc='remove voxels with zero variance')
     save_std = traits.Bool(True, usedefault=True,
                            desc='save standardized DVARS')
@@ -255,7 +255,7 @@ def _run_interface(self, runtime):
             'out_file': op.abspath(self.inputs.out_file),
             'fd_average': float(fd_res.mean())
         }
-        np.savetxt(self.inputs.out_file, fd_res)
+        np.savetxt(self.inputs.out_file, fd_res, header='framewise_displacement')
 
         if self.inputs.save_plot:
             tr = None
@@ -291,6 +291,8 @@ class CompCorInputSpec(BaseInterfaceInputSpec):
                                    '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"')
 
 class CompCorOutputSpec(TraitedSpec):
     components_file = File(exists=True,
@@ -329,6 +331,13 @@ class CompCor(BaseInterface):
     def _run_interface(self, runtime):
         imgseries = nb.load(self.inputs.realigned_file).get_data()
         mask = nb.load(self.inputs.mask_file).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))
+
         voxel_timecourses = imgseries[mask > 0]
         # Zero-out any bad values
         voxel_timecourses[np.isnan(np.sum(voxel_timecourses, axis=1)), :] = 0
@@ -352,7 +361,10 @@ def _run_interface(self, runtime):
         u, _, _ = linalg.svd(M, full_matrices=False)
         components = u[:, :self.inputs.num_components]
         components_file = os.path.join(os.getcwd(), self.inputs.components_file)
-        np.savetxt(components_file, components, fmt=b"%.10f")
+
+        self._set_header()
+        np.savetxt(components_file, components, fmt=b"%.10f", delimiter='\t',
+                   header=self._make_headers(components.shape[1]))
         return runtime
 
     def _list_outputs(self):
@@ -367,6 +379,26 @@ def _compute_tSTD(self, M, x):
         stdM[np.isnan(stdM)] = x
         return stdM
 
+    def _set_header(self, header='CompCor'):
+        self.inputs.header = self.inputs.header if isdefined(self.inputs.header) else header
+
+    def _make_headers(self, num_col):
+        headers = []
+        for i in range(num_col):
+            headers.append(self.inputs.header + str(i))
+        return '\t'.join(headers)
+
+
+class ACompCor(CompCor):
+    ''' Anatomical compcor; for input/output, see CompCor.
+    If the mask provided is an anatomical mask, CompCor == ACompCor '''
+
+    def __init__(self, *args, **kwargs):
+        ''' exactly the same as compcor except the header '''
+        super(ACompCor, self).__init__(*args, **kwargs)
+        self._set_header('aCompCor')
+
+
 class TCompCorInputSpec(CompCorInputSpec):
     # and all the fields in CompCorInputSpec
     percentile_threshold = traits.Range(low=0., high=1., value=.02,
@@ -401,6 +433,11 @@ class TCompCor(CompCor):
     def _run_interface(self, runtime):
         imgseries = nb.load(self.inputs.realigned_file).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))
+
         # 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
@@ -419,18 +456,19 @@ def _run_interface(self, runtime):
         threshold_index = int(num_voxels * (1. - self.inputs.percentile_threshold))
         threshold_std = sortSTD[threshold_index]
         mask = tSTD >= threshold_std
-        mask = mask.astype(int)
+        mask = mask.astype(int).T
 
         # save mask
-        mask_file = 'mask.nii'
+        mask_file = os.path.abspath('mask.nii')
         nb.nifti1.save(nb.Nifti1Image(mask, np.eye(4)), 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')
 
         super(TCompCor, self)._run_interface(runtime)
         return runtime
 
-ACompCor = CompCor
-
 class TSNRInputSpec(BaseInterfaceInputSpec):
     in_file = InputMultiPath(File(exists=True), mandatory=True,
                              desc='realigned 4D file or a list of 3D files')
@@ -512,6 +550,8 @@ def regress_poly(degree, data, remove_mean=True, axis=-1):
     If remove_mean is True (default), the data is demeaned (i.e. degree 0).
     If remove_mean is false, the data is not.
     '''
+    IFLOG.debug('Performing polynomial regression on data of shape ' + str(data.shape))
+
     datashape = data.shape
     timepoints = datashape[axis]
 
@@ -570,6 +610,7 @@ def compute_dvars(in_file, in_mask, remove_zerovariance=False):
     import numpy as np
     import nibabel as nb
     from nitime.algorithms import AR_est_YW
+    import warnings
 
     func = nb.load(in_file).get_data().astype(np.float32)
     mask = nb.load(in_mask).get_data().astype(np.uint8)
@@ -585,7 +626,7 @@ def compute_dvars(in_file, in_mask, remove_zerovariance=False):
 
     if remove_zerovariance:
         # Remove zero-variance voxels across time axis
-        mask = zero_variance(func, mask)
+        mask = zero_remove(func_sd, mask)
 
     idx = np.where(mask > 0)
     mfunc = func[idx[0], idx[1], idx[2], :]
@@ -609,31 +650,28 @@ def compute_dvars(in_file, in_mask, remove_zerovariance=False):
     # standardization
     dvars_stdz = dvars_nstd / diff_sd_mean
 
-    # voxelwise standardization
-    diff_vx_stdz = func_diff / np.array([diff_sdhat] * func_diff.shape[-1]).T
-    dvars_vx_stdz = diff_vx_stdz.std(axis=0, ddof=1)
+    with warnings.catch_warnings(): # catch, e.g., divide by zero errors
+        warnings.filterwarnings('error')
+
+        # voxelwise standardization
+        diff_vx_stdz = func_diff / np.array([diff_sdhat] * func_diff.shape[-1]).T            
+        dvars_vx_stdz = diff_vx_stdz.std(axis=0, ddof=1)
 
     return (dvars_stdz, dvars_nstd, dvars_vx_stdz)
 
-def zero_variance(func, mask):
+def zero_remove(data, mask):
     """
-    Mask out voxels with zero variance across t-axis
+    Modify inputted mask to also mask out zero values
 
-    :param numpy.ndarray func: input fMRI dataset, after motion correction
-    :param numpy.ndarray mask: 3D brain mask
-    :return: the 3D mask of voxels with nonzero variance across :math:`t`.
+    :param numpy.ndarray data: e.g. voxelwise stddev of fMRI dataset, after motion correction
+    :param numpy.ndarray mask: brain mask (same dimensions as data)
+    :return: the mask with any additional zero voxels removed (same dimensions as inputs)
     :rtype: numpy.ndarray
 
     """
-    idx = np.where(mask > 0)
-    func = func[idx[0], idx[1], idx[2], :]
-    tvariance = func.var(axis=1)
-    tv_mask = np.zeros_like(tvariance, dtype=np.uint8)
-    tv_mask[tvariance > 0] = 1
-
-    newmask = np.zeros_like(mask, dtype=np.uint8)
-    newmask[idx] = tv_mask
-    return newmask
+    new_mask = mask.copy()
+    new_mask[data == 0] = 0
+    return new_mask
 
 def plot_confound(tseries, figsize, name, units=None,
                   series_tr=None, normalize=False):

nipype/algorithms/tests/test_compcor.py

@@ -8,7 +8,7 @@
 import nibabel as nb
 import numpy as np
 
-from ...testing import assert_equal, assert_true, utils
+from ...testing import assert_equal, assert_true, utils, assert_in
 from ..confounds import CompCor, TCompCor, ACompCor
 
 class TestCompCor(unittest.TestCase):
@@ -38,25 +38,20 @@ def test_compcor(self):
                                ['0.4206466244', '-0.3361270124'],
                                ['-0.1246655485', '-0.1235705610']]
 
-        ccresult = self.run_cc(CompCor(realigned_file=self.realigned_file,
-                                       mask_file=self.mask_file),
-                               expected_components)
+        self.run_cc(CompCor(realigned_file=self.realigned_file, mask_file=self.mask_file),
+                    expected_components)
 
-        accresult = self.run_cc(ACompCor(realigned_file=self.realigned_file,
-                                         mask_file=self.mask_file,
-                                         components_file='acc_components_file'),
-                                expected_components)
-
-        assert_equal(os.path.getsize(ccresult.outputs.components_file),
-                     os.path.getsize(accresult.outputs.components_file))
+        self.run_cc(ACompCor(realigned_file=self.realigned_file, mask_file=self.mask_file,
+                             components_file='acc_components_file'),
+                    expected_components, 'aCompCor')
 
     def test_tcompcor(self):
         ccinterface = TCompCor(realigned_file=self.realigned_file, percentile_threshold=0.75)
         self.run_cc(ccinterface, [['-0.1114536190', '-0.4632908609'],
                                   ['0.4566907310', '0.6983205193'],
                                   ['-0.7132557407', '0.1340170559'],
                                   ['0.5022537643', '-0.5098322262'],
-                                  ['-0.1342351356', '0.1407855119']])
+                                  ['-0.1342351356', '0.1407855119']], 'tCompCor')
 
     def test_tcompcor_no_percentile(self):
         ccinterface = TCompCor(realigned_file=self.realigned_file)
@@ -74,7 +69,29 @@ def test_compcor_no_regress_poly(self):
                                                       ['-0.5367548139', '0.0059943226'],
                                                       ['-0.0520809054', '0.2940637551']])
 
-    def run_cc(self, ccinterface, expected_components):
+    def test_tcompcor_asymmetric_dim(self):
+        asymmetric_shape = (2, 3, 4, 5)
+        asymmetric_data = utils.save_toy_nii(np.zeros(asymmetric_shape), 'asymmetric.nii')
+
+        TCompCor(realigned_file=asymmetric_data).run()
+        self.assertEqual(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)
+
+        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)
+            self.assertRaisesRegexp(ValueError, "dimensions", interface.run)
+
+    def test_tcompcor_bad_input_dim(self):
+        bad_dims = (2, 2, 2)
+        data_file = utils.save_toy_nii(np.zeros(bad_dims), 'temp.nii')
+        interface = TCompCor(realigned_file=data_file)
+        self.assertRaisesRegexp(ValueError, '4-D', interface.run)
+
+    def run_cc(self, ccinterface, expected_components, expected_header='CompCor'):
         # run
         ccresult = ccinterface.run()
 
@@ -86,12 +103,21 @@ def run_cc(self, ccinterface, expected_components):
         assert_equal(ccinterface.inputs.num_components, 6)
 
         with open(ccresult.outputs.components_file, 'r') as components_file:
-            components_data = [line.split() for line in components_file]
-            num_got_components = len(components_data)
-            assert_true(num_got_components == ccinterface.inputs.num_components
-                        or num_got_components == self.fake_data.shape[3])
-            first_two = [row[:2] for row in components_data]
-            assert_equal(first_two, expected_components)
+            expected_n_components = min(ccinterface.inputs.num_components, self.fake_data.shape[3])
+
+            components_data = [line.split('\t') for line in components_file]
+
+            header = components_data.pop(0) # the first item will be '#', we can throw it out
+            expected_header = [expected_header + str(i) for i in range(expected_n_components)]
+            for i, heading in enumerate(header):
+                assert_in(expected_header[i], heading)
+
+            num_got_timepoints = len(components_data)
+            assert_equal(num_got_timepoints, self.fake_data.shape[3])
+            for index, timepoint in enumerate(components_data):
+                assert_true(len(timepoint) == ccinterface.inputs.num_components
+                            or len(timepoint) == self.fake_data.shape[3])
+                assert_equal(timepoint[:2], expected_components[index])
         return ccresult
 
     def tearDown(self):

nipype/algorithms/tests/test_confounds.py

@@ -4,7 +4,9 @@
 from tempfile import mkdtemp
 from shutil import rmtree
 
-from nipype.testing import (assert_equal, example_data, skipif, assert_true)
+from io import open
+
+from nipype.testing import (assert_equal, example_data, skipif, assert_true, assert_in)
 from nipype.algorithms.confounds import FramewiseDisplacement, ComputeDVARS
 import numpy as np
 
@@ -24,8 +26,14 @@ def test_fd():
                                           out_file=tempdir + '/fd.txt')
     res = fdisplacement.run()
 
+    with open(res.outputs.out_file) as all_lines:
+        for line in all_lines:
+            yield assert_in, 'framewise_displacement', line
+            break
+
     yield assert_true, np.allclose(ground_truth, np.loadtxt(res.outputs.out_file), atol=.16)
     yield assert_true, np.abs(ground_truth.mean() - res.outputs.fd_average) < 1e-2
+
     rmtree(tempdir)
 
 @skipif(nonitime)
@@ -35,8 +43,14 @@ def test_dvars():
     dvars = ComputeDVARS(in_file=example_data('ds003_sub-01_mc.nii.gz'),
                          in_mask=example_data('ds003_sub-01_mc_brainmask.nii.gz'),
                          save_all=True)
+
+    origdir = os.getcwd()
     os.chdir(tempdir)
+
     res = dvars.run()
 
     dv1 = np.loadtxt(res.outputs.out_std)
     yield assert_equal, (np.abs(dv1 - ground_truth).sum()/ len(dv1)) < 0.05, True
+
+    os.chdir(origdir)
+    rmtree(tempdir)

nipype/interfaces/ants/resampling.py

@@ -246,9 +246,10 @@ class ApplyTransformsInputSpec(ANTSCommandInputSpec):
     interpolation_parameters = traits.Either(traits.Tuple(traits.Int()),  # BSpline (order)
                                              traits.Tuple(traits.Float(),  # Gaussian/MultiLabel (sigma, alpha)
                                                           traits.Float())
-                                             )
-    transforms = InputMultiPath(
-        File(exists=True), argstr='%s', mandatory=True, desc='transform files: will be applied in reverse order. For example, the last specified transform will be applied first')
+                                            )
+    transforms = InputMultiPath(File(exists=True), argstr='%s', mandatory=True,
+                                desc='transform files: will be applied in reverse order. For '
+                                'example, the last specified transform will be applied first.')
     invert_transform_flags = InputMultiPath(traits.Bool())
     default_value = traits.Float(0.0, argstr='--default-value %g', usedefault=True)
     print_out_composite_warp_file = traits.Bool(False, requires=["output_image"],
@@ -296,8 +297,6 @@ class ApplyTransforms(ANTSCommand):
     'antsApplyTransforms --default-value 0 --dimensionality 3 --input moving1.nii --interpolation BSpline[ 5 ] \
 --output deformed_moving1.nii --reference-image fixed1.nii --transform [ ants_Warp.nii.gz, 0 ] \
 --transform [ trans.mat, 0 ]'
-
-
     """
     _cmd = 'antsApplyTransforms'
     input_spec = ApplyTransformsInputSpec

nipype/interfaces/fsl/epi.py

@@ -307,7 +307,7 @@ def _overload_extension(self, value, name=None):
 
 class ApplyTOPUPInputSpec(FSLCommandInputSpec):
     in_files = InputMultiPath(File(exists=True), mandatory=True,
-                              desc='name of 4D file with images',
+                              desc='name of file with images',
                               argstr='--imain=%s', sep=',')
     encoding_file = File(exists=True, mandatory=True,
                          desc='name of text file with PE directions/times',