Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

ENH Expand spm interface thresholding options #3333

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open

ENH Expand spm interface thresholding options #3333

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
d78ba24
Allow for choosing between FDR, FWE and uncorrected threshold
May 3, 2021
662ec09
Add FWEc clustersize correction (use_topo_fwe) and turn off default u…
May 3, 2021
79285f7
Added my name in the .zenodo.json file
May 3, 2021
11f13ae
Merge branch 'master' into expand_spm_functionalities
ldeangelisphys May 4, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions .zenodo.json
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -810,6 +810,11 @@
"affiliation": "MIT, HMS",
"name": "Ghosh, Satrajit",
"orcid": "0000-0002-5312-6729"
},
{
"affiliation": "Netherlands Institute for Neuroscience",
"name": "De Angelis, Lorenzo",
"orcid": "0000-0002-1798-9383"
}
],
"keywords": [
Expand Down
81 changes: 50 additions & 31 deletions nipype/interfaces/spm/model.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -574,23 +574,26 @@ class ThresholdInputSpec(SPMCommandInputSpec):
)
stat_image = File(exists=True, desc="stat image", copyfile=False, mandatory=True)
contrast_index = traits.Int(
mandatory=True, desc="which contrast in the SPM.mat to use"
)
use_fwe_correction = traits.Bool(
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We have to keep the old options around so that we don't break existing workflows. You can deprecate it though (see https://nipype.readthedocs.io/en/latest/devel/interface_specs.html).

True,
mandatory=True, desc='which contrast in the SPM.mat to use')
multitest_correction = traits.Enum(
'none',
'FWE',
'FDR',
usedefault=True,
desc=(
"whether to use FWE (Bonferroni) "
"correction for initial threshold "
"(height_threshold_type has to be "
"set to p-value)"
),
)
desc=('Whether to use a correction '
'for multiple test. '
'Possible choices are FWE, FDR '
'or none'))
use_topo_fdr = traits.Bool(
True,
usedefault=True,
desc=("whether to use FDR over cluster extent probabilities"),
)
False,
usedefault=False,
desc=('whether to use FDR over cluster extent '
'probabilities'))
use_topo_fwe = traits.Bool(
False,
usedefault=False,
desc=('whether to use FWE over cluster extent '
'probabilities'))
height_threshold = traits.Float(
0.05,
usedefault=True,
Expand Down Expand Up @@ -629,6 +632,8 @@ class ThresholdOutputSpec(TraitedSpec):
pre_topo_n_clusters = traits.Int()
activation_forced = traits.Bool()
cluster_forming_thr = traits.Float()
FWEc = traits.Float()
FDRc = traits.Float()


class Threshold(SPMCommand):
Expand All @@ -638,7 +643,6 @@ class Threshold(SPMCommand):

Examples
--------

>>> thresh = Threshold()
>>> thresh.inputs.spm_mat_file = 'SPM.mat'
>>> thresh.inputs.stat_image = 'spmT_0001.img'
Expand All @@ -661,15 +665,19 @@ def _gen_pre_topo_map_filename(self):
def _make_matlab_command(self, _):
script = "con_index = %d;\n" % self.inputs.contrast_index
script += "cluster_forming_thr = %f;\n" % self.inputs.height_threshold
if self.inputs.use_fwe_correction:
script += "thresDesc = 'FWE';\n"
else:
script += "thresDesc = 'none';\n"

script += "thresDesc = '%s';\n" % self.inputs.multitest_correction
# if self.inputs.use_fwe_correction:
# script += "thresDesc = 'FWE';\n"
# else:
# script += "thresDesc = 'none';\n"
if self.inputs.use_topo_fdr:
script += "use_topo_fdr = 1;\n"
else:
script += "use_topo_fdr = 0;\n"
if self.inputs.use_topo_fwe:
script += "use_topo_fwe = 1;\n"
else:
script += "use_topo_fwe = 0;\n"

if self.inputs.force_activation:
script += "force_activation = 1;\n"
Expand All @@ -690,11 +698,15 @@ def _make_matlab_command(self, _):
R = SPM.xVol.R;
S = SPM.xVol.S;
n = 1;
VspmSv = cat(1,SPM.xCon(con_index).Vspm);

switch thresDesc
case 'FWE'
cluster_forming_thr = spm_uc(cluster_forming_thr,df,STAT,R,n,S);

case 'FDR'
cluster_forming_thr = spm_uc_FDR(cluster_forming_thr,df,STAT,n,VspmSv,0);

case 'none'
if strcmp(height_threshold_type, 'p-value')
cluster_forming_thr = spm_u(cluster_forming_thr^(1/n),df,STAT);
Expand All @@ -721,14 +733,15 @@ def _make_matlab_command(self, _):
max_size_index = 0;
th_nclusters = 0;
nclusters = 0;
uc = 0;
ue = 0;
if isempty(XYZth)
thresholded_XYZ = [];
thresholded_Z = [];
else
if use_topo_fdr
V2R = 1/prod(FWHM(stat_map_vol.dim > 1));
[uc,Pc,ue] = spm_uc_clusterFDR(cluster_extent_p_fdr_thr,df,STAT,R,n,Z,XYZ,V2R,cluster_forming_thr);
end

V2R = 1/prod(FWHM(stat_map_vol.dim > 1));
[uc,Pc,ue] = spm_uc_clusterFDR(cluster_extent_p_fdr_thr,df,STAT,R,n,Z,XYZ,V2R,cluster_forming_thr);

voxel_labels = spm_clusters(XYZth);
nclusters = max(voxel_labels);
Expand All @@ -738,7 +751,7 @@ def _make_matlab_command(self, _):

for i = 1:nclusters
cluster_size = sum(voxel_labels==i);
if cluster_size > extent_threshold && (~use_topo_fdr || (cluster_size - uc) > -1)
if cluster_size > extent_threshold && (~use_topo_fdr || (cluster_size - uc) > -1) && (~use_topo_fwe || (cluster_size - ue) > -1)
thresholded_XYZ = cat(2, thresholded_XYZ, XYZth(:,voxel_labels == i));
thresholded_Z = cat(2, thresholded_Z, Zth(voxel_labels == i));
th_nclusters = th_nclusters + 1;
Expand Down Expand Up @@ -771,6 +784,8 @@ def _make_matlab_command(self, _):
fprintf('pre_topo_n_clusters = %d\\n',nclusters);
fprintf('n_clusters = %d\\n',th_nclusters);
fprintf('cluster_forming_thr = %f\\n',cluster_forming_thr);
fprintf('FWEc = %d\\n',ue);
fprintf('FDRc = %d\\n',uc);

"""
script += (
Expand Down Expand Up @@ -803,11 +818,15 @@ def aggregate_outputs(self, runtime=None):
int(line[len("pre_topo_n_clusters = ") :].strip()),
)
elif line.startswith("cluster_forming_thr = "):
setattr(
outputs,
"cluster_forming_thr",
float(line[len("cluster_forming_thr = ") :].strip()),
)
setattr(outputs, 'cluster_forming_thr',
float(line[len("cluster_forming_thr = "):].strip()))
elif line.startswith("FWEc = "):
setattr(outputs, 'FWEc',
float(line[len("FWEc = "):].strip()))
elif line.startswith("FDRc = "):
setattr(outputs, 'FDRc',
float(line[len("FDRc = "):].strip()))

return outputs

def _list_outputs(self):
Expand Down