diff --git a/nipype/interfaces/nilearn.py b/nipype/interfaces/nilearn.py new file mode 100644 index 0000000000..bc662826d0 --- /dev/null +++ b/nipype/interfaces/nilearn.py @@ -0,0 +1,141 @@ +# -*- coding: utf-8 -*- +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +''' +Algorithms to compute statistics on :abbr:`fMRI (functional MRI)` + + Change directory to provide relative paths for doctests + >>> import os + >>> filepath = os.path.dirname(os.path.realpath(__file__)) + >>> datadir = os.path.realpath(os.path.join(filepath, '../testing/data')) + >>> os.chdir(datadir) + +''' +from __future__ import (print_function, division, unicode_literals, + absolute_import) + +import numpy as np +import nibabel as nb + +from .. import logging +from ..interfaces.base import (traits, TraitedSpec, BaseInterface, + BaseInterfaceInputSpec, File, InputMultiPath) +IFLOG = logging.getLogger('interface') + +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, + mandatory=False, desc='The name of the file to output to. ' + 'signals.tsv by default') + incl_shared_variance = traits.Bool(True, usedefault=True, mandatory=False, desc='By default ' + '(True), returns simple time series calculated from each ' + 'region independently (e.g., for noise regression). If ' + 'False, returns unique signals for each region, discarding ' + 'shared variance (e.g., for connectivity. Only has effect ' + 'with 4D probability maps.') + include_global = traits.Bool(False, usedefault=True, mandatory=False, + desc='If True, include an extra column ' + 'labeled "global", with values calculated from the entire brain ' + '(instead of just regions).') + detrend = traits.Bool(False, usedefault=True, mandatory=False, + desc='If True, perform detrending using nilearn.') + +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(BaseInterface): + ''' + Extracts signals over tissue classes or brain regions + + >>> seinterface = SignalExtraction() + >>> seinterface.inputs.in_file = 'functional.nii' + >>> seinterface.inputs.label_files = 'segmentation0.nii.gz' + >>> seinterface.inputs.out_file = 'means.tsv' + >>> segments = ['CSF', 'gray', 'white'] + >>> seinterface.inputs.class_labels = segments + >>> seinterface.inputs.detrend = True + >>> seinterface.inputs.include_global = True + ''' + input_spec = SignalExtractionInputSpec + output_spec = SignalExtractionOutputSpec + + def _run_interface(self, runtime): + maskers = self._process_inputs() + + signals = [] + for masker in maskers: + signals.append(masker.fit_transform(self.inputs.in_file)) + region_signals = np.hstack(signals) + + output = np.vstack((self.inputs.class_labels, region_signals.astype(str))) + + # save output + np.savetxt(self.inputs.out_file, output, fmt=b'%s', delimiter='\t') + return runtime + + def _process_inputs(self): + ''' validate and process inputs into useful form. + Returns a list of nilearn maskers and the list of corresponding label names.''' + import nilearn.input_data as nl + import nilearn.image as nli + + label_data = nli.concat_imgs(self.inputs.label_files) + maskers = [] + + # determine form of label files, choose appropriate nilearn masker + if np.amax(label_data.get_data()) > 1: # 3d label file + n_labels = np.amax(label_data.get_data()) + maskers.append(nl.NiftiLabelsMasker(label_data)) + else: # 4d labels + n_labels = label_data.get_data().shape[3] + if self.inputs.incl_shared_variance: # 4d labels, independent computation + for img in nli.iter_img(label_data): + maskers.append(nl.NiftiMapsMasker(self._4d(img.get_data(), img.affine))) + else: # 4d labels, one computation fitting all + maskers.append(nl.NiftiMapsMasker(label_data)) + + # check label list size + if len(self.inputs.class_labels) != n_labels: + raise ValueError('The length of class_labels {} does not ' + 'match the number of regions {} found in ' + 'label_files {}'.format(self.inputs.class_labels, + n_labels, + self.inputs.label_files)) + + if self.inputs.include_global: + global_label_data = label_data.get_data().sum(axis=3) # sum across all regions + global_label_data = np.rint(global_label_data).astype(int).clip(0, 1) # binarize + global_label_data = self._4d(global_label_data, label_data.affine) + global_masker = nl.NiftiLabelsMasker(global_label_data, detrend=self.inputs.detrend) + maskers.insert(0, global_masker) + self.inputs.class_labels.insert(0, 'global') + + for masker in maskers: + masker.set_params(detrend=self.inputs.detrend) + + return maskers + + def _4d(self, array, affine): + ''' takes a 3-dimensional numpy array and an affine, + returns the equivalent 4th dimensional nifti file ''' + return nb.Nifti1Image(array[:, :, :, np.newaxis], affine) + + def _list_outputs(self): + outputs = self._outputs().get() + outputs['out_file'] = self.inputs.out_file + return outputs diff --git a/nipype/interfaces/tests/test_nilearn.py b/nipype/interfaces/tests/test_nilearn.py new file mode 100644 index 0000000000..074e71fee9 --- /dev/null +++ b/nipype/interfaces/tests/test_nilearn.py @@ -0,0 +1,199 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +import unittest +import os +import tempfile +import shutil + +import numpy as np + +from ...testing import (assert_equal, utils, assert_almost_equal, raises, + skipif) +from .. import nilearn as iface + +no_nilearn = True +try: + __import__('nilearn') + no_nilearn = False +except ImportError: + pass + +class TestSignalExtraction(unittest.TestCase): + + filenames = { + 'in_file': 'fmri.nii', + 'label_files': 'labels.nii', + '4d_label_file': '4dlabels.nii', + 'out_file': 'signals.tsv' + } + labels = ['csf', 'gray', 'white'] + global_labels = ['global'] + labels + + def setUp(self): + self.orig_dir = os.getcwd() + self.temp_dir = tempfile.mkdtemp() + os.chdir(self.temp_dir) + + utils.save_toy_nii(self.fake_fmri_data, self.filenames['in_file']) + utils.save_toy_nii(self.fake_label_data, self.filenames['label_files']) + + @skipif(no_nilearn) + def test_signal_extract_no_shared(self): + # run + iface.SignalExtraction(in_file=self.filenames['in_file'], + label_files=self.filenames['label_files'], + class_labels=self.labels, + incl_shared_variance=False).run() + # assert + self.assert_expected_output(self.labels, self.base_wanted) + + + @skipif(no_nilearn) + @raises(ValueError) + def test_signal_extr_bad_label_list(self): + # run + iface.SignalExtraction(in_file=self.filenames['in_file'], + label_files=self.filenames['label_files'], + class_labels=['bad'], + incl_shared_variance=False).run() + + @skipif(no_nilearn) + def test_signal_extr_equiv_4d_no_shared(self): + self._test_4d_label(self.base_wanted, self.fake_equiv_4d_label_data, + incl_shared_variance=False) + + @skipif(no_nilearn) + def test_signal_extr_4d_no_shared(self): + # set up & run & assert + self._test_4d_label(self.fourd_wanted, self.fake_4d_label_data, incl_shared_variance=False) + + @skipif(no_nilearn) + def test_signal_extr_global_no_shared(self): + # set up + wanted_global = [[-4./6], [-1./6], [3./6], [-1./6], [-7./6]] + for i, vals in enumerate(self.base_wanted): + wanted_global[i].extend(vals) + + # run + iface.SignalExtraction(in_file=self.filenames['in_file'], + label_files=self.filenames['label_files'], + class_labels=self.labels, + include_global=True, + incl_shared_variance=False).run() + + # assert + self.assert_expected_output(self.global_labels, wanted_global) + + @skipif(no_nilearn) + def test_signal_extr_4d_global_no_shared(self): + # set up + wanted_global = [[3./8], [-3./8], [1./8], [-7./8], [-9./8]] + for i, vals in enumerate(self.fourd_wanted): + wanted_global[i].extend(vals) + + # run & assert + self._test_4d_label(wanted_global, self.fake_4d_label_data, + include_global=True, incl_shared_variance=False) + + @skipif(no_nilearn) + def test_signal_extr_shared(self): + # set up + wanted = [] + for vol in range(self.fake_fmri_data.shape[3]): + volume = self.fake_fmri_data[:, :, :, vol].flatten() + wanted_row = [] + for reg in range(self.fake_4d_label_data.shape[3]): + region = self.fake_4d_label_data[:, :, :, reg].flatten() + wanted_row.append((volume*region).sum()/(region*region).sum()) + + wanted.append(wanted_row) + # run & assert + self._test_4d_label(wanted, self.fake_4d_label_data) + + def _test_4d_label(self, wanted, fake_labels, include_global=False, incl_shared_variance=True): + # set up + utils.save_toy_nii(fake_labels, self.filenames['4d_label_file']) + + # run + iface.SignalExtraction(in_file=self.filenames['in_file'], + label_files=self.filenames['4d_label_file'], + class_labels=self.labels, + incl_shared_variance=incl_shared_variance, + include_global=include_global).run() + + wanted_labels = self.global_labels if include_global else self.labels + + # assert + self.assert_expected_output(wanted_labels, wanted) + + def assert_expected_output(self, labels, wanted): + with open(self.filenames['out_file'], 'r') as output: + got = [line.split() for line in output] + labels_got = got.pop(0) # remove header + assert_equal(labels_got, labels) + assert_equal(len(got), self.fake_fmri_data.shape[3], + 'num rows and num volumes') + # convert from string to float + got = [[float(num) for num in row] for row in got] + for i, time in enumerate(got): + assert_equal(len(labels), len(time)) + for j, segment in enumerate(time): + assert_almost_equal(segment, wanted[i][j], decimal=1) + + + def tearDown(self): + os.chdir(self.orig_dir) + shutil.rmtree(self.temp_dir) + + fake_fmri_data = np.array([[[[2, -1, 4, -2, 3], + [4, -2, -5, -1, 0]], + + [[-2, 0, 1, 4, 4], + [-5, 3, -3, 1, -5]]], + + + [[[2, -2, -1, -2, -5], + [3, 0, 3, -5, -2]], + + [[-4, -2, -2, 1, -2], + [3, 1, 4, -3, -2]]]]) + + fake_label_data = np.array([[[1, 0], + [3, 1]], + + [[2, 0], + [1, 3]]]) + + fake_equiv_4d_label_data = np.array([[[[1., 0., 0.], + [0., 0., 0.]], + [[0., 0., 1.], + [1., 0., 0.]]], + [[[0., 1., 0.], + [0., 0., 0.]], + [[1., 0., 0.], + [0., 0., 1.]]]]) + + base_wanted = [[-2.33333, 2, .5], + [0, -2, .5], + [-.3333333, -1, 2.5], + [0, -2, .5], + [-1.3333333, -5, 1]] + + fake_4d_label_data = np.array([[[[0.2, 0.3, 0.5], + [0.1, 0.1, 0.8]], + + [[0.1, 0.3, 0.6], + [0.3, 0.4, 0.3]]], + + [[[0.2, 0.2, 0.6], + [0., 0.3, 0.7]], + + [[0.3, 0.3, 0.4], + [0.3, 0.4, 0.3]]]]) + + + fourd_wanted = [[-5.0652173913, -5.44565217391, 5.50543478261], + [-7.02173913043, 11.1847826087, -4.33152173913], + [-19.0869565217, 21.2391304348, -4.57608695652], + [5.19565217391, -3.66304347826, -1.51630434783], + [-12.0, 3., 0.5]]