@@ -372,8 +372,8 @@ def _list_outputs(self):
372372 outputs ['components_file' ] = os .path .abspath (self .inputs .components_file )
373373 return outputs
374374
375- def _compute_tSTD (self , M , x ):
376- stdM = np .std (M , axis = 0 )
375+ def _compute_tSTD (self , M , x , axis = 0 ):
376+ stdM = np .std (M , axis = axis )
377377 # set bad values to x
378378 stdM [stdM == 0 ] = x
379379 stdM [np .isnan (stdM )] = x
@@ -411,6 +411,10 @@ class TCompCorInputSpec(CompCorInputSpec):
411411 'That is, the 2% of voxels '
412412 'with the highest variance are used.' )
413413
414+ class TCompCorOutputSpec (CompCorInputSpec ):
415+ # and all the fields in CompCorInputSpec
416+ high_variance_mask = File (exists = True , desc = "voxels excedding the variance threshold" )
417+
414418class TCompCor (CompCor ):
415419 '''
416420 Interface for tCompCor. Computes a ROI mask based on variance of voxels.
@@ -428,7 +432,7 @@ class TCompCor(CompCor):
428432 '''
429433
430434 input_spec = TCompCorInputSpec
431- output_spec = CompCorOutputSpec
435+ output_spec = TCompCorOutputSpec
432436
433437 def _run_interface (self , runtime ):
434438 imgseries = nb .load (self .inputs .realigned_file ).get_data ()
@@ -438,29 +442,34 @@ def _run_interface(self, runtime):
438442 '(shape {})'
439443 .format (self .inputs .realigned_file , imgseries .ndim , imgseries .shape ))
440444
445+ if isdefined (self .inputs .mask_file ):
446+ in_mask_data = nb .load (self .inputs .mask_file ).get_data ()
447+ imgseries = imgseries [in_mask_data != 0 , :]
448+
441449 # From the paper:
442450 # "For each voxel time series, the temporal standard deviation is
443451 # defined as the standard deviation of the time series after the removal
444452 # of low-frequency nuisance terms (e.g., linear and quadratic drift)."
445453 imgseries = regress_poly (2 , imgseries )
446454
447- time_voxels = imgseries .T
448- num_voxels = np .prod (time_voxels .shape [1 :])
449-
450455 # "To construct the tSTD noise ROI, we sorted the voxels by their
451456 # temporal standard deviation ..."
452- tSTD = self ._compute_tSTD (time_voxels , 0 )
453- sortSTD = np .sort (tSTD , axis = None ) # flattened sorted matrix
457+ tSTD = self ._compute_tSTD (imgseries , 0 , axis = - 1 )
454458
455459 # use percentile_threshold to pick voxels
456- threshold_index = int (num_voxels * (1. - self .inputs .percentile_threshold ))
457- threshold_std = sortSTD [threshold_index ]
460+ threshold_std = np .percentile (tSTD , 100. * (1. - self .inputs .percentile_threshold ))
458461 mask = tSTD >= threshold_std
459- mask = mask .astype (int ).T
462+
463+ if isdefined (self .inputs .mask_file ):
464+ mask_data = np .zeros_like (in_mask_data )
465+ mask_data [in_mask_data != 0 ] = mask
466+ else :
467+ mask_data = mask .astype (int )
460468
461469 # save mask
462470 mask_file = os .path .abspath ('mask.nii' )
463- nb .nifti1 .save (nb .Nifti1Image (mask , np .eye (4 )), mask_file )
471+ nb .Nifti1Image (mask_data ,
472+ nb .load (self .inputs .realigned_file ).affine ).to_filename (mask_file )
464473 IFLOG .debug ('tCompcor computed and saved mask of shape {} to mask_file {}'
465474 .format (mask .shape , mask_file ))
466475 self .inputs .mask_file = mask_file
@@ -469,6 +478,11 @@ def _run_interface(self, runtime):
469478 super (TCompCor , self )._run_interface (runtime )
470479 return runtime
471480
481+ def _list_outputs (self ):
482+ outputs = super (TCompCor , self )._list_outputs ()
483+ outputs ['high_variance_mask' ] = self .inputs .mask_file
484+ return outputs
485+
472486class TSNRInputSpec (BaseInterfaceInputSpec ):
473487 in_file = InputMultiPath (File (exists = True ), mandatory = True ,
474488 desc = 'realigned 4D file or a list of 3D files' )
@@ -654,7 +668,7 @@ def compute_dvars(in_file, in_mask, remove_zerovariance=False):
654668 warnings .filterwarnings ('error' )
655669
656670 # voxelwise standardization
657- diff_vx_stdz = func_diff / np .array ([diff_sdhat ] * func_diff .shape [- 1 ]).T
671+ diff_vx_stdz = func_diff / np .array ([diff_sdhat ] * func_diff .shape [- 1 ]).T
658672 dvars_vx_stdz = diff_vx_stdz .std (axis = 0 , ddof = 1 )
659673
660674 return (dvars_stdz , dvars_nstd , dvars_vx_stdz )
0 commit comments