ENH CompCor

nipype/algorithms/compcor.py

@@ -0,0 +1,81 @@
+from ..interfaces.base import (BaseInterfaceInputSpec, TraitedSpec,
+                               BaseInterface, traits, File)
+import nibabel as nb
+import numpy as np
+from scipy import linalg, stats
+import os
+
+from nipype.pipeline.engine import Workflow
+
+class CompCoreInputSpec(BaseInterfaceInputSpec):
+    realigned_file = File(exists=True, mandatory=True, desc='already realigned brain image (4D)')
+    mask_file = File(exists=True, mandatory=True, desc='mask file that determines ROI (3D)')
+    num_components = traits.Int(6, usedefault=True) # 6 for BOLD, 4 for ASL
+    # additional_regressors??
+
+class CompCoreOutputSpec(TraitedSpec):
+    components_file = File(desc='text file containing the noise components', exists=True)
+
+class CompCore(BaseInterface):
+    '''
+    Interface with core CompCor computation, used in aCompCor and tCompCor
+
+    Example
+    -------
+
+    >>> ccinterface = CompCore()
+    >>> ccinterface.inputs.realigned_file = 'nipype/testing/data/functional.nii'
+    >>> ccinterface.inputs.mask_file = 'nipype/testing/data/mask.nii'
+    >>> ccinterface.inputs.num_components = 1
+    '''
+    input_spec = CompCoreInputSpec
+    output_spec = CompCoreOutputSpec
+
+    def _run_interface(self, runtime):
+        imgseries = nb.load(self.inputs.realigned_file).get_data()
+        mask = nb.load(self.inputs.mask_file).get_data()
+        voxel_timecourses = imgseries[mask > 0]
+        # Zero-out any bad values
+        voxel_timecourses[np.isnan(np.sum(voxel_timecourses, axis=1)), :] = 0
+
+        # from paper:
+        # "Voxel time series from the noise ROI (either anatomical or tSTD) were
+        # placed in a matrix M of size Nxm, with time along the row dimension
+        # and voxels along the column dimension."
+        # voxel_timecourses.shape == [nvoxels, time]
+        M = voxel_timecourses.T
+        numvols = M.shape[0]
+        numvoxels = M.shape[1]
+
+        # "The constant and linear trends of the columns in the matrix M were removed ..."
+        timesteps = range(numvols)
+        for voxel in range(numvoxels):
+            m, b, _, _, _ = stats.linregress(M[:, voxel], timesteps)
+            M[:, voxel] = M[:, voxel] - [m*t + b for t in timesteps]
+
+        # "... prior to column-wise variance normalization."
+        stdM = np.std(M, axis=0)
+        # set bad values to division identity
+        stdM[stdM == 0] = 1.
+        stdM[np.isnan(stdM)] = 1.
+        stdM[np.isinf(stdM)] = 1.
+        M = M / stdM
+
+        # "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(), "components_file.txt")
+        np.savetxt(components_file, components, fmt="%.10f")
+        return runtime
+
+    def _list_outputs(self):
+        outputs = self._outputs().get()
+        outputs['components_file'] = os.path.abspath("components_file.txt")
+        return outputs
+
+class aCompCor(Workflow):
+    pass
+
+class tCompCor(Workflow):
+    pass

nipype/algorithms/tests/test_compcor.py

@@ -0,0 +1,60 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+import nipype
+from nipype.testing import assert_equal, assert_true, assert_false, skipif
+from nipype.algorithms.compcor import CompCore
+
+import nibabel as nb
+import numpy as np
+import os
+
+functionalnii = 'func.nii'
+masknii = 'mask.nii'
+
+def test_compcore():
+    # setup
+    noise = np.fromfunction(fake_noise_fun, fake_data.shape)
+    realigned_file = make_toy(fake_data + noise, functionalnii)
+
+    mask = np.ones(fake_data.shape[:3])
+    mask[0,0,0] = 0
+    mask[0,0,1] = 0
+    mask_file = make_toy(mask, masknii)
+
+    # run
+    ccinterface = CompCore(realigned_file=realigned_file, mask_file=mask_file)
+    ccresult = ccinterface.run()
+
+    # asserts
+    components = ccinterface._list_outputs()['components_file']
+    assert_true(os.path.exists(components))
+    assert_true(os.path.getsize(components) > 0)
+    assert_equal(ccinterface.inputs.num_components, 6)
+
+    #  apply components_file to realigned_file, is it better? the same as before adding noise?
+
+    # remove temporary nifti files
+    os.remove(functionalnii)
+    os.remove(masknii)
+    os.remove(components)
+
+def make_toy(ndarray, filename):
+    toy = nb.Nifti1Image(ndarray, np.eye(4))
+    nb.nifti1.save(toy, filename)
+    return filename
+
+fake_data = np.array([[[[8, 5, 3, 8, 0],
+                        [6, 7, 4, 7, 1]],
+
+                       [[7, 9, 1, 6, 5],
+                        [0, 7, 4, 7, 7]]],
+
+
+                      [[[2, 4, 5, 7, 0],
+                        [1, 7, 0, 5, 4]],
+
+                       [[7, 3, 9, 0, 4],
+                        [9, 4, 1, 5, 0]]]])
+
+def fake_noise_fun(i, j, l, m):
+    return m*i + l - j