From 8b6b2636cac650d858cc2fc45a1d33295d9a213e Mon Sep 17 00:00:00 2001 From: Chris Gorgolewski Date: Mon, 20 Feb 2017 15:50:56 -0800 Subject: [PATCH 1/9] Fix DVARS calculation. --- nipype/algorithms/confounds.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nipype/algorithms/confounds.py b/nipype/algorithms/confounds.py index c6503c703a..e396b8da2a 100644 --- a/nipype/algorithms/confounds.py +++ b/nipype/algorithms/confounds.py @@ -662,7 +662,7 @@ def compute_dvars(in_file, in_mask, remove_zerovariance=False): 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 From 0ce153d3bb5f0ed4433e0f93197e0421b4bde819 Mon Sep 17 00:00:00 2001 From: Chris Gorgolewski Date: Mon, 20 Feb 2017 17:26:53 -0800 Subject: [PATCH 2/9] added intensity normalization --- nipype/algorithms/confounds.py | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/nipype/algorithms/confounds.py b/nipype/algorithms/confounds.py index e396b8da2a..333949a641 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)` @@ -651,6 +663,9 @@ def compute_dvars(in_file, in_mask, remove_zerovariance=False): # Demean mfunc = regress_poly(0, mfunc, remove_mean=True).astype(np.float32) + if intensity_normalization != False: + mfunc = (mfunc / np.median(mfunc)) * intensity_normalization + # Compute (non-robust) estimate of lag-1 autocorrelation ar1 = np.apply_along_axis(AR_est_YW, 1, mfunc, 1)[:, 0] From 36caef9e4fdb76e040375eec3748d28a705fe779 Mon Sep 17 00:00:00 2001 From: Chris Gorgolewski Date: Mon, 20 Feb 2017 17:39:03 -0800 Subject: [PATCH 3/9] cleanup --- nipype/algorithms/confounds.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nipype/algorithms/confounds.py b/nipype/algorithms/confounds.py index 333949a641..621f18eb28 100644 --- a/nipype/algorithms/confounds.py +++ b/nipype/algorithms/confounds.py @@ -663,7 +663,7 @@ def compute_dvars(in_file, in_mask, remove_zerovariance=False, # Demean mfunc = regress_poly(0, mfunc, remove_mean=True).astype(np.float32) - if intensity_normalization != False: + if intensity_normalization != 0: mfunc = (mfunc / np.median(mfunc)) * intensity_normalization # Compute (non-robust) estimate of lag-1 autocorrelation From bab66ef48ef1a08408dd80ead792214fcc94a8cb Mon Sep 17 00:00:00 2001 From: Chris Gorgolewski Date: Mon, 20 Feb 2017 17:39:19 -0800 Subject: [PATCH 4/9] updated test data --- nipype/testing/data/ds003_sub-01_mc.DVARS | 38 +++++++++++------------ 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/nipype/testing/data/ds003_sub-01_mc.DVARS b/nipype/testing/data/ds003_sub-01_mc.DVARS index 9622f82fff..a8e3f6b11d 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 +1.54871 +0.921419 +1.26058 +1.00079 +0.929074 +0.741207 +1.07913 +1.2969 +0.767387 +0.847059 +0.984061 +0.852897 +0.927778 +0.857544 +0.780098 +1.05496 +1.32099 +0.691529 From 7295e6717c7ce2372792131dc0137f089da0f627 Mon Sep 17 00:00:00 2001 From: Chris Gorgolewski Date: Mon, 20 Feb 2017 17:45:17 -0800 Subject: [PATCH 5/9] changelog --- CHANGES | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGES b/CHANGES index 995e324966..af8f71e455 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) * FIX: Issues in Docker image permissions, and docker documentation (https://github.com/nipy/nipype/pull/1825) * ENH: Revised all Dockerfiles and automated deployment to Docker Hub from CircleCI (https://github.com/nipy/nipype/pull/1815) From 8c0ab545441bf8b48213ee55f7123cc10ac7f8cf Mon Sep 17 00:00:00 2001 From: Chris Gorgolewski Date: Mon, 20 Feb 2017 19:03:06 -0800 Subject: [PATCH 6/9] refactor --- nipype/algorithms/confounds.py | 32 +++++++++----------------------- 1 file changed, 9 insertions(+), 23 deletions(-) diff --git a/nipype/algorithms/confounds.py b/nipype/algorithms/confounds.py index 621f18eb28..ba73b2c624 100644 --- a/nipype/algorithms/confounds.py +++ b/nipype/algorithms/confounds.py @@ -648,15 +648,6 @@ 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], :] @@ -666,11 +657,19 @@ def compute_dvars(in_file, in_mask, remove_zerovariance=False, if intensity_normalization != 0: mfunc = (mfunc / np.median(mfunc)) * intensity_normalization + # Robust standard deviation + func_sd = (np.percentile(mfunc, 75, axis=1) - + np.percentile(mfunc, 25, axis=1)) / 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] # 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 @@ -691,19 +690,6 @@ def compute_dvars(in_file, in_mask, remove_zerovariance=False, 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): From 078037f7f97ed543567ea31688b341ef957183f2 Mon Sep 17 00:00:00 2001 From: Chris Gorgolewski Date: Mon, 20 Feb 2017 19:40:08 -0800 Subject: [PATCH 7/9] try without intensity normalization --- nipype/algorithms/tests/test_confounds.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/nipype/algorithms/tests/test_confounds.py b/nipype/algorithms/tests/test_confounds.py index a1d6ec9557..778c4c88c3 100644 --- a/nipype/algorithms/tests/test_confounds.py +++ b/nipype/algorithms/tests/test_confounds.py @@ -39,7 +39,8 @@ 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=False) os.chdir(str(tmpdir)) res = dvars.run() From ebdbf21ab8d7d7b22c753837bfada1bceb86b153 Mon Sep 17 00:00:00 2001 From: Chris Gorgolewski Date: Mon, 20 Feb 2017 19:52:26 -0800 Subject: [PATCH 8/9] fixed voxel wise standardization (see https://github.com/nicholst/DVARS/compare/e0db49d15c2b53359a5060404b82ff6e037ca0ea...master#diff-af96b24619b538d0e8b556069fd0883aR150) --- nipype/algorithms/confounds.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/nipype/algorithms/confounds.py b/nipype/algorithms/confounds.py index ba73b2c624..84837e1c27 100644 --- a/nipype/algorithms/confounds.py +++ b/nipype/algorithms/confounds.py @@ -685,8 +685,8 @@ def compute_dvars(in_file, in_mask, remove_zerovariance=False, 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) From 948ac54d7ffdf57bb046a76adbb4fb2092fcba9f Mon Sep 17 00:00:00 2001 From: Chris Gorgolewski Date: Mon, 20 Feb 2017 22:03:40 -0800 Subject: [PATCH 9/9] Fixed tests. --- nipype/algorithms/confounds.py | 19 ++++++------ nipype/algorithms/tests/test_confounds.py | 23 ++++++++++++-- nipype/testing/data/ds003_sub-01_mc.DVARS | 38 +++++++++++------------ 3 files changed, 49 insertions(+), 31 deletions(-) diff --git a/nipype/algorithms/confounds.py b/nipype/algorithms/confounds.py index 84837e1c27..db616f1f04 100644 --- a/nipype/algorithms/confounds.py +++ b/nipype/algorithms/confounds.py @@ -651,22 +651,22 @@ def compute_dvars(in_file, in_mask, remove_zerovariance=False, 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 - func_sd = (np.percentile(mfunc, 75, axis=1) - - np.percentile(mfunc, 25, axis=1)) / 1.349 + # 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 @@ -681,11 +681,12 @@ def compute_dvars(in_file, in_mask, remove_zerovariance=False, # 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 = np.square(func_diff) / np.array([diff_sdhat] * func_diff.shape[-1]).T + 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) diff --git a/nipype/algorithms/tests/test_confounds.py b/nipype/algorithms/tests/test_confounds.py index 778c4c88c3..7b1fe412d9 100644 --- a/nipype/algorithms/tests/test_confounds.py +++ b/nipype/algorithms/tests/test_confounds.py @@ -40,9 +40,26 @@ def test_dvars(tmpdir): 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, - intensity_normalization=False) + 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 a8e3f6b11d..bad890e227 100644 --- a/nipype/testing/data/ds003_sub-01_mc.DVARS +++ b/nipype/testing/data/ds003_sub-01_mc.DVARS @@ -1,19 +1,19 @@ -2.02915 -1.54871 -0.921419 -1.26058 -1.00079 -0.929074 -0.741207 -1.07913 -1.2969 -0.767387 -0.847059 -0.984061 -0.852897 -0.927778 -0.857544 -0.780098 -1.05496 -1.32099 -0.691529 +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