From 50da63f9a8c74e60ef43b48349d9e5e19f5ad719 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Thu, 1 Sep 2016 00:00:19 +0000 Subject: [PATCH 01/61] compcor skeleton --- nipype/algorithms/compcor.py | 56 ++++++++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) create mode 100644 nipype/algorithms/compcor.py diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py new file mode 100644 index 0000000000..faddc3ae54 --- /dev/null +++ b/nipype/algorithms/compcor.py @@ -0,0 +1,56 @@ +from __future__ import division +from builtins import range +from numpy import ones, kron, mean, eye, hstack, dot, tile +from scipy.linalg import pinv +from ..interfaces.base import BaseInterfaceInputSpec, TraitedSpec, \ + BaseInterface, traits, File +import nibabel as nb +import numpy as np +import os + +class CompCoreInputSpec(BaseInterfaceInputSpec): + realigned_file = File(exists=True, mandatory=True, desc='already realigned brain image') + mask_file = File(exists=True, mandatory=True, desc='mask file that determines ROI') + num_components = traits.Int(default=6, usedefault=True) # 6 for BOLD, 4 for ASL + +class CompCoreOutputSpec(TraitedSpec): + components_file = File(desc='text file containing the noise components') + +class CompCore(BaseInterface): + ''' + Node with core CompCor computation, used in aCompCor and tCompCor + ''' + input_spec = CompCoreInputSpec + output_spec = CompCoreOutputSpec + + def _run_interface(self, runtime): + # just copied for now + imgseries = nb.load(realigned_file) + components = None + mask = nb.load(mask_file).get_data() + + voxel_timecourses = imgseries.get_data()[mask > 0] + voxel_timecourses[np.isnan(np.sum(voxel_timecourses, axis=1)), :] = 0 + # remove mean and normalize by variance + # voxel_timecourses.shape == [nvoxels, time] + X = voxel_timecourses.T + stdX = np.std(X, axis=0) + stdX[stdX == 0] = 1. + stdX[np.isnan(stdX)] = 1. + stdX[np.isinf(stdX)] = 1. + X = (X - np.mean(X, axis=0)) / stdX + u, _, _ = sp.linalg.svd(X, full_matrices=False) + if components is None: + components = u[:, :num_components] + else: + components = np.hstack((components, u[:, :num_components])) + + components_file = os.path.join(os.getcwd(), 'noise_components.txt') + np.savetxt(components_file, components, fmt="%.10f") + return components_file + +class aCompCor(Node): + pass + +class tCompCor(Node): + pass From 9e0d32796cec4ef237e2ccb2b7b74554f01c6ad8 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Fri, 2 Sep 2016 04:36:14 +0000 Subject: [PATCH 02/61] incremental changes, baby test --- nipype/algorithms/compcor.py | 18 ++++++++++-------- nipype/algorithms/tests/test_compcor.py | 22 ++++++++++++++++++++++ 2 files changed, 32 insertions(+), 8 deletions(-) create mode 100644 nipype/algorithms/tests/test_compcor.py diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index faddc3ae54..1411106f9f 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -8,10 +8,13 @@ import numpy as np import os +from nipype.pipeline.engine import Workflow + class CompCoreInputSpec(BaseInterfaceInputSpec): realigned_file = File(exists=True, mandatory=True, desc='already realigned brain image') mask_file = File(exists=True, mandatory=True, desc='mask file that determines ROI') num_components = traits.Int(default=6, usedefault=True) # 6 for BOLD, 4 for ASL + # additional_regressors?? class CompCoreOutputSpec(TraitedSpec): components_file = File(desc='text file containing the noise components') @@ -24,10 +27,9 @@ class CompCore(BaseInterface): output_spec = CompCoreOutputSpec def _run_interface(self, runtime): - # just copied for now - imgseries = nb.load(realigned_file) + imgseries = nb.load(self.inputs.realigned_file) components = None - mask = nb.load(mask_file).get_data() + mask = nb.load(self.inputs.realigned_file, self.inputs.mask_file).get_data() voxel_timecourses = imgseries.get_data()[mask > 0] voxel_timecourses[np.isnan(np.sum(voxel_timecourses, axis=1)), :] = 0 @@ -41,16 +43,16 @@ def _run_interface(self, runtime): X = (X - np.mean(X, axis=0)) / stdX u, _, _ = sp.linalg.svd(X, full_matrices=False) if components is None: - components = u[:, :num_components] + components = u[:, :self.inputs.num_components] else: - components = np.hstack((components, u[:, :num_components])) + components = np.hstack((components, u[:, :self.inputs.num_components])) components_file = os.path.join(os.getcwd(), 'noise_components.txt') np.savetxt(components_file, components, fmt="%.10f") - return components_file + return runtime -class aCompCor(Node): +class aCompCor(Workflow): pass -class tCompCor(Node): +class tCompCor(Workflow): pass diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py new file mode 100644 index 0000000000..f758c5d7ab --- /dev/null +++ b/nipype/algorithms/tests/test_compcor.py @@ -0,0 +1,22 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +from __future__ import print_function +from builtins import zip +from builtins import range +from builtins import open + +import os +import glob +import shutil +import os.path as op +from tempfile import mkstemp, mkdtemp +from subprocess import Popen + +from nose.tools import assert_raises +import nipype +from nipype.testing import assert_equal, assert_true, assert_false, skipif +from nipype.algorithms.compcor import CompCore + +def test_compcore(): + ccinterface = CompCore(realigned_file='../../testing/data/functional.nii', mask_file='../../testing/data/mask.nii') + ccresult = ccinterface.run() From 9f7886c0d928ee0838cb3218b0a576be982f4f61 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Fri, 2 Sep 2016 19:34:14 +0000 Subject: [PATCH 03/61] add doctest --- nipype/algorithms/compcor.py | 9 +++++++++ nipype/algorithms/tests/test_compcor.py | 1 + 2 files changed, 10 insertions(+) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index 1411106f9f..4aee03faca 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -22,6 +22,15 @@ class CompCoreOutputSpec(TraitedSpec): class CompCore(BaseInterface): ''' Node with core CompCor computation, used in aCompCor and tCompCor + + Example + ------- + + >>> ccinterface = CompCore() + >>> ccinterface.inputs.realigned_file = '../../testing/data/functional.nii' + >>> ccinterface.inputs.mask_file = '../../testing/data/mask.nii' + >>> ccinterface.inputs.num_components = 1 + ''' input_spec = CompCoreInputSpec output_spec = CompCoreOutputSpec diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index f758c5d7ab..34535180ac 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -17,6 +17,7 @@ from nipype.testing import assert_equal, assert_true, assert_false, skipif from nipype.algorithms.compcor import CompCore +@skipif(True) def test_compcore(): ccinterface = CompCore(realigned_file='../../testing/data/functional.nii', mask_file='../../testing/data/mask.nii') ccresult = ccinterface.run() From 55aec65e303d8b6311c48d684bdf4f555884c4cc Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Fri, 2 Sep 2016 19:34:14 +0000 Subject: [PATCH 04/61] add doctest --- nipype/algorithms/compcor.py | 9 +++++++++ nipype/algorithms/tests/test_compcor.py | 1 + 2 files changed, 10 insertions(+) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index 1411106f9f..54c5d8f0ce 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -22,6 +22,15 @@ class CompCoreOutputSpec(TraitedSpec): class CompCore(BaseInterface): ''' Node 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 diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index f758c5d7ab..34535180ac 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -17,6 +17,7 @@ from nipype.testing import assert_equal, assert_true, assert_false, skipif from nipype.algorithms.compcor import CompCore +@skipif(True) def test_compcore(): ccinterface = CompCore(realigned_file='../../testing/data/functional.nii', mask_file='../../testing/data/mask.nii') ccresult = ccinterface.run() From f454cbed6a7bb95b984fa0124dd19cab87ac7f50 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Sat, 3 Sep 2016 06:22:43 +0000 Subject: [PATCH 05/61] fix imports --- nipype/algorithms/compcor.py | 5 +---- nipype/algorithms/tests/test_compcor.py | 13 ------------- 2 files changed, 1 insertion(+), 17 deletions(-) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index 54c5d8f0ce..84fabcdfa3 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -1,7 +1,4 @@ -from __future__ import division -from builtins import range -from numpy import ones, kron, mean, eye, hstack, dot, tile -from scipy.linalg import pinv +from scipy.linalg import svd # singular value decomposition from ..interfaces.base import BaseInterfaceInputSpec, TraitedSpec, \ BaseInterface, traits, File import nibabel as nb diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index 34535180ac..d4edd1991f 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -1,18 +1,5 @@ # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: -from __future__ import print_function -from builtins import zip -from builtins import range -from builtins import open - -import os -import glob -import shutil -import os.path as op -from tempfile import mkstemp, mkdtemp -from subprocess import Popen - -from nose.tools import assert_raises import nipype from nipype.testing import assert_equal, assert_true, assert_false, skipif from nipype.algorithms.compcor import CompCore From ca930d364813fa0249856cf49921ca16dcdfff55 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Sat, 3 Sep 2016 06:31:22 +0000 Subject: [PATCH 06/61] better comments, minor re-arranging --- nipype/algorithms/compcor.py | 56 +++++++++++++++---------- nipype/algorithms/tests/test_compcor.py | 15 ++++++- 2 files changed, 47 insertions(+), 24 deletions(-) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index 84fabcdfa3..b4bfabd7fa 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -8,17 +8,17 @@ from nipype.pipeline.engine import Workflow class CompCoreInputSpec(BaseInterfaceInputSpec): - realigned_file = File(exists=True, mandatory=True, desc='already realigned brain image') - mask_file = File(exists=True, mandatory=True, desc='mask file that determines ROI') + 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(default=6, usedefault=True) # 6 for BOLD, 4 for ASL # additional_regressors?? class CompCoreOutputSpec(TraitedSpec): - components_file = File(desc='text file containing the noise components') + components_file = File(desc='text file containing the noise components', exists=True) class CompCore(BaseInterface): ''' - Node with core CompCor computation, used in aCompCor and tCompCor + Interface with core CompCor computation, used in aCompCor and tCompCor Example ------- @@ -27,33 +27,43 @@ class CompCore(BaseInterface): >>> 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) - components = None - mask = nb.load(self.inputs.realigned_file, self.inputs.mask_file).get_data() - - voxel_timecourses = imgseries.get_data()[mask > 0] + 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 - # remove mean and normalize by variance + + # 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] - X = voxel_timecourses.T - stdX = np.std(X, axis=0) - stdX[stdX == 0] = 1. - stdX[np.isnan(stdX)] = 1. - stdX[np.isinf(stdX)] = 1. - X = (X - np.mean(X, axis=0)) / stdX - u, _, _ = sp.linalg.svd(X, full_matrices=False) - if components is None: - components = u[:, :self.inputs.num_components] - else: - components = np.hstack((components, u[:, :self.inputs.num_components])) + M = voxel_timecourses.T + + 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. + + # "The constant and linear trends of the columns in the matrix M were + # removed prior to column-wise variance normalization." + M = (M - np.mean(M, axis=0)) / stdM + print(M) + # "The covariance matrix C = MMT was constructed and decomposed into its + # principal components using a singular value decomposition." + u, _, _ = svd(M, full_matrices=False) + + components = u[:, :self.inputs.num_components] - components_file = os.path.join(os.getcwd(), 'noise_components.txt') + components_file = os.path.join(os.getcwd(), + self._outputs().get()["components_file"]) np.savetxt(components_file, components, fmt="%.10f") return runtime diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index d4edd1991f..26df22a0e2 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -4,7 +4,20 @@ 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 + +# first 3 = spatial; last = temporal +dim = 2, 3, 4, 5 + @skipif(True) def test_compcore(): - ccinterface = CompCore(realigned_file='../../testing/data/functional.nii', mask_file='../../testing/data/mask.nii') + ccinterface = CompCore(realigned_file=make_toy(np.random.rand(*dim), 'func.nii'), + mask_file=make_toy(np.random.randint(0, 2, dim[:2]), + 'mask.nii')) ccresult = ccinterface.run() + +def make_toy(array, filename): + toy = nb.Nifti1Image(array, np.eye(4)) + nb.nifti1.save(toy, filename) + return filename From 414f11a84b9861a81ab122dc20e2146c0dcb98cb Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Sun, 4 Sep 2016 03:38:09 +0000 Subject: [PATCH 07/61] set up test prereq's --- nipype/algorithms/compcor.py | 16 ++++++------- nipype/algorithms/tests/test_compcor.py | 30 ++++++++++++++++++++++--- 2 files changed, 34 insertions(+), 12 deletions(-) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index b4bfabd7fa..57775b8858 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -10,7 +10,7 @@ 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(default=6, usedefault=True) # 6 for BOLD, 4 for ASL + num_components = traits.Int(default=6) # 6 for BOLD, 4 for ASL # additional_regressors?? class CompCoreOutputSpec(TraitedSpec): @@ -45,25 +45,23 @@ def _run_interface(self, runtime): # voxel_timecourses.shape == [nvoxels, time] M = voxel_timecourses.T - stdM = np.std(M, axis=0) + # "The constant and linear trends of the columns in the matrix M were removed ..." + M = (M - np.mean(M, axis=0)) + # "... 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 constant and linear trends of the columns in the matrix M were - # removed prior to column-wise variance normalization." - M = (M - np.mean(M, axis=0)) / stdM - print(M) # "The covariance matrix C = MMT was constructed and decomposed into its # principal components using a singular value decomposition." u, _, _ = svd(M, full_matrices=False) components = u[:, :self.inputs.num_components] - - components_file = os.path.join(os.getcwd(), - self._outputs().get()["components_file"]) + components_file = os.path.join(os.getcwd(), "components_file.txt") np.savetxt(components_file, components, fmt="%.10f") return runtime diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index 26df22a0e2..ae0d89b3d9 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -12,12 +12,36 @@ @skipif(True) def test_compcore(): - ccinterface = CompCore(realigned_file=make_toy(np.random.rand(*dim), 'func.nii'), - mask_file=make_toy(np.random.randint(0, 2, dim[:2]), - 'mask.nii')) + # setup + noise = np.fromfunction(fake_noise_fun, fake_data.shape) + realigned_file = make_toy(fake_data + noise, 'func.nii') + + mask = np.ones(fake_data.shape[:3]) + mask[0,0,0] = 0 + mask[0,0,1] = 0 + mask_file = make_toy(mask, 'mask.nii') + + # run + ccinterface = CompCore(realigned_file=realigned_file, mask_file=mask_file) ccresult = ccinterface.run() def make_toy(array, filename): toy = nb.Nifti1Image(array, 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 From bca197e07139f43fed0161fc2cde8affcc6010a8 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Sun, 4 Sep 2016 22:35:39 +0000 Subject: [PATCH 08/61] complete _run_interface implementation --- nipype/algorithms/compcor.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index 57775b8858..09a999e19f 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -1,8 +1,8 @@ -from scipy.linalg import svd # singular value decomposition -from ..interfaces.base import BaseInterfaceInputSpec, TraitedSpec, \ - BaseInterface, traits, File +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 @@ -10,7 +10,7 @@ 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(default=6) # 6 for BOLD, 4 for ASL + num_components = traits.Int(6, usedefault=True) # 6 for BOLD, 4 for ASL # additional_regressors?? class CompCoreOutputSpec(TraitedSpec): @@ -44,9 +44,14 @@ def _run_interface(self, runtime): # 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 ..." - M = (M - np.mean(M, axis=0)) + 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) @@ -58,8 +63,7 @@ def _run_interface(self, runtime): # "The covariance matrix C = MMT was constructed and decomposed into its # principal components using a singular value decomposition." - u, _, _ = svd(M, full_matrices=False) - + 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") From 124d3461fd1e6cbe1886948f3a9ca44feb90a6c0 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Sun, 4 Sep 2016 22:39:30 +0000 Subject: [PATCH 09/61] implement _list_outputs --- nipype/algorithms/compcor.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index 09a999e19f..23d7c28240 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -69,6 +69,11 @@ def _run_interface(self, runtime): 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 From ef50858d3b3016018f0078af6b7da29412803c00 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Sun, 4 Sep 2016 23:57:33 +0000 Subject: [PATCH 10/61] pretty ok test --- nipype/algorithms/tests/test_compcor.py | 27 ++++++++++++++++++------- 1 file changed, 20 insertions(+), 7 deletions(-) diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index ae0d89b3d9..cf54d62c58 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -6,27 +6,40 @@ import nibabel as nb import numpy as np +import os -# first 3 = spatial; last = temporal -dim = 2, 3, 4, 5 +functionalnii = 'func.nii' +masknii = 'mask.nii' -@skipif(True) def test_compcore(): # setup noise = np.fromfunction(fake_noise_fun, fake_data.shape) - realigned_file = make_toy(fake_data + noise, 'func.nii') + 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, 'mask.nii') + mask_file = make_toy(mask, masknii) # run ccinterface = CompCore(realigned_file=realigned_file, mask_file=mask_file) ccresult = ccinterface.run() -def make_toy(array, filename): - toy = nb.Nifti1Image(array, np.eye(4)) + # 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 From 5e8b99775cc400f978c0c7d471439145b0c36bee Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Tue, 6 Sep 2016 21:39:35 +0000 Subject: [PATCH 11/61] compcore -> compcor --- nipype/algorithms/compcor.py | 18 ++++++++++-------- nipype/algorithms/tests/test_compcor.py | 8 +++++--- 2 files changed, 15 insertions(+), 11 deletions(-) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index 23d7c28240..a926b433b1 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -1,35 +1,37 @@ from ..interfaces.base import (BaseInterfaceInputSpec, TraitedSpec, BaseInterface, traits, File) + +from nipype.pipeline import engine as pe +from nipype.interfaces.utility import IdentityInterface + import nibabel as nb import numpy as np from scipy import linalg, stats import os -from nipype.pipeline.engine import Workflow - -class CompCoreInputSpec(BaseInterfaceInputSpec): +class CompCorInputSpec(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): +class CompCorOutputSpec(TraitedSpec): components_file = File(desc='text file containing the noise components', exists=True) -class CompCore(BaseInterface): +class CompCor(BaseInterface): ''' Interface with core CompCor computation, used in aCompCor and tCompCor Example ------- - >>> ccinterface = CompCore() + >>> ccinterface = CompCor() >>> 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 + input_spec = CompCorInputSpec + output_spec = CompCorOutputSpec def _run_interface(self, runtime): imgseries = nb.load(self.inputs.realigned_file).get_data() diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index cf54d62c58..ad5aca8fc6 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -2,8 +2,10 @@ # 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 +from nipype.algorithms.compcor import (CompCor, CompCorWorkflow, ACompCor, + TCompCor) +import mock import nibabel as nb import numpy as np import os @@ -11,7 +13,7 @@ functionalnii = 'func.nii' masknii = 'mask.nii' -def test_compcore(): +def test_compcor(): # setup noise = np.fromfunction(fake_noise_fun, fake_data.shape) realigned_file = make_toy(fake_data + noise, functionalnii) @@ -22,7 +24,7 @@ def test_compcore(): mask_file = make_toy(mask, masknii) # run - ccinterface = CompCore(realigned_file=realigned_file, mask_file=mask_file) + ccinterface = CompCor(realigned_file=realigned_file, mask_file=mask_file) ccresult = ccinterface.run() # asserts From ac8dbfe1f8dbe1fffbecefebd09f14b0d4ce3c4b Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Wed, 7 Sep 2016 03:44:44 +0000 Subject: [PATCH 12/61] TCompCor skeleton, tests --- nipype/algorithms/compcor.py | 12 ++-- nipype/algorithms/tests/test_compcor.py | 96 +++++++++++++++---------- 2 files changed, 64 insertions(+), 44 deletions(-) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index a926b433b1..e213d612a1 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -11,7 +11,7 @@ class CompCorInputSpec(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)') + mask_file = File(exists=True, mandatory=False, desc='mask file that determines ROI (3D)') num_components = traits.Int(6, usedefault=True) # 6 for BOLD, 4 for ASL # additional_regressors?? @@ -76,8 +76,10 @@ def _list_outputs(self): outputs['components_file'] = os.path.abspath("components_file.txt") return outputs -class aCompCor(Workflow): - pass +class TCompCor(CompCor): -class tCompCor(Workflow): - pass + def _run_interface(self, runtime): + # create mask here + + super(TCompCor, self)._run_interface(runtime) + return runtime diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index ad5aca8fc6..0862ac1205 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -2,61 +2,79 @@ # 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 (CompCor, CompCorWorkflow, ACompCor, - TCompCor) +from nipype.algorithms.compcor import CompCor, TCompCor +import unittest import mock import nibabel as nb import numpy as np import os -functionalnii = 'func.nii' -masknii = 'mask.nii' +class TestCompCor(unittest.TestCase): + ''' Note: Tests currently do a poor job of testing functionality ''' -def test_compcor(): - # setup - noise = np.fromfunction(fake_noise_fun, fake_data.shape) - realigned_file = make_toy(fake_data + noise, functionalnii) + functionalnii = 'func.nii' + masknii = 'mask.nii' + components_file = None - mask = np.ones(fake_data.shape[:3]) - mask[0,0,0] = 0 - mask[0,0,1] = 0 - mask_file = make_toy(mask, masknii) + def setUp(self): + # setup + noise = np.fromfunction(self.fake_noise_fun, self.fake_data.shape) + self.realigned_file = self.make_toy(self.fake_data + noise, + self.functionalnii) - # run - ccinterface = CompCor(realigned_file=realigned_file, mask_file=mask_file) - ccresult = ccinterface.run() + def test_compcor(self): + mask = np.ones(self.fake_data.shape[:3]) + mask[0,0,0] = 0 + mask[0,0,1] = 0 + mask_file = self.make_toy(mask, self.masknii) - # 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) + ccinterface = CompCor(realigned_file=self.realigned_file, + mask_file=mask_file) + self.meat(ccinterface) - # apply components_file to realigned_file, is it better? the same as before adding noise? + @skipif(True) + def test_tcompcor(self): + ccinterface = TCompCor(realigned_file=self.realigned_file) + self.meat(ccinterface) - # remove temporary nifti files - os.remove(functionalnii) - os.remove(masknii) - os.remove(components) + def meat(self, ccinterface): + # run + ccresult = ccinterface.run() -def make_toy(ndarray, filename): - toy = nb.Nifti1Image(ndarray, np.eye(4)) - nb.nifti1.save(toy, filename) - return filename + # assert + print(ccresult.outputs.components_file) + self.components_file = ccinterface._list_outputs()['components_file'] + assert_equal(ccresult.outputs.components_file, self.components_file) + assert_true(os.path.exists(self.components_file)) + assert_true(os.path.getsize(self.components_file) > 0) + assert_equal(ccinterface.inputs.num_components, 6) -fake_data = np.array([[[[8, 5, 3, 8, 0], - [6, 7, 4, 7, 1]], + def tearDown(self): + # remove temporary nifti files + try: + os.remove(self.functionalnii) + os.remove(self.components_file) + os.remove(self.masknii) + except (FileNotFoundError, TypeError) as e: + print(e) - [[7, 9, 1, 6, 5], - [0, 7, 4, 7, 7]]], + def make_toy(self, ndarray, filename): + toy = nb.Nifti1Image(ndarray, np.eye(4)) + nb.nifti1.save(toy, filename) + return filename + def fake_noise_fun(self, i, j, l, m): + return m*i + l - j - [[[2, 4, 5, 7, 0], - [1, 7, 0, 5, 4]], + fake_data = np.array([[[[8, 5, 3, 8, 0], + [6, 7, 4, 7, 1]], - [[7, 3, 9, 0, 4], - [9, 4, 1, 5, 0]]]]) + [[7, 9, 1, 6, 5], + [0, 7, 4, 7, 7]]], -def fake_noise_fun(i, j, l, m): - return m*i + l - j + [[[2, 4, 5, 7, 0], + [1, 7, 0, 5, 4]], + + [[7, 3, 9, 0, 4], + [9, 4, 1, 5, 0]]]]) From a5402f7f5a2ac0ec42becd627b21066c04783625 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Wed, 7 Sep 2016 03:54:52 +0000 Subject: [PATCH 13/61] make python2 happy --- nipype/algorithms/tests/test_compcor.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index 0862ac1205..4c03fd775b 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -56,7 +56,7 @@ def tearDown(self): os.remove(self.functionalnii) os.remove(self.components_file) os.remove(self.masknii) - except (FileNotFoundError, TypeError) as e: + except (OSError, TypeError) as e: print(e) def make_toy(self, ndarray, filename): From 8f5eabcfc98d58003332c4ee01175660aa95be6a Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Wed, 7 Sep 2016 04:18:55 +0000 Subject: [PATCH 14/61] factor out tSTD calculation --- nipype/algorithms/compcor.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index e213d612a1..d3cb5d1808 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -56,12 +56,7 @@ def _run_interface(self, runtime): 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 + M = M / self._compute_tSTD(M) # "The covariance matrix C = MMT was constructed and decomposed into its # principal components using a singular value decomposition." @@ -76,6 +71,14 @@ def _list_outputs(self): outputs['components_file'] = os.path.abspath("components_file.txt") return outputs + def _compute_tSTD(self, matrix): + 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. + return stdM + class TCompCor(CompCor): def _run_interface(self, runtime): From d87ac7160fedffccf2b975f6a7f75557d8222266 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Wed, 7 Sep 2016 05:50:59 +0000 Subject: [PATCH 15/61] tCompCor first pass --- nipype/algorithms/compcor.py | 38 +++++++++++++++++++------ nipype/algorithms/tests/test_compcor.py | 1 - 2 files changed, 30 insertions(+), 9 deletions(-) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index d3cb5d1808..4fb557e382 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -56,7 +56,7 @@ def _run_interface(self, runtime): M[:, voxel] = M[:, voxel] - [m*t + b for t in timesteps] # "... prior to column-wise variance normalization." - M = M / self._compute_tSTD(M) + M = M / self._compute_tSTD(M, 1.) # "The covariance matrix C = MMT was constructed and decomposed into its # principal components using a singular value decomposition." @@ -71,18 +71,40 @@ def _list_outputs(self): outputs['components_file'] = os.path.abspath("components_file.txt") return outputs - def _compute_tSTD(self, matrix): + def _compute_tSTD(self, M, x): 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. + # set bad values to x + stdM[stdM == 0] = x + stdM[np.isnan(stdM)] = x + stdM[np.isinf(stdM)] = x return stdM class TCompCor(CompCor): def _run_interface(self, runtime): - # create mask here - + imgseries = nb.load(self.inputs.realigned_file).get_data() + time_voxels = imgseries.T + num_voxels = np.prod(time_voxels.shape[1:]) + + # 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)." + + # "To construct the tSTD noise ROI, we sorted the voxels by their + # temporal standard deviation ..." + tSTD = self._compute_tSTD(time_voxels, 0) + sortSTD = np.sort(tSTD, axis=None) # flattened sorted matrix + + # "... and retained a pre-specified upper fraction of the sorted voxels + # within each slice ... we chose a 2% threshold" + threshold = sortSTD[int(num_voxels * .98)] + mask = tSTD >= threshold + mask = mask.astype(int) + + # save mask + mask_file = 'mask.nii' + nb.nifti1.save(nb.Nifti1Image(mask, np.eye(4)), mask_file) + self.inputs.mask_file = 'mask.nii' super(TCompCor, self)._run_interface(runtime) return runtime diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index 4c03fd775b..1f2cd313d7 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -33,7 +33,6 @@ def test_compcor(self): mask_file=mask_file) self.meat(ccinterface) - @skipif(True) def test_tcompcor(self): ccinterface = TCompCor(realigned_file=self.realigned_file) self.meat(ccinterface) From 757bd0e011762b15aa18e52f79a5431076b1369d Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Wed, 7 Sep 2016 06:23:47 +0000 Subject: [PATCH 16/61] alias ACompCor => CompCor --- nipype/algorithms/compcor.py | 15 ++++++++----- nipype/algorithms/tests/test_compcor.py | 29 +++++++++++++++---------- 2 files changed, 28 insertions(+), 16 deletions(-) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index 4fb557e382..c2eaced2c4 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -10,10 +10,13 @@ import os class CompCorInputSpec(BaseInterfaceInputSpec): - realigned_file = File(exists=True, mandatory=True, desc='already realigned brain image (4D)') - mask_file = File(exists=True, mandatory=False, desc='mask file that determines ROI (3D)') + realigned_file = File(exists=True, mandatory=True, + desc='already realigned brain image (4D)') + mask_file = File(exists=True, mandatory=False, + desc='mask file that determines ROI (3D)') + components_file = File('components_file.txt', exists=False, mandatory=False, usedefault=True, + desc='filename to store physiological components in') num_components = traits.Int(6, usedefault=True) # 6 for BOLD, 4 for ASL - # additional_regressors?? class CompCorOutputSpec(TraitedSpec): components_file = File(desc='text file containing the noise components', exists=True) @@ -62,13 +65,13 @@ def _run_interface(self, runtime): # 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") + components_file = os.path.join(os.getcwd(), self.inputs.components_file) 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") + outputs['components_file'] = os.path.abspath(self.inputs.components_file) return outputs def _compute_tSTD(self, M, x): @@ -108,3 +111,5 @@ def _run_interface(self, runtime): self.inputs.mask_file = 'mask.nii' super(TCompCor, self)._run_interface(runtime) return runtime + +ACompCor = CompCor diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index 1f2cd313d7..6c0f922dd5 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -2,7 +2,7 @@ # 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 CompCor, TCompCor +from nipype.algorithms.compcor import CompCor, TCompCor, ACompCor import unittest import mock @@ -29,26 +29,33 @@ def test_compcor(self): mask[0,0,1] = 0 mask_file = self.make_toy(mask, self.masknii) - ccinterface = CompCor(realigned_file=self.realigned_file, - mask_file=mask_file) - self.meat(ccinterface) + ccresult = self.run_cc(CompCor(realigned_file=self.realigned_file, + mask_file=mask_file)) + + accresult = self.run_cc(ACompCor(realigned_file=self.realigned_file, + mask_file=mask_file, + components_file='acc_components_file')) + + assert_equal(os.path.getsize(ccresult.outputs.components_file), + os.path.getsize(accresult.outputs.components_file)) def test_tcompcor(self): ccinterface = TCompCor(realigned_file=self.realigned_file) - self.meat(ccinterface) + self.run_cc(ccinterface) - def meat(self, ccinterface): + def run_cc(self, ccinterface): # run ccresult = ccinterface.run() # assert - print(ccresult.outputs.components_file) - self.components_file = ccinterface._list_outputs()['components_file'] - assert_equal(ccresult.outputs.components_file, self.components_file) - assert_true(os.path.exists(self.components_file)) - assert_true(os.path.getsize(self.components_file) > 0) + expected_file = ccinterface._list_outputs()['components_file'] + assert_equal(ccresult.outputs.components_file, expected_file) + assert_true(os.path.exists(expected_file)) + assert_true(os.path.getsize(expected_file) > 0) assert_equal(ccinterface.inputs.num_components, 6) + return ccresult + def tearDown(self): # remove temporary nifti files try: From a54120f65139af03ef39f88f3d47bef7259ce115 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Wed, 7 Sep 2016 06:39:52 +0000 Subject: [PATCH 17/61] prevent name collision in tests --- nipype/algorithms/tests/test_compcor.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index 6c0f922dd5..e2f5df59cd 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -14,7 +14,7 @@ class TestCompCor(unittest.TestCase): ''' Note: Tests currently do a poor job of testing functionality ''' functionalnii = 'func.nii' - masknii = 'mask.nii' + masknii = 'amask.nii' components_file = None def setUp(self): From c1b0a8705acb652c4f81bd84e527011ddaa2148b Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Wed, 7 Sep 2016 21:59:35 +0000 Subject: [PATCH 18/61] better unique file names --- nipype/algorithms/tests/test_compcor.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index e2f5df59cd..75bf822706 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -13,8 +13,8 @@ class TestCompCor(unittest.TestCase): ''' Note: Tests currently do a poor job of testing functionality ''' - functionalnii = 'func.nii' - masknii = 'amask.nii' + functionalnii = 'compcorfunc.nii' + masknii = 'compcormask.nii' components_file = None def setUp(self): From 9a88a118aa8f67663529f8401631b1d0b65924a2 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Wed, 7 Sep 2016 23:08:35 +0000 Subject: [PATCH 19/61] refactor toy maker function --- nipype/algorithms/tests/test_compcor.py | 14 ++-- nipype/algorithms/tests/test_tsnr.py | 86 +++++++++++++++++++++++++ nipype/testing/utils.py | 8 +++ 3 files changed, 99 insertions(+), 9 deletions(-) create mode 100644 nipype/algorithms/tests/test_tsnr.py diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index 75bf822706..a5165e0adf 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -1,7 +1,8 @@ # 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.testing import (assert_equal, assert_true, assert_false, skipif, + utils) from nipype.algorithms.compcor import CompCor, TCompCor, ACompCor import unittest @@ -20,14 +21,14 @@ class TestCompCor(unittest.TestCase): def setUp(self): # setup noise = np.fromfunction(self.fake_noise_fun, self.fake_data.shape) - self.realigned_file = self.make_toy(self.fake_data + noise, - self.functionalnii) + self.realigned_file = utils.save_toy_nii(self.fake_data + noise, + self.functionalnii) def test_compcor(self): mask = np.ones(self.fake_data.shape[:3]) mask[0,0,0] = 0 mask[0,0,1] = 0 - mask_file = self.make_toy(mask, self.masknii) + mask_file = utils.save_toy_nii(mask, self.masknii) ccresult = self.run_cc(CompCor(realigned_file=self.realigned_file, mask_file=mask_file)) @@ -65,11 +66,6 @@ def tearDown(self): except (OSError, TypeError) as e: print(e) - def make_toy(self, ndarray, filename): - toy = nb.Nifti1Image(ndarray, np.eye(4)) - nb.nifti1.save(toy, filename) - return filename - def fake_noise_fun(self, i, j, l, m): return m*i + l - j diff --git a/nipype/algorithms/tests/test_tsnr.py b/nipype/algorithms/tests/test_tsnr.py new file mode 100644 index 0000000000..75bf822706 --- /dev/null +++ b/nipype/algorithms/tests/test_tsnr.py @@ -0,0 +1,86 @@ +# 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 CompCor, TCompCor, ACompCor + +import unittest +import mock +import nibabel as nb +import numpy as np +import os + +class TestCompCor(unittest.TestCase): + ''' Note: Tests currently do a poor job of testing functionality ''' + + functionalnii = 'compcorfunc.nii' + masknii = 'compcormask.nii' + components_file = None + + def setUp(self): + # setup + noise = np.fromfunction(self.fake_noise_fun, self.fake_data.shape) + self.realigned_file = self.make_toy(self.fake_data + noise, + self.functionalnii) + + def test_compcor(self): + mask = np.ones(self.fake_data.shape[:3]) + mask[0,0,0] = 0 + mask[0,0,1] = 0 + mask_file = self.make_toy(mask, self.masknii) + + ccresult = self.run_cc(CompCor(realigned_file=self.realigned_file, + mask_file=mask_file)) + + accresult = self.run_cc(ACompCor(realigned_file=self.realigned_file, + mask_file=mask_file, + components_file='acc_components_file')) + + assert_equal(os.path.getsize(ccresult.outputs.components_file), + os.path.getsize(accresult.outputs.components_file)) + + def test_tcompcor(self): + ccinterface = TCompCor(realigned_file=self.realigned_file) + self.run_cc(ccinterface) + + def run_cc(self, ccinterface): + # run + ccresult = ccinterface.run() + + # assert + expected_file = ccinterface._list_outputs()['components_file'] + assert_equal(ccresult.outputs.components_file, expected_file) + assert_true(os.path.exists(expected_file)) + assert_true(os.path.getsize(expected_file) > 0) + assert_equal(ccinterface.inputs.num_components, 6) + + return ccresult + + def tearDown(self): + # remove temporary nifti files + try: + os.remove(self.functionalnii) + os.remove(self.components_file) + os.remove(self.masknii) + except (OSError, TypeError) as e: + print(e) + + def make_toy(self, ndarray, filename): + toy = nb.Nifti1Image(ndarray, np.eye(4)) + nb.nifti1.save(toy, filename) + return filename + + def fake_noise_fun(self, i, j, l, m): + return m*i + l - j + + 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]]]]) diff --git a/nipype/testing/utils.py b/nipype/testing/utils.py index d98bda9de6..40c2d21b27 100644 --- a/nipype/testing/utils.py +++ b/nipype/testing/utils.py @@ -15,6 +15,9 @@ from nose import SkipTest from future.utils import raise_from +import numpy as np +import nibabel as nb + def skip_if_no_package(*args, **kwargs): """Raise SkipTest if package_check fails @@ -98,3 +101,8 @@ def __exit__(self, exc_type, exc_val, exc_tb): assert not os.path.exists(self.canary) self.dev_null.close() shutil.rmtree(self.tmpdir) + +def save_toy_nii(ndarray, filename): + toy = nb.Nifti1Image(ndarray, np.eye(4)) + nb.nifti1.save(toy, filename) + return filename From 61da367f3938d810b67eef329bc5bbf3076753dd Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Thu, 8 Sep 2016 01:04:55 +0000 Subject: [PATCH 20/61] testing skeleton for tsnr --- nipype/algorithms/tests/test_tsnr.py | 81 ++++++++-------------------- 1 file changed, 22 insertions(+), 59 deletions(-) diff --git a/nipype/algorithms/tests/test_tsnr.py b/nipype/algorithms/tests/test_tsnr.py index 75bf822706..d90f720ebb 100644 --- a/nipype/algorithms/tests/test_tsnr.py +++ b/nipype/algorithms/tests/test_tsnr.py @@ -1,8 +1,8 @@ # 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 CompCor, TCompCor, ACompCor +from nipype.testing import assert_equal, assert_true, skipif, utils +from nipype.algorithms.misc import TSNR import unittest import mock @@ -10,77 +10,40 @@ import numpy as np import os -class TestCompCor(unittest.TestCase): +class TestTSNR(unittest.TestCase): ''' Note: Tests currently do a poor job of testing functionality ''' - - functionalnii = 'compcorfunc.nii' - masknii = 'compcormask.nii' - components_file = None + ''' + in_file = InputMultiPath(File(exists=True), mandatory=True, + regress_poly = traits.Range(low=1, desc='Remove polynomials') + ''' + in_file_name = 'tsnrinfile.nii' def setUp(self): # setup - noise = np.fromfunction(self.fake_noise_fun, self.fake_data.shape) - self.realigned_file = self.make_toy(self.fake_data + noise, - self.functionalnii) - - def test_compcor(self): - mask = np.ones(self.fake_data.shape[:3]) - mask[0,0,0] = 0 - mask[0,0,1] = 0 - mask_file = self.make_toy(mask, self.masknii) - - ccresult = self.run_cc(CompCor(realigned_file=self.realigned_file, - mask_file=mask_file)) - - accresult = self.run_cc(ACompCor(realigned_file=self.realigned_file, - mask_file=mask_file, - components_file='acc_components_file')) - - assert_equal(os.path.getsize(ccresult.outputs.components_file), - os.path.getsize(accresult.outputs.components_file)) - - def test_tcompcor(self): - ccinterface = TCompCor(realigned_file=self.realigned_file) - self.run_cc(ccinterface) + utils.save_toy_nii(self.fake_data, self.in_file_name) - def run_cc(self, ccinterface): + def test_tsnr(self): + # setup # run - ccresult = ccinterface.run() - + tsnrresult = TSNR(in_file=self.in_file_name, regress_poly=1) # assert - expected_file = ccinterface._list_outputs()['components_file'] - assert_equal(ccresult.outputs.components_file, expected_file) - assert_true(os.path.exists(expected_file)) - assert_true(os.path.getsize(expected_file) > 0) - assert_equal(ccinterface.inputs.num_components, 6) - - return ccresult + # cleanup def tearDown(self): # remove temporary nifti files try: - os.remove(self.functionalnii) - os.remove(self.components_file) - os.remove(self.masknii) + os.remove(self.in_file_name) except (OSError, TypeError) as e: print(e) - def make_toy(self, ndarray, filename): - toy = nb.Nifti1Image(ndarray, np.eye(4)) - nb.nifti1.save(toy, filename) - return filename - - def fake_noise_fun(self, i, j, l, m): - return m*i + l - j - - fake_data = np.array([[[[8, 5, 3, 8, 0], - [6, 7, 4, 7, 1]], + fake_data = np.array([[[[2, 4, 3, 9, 1], + [3, 6, 4, 7, 4]], - [[7, 9, 1, 6, 5], - [0, 7, 4, 7, 7]]], + [[8, 3, 4, 6, 2], + [4, 0, 4, 4, 2]]], - [[[2, 4, 5, 7, 0], - [1, 7, 0, 5, 4]], + [[[9, 7, 5, 5, 7], + [7, 8, 4, 8, 4]], - [[7, 3, 9, 0, 4], - [9, 4, 1, 5, 0]]]]) + [[0, 4, 7, 1, 7], + [6, 8, 8, 8, 7]]]]) From 89b530a0bc16aab1014f0303117e14938b5bf78d Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Thu, 8 Sep 2016 01:19:59 +0000 Subject: [PATCH 21/61] factor out deleting temp files --- nipype/algorithms/tests/test_compcor.py | 18 ++++++------------ nipype/algorithms/tests/test_tsnr.py | 19 ++++++++++++------- nipype/testing/utils.py | 8 ++++++++ 3 files changed, 26 insertions(+), 19 deletions(-) diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index a5165e0adf..ffb5b198b5 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -14,21 +14,21 @@ class TestCompCor(unittest.TestCase): ''' Note: Tests currently do a poor job of testing functionality ''' - functionalnii = 'compcorfunc.nii' - masknii = 'compcormask.nii' - components_file = None + filenames = {'functionalnii': 'compcorfunc.nii', + 'masknii': 'compcormask.nii', + 'components_file': None} def setUp(self): # setup noise = np.fromfunction(self.fake_noise_fun, self.fake_data.shape) self.realigned_file = utils.save_toy_nii(self.fake_data + noise, - self.functionalnii) + self.filenames['functionalnii']) def test_compcor(self): mask = np.ones(self.fake_data.shape[:3]) mask[0,0,0] = 0 mask[0,0,1] = 0 - mask_file = utils.save_toy_nii(mask, self.masknii) + mask_file = utils.save_toy_nii(mask, self.filenames['masknii']) ccresult = self.run_cc(CompCor(realigned_file=self.realigned_file, mask_file=mask_file)) @@ -58,13 +58,7 @@ def run_cc(self, ccinterface): return ccresult def tearDown(self): - # remove temporary nifti files - try: - os.remove(self.functionalnii) - os.remove(self.components_file) - os.remove(self.masknii) - except (OSError, TypeError) as e: - print(e) + utils.remove_nii(self.filenames.values()) def fake_noise_fun(self, i, j, l, m): return m*i + l - j diff --git a/nipype/algorithms/tests/test_tsnr.py b/nipype/algorithms/tests/test_tsnr.py index d90f720ebb..6448d70d83 100644 --- a/nipype/algorithms/tests/test_tsnr.py +++ b/nipype/algorithms/tests/test_tsnr.py @@ -16,11 +16,20 @@ class TestTSNR(unittest.TestCase): in_file = InputMultiPath(File(exists=True), mandatory=True, regress_poly = traits.Range(low=1, desc='Remove polynomials') ''' - in_file_name = 'tsnrinfile.nii' + + filenames = { + 'in_file': 'tsnrinfile.nii', + + # default output file names + 'detrended_file': '/home/ubuntu/nipype/detrend.nii.gz', + 'mean_file': '/home/ubuntu/nipype/mean.nii.gz', + 'stddev_file': '/home/ubuntu/nipype/stdev.nii.gz', + 'tsnr_file': '/home/ubuntu/nipype/tsnr.nii.gz' + } def setUp(self): # setup - utils.save_toy_nii(self.fake_data, self.in_file_name) + utils.save_toy_nii(self.fake_data, self.filenames['in_file']) def test_tsnr(self): # setup @@ -30,11 +39,7 @@ def test_tsnr(self): # cleanup def tearDown(self): - # remove temporary nifti files - try: - os.remove(self.in_file_name) - except (OSError, TypeError) as e: - print(e) + utils.remove_nii(self.filenames.values()) fake_data = np.array([[[[2, 4, 3, 9, 1], [3, 6, 4, 7, 4]], diff --git a/nipype/testing/utils.py b/nipype/testing/utils.py index 40c2d21b27..423516d452 100644 --- a/nipype/testing/utils.py +++ b/nipype/testing/utils.py @@ -106,3 +106,11 @@ def save_toy_nii(ndarray, filename): toy = nb.Nifti1Image(ndarray, np.eye(4)) nb.nifti1.save(toy, filename) return filename + +def remove_nii(filenames): + ''' remove temporary nifti files''' + for filename in filenames: + try: + os.remove(filename) + except (OSError, TypeError) as e: + print(e) From 87108d9e6dda39c0b3447bbe9ac674c2de458aba Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Thu, 8 Sep 2016 04:24:26 +0000 Subject: [PATCH 22/61] pre-refactor tests --- nipype/algorithms/tests/test_tsnr.py | 78 +++++++++++++++++++++++++--- 1 file changed, 72 insertions(+), 6 deletions(-) diff --git a/nipype/algorithms/tests/test_tsnr.py b/nipype/algorithms/tests/test_tsnr.py index 6448d70d83..57e2a5ffc5 100644 --- a/nipype/algorithms/tests/test_tsnr.py +++ b/nipype/algorithms/tests/test_tsnr.py @@ -4,6 +4,7 @@ from nipype.testing import assert_equal, assert_true, skipif, utils from nipype.algorithms.misc import TSNR +from hashlib import sha1 import unittest import mock import nibabel as nb @@ -17,9 +18,11 @@ class TestTSNR(unittest.TestCase): regress_poly = traits.Range(low=1, desc='Remove polynomials') ''' - filenames = { + in_filenames = { 'in_file': 'tsnrinfile.nii', + } + out_filenames = { # default output file names 'detrended_file': '/home/ubuntu/nipype/detrend.nii.gz', 'mean_file': '/home/ubuntu/nipype/mean.nii.gz', @@ -29,17 +32,80 @@ class TestTSNR(unittest.TestCase): def setUp(self): # setup - utils.save_toy_nii(self.fake_data, self.filenames['in_file']) + utils.save_toy_nii(self.fake_data, self.in_filenames['in_file']) def test_tsnr(self): - # setup # run - tsnrresult = TSNR(in_file=self.in_file_name, regress_poly=1) + tsnrresult = TSNR(in_file=self.in_filenames['in_file']).run() + + # assert + self.assert_expected_outputs(tsnrresult, { + 'mean_file': '1a55bcdf49901f25a2a838c90769989b9e4f2f19', + 'stddev_file': '0ba52a51fae90a9db6090c735432df5b742d663a', + 'tsnr_file': 'a794fc55766c8ad725437d3ff8b1153bd5f6e9b0' + }) + + def test_tsnr_withpoly1(self): + # run + tsnrresult = TSNR(in_file=self.in_filenames['in_file'], + regress_poly=1).run() + # assert - # cleanup + self.assert_expected_outputs_poly(tsnrresult, { + 'detrended_file': 'ee4f6c0b0e8c547617fc11aa50cf659436f9ccf0', + 'mean_file': '1a55bcdf49901f25a2a838c90769989b9e4f2f19', + 'stddev_file': 'e61d94e3cfea20b0c86c81bfdf80d82c55e9203b', + 'tsnr_file': 'a49f1cbd88f2aa71183dcd7aa4b86b17df622f0c' + }) + + def test_tsnr_withpoly2(self): + # run + tsnrresult = TSNR(in_file=self.in_filenames['in_file'], + regress_poly=2).run() + + # assert + self.assert_expected_outputs_poly(tsnrresult, { + 'detrended_file': '22cb7f058d61cc090eb1a9dd7d31550bd4362a61', + 'mean_file': '4cee6776461f6bc238d11a55c0a8d1947a5732df', + 'stddev_file': '7267de4d9b63fcc553115c0198f1fa3bbb6a5a13', + 'tsnr_file': '1c6ed05f94806838f7b563df65900f37e60f8a6d' + }) + + def test_tsnr_withpoly3(self): + # run + tsnrresult = TSNR(in_file=self.in_filenames['in_file'], + regress_poly=3).run() + + # assert + self.assert_expected_outputs_poly(tsnrresult, { + 'detrended_file': '3f2c1c7da233f128a7020b6fed79d6be2ec59fca', + 'mean_file': '4cee6776461f6bc238d11a55c0a8d1947a5732df', + 'stddev_file': '82bb793b76fab503d1d6b3e2d1b20a1bdebd7a2a', + 'tsnr_file': 'e004bd6096a0077b15058aabd4b0339bf6e21f64' + }) + + def assert_expected_outputs_poly(self, tsnrresult, hash_dict): + assert_equal(tsnrresult.outputs.detrended_file, + self.out_filenames['detrended_file']) + self.assert_expected_outputs(tsnrresult, hash_dict) + + def assert_expected_outputs(self, tsnrresult, hash_dict): + self.assert_default_outputs(tsnrresult.outputs) + self.assert_unchanged(hash_dict) + + def assert_default_outputs(self, outputs): + assert_equal(outputs.mean_file, self.out_filenames['mean_file']) + assert_equal(outputs.stddev_file, self.out_filenames['stddev_file']) + assert_equal(outputs.tsnr_file, self.out_filenames['tsnr_file']) + + def assert_unchanged(self, expected_hashes): + for key, hexhash in expected_hashes.iteritems(): + data = np.asanyarray(nb.load(self.out_filenames[key])._data) + assert_equal(sha1(str(data)).hexdigest(), hexhash) def tearDown(self): - utils.remove_nii(self.filenames.values()) + utils.remove_nii(self.in_filenames.values()) + utils.remove_nii(self.out_filenames.values()) fake_data = np.array([[[[2, 4, 3, 9, 1], [3, 6, 4, 7, 4]], From 33d2f309baf74303cc5e8f33061707f444a61c2b Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Thu, 8 Sep 2016 04:51:12 +0000 Subject: [PATCH 23/61] move tsnr to own file --- nipype/algorithms/misc.py | 91 ++------------------------------- nipype/algorithms/tsnr.py | 105 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 109 insertions(+), 87 deletions(-) create mode 100644 nipype/algorithms/tsnr.py diff --git a/nipype/algorithms/misc.py b/nipype/algorithms/misc.py index c476ed003a..911340eb97 100644 --- a/nipype/algorithms/misc.py +++ b/nipype/algorithms/misc.py @@ -39,6 +39,10 @@ BaseInterfaceInputSpec, isdefined, DynamicTraitedSpec, Undefined) from nipype.utils.filemanip import fname_presuffix, split_filename + +# for backwards compatibility +from nipype.algorithms.tsnr import TSNR + iflogger = logging.getLogger('interface') @@ -259,93 +263,6 @@ def _list_outputs(self): outputs['nifti_file'] = self._gen_output_file_name() return outputs - -class TSNRInputSpec(BaseInterfaceInputSpec): - in_file = InputMultiPath(File(exists=True), mandatory=True, - desc='realigned 4D file or a list of 3D files') - regress_poly = traits.Range(low=1, desc='Remove polynomials') - tsnr_file = File('tsnr.nii.gz', usedefault=True, hash_files=False, - desc='output tSNR file') - mean_file = File('mean.nii.gz', usedefault=True, hash_files=False, - desc='output mean file') - stddev_file = File('stdev.nii.gz', usedefault=True, hash_files=False, - desc='output tSNR file') - detrended_file = File('detrend.nii.gz', usedefault=True, hash_files=False, - desc='input file after detrending') - - -class TSNROutputSpec(TraitedSpec): - tsnr_file = File(exists=True, desc='tsnr image file') - mean_file = File(exists=True, desc='mean image file') - stddev_file = File(exists=True, desc='std dev image file') - detrended_file = File(desc='detrended input file') - - -class TSNR(BaseInterface): - """Computes the time-course SNR for a time series - - Typically you want to run this on a realigned time-series. - - Example - ------- - - >>> tsnr = TSNR() - >>> tsnr.inputs.in_file = 'functional.nii' - >>> res = tsnr.run() # doctest: +SKIP - - """ - input_spec = TSNRInputSpec - output_spec = TSNROutputSpec - - def _run_interface(self, runtime): - img = nb.load(self.inputs.in_file[0]) - header = img.header.copy() - vollist = [nb.load(filename) for filename in self.inputs.in_file] - data = np.concatenate([vol.get_data().reshape( - vol.get_shape()[:3] + (-1,)) for vol in vollist], axis=3) - data = np.nan_to_num(data) - - if data.dtype.kind == 'i': - header.set_data_dtype(np.float32) - data = data.astype(np.float32) - - if isdefined(self.inputs.regress_poly): - timepoints = img.shape[-1] - X = np.ones((timepoints, 1)) - for i in range(self.inputs.regress_poly): - X = np.hstack((X, legendre( - i + 1)(np.linspace(-1, 1, timepoints))[:, None])) - betas = np.dot(np.linalg.pinv(X), np.rollaxis(data, 3, 2)) - datahat = np.rollaxis(np.dot(X[:, 1:], - np.rollaxis( - betas[1:, :, :, :], 0, 3)), - 0, 4) - data = data - datahat - img = nb.Nifti1Image(data, img.get_affine(), header) - nb.save(img, op.abspath(self.inputs.detrended_file)) - - meanimg = np.mean(data, axis=3) - stddevimg = np.std(data, axis=3) - tsnr = np.zeros_like(meanimg) - tsnr[stddevimg > 1.e-3] = meanimg[stddevimg > 1.e-3] / stddevimg[stddevimg > 1.e-3] - img = nb.Nifti1Image(tsnr, img.get_affine(), header) - nb.save(img, op.abspath(self.inputs.tsnr_file)) - img = nb.Nifti1Image(meanimg, img.get_affine(), header) - nb.save(img, op.abspath(self.inputs.mean_file)) - img = nb.Nifti1Image(stddevimg, img.get_affine(), header) - nb.save(img, op.abspath(self.inputs.stddev_file)) - return runtime - - def _list_outputs(self): - outputs = self._outputs().get() - for k in ['tsnr_file', 'mean_file', 'stddev_file']: - outputs[k] = op.abspath(getattr(self.inputs, k)) - - if isdefined(self.inputs.regress_poly): - outputs['detrended_file'] = op.abspath(self.inputs.detrended_file) - return outputs - - class GunzipInputSpec(BaseInterfaceInputSpec): in_file = File(exists=True, mandatory=True) diff --git a/nipype/algorithms/tsnr.py b/nipype/algorithms/tsnr.py new file mode 100644 index 0000000000..b0cd3cb38c --- /dev/null +++ b/nipype/algorithms/tsnr.py @@ -0,0 +1,105 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +''' +Miscellaneous algorithms + + 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 nipype.interfaces.base import (BaseInterface, traits, TraitedSpec, File, + InputMultiPath, BaseInterfaceInputSpec, + isdefined) + +import nibabel as nb +import numpy as np +import os.path as op +from scipy.special import legendre + +class TSNRInputSpec(BaseInterfaceInputSpec): + in_file = InputMultiPath(File(exists=True), mandatory=True, + desc='realigned 4D file or a list of 3D files') + regress_poly = traits.Range(low=1, desc='Remove polynomials') + tsnr_file = File('tsnr.nii.gz', usedefault=True, hash_files=False, + desc='output tSNR file') + mean_file = File('mean.nii.gz', usedefault=True, hash_files=False, + desc='output mean file') + stddev_file = File('stdev.nii.gz', usedefault=True, hash_files=False, + desc='output tSNR file') + detrended_file = File('detrend.nii.gz', usedefault=True, hash_files=False, + desc='input file after detrending') + + +class TSNROutputSpec(TraitedSpec): + tsnr_file = File(exists=True, desc='tsnr image file') + mean_file = File(exists=True, desc='mean image file') + stddev_file = File(exists=True, desc='std dev image file') + detrended_file = File(desc='detrended input file') + + +class TSNR(BaseInterface): + """Computes the time-course SNR for a time series + + Typically you want to run this on a realigned time-series. + + Example + ------- + + >>> tsnr = TSNR() + >>> tsnr.inputs.in_file = 'functional.nii' + >>> res = tsnr.run() # doctest: +SKIP + + """ + input_spec = TSNRInputSpec + output_spec = TSNROutputSpec + + def _run_interface(self, runtime): + img = nb.load(self.inputs.in_file[0]) + header = img.header.copy() + vollist = [nb.load(filename) for filename in self.inputs.in_file] + data = np.concatenate([vol.get_data().reshape( + vol.get_shape()[:3] + (-1,)) for vol in vollist], axis=3) + data = np.nan_to_num(data) + + if data.dtype.kind == 'i': + header.set_data_dtype(np.float32) + data = data.astype(np.float32) + + if isdefined(self.inputs.regress_poly): + timepoints = img.shape[-1] + X = np.ones((timepoints, 1)) + for i in range(self.inputs.regress_poly): + X = np.hstack((X, legendre( + i + 1)(np.linspace(-1, 1, timepoints))[:, None])) + betas = np.dot(np.linalg.pinv(X), np.rollaxis(data, 3, 2)) + datahat = np.rollaxis(np.dot(X[:, 1:], + np.rollaxis( + betas[1:, :, :, :], 0, 3)), + 0, 4) + data = data - datahat + img = nb.Nifti1Image(data, img.get_affine(), header) + nb.save(img, op.abspath(self.inputs.detrended_file)) + + meanimg = np.mean(data, axis=3) + stddevimg = np.std(data, axis=3) + tsnr = np.zeros_like(meanimg) + tsnr[stddevimg > 1.e-3] = meanimg[stddevimg > 1.e-3] / stddevimg[stddevimg > 1.e-3] + img = nb.Nifti1Image(tsnr, img.get_affine(), header) + nb.save(img, op.abspath(self.inputs.tsnr_file)) + img = nb.Nifti1Image(meanimg, img.get_affine(), header) + nb.save(img, op.abspath(self.inputs.mean_file)) + img = nb.Nifti1Image(stddevimg, img.get_affine(), header) + nb.save(img, op.abspath(self.inputs.stddev_file)) + return runtime + + def _list_outputs(self): + outputs = self._outputs().get() + for k in ['tsnr_file', 'mean_file', 'stddev_file']: + outputs[k] = op.abspath(getattr(self.inputs, k)) + + if isdefined(self.inputs.regress_poly): + outputs['detrended_file'] = op.abspath(self.inputs.detrended_file) + return outputs From 1511beac8b7220e4cbc8e0d29d94ea3b2b071653 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Thu, 8 Sep 2016 05:13:03 +0000 Subject: [PATCH 24/61] naive refactor --- nipype/algorithms/tsnr.py | 31 ++++++++++++++++++------------- 1 file changed, 18 insertions(+), 13 deletions(-) diff --git a/nipype/algorithms/tsnr.py b/nipype/algorithms/tsnr.py index b0cd3cb38c..b3d9be8385 100644 --- a/nipype/algorithms/tsnr.py +++ b/nipype/algorithms/tsnr.py @@ -69,19 +69,8 @@ def _run_interface(self, runtime): data = data.astype(np.float32) if isdefined(self.inputs.regress_poly): - timepoints = img.shape[-1] - X = np.ones((timepoints, 1)) - for i in range(self.inputs.regress_poly): - X = np.hstack((X, legendre( - i + 1)(np.linspace(-1, 1, timepoints))[:, None])) - betas = np.dot(np.linalg.pinv(X), np.rollaxis(data, 3, 2)) - datahat = np.rollaxis(np.dot(X[:, 1:], - np.rollaxis( - betas[1:, :, :, :], 0, 3)), - 0, 4) - data = data - datahat - img = nb.Nifti1Image(data, img.get_affine(), header) - nb.save(img, op.abspath(self.inputs.detrended_file)) + img, data = regress_poly(self.inputs.regress_poly, header, img, data, + self.inputs.detrended_file) meanimg = np.mean(data, axis=3) stddevimg = np.std(data, axis=3) @@ -103,3 +92,19 @@ def _list_outputs(self): if isdefined(self.inputs.regress_poly): outputs['detrended_file'] = op.abspath(self.inputs.detrended_file) return outputs + +def regress_poly(degree, header, img, data, filename): + timepoints = img.shape[-1] + X = np.ones((timepoints, 1)) + for i in range(degree): + X = np.hstack((X, legendre( + i + 1)(np.linspace(-1, 1, timepoints))[:, None])) + betas = np.dot(np.linalg.pinv(X), np.rollaxis(data, 3, 2)) + datahat = np.rollaxis(np.dot(X[:, 1:], + np.rollaxis( + betas[1:, :, :, :], 0, 3)), + 0, 4) + regressed_data = data - datahat + regressed_img = nb.Nifti1Image(regressed_data, img.get_affine(), header) + nb.save(regressed_img, op.abspath(filename)) + return regressed_img, regressed_data From b66db5ff5c158abf2ef5889a8aec18a72c90eb1e Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Thu, 8 Sep 2016 06:03:57 +0000 Subject: [PATCH 25/61] make compcor tests more sensitive to changes, simplify regress_poly() --- nipype/algorithms/compcor.py | 8 ++++++-- nipype/algorithms/tests/test_compcor.py | 22 +++++++++++++++++----- nipype/algorithms/tsnr.py | 13 ++++++------- 3 files changed, 29 insertions(+), 14 deletions(-) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index c2eaced2c4..cdbcf1b982 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -3,6 +3,7 @@ from nipype.pipeline import engine as pe from nipype.interfaces.utility import IdentityInterface +from nipype.algorithms.tsnr import regress_poly import nibabel as nb import numpy as np @@ -14,12 +15,15 @@ class CompCorInputSpec(BaseInterfaceInputSpec): desc='already realigned brain image (4D)') mask_file = File(exists=True, mandatory=False, desc='mask file that determines ROI (3D)') - components_file = File('components_file.txt', exists=False, mandatory=False, usedefault=True, + components_file = File('components_file.txt', exists=False, mandatory=False, + usedefault=True, desc='filename to store physiological components in') num_components = traits.Int(6, usedefault=True) # 6 for BOLD, 4 for ASL + regress_poly = traits.Range(low=1, default=1, usedefault=True) class CompCorOutputSpec(TraitedSpec): - components_file = File(desc='text file containing the noise components', exists=True) + components_file = File(exists=True, + desc='text file containing the noise components') class CompCor(BaseInterface): ''' diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index ffb5b198b5..8dba5ece7c 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -10,6 +10,7 @@ import nibabel as nb import numpy as np import os +from hashlib import sha1 class TestCompCor(unittest.TestCase): ''' Note: Tests currently do a poor job of testing functionality ''' @@ -30,21 +31,25 @@ def test_compcor(self): mask[0,0,1] = 0 mask_file = utils.save_toy_nii(mask, self.filenames['masknii']) + expected_sha1 = 'b0dd7f9ab7ba8f516712eb0204dacc9e397fc6aa' + ccresult = self.run_cc(CompCor(realigned_file=self.realigned_file, - mask_file=mask_file)) + mask_file=mask_file), + expected_sha1) accresult = self.run_cc(ACompCor(realigned_file=self.realigned_file, - mask_file=mask_file, - components_file='acc_components_file')) + mask_file=mask_file, + components_file='acc_components_file'), + expected_sha1) assert_equal(os.path.getsize(ccresult.outputs.components_file), os.path.getsize(accresult.outputs.components_file)) def test_tcompcor(self): ccinterface = TCompCor(realigned_file=self.realigned_file) - self.run_cc(ccinterface) + self.run_cc(ccinterface, '12e54c07281a28ac0da3b934dce5c9d27626848a') - def run_cc(self, ccinterface): + def run_cc(self, ccinterface, expected_components_data_sha1): # run ccresult = ccinterface.run() @@ -55,6 +60,13 @@ def run_cc(self, ccinterface): assert_true(os.path.getsize(expected_file) > 0) assert_equal(ccinterface.inputs.num_components, 6) + with open(ccresult.outputs.components_file, 'r') as components_file: + components_data = [line 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]) + print(str(components_data), "comdata") + assert_equal(sha1(str(components_data)).hexdigest(), expected_components_data_sha1) return ccresult def tearDown(self): diff --git a/nipype/algorithms/tsnr.py b/nipype/algorithms/tsnr.py index b3d9be8385..fbc39c6d3e 100644 --- a/nipype/algorithms/tsnr.py +++ b/nipype/algorithms/tsnr.py @@ -69,8 +69,9 @@ def _run_interface(self, runtime): data = data.astype(np.float32) if isdefined(self.inputs.regress_poly): - img, data = regress_poly(self.inputs.regress_poly, header, img, data, - self.inputs.detrended_file) + data = regress_poly(self.inputs.regress_poly, data) + img = nb.Nifti1Image(data, img.get_affine(), header) + nb.save(img, op.abspath(self.inputs.detrended_file)) meanimg = np.mean(data, axis=3) stddevimg = np.std(data, axis=3) @@ -93,8 +94,8 @@ def _list_outputs(self): outputs['detrended_file'] = op.abspath(self.inputs.detrended_file) return outputs -def regress_poly(degree, header, img, data, filename): - timepoints = img.shape[-1] +def regress_poly(degree, data): + timepoints = data.shape[-1] X = np.ones((timepoints, 1)) for i in range(degree): X = np.hstack((X, legendre( @@ -105,6 +106,4 @@ def regress_poly(degree, header, img, data, filename): betas[1:, :, :, :], 0, 3)), 0, 4) regressed_data = data - datahat - regressed_img = nb.Nifti1Image(regressed_data, img.get_affine(), header) - nb.save(regressed_img, op.abspath(filename)) - return regressed_img, regressed_data + return regressed_data From 4b9dbe495bfacae7dd4951bfc883758225eb677a Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Thu, 8 Sep 2016 19:26:04 +0000 Subject: [PATCH 26/61] fix test oops --- nipype/algorithms/tests/test_tsnr.py | 26 +++++++++++++++----------- 1 file changed, 15 insertions(+), 11 deletions(-) diff --git a/nipype/algorithms/tests/test_tsnr.py b/nipype/algorithms/tests/test_tsnr.py index 57e2a5ffc5..0cdfe8b5d0 100644 --- a/nipype/algorithms/tests/test_tsnr.py +++ b/nipype/algorithms/tests/test_tsnr.py @@ -22,13 +22,14 @@ class TestTSNR(unittest.TestCase): 'in_file': 'tsnrinfile.nii', } - out_filenames = { - # default output file names - 'detrended_file': '/home/ubuntu/nipype/detrend.nii.gz', - 'mean_file': '/home/ubuntu/nipype/mean.nii.gz', - 'stddev_file': '/home/ubuntu/nipype/stdev.nii.gz', - 'tsnr_file': '/home/ubuntu/nipype/tsnr.nii.gz' - } + out_filenames = {key: os.path.abspath(filename) for key, filename in + { + # default output file names + 'detrended_file': 'detrend.nii.gz', + 'mean_file': 'mean.nii.gz', + 'stddev_file': 'stdev.nii.gz', + 'tsnr_file': 'tsnr.nii.gz' + }.items()} def setUp(self): # setup @@ -85,7 +86,7 @@ def test_tsnr_withpoly3(self): }) def assert_expected_outputs_poly(self, tsnrresult, hash_dict): - assert_equal(tsnrresult.outputs.detrended_file, + assert_equal(os.path.abspath(tsnrresult.outputs.detrended_file), self.out_filenames['detrended_file']) self.assert_expected_outputs(tsnrresult, hash_dict) @@ -94,9 +95,12 @@ def assert_expected_outputs(self, tsnrresult, hash_dict): self.assert_unchanged(hash_dict) def assert_default_outputs(self, outputs): - assert_equal(outputs.mean_file, self.out_filenames['mean_file']) - assert_equal(outputs.stddev_file, self.out_filenames['stddev_file']) - assert_equal(outputs.tsnr_file, self.out_filenames['tsnr_file']) + assert_equal(os.path.abspath(outputs.mean_file), + self.out_filenames['mean_file']) + assert_equal(os.path.abspath(outputs.stddev_file), + self.out_filenames['stddev_file']) + assert_equal(os.path.abspath(outputs.tsnr_file), + self.out_filenames['tsnr_file']) def assert_unchanged(self, expected_hashes): for key, hexhash in expected_hashes.iteritems(): From bd9113a41582e85e50573ecb5221a6bff0727b4e Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Thu, 8 Sep 2016 19:44:38 +0000 Subject: [PATCH 27/61] m rm unnecessary/confusing test msgs --- nipype/testing/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nipype/testing/utils.py b/nipype/testing/utils.py index 423516d452..c99f3a2f61 100644 --- a/nipype/testing/utils.py +++ b/nipype/testing/utils.py @@ -113,4 +113,4 @@ def remove_nii(filenames): try: os.remove(filename) except (OSError, TypeError) as e: - print(e) + pass From e78b0ae025933e3b540be57be4a429ce1a79810b Mon Sep 17 00:00:00 2001 From: oesteban Date: Thu, 8 Sep 2016 16:42:13 -0700 Subject: [PATCH 28/61] make regress_poly take any kind of input data --- nipype/algorithms/tsnr.py | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/nipype/algorithms/tsnr.py b/nipype/algorithms/tsnr.py index fbc39c6d3e..762ff00ca1 100644 --- a/nipype/algorithms/tsnr.py +++ b/nipype/algorithms/tsnr.py @@ -95,15 +95,24 @@ def _list_outputs(self): return outputs def regress_poly(degree, data): - timepoints = data.shape[-1] + datashape = data.shape + timepoints = datashape[-1] + + # Rearrange all voxel-wise time-series in rows + data = data.reshape((-1, timepoints)) + + # Generate design matrix X = np.ones((timepoints, 1)) for i in range(degree): X = np.hstack((X, legendre( i + 1)(np.linspace(-1, 1, timepoints))[:, None])) - betas = np.dot(np.linalg.pinv(X), np.rollaxis(data, 3, 2)) - datahat = np.rollaxis(np.dot(X[:, 1:], - np.rollaxis( - betas[1:, :, :, :], 0, 3)), - 0, 4) + + # Calculate coefficients + betas = np.linalg.pinv(X).dot(data.T) + + # Estimation + datahat = X[:, 1:].dot(betas[1:, ...]).T regressed_data = data - datahat - return regressed_data + + # Back to original shape + return regressed_data.reshape(datashape) From 5bf2d540811035369a131422128f01de22f8b518 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Fri, 9 Sep 2016 00:30:31 +0000 Subject: [PATCH 29/61] improve code readability --- nipype/algorithms/tsnr.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/nipype/algorithms/tsnr.py b/nipype/algorithms/tsnr.py index fbc39c6d3e..c2644c5de0 100644 --- a/nipype/algorithms/tsnr.py +++ b/nipype/algorithms/tsnr.py @@ -13,7 +13,7 @@ from nipype.interfaces.base import (BaseInterface, traits, TraitedSpec, File, InputMultiPath, BaseInterfaceInputSpec, isdefined) - +from numpy import newaxis import nibabel as nb import numpy as np import os.path as op @@ -95,11 +95,16 @@ def _list_outputs(self): return outputs def regress_poly(degree, data): + ''' returns data with degree polynomial regressed out. + The last dimension (i.e. data.shape[-1]) should be time. + ''' timepoints = data.shape[-1] X = np.ones((timepoints, 1)) for i in range(degree): - X = np.hstack((X, legendre( - i + 1)(np.linspace(-1, 1, timepoints))[:, None])) + polynomial_func = legendre(i+1) + value_array = np.linspace(-1, 1, timepoints) + X = np.hstack((X, polynomial_func(value_array)[:, newaxis])) + betas = np.dot(np.linalg.pinv(X), np.rollaxis(data, 3, 2)) datahat = np.rollaxis(np.dot(X[:, 1:], np.rollaxis( From fd45c05d09997841baccdc0ff6980be0cf5108cf Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Fri, 9 Sep 2016 00:47:43 +0000 Subject: [PATCH 30/61] fixed for real this time? --- nipype/algorithms/tests/test_tsnr.py | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/nipype/algorithms/tests/test_tsnr.py b/nipype/algorithms/tests/test_tsnr.py index 0cdfe8b5d0..512edda560 100644 --- a/nipype/algorithms/tests/test_tsnr.py +++ b/nipype/algorithms/tests/test_tsnr.py @@ -22,14 +22,12 @@ class TestTSNR(unittest.TestCase): 'in_file': 'tsnrinfile.nii', } - out_filenames = {key: os.path.abspath(filename) for key, filename in - { - # default output file names - 'detrended_file': 'detrend.nii.gz', - 'mean_file': 'mean.nii.gz', - 'stddev_file': 'stdev.nii.gz', - 'tsnr_file': 'tsnr.nii.gz' - }.items()} + out_filenames = {# default output file names + 'detrended_file': 'detrend.nii.gz', + 'mean_file': 'mean.nii.gz', + 'stddev_file': 'stdev.nii.gz', + 'tsnr_file': 'tsnr.nii.gz' + } def setUp(self): # setup @@ -86,7 +84,7 @@ def test_tsnr_withpoly3(self): }) def assert_expected_outputs_poly(self, tsnrresult, hash_dict): - assert_equal(os.path.abspath(tsnrresult.outputs.detrended_file), + assert_equal(os.path.basename(tsnrresult.outputs.detrended_file), self.out_filenames['detrended_file']) self.assert_expected_outputs(tsnrresult, hash_dict) @@ -95,11 +93,11 @@ def assert_expected_outputs(self, tsnrresult, hash_dict): self.assert_unchanged(hash_dict) def assert_default_outputs(self, outputs): - assert_equal(os.path.abspath(outputs.mean_file), + assert_equal(os.path.basename(outputs.mean_file), self.out_filenames['mean_file']) - assert_equal(os.path.abspath(outputs.stddev_file), + assert_equal(os.path.basename(outputs.stddev_file), self.out_filenames['stddev_file']) - assert_equal(os.path.abspath(outputs.tsnr_file), + assert_equal(os.path.basename(outputs.tsnr_file), self.out_filenames['tsnr_file']) def assert_unchanged(self, expected_hashes): From 62bb3bd06250d1a410348d6c1f9ebc5661b253ae Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Fri, 9 Sep 2016 02:05:27 +0000 Subject: [PATCH 31/61] change absolute imports to explicit relative imports --- nipype/algorithms/compcor.py | 9 ++++----- nipype/algorithms/tests/test_compcor.py | 5 ++--- nipype/algorithms/tests/test_tsnr.py | 4 ++-- nipype/algorithms/tsnr.py | 6 +++--- 4 files changed, 11 insertions(+), 13 deletions(-) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index cdbcf1b982..2fce73d17f 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -1,9 +1,8 @@ from ..interfaces.base import (BaseInterfaceInputSpec, TraitedSpec, BaseInterface, traits, File) - -from nipype.pipeline import engine as pe -from nipype.interfaces.utility import IdentityInterface -from nipype.algorithms.tsnr import regress_poly +from ..pipeline import engine as pe +from ..interfaces.utility import IdentityInterface +from .tsnr import regress_poly import nibabel as nb import numpy as np @@ -16,7 +15,7 @@ class CompCorInputSpec(BaseInterfaceInputSpec): mask_file = File(exists=True, mandatory=False, desc='mask file that determines ROI (3D)') components_file = File('components_file.txt', exists=False, mandatory=False, - usedefault=True, + usedefault=True, desc='filename to store physiological components in') num_components = traits.Int(6, usedefault=True) # 6 for BOLD, 4 for ASL regress_poly = traits.Range(low=1, default=1, usedefault=True) diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index 8dba5ece7c..b112a5f697 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -1,9 +1,8 @@ # 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, - utils) -from nipype.algorithms.compcor import CompCor, TCompCor, ACompCor +from ...testing import assert_equal, assert_true, assert_false, skipif, utils +from ..compcor import CompCor, TCompCor, ACompCor import unittest import mock diff --git a/nipype/algorithms/tests/test_tsnr.py b/nipype/algorithms/tests/test_tsnr.py index 512edda560..c3a2f55e5d 100644 --- a/nipype/algorithms/tests/test_tsnr.py +++ b/nipype/algorithms/tests/test_tsnr.py @@ -1,8 +1,8 @@ # 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, skipif, utils -from nipype.algorithms.misc import TSNR +from ...testing import assert_equal, assert_true, skipif, utils +from ..misc import TSNR from hashlib import sha1 import unittest diff --git a/nipype/algorithms/tsnr.py b/nipype/algorithms/tsnr.py index c2644c5de0..ee188495a2 100644 --- a/nipype/algorithms/tsnr.py +++ b/nipype/algorithms/tsnr.py @@ -10,9 +10,9 @@ >>> os.chdir(datadir) ''' -from nipype.interfaces.base import (BaseInterface, traits, TraitedSpec, File, - InputMultiPath, BaseInterfaceInputSpec, - isdefined) +from ..interfaces.base import (BaseInterface, traits, TraitedSpec, File, + InputMultiPath, BaseInterfaceInputSpec, + isdefined) from numpy import newaxis import nibabel as nb import numpy as np From da1e59a068518625b95521844309f42d0c84a307 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Fri, 9 Sep 2016 02:19:54 +0000 Subject: [PATCH 32/61] add use_regress_poly Boolean to compcor --- nipype/algorithms/compcor.py | 23 +++++++++++++++++------ 1 file changed, 17 insertions(+), 6 deletions(-) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index 2fce73d17f..5e60e6a35a 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -14,11 +14,15 @@ class CompCorInputSpec(BaseInterfaceInputSpec): desc='already realigned brain image (4D)') mask_file = File(exists=True, mandatory=False, desc='mask file that determines ROI (3D)') - components_file = File('components_file.txt', exists=False, mandatory=False, - usedefault=True, - desc='filename to store physiological components in') + components_file = File('components_file.txt', exists=False, + mandatory=False, usedefault=True, + desc='filename to store physiological components') num_components = traits.Int(6, usedefault=True) # 6 for BOLD, 4 for ASL - regress_poly = traits.Range(low=1, default=1, usedefault=True) + use_regress_poly = traits.Bool(True, usedefault=True, + desc='use polynomial regression' + 'pre-component extraction') + regress_poly_degree = traits.Range(low=1, default=1, usedefault=True, + desc='the degree polynomial to use') class CompCorOutputSpec(TraitedSpec): components_file = File(exists=True, @@ -35,6 +39,8 @@ class CompCor(BaseInterface): >>> ccinterface.inputs.realigned_file = 'nipype/testing/data/functional.nii' >>> ccinterface.inputs.mask_file = 'nipype/testing/data/mask.nii' >>> ccinterface.inputs.num_components = 1 + >>> ccinterface.inputs.use_regress_poly = True + >>> ccinterface.inputs.regress_poly_degree = 2 ''' input_spec = CompCorInputSpec output_spec = CompCorOutputSpec @@ -51,16 +57,21 @@ def _run_interface(self, runtime): # 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 ..." + if self.inputs.use_regress_poly: + # "The constant and linear trends of the columns in the matrix M were removed ..." + regress_poly(self.inputs.regress_poly_degree, voxel_timecourses) + + ''' 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." M = M / self._compute_tSTD(M, 1.) From ee726a17a9fccf626904ee210d598a11a4afdafe Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Fri, 9 Sep 2016 03:36:00 +0000 Subject: [PATCH 33/61] undo accidental changes in merge --- Makefile | 5 +---- nipype/algorithms/tsnr.py | 7 ++----- 2 files changed, 3 insertions(+), 9 deletions(-) diff --git a/Makefile b/Makefile index e36e222b8a..8c7856e7ca 100644 --- a/Makefile +++ b/Makefile @@ -53,10 +53,7 @@ inplace: $(PYTHON) setup.py build_ext -i test-code: in -# $(NOSETESTS) -s nipype/algorithms/compcor.py --with-doctest -v - $(NOSETESTS) -s nipype/algorithms/tests/test_tsnr.py #--with-coverage - $(NOSETESTS) -s nipype/algorithms/tests/test_compcor.py #--with-coverage -# $(NOSETESTS) -s nipype/algorithms/misc.py --with-doctest -v + $(NOSETESTS) -s nipype --with-doctest test-doc: $(NOSETESTS) -s --with-doctest --doctest-tests --doctest-extension=rst \ diff --git a/nipype/algorithms/tsnr.py b/nipype/algorithms/tsnr.py index ff69e13dd6..8a69006759 100644 --- a/nipype/algorithms/tsnr.py +++ b/nipype/algorithms/tsnr.py @@ -60,11 +60,8 @@ def _run_interface(self, runtime): img = nb.load(self.inputs.in_file[0]) header = img.header.copy() vollist = [nb.load(filename) for filename in self.inputs.in_file] - data = np.concatenate( - [vol.get_data().reshape( - vol.get_shape()[:3] + (-1,)) - for vol in vollist], - axis=3) + data = np.concatenate([vol.get_data().reshape( + vol.get_shape()[:3] + (-1,)) for vol in vollist], axis=3) data = np.nan_to_num(data) if data.dtype.kind == 'i': From 0a5bc7a612a3fc942c50ffcf5d4f9b284744168e Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Fri, 9 Sep 2016 03:46:14 +0000 Subject: [PATCH 34/61] move tsnr back into misc --- nipype/algorithms/compcor.py | 2 +- nipype/algorithms/misc.py | 105 +++++++++++++++++++++++++++++- nipype/algorithms/tsnr.py | 122 ----------------------------------- 3 files changed, 103 insertions(+), 126 deletions(-) delete mode 100644 nipype/algorithms/tsnr.py diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index 4fc94f1c94..fac35d1e7d 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -2,7 +2,7 @@ BaseInterface, traits, File) from ..pipeline import engine as pe from ..interfaces.utility import IdentityInterface -from .tsnr import regress_poly +from .misc import regress_poly import nibabel as nb import numpy as np diff --git a/nipype/algorithms/misc.py b/nipype/algorithms/misc.py index 911340eb97..1c3f8c9354 100644 --- a/nipype/algorithms/misc.py +++ b/nipype/algorithms/misc.py @@ -40,9 +40,6 @@ DynamicTraitedSpec, Undefined) from nipype.utils.filemanip import fname_presuffix, split_filename -# for backwards compatibility -from nipype.algorithms.tsnr import TSNR - iflogger = logging.getLogger('interface') @@ -263,6 +260,108 @@ def _list_outputs(self): outputs['nifti_file'] = self._gen_output_file_name() return outputs +class TSNRInputSpec(BaseInterfaceInputSpec): + in_file = InputMultiPath(File(exists=True), mandatory=True, + desc='realigned 4D file or a list of 3D files') + regress_poly = traits.Range(low=1, desc='Remove polynomials') + tsnr_file = File('tsnr.nii.gz', usedefault=True, hash_files=False, + desc='output tSNR file') + mean_file = File('mean.nii.gz', usedefault=True, hash_files=False, + desc='output mean file') + stddev_file = File('stdev.nii.gz', usedefault=True, hash_files=False, + desc='output tSNR file') + detrended_file = File('detrend.nii.gz', usedefault=True, hash_files=False, + desc='input file after detrending') + + +class TSNROutputSpec(TraitedSpec): + tsnr_file = File(exists=True, desc='tsnr image file') + mean_file = File(exists=True, desc='mean image file') + stddev_file = File(exists=True, desc='std dev image file') + detrended_file = File(desc='detrended input file') + + +class TSNR(BaseInterface): + """Computes the time-course SNR for a time series + + Typically you want to run this on a realigned time-series. + + Example + ------- + + >>> tsnr = TSNR() + >>> tsnr.inputs.in_file = 'functional.nii' + >>> res = tsnr.run() # doctest: +SKIP + + """ + input_spec = TSNRInputSpec + output_spec = TSNROutputSpec + + def _run_interface(self, runtime): + img = nb.load(self.inputs.in_file[0]) + header = img.header.copy() + vollist = [nb.load(filename) for filename in self.inputs.in_file] + data = np.concatenate([vol.get_data().reshape( + vol.get_shape()[:3] + (-1,)) for vol in vollist], axis=3) + data = np.nan_to_num(data) + + if data.dtype.kind == 'i': + header.set_data_dtype(np.float32) + data = data.astype(np.float32) + + if isdefined(self.inputs.regress_poly): + data = regress_poly(self.inputs.regress_poly, data) + img = nb.Nifti1Image(data, img.get_affine(), header) + nb.save(img, op.abspath(self.inputs.detrended_file)) + + meanimg = np.mean(data, axis=3) + stddevimg = np.std(data, axis=3) + tsnr = np.zeros_like(meanimg) + tsnr[stddevimg > 1.e-3] = meanimg[stddevimg > 1.e-3] / stddevimg[stddevimg > 1.e-3] + img = nb.Nifti1Image(tsnr, img.get_affine(), header) + nb.save(img, op.abspath(self.inputs.tsnr_file)) + img = nb.Nifti1Image(meanimg, img.get_affine(), header) + nb.save(img, op.abspath(self.inputs.mean_file)) + img = nb.Nifti1Image(stddevimg, img.get_affine(), header) + nb.save(img, op.abspath(self.inputs.stddev_file)) + return runtime + + def _list_outputs(self): + outputs = self._outputs().get() + for k in ['tsnr_file', 'mean_file', 'stddev_file']: + outputs[k] = op.abspath(getattr(self.inputs, k)) + + if isdefined(self.inputs.regress_poly): + outputs['detrended_file'] = op.abspath(self.inputs.detrended_file) + return outputs + +def regress_poly(degree, data): + ''' returns data with degree polynomial regressed out. + The last dimension (i.e. data.shape[-1]) should be time. + ''' + datashape = data.shape + timepoints = datashape[-1] + + # Rearrange all voxel-wise time-series in rows + data = data.reshape((-1, timepoints)) + + # Generate design matrix + X = np.ones((timepoints, 1)) + for i in range(degree): + polynomial_func = legendre(i+1) + value_array = np.linspace(-1, 1, timepoints) + X = np.hstack((X, polynomial_func(value_array)[:, np.newaxis])) + + # Calculate coefficients + betas = np.linalg.pinv(X).dot(data.T) + + # Estimation + datahat = X[:, 1:].dot(betas[1:, ...]).T + regressed_data = data - datahat + + # Back to original shape + return regressed_data.reshape(datashape) + class GunzipInputSpec(BaseInterfaceInputSpec): in_file = File(exists=True, mandatory=True) diff --git a/nipype/algorithms/tsnr.py b/nipype/algorithms/tsnr.py deleted file mode 100644 index 8a69006759..0000000000 --- a/nipype/algorithms/tsnr.py +++ /dev/null @@ -1,122 +0,0 @@ -# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- -# vi: set ft=python sts=4 ts=4 sw=4 et: -''' -Miscellaneous algorithms - - 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 ..interfaces.base import (BaseInterface, traits, TraitedSpec, File, - InputMultiPath, BaseInterfaceInputSpec, - isdefined) -from numpy import newaxis -import nibabel as nb -import numpy as np -import os.path as op -from scipy.special import legendre - -class TSNRInputSpec(BaseInterfaceInputSpec): - in_file = InputMultiPath(File(exists=True), mandatory=True, - desc='realigned 4D file or a list of 3D files') - regress_poly = traits.Range(low=1, desc='Remove polynomials') - tsnr_file = File('tsnr.nii.gz', usedefault=True, hash_files=False, - desc='output tSNR file') - mean_file = File('mean.nii.gz', usedefault=True, hash_files=False, - desc='output mean file') - stddev_file = File('stdev.nii.gz', usedefault=True, hash_files=False, - desc='output tSNR file') - detrended_file = File('detrend.nii.gz', usedefault=True, hash_files=False, - desc='input file after detrending') - - -class TSNROutputSpec(TraitedSpec): - tsnr_file = File(exists=True, desc='tsnr image file') - mean_file = File(exists=True, desc='mean image file') - stddev_file = File(exists=True, desc='std dev image file') - detrended_file = File(desc='detrended input file') - - -class TSNR(BaseInterface): - """Computes the time-course SNR for a time series - - Typically you want to run this on a realigned time-series. - - Example - ------- - - >>> tsnr = TSNR() - >>> tsnr.inputs.in_file = 'functional.nii' - >>> res = tsnr.run() # doctest: +SKIP - - """ - input_spec = TSNRInputSpec - output_spec = TSNROutputSpec - - def _run_interface(self, runtime): - img = nb.load(self.inputs.in_file[0]) - header = img.header.copy() - vollist = [nb.load(filename) for filename in self.inputs.in_file] - data = np.concatenate([vol.get_data().reshape( - vol.get_shape()[:3] + (-1,)) for vol in vollist], axis=3) - data = np.nan_to_num(data) - - if data.dtype.kind == 'i': - header.set_data_dtype(np.float32) - data = data.astype(np.float32) - - if isdefined(self.inputs.regress_poly): - data = regress_poly(self.inputs.regress_poly, data) - img = nb.Nifti1Image(data, img.get_affine(), header) - nb.save(img, op.abspath(self.inputs.detrended_file)) - - meanimg = np.mean(data, axis=3) - stddevimg = np.std(data, axis=3) - tsnr = np.zeros_like(meanimg) - tsnr[stddevimg > 1.e-3] = meanimg[stddevimg > 1.e-3] / stddevimg[stddevimg > 1.e-3] - img = nb.Nifti1Image(tsnr, img.get_affine(), header) - nb.save(img, op.abspath(self.inputs.tsnr_file)) - img = nb.Nifti1Image(meanimg, img.get_affine(), header) - nb.save(img, op.abspath(self.inputs.mean_file)) - img = nb.Nifti1Image(stddevimg, img.get_affine(), header) - nb.save(img, op.abspath(self.inputs.stddev_file)) - return runtime - - def _list_outputs(self): - outputs = self._outputs().get() - for k in ['tsnr_file', 'mean_file', 'stddev_file']: - outputs[k] = op.abspath(getattr(self.inputs, k)) - - if isdefined(self.inputs.regress_poly): - outputs['detrended_file'] = op.abspath(self.inputs.detrended_file) - return outputs - -def regress_poly(degree, data): - ''' returns data with degree polynomial regressed out. - The last dimension (i.e. data.shape[-1]) should be time. - ''' - datashape = data.shape - timepoints = datashape[-1] - - # Rearrange all voxel-wise time-series in rows - data = data.reshape((-1, timepoints)) - - # Generate design matrix - X = np.ones((timepoints, 1)) - for i in range(degree): - polynomial_func = legendre(i+1) - value_array = np.linspace(-1, 1, timepoints) - X = np.hstack((X, polynomial_func(value_array)[:, newaxis])) - - # Calculate coefficients - betas = np.linalg.pinv(X).dot(data.T) - - # Estimation - datahat = X[:, 1:].dot(betas[1:, ...]).T - regressed_data = data - datahat - - # Back to original shape - return regressed_data.reshape(datashape) From 61c96ad1d683000382faf176385aca3c9a2dec4d Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Fri, 9 Sep 2016 03:47:57 +0000 Subject: [PATCH 35/61] remove rogue file --- nipype/testing/data/mask.nii | Bin 416 -> 0 bytes 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 nipype/testing/data/mask.nii diff --git a/nipype/testing/data/mask.nii b/nipype/testing/data/mask.nii deleted file mode 100644 index 940478111f2eeb923e22f549c6441eeb0c83a8ba..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 416 zcma!HWFQJKGcbW6BLf7YX<~5z3pCg>FyO-*oFVdPx(AHQgvVx(KhOcOVhjv<+J@+w J@$zv+0|3s}3u6EP From 4fd169a904ad6415123abd4f98ff4b7ee8631211 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Fri, 9 Sep 2016 06:01:18 +0000 Subject: [PATCH 36/61] fix lingregress (swapped z/y! ) --- nipype/algorithms/compcor.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index fac35d1e7d..e8e1a96271 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -69,7 +69,7 @@ def _run_interface(self, runtime): ''' timesteps = range(numvols) for voxel in range(numvoxels): - m, b, _, _, _ = stats.linregress(M[:, voxel], timesteps) + m, b, _, _, _ = stats.linregress(timesteps, M[:, voxel]) M[:, voxel] = M[:, voxel] - [m*t + b for t in timesteps] ''' # "... prior to column-wise variance normalization." From dfdd6b0fa669cc7d8a98af6052b3090b86ed676c Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Fri, 9 Sep 2016 06:36:45 +0000 Subject: [PATCH 37/61] more intelligent testing --- nipype/algorithms/compcor.py | 3 +-- nipype/algorithms/tests/test_compcor.py | 22 ++++++++++++++-------- 2 files changed, 15 insertions(+), 10 deletions(-) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index e8e1a96271..99dcd64ca2 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -66,12 +66,11 @@ def _run_interface(self, runtime): # "The constant and linear trends of the columns in the matrix M were removed ..." regress_poly(self.inputs.regress_poly_degree, voxel_timecourses) - ''' timesteps = range(numvols) for voxel in range(numvoxels): m, b, _, _, _ = stats.linregress(timesteps, M[:, voxel]) M[:, voxel] = M[:, voxel] - [m*t + b for t in timesteps] - ''' + # "... prior to column-wise variance normalization." M = M / self._compute_tSTD(M, 1.) diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index b112a5f697..b8a0a242a0 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -30,25 +30,31 @@ def test_compcor(self): mask[0,0,1] = 0 mask_file = utils.save_toy_nii(mask, self.filenames['masknii']) - expected_sha1 = 'b0dd7f9ab7ba8f516712eb0204dacc9e397fc6aa' + expected_components = [['-0.1989607212', '-0.5753813646'], + ['0.5692369697', '0.5674945949'], + ['-0.6662573243', '0.4675843432'], + ['0.4206466244', '-0.3361270124'], + ['-0.1246655485', '-0.1235705610']] ccresult = self.run_cc(CompCor(realigned_file=self.realigned_file, mask_file=mask_file), - expected_sha1) + expected_components) accresult = self.run_cc(ACompCor(realigned_file=self.realigned_file, mask_file=mask_file, components_file='acc_components_file'), - expected_sha1) + expected_components) assert_equal(os.path.getsize(ccresult.outputs.components_file), os.path.getsize(accresult.outputs.components_file)) def test_tcompcor(self): ccinterface = TCompCor(realigned_file=self.realigned_file) - self.run_cc(ccinterface, '12e54c07281a28ac0da3b934dce5c9d27626848a') + self.run_cc(ccinterface, [['-0.2846272268'], ['0.7115680670'], + ['-0.6048328569'], ['0.2134704201'], + ['-0.0355784033']]) - def run_cc(self, ccinterface, expected_components_data_sha1): + def run_cc(self, ccinterface, expected_components): # run ccresult = ccinterface.run() @@ -60,12 +66,12 @@ def run_cc(self, ccinterface, expected_components_data_sha1): assert_equal(ccinterface.inputs.num_components, 6) with open(ccresult.outputs.components_file, 'r') as components_file: - components_data = [line for line in 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]) - print(str(components_data), "comdata") - assert_equal(sha1(str(components_data)).hexdigest(), expected_components_data_sha1) + first_two = [row[:2] for row in components_data] + assert_equal(first_two, expected_components) return ccresult def tearDown(self): From 1995a0187fa52605cb12deedf70cb3c9adbfc094 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Fri, 9 Sep 2016 06:45:43 +0000 Subject: [PATCH 38/61] use regress_poly from TSNR in compcor --- nipype/algorithms/compcor.py | 22 +++++++++------------- 1 file changed, 9 insertions(+), 13 deletions(-) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index 99dcd64ca2..efec86071d 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -53,25 +53,21 @@ def _run_interface(self, runtime): voxel_timecourses[np.isnan(np.sum(voxel_timecourses, axis=1)), :] = 0 # from paper: + # "The constant and linear trends of the columns in the matrix M were + # removed [prior to ...]" + if self.inputs.use_regress_poly: + regressed = regress_poly(self.inputs.regress_poly_degree, + voxel_timecourses) + regressed = regressed - np.mean(regressed, axis=1)[:,None] + # "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 + M = regressed.T numvols = M.shape[0] numvoxels = M.shape[1] - if self.inputs.use_regress_poly: - # "The constant and linear trends of the columns in the matrix M were removed ..." - regress_poly(self.inputs.regress_poly_degree, voxel_timecourses) - - timesteps = range(numvols) - for voxel in range(numvoxels): - m, b, _, _, _ = stats.linregress(timesteps, M[:, voxel]) - M[:, voxel] = M[:, voxel] - [m*t + b for t in timesteps] - - # "... prior to column-wise variance normalization." + # "[... were removed] prior to column-wise variance normalization." M = M / self._compute_tSTD(M, 1.) # "The covariance matrix C = MMT was constructed and decomposed into its From 13838bda592e7ee7858c90f336713b5dc9353e8b Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Fri, 9 Sep 2016 20:03:19 +0000 Subject: [PATCH 39/61] add mask.nii back --- nipype/testing/data/mask.nii | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 nipype/testing/data/mask.nii diff --git a/nipype/testing/data/mask.nii b/nipype/testing/data/mask.nii new file mode 100644 index 0000000000..e69de29bb2 From 027d47eda08b18948be7eaede4dc3c6167ba15e1 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Fri, 9 Sep 2016 20:53:59 +0000 Subject: [PATCH 40/61] smarter tests --- nipype/algorithms/tests/test_compcor.py | 1 - nipype/algorithms/tests/test_tsnr.py | 53 +++++++++++++------------ 2 files changed, 27 insertions(+), 27 deletions(-) diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index b8a0a242a0..93f4ee79ec 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -9,7 +9,6 @@ import nibabel as nb import numpy as np import os -from hashlib import sha1 class TestCompCor(unittest.TestCase): ''' Note: Tests currently do a poor job of testing functionality ''' diff --git a/nipype/algorithms/tests/test_tsnr.py b/nipype/algorithms/tests/test_tsnr.py index c3a2f55e5d..b366e82a22 100644 --- a/nipype/algorithms/tests/test_tsnr.py +++ b/nipype/algorithms/tests/test_tsnr.py @@ -1,10 +1,10 @@ # 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 ...testing import assert_equal, assert_true, skipif, utils + +from ...testing import (assert_equal, assert_true, assert_almost_equal, + skipif, utils) from ..misc import TSNR -from hashlib import sha1 import unittest import mock import nibabel as nb @@ -39,9 +39,9 @@ def test_tsnr(self): # assert self.assert_expected_outputs(tsnrresult, { - 'mean_file': '1a55bcdf49901f25a2a838c90769989b9e4f2f19', - 'stddev_file': '0ba52a51fae90a9db6090c735432df5b742d663a', - 'tsnr_file': 'a794fc55766c8ad725437d3ff8b1153bd5f6e9b0' + 'mean_file': (2.8, 7.4), + 'stddev_file': (0.8, 2.9), + 'tsnr_file': (1.3, 9.25) }) def test_tsnr_withpoly1(self): @@ -51,10 +51,10 @@ def test_tsnr_withpoly1(self): # assert self.assert_expected_outputs_poly(tsnrresult, { - 'detrended_file': 'ee4f6c0b0e8c547617fc11aa50cf659436f9ccf0', - 'mean_file': '1a55bcdf49901f25a2a838c90769989b9e4f2f19', - 'stddev_file': 'e61d94e3cfea20b0c86c81bfdf80d82c55e9203b', - 'tsnr_file': 'a49f1cbd88f2aa71183dcd7aa4b86b17df622f0c' + 'detrended_file': (-0.1, 8.7), + 'mean_file': (2.8, 7.4), + 'stddev_file': (0.75, 2.75), + 'tsnr_file': (1.4, 9.9) }) def test_tsnr_withpoly2(self): @@ -64,10 +64,10 @@ def test_tsnr_withpoly2(self): # assert self.assert_expected_outputs_poly(tsnrresult, { - 'detrended_file': '22cb7f058d61cc090eb1a9dd7d31550bd4362a61', - 'mean_file': '4cee6776461f6bc238d11a55c0a8d1947a5732df', - 'stddev_file': '7267de4d9b63fcc553115c0198f1fa3bbb6a5a13', - 'tsnr_file': '1c6ed05f94806838f7b563df65900f37e60f8a6d' + 'detrended_file': (-0.22, 8.55), + 'mean_file': (2.8, 7.7), + 'stddev_file': (0.21, 2.4), + 'tsnr_file': (1.7, 35.9) }) def test_tsnr_withpoly3(self): @@ -77,20 +77,20 @@ def test_tsnr_withpoly3(self): # assert self.assert_expected_outputs_poly(tsnrresult, { - 'detrended_file': '3f2c1c7da233f128a7020b6fed79d6be2ec59fca', - 'mean_file': '4cee6776461f6bc238d11a55c0a8d1947a5732df', - 'stddev_file': '82bb793b76fab503d1d6b3e2d1b20a1bdebd7a2a', - 'tsnr_file': 'e004bd6096a0077b15058aabd4b0339bf6e21f64' + 'detrended_file': (1.8, 7.95), + 'mean_file': (2.8, 7.7), + 'stddev_file': (0.1, 1.7), + 'tsnr_file': (2.6, 57.3) }) - def assert_expected_outputs_poly(self, tsnrresult, hash_dict): + def assert_expected_outputs_poly(self, tsnrresult, expected_ranges): assert_equal(os.path.basename(tsnrresult.outputs.detrended_file), self.out_filenames['detrended_file']) - self.assert_expected_outputs(tsnrresult, hash_dict) + self.assert_expected_outputs(tsnrresult, expected_ranges) - def assert_expected_outputs(self, tsnrresult, hash_dict): + def assert_expected_outputs(self, tsnrresult, expected_ranges): self.assert_default_outputs(tsnrresult.outputs) - self.assert_unchanged(hash_dict) + self.assert_unchanged(expected_ranges) def assert_default_outputs(self, outputs): assert_equal(os.path.basename(outputs.mean_file), @@ -100,10 +100,11 @@ def assert_default_outputs(self, outputs): assert_equal(os.path.basename(outputs.tsnr_file), self.out_filenames['tsnr_file']) - def assert_unchanged(self, expected_hashes): - for key, hexhash in expected_hashes.iteritems(): - data = np.asanyarray(nb.load(self.out_filenames[key])._data) - assert_equal(sha1(str(data)).hexdigest(), hexhash) + def assert_unchanged(self, expected_ranges): + for key, (min_, max_) in expected_ranges.iteritems(): + data = np.asarray(nb.load(self.out_filenames[key])._data) + assert_almost_equal(np.amin(data), min_, decimal=1) + assert_almost_equal(np.amax(data), max_, decimal=1) def tearDown(self): utils.remove_nii(self.in_filenames.values()) From 54649023aae7d872b3b7f2842b03f66f45dc94b4 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Fri, 9 Sep 2016 21:33:25 +0000 Subject: [PATCH 41/61] rsfmri workflow test skeleton --- nipype/algorithms/tests/test_tsnr.py | 4 -- .../rsfmri/fsl/tests/test_resting.py | 53 +++++++++++++++++++ 2 files changed, 53 insertions(+), 4 deletions(-) create mode 100644 nipype/workflows/rsfmri/fsl/tests/test_resting.py diff --git a/nipype/algorithms/tests/test_tsnr.py b/nipype/algorithms/tests/test_tsnr.py index b366e82a22..391834e5db 100644 --- a/nipype/algorithms/tests/test_tsnr.py +++ b/nipype/algorithms/tests/test_tsnr.py @@ -13,10 +13,6 @@ class TestTSNR(unittest.TestCase): ''' Note: Tests currently do a poor job of testing functionality ''' - ''' - in_file = InputMultiPath(File(exists=True), mandatory=True, - regress_poly = traits.Range(low=1, desc='Remove polynomials') - ''' in_filenames = { 'in_file': 'tsnrinfile.nii', diff --git a/nipype/workflows/rsfmri/fsl/tests/test_resting.py b/nipype/workflows/rsfmri/fsl/tests/test_resting.py new file mode 100644 index 0000000000..2578461237 --- /dev/null +++ b/nipype/workflows/rsfmri/fsl/tests/test_resting.py @@ -0,0 +1,53 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: + +from nipype.testing import (assert_equal, assert_true, assert_almost_equal, + skipif, utils) +from nipype.workflows.rsfmri.fsl import resting + +import unittest +import mock +import nibabel as nb +import numpy as np +import os + +class TestResting(unittest.TestCase): + + in_filenames = { + 'in_file': 'rsfmrifunc.nii', + } + + out_filenames = { + 'noise_mask_file': '', + 'filtered_file': '' + } + + def setUp(self): + # setup + utils.save_toy_nii(self.fake_data, self.in_filenames['in_file']) + + @skipif(True) + def test_create_resting_preproc(self): + # setup + print(np.random.randint(0, 10, (2, 2, 2, 5))) + # run + wf = resting.create_resting_preproc() + wf.run() + + # assert + + def tearDown(self): + utils.remove_nii(self.in_filenames.values()) + utils.remove_nii(self.out_filenames.values()) + + fake_data = np.array([[[[2, 4, 3, 9, 1], + [3, 6, 4, 7, 4]], + + [[8, 3, 4, 6, 2], + [4, 0, 4, 4, 2]]], + + [[[9, 7, 5, 5, 7], + [7, 8, 4, 8, 4]], + + [[0, 4, 7, 1, 7], + [6, 8, 8, 8, 7]]]]) From f8cea6145666b64ef240924abdec362b014694e2 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Sat, 10 Sep 2016 03:31:46 +0000 Subject: [PATCH 42/61] set up mocks --- .../rsfmri/fsl/tests/test_resting.py | 25 +++++++++++++++---- 1 file changed, 20 insertions(+), 5 deletions(-) diff --git a/nipype/workflows/rsfmri/fsl/tests/test_resting.py b/nipype/workflows/rsfmri/fsl/tests/test_resting.py index 2578461237..c1ad16fa55 100644 --- a/nipype/workflows/rsfmri/fsl/tests/test_resting.py +++ b/nipype/workflows/rsfmri/fsl/tests/test_resting.py @@ -1,16 +1,27 @@ # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- # vi: set ft=python sts=4 ts=4 sw=4 et: -from nipype.testing import (assert_equal, assert_true, assert_almost_equal, - skipif, utils) -from nipype.workflows.rsfmri.fsl import resting +from .....testing import (assert_equal, assert_true, assert_almost_equal, + skipif, utils) +from .....interfaces import fsl + +from ..resting import create_resting_preproc import unittest import mock +from mock import MagicMock import nibabel as nb import numpy as np import os +def mock_node_factory(*args, **kwargs): + mock = MagicMock() + mock.name = kwargs['name'] if 'name' in kwargs.keys() else '' + if 'interface' in kwargs.keys(): + mock = mock.create_autospec(kwargs['interface'], instance=True) + mock.iterables = None + return mock + class TestResting(unittest.TestCase): in_filenames = { @@ -27,11 +38,15 @@ def setUp(self): utils.save_toy_nii(self.fake_data, self.in_filenames['in_file']) @skipif(True) - def test_create_resting_preproc(self): + @mock.patch('nipype.pipeline.engine.Workflow._write_report_info') + @mock.patch('nipype.workflows.rsfmri.fsl.resting.create_realign_flow', + side_effect=mock_node_factory) + @mock.patch('nipype.pipeline.engine.Node', side_effect=mock_node_factory) + def test_create_resting_preproc(self, mock_Node, mock_realign_wf, nothing): # setup print(np.random.randint(0, 10, (2, 2, 2, 5))) # run - wf = resting.create_resting_preproc() + wf = create_resting_preproc() wf.run() # assert From 1c655f84dd1538a1a7f825c105b99e4a983f83b2 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Sat, 10 Sep 2016 20:36:47 +0000 Subject: [PATCH 43/61] mock only what isn't being tested --- nipype/workflows/rsfmri/fsl/resting.py | 2 +- .../rsfmri/fsl/tests/test_resting.py | 25 ++++++++++++++++--- 2 files changed, 22 insertions(+), 5 deletions(-) diff --git a/nipype/workflows/rsfmri/fsl/resting.py b/nipype/workflows/rsfmri/fsl/resting.py index 2bf8f28bd8..6906504b74 100644 --- a/nipype/workflows/rsfmri/fsl/resting.py +++ b/nipype/workflows/rsfmri/fsl/resting.py @@ -151,7 +151,7 @@ def create_resting_preproc(name='restpreproc'): 'num_components'], output_names=['noise_components'], function=extract_noise_components), - name='compcorr') + name='compcor') remove_noise = pe.Node(fsl.FilterRegressor(filter_all=True), name='remove_noise') bandpass_filter = pe.Node(fsl.TemporalFilter(), diff --git a/nipype/workflows/rsfmri/fsl/tests/test_resting.py b/nipype/workflows/rsfmri/fsl/tests/test_resting.py index c1ad16fa55..5ff4902360 100644 --- a/nipype/workflows/rsfmri/fsl/tests/test_resting.py +++ b/nipype/workflows/rsfmri/fsl/tests/test_resting.py @@ -3,7 +3,8 @@ from .....testing import (assert_equal, assert_true, assert_almost_equal, skipif, utils) -from .....interfaces import fsl +from .....interfaces import fsl, IdentityInterface +from .....pipeline.engine import Node from ..resting import create_resting_preproc @@ -15,8 +16,24 @@ import os def mock_node_factory(*args, **kwargs): + ''' return mocks for all the nodes except compcor and compcor's neighbors''' mock = MagicMock() - mock.name = kwargs['name'] if 'name' in kwargs.keys() else '' + if 'name' in kwargs.keys(): + name = kwargs['name'] + if name == 'compcor': + return Node(*args, **kwargs) + if name in ['threshold', 'inputspec']: + # nodes that provide inputs for compcor + # just give all of them all of the fields + return Node(IdentityInterface(fields=['out_file', 'lowpass_sigma', + 'num_noise_components', + 'func', 'highpass_sigma']), + name='fake'+name) + if name in ('remove_noise'): + # node that takes output from compcor + return Node(IdentityInterface(fields=['design_file', 'out_file']), + name=name) + mock.name = kwargs['name'] if 'interface' in kwargs.keys(): mock = mock.create_autospec(kwargs['interface'], instance=True) mock.iterables = None @@ -40,11 +57,11 @@ def setUp(self): @skipif(True) @mock.patch('nipype.pipeline.engine.Workflow._write_report_info') @mock.patch('nipype.workflows.rsfmri.fsl.resting.create_realign_flow', - side_effect=mock_node_factory) + return_value=Node(name='realigner', interface=IdentityInterface( + fields=['outputspec.realigned_file']))) @mock.patch('nipype.pipeline.engine.Node', side_effect=mock_node_factory) def test_create_resting_preproc(self, mock_Node, mock_realign_wf, nothing): # setup - print(np.random.randint(0, 10, (2, 2, 2, 5))) # run wf = create_resting_preproc() wf.run() From 08bccfacd7885905b0873abb998f486577eaee63 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Sat, 10 Sep 2016 21:13:37 +0000 Subject: [PATCH 44/61] run compcor tests in tempfile --- nipype/algorithms/tests/test_compcor.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index 93f4ee79ec..e06408d554 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -9,6 +9,8 @@ import nibabel as nb import numpy as np import os +import tempfile +import shutil class TestCompCor(unittest.TestCase): ''' Note: Tests currently do a poor job of testing functionality ''' @@ -19,6 +21,9 @@ class TestCompCor(unittest.TestCase): def setUp(self): # setup + self.orig_dir = os.getcwd() + self.temp_dir = tempfile.mkdtemp() + os.chdir(self.temp_dir) noise = np.fromfunction(self.fake_noise_fun, self.fake_data.shape) self.realigned_file = utils.save_toy_nii(self.fake_data + noise, self.filenames['functionalnii']) @@ -74,7 +79,8 @@ def run_cc(self, ccinterface, expected_components): return ccresult def tearDown(self): - utils.remove_nii(self.filenames.values()) + os.chdir(self.orig_dir) + shutil.rmtree(self.temp_dir) def fake_noise_fun(self, i, j, l, m): return m*i + l - j From 4401bb84b609d443849f0c18f5415769fe98cff3 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Sat, 10 Sep 2016 22:48:32 +0000 Subject: [PATCH 45/61] set up 2 out of 3 inputs for compcor node --- nipype/workflows/rsfmri/fsl/tests/test_resting.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/nipype/workflows/rsfmri/fsl/tests/test_resting.py b/nipype/workflows/rsfmri/fsl/tests/test_resting.py index 5ff4902360..1fdba0590c 100644 --- a/nipype/workflows/rsfmri/fsl/tests/test_resting.py +++ b/nipype/workflows/rsfmri/fsl/tests/test_resting.py @@ -15,9 +15,15 @@ import numpy as np import os + + + + def mock_node_factory(*args, **kwargs): ''' return mocks for all the nodes except compcor and compcor's neighbors''' mock = MagicMock() + if 'interface' in kwargs.keys(): + mock = mock.create_autospec(kwargs['interface'], instance=True) if 'name' in kwargs.keys(): name = kwargs['name'] if name == 'compcor': @@ -28,14 +34,12 @@ def mock_node_factory(*args, **kwargs): return Node(IdentityInterface(fields=['out_file', 'lowpass_sigma', 'num_noise_components', 'func', 'highpass_sigma']), - name='fake'+name) + name=name) if name in ('remove_noise'): # node that takes output from compcor return Node(IdentityInterface(fields=['design_file', 'out_file']), name=name) mock.name = kwargs['name'] - if 'interface' in kwargs.keys(): - mock = mock.create_autospec(kwargs['interface'], instance=True) mock.iterables = None return mock @@ -43,6 +47,7 @@ class TestResting(unittest.TestCase): in_filenames = { 'in_file': 'rsfmrifunc.nii', + 'noise_mask_file': 'rsfmrimask.nii' } out_filenames = { @@ -64,6 +69,8 @@ def test_create_resting_preproc(self, mock_Node, mock_realign_wf, nothing): # setup # run wf = create_resting_preproc() + wf.inputs.inputspec.num_noise_components = 6 + wf.get_node('threshold').inputs.out_file = self.in_filenames['noise_mask_file'] wf.run() # assert From 402b6b1ccaea130cfb838387b6d1efa5a9ca486b Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Sun, 11 Sep 2016 04:34:42 +0000 Subject: [PATCH 46/61] add __init__.py file --- nipype/workflows/rsfmri/fsl/tests/__init__.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 nipype/workflows/rsfmri/fsl/tests/__init__.py diff --git a/nipype/workflows/rsfmri/fsl/tests/__init__.py b/nipype/workflows/rsfmri/fsl/tests/__init__.py new file mode 100644 index 0000000000..e69de29bb2 From c18a780fa8fec94151c4bb3040415798589bef1a Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Sun, 11 Sep 2016 04:48:34 +0000 Subject: [PATCH 47/61] mocked-up rsfmri workflow runs; use tempfiles --- .../rsfmri/fsl/tests/test_resting.py | 89 ++++++++++--------- 1 file changed, 45 insertions(+), 44 deletions(-) diff --git a/nipype/workflows/rsfmri/fsl/tests/test_resting.py b/nipype/workflows/rsfmri/fsl/tests/test_resting.py index 1fdba0590c..62b714a0d0 100644 --- a/nipype/workflows/rsfmri/fsl/tests/test_resting.py +++ b/nipype/workflows/rsfmri/fsl/tests/test_resting.py @@ -4,7 +4,7 @@ from .....testing import (assert_equal, assert_true, assert_almost_equal, skipif, utils) from .....interfaces import fsl, IdentityInterface -from .....pipeline.engine import Node +from .....pipeline.engine import Node, Workflow from ..resting import create_resting_preproc @@ -14,41 +14,35 @@ import nibabel as nb import numpy as np import os - - - - - -def mock_node_factory(*args, **kwargs): - ''' return mocks for all the nodes except compcor and compcor's neighbors''' - mock = MagicMock() - if 'interface' in kwargs.keys(): - mock = mock.create_autospec(kwargs['interface'], instance=True) - if 'name' in kwargs.keys(): - name = kwargs['name'] - if name == 'compcor': - return Node(*args, **kwargs) - if name in ['threshold', 'inputspec']: - # nodes that provide inputs for compcor - # just give all of them all of the fields - return Node(IdentityInterface(fields=['out_file', 'lowpass_sigma', - 'num_noise_components', - 'func', 'highpass_sigma']), - name=name) - if name in ('remove_noise'): - # node that takes output from compcor - return Node(IdentityInterface(fields=['design_file', 'out_file']), - name=name) - mock.name = kwargs['name'] - mock.iterables = None - return mock +import tempfile +import shutil + +all_fields = ['func', 'in_file', 'slice_time_corrected_file', 'stddev_file', + 'out_stat', 'thresh', 'num_noise_components', 'detrended_file', + 'design_file', 'highpass_sigma', 'lowpass_sigma', 'out_file', + 'noise_mask_file', 'filtered_file'] + +def stub_node_factory(*args, **kwargs): + if 'name' not in kwargs.keys(): + raise Exception() + name = kwargs['name'] + if name == 'compcor': + return Node(*args, **kwargs) + else: # replace with an IdentityInterface + return Node(IdentityInterface(fields=all_fields), + name=name) + +def stub_wf(*args, **kwargs): + wf = Workflow(name='realigner') + inputnode = Node(IdentityInterface(fields=['func']), name='inputspec') + outputnode = Node(interface=IdentityInterface(fields=['realigned_file']), + name='outputspec') + wf.connect(inputnode, 'func', outputnode, 'realigned_file') + return wf class TestResting(unittest.TestCase): - in_filenames = { - 'in_file': 'rsfmrifunc.nii', - 'noise_mask_file': 'rsfmrimask.nii' - } + in_filenames = {} out_filenames = { 'noise_mask_file': '', @@ -57,27 +51,34 @@ class TestResting(unittest.TestCase): def setUp(self): # setup - utils.save_toy_nii(self.fake_data, self.in_filenames['in_file']) + self.orig_dir = os.getcwd() + self.temp_dir = tempfile.mkdtemp() + os.chdir(self.temp_dir) + self.in_filenames['realigned_file'] = utils.save_toy_nii(self.fake_data, os.path.abspath('rsfmrifunc.nii')) + mask = np.zeros(self.fake_data.shape[:3]) + for i in range(mask.shape[0]): + for j in range(mask.shape[1]): + if i==j: + mask[i,j] = 1 + self.in_filenames['noise_mask_file'] = utils.save_toy_nii(mask, os.path.abspath('rsfmrimask.nii')) - @skipif(True) - @mock.patch('nipype.pipeline.engine.Workflow._write_report_info') @mock.patch('nipype.workflows.rsfmri.fsl.resting.create_realign_flow', - return_value=Node(name='realigner', interface=IdentityInterface( - fields=['outputspec.realigned_file']))) - @mock.patch('nipype.pipeline.engine.Node', side_effect=mock_node_factory) - def test_create_resting_preproc(self, mock_Node, mock_realign_wf, nothing): - # setup - # run + side_effect=stub_wf) + @mock.patch('nipype.pipeline.engine.Node', side_effect=stub_node_factory) + def test_create_resting_preproc(self, mock_Node, mock_realign_wf): wf = create_resting_preproc() + wf.inputs.inputspec.num_noise_components = 6 + wf.get_node('slicetimer').inputs.slice_time_corrected_file = self.in_filenames['realigned_file'] wf.get_node('threshold').inputs.out_file = self.in_filenames['noise_mask_file'] + wf.run() # assert def tearDown(self): - utils.remove_nii(self.in_filenames.values()) - utils.remove_nii(self.out_filenames.values()) + os.chdir(self.orig_dir) + shutil.rmtree(self.temp_dir) fake_data = np.array([[[[2, 4, 3, 9, 1], [3, 6, 4, 7, 4]], From 083c16f53f105752e7e48550487054e6efe23bd7 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Sun, 11 Sep 2016 06:58:26 +0000 Subject: [PATCH 48/61] completed rsfmri wf test --- nipype/workflows/rsfmri/fsl/resting.py | 4 +- .../rsfmri/fsl/tests/test_resting.py | 44 ++++++++++++++----- 2 files changed, 35 insertions(+), 13 deletions(-) diff --git a/nipype/workflows/rsfmri/fsl/resting.py b/nipype/workflows/rsfmri/fsl/resting.py index 8d0befff7f..a21022f41e 100644 --- a/nipype/workflows/rsfmri/fsl/resting.py +++ b/nipype/workflows/rsfmri/fsl/resting.py @@ -96,7 +96,7 @@ def create_realign_flow(name='realign'): return realignflow -def create_resting_preproc(name='restpreproc'): +def create_resting_preproc(name='restpreproc', base_dir=None): """Create a "resting" time series preprocessing workflow The noise removal is based on Behzadi et al. (2007) @@ -128,7 +128,7 @@ def create_resting_preproc(name='restpreproc'): """ - restpreproc = pe.Workflow(name=name) + restpreproc = pe.Workflow(name=name, base_dir=base_dir) # Define nodes inputnode = pe.Node(interface=util.IdentityInterface(fields=['func', diff --git a/nipype/workflows/rsfmri/fsl/tests/test_resting.py b/nipype/workflows/rsfmri/fsl/tests/test_resting.py index 62b714a0d0..6590283748 100644 --- a/nipype/workflows/rsfmri/fsl/tests/test_resting.py +++ b/nipype/workflows/rsfmri/fsl/tests/test_resting.py @@ -3,7 +3,7 @@ from .....testing import (assert_equal, assert_true, assert_almost_equal, skipif, utils) -from .....interfaces import fsl, IdentityInterface +from .....interfaces import fsl, IdentityInterface, utility from .....pipeline.engine import Node, Workflow from ..resting import create_resting_preproc @@ -42,39 +42,61 @@ def stub_wf(*args, **kwargs): class TestResting(unittest.TestCase): - in_filenames = {} + in_filenames = { + 'realigned_file': 'rsfmrifunc.nii', + 'mask_file': 'rsfmrimask.nii' + } out_filenames = { - 'noise_mask_file': '', - 'filtered_file': '' + 'components_file': 'restpreproc/compcor/noise_components.txt' } + num_noise_components = 6 + def setUp(self): - # setup + # setup temp folder self.orig_dir = os.getcwd() self.temp_dir = tempfile.mkdtemp() os.chdir(self.temp_dir) - self.in_filenames['realigned_file'] = utils.save_toy_nii(self.fake_data, os.path.abspath('rsfmrifunc.nii')) + self.in_filenames = {key: os.path.abspath(value) + for key, value in self.in_filenames.items()} + + # create&save input files + utils.save_toy_nii(self.fake_data, self.in_filenames['realigned_file']) mask = np.zeros(self.fake_data.shape[:3]) for i in range(mask.shape[0]): for j in range(mask.shape[1]): if i==j: mask[i,j] = 1 - self.in_filenames['noise_mask_file'] = utils.save_toy_nii(mask, os.path.abspath('rsfmrimask.nii')) + utils.save_toy_nii(mask, self.in_filenames['mask_file']) @mock.patch('nipype.workflows.rsfmri.fsl.resting.create_realign_flow', side_effect=stub_wf) @mock.patch('nipype.pipeline.engine.Node', side_effect=stub_node_factory) def test_create_resting_preproc(self, mock_Node, mock_realign_wf): - wf = create_resting_preproc() + wf = create_resting_preproc(base_dir=os.getcwd()) - wf.inputs.inputspec.num_noise_components = 6 - wf.get_node('slicetimer').inputs.slice_time_corrected_file = self.in_filenames['realigned_file'] - wf.get_node('threshold').inputs.out_file = self.in_filenames['noise_mask_file'] + wf.inputs.inputspec.num_noise_components = self.num_noise_components + mask_in = wf.get_node('threshold').inputs + mask_in.out_file = self.in_filenames['mask_file'] + func_in = wf.get_node('slicetimer').inputs + func_in.slice_time_corrected_file = self.in_filenames['realigned_file'] wf.run() # assert + expected_file = os.path.abspath(self.out_filenames['components_file']) + with open(expected_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 == self.num_noise_components + or num_got_components == self.fake_data.shape[3]) + first_two = [row[:2] for row in components_data] + assert_equal(first_two, [['-0.5172356654', '-0.6973053243'], + ['0.2574722644', '0.1645270737'], + ['-0.0806469590', '0.5156853779'], + ['0.7187176051', '-0.3235820287'], + ['-0.3783072450', '0.3406749013']]) def tearDown(self): os.chdir(self.orig_dir) From be215f5f46622bf4ab12a54236bbf04a953404cc Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Sun, 11 Sep 2016 07:20:53 +0000 Subject: [PATCH 49/61] refactor rsfmri fsl workflow --- nipype/algorithms/compcor.py | 7 +++-- nipype/workflows/rsfmri/fsl/resting.py | 43 +++----------------------- 2 files changed, 9 insertions(+), 41 deletions(-) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index efec86071d..af8e65948f 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -56,14 +56,15 @@ def _run_interface(self, runtime): # "The constant and linear trends of the columns in the matrix M were # removed [prior to ...]" if self.inputs.use_regress_poly: - regressed = regress_poly(self.inputs.regress_poly_degree, + voxel_timecourses = regress_poly(self.inputs.regress_poly_degree, voxel_timecourses) - regressed = regressed - np.mean(regressed, axis=1)[:,None] + voxel_timecourses = voxel_timecourses - np.mean(voxel_timecourses, + axis=1)[:, np.newaxis] # "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." - M = regressed.T + M = voxel_timecourses.T numvols = M.shape[0] numvoxels = M.shape[1] diff --git a/nipype/workflows/rsfmri/fsl/resting.py b/nipype/workflows/rsfmri/fsl/resting.py index a21022f41e..a2778fe2d9 100644 --- a/nipype/workflows/rsfmri/fsl/resting.py +++ b/nipype/workflows/rsfmri/fsl/resting.py @@ -8,37 +8,7 @@ from ....interfaces import utility as util # utility from ....pipeline import engine as pe # pypeline engine from ....algorithms.misc import TSNR - - -def extract_noise_components(realigned_file, noise_mask_file, num_components): - """Derive components most reflective of physiological noise - """ - import os - from nibabel import load - import numpy as np - import scipy as sp - imgseries = load(realigned_file) - components = None - mask = load(noise_mask_file).get_data() - voxel_timecourses = imgseries.get_data()[mask > 0] - voxel_timecourses[np.isnan(np.sum(voxel_timecourses, axis=1)), :] = 0 - # remove mean and normalize by variance - # voxel_timecourses.shape == [nvoxels, time] - X = voxel_timecourses.T - stdX = np.std(X, axis=0) - stdX[stdX == 0] = 1. - stdX[np.isnan(stdX)] = 1. - stdX[np.isinf(stdX)] = 1. - X = (X - np.mean(X, axis=0)) / stdX - u, _, _ = sp.linalg.svd(X, full_matrices=False) - if components is None: - components = u[:, :num_components] - else: - components = np.hstack((components, u[:, :num_components])) - components_file = os.path.join(os.getcwd(), 'noise_components.txt') - np.savetxt(components_file, components, fmt=b"%.10f") - return components_file - +from ....algorithms import compcor as cc def select_volume(filename, which): """Return the middle index of a file @@ -148,11 +118,8 @@ def create_resting_preproc(name='restpreproc', base_dir=None): getthresh = pe.Node(interface=fsl.ImageStats(op_string='-p 98'), name='getthreshold') threshold_stddev = pe.Node(fsl.Threshold(), name='threshold') - compcor = pe.Node(util.Function(input_names=['realigned_file', - 'noise_mask_file', - 'num_components'], - output_names=['noise_components'], - function=extract_noise_components), + compcor = pe.Node(cc.ACompCor(components_file="noise_components.txt", + use_regress_poly=False), name='compcor') remove_noise = pe.Node(fsl.FilterRegressor(filter_all=True), name='remove_noise') @@ -170,12 +137,12 @@ def create_resting_preproc(name='restpreproc', base_dir=None): restpreproc.connect(realigner, 'outputspec.realigned_file', compcor, 'realigned_file') restpreproc.connect(threshold_stddev, 'out_file', - compcor, 'noise_mask_file') + compcor, 'mask_file') restpreproc.connect(inputnode, 'num_noise_components', compcor, 'num_components') restpreproc.connect(tsnr, 'detrended_file', remove_noise, 'in_file') - restpreproc.connect(compcor, 'noise_components', + restpreproc.connect(compcor, 'components_file', remove_noise, 'design_file') restpreproc.connect(inputnode, 'highpass_sigma', bandpass_filter, 'highpass_sigma') From 5f39bfcfb027410d5cd214c72382e7755c48c084 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Sun, 11 Sep 2016 07:52:45 +0000 Subject: [PATCH 50/61] make python 3 happy --- nipype/algorithms/tests/test_tsnr.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nipype/algorithms/tests/test_tsnr.py b/nipype/algorithms/tests/test_tsnr.py index 391834e5db..fcf3b19b9b 100644 --- a/nipype/algorithms/tests/test_tsnr.py +++ b/nipype/algorithms/tests/test_tsnr.py @@ -97,7 +97,7 @@ def assert_default_outputs(self, outputs): self.out_filenames['tsnr_file']) def assert_unchanged(self, expected_ranges): - for key, (min_, max_) in expected_ranges.iteritems(): + for key, (min_, max_) in expected_ranges.items(): data = np.asarray(nb.load(self.out_filenames[key])._data) assert_almost_equal(np.amin(data), min_, decimal=1) assert_almost_equal(np.amax(data), max_, decimal=1) From 6869a5f9125f43790c7f5b885c6c6b9d3655d522 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Sun, 11 Sep 2016 07:57:31 +0000 Subject: [PATCH 51/61] use mkdtemp for test_tsnr --- nipype/algorithms/tests/test_tsnr.py | 12 +++++++++--- nipype/testing/utils.py | 8 -------- 2 files changed, 9 insertions(+), 11 deletions(-) diff --git a/nipype/algorithms/tests/test_tsnr.py b/nipype/algorithms/tests/test_tsnr.py index fcf3b19b9b..3833d1f27f 100644 --- a/nipype/algorithms/tests/test_tsnr.py +++ b/nipype/algorithms/tests/test_tsnr.py @@ -10,6 +10,8 @@ import nibabel as nb import numpy as np import os +import tempfile +import shutil class TestTSNR(unittest.TestCase): ''' Note: Tests currently do a poor job of testing functionality ''' @@ -26,7 +28,11 @@ class TestTSNR(unittest.TestCase): } def setUp(self): - # setup + # setup temp folder + self.orig_dir = os.getcwd() + self.temp_dir = tempfile.mkdtemp() + os.chdir(self.temp_dir) + utils.save_toy_nii(self.fake_data, self.in_filenames['in_file']) def test_tsnr(self): @@ -103,8 +109,8 @@ def assert_unchanged(self, expected_ranges): assert_almost_equal(np.amax(data), max_, decimal=1) def tearDown(self): - utils.remove_nii(self.in_filenames.values()) - utils.remove_nii(self.out_filenames.values()) + os.chdir(self.orig_dir) + shutil.rmtree(self.temp_dir) fake_data = np.array([[[[2, 4, 3, 9, 1], [3, 6, 4, 7, 4]], diff --git a/nipype/testing/utils.py b/nipype/testing/utils.py index be924b36fb..092f53ea8e 100644 --- a/nipype/testing/utils.py +++ b/nipype/testing/utils.py @@ -110,11 +110,3 @@ def save_toy_nii(ndarray, filename): toy = nb.Nifti1Image(ndarray, np.eye(4)) nb.nifti1.save(toy, filename) return filename - -def remove_nii(filenames): - ''' remove temporary nifti files''' - for filename in filenames: - try: - os.remove(filename) - except (OSError, TypeError) as e: - pass From c9aa6d6d06431d9b23a9966b10716c11063f702d Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Mon, 12 Sep 2016 01:32:20 +0000 Subject: [PATCH 52/61] enable extra_regressors (needed by rsfmri_surface_vol...) --- nipype/algorithms/compcor.py | 10 ++++++++++ nipype/algorithms/tests/test_compcor.py | 17 ++++++++++++----- 2 files changed, 22 insertions(+), 5 deletions(-) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index af8e65948f..291d2af2e1 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -23,6 +23,8 @@ class CompCorInputSpec(BaseInterfaceInputSpec): 'pre-component extraction') regress_poly_degree = traits.Range(low=1, default=1, usedefault=True, desc='the degree polynomial to use') + extra_regressors = File(exists=True, mandatory=False, + desc='additional regressors to add') class CompCorOutputSpec(TraitedSpec): components_file = File(exists=True, @@ -75,6 +77,10 @@ def _run_interface(self, runtime): # principal components using a singular value decomposition." u, _, _ = linalg.svd(M, full_matrices=False) components = u[:, :self.inputs.num_components] + if self.inputs.extra_regressors: + components = self._add_extras(components, + self.inputs.extra_regressors) + components_file = os.path.join(os.getcwd(), self.inputs.components_file) np.savetxt(components_file, components, fmt="%.10f") return runtime @@ -92,6 +98,10 @@ def _compute_tSTD(self, M, x): stdM[np.isinf(stdM)] = x return stdM + def _add_extras(self, components, extra_regressors): + regressors = np.genfromtxt(self.inputs.extra_regressors) + return np.hstack((components, regressors)) + class TCompCor(CompCor): def _run_interface(self, runtime): diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index e06408d554..8e30ce2e8b 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -27,13 +27,12 @@ def setUp(self): noise = np.fromfunction(self.fake_noise_fun, self.fake_data.shape) self.realigned_file = utils.save_toy_nii(self.fake_data + noise, self.filenames['functionalnii']) - - def test_compcor(self): mask = np.ones(self.fake_data.shape[:3]) mask[0,0,0] = 0 mask[0,0,1] = 0 - mask_file = utils.save_toy_nii(mask, self.filenames['masknii']) + self.mask_file = utils.save_toy_nii(mask, self.filenames['masknii']) + def test_compcor(self): expected_components = [['-0.1989607212', '-0.5753813646'], ['0.5692369697', '0.5674945949'], ['-0.6662573243', '0.4675843432'], @@ -41,17 +40,25 @@ def test_compcor(self): ['-0.1246655485', '-0.1235705610']] ccresult = self.run_cc(CompCor(realigned_file=self.realigned_file, - mask_file=mask_file), + mask_file=self.mask_file), expected_components) accresult = self.run_cc(ACompCor(realigned_file=self.realigned_file, - mask_file=mask_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)) + @mock.patch('nipype.algorithms.compcor.CompCor._add_extras') + def test_compcor_with_extra_regressors(self, mock_add_extras): + regressors_file ='regress.txt' + open(regressors_file, 'a').close() # make sure file exists + CompCor(realigned_file=self.realigned_file, mask_file=self.mask_file, + extra_regressors=regressors_file).run() + assert_true(mock_add_extras.called) + def test_tcompcor(self): ccinterface = TCompCor(realigned_file=self.realigned_file) self.run_cc(ccinterface, [['-0.2846272268'], ['0.7115680670'], From 9aa2300457f9c65b4cff206e86dc63805342b62a Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Mon, 12 Sep 2016 02:43:11 +0000 Subject: [PATCH 53/61] make percentile threshold a parameter for tCompCor --- nipype/algorithms/compcor.py | 22 +++++++++++++++++----- nipype/algorithms/tests/test_compcor.py | 8 ++++++++ 2 files changed, 25 insertions(+), 5 deletions(-) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index 291d2af2e1..601093f4fa 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -102,8 +102,20 @@ def _add_extras(self, components, extra_regressors): regressors = np.genfromtxt(self.inputs.extra_regressors) return np.hstack((components, regressors)) +class TCompCorInputSpec(CompCorInputSpec): + # and all the fields in CompCorInputSpec + percentile_threshold = traits.Range(low=0., high=1., value=.02, + exclude_low=True, exclude_high=True, + usedefault=True, desc='the percentile ' + 'used to select highest-variance ' + 'voxels. By default the 2% of voxels ' + 'with the highest variance are used.') + class TCompCor(CompCor): + input_spec = TCompCorInputSpec + output_spec = CompCorOutputSpec + def _run_interface(self, runtime): imgseries = nb.load(self.inputs.realigned_file).get_data() time_voxels = imgseries.T @@ -119,16 +131,16 @@ def _run_interface(self, runtime): tSTD = self._compute_tSTD(time_voxels, 0) sortSTD = np.sort(tSTD, axis=None) # flattened sorted matrix - # "... and retained a pre-specified upper fraction of the sorted voxels - # within each slice ... we chose a 2% threshold" - threshold = sortSTD[int(num_voxels * .98)] - mask = tSTD >= threshold + # use percentile_threshold to pick voxels + threshold_index = int(num_voxels * (1. - self.inputs.percentile_threshold)) + threshold_std = sortSTD[threshold_index] + mask = tSTD >= threshold_std mask = mask.astype(int) # save mask mask_file = 'mask.nii' nb.nifti1.save(nb.Nifti1Image(mask, np.eye(4)), mask_file) - self.inputs.mask_file = 'mask.nii' + self.inputs.mask_file = mask_file super(TCompCor, self)._run_interface(runtime) return runtime diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index 8e30ce2e8b..d8a6e74c02 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -65,6 +65,14 @@ def test_tcompcor(self): ['-0.6048328569'], ['0.2134704201'], ['-0.0355784033']]) + def test_tcompcor_with_percentile(self): + ccinterface = TCompCor(realigned_file=self.realigned_file, percentile_threshold=0.2) + ccinterface.run() + + mask = nb.load('mask.nii').get_data() + num_nonmasked_voxels = np.count_nonzero(mask) + assert_equal(num_nonmasked_voxels, 2) + def run_cc(self, ccinterface, expected_components): # run ccresult = ccinterface.run() From 276fb576e04b8a3371b24c7fc69c7ae4c612c62b Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Mon, 12 Sep 2016 02:58:22 +0000 Subject: [PATCH 54/61] use CompCor interface in example --- .../rsfmri_vol_surface_preprocessing_nipy.py | 60 ++----------------- 1 file changed, 6 insertions(+), 54 deletions(-) diff --git a/examples/rsfmri_vol_surface_preprocessing_nipy.py b/examples/rsfmri_vol_surface_preprocessing_nipy.py index b46ae28ab6..4f5bae4c8d 100644 --- a/examples/rsfmri_vol_surface_preprocessing_nipy.py +++ b/examples/rsfmri_vol_surface_preprocessing_nipy.py @@ -64,6 +64,7 @@ from nipype.algorithms.rapidart import ArtifactDetect from nipype.algorithms.misc import TSNR +from nipype.algorithms.compcor import ACompCor from nipype.interfaces.utility import Rename, Merge, IdentityInterface from nipype.utils.filemanip import filename_to_list from nipype.interfaces.io import DataSink, FreeSurferSource @@ -228,51 +229,6 @@ def build_filter1(motion_params, comp_norm, outliers, detrend_poly=None): out_files.append(filename) return out_files - -def extract_noise_components(realigned_file, mask_file, num_components=5, - extra_regressors=None): - """Derive components most reflective of physiological noise - - Parameters - ---------- - realigned_file: a 4D Nifti file containing realigned volumes - mask_file: a 3D Nifti file containing white matter + ventricular masks - num_components: number of components to use for noise decomposition - extra_regressors: additional regressors to add - - Returns - ------- - components_file: a text file containing the noise components - """ - imgseries = nb.load(realigned_file) - components = None - for filename in filename_to_list(mask_file): - mask = nb.load(filename).get_data() - if len(np.nonzero(mask > 0)[0]) == 0: - continue - voxel_timecourses = imgseries.get_data()[mask > 0] - voxel_timecourses[np.isnan(np.sum(voxel_timecourses, axis=1)), :] = 0 - # remove mean and normalize by variance - # voxel_timecourses.shape == [nvoxels, time] - X = voxel_timecourses.T - stdX = np.std(X, axis=0) - stdX[stdX == 0] = 1. - stdX[np.isnan(stdX)] = 1. - stdX[np.isinf(stdX)] = 1. - X = (X - np.mean(X, axis=0)) / stdX - u, _, _ = sp.linalg.svd(X, full_matrices=False) - if components is None: - components = u[:, :num_components] - else: - components = np.hstack((components, u[:, :num_components])) - if extra_regressors: - regressors = np.genfromtxt(extra_regressors) - components = np.hstack((components, regressors)) - components_file = os.path.join(os.getcwd(), 'noise_components.txt') - np.savetxt(components_file, components, fmt=b"%.10f") - return components_file - - def rename(in_files, suffix=None): from nipype.utils.filemanip import (filename_to_list, split_filename, list_to_filename) @@ -592,7 +548,7 @@ def create_workflow(files, realign.inputs.slice_info = 2 realign.plugin_args = {'sbatch_args': '-c%d' % 4} - # Comute TSNR on realigned data regressing polynomials upto order 2 + # Compute TSNR on realigned data regressing polynomials up to order 2 tsnr = MapNode(TSNR(regress_poly=2), iterfield=['in_file'], name='tsnr') wf.connect(realign, "out_file", tsnr, "in_file") @@ -694,14 +650,10 @@ def merge_files(in1, in2): filter1, 'out_res_name') wf.connect(createfilter1, 'out_files', filter1, 'design') - createfilter2 = MapNode(Function(input_names=['realigned_file', 'mask_file', - 'num_components', - 'extra_regressors'], - output_names=['out_files'], - function=extract_noise_components, - imports=imports), + createfilter2 = MapNode(ACompCor(), iterfield=['realigned_file', 'extra_regressors'], name='makecompcorrfilter') + createfilter2.inputs.components_file('noise_components.txt') createfilter2.inputs.num_components = num_components wf.connect(createfilter1, 'out_files', createfilter2, 'extra_regressors') @@ -717,7 +669,7 @@ def merge_files(in1, in2): wf.connect(filter1, 'out_res', filter2, 'in_file') wf.connect(filter1, ('out_res', rename, '_cleaned'), filter2, 'out_res_name') - wf.connect(createfilter2, 'out_files', filter2, 'design') + wf.connect(createfilter2, 'components_file', filter2, 'design') wf.connect(mask, 'mask_file', filter2, 'mask') bandpass = Node(Function(input_names=['files', 'lowpass_freq', @@ -923,7 +875,7 @@ def get_names(files, suffix): wf.connect(smooth, 'out_file', datasink, 'resting.timeseries.@smoothed') wf.connect(createfilter1, 'out_files', datasink, 'resting.regress.@regressors') - wf.connect(createfilter2, 'out_files', + wf.connect(createfilter2, 'components_file', datasink, 'resting.regress.@compcorr') wf.connect(maskts, 'out_file', datasink, 'resting.timeseries.target') wf.connect(sampleaparc, 'summary_file', From ea607824166da3ebbf00a3243be5410940e07073 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Mon, 12 Sep 2016 03:18:54 +0000 Subject: [PATCH 55/61] forgot to regress pre-computing std --- nipype/algorithms/compcor.py | 7 +++++-- nipype/algorithms/tests/test_compcor.py | 6 +++--- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index 601093f4fa..4687e533b0 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -118,13 +118,16 @@ class TCompCor(CompCor): def _run_interface(self, runtime): imgseries = nb.load(self.inputs.realigned_file).get_data() - time_voxels = imgseries.T - num_voxels = np.prod(time_voxels.shape[1:]) # 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) + imgseries = imgseries - np.mean(imgseries, axis=1)[:, np.newaxis] + + time_voxels = imgseries.T + num_voxels = np.prod(time_voxels.shape[1:]) # "To construct the tSTD noise ROI, we sorted the voxels by their # temporal standard deviation ..." diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index d8a6e74c02..38ecde900c 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -61,9 +61,9 @@ def test_compcor_with_extra_regressors(self, mock_add_extras): def test_tcompcor(self): ccinterface = TCompCor(realigned_file=self.realigned_file) - self.run_cc(ccinterface, [['-0.2846272268'], ['0.7115680670'], - ['-0.6048328569'], ['0.2134704201'], - ['-0.0355784033']]) + self.run_cc(ccinterface, [['-0.2500000000'], ['0.2500000000'], + ['-0.2500000000'], ['0.7500000000'], + ['-0.5000000000']]) def test_tcompcor_with_percentile(self): ccinterface = TCompCor(realigned_file=self.realigned_file, percentile_threshold=0.2) From 6a1497de3af995999ba6bcde9cdedb18becb06fa Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Mon, 12 Sep 2016 03:25:06 +0000 Subject: [PATCH 56/61] use relative paths in doctests --- nipype/algorithms/compcor.py | 31 +++++++++++++++++++++++++++++-- 1 file changed, 29 insertions(+), 2 deletions(-) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index 4687e533b0..5be6a72c7a 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -1,3 +1,16 @@ +# -*- 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: +''' +Miscellaneous algorithms + + 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 ..interfaces.base import (BaseInterfaceInputSpec, TraitedSpec, BaseInterface, traits, File) from ..pipeline import engine as pe @@ -38,8 +51,8 @@ class CompCor(BaseInterface): ------- >>> ccinterface = CompCor() - >>> ccinterface.inputs.realigned_file = 'nipype/testing/data/functional.nii' - >>> ccinterface.inputs.mask_file = 'nipype/testing/data/mask.nii' + >>> ccinterface.inputs.realigned_file = 'functional.nii' + >>> ccinterface.inputs.mask_file = 'mask.nii' >>> ccinterface.inputs.num_components = 1 >>> ccinterface.inputs.use_regress_poly = True >>> ccinterface.inputs.regress_poly_degree = 2 @@ -112,6 +125,20 @@ class TCompCorInputSpec(CompCorInputSpec): 'with the highest variance are used.') class TCompCor(CompCor): + ''' + Interface for tCompCor. Computes a ROI mask based on variance of voxels. + + Example + ------- + + >>> ccinterface = TCompCor() + >>> ccinterface.inputs.realigned_file = 'functional.nii' + >>> ccinterface.inputs.mask_file = 'mask.nii' + >>> ccinterface.inputs.num_components = 1 + >>> ccinterface.inputs.use_regress_poly = True + >>> ccinterface.inputs.regress_poly_degree = 2 + >>> ccinterface.inputs.percentile_threshold = .03 + ''' input_spec = TCompCorInputSpec output_spec = CompCorOutputSpec From e245073b6c818c7508b6e236f6aba8ad4ce2db6a Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Mon, 12 Sep 2016 03:48:27 +0000 Subject: [PATCH 57/61] update documentation --- examples/README | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/README b/examples/README index 0f0159013e..3e482bb975 100644 --- a/examples/README +++ b/examples/README @@ -1,4 +1,4 @@ A dataset for use with these scripts can be downloaded from the nipype website. At the time of writing, it's at: -http://nipy.org/nipype/users/pipeline_tutorial.html \ No newline at end of file +http://nipype.readthedocs.io/en/0.12.0/users/pipeline_tutorial.html From 6a88b3a8fa8a988d6ce99195bf249622979cbe10 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Mon, 12 Sep 2016 05:55:17 +0000 Subject: [PATCH 58/61] robuster tests --- nipype/algorithms/tests/test_compcor.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index 38ecde900c..4beb756bba 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -60,18 +60,20 @@ def test_compcor_with_extra_regressors(self, mock_add_extras): assert_true(mock_add_extras.called) def test_tcompcor(self): + ccinterface = TCompCor(realigned_file=self.realigned_file, percentile_threshold=0.75) + self.run_cc(ccinterface, [['-0.1535587949', '-0.4318584065'], + ['0.4286265207', '0.7162580102'], + ['-0.6288093202', '0.0465452630'], + ['0.5859742580', '-0.5144309306'], + ['-0.2322326636', '0.1834860639']]) + + def test_tcompcor_with_default_percentile(self): ccinterface = TCompCor(realigned_file=self.realigned_file) - self.run_cc(ccinterface, [['-0.2500000000'], ['0.2500000000'], - ['-0.2500000000'], ['0.7500000000'], - ['-0.5000000000']]) - - def test_tcompcor_with_percentile(self): - ccinterface = TCompCor(realigned_file=self.realigned_file, percentile_threshold=0.2) ccinterface.run() mask = nb.load('mask.nii').get_data() num_nonmasked_voxels = np.count_nonzero(mask) - assert_equal(num_nonmasked_voxels, 2) + assert_equal(num_nonmasked_voxels, 1) def run_cc(self, ccinterface, expected_components): # run From 53478811b8dfbad2f93a36fdcc788e8029173f0a Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Mon, 12 Sep 2016 06:28:08 +0000 Subject: [PATCH 59/61] fix dumb error, comments --- examples/rsfmri_vol_surface_preprocessing_nipy.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/rsfmri_vol_surface_preprocessing_nipy.py b/examples/rsfmri_vol_surface_preprocessing_nipy.py index 4f5bae4c8d..b745704023 100644 --- a/examples/rsfmri_vol_surface_preprocessing_nipy.py +++ b/examples/rsfmri_vol_surface_preprocessing_nipy.py @@ -52,7 +52,9 @@ from nipype.interfaces.base import CommandLine CommandLine.set_default_terminal_output('allatonce') +# https://github.com/moloney/dcmstack from dcmstack.extract import default_extractor +# pip install pydicom from dicom import read_file from nipype.interfaces import (fsl, Function, ants, freesurfer, nipy) @@ -653,7 +655,7 @@ def merge_files(in1, in2): createfilter2 = MapNode(ACompCor(), iterfield=['realigned_file', 'extra_regressors'], name='makecompcorrfilter') - createfilter2.inputs.components_file('noise_components.txt') + createfilter2.inputs.components_file = 'noise_components.txt' createfilter2.inputs.num_components = num_components wf.connect(createfilter1, 'out_files', createfilter2, 'extra_regressors') From cb2faa2ec6d7ef212d9bedacb104c05080c109d1 Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Tue, 13 Sep 2016 05:25:19 +0000 Subject: [PATCH 60/61] clarify documentation --- nipype/algorithms/compcor.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index 5be6a72c7a..68201f234c 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -121,7 +121,10 @@ class TCompCorInputSpec(CompCorInputSpec): exclude_low=True, exclude_high=True, usedefault=True, desc='the percentile ' 'used to select highest-variance ' - 'voxels. By default the 2% of voxels ' + 'voxels, represented by a number ' + 'between 0 and 1, exclusive. By ' + 'default, this value is set to .02. ' + 'That is, the 2% of voxels ' 'with the highest variance are used.') class TCompCor(CompCor): From bb45785c506c801f9c9c83548efb086a3c8349ee Mon Sep 17 00:00:00 2001 From: Shoshana Berleant Date: Tue, 13 Sep 2016 18:41:13 +0000 Subject: [PATCH 61/61] remove stuff --- nipype/algorithms/compcor.py | 11 ----------- nipype/algorithms/tests/test_compcor.py | 9 --------- 2 files changed, 20 deletions(-) diff --git a/nipype/algorithms/compcor.py b/nipype/algorithms/compcor.py index 68201f234c..a2f2f434b3 100644 --- a/nipype/algorithms/compcor.py +++ b/nipype/algorithms/compcor.py @@ -36,8 +36,6 @@ class CompCorInputSpec(BaseInterfaceInputSpec): 'pre-component extraction') regress_poly_degree = traits.Range(low=1, default=1, usedefault=True, desc='the degree polynomial to use') - extra_regressors = File(exists=True, mandatory=False, - desc='additional regressors to add') class CompCorOutputSpec(TraitedSpec): components_file = File(exists=True, @@ -90,10 +88,6 @@ def _run_interface(self, runtime): # principal components using a singular value decomposition." u, _, _ = linalg.svd(M, full_matrices=False) components = u[:, :self.inputs.num_components] - if self.inputs.extra_regressors: - components = self._add_extras(components, - self.inputs.extra_regressors) - components_file = os.path.join(os.getcwd(), self.inputs.components_file) np.savetxt(components_file, components, fmt="%.10f") return runtime @@ -108,13 +102,8 @@ def _compute_tSTD(self, M, x): # set bad values to x stdM[stdM == 0] = x stdM[np.isnan(stdM)] = x - stdM[np.isinf(stdM)] = x return stdM - def _add_extras(self, components, extra_regressors): - regressors = np.genfromtxt(self.inputs.extra_regressors) - return np.hstack((components, regressors)) - class TCompCorInputSpec(CompCorInputSpec): # and all the fields in CompCorInputSpec percentile_threshold = traits.Range(low=0., high=1., value=.02, diff --git a/nipype/algorithms/tests/test_compcor.py b/nipype/algorithms/tests/test_compcor.py index 4beb756bba..256c258a17 100644 --- a/nipype/algorithms/tests/test_compcor.py +++ b/nipype/algorithms/tests/test_compcor.py @@ -5,7 +5,6 @@ from ..compcor import CompCor, TCompCor, ACompCor import unittest -import mock import nibabel as nb import numpy as np import os @@ -51,14 +50,6 @@ def test_compcor(self): assert_equal(os.path.getsize(ccresult.outputs.components_file), os.path.getsize(accresult.outputs.components_file)) - @mock.patch('nipype.algorithms.compcor.CompCor._add_extras') - def test_compcor_with_extra_regressors(self, mock_add_extras): - regressors_file ='regress.txt' - open(regressors_file, 'a').close() # make sure file exists - CompCor(realigned_file=self.realigned_file, mask_file=self.mask_file, - extra_regressors=regressors_file).run() - assert_true(mock_add_extras.called) - def test_tcompcor(self): ccinterface = TCompCor(realigned_file=self.realigned_file, percentile_threshold=0.75) self.run_cc(ccinterface, [['-0.1535587949', '-0.4318584065'],