refactor rsfmri fsl workflow

Author: Shoshana Berleant

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]
 

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')