Merge pull request #1839 from chrisfilo/enh/nstdstate

Non steady state detector

Author: satra

CHANGES

@@ -1,6 +1,7 @@
 Upcoming Release
 =====================
 
+* ENH: Added non-steady state detector for EPI data (https://github.com/nipy/nipype/pull/1839)
 * ENH: Enable new BBRegister init options for FSv6+ (https://github.com/nipy/nipype/pull/1811)
 * REF: Splits nipype.interfaces.utility into base, csv, and wrappers (https://github.com/nipy/nipype/pull/1828)
 * FIX: Makespec now runs with nipype in current directory (https://github.com/nipy/nipype/pull/1813)

nipype/algorithms/confounds.py

@@ -579,6 +579,79 @@ def _list_outputs(self):
             outputs['detrended_file'] = op.abspath(self.inputs.detrended_file)
         return outputs
 
+
+class NonSteadyStateDetectorInputSpec(BaseInterfaceInputSpec):
+    in_file = File(exists=True, mandatory=True, desc='4D NIFTI EPI file')
+
+
+class NonSteadyStateDetectorOutputSpec(TraitedSpec):
+    n_volumes_to_discard = traits.Int(desc='Number of non-steady state volumes'
+                                           'detected in the beginning of the scan.')
+
+
+class NonSteadyStateDetector(BaseInterface):
+    """
+    Returns the number of non-steady state volumes detected at the beginning
+    of the scan.
+    """
+
+    input_spec = NonSteadyStateDetectorInputSpec
+    output_spec = NonSteadyStateDetectorOutputSpec
+
+    def _run_interface(self, runtime):
+        in_nii = nb.load(self.inputs.in_plots)
+        global_signal = in_nii.get_data()[:,:,:,:50].mean(axis=0).mean(axis=0).mean(axis=0)
+
+        self._results = {
+            'out_file': _is_outlier(global_signal)
+        }
+
+        return runtime
+
+    def _list_outputs(self):
+        return self._results
+
+def _is_outlier(points, thresh=3.5):
+    """
+    Returns a boolean array with True if points are outliers and False
+    otherwise.
+
+    Parameters:
+    -----------
+        points : An numobservations by numdimensions array of observations
+        thresh : The modified z-score to use as a threshold. Observations with
+            a modified z-score (based on the median absolute deviation) greater
+            than this value will be classified as outliers.
+
+    Returns:
+    --------
+        mask : A numobservations-length boolean array.
+
+    References:
+    ----------
+        Boris Iglewicz and David Hoaglin (1993), "Volume 16: How to Detect and
+        Handle Outliers", The ASQC Basic References in Quality Control:
+        Statistical Techniques, Edward F. Mykytka, Ph.D., Editor.
+    """
+    if len(points.shape) == 1:
+        points = points[:, None]
+    median = np.median(points, axis=0)
+    diff = np.sum((points - median) ** 2, axis=-1)
+    diff = np.sqrt(diff)
+    med_abs_deviation = np.median(diff)
+
+    modified_z_score = 0.6745 * diff / med_abs_deviation
+
+    timepoints_to_discard = 0
+    for i in range(len(modified_z_score)):
+        if modified_z_score[i] <= thresh:
+            break
+        else:
+            timepoints_to_discard += 1
+
+    return timepoints_to_discard
+
+
 def regress_poly(degree, data, remove_mean=True, axis=-1):
     ''' returns data with degree polynomial regressed out.
     Be default it is calculated along the last axis (usu. time).

nipype/algorithms/tests/test_confounds.py

@@ -6,7 +6,8 @@
 
 import pytest
 from nipype.testing import example_data
-from nipype.algorithms.confounds import FramewiseDisplacement, ComputeDVARS
+from nipype.algorithms.confounds import FramewiseDisplacement, ComputeDVARS, \
+    _is_outlier
 import numpy as np
 
 
@@ -63,4 +64,11 @@ def test_dvars(tmpdir):
 
     assert (np.abs(dv1[:, 1] - ground_truth[:, 1]).sum() / len(dv1)) > 0.05
 
-    assert (np.abs(dv1[:, 2] - ground_truth[:, 2]).sum() / len(dv1)) < 0.05
+    assert (np.abs(dv1[:, 2] - ground_truth[:, 2]).sum() / len(dv1)) < 0.05
+
+def test_outliers(tmpdir):
+    in_data = np.random.randn(100)
+    in_data[0] += 10
+
+    assert _is_outlier(in_data) == 1
+