ENH CompCor

examples/README

@@ -1,4 +1,4 @@
 A dataset for use with these scripts can be downloaded from the nipype
 website. At the time of writing, it's at:
 
-http://nipy.org/nipype/users/pipeline_tutorial.html
+http://nipype.readthedocs.io/en/0.12.0/users/pipeline_tutorial.html

examples/rsfmri_vol_surface_preprocessing_nipy.py

@@ -52,7 +52,9 @@
 from nipype.interfaces.base import CommandLine
 CommandLine.set_default_terminal_output('allatonce')
 
+# https://github.com/moloney/dcmstack
 from dcmstack.extract import default_extractor
+# pip install pydicom
 from dicom import read_file
 
 from nipype.interfaces import (fsl, Function, ants, freesurfer, nipy)
@@ -64,6 +66,7 @@
 
 from nipype.algorithms.rapidart import ArtifactDetect
 from nipype.algorithms.misc import TSNR
+from nipype.algorithms.compcor import ACompCor
 from nipype.interfaces.utility import Rename, Merge, IdentityInterface
 from nipype.utils.filemanip import filename_to_list
 from nipype.interfaces.io import DataSink, FreeSurferSource
@@ -228,51 +231,6 @@ def build_filter1(motion_params, comp_norm, outliers, detrend_poly=None):
         out_files.append(filename)
     return out_files
 
-
-def extract_noise_components(realigned_file, mask_file, num_components=5,
-                             extra_regressors=None):
-    """Derive components most reflective of physiological noise
-
-    Parameters
-    ----------
-    realigned_file: a 4D Nifti file containing realigned volumes
-    mask_file: a 3D Nifti file containing white matter + ventricular masks
-    num_components: number of components to use for noise decomposition
-    extra_regressors: additional regressors to add
-
-    Returns
-    -------
-    components_file: a text file containing the noise components
-    """
-    imgseries = nb.load(realigned_file)
-    components = None
-    for filename in filename_to_list(mask_file):
-        mask = nb.load(filename).get_data()
-        if len(np.nonzero(mask > 0)[0]) == 0:
-            continue
-        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]))
-    if extra_regressors:
-        regressors = np.genfromtxt(extra_regressors)
-        components = np.hstack((components, regressors))
-    components_file = os.path.join(os.getcwd(), 'noise_components.txt')
-    np.savetxt(components_file, components, fmt=b"%.10f")
-    return components_file
-
-
 def rename(in_files, suffix=None):
     from nipype.utils.filemanip import (filename_to_list, split_filename,
                                         list_to_filename)
@@ -592,7 +550,7 @@ def create_workflow(files,
     realign.inputs.slice_info = 2
     realign.plugin_args = {'sbatch_args': '-c%d' % 4}
 
-    # Comute TSNR on realigned data regressing polynomials upto order 2
+    # Compute TSNR on realigned data regressing polynomials up to order 2
     tsnr = MapNode(TSNR(regress_poly=2), iterfield=['in_file'], name='tsnr')
     wf.connect(realign, "out_file", tsnr, "in_file")
 
@@ -694,14 +652,10 @@ def merge_files(in1, in2):
                filter1, 'out_res_name')
     wf.connect(createfilter1, 'out_files', filter1, 'design')
 
