1212
1313'''
1414from __future__ import print_function , division , unicode_literals , absolute_import
15- from builtins import str , zip , range , open
15+ from builtins import range
1616
1717import os
1818import os .path as op
1919
2020import nibabel as nb
2121import numpy as np
22- from scipy import linalg
23- from scipy . special import legendre
22+ from numpy . polynomial import Legendre
23+ from scipy import linalg , signal
2424
2525from .. import logging
26- from ..external .due import due , Doi , BibTeX
26+ from ..external .due import BibTeX
2727from ..interfaces .base import (traits , TraitedSpec , BaseInterface ,
2828 BaseInterfaceInputSpec , File , isdefined ,
2929 InputMultiPath )
@@ -160,8 +160,8 @@ def _run_interface(self, runtime):
160160 'dvars_nstd' , ext = self .inputs .figformat )
161161 fig = plot_confound (dvars [1 ], self .inputs .figsize , 'DVARS' , series_tr = tr )
162162 fig .savefig (self ._results ['fig_nstd' ], dpi = float (self .inputs .figdpi ),
163- format = self .inputs .figformat ,
164- bbox_inches = 'tight' )
163+ format = self .inputs .figformat ,
164+ bbox_inches = 'tight' )
165165 fig .clf ()
166166
167167 if self .inputs .save_vxstd :
@@ -175,8 +175,8 @@ def _run_interface(self, runtime):
175175 fig = plot_confound (dvars [2 ], self .inputs .figsize , 'Voxelwise std DVARS' ,
176176 series_tr = tr )
177177 fig .savefig (self ._results ['fig_vxstd' ], dpi = float (self .inputs .figdpi ),
178- format = self .inputs .figformat ,
179- bbox_inches = 'tight' )
178+ format = self .inputs .figformat ,
179+ bbox_inches = 'tight' )
180180 fig .clf ()
181181
182182 if self .inputs .save_all :
@@ -192,7 +192,7 @@ def _list_outputs(self):
192192
193193
194194class FramewiseDisplacementInputSpec (BaseInterfaceInputSpec ):
195- in_plots = File (exists = True , desc = 'motion parameters as written by FSL MCFLIRT' )
195+ in_plots = File (exists = True , mandatory = True , desc = 'motion parameters as written by FSL MCFLIRT' )
196196 radius = traits .Float (50 , usedefault = True ,
197197 desc = 'radius in mm to calculate angular FDs, 50mm is the '
198198 'default since it is used in Power et al. 2012' )
@@ -313,6 +313,19 @@ class CompCor(BaseInterface):
313313 '''
314314 input_spec = CompCorInputSpec
315315 output_spec = CompCorOutputSpec
316+ references_ = [{'entry' : BibTeX ("@article{compcor_2007,"
317+ "title = {A component based noise correction method (CompCor) for BOLD and perfusion based},"
318+ "volume = {37},"
319+ "number = {1},"
320+ "doi = {10.1016/j.neuroimage.2007.04.042},"
321+ "urldate = {2016-08-13},"
322+ "journal = {NeuroImage},"
323+ "author = {Behzadi, Yashar and Restom, Khaled and Liau, Joy and Liu, Thomas T.},"
324+ "year = {2007},"
325+ "pages = {90-101},}"
326+ ),
327+ 'tags' : ['method' , 'implementation' ]
328+ }]
316329
317330 def _run_interface (self , runtime ):
318331 imgseries = nb .load (self .inputs .realigned_file ).get_data ()
@@ -324,18 +337,13 @@ def _run_interface(self, runtime):
324337 # from paper:
325338 # "The constant and linear trends of the columns in the matrix M were
326339 # removed [prior to ...]"
327- if self .inputs .use_regress_poly :
328- voxel_timecourses = regress_poly (self .inputs .regress_poly_degree ,
329- voxel_timecourses )
330- voxel_timecourses = voxel_timecourses - np .mean (voxel_timecourses ,
331- axis = 1 )[:, np .newaxis ]
340+ degree = self .inputs .regress_poly_degree if self .inputs .use_regress_poly else 0
341+ voxel_timecourses = regress_poly (degree , voxel_timecourses )
332342
333343 # "Voxel time series from the noise ROI (either anatomical or tSTD) were
334344 # placed in a matrix M of size Nxm, with time along the row dimension
335345 # and voxels along the column dimension."
336346 M = voxel_timecourses .T
337- numvols = M .shape [0 ]
338- numvoxels = M .shape [1 ]
339347
340348 # "[... were removed] prior to column-wise variance normalization."
341349 M = M / self ._compute_tSTD (M , 1. )
@@ -399,7 +407,6 @@ def _run_interface(self, runtime):
399407 # defined as the standard deviation of the time series after the removal
400408 # of low-frequency nuisance terms (e.g., linear and quadratic drift)."
401409 imgseries = regress_poly (2 , imgseries )
402- imgseries = imgseries - np .mean (imgseries , axis = 1 )[:, np .newaxis ]
403410
404411 time_voxels = imgseries .T
405412 num_voxels = np .prod (time_voxels .shape [1 :])
@@ -475,7 +482,7 @@ def _run_interface(self, runtime):
475482 data = data .astype (np .float32 )
476483
477484 if isdefined (self .inputs .regress_poly ):
478- data = regress_poly (self .inputs .regress_poly , data )
485+ data = regress_poly (self .inputs .regress_poly , data , remove_mean = False )
479486 img = nb .Nifti1Image (data , img .get_affine (), header )
480487 nb .save (img , op .abspath (self .inputs .detrended_file ))
481488
@@ -500,28 +507,33 @@ def _list_outputs(self):
500507 outputs ['detrended_file' ] = op .abspath (self .inputs .detrended_file )
501508 return outputs
502509
503- def regress_poly (degree , data ):
510+ def regress_poly (degree , data , remove_mean = True , axis = - 1 ):
504511 ''' returns data with degree polynomial regressed out.
505- The last dimension (i.e. data.shape[-1]) should be time.
512+ Be default it is calculated along the last axis (usu. time).
513+ If remove_mean is True (default), the data is demeaned (i.e. degree 0).
514+ If remove_mean is false, the data is not.
506515 '''
507516 datashape = data .shape
508- timepoints = datashape [- 1 ]
517+ timepoints = datashape [axis ]
509518
510519 # Rearrange all voxel-wise time-series in rows
511520 data = data .reshape ((- 1 , timepoints ))
512521
513522 # Generate design matrix
514- X = np .ones ((timepoints , 1 ))
523+ X = np .ones ((timepoints , 1 )) # quick way to calc degree 0
515524 for i in range (degree ):
516- polynomial_func = legendre ( i + 1 )
525+ polynomial_func = Legendre . basis ( i + 1 )
517526 value_array = np .linspace (- 1 , 1 , timepoints )
518527 X = np .hstack ((X , polynomial_func (value_array )[:, np .newaxis ]))
519528
520529 # Calculate coefficients
521530 betas = np .linalg .pinv (X ).dot (data .T )
522531
523532 # Estimation
524- datahat = X [:, 1 :].dot (betas [1 :, ...]).T
533+ if remove_mean :
534+ datahat = X .dot (betas ).T
535+ else : # disregard the first layer of X, which is degree 0
536+ datahat = X [:, 1 :].dot (betas [1 :, ...]).T
525537 regressed_data = data - datahat
526538
527539 # Back to original shape
@@ -556,7 +568,6 @@ def compute_dvars(in_file, in_mask, remove_zerovariance=False):
556568 :return: the standardized DVARS
557569
558570 """
559- import os .path as op
560571 import numpy as np
561572 import nibabel as nb
562573 from nitime .algorithms import AR_est_YW
@@ -581,7 +592,7 @@ def compute_dvars(in_file, in_mask, remove_zerovariance=False):
581592 mfunc = func [idx [0 ], idx [1 ], idx [2 ], :]
582593
583594 # Demean
584- mfunc -= mfunc . mean ( axis = 1 ).astype (np .float32 )[..., np . newaxis ]
595+ mfunc = regress_poly ( 0 , mfunc , remove_mean = True ).astype (np .float32 )
585596
586597 # Compute (non-robust) estimate of lag-1 autocorrelation
587598 ar1 = np .apply_along_axis (AR_est_YW , 1 , mfunc , 1 )[:, 0 ]
0 commit comments