Merge pull request #1651 from shoshber/regresspoly

ENH Make regresspoly more universal

Author: satra

nipype/algorithms/confounds.py

@@ -12,18 +12,18 @@
 
 '''
 from __future__ import print_function, division, unicode_literals, absolute_import
-from builtins import str, zip, range, open
+from builtins import range
 
 import os
 import os.path as op
 
 import nibabel as nb
 import numpy as np
-from scipy import linalg
-from scipy.special import legendre
+from numpy.polynomial import Legendre
+from scipy import linalg, signal
 
 from .. import logging
-from ..external.due import due, Doi, BibTeX
+from ..external.due import BibTeX
 from ..interfaces.base import (traits, TraitedSpec, BaseInterface,
                                BaseInterfaceInputSpec, File, isdefined,
                                InputMultiPath)
@@ -160,8 +160,8 @@ def _run_interface(self, runtime):
                     'dvars_nstd', ext=self.inputs.figformat)
                 fig = plot_confound(dvars[1], self.inputs.figsize, 'DVARS', series_tr=tr)
                 fig.savefig(self._results['fig_nstd'], dpi=float(self.inputs.figdpi),
-                        format=self.inputs.figformat,
-                        bbox_inches='tight')
+                            format=self.inputs.figformat,
+                            bbox_inches='tight')
                 fig.clf()
 
         if self.inputs.save_vxstd:
@@ -175,8 +175,8 @@ def _run_interface(self, runtime):
                 fig = plot_confound(dvars[2], self.inputs.figsize, 'Voxelwise std DVARS',
                                     series_tr=tr)
                 fig.savefig(self._results['fig_vxstd'], dpi=float(self.inputs.figdpi),
-                        format=self.inputs.figformat,
-                        bbox_inches='tight')
+                            format=self.inputs.figformat,
+                            bbox_inches='tight')
                 fig.clf()
 
         if self.inputs.save_all:
@@ -323,9 +323,9 @@ class CompCor(BaseInterface):
                                     "author = {Behzadi, Yashar and Restom, Khaled and Liau, Joy and Liu, Thomas T.},"
                                     "year = {2007},"
                                     "pages = {90-101},}"
-                                ),
+                                   ),
                     'tags': ['method', 'implementation']
-                }]
+                   }]
 
     def _run_interface(self, runtime):
         imgseries = nb.load(self.inputs.realigned_file).get_data()
@@ -337,18 +337,13 @@ def _run_interface(self, runtime):
         # from paper:
         # "The constant and linear trends of the columns in the matrix M were
         # removed [prior to ...]"
-        if self.inputs.use_regress_poly:
-            voxel_timecourses = regress_poly(self.inputs.regress_poly_degree,
-                                             voxel_timecourses)
-        voxel_timecourses = voxel_timecourses - np.mean(voxel_timecourses,
-                                                        axis=1)[:, np.newaxis]
+        degree = self.inputs.regress_poly_degree if self.inputs.use_regress_poly else 0
+        voxel_timecourses = regress_poly(degree, voxel_timecourses)
 
         # "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 = voxel_timecourses.T
-        numvols = M.shape[0]
-        numvoxels = M.shape[1]
 
         # "[... were removed] prior to column-wise variance normalization."
         M = M / self._compute_tSTD(M, 1.)
@@ -412,7 +407,6 @@ def _run_interface(self, runtime):
         # 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:])
@@ -488,7 +482,7 @@ def _run_interface(self, runtime):
             data = data.astype(np.float32)
 
         if isdefined(self.inputs.regress_poly):
-            data = regress_poly(self.inputs.regress_poly, data)
+            data = regress_poly(self.inputs.regress_poly, data, remove_mean=False)
             img = nb.Nifti1Image(data, img.get_affine(), header)
             nb.save(img, op.abspath(self.inputs.detrended_file))
 
@@ -513,28 +507,33 @@ def _list_outputs(self):
             outputs['detrended_file'] = op.abspath(self.inputs.detrended_file)
         return outputs
 
-def regress_poly(degree, data):
+def regress_poly(degree, data, remove_mean=True, axis=-1):
     ''' returns data with degree polynomial regressed out.
