Merge pull request #1170 from effigies/rf/signal_extract_lowmem

RF: Reimplement SignalExtraction

Author: oesteban

fmriprep/interfaces/images.py

@@ -572,3 +572,77 @@ def nii_ones_like(in_file, value, dtype, newpath=None):
     nii.to_filename(out_file)
 
     return out_file
+
+
+class SignalExtractionInputSpec(BaseInterfaceInputSpec):
+    in_file = File(exists=True, mandatory=True, desc='4-D fMRI nii file')
+    label_files = InputMultiPath(
+        File(exists=True),
+        mandatory=True,
+        desc='a 3-D label image, with 0 denoting '
+        'background, or a list of 3-D probability '
+        'maps (one per label) or the equivalent 4D '
+        'file.')
+    class_labels = traits.List(
+        mandatory=True,
+        desc='Human-readable labels for each segment '
+        'in the label file, in order. The length of '
+        'class_labels must be equal to the number of '
+        'segments (background excluded). This list '
+        'corresponds to the class labels in label_file '
+        'in ascending order')
+    out_file = File(
+        'signals.tsv',
+        usedefault=True,
+        exists=False,
+        desc='The name of the file to output to. '
+        'signals.tsv by default')
+
+
+class SignalExtractionOutputSpec(TraitedSpec):
+    out_file = File(
+        exists=True,
+        desc='tsv file containing the computed '
+        'signals, with as many columns as there are labels and as '
+        'many rows as there are timepoints in in_file, plus a '
+        'header row with values from class_labels')
+
+
+class SignalExtraction(SimpleInterface):
+    """ Extract mean signals from a time series within a set of ROIs
+
+    This interface is intended to be a memory-efficient alternative to
+    nipype.interfaces.nilearn.SignalExtraction.
+    Not all features of nilearn.SignalExtraction are implemented at
+    this time.
+    """
+    input_spec = SignalExtractionInputSpec
+    output_spec = SignalExtractionOutputSpec
+
+    def _run_interface(self, runtime):
+        mask_imgs = [nb.load(fname) for fname in self.inputs.label_files]
+        if len(mask_imgs) == 1:
+            mask_imgs = nb.four_to_three(mask_imgs[0])
+
+        masks = [mask_img.get_data().astype(np.bool) for mask_img in mask_imgs]
+
+        n_masks = len(masks)
+
+        if n_masks != len(self.inputs.class_labels):
+            raise ValueError("Number of masks must match number of labels")
+
+        img = nb.load(self.inputs.in_file)
+
+        series = np.zeros((img.shape[3], n_masks))
+
+        data = img.get_data()
+        for j in range(n_masks):
+            series[:, j] = data[masks[j], :].mean(axis=0)
+
+        output = np.vstack((self.inputs.class_labels, series.astype(str)))
+        self._results['out_file'] = os.path.join(runtime.cwd,
+                                                 self.inputs.out_file)
+        np.savetxt(
+            self._results['out_file'], output, fmt=b'%s', delimiter='\t')
+
+        return runtime

fmriprep/interfaces/tests/__init__.py

@@ -0,0 +1,53 @@
+import time
+import numpy as np
+import nibabel as nb
+from nipype.interfaces import nilearn as nl
+from .. import images as im
+
+import pytest
+
+
+@pytest.mark.parametrize('nvols, nmasks, ext, factor', [
+    (500, 10, '.nii', 2),
+    (500, 10, '.nii.gz', 5),
+    (200, 3, '.nii', 1.1),
+    (200, 3, '.nii.gz', 2),
+    (200, 10, '.nii', 1.1),
+    (200, 10, '.nii.gz', 2),
+    ])
+def test_signal_extraction_equivalence(tmpdir, nvols, nmasks, ext, factor):
+    tmpdir.chdir()
+
+    vol_shape = (64, 64, 40)
+
+    img_fname = 'img' + ext
+    masks_fname = 'masks' + ext
+
+    random_data = np.random.random(size=vol_shape + (nvols,)) * 2000
+    random_mask_data = np.random.random(size=vol_shape + (nmasks,)) < 0.2
+
+    nb.Nifti1Image(random_data, np.eye(4)).to_filename(img_fname)
+    nb.Nifti1Image(random_mask_data.astype(np.uint8), np.eye(4)).to_filename(masks_fname)
+
+    se1 = nl.SignalExtraction(in_file=img_fname, label_files=masks_fname,
+                              class_labels=['a%d' % i for i in range(nmasks)],
+                              out_file='nlsignals.tsv')
+    se2 = im.SignalExtraction(in_file=img_fname, label_files=masks_fname,
+                              class_labels=['a%d' % i for i in range(nmasks)],
+                              out_file='imsignals.tsv')
+
+    tic = time.time()
+    se1.run()
+    toc = time.time()
+    se2.run()
+    toc2 = time.time()
+
+    tab1 = np.loadtxt('nlsignals.tsv', skiprows=1)
+    tab2 = np.loadtxt('imsignals.tsv', skiprows=1)
+
+    assert np.allclose(tab1, tab2)
+
+    t1 = toc - tic
+    t2 = toc2 - toc
+
+    assert t2 < t1 / factor