diff --git a/CHANGES b/CHANGES index 52b2c8f00b..1583854034 100644 --- a/CHANGES +++ b/CHANGES @@ -1,6 +1,8 @@ Upcoming release 0.13 ===================== +* ENH: DVARS includes intensity normalization feature - turned on by default (https://github.com/nipy/nipype/pull/1827) +* FIX: DVARS is correctly using sum of squares instead of standard deviation (https://github.com/nipy/nipype/pull/1827) * ENH: Refactoring of nipype.interfaces.utility (https://github.com/nipy/nipype/pull/1828) * FIX: CircleCI were failing silently. Some fixes to tests (https://github.com/nipy/nipype/pull/1833) * FIX: Issues in Docker image permissions, and docker documentation (https://github.com/nipy/nipype/pull/1825) diff --git a/nipype/algorithms/confounds.py b/nipype/algorithms/confounds.py index c6503c703a..db616f1f04 100644 --- a/nipype/algorithms/confounds.py +++ b/nipype/algorithms/confounds.py @@ -52,6 +52,16 @@ class ComputeDVARSInputSpec(BaseInterfaceInputSpec): desc='output figure size') figformat = traits.Enum('png', 'pdf', 'svg', usedefault=True, desc='output format for figures') + intensity_normalization = traits.Float(1000.0, usedefault=True, + desc='Divide value in each voxel at each timepoint ' + 'by the median calculated across all voxels' + 'and timepoints within the mask (if specified)' + 'and then multiply by the value specified by' + 'this parameter. By using the default (1000)' \ + 'output DVARS will be expressed in ' \ + 'x10 % BOLD units compatible with Power et al.' \ + '2012. Set this to 0 to disable intensity' \ + 'normalization altogether.') @@ -128,7 +138,8 @@ def _gen_fname(self, suffix, ext=None): def _run_interface(self, runtime): dvars = compute_dvars(self.inputs.in_file, self.inputs.in_mask, - remove_zerovariance=self.inputs.remove_zerovariance) + remove_zerovariance=self.inputs.remove_zerovariance, + intensity_normalization=self.inputs.intensity_normalization) (self._results['avg_std'], self._results['avg_nstd'], @@ -595,7 +606,8 @@ def regress_poly(degree, data, remove_mean=True, axis=-1): # Back to original shape return regressed_data.reshape(datashape) -def compute_dvars(in_file, in_mask, remove_zerovariance=False): +def compute_dvars(in_file, in_mask, remove_zerovariance=False, + intensity_normalization=1000): """ Compute the :abbr:`DVARS (D referring to temporal derivative of timecourses, VARS referring to RMS variance over voxels)` @@ -636,59 +648,49 @@ def compute_dvars(in_file, in_mask, remove_zerovariance=False): raise RuntimeError( "Input fMRI dataset should be 4-dimensional") - # Robust standard deviation - func_sd = (np.percentile(func, 75, axis=3) - - np.percentile(func, 25, axis=3)) / 1.349 - func_sd[mask <= 0] = 0 - - if remove_zerovariance: - # Remove zero-variance voxels across time axis - mask = zero_remove(func_sd, mask) - idx = np.where(mask > 0) mfunc = func[idx[0], idx[1], idx[2], :] - # Demean - mfunc = regress_poly(0, mfunc, remove_mean=True).astype(np.float32) + if intensity_normalization != 0: + mfunc = (mfunc / np.median(mfunc)) * intensity_normalization + + # Robust standard deviation (we are using "lower" interpolation + # because this is what FSL is doing + func_sd = (np.percentile(mfunc, 75, axis=1, interpolation="lower") - + np.percentile(mfunc, 25, axis=1, interpolation="lower")) / 1.349 + + if remove_zerovariance: + mfunc = mfunc[func_sd != 0, :] + func_sd = func_sd[func_sd != 0] # Compute (non-robust) estimate of lag-1 autocorrelation - ar1 = np.apply_along_axis(AR_est_YW, 1, mfunc, 1)[:, 0] + ar1 = np.apply_along_axis(AR_est_YW, 1, + regress_poly(0, mfunc, remove_mean=True).astype( + np.float32), 1)[:, 0] # Compute (predicted) standard deviation of temporal difference time series - diff_sdhat = np.squeeze(np.sqrt(((1 - ar1) * 2).tolist())) * func_sd[mask > 0].reshape(-1) + diff_sdhat = np.squeeze(np.sqrt(((1 - ar1) * 2).tolist())) * func_sd diff_sd_mean = diff_sdhat.mean() # Compute temporal difference time series func_diff = np.diff(mfunc, axis=1) # DVARS (no standardization) - dvars_nstd = func_diff.std(axis=0) + dvars_nstd = np.sqrt(np.square(func_diff).mean(axis=0)) # standardization dvars_stdz = dvars_nstd / diff_sd_mean - with warnings.catch_warnings(): # catch, e.g., divide by zero errors + with warnings.catch_warnings(): # catch, e.g., divide by zero errors warnings.filterwarnings('error') # voxelwise standardization - diff_vx_stdz = func_diff / np.array([diff_sdhat] * func_diff.shape[-1]).T - dvars_vx_stdz = diff_vx_stdz.std(axis=0, ddof=1) + diff_vx_stdz = np.square( + func_diff / np.array([diff_sdhat] * func_diff.shape[-1]).T) + dvars_vx_stdz = np.sqrt(diff_vx_stdz.mean(axis=0)) return (dvars_stdz, dvars_nstd, dvars_vx_stdz) -def zero_remove(data, mask): - """ - Modify inputted mask to also mask out zero values - - :param numpy.ndarray data: e.g. voxelwise stddev of fMRI dataset, after motion correction - :param numpy.ndarray mask: brain mask (same dimensions as data) - :return: the mask with any additional zero voxels removed (same dimensions as inputs) - :rtype: numpy.ndarray - - """ - new_mask = mask.copy() - new_mask[data == 0] = 0 - return new_mask def plot_confound(tseries, figsize, name, units=None, series_tr=None, normalize=False): diff --git a/nipype/algorithms/tests/test_confounds.py b/nipype/algorithms/tests/test_confounds.py index a1d6ec9557..7b1fe412d9 100644 --- a/nipype/algorithms/tests/test_confounds.py +++ b/nipype/algorithms/tests/test_confounds.py @@ -39,9 +39,27 @@ def test_dvars(tmpdir): ground_truth = np.loadtxt(example_data('ds003_sub-01_mc.DVARS')) dvars = ComputeDVARS(in_file=example_data('ds003_sub-01_mc.nii.gz'), in_mask=example_data('ds003_sub-01_mc_brainmask.nii.gz'), - save_all=True) + save_all=True, + intensity_normalization=0) os.chdir(str(tmpdir)) res = dvars.run() - dv1 = np.loadtxt(res.outputs.out_std) - assert (np.abs(dv1 - ground_truth).sum()/ len(dv1)) < 0.05 + dv1 = np.loadtxt(res.outputs.out_all, skiprows=1) + assert (np.abs(dv1[:, 0] - ground_truth[:, 0]).sum()/ len(dv1)) < 0.05 + + assert (np.abs(dv1[:, 1] - ground_truth[:, 1]).sum() / len(dv1)) < 0.05 + + assert (np.abs(dv1[:, 2] - ground_truth[:, 2]).sum() / len(dv1)) < 0.05 + + dvars = ComputeDVARS(in_file=example_data('ds003_sub-01_mc.nii.gz'), + in_mask=example_data( + 'ds003_sub-01_mc_brainmask.nii.gz'), + save_all=True) + res = dvars.run() + + dv1 = np.loadtxt(res.outputs.out_all, skiprows=1) + assert (np.abs(dv1[:, 0] - ground_truth[:, 0]).sum() / len(dv1)) < 0.05 + + assert (np.abs(dv1[:, 1] - ground_truth[:, 1]).sum() / len(dv1)) > 0.05 + + assert (np.abs(dv1[:, 2] - ground_truth[:, 2]).sum() / len(dv1)) < 0.05 \ No newline at end of file diff --git a/nipype/testing/data/ds003_sub-01_mc.DVARS b/nipype/testing/data/ds003_sub-01_mc.DVARS index 9622f82fff..bad890e227 100644 --- a/nipype/testing/data/ds003_sub-01_mc.DVARS +++ b/nipype/testing/data/ds003_sub-01_mc.DVARS @@ -1,19 +1,19 @@ -1.54062 -1.31972 -0.921541 -1.26107 -0.99986 -0.929237 -0.715096 -1.05153 -1.29109 -0.700641 -0.844657 -0.884972 -0.807096 -0.881976 -0.843652 -0.780457 -1.05401 -1.32161 -0.686738 +2.02915 5.2016 1.74221 +1.54871 3.97002 1.18108 +0.921419 2.362 0.784497 +1.26058 3.23142 0.734119 +1.00079 2.56548 0.787452 +0.929074 2.38163 0.828835 +0.741207 1.90004 0.746263 +1.07913 2.7663 0.779829 +1.2969 3.32452 0.73856 +0.767387 1.96715 0.772047 +0.847059 2.17138 0.774103 +0.984061 2.52258 0.88097 +0.852897 2.18635 0.794655 +0.927778 2.3783 0.756786 +0.857544 2.19826 0.796125 +0.780098 1.99973 0.731265 +1.05496 2.70434 0.788584 +1.32099 3.38628 0.831803 +0.691529 1.77269 0.738788