move TSNR to confounds

Author: Shoshana Berleant

nipype/algorithms/confounds.py

@@ -20,12 +20,13 @@
 import nibabel as nb
 import numpy as np
 from scipy import linalg
+from scipy.special import legendre
 
 from .. import logging
 from ..external.due import due, Doi, BibTeX
 from ..interfaces.base import (traits, TraitedSpec, BaseInterface,
-                               BaseInterfaceInputSpec, File, isdefined)
-from .misc import regress_poly
+                               BaseInterfaceInputSpec, File, isdefined,
+                               InputMultiPath)
 IFLOG = logging.getLogger('interface')
 
 class ComputeDVARSInputSpec(BaseInterfaceInputSpec):
@@ -424,6 +425,108 @@ def _run_interface(self, runtime):
 
 ACompCor = CompCor
 
+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)
+
 def compute_dvars(in_file, in_mask, remove_zerovariance=False):
     """
     Compute the :abbr:`DVARS (D referring to temporal

nipype/algorithms/misc.py

@@ -22,7 +22,6 @@
 import numpy as np
 from math import floor, ceil
 from scipy.ndimage.morphology import grey_dilation
-from scipy.special import legendre
 import scipy.io as sio
 import itertools
 import scipy.stats as stats
@@ -35,6 +34,7 @@
                                BaseInterfaceInputSpec, isdefined,
                                DynamicTraitedSpec, Undefined)
 from ..utils.filemanip import fname_presuffix, split_filename
+from . import confounds
 
 iflogger = logging.getLogger('interface')
 
@@ -258,108 +258,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):
-            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)
 
@@ -1520,3 +1418,14 @@ def __init__(self, **inputs):
         warnings.warn(("This interface has been deprecated since 0.10.0,"
                        " please use nipype.algorithms.metrics.FuzzyOverlap"),
                       DeprecationWarning)
+
+class TSNR(confounds.TSNR):
+    """
+    .. deprecated:: 0.12.1
+       Use :py:class:`nipype.algorithms.confounds.TSNR` instead
+    """
+    def __init__(self, **inputs):
+        super(confounds.TSNR, self).__init__(**inputs)
+        warnings.warn(("This interface has been deprecated since 0.12.0,"
+                       " please use nipype.algorithms.confounds.TSNR"),
+                      DeprecationWarning)

nipype/algorithms/tests/test_tsnr.py

@@ -3,7 +3,8 @@
 
 from ...testing import (assert_equal, assert_true, assert_almost_equal,
                         skipif, utils)
-from ..misc import TSNR
+from ..confounds import TSNR
+from .. import misc
 
 import unittest
 import mock
@@ -85,6 +86,14 @@ def test_tsnr_withpoly3(self):
             'tsnr_file': (2.6, 57.3)
         })
 
+    @mock.patch('warnings.warn')
+    def test_deprecation_warning(self, mock_warn):
+        # run
+        misc.TSNR(in_file=self.in_filenames['in_file'])
+
+        # assert
+        mock_warn.assert_called_once_with(mock.ANY, DeprecationWarning)
+
     def assert_expected_outputs_poly(self, tsnrresult, expected_ranges):
         assert_equal(os.path.basename(tsnrresult.outputs.detrended_file),
                      self.out_filenames['detrended_file'])