@@ -50,6 +50,11 @@ class ComputeDVARSInputSpec(BaseInterfaceInputSpec):
5050 remove_zerovariance = traits .Bool (
5151 True , usedefault = True , desc = "remove voxels with zero variance"
5252 )
53+ variance_tol = traits .Float (
54+ 1e-7 ,
55+ usedefault = True ,
56+ desc = "maximum variance to consider \" close to\" zero for the purposes of removal" ,
57+ )
5358 save_std = traits .Bool (True , usedefault = True , desc = "save standardized DVARS" )
5459 save_nstd = traits .Bool (False , usedefault = True , desc = "save non-standardized DVARS" )
5560 save_vxstd = traits .Bool (
@@ -167,6 +172,7 @@ def _run_interface(self, runtime):
167172 self .inputs .in_file ,
168173 self .inputs .in_mask ,
169174 remove_zerovariance = self .inputs .remove_zerovariance ,
175+ variance_tol = self .inputs .variance_tol ,
170176 intensity_normalization = self .inputs .intensity_normalization ,
171177 )
172178
@@ -994,8 +1000,19 @@ def _list_outputs(self):
9941000 return self ._results
9951001
9961002
1003+ def _AR_est_YW (x , order , rxx = None ):
1004+ """Retrieve AR coefficients while dropping the sig_sq return value"""
1005+ from nitime .algorithms import AR_est_YW
1006+ 1007+ return AR_est_YW (x , order , rxx = rxx )[0 ]
1008+ 1009+ 9971010def compute_dvars (
998- in_file , in_mask , remove_zerovariance = False , intensity_normalization = 1000
1011+ in_file ,
1012+ in_mask ,
1013+ remove_zerovariance = False ,
1014+ intensity_normalization = 1000 ,
1015+ variance_tol = 0.0 ,
9991016):
10001017 """
10011018 Compute the :abbr:`DVARS (D referring to temporal
@@ -1027,36 +1044,35 @@ def compute_dvars(
10271044 """
10281045 import numpy as np
10291046 import nibabel as nb
1030- from nitime .algorithms import AR_est_YW
10311047 import warnings
10321048
1033- func = nb .load (in_file ).get_fdata ( dtype = np . float32 )
1034- mask = np .asanyarray (nb .load (in_mask ).dataobj ). astype ( np . uint8 )
1049+ func = np . float32 ( nb .load (in_file ).dataobj )
1050+ mask = np .bool_ (nb .load (in_mask ).dataobj )
10351051
10361052 if len (func .shape ) != 4 :
10371053 raise RuntimeError ("Input fMRI dataset should be 4-dimensional" )
10381054
1039- idx = np .where (mask > 0 )
1040- mfunc = func [idx [0 ], idx [1 ], idx [2 ], :]
1055+ mfunc = func [mask ]
10411056
10421057 if intensity_normalization != 0 :
10431058 mfunc = (mfunc / np .median (mfunc )) * intensity_normalization
10441059
10451060 # Robust standard deviation (we are using "lower" interpolation
10461061 # because this is what FSL is doing
10471062 func_sd = (
1048- np .percentile (mfunc , 75 , axis = 1 , interpolation = "lower" )
1049- - np .percentile (mfunc , 25 , axis = 1 , interpolation = "lower" )
1063+ np .percentile (mfunc , 75 , axis = 1 , method = "lower" )
1064+ - np .percentile (mfunc , 25 , axis = 1 , method = "lower" )
10501065 ) / 1.349
10511066
10521067 if remove_zerovariance :
1053- mfunc = mfunc [func_sd != 0 , :]
1054- func_sd = func_sd [func_sd != 0 ]
1068+ zero_variance_voxels = func_sd > variance_tol
1069+ mfunc = mfunc [zero_variance_voxels , :]
1070+ func_sd = func_sd [zero_variance_voxels ]
10551071
10561072 # Compute (non-robust) estimate of lag-1 autocorrelation
10571073 ar1 = np .apply_along_axis (
1058- AR_est_YW , 1 , regress_poly (0 , mfunc , remove_mean = True )[0 ].astype (np .float32 ), 1
1059- )[:, 0 ]
1074+ _AR_est_YW , 1 , regress_poly (0 , mfunc , remove_mean = True )[0 ].astype (np .float32 ), 1
1075+ )
10601076
10611077 # Compute (predicted) standard deviation of temporal difference time series
10621078 diff_sdhat = np .squeeze (np .sqrt (((1 - ar1 ) * 2 ).tolist ())) * func_sd
0 commit comments