-    The last dimension (i.e. data.shape[-1]) should be time.
+    Be default it is calculated along the last axis (usu. time).
+    If remove_mean is True (default), the data is demeaned (i.e. degree 0).
+    If remove_mean is false, the data is not.
     '''
     datashape = data.shape
-    timepoints = datashape[-1]
+    timepoints = datashape[axis]
 
     # Rearrange all voxel-wise time-series in rows
     data = data.reshape((-1, timepoints))
 
     # Generate design matrix
-    X = np.ones((timepoints, 1))
+    X = np.ones((timepoints, 1)) # quick way to calc degree 0
     for i in range(degree):
-        polynomial_func = legendre(i+1)
+        polynomial_func = Legendre.basis(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
+    if remove_mean:
+        datahat = X.dot(betas).T
+    else: # disregard the first layer of X, which is degree 0
+        datahat = X[:, 1:].dot(betas[1:, ...]).T
     regressed_data = data - datahat
 
     # Back to original shape
@@ -569,7 +568,6 @@ def compute_dvars(in_file, in_mask, remove_zerovariance=False):
     :return: the standardized DVARS
 
     """
-    import os.path as op
     import numpy as np
     import nibabel as nb
     from nitime.algorithms import AR_est_YW
@@ -594,7 +592,7 @@ def compute_dvars(in_file, in_mask, remove_zerovariance=False):
     mfunc = func[idx[0], idx[1], idx[2], :]
 
     # Demean
-    mfunc -= mfunc.mean(axis=1).astype(np.float32)[..., np.newaxis]
+    mfunc = regress_poly(0, mfunc, remove_mean=True).astype(np.float32)
 
     # Compute (non-robust) estimate of lag-1 autocorrelation
     ar1 = np.apply_along_axis(AR_est_YW, 1, mfunc, 1)[:, 0]

nipype/algorithms/tests/test_compcor.py

@@ -1,15 +1,15 @@
 # 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, assert_false, skipif, utils
-from ..confounds import CompCor, TCompCor, ACompCor
-
-import unittest
-import nibabel as nb
-import numpy as np
 import os
 import tempfile
 import shutil
+import unittest
+
+import nibabel as nb
+import numpy as np
+
+from ...testing import assert_equal, assert_true, utils
+from ..confounds import CompCor, TCompCor, ACompCor
 
 class TestCompCor(unittest.TestCase):
     ''' Note: Tests currently do a poor job of testing functionality '''
@@ -27,8 +27,8 @@ def setUp(self):
         self.realigned_file = utils.save_toy_nii(self.fake_data + noise,
                                                  self.filenames['functionalnii'])
         mask = np.ones(self.fake_data.shape[:3])
-        mask[0,0,0] = 0
-        mask[0,0,1] = 0
+        mask[0, 0, 0] = 0
+        mask[0, 0, 1] = 0
         self.mask_file = utils.save_toy_nii(mask, self.filenames['masknii'])
 
     def test_compcor(self):
@@ -52,20 +52,28 @@ def test_compcor(self):
 
     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']])
+        self.run_cc(ccinterface, [['-0.1114536190', '-0.4632908609'],
+                                  ['0.4566907310', '0.6983205193'],
+                                  ['-0.7132557407', '0.1340170559'],
+                                  ['0.5022537643', '-0.5098322262'],
+                                  ['-0.1342351356', '0.1407855119']])
 
-    def test_tcompcor_with_default_percentile(self):
+    def test_tcompcor_no_percentile(self):
         ccinterface = TCompCor(realigned_file=self.realigned_file)
         ccinterface.run()
 
         mask = nb.load('mask.nii').get_data()
         num_nonmasked_voxels = np.count_nonzero(mask)
         assert_equal(num_nonmasked_voxels, 1)
 
+    def test_compcor_no_regress_poly(self):
+        self.run_cc(CompCor(realigned_file=self.realigned_file, mask_file=self.mask_file,
+                            use_regress_poly=False), [['0.4451946442', '-0.7683311482'],
+                                                      ['-0.4285129505', '-0.0926034137'],
+                                                      ['0.5721540256', '0.5608764842'],
+                                                      ['-0.5367548139', '0.0059943226'],
+                                                      ['-0.0520809054', '0.2940637551']])
+
     def run_cc(self, ccinterface, expected_components):
         # run
         ccresult = ccinterface.run()