-    createfilter2 = MapNode(Function(input_names=['realigned_file', 'mask_file',
-                                                  'num_components',
-                                                  'extra_regressors'],
-                                     output_names=['out_files'],
-                                     function=extract_noise_components,
-                                     imports=imports),
+    createfilter2 = MapNode(ACompCor(),
                             iterfield=['realigned_file', 'extra_regressors'],
                             name='makecompcorrfilter')
+    createfilter2.inputs.components_file = 'noise_components.txt'
     createfilter2.inputs.num_components = num_components
 
     wf.connect(createfilter1, 'out_files', createfilter2, 'extra_regressors')
@@ -717,7 +671,7 @@ def merge_files(in1, in2):
     wf.connect(filter1, 'out_res', filter2, 'in_file')
     wf.connect(filter1, ('out_res', rename, '_cleaned'),
                filter2, 'out_res_name')
-    wf.connect(createfilter2, 'out_files', filter2, 'design')
+    wf.connect(createfilter2, 'components_file', filter2, 'design')
     wf.connect(mask, 'mask_file', filter2, 'mask')
 
     bandpass = Node(Function(input_names=['files', 'lowpass_freq',
@@ -923,7 +877,7 @@ def get_names(files, suffix):
     wf.connect(smooth, 'out_file', datasink, 'resting.timeseries.@smoothed')
     wf.connect(createfilter1, 'out_files',
                datasink, 'resting.regress.@regressors')
-    wf.connect(createfilter2, 'out_files',
+    wf.connect(createfilter2, 'components_file',
                datasink, 'resting.regress.@compcorr')
     wf.connect(maskts, 'out_file', datasink, 'resting.timeseries.target')
     wf.connect(sampleaparc, 'summary_file',

nipype/algorithms/compcor.py

@@ -0,0 +1,178 @@
+# -*- coding: utf-8 -*-
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+'''
+Miscellaneous algorithms
+
+    Change directory to provide relative paths for doctests
+    >>> import os
+    >>> filepath = os.path.dirname(os.path.realpath(__file__))
+    >>> datadir = os.path.realpath(os.path.join(filepath, '../testing/data'))
+    >>> os.chdir(datadir)
+
+'''
+from ..interfaces.base import (BaseInterfaceInputSpec, TraitedSpec,
+                               BaseInterface, traits, File)
+from ..pipeline import engine as pe
+from ..interfaces.utility import IdentityInterface
+from .misc import regress_poly
+
+import nibabel as nb
+import numpy as np
+from scipy import linalg, stats
+import os
+
+class CompCorInputSpec(BaseInterfaceInputSpec):
+    realigned_file = File(exists=True, mandatory=True,
+                          desc='already realigned brain image (4D)')
+    mask_file = File(exists=True, mandatory=False,
+                     desc='mask file that determines ROI (3D)')
+    components_file = File('components_file.txt', exists=False,
+                           mandatory=False, usedefault=True,
+                           desc='filename to store physiological components')
+    num_components = traits.Int(6, usedefault=True) # 6 for BOLD, 4 for ASL
+    use_regress_poly = traits.Bool(True, usedefault=True,
+                                   desc='use polynomial regression'
+                                   'pre-component extraction')
+    regress_poly_degree = traits.Range(low=1, default=1, usedefault=True,
+                                       desc='the degree polynomial to use')
+    extra_regressors = File(exists=True, mandatory=False,
+                            desc='additional regressors to add')
+
+class CompCorOutputSpec(TraitedSpec):
+    components_file = File(exists=True,
+                           desc='text file containing the noise components')
+
+class CompCor(BaseInterface):
+    '''
+    Interface with core CompCor computation, used in aCompCor and tCompCor
+
+    Example
+    -------
+
+    >>> ccinterface = CompCor()
+    >>> ccinterface.inputs.realigned_file = 'functional.nii'
+    >>> ccinterface.inputs.mask_file = 'mask.nii'
+    >>> ccinterface.inputs.num_components = 1
+    >>> ccinterface.inputs.use_regress_poly = True
+    >>> ccinterface.inputs.regress_poly_degree = 2
+    '''
+    input_spec = CompCorInputSpec
+    output_spec = CompCorOutputSpec
+
+    def _run_interface(self, runtime):
+        imgseries = nb.load(self.inputs.realigned_file).get_data()
+        mask = nb.load(self.inputs.mask_file).get_data()
+        voxel_timecourses = imgseries[mask > 0]
+        # Zero-out any bad values
+        voxel_timecourses[np.isnan(np.sum(voxel_timecourses, axis=1)), :] = 0
+
+        # 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]
+
+        # "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.)
+
+        # "The covariance matrix C = MMT was constructed and decomposed into its
+        # principal components using a singular value decomposition."
+        u, _, _ = linalg.svd(M, full_matrices=False)
+        components = u[:, :self.inputs.num_components]
+        if self.inputs.extra_regressors:
+            components = self._add_extras(components,
+                                          self.inputs.extra_regressors)
+
+        components_file = os.path.join(os.getcwd(), self.inputs.components_file)
+        np.savetxt(components_file, components, fmt="%.10f")
+        return runtime
+
+    def _list_outputs(self):
+        outputs = self._outputs().get()
+        outputs['components_file'] = os.path.abspath(self.inputs.components_file)
+        return outputs
+
+    def _compute_tSTD(self, M, x):
+        stdM = np.std(M, axis=0)
+        # set bad values to x
+        stdM[stdM == 0] = x
+        stdM[np.isnan(stdM)] = x
+        stdM[np.isinf(stdM)] = x
+        return stdM
+
+    def _add_extras(self, components, extra_regressors):
+            regressors = np.genfromtxt(self.inputs.extra_regressors)
+            return np.hstack((components, regressors))
+
+class TCompCorInputSpec(CompCorInputSpec):
+    # and all the fields in CompCorInputSpec
+    percentile_threshold = traits.Range(low=0., high=1., value=.02,
+                                        exclude_low=True, exclude_high=True,
+                                        usedefault=True, desc='the percentile '
+                                        'used to select highest-variance '
+                                        'voxels. By default the 2% of voxels '
+                                        'with the highest variance are used.')
+
+class TCompCor(CompCor):
+    '''
+    Interface for tCompCor. Computes a ROI mask based on variance of voxels.
+
+    Example
+    -------
+
+    >>> ccinterface = TCompCor()
+    >>> ccinterface.inputs.realigned_file = 'functional.nii'
+    >>> ccinterface.inputs.mask_file = 'mask.nii'
+    >>> ccinterface.inputs.num_components = 1
+    >>> ccinterface.inputs.use_regress_poly = True
+    >>> ccinterface.inputs.regress_poly_degree = 2
+    >>> ccinterface.inputs.percentile_threshold = .03
+    '''
+
+    input_spec = TCompCorInputSpec
+    output_spec = CompCorOutputSpec
+
+    def _run_interface(self, runtime):
+        imgseries = nb.load(self.inputs.realigned_file).get_data()
+
+        # From the paper:
+        # "For each voxel time series, the temporal standard deviation is
+        # 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:])
+
+        # "To construct the tSTD noise ROI, we sorted the voxels by their
+        # temporal standard deviation ..."
+        tSTD = self._compute_tSTD(time_voxels, 0)
+        sortSTD = np.sort(tSTD, axis=None) # flattened sorted matrix
+
+        # use percentile_threshold to pick voxels
+        threshold_index = int(num_voxels * (1. - self.inputs.percentile_threshold))
+        threshold_std = sortSTD[threshold_index]
+        mask = tSTD >= threshold_std
+        mask = mask.astype(int)
+
+        # save mask
+        mask_file = 'mask.nii'
+        nb.nifti1.save(nb.Nifti1Image(mask, np.eye(4)), mask_file)
+        self.inputs.mask_file = mask_file
+
+        super(TCompCor, self)._run_interface(runtime)
+        return runtime
+
+ACompCor = CompCor

nipype/algorithms/misc.py

@@ -35,6 +35,7 @@
                                BaseInterfaceInputSpec, isdefined,
                                DynamicTraitedSpec, Undefined)
 from ..utils.filemanip import fname_presuffix, split_filename
+
 iflogger = logging.getLogger('interface')
 
 
@@ -257,7 +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')
@@ -308,17 +308,7 @@ def _run_interface(self, runtime):
             data = data.astype(np.float32)
 
         if isdefined(self.inputs.regress_poly):
-            timepoints = img.shape[-1]
-            X = np.ones((timepoints, 1))
-            for i in range(self.inputs.regress_poly):
-                X = np.hstack((X, legendre(
-                    i + 1)(np.linspace(-1, 1, timepoints))[:, None]))
-            betas = np.dot(np.linalg.pinv(X), np.rollaxis(data, 3, 2))
-            datahat = np.rollaxis(np.dot(X[:, 1:],
-                                         np.rollaxis(
-                                             betas[1:, :, :, :], 0, 3)),
-                                  0, 4)
-            data = data - datahat
+            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))
 
@@ -343,6 +333,32 @@ def _list_outputs(self):
             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)