move tsnr back into misc

Author: Shoshana Berleant

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

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)