From 1f5687602a3a06834e9cfd5213b195a380fe860b Mon Sep 17 00:00:00 2001 From: cmaumet Date: Tue, 14 Mar 2017 13:50:55 +0000 Subject: [PATCH 01/29] Create json structure for nidm export --- snpm_cp.m | 23 ++++++++++++++++++++++ snpm_pp.m | 58 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ snpm_ui.m | 12 ++++++++++-- 3 files changed, 91 insertions(+), 2 deletions(-) diff --git a/snpm_cp.m b/snpm_cp.m index e322346..533aeb2 100644 --- a/snpm_cp.m +++ b/snpm_cp.m @@ -223,15 +223,27 @@ function snpm_cp(CWD) error('SnPM:InvalidContrast','Contrast problem; wrong number of columns'); end if size(CONT,1) > 1 + nidm_json.('Contrasts') = struct(); + contrast = struct(); + contrast.('obo_contrastweightmatrix/prov:value') = CONT; warning('SnPM:FContrast', ... 'F contrast! F statistic images are being created.'); STAT = 'F'; + contrast.('nidm_StatisticMap/nidm_statisticType') = 'obo_Fstatistic'; if (CONT(1,:) == -CONT(2,:)) CONT = CONT(1,:); end + con_name = [mat2str(CONT) ' (F)']; + nidm_json.('Contrasts').(con_name) = contrast; else + con_name = [mat2str(CONT) ' (T)']; + con_neg_name = [mat2str(-CONT) ' (T)']; STAT = 'T'; + contrast.('nidm_StatisticMap/nidm_statisticType') = 'obo_tstatistic'; + nidm_json.('Contrasts').(con_name) = contrast; + nidm_json.('Contrasts').(con_neg_name) = contrast; end + if rank(CONT) u -tol) / nPermReal; + nidm_inference.('nidm_HeightThreshold/prov:type') = 'obo_statistic'; + nidm_inference.('nidm_HeightThreshold/prov:value') = u; case 'FDRth' alph_FDR = BoundCheck(job.Thr.Vox.VoxSig.FDRth,[0 1],'Invalid FDR level'); alpha_ucp = snpm_uc_FDR(alph_FDR,[],'P',[],sSnPMucp'); + nidm_inference.('nidm_HeightThreshold/prov:type') = 'obo_qvalue'; + nidm_inference.('nidm_HeightThreshold/prov:value') = alph_FDR; case 'FWEth' alph_FWE = BoundCheck(job.Thr.Vox.VoxSig.FWEth,[0 1],'Invalid FWE level'); iFWE = ceil((1-alph_FWE)*nPermReal); @@ -508,7 +523,12 @@ function snpm_pp(CWD,varargin) C_MaxT = 0; end u = C_MaxT; + nidm_inference.('nidm_HeightThreshold/prov:type') = 'obo_FWERadjustedpvalue'; + nidm_inference.('nidm_HeightThreshold/prov:value') = alph_FWE; end + % No extent thresholding when voxelwise threshold is requested + nidm_inference.('nidm_ExtentThreshold/prov:type') = 'obo_statistic'; + nidm_inference.('nidm_ExtentThreshold/prov:value') = 0; else % Cluster-wise inference if exist(fullfile(CWD,'SnPM_ST.mat'))~=2 & exist(fullfile(CWD,'STCS.mat'))~=2 @@ -523,6 +543,14 @@ function snpm_pp(CWD,varargin) % Save original ST_Ut ST_Ut_0 = ST_Ut; CFth=job.Thr.Clus.ClusSize.CFth; + if (CFth<1) + nidm_inference.('nidm_ExtentThreshold/prov:type') = 'nidm_PValueUncorrected'; + nidm_inference.('nidm_ExtentThreshold/prov:value') = CFth; + else + nidm_inference.('nidm_ExtentThreshold/prov:type') = 'obo_statistic'; + nidm_inference.('nidm_ExtentThreshold/prov:value') = CFth; + end + if (CFth<=0) error('SnPM:InvalidClusterFormingThresh', 'ERROR: Cluster-forming threshold must be strictly positive.\nRe-run results with a cluster-forming threshold greater than 0.\n') end @@ -585,11 +613,17 @@ function snpm_pp(CWD,varargin) switch tmp{1} case 'Cth' C_STCS = job.Thr.Clus.ClusSize.ClusSig.Cth; + nidm_inference.('nidm_ExtentThreshold/prov:type') = 'obo_statistic'; + nidm_inference.('nidm_ExtentThreshold/prov:value') = C_STCS; case 'PthC' alpha_ucp = BoundCheck(job.Thr.Clus.ClusSize.ClusSig.PthC,[0 1],'Invalid uncorrected P(k)'); + nidm_inference.('nidm_ExtentThreshold/prov:type') = 'nidm_PValueUncorrected'; + nidm_inference.('nidm_ExtentThreshold/prov:value') = alpha_ucp; case 'FWEthC' alph_FWE = BoundCheck(job.Thr.Clus.ClusSize.ClusSig.FWEthC,[0 1],'Invalid FWE level (cluster-level inference)'); iFWE = ceil((1-alph_FWE)*nPermReal); + nidm_inference.('nidm_ExtentThreshold/prov:type') = 'obo_FWERadjustedpvalue'; + nidm_inference.('nidm_ExtentThreshold/prov:value') = alph_FWE; end end % END: Cluster-wise inference @@ -896,6 +930,8 @@ function snpm_pp(CWD,varargin) if i==1 %-Save perm 1 stats for use later - [X;Y;Z;T;perm;STCno] tmp = spm_clusters(Locs_vox(1:3,:)); + + nidm_json.('nidm_ClusterDefinitionCriteria/nidm_hasConnectivityCriterion') = 'nidm_voxel18connected'; if isPos==1 STCstats_Pos = [ SnPM_ST(:,tQ); tmp]; if bNeg==0 @@ -1311,7 +1347,11 @@ function snpm_pp(CWD,varargin) %----------------------------------------------------------------------- r = 1; bUsed = zeros(size(STC_SnPMt)); + nidm_clusters = struct(); + while max(STC_SnPMt.*(~bUsed)) & (y > 3) + nidm_cluster = struct(); + nidm_peaks = struct(); [null, i] = max(STC_SnPMt.*(~bUsed)); % Largest t value j = find(STC_r == STC_r(i)); % Maxima in same region @@ -1341,10 +1381,19 @@ function snpm_pp(CWD,varargin) text(tCol(8),y,sprintf(Fmtst{8},STC_XYZ(1,i)),'UserData',STC_XYZ(:,i),StrAttrB{:}) text(tCol(9),y,sprintf(Fmtst{9},STC_XYZ(2,i)),'UserData',STC_XYZ(:,i),StrAttrB{:}) text(tCol(10),y,sprintf(Fmtst{10},STC_XYZ(3,i)),'UserData',STC_XYZ(:,i),StrAttrB{:}) + + nidm_cluster.('nidm_SupraThresholdCluster/nidm_clusterSizeInVoxels') = STC_N(i); + nidm_peaks.('1').('nidm_Peak/nidm_pValueFWER') = Pt(i); + nidm_peaks.('1').('nidm_Peak/nidm_qValueFDR') = Pfdr(i); + nidm_peaks.('1').('nidm_Peak/prov:value') = STC_SnPMt(i); + nidm_peaks.('1').('nidm_Peak/nidm_pValueUncorrected') = Pu(i); + nidm_peaks.('1').('nidm_Coordinate/nidm_coordinateVector') = STC_XYZ(1:3,i); + y = y -1; %-Print up to 3 secondary maxima (>8mm apart) %------------------------------------------------------------------- + nidm_inference.('nidm_PeakDefinitionCriteria/nidm_minDistanceBetweenPeaks') = 8; [null, k] = sort(-STC_SnPMt(j)); % Sort on t value D = i; for i = 1:length(k) @@ -1365,11 +1414,20 @@ function snpm_pp(CWD,varargin) y = y -1; end end + nidm_peaks.(num2str(i)).('nidm_Peak/nidm_pValueFWER') = Pt(d); + nidm_peaks.(num2str(i)).('nidm_Peak/nidm_qValueFDR') = Pfdr(d); + nidm_peaks.(num2str(i)).('nidm_Peak/prov:value') = STC_SnPMt(d); + nidm_peaks.(num2str(i)).('nidm_Peak/nidm_pValueUncorrected') = Pu(d); + nidm_peaks.(num2str(i)).('nidm_Coordinate/nidm_coordinateVector') = STC_XYZ(1:3,d); end bUsed(j) = (bUsed(j) | 1 ); %-Mark maxima as "used" r = r + 1; % Next region + + nidm_cluster.('Peaks') = nidm_peaks; + nidm_clusters(num2str(r)) = nidm_cluster; end + nidm_inference.('Clusters') = nidm_clusters; clear i j k D d r diff --git a/snpm_ui.m b/snpm_ui.m index 7f146bf..a7c76e0 100644 --- a/snpm_ui.m +++ b/snpm_ui.m @@ -123,6 +123,8 @@ function snpm_ui(varargin) % V Memory mapping handles % MASK Filename of explicit mask image % ImMASK Implicit masking; 0=none; 1=zeros are equivalent to NaN +% nidm_json json structure storing minimal information required for a +% NIDM-Results export % % df degrees of freedom due to error % sDesSave String of PlugIn variables to save to cfg file @@ -190,7 +192,9 @@ function snpm_ui(varargin) end cd(job.dir{1}) end - + +nidm_json = struct(); + %-Definitions & Design parameters %======================================================================= sDesigns=str2mat(... @@ -368,6 +372,7 @@ function snpm_ui(varargin) end if (iGMsca==2) % CHANGED from 1 to 2 as should not ask for a value if grand mean scaling is not required. + nidm_json.('nidm_Data/nidm_grandMeanScaling') = true; if (iGloNorm==2) % Proportional scaling str = 'PropSca global mean to'; else @@ -381,8 +386,10 @@ function snpm_ui(varargin) case 'gmsca_no', GM = 50; end + nidm_json.('nidm_Data/nidm_targetIntensity') = GM; elseif (iGMsca==1) % No grand mean scaling GM = 0; + nidm_json.('nidm_Data/nidm_grandMeanScaling') = false; end @@ -572,13 +579,14 @@ function snpm_ui(varargin) %-Construct full design matrix and name matrices for display %----------------------------------------------------------------------- [nHCBG,HCBGnames] = spm_DesMtx('Sca',H,Hnames,C,Cnames,B,Bnames,G,Gnames); +nidm_json.('nidm_DesignMatrix/prov:value') = nHCBG; %-Setup is complete - save SnPMcfg Mat file %----------------------------------------------------------------------- s_SnPMcfg_save = ['s_SnPMcfg_save H C B G HCBGnames P PiCond ',... 'sPiCond bhPerms sHCform iGloNorm sGloNorm GM rg GX GMscale CONT ',... 'THRESH MASK ImMASK TH bVarSm vFWHM sVarSm bVolm bST sDesFile sDesign ',... - 'V pU_ST_Ut df1 ', ... + 'V pU_ST_Ut df1 nidm_json', ... 'sDesSave ',sDesSave]; eval(['save SnPMcfg ',s_SnPMcfg_save]) From 4530ddb9288385180fbfc3c3db253ff7ea5ef098 Mon Sep 17 00:00:00 2001 From: cmaumet Date: Tue, 14 Mar 2017 14:10:44 +0000 Subject: [PATCH 02/29] Save and load NIDM json from file --- snpm_cp.m | 4 ++++ snpm_pp.m | 6 ++++-- snpm_ui.m | 2 +- 3 files changed, 9 insertions(+), 3 deletions(-) diff --git a/snpm_cp.m b/snpm_cp.m index 533aeb2..ed31ff7 100644 --- a/snpm_cp.m +++ b/snpm_cp.m @@ -209,10 +209,13 @@ function snpm_cp(CWD) fprintf('Changing directory to %s\n',CWD); end CfgFile = fullfile(CWD,'SnPMcfg.mat'); +NIDMFile = fullfile(CWD,'snpm_nidm.json'); %-Load config file & catch all problem cases now %----------------------------------------------------------------------- load(CfgFile); +nidm_json = spm_jsonread(NIDMFile); + if isempty([H C]) error('SnPM:NoModel', 'No model specified; [H C] empty'); end @@ -1194,6 +1197,7 @@ function snpm_cp(CWD) %-Save key variables %======================================================================= eval(['save SnPM ',s_SnPM_save]) +spm_jsonwrite(NIDMFile, nidm_json) %-Print quick summary info (allowing for STOPing) %======================================================================= diff --git a/snpm_pp.m b/snpm_pp.m index 037666d..6e98277 100644 --- a/snpm_pp.m +++ b/snpm_pp.m @@ -369,7 +369,7 @@ function snpm_pp(CWD,varargin) load(fullfile(CWD,'SnPMcfg')) load(fullfile(CWD,'SnPM')) load(fullfile(CWD,'SnPMucp')) - +nidm_json = spm_jsonread(fullfile(CWD,'snpm_nidm.json')); %-Ask whether positive or negative effects be analysed %----------------------------------------------------------------------- @@ -1430,7 +1430,8 @@ function snpm_pp(CWD,varargin) nidm_inference.('Clusters') = nidm_clusters; clear i j k D d r - + nidm_json.('Inferences') = nidm_inference; + spm_jsonwrite('snpm_nidm.json', nidm_json) %-Footnote with SnPM parameters %======================================================================= line([0,1],[0.5,0.5],'LineWidth',1,'Color','r') @@ -1689,6 +1690,7 @@ function ShowDist(T,cT,aT,C,cC,aC,Typ) end + return diff --git a/snpm_ui.m b/snpm_ui.m index a7c76e0..06bc80b 100644 --- a/snpm_ui.m +++ b/snpm_ui.m @@ -589,7 +589,7 @@ function snpm_ui(varargin) 'V pU_ST_Ut df1 nidm_json', ... 'sDesSave ',sDesSave]; eval(['save SnPMcfg ',s_SnPMcfg_save]) - +spm_jsonwrite('snpm_nidm.json',nidm_json) if ~spm_get_defaults('cmdline') From e342f9420b812a96d22f443e14b0d89b996c7040 Mon Sep 17 00:00:00 2001 From: cmaumet Date: Tue, 14 Mar 2017 14:22:28 +0000 Subject: [PATCH 03/29] Remove special characters from field names MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit To avoid « Invalid field name » error in Matlab --- snpm_cp.m | 16 +++++----- snpm_pp.m | 68 ++++++++++++++++++++--------------------- snpm_ui.m | 8 ++--- test/snpm_test_config.m | 2 +- 4 files changed, 47 insertions(+), 47 deletions(-) diff --git a/snpm_cp.m b/snpm_cp.m index ed31ff7..9c72594 100644 --- a/snpm_cp.m +++ b/snpm_cp.m @@ -228,11 +228,11 @@ function snpm_cp(CWD) if size(CONT,1) > 1 nidm_json.('Contrasts') = struct(); contrast = struct(); - contrast.('obo_contrastweightmatrix/prov:value') = CONT; + contrast.('obo_contrastweightmatrix__prov_value') = CONT; warning('SnPM:FContrast', ... 'F contrast! F statistic images are being created.'); STAT = 'F'; - contrast.('nidm_StatisticMap/nidm_statisticType') = 'obo_Fstatistic'; + contrast.('nidm_StatisticMap__nidm_statisticType') = 'obo_Fstatistic'; if (CONT(1,:) == -CONT(2,:)) CONT = CONT(1,:); end @@ -242,7 +242,7 @@ function snpm_cp(CWD) con_name = [mat2str(CONT) ' (T)']; con_neg_name = [mat2str(-CONT) ' (T)']; STAT = 'T'; - contrast.('nidm_StatisticMap/nidm_statisticType') = 'obo_tstatistic'; + contrast.('nidm_StatisticMap__nidm_statisticType') = 'obo_tstatistic'; nidm_json.('Contrasts').(con_name) = contrast; nidm_json.('Contrasts').(con_neg_name) = contrast; end @@ -307,7 +307,7 @@ function snpm_cp(CWD) df = q - r; %-Residual degrees of freedom nPerm = size(PiCond,1); %-# permutations -nidm_json.('Contrasts').(con_name).('nidm_StatisticMap/nidm_errorDegreesOfFreedom') = df; +nidm_json.('Contrasts').(con_name).('nidm_StatisticMap__nidm_errorDegreesOfFreedom') = df; %-Get ORIGIN, etc %----------------------------------------------------------------------- @@ -390,13 +390,13 @@ function snpm_cp(CWD) Vbeta(ii)=snpm_clone_vol(Vt,fname,descrip); nidm_json.('ParameterEstimateMaps').(descrip).(... - 'nidm_ParameterEstimateMap/prov:atLocation') = fname; + 'nidm_ParameterEstimateMap__prov_atLocation') = fname; end Vbeta = spm_create_vol(Vbeta); VResMS=snpm_clone_vol(Vt,'ResMS.img','Residual sum-of-squares'); VResMS=spm_create_vol(VResMS); -nidm_json.('nidm_ResidualMeanSquaresMap/prov:atLocation') = 'ResMS.img'; +nidm_json.('nidm_ResidualMeanSquaresMap__prov_atLocation') = 'ResMS.img'; if bVarSm==0 str = sprintf('%c_{%d} statistic',STAT,df); @@ -411,10 +411,10 @@ function snpm_cp(CWD) if STAT=='T' VT_pos=snpm_clone_vol(Vt,'snpmT+.img',[str,' (+ve)']); VT_pos=spm_create_vol(VT_pos); - nidm_json.(con_name).('nidm_StatisticMap/prov:atLocation') = 'snpmT+.img'; + nidm_json.(con_name).('nidm_StatisticMap__prov_atLocation') = 'snpmT+.img'; VT_neg=snpm_clone_vol(Vt,'snpmT-.img',[str,' (-ve)']); VT_neg=spm_create_vol(VT_neg); - nidm_json.(con_neg_name).('nidm_StatisticMap/prov:atLocation') = 'snpmT-.img'; + nidm_json.(con_neg_name).('nidm_StatisticMap__prov_atLocation') = 'snpmT-.img'; elseif STAT=='F' VF=snpm_clone_vol(Vt,'snpmF.img',str); VF=spm_create_vol(VF); diff --git a/snpm_pp.m b/snpm_pp.m index 6e98277..9f33a79 100644 --- a/snpm_pp.m +++ b/snpm_pp.m @@ -492,7 +492,7 @@ function snpm_pp(CWD,varargin) alph_FDR = NaN; % FDR rate of a specified alpha_ucp nidm_inference = struct(); % NIDM structre storing information about inference -nidm_inference.('nidm_Inference/nidm_hasAlternativeHypothesis') = 'nidm_OneTailedTest'; +nidm_inference.('nidm_Inference__nidm_hasAlternativeHypothesis') = 'nidm_OneTailedTest'; if BATCH bSpatEx = isfield(job.Thr,'Clus'); if ~bSpatEx @@ -502,18 +502,18 @@ function snpm_pp(CWD,varargin) case 'Pth' alpha_ucp = BoundCheck(job.Thr.Vox.VoxSig.Pth,[0 1],'Invalid Uncorrected P'); alph_FDR = snpm_P_FDR(alpha_ucp,[],'P',[],sSnPMucp'); - nidm_inference.('nidm_HeightThreshold/prov:type') = 'nidm_PValueUncorrected'; - nidm_inference.('nidm_HeightThreshold/prov:value') = alpha_ucp; + nidm_inference.('nidm_HeightThreshold__prov_type') = 'nidm_PValueUncorrected'; + nidm_inference.('nidm_HeightThreshold__prov_value') = alpha_ucp; case 'TFth' u = BoundCheck(job.Thr.Vox.VoxSig.TFth,[0 Inf],'Negative Threshold!'); alph_FWE = sum(MaxT > u -tol) / nPermReal; - nidm_inference.('nidm_HeightThreshold/prov:type') = 'obo_statistic'; - nidm_inference.('nidm_HeightThreshold/prov:value') = u; + nidm_inference.('nidm_HeightThreshold__prov_type') = 'obo_statistic'; + nidm_inference.('nidm_HeightThreshold__prov_value') = u; case 'FDRth' alph_FDR = BoundCheck(job.Thr.Vox.VoxSig.FDRth,[0 1],'Invalid FDR level'); alpha_ucp = snpm_uc_FDR(alph_FDR,[],'P',[],sSnPMucp'); - nidm_inference.('nidm_HeightThreshold/prov:type') = 'obo_qvalue'; - nidm_inference.('nidm_HeightThreshold/prov:value') = alph_FDR; + nidm_inference.('nidm_HeightThreshold__prov_type') = 'obo_qvalue'; + nidm_inference.('nidm_HeightThreshold__prov_value') = alph_FDR; case 'FWEth' alph_FWE = BoundCheck(job.Thr.Vox.VoxSig.FWEth,[0 1],'Invalid FWE level'); iFWE = ceil((1-alph_FWE)*nPermReal); @@ -523,12 +523,12 @@ function snpm_pp(CWD,varargin) C_MaxT = 0; end u = C_MaxT; - nidm_inference.('nidm_HeightThreshold/prov:type') = 'obo_FWERadjustedpvalue'; - nidm_inference.('nidm_HeightThreshold/prov:value') = alph_FWE; + nidm_inference.('nidm_HeightThreshold__prov_type') = 'obo_FWERadjustedpvalue'; + nidm_inference.('nidm_HeightThreshold__prov_value') = alph_FWE; end % No extent thresholding when voxelwise threshold is requested - nidm_inference.('nidm_ExtentThreshold/prov:type') = 'obo_statistic'; - nidm_inference.('nidm_ExtentThreshold/prov:value') = 0; + nidm_inference.('nidm_ExtentThreshold__prov_type') = 'obo_statistic'; + nidm_inference.('nidm_ExtentThreshold__prov_value') = 0; else % Cluster-wise inference if exist(fullfile(CWD,'SnPM_ST.mat'))~=2 & exist(fullfile(CWD,'STCS.mat'))~=2 @@ -544,11 +544,11 @@ function snpm_pp(CWD,varargin) ST_Ut_0 = ST_Ut; CFth=job.Thr.Clus.ClusSize.CFth; if (CFth<1) - nidm_inference.('nidm_ExtentThreshold/prov:type') = 'nidm_PValueUncorrected'; - nidm_inference.('nidm_ExtentThreshold/prov:value') = CFth; + nidm_inference.('nidm_ExtentThreshold__prov_type') = 'nidm_PValueUncorrected'; + nidm_inference.('nidm_ExtentThreshold__prov_value') = CFth; else - nidm_inference.('nidm_ExtentThreshold/prov:type') = 'obo_statistic'; - nidm_inference.('nidm_ExtentThreshold/prov:value') = CFth; + nidm_inference.('nidm_ExtentThreshold__prov_type') = 'obo_statistic'; + nidm_inference.('nidm_ExtentThreshold__prov_value') = CFth; end if (CFth<=0) @@ -613,17 +613,17 @@ function snpm_pp(CWD,varargin) switch tmp{1} case 'Cth' C_STCS = job.Thr.Clus.ClusSize.ClusSig.Cth; - nidm_inference.('nidm_ExtentThreshold/prov:type') = 'obo_statistic'; - nidm_inference.('nidm_ExtentThreshold/prov:value') = C_STCS; + nidm_inference.('nidm_ExtentThreshold__prov_type') = 'obo_statistic'; + nidm_inference.('nidm_ExtentThreshold__prov_value') = C_STCS; case 'PthC' alpha_ucp = BoundCheck(job.Thr.Clus.ClusSize.ClusSig.PthC,[0 1],'Invalid uncorrected P(k)'); - nidm_inference.('nidm_ExtentThreshold/prov:type') = 'nidm_PValueUncorrected'; - nidm_inference.('nidm_ExtentThreshold/prov:value') = alpha_ucp; + nidm_inference.('nidm_ExtentThreshold__prov_type') = 'nidm_PValueUncorrected'; + nidm_inference.('nidm_ExtentThreshold__prov_value') = alpha_ucp; case 'FWEthC' alph_FWE = BoundCheck(job.Thr.Clus.ClusSize.ClusSig.FWEthC,[0 1],'Invalid FWE level (cluster-level inference)'); iFWE = ceil((1-alph_FWE)*nPermReal); - nidm_inference.('nidm_ExtentThreshold/prov:type') = 'obo_FWERadjustedpvalue'; - nidm_inference.('nidm_ExtentThreshold/prov:value') = alph_FWE; + nidm_inference.('nidm_ExtentThreshold__prov_type') = 'obo_FWERadjustedpvalue'; + nidm_inference.('nidm_ExtentThreshold__prov_value') = alph_FWE; end end % END: Cluster-wise inference @@ -931,7 +931,7 @@ function snpm_pp(CWD,varargin) %-Save perm 1 stats for use later - [X;Y;Z;T;perm;STCno] tmp = spm_clusters(Locs_vox(1:3,:)); - nidm_json.('nidm_ClusterDefinitionCriteria/nidm_hasConnectivityCriterion') = 'nidm_voxel18connected'; + nidm_json.('nidm_ClusterDefinitionCriteria__nidm_hasConnectivityCriterion') = 'nidm_voxel18connected'; if isPos==1 STCstats_Pos = [ SnPM_ST(:,tQ); tmp]; if bNeg==0 @@ -1382,18 +1382,18 @@ function snpm_pp(CWD,varargin) text(tCol(9),y,sprintf(Fmtst{9},STC_XYZ(2,i)),'UserData',STC_XYZ(:,i),StrAttrB{:}) text(tCol(10),y,sprintf(Fmtst{10},STC_XYZ(3,i)),'UserData',STC_XYZ(:,i),StrAttrB{:}) - nidm_cluster.('nidm_SupraThresholdCluster/nidm_clusterSizeInVoxels') = STC_N(i); - nidm_peaks.('1').('nidm_Peak/nidm_pValueFWER') = Pt(i); - nidm_peaks.('1').('nidm_Peak/nidm_qValueFDR') = Pfdr(i); - nidm_peaks.('1').('nidm_Peak/prov:value') = STC_SnPMt(i); - nidm_peaks.('1').('nidm_Peak/nidm_pValueUncorrected') = Pu(i); - nidm_peaks.('1').('nidm_Coordinate/nidm_coordinateVector') = STC_XYZ(1:3,i); + nidm_cluster.('nidm_SupraThresholdCluster__nidm_clusterSizeInVoxels') = STC_N(i); + nidm_peaks.('1').('nidm_Peak__nidm_pValueFWER') = Pt(i); + nidm_peaks.('1').('nidm_Peak__nidm_qValueFDR') = Pfdr(i); + nidm_peaks.('1').('nidm_Peak__prov_value') = STC_SnPMt(i); + nidm_peaks.('1').('nidm_Peak__nidm_pValueUncorrected') = Pu(i); + nidm_peaks.('1').('nidm_Coordinate__nidm_coordinateVector') = STC_XYZ(1:3,i); y = y -1; %-Print up to 3 secondary maxima (>8mm apart) %------------------------------------------------------------------- - nidm_inference.('nidm_PeakDefinitionCriteria/nidm_minDistanceBetweenPeaks') = 8; + nidm_inference.('nidm_PeakDefinitionCriteria__nidm_minDistanceBetweenPeaks') = 8; [null, k] = sort(-STC_SnPMt(j)); % Sort on t value D = i; for i = 1:length(k) @@ -1414,11 +1414,11 @@ function snpm_pp(CWD,varargin) y = y -1; end end - nidm_peaks.(num2str(i)).('nidm_Peak/nidm_pValueFWER') = Pt(d); - nidm_peaks.(num2str(i)).('nidm_Peak/nidm_qValueFDR') = Pfdr(d); - nidm_peaks.(num2str(i)).('nidm_Peak/prov:value') = STC_SnPMt(d); - nidm_peaks.(num2str(i)).('nidm_Peak/nidm_pValueUncorrected') = Pu(d); - nidm_peaks.(num2str(i)).('nidm_Coordinate/nidm_coordinateVector') = STC_XYZ(1:3,d); + nidm_peaks.(num2str(i)).('nidm_Peak__nidm_pValueFWER') = Pt(d); + nidm_peaks.(num2str(i)).('nidm_Peak__nidm_qValueFDR') = Pfdr(d); + nidm_peaks.(num2str(i)).('nidm_Peak__prov_value') = STC_SnPMt(d); + nidm_peaks.(num2str(i)).('nidm_Peak__nidm_pValueUncorrected') = Pu(d); + nidm_peaks.(num2str(i)).('nidm_Coordinate__nidm_coordinateVector') = STC_XYZ(1:3,d); end bUsed(j) = (bUsed(j) | 1 ); %-Mark maxima as "used" diff --git a/snpm_ui.m b/snpm_ui.m index 06bc80b..dfbab63 100644 --- a/snpm_ui.m +++ b/snpm_ui.m @@ -372,7 +372,7 @@ function snpm_ui(varargin) end if (iGMsca==2) % CHANGED from 1 to 2 as should not ask for a value if grand mean scaling is not required. - nidm_json.('nidm_Data/nidm_grandMeanScaling') = true; + nidm_json.('nidm_Data__nidm_grandMeanScaling') = true; if (iGloNorm==2) % Proportional scaling str = 'PropSca global mean to'; else @@ -386,10 +386,10 @@ function snpm_ui(varargin) case 'gmsca_no', GM = 50; end - nidm_json.('nidm_Data/nidm_targetIntensity') = GM; + nidm_json.('nidm_Data__nidm_targetIntensity') = GM; elseif (iGMsca==1) % No grand mean scaling GM = 0; - nidm_json.('nidm_Data/nidm_grandMeanScaling') = false; + nidm_json.('nidm_Data__nidm_grandMeanScaling') = false; end @@ -579,7 +579,7 @@ function snpm_ui(varargin) %-Construct full design matrix and name matrices for display %----------------------------------------------------------------------- [nHCBG,HCBGnames] = spm_DesMtx('Sca',H,Hnames,C,Cnames,B,Bnames,G,Gnames); -nidm_json.('nidm_DesignMatrix/prov:value') = nHCBG; +nidm_json.('nidm_DesignMatrix__prov_value') = nHCBG; %-Setup is complete - save SnPMcfg Mat file %----------------------------------------------------------------------- diff --git a/test/snpm_test_config.m b/test/snpm_test_config.m index 0f781d2..97b26be 100644 --- a/test/snpm_test_config.m +++ b/test/snpm_test_config.m @@ -8,6 +8,6 @@ % Camille Maumet global testDataDir; -testDataDir = ''; % Specify directory containing test data +testDataDir = '/Users/cmaumet/Desktop/SnPM_test_data'; % Specify directory containing test data global SnPMrefVersion SnPMrefVersion = 'SnPM8'; \ No newline at end of file From 02d3cd9f363c829a5d7fe2e312dffbece2945225 Mon Sep 17 00:00:00 2001 From: cmaumet Date: Tue, 14 Mar 2017 15:09:23 +0000 Subject: [PATCH 04/29] first run on test one sample --- snpm_cp.m | 27 +++++++++++++++++++-------- snpm_pp.m | 36 ++++++++++++++++++------------------ snpm_ui.m | 2 +- 3 files changed, 38 insertions(+), 27 deletions(-) diff --git a/snpm_cp.m b/snpm_cp.m index 9c72594..bfa04a6 100644 --- a/snpm_cp.m +++ b/snpm_cp.m @@ -236,15 +236,23 @@ function snpm_cp(CWD) if (CONT(1,:) == -CONT(2,:)) CONT = CONT(1,:); end - con_name = [mat2str(CONT) ' (F)']; + con_name = 'Positive'; nidm_json.('Contrasts').(con_name) = contrast; else - con_name = [mat2str(CONT) ' (T)']; - con_neg_name = [mat2str(-CONT) ' (T)']; + con_name = 'Positive'; + con_neg_name = 'Negative'; STAT = 'T'; + contrast.('nidm_StatisticMap__nidm_statisticType') = 'obo_tstatistic'; - nidm_json.('Contrasts').(con_name) = contrast; - nidm_json.('Contrasts').(con_neg_name) = contrast; + contrast_pos = contrast; + contrast_neg = contrast; + contrast_pos.('nidm_ContrastMap__nidm_contrastName') = ['T: ' mat2str(CONT)]; + contrast_pos.('obo_contrastweightmatrix__prov_value') = CONT; + contrast_neg.('nidm_ContrastMap__nidm_contrastName') = ['T: ' mat2str(-CONT)]; + contrast_neg.('obo_contrastweightmatrix__prov_value') = -CONT; + + nidm_json.('Contrasts').(con_name) = contrast_pos; + nidm_json.('Contrasts').(con_neg_name) = contrast_neg; end if rank(CONT)8mm apart) %------------------------------------------------------------------- - nidm_inference.('nidm_PeakDefinitionCriteria__nidm_minDistanceBetweenPeaks') = 8; + nidm_json.('Inferences').('nidm_PeakDefinitionCriteria__nidm_minDistanceBetweenPeaks') = 8; [null, k] = sort(-STC_SnPMt(j)); % Sort on t value D = i; for i = 1:length(k) @@ -1414,24 +1414,24 @@ function snpm_pp(CWD,varargin) y = y -1; end end - nidm_peaks.(num2str(i)).('nidm_Peak__nidm_pValueFWER') = Pt(d); - nidm_peaks.(num2str(i)).('nidm_Peak__nidm_qValueFDR') = Pfdr(d); - nidm_peaks.(num2str(i)).('nidm_Peak__prov_value') = STC_SnPMt(d); - nidm_peaks.(num2str(i)).('nidm_Peak__nidm_pValueUncorrected') = Pu(d); - nidm_peaks.(num2str(i)).('nidm_Coordinate__nidm_coordinateVector') = STC_XYZ(1:3,d); + nidm_peaks.(['Peak_' num2str(i)]).('nidm_Peak__nidm_pValueFWER') = Pt(d); + nidm_peaks.(['Peak_' num2str(i)]).('nidm_Peak__nidm_qValueFDR') = Pfdr(d); + nidm_peaks.(['Peak_' num2str(i)]).('nidm_Peak__prov_value') = STC_SnPMt(d); + nidm_peaks.(['Peak_' num2str(i)]).('nidm_Peak__nidm_pValueUncorrected') = Pu(d); + nidm_peaks.(['Peak_' num2str(i)]).('nidm_Coordinate__nidm_coordinateVector') = STC_XYZ(1:3,d); end + nidm_cluster.('Peaks') = nidm_peaks; + nidm_clusters.(['Cluster_' num2str(r)]) = nidm_cluster; + bUsed(j) = (bUsed(j) | 1 ); %-Mark maxima as "used" r = r + 1; % Next region - - nidm_cluster.('Peaks') = nidm_peaks; - nidm_clusters(num2str(r)) = nidm_cluster; end nidm_inference.('Clusters') = nidm_clusters; clear i j k D d r - nidm_json.('Inferences') = nidm_inference; - spm_jsonwrite('snpm_nidm.json', nidm_json) + nidm_json.('Inferences').(contrast_id) = nidm_inference; + spm_jsonwrite('snpm_nidm_thresh.json', nidm_json) %-Footnote with SnPM parameters %======================================================================= line([0,1],[0.5,0.5],'LineWidth',1,'Color','r') diff --git a/snpm_ui.m b/snpm_ui.m index dfbab63..bed625f 100644 --- a/snpm_ui.m +++ b/snpm_ui.m @@ -586,7 +586,7 @@ function snpm_ui(varargin) s_SnPMcfg_save = ['s_SnPMcfg_save H C B G HCBGnames P PiCond ',... 'sPiCond bhPerms sHCform iGloNorm sGloNorm GM rg GX GMscale CONT ',... 'THRESH MASK ImMASK TH bVarSm vFWHM sVarSm bVolm bST sDesFile sDesign ',... - 'V pU_ST_Ut df1 nidm_json', ... + 'V pU_ST_Ut df1 nidm_json ', ... 'sDesSave ',sDesSave]; eval(['save SnPMcfg ',s_SnPMcfg_save]) spm_jsonwrite('snpm_nidm.json',nidm_json) From ca62a2a51f6475f58741dce86ccbe57241a67b3f Mon Sep 17 00:00:00 2001 From: cmaumet Date: Tue, 14 Mar 2017 15:32:46 +0000 Subject: [PATCH 05/29] Add permutation terms --- snpm_cp.m | 12 ++++++++++-- snpm_ui.m | 2 +- 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/snpm_cp.m b/snpm_cp.m index bfa04a6..5ccd6a0 100644 --- a/snpm_cp.m +++ b/snpm_cp.m @@ -242,8 +242,11 @@ function snpm_cp(CWD) con_name = 'Positive'; con_neg_name = 'Negative'; STAT = 'T'; - - contrast.('nidm_StatisticMap__nidm_statisticType') = 'obo_tstatistic'; + if bVarSm + contrast.('nidm_StatisticMap__nidm_statisticType') = 'nidm_smoothedtstatistic'; + else + contrast.('nidm_StatisticMap__nidm_statisticType') = 'obo_tstatistic'; + end contrast_pos = contrast; contrast_neg = contrast; contrast_pos.('nidm_ContrastMap__nidm_contrastName') = ['T: ' mat2str(CONT)]; @@ -315,6 +318,11 @@ function snpm_cp(CWD) df = q - r; %-Residual degrees of freedom nPerm = size(PiCond,1); %-# permutations +nidm_json.('Inferences').('nidm_NonParametricNullDistribution__nidm_numberOfPermutations') = nPerm; +nidm_json.('Inferences').('nidm_NonParametricNullDistribution__nidm_hasResamplingScheme') = 'nidm_Permutation'; +nidm_json.('Inferences').('nidm_NonParametricNullDistribution__nidm_hasApproximationMethod') = 'nidm_MonteCarlo'; +nidm_json.('Inferences').('nidm_NonParametricNullDistribution__nidm_maximumNumberOfPermutations') = nPiCond_mx; + con_names = fieldnames(nidm_json.('Contrasts')); for i = 1:numel(con_names) nidm_json.('Contrasts').(con_names{i}).('nidm_StatisticMap__nidm_errorDegreesOfFreedom') = df; diff --git a/snpm_ui.m b/snpm_ui.m index bed625f..23a18d7 100644 --- a/snpm_ui.m +++ b/snpm_ui.m @@ -586,7 +586,7 @@ function snpm_ui(varargin) s_SnPMcfg_save = ['s_SnPMcfg_save H C B G HCBGnames P PiCond ',... 'sPiCond bhPerms sHCform iGloNorm sGloNorm GM rg GX GMscale CONT ',... 'THRESH MASK ImMASK TH bVarSm vFWHM sVarSm bVolm bST sDesFile sDesign ',... - 'V pU_ST_Ut df1 nidm_json ', ... + 'V pU_ST_Ut df1 nidm_json nPiCond_mx ', ... 'sDesSave ',sDesSave]; eval(['save SnPMcfg ',s_SnPMcfg_save]) spm_jsonwrite('snpm_nidm.json',nidm_json) From 6f0666726aad6a4cae08d3fafd7d8306fe37f733 Mon Sep 17 00:00:00 2001 From: cmaumet Date: Tue, 14 Mar 2017 15:34:43 +0000 Subject: [PATCH 06/29] fix: statmap inside contrasts --- snpm_cp.m | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/snpm_cp.m b/snpm_cp.m index 5ccd6a0..fbdfcf7 100644 --- a/snpm_cp.m +++ b/snpm_cp.m @@ -430,13 +430,14 @@ function snpm_cp(CWD) if STAT=='T' VT_pos=snpm_clone_vol(Vt,'snpmT+.img',[str,' (+ve)']); VT_pos=spm_create_vol(VT_pos); - nidm_json.(con_name).('nidm_StatisticMap__prov_atLocation') = 'snpmT+.img'; + nidm_json.('Contrasts').(con_name).('nidm_StatisticMap__prov_atLocation') = 'snpmT+.img'; VT_neg=snpm_clone_vol(Vt,'snpmT-.img',[str,' (-ve)']); VT_neg=spm_create_vol(VT_neg); - nidm_json.(con_neg_name).('nidm_StatisticMap__prov_atLocation') = 'snpmT-.img'; + nidm_json.('Contrasts').(con_neg_name).('nidm_StatisticMap__prov_atLocation') = 'snpmT-.img'; elseif STAT=='F' VF=snpm_clone_vol(Vt,'snpmF.img',str); VF=spm_create_vol(VF); + nidm_json.('Contrasts').(con_name).('nidm_StatisticMap__prov_atLocation') = 'snpmF.img'; end VlP_pos=snpm_clone_vol(Vt, 'lP+.img', '-log10(uncor. non-para. P, +ve)'); From bbfeef16ba06c4111384cae169daa239668f5332 Mon Sep 17 00:00:00 2001 From: cmaumet Date: Thu, 16 Mar 2017 14:09:45 +0000 Subject: [PATCH 07/29] Add error models --- snpm_cp.m | 18 ++++++++++++++++++ snpm_pp.m | 11 ++++++++--- 2 files changed, 26 insertions(+), 3 deletions(-) diff --git a/snpm_cp.m b/snpm_cp.m index fbdfcf7..4b87bf2 100644 --- a/snpm_cp.m +++ b/snpm_cp.m @@ -216,6 +216,24 @@ function snpm_cp(CWD) load(CfgFile); nidm_json = spm_jsonread(NIDMFile); +if strcmp(sDesFile, 'snpm_pi_OneSampT') || ... + strcmp(sDesFile, 'snpm_pi_ANOVAwithinS') + % Sign flipping + nidm_json.('nidm_ErrorModel__nidm_hasErrorDistribution') = {'obo_nonparametricdistribution', 'obo_symmetricdistribution'}; + nidm_json.('nidm_ErrorModel__nidm_errorVarianceHomogeneous') = false; + nidm_json.('nidm_ErrorModel__nidm_varianceMapWiseDependence') = 'nidm_IndependentParameter'; + nidm_json.('nidm_ErrorModel__nidm_hasErrorDependence') = 'nidm_IndependentError'; +else + % Permutation + nidm_json.('nidm_ErrorModel__nidm_hasErrorDistribution') = 'obo_nonparametricdistribution'; + nidm_json.('nidm_ErrorModel__nidm_errorVarianceHomogeneous') = true; + nidm_json.('nidm_ErrorModel__nidm_varianceMapWiseDependence') = 'nidm_IndependentParameter'; + % TODO: the 'obo_exchangeable' term is not yet in STATO + nidm_json.('nidm_ErrorModel__nidm_hasErrorDependence') = 'obo_exchangeable'; + nidm_json.('nidm_ErrorModel__nidm_dependenceMapWiseDependence') = 'nidm_IndependentParameter'; +end + + if isempty([H C]) error('SnPM:NoModel', 'No model specified; [H C] empty'); end diff --git a/snpm_pp.m b/snpm_pp.m index 1d1789a..5ce5e05 100644 --- a/snpm_pp.m +++ b/snpm_pp.m @@ -1178,7 +1178,8 @@ function snpm_pp(CWD,varargin) end % Display only if *not* in command line mode -if true%~spm_get_defaults('cmdline') +nidm=true; +if ~spm_get_defaults('cmdline') || nidm %======================================================================= %-D I S P L A Y : Max report @@ -1431,7 +1432,7 @@ function snpm_pp(CWD,varargin) clear i j k D d r nidm_json.('Inferences').(contrast_id) = nidm_inference; - spm_jsonwrite('snpm_nidm_thresh.json', nidm_json) + %-Footnote with SnPM parameters %======================================================================= line([0,1],[0.5,0.5],'LineWidth',1,'Color','r') @@ -1504,7 +1505,7 @@ function snpm_pp(CWD,varargin) %- Image output? %======================================================================= %-Write out filtered SnPMt? -if WrtFlt +if WrtFlt || nidm Fname = WrtFltFn; @@ -1549,6 +1550,10 @@ function snpm_pp(CWD,varargin) end Vs = sf_close_vol(Vs); clear t + + nidm_json.('Inferences').(contrast_id).(... + 'nidm_ExcursionSetMap__prov_atLocation') = Fname; + spm_jsonwrite('snpm_nidm_thresh.json', nidm_json) end %-Reset Interactive Window From 628d2d9685798fe076d5145f82acd974b4225fef Mon Sep 17 00:00:00 2001 From: cmaumet Date: Fri, 24 Mar 2017 16:55:57 +0000 Subject: [PATCH 08/29] testDataDir should be empty --- test/snpm_test_config.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/snpm_test_config.m b/test/snpm_test_config.m index 97b26be..0f781d2 100644 --- a/test/snpm_test_config.m +++ b/test/snpm_test_config.m @@ -8,6 +8,6 @@ % Camille Maumet global testDataDir; -testDataDir = '/Users/cmaumet/Desktop/SnPM_test_data'; % Specify directory containing test data +testDataDir = ''; % Specify directory containing test data global SnPMrefVersion SnPMrefVersion = 'SnPM8'; \ No newline at end of file From 2a6a7776be448463b349eb9c17e10492440f6228 Mon Sep 17 00:00:00 2001 From: cmaumet Date: Wed, 29 Mar 2017 14:02:31 +0100 Subject: [PATCH 09/29] Add test functions to path when loading SnPM toolbox --- tbx_cfg_snpm13.m | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tbx_cfg_snpm13.m b/tbx_cfg_snpm13.m index 5ab8c24..fbd05f0 100644 --- a/tbx_cfg_snpm13.m +++ b/tbx_cfg_snpm13.m @@ -15,6 +15,8 @@ addpath(toolboxDir); addpath(fullfile(toolboxDir, 'config')); +addpath(fullfile(toolboxDir, 'test')); +addpath(fullfile(toolboxDir, 'test', 'common')); snpmBatch = snpm_cfg_master; end \ No newline at end of file From 4c0d1c1b4511daf7d96e71280c03d39912a35357 Mon Sep 17 00:00:00 2001 From: cmaumet Date: Fri, 31 Mar 2017 16:02:11 +0100 Subject: [PATCH 10/29] Use map rather than struct --- snpm_cp.m | 113 +++++++++++++++++++++++++++++++++--------------------- snpm_pp.m | 77 ++++++++++++++++++++++++------------- snpm_ui.m | 13 +++---- 3 files changed, 126 insertions(+), 77 deletions(-) diff --git a/snpm_cp.m b/snpm_cp.m index 4b87bf2..fd5cbe3 100644 --- a/snpm_cp.m +++ b/snpm_cp.m @@ -209,28 +209,26 @@ function snpm_cp(CWD) fprintf('Changing directory to %s\n',CWD); end CfgFile = fullfile(CWD,'SnPMcfg.mat'); -NIDMFile = fullfile(CWD,'snpm_nidm.json'); %-Load config file & catch all problem cases now %----------------------------------------------------------------------- load(CfgFile); -nidm_json = spm_jsonread(NIDMFile); if strcmp(sDesFile, 'snpm_pi_OneSampT') || ... strcmp(sDesFile, 'snpm_pi_ANOVAwithinS') % Sign flipping - nidm_json.('nidm_ErrorModel__nidm_hasErrorDistribution') = {'obo_nonparametricdistribution', 'obo_symmetricdistribution'}; - nidm_json.('nidm_ErrorModel__nidm_errorVarianceHomogeneous') = false; - nidm_json.('nidm_ErrorModel__nidm_varianceMapWiseDependence') = 'nidm_IndependentParameter'; - nidm_json.('nidm_ErrorModel__nidm_hasErrorDependence') = 'nidm_IndependentError'; + nidm_json('nidm_ErrorModel/nidm_hasErrorDistribution') = {'obo_nonparametricdistribution', 'obo_symmetricdistribution'}; + nidm_json('nidm_ErrorModel/nidm_errorVarianceHomogeneous') = false; + nidm_json('nidm_ErrorModel/nidm_varianceMapWiseDependence') = 'nidm_IndependentParameter'; + nidm_json('nidm_ErrorModel/nidm_hasErrorDependence') = 'nidm_IndependentError'; else % Permutation - nidm_json.('nidm_ErrorModel__nidm_hasErrorDistribution') = 'obo_nonparametricdistribution'; - nidm_json.('nidm_ErrorModel__nidm_errorVarianceHomogeneous') = true; - nidm_json.('nidm_ErrorModel__nidm_varianceMapWiseDependence') = 'nidm_IndependentParameter'; + nidm_json('nidm_ErrorModel/nidm_hasErrorDistribution') = 'obo_nonparametricdistribution'; + nidm_json('nidm_ErrorModel/nidm_errorVarianceHomogeneous') = true; + nidm_json('nidm_ErrorModel/nidm_varianceMapWiseDependence') = 'nidm_IndependentParameter'; % TODO: the 'obo_exchangeable' term is not yet in STATO - nidm_json.('nidm_ErrorModel__nidm_hasErrorDependence') = 'obo_exchangeable'; - nidm_json.('nidm_ErrorModel__nidm_dependenceMapWiseDependence') = 'nidm_IndependentParameter'; + nidm_json('nidm_ErrorModel/nidm_hasErrorDependence') = 'obo_exchangeable'; + nidm_json('nidm_ErrorModel/nidm_dependenceMapWiseDependence') = 'nidm_IndependentParameter'; end @@ -243,37 +241,44 @@ function snpm_cp(CWD) if size(CONT,2) ~= size([H C B G],2) error('SnPM:InvalidContrast','Contrast problem; wrong number of columns'); end + if size(CONT,1) > 1 - nidm_json.('Contrasts') = struct(); - contrast = struct(); - contrast.('obo_contrastweightmatrix__prov_value') = CONT; + contrastMap = containers.Map(); + contrastMap('obo_contrastweightmatrix/prov:value') = CONT; warning('SnPM:FContrast', ... 'F contrast! F statistic images are being created.'); STAT = 'F'; - contrast.('nidm_StatisticMap__nidm_statisticType') = 'obo_Fstatistic'; + contrastMap('nidm_StatisticMap/nidm_statisticType') = 'obo_Fstatistic'; if (CONT(1,:) == -CONT(2,:)) CONT = CONT(1,:); end con_name = 'Positive'; - nidm_json.('Contrasts').(con_name) = contrast; + + nidm_json('Contrasts') = containers.Map({con_name}, {contrastMap}); else + contrast_pos = containers.Map(); + contrast_neg = containers.Map(); % There is no deep copy of + % containers.Map so we need a separate object from the start con_name = 'Positive'; con_neg_name = 'Negative'; STAT = 'T'; if bVarSm - contrast.('nidm_StatisticMap__nidm_statisticType') = 'nidm_smoothedtstatistic'; + contrast_pos('nidm_StatisticMap/nidm_statisticType') = 'nidm_smoothedtstatistic'; + contrast_neg('nidm_StatisticMap/nidm_statisticType') = 'nidm_smoothedtstatistic'; else - contrast.('nidm_StatisticMap__nidm_statisticType') = 'obo_tstatistic'; + contrast_pos('nidm_StatisticMap/nidm_statisticType') = 'obo_tstatistic'; + contrast_neg('nidm_StatisticMap/nidm_statisticType') = 'obo_tstatistic'; end - contrast_pos = contrast; - contrast_neg = contrast; - contrast_pos.('nidm_ContrastMap__nidm_contrastName') = ['T: ' mat2str(CONT)]; - contrast_pos.('obo_contrastweightmatrix__prov_value') = CONT; - contrast_neg.('nidm_ContrastMap__nidm_contrastName') = ['T: ' mat2str(-CONT)]; - contrast_neg.('obo_contrastweightmatrix__prov_value') = -CONT; + contrast_pos('nidm_ContrastMap/nidm_contrastName') = ['Positive T-Contrast: [' mat2str(CONT) ']']; + contrast_pos('obo_contrastweightmatrix/prov:value') = CONT; + contrast_neg('nidm_ContrastMap/nidm_contrastName') = ['Negative T-Contrast: [' mat2str(-CONT) ']']; + contrast_neg('obo_contrastweightmatrix/prov:value') = -CONT; + +% nidm_json('Contrasts').(con_name) = contrast_pos; +% nidm_json('Contrasts').(con_neg_name) = contrast_neg; - nidm_json.('Contrasts').(con_name) = contrast_pos; - nidm_json.('Contrasts').(con_neg_name) = contrast_neg; + nidm_json('Contrasts') = containers.Map({con_name, con_neg_name}, ... + {contrast_pos, contrast_neg}); end if rank(CONT) 3) - nidm_cluster = struct(); - nidm_peaks = struct(); + nidm_cluster = containers.Map(); + nidm_peaks = containers.Map(); [null, i] = max(STC_SnPMt.*(~bUsed)); % Largest t value j = find(STC_r == STC_r(i)); % Maxima in same region @@ -1383,18 +1384,25 @@ function snpm_pp(CWD,varargin) text(tCol(9),y,sprintf(Fmtst{9},STC_XYZ(2,i)),'UserData',STC_XYZ(:,i),StrAttrB{:}) text(tCol(10),y,sprintf(Fmtst{10},STC_XYZ(3,i)),'UserData',STC_XYZ(:,i),StrAttrB{:}) - nidm_cluster.('nidm_SupraThresholdCluster__nidm_clusterSizeInVoxels') = STC_N(i); - nidm_peaks.('Peak_1').('nidm_Peak__nidm_pValueFWER') = Pt(i); - nidm_peaks.('Peak_1').('nidm_Peak__nidm_qValueFDR') = Pfdr(i); - nidm_peaks.('Peak_1').('nidm_Peak__prov_value') = STC_SnPMt(i); - nidm_peaks.('Peak_1').('nidm_Peak__nidm_pValueUncorrected') = Pu(i); - nidm_peaks.('Peak_1').('nidm_Coordinate__nidm_coordinateVector') = STC_XYZ(1:3,i); + nidm_cluster('nidm_SupraThresholdCluster/nidm_clusterSizeInVoxels') = STC_N(i); + nidm_peaks('Peak_1') = containers.Map(... + { ... + 'nidm_Peak/nidm_pValueFWER', ... + 'nidm_Peak/nidm_qValueFDR', ... + 'nidm_Peak/prov:value', ... + 'nidm_Peak/nidm_pValueUncorrected', ... + 'nidm_Coordinate/nidm_coordinateVector', ... + }, ... + {Pt(i), Pfdr(i), STC_SnPMt(i), Pu(i), STC_XYZ(1:3,i)}); y = y -1; %-Print up to 3 secondary maxima (>8mm apart) - %------------------------------------------------------------------- - nidm_json.('Inferences').('nidm_PeakDefinitionCriteria__nidm_minDistanceBetweenPeaks') = 8; + %------------------------------------------------------------------- + nidm_inferences = nidm_json('Inferences'); + nidm_inferences('nidm_PeakDefinitionCriteria/nidm_minDistanceBetweenPeaks') = 8; + nidm_json('Inferences') = nidm_inferences; + [null, k] = sort(-STC_SnPMt(j)); % Sort on t value D = i; for i = 1:length(k) @@ -1415,23 +1423,33 @@ function snpm_pp(CWD,varargin) y = y -1; end end - nidm_peaks.(['Peak_' num2str(i)]).('nidm_Peak__nidm_pValueFWER') = Pt(d); - nidm_peaks.(['Peak_' num2str(i)]).('nidm_Peak__nidm_qValueFDR') = Pfdr(d); - nidm_peaks.(['Peak_' num2str(i)]).('nidm_Peak__prov_value') = STC_SnPMt(d); - nidm_peaks.(['Peak_' num2str(i)]).('nidm_Peak__nidm_pValueUncorrected') = Pu(d); - nidm_peaks.(['Peak_' num2str(i)]).('nidm_Coordinate__nidm_coordinateVector') = STC_XYZ(1:3,d); + + nidm_peaks(['Peak_' num2str(i)]) = containers.Map(... + { ... + 'nidm_Peak/nidm_pValueFWER', ... + 'nidm_Peak/nidm_qValueFDR', ... + 'nidm_Peak/prov:value', ... + 'nidm_Peak/nidm_pValueUncorrected', ... + 'nidm_Coordinate/nidm_coordinateVector', ... + }, ... + {Pt(d), Pfdr(d), STC_SnPMt(d), Pu(d), STC_XYZ(1:3,d)}); + end - nidm_cluster.('Peaks') = nidm_peaks; - nidm_clusters.(['Cluster_' num2str(r)]) = nidm_cluster; + nidm_cluster('Peaks') = nidm_peaks; + nidm_clusters(['Cluster_' num2str(r)]) = nidm_cluster; bUsed(j) = (bUsed(j) | 1 ); %-Mark maxima as "used" r = r + 1; % Next region end - nidm_inference.('Clusters') = nidm_clusters; - clear i j k D d r - nidm_json.('Inferences').(contrast_id) = nidm_inference; + nidm_inferences = nidm_json('Inferences'); + nidm_inference = containers.Map(); + nidm_inference('Clusters') = nidm_clusters; + nidm_inferences(contrast_id) = nidm_inference; + nidm_json('Inferences') = nidm_inferences; + + clear i j k D d r %-Footnote with SnPM parameters %======================================================================= @@ -1551,9 +1569,16 @@ function snpm_pp(CWD,varargin) Vs = sf_close_vol(Vs); clear t - nidm_json.('Inferences').(contrast_id).(... - 'nidm_ExcursionSetMap__prov_atLocation') = Fname; - spm_jsonwrite('snpm_nidm_thresh.json', nidm_json) + + nidm_inferences = nidm_json('Inferences'); + nidm_inference = nidm_inferences(contrast_id); + nidm_inference('nidm_ExcursionSetMap/prov:atLocation') = Fname; + nidm_inferences(contrast_id) = nidm_inference; + nidm_json('Inferences') = nidm_inferences; + + % TODO: This temp file should only be produced if NIDM export is requested + jsonwrite('snpm_nidm_thresh.json', nidm_json, ... + struct('indent',' ', 'escape', false)); end %-Reset Interactive Window diff --git a/snpm_ui.m b/snpm_ui.m index 23a18d7..2ce66c8 100644 --- a/snpm_ui.m +++ b/snpm_ui.m @@ -193,7 +193,7 @@ function snpm_ui(varargin) cd(job.dir{1}) end -nidm_json = struct(); +nidm_json = containers.Map(); %-Definitions & Design parameters %======================================================================= @@ -372,7 +372,7 @@ function snpm_ui(varargin) end if (iGMsca==2) % CHANGED from 1 to 2 as should not ask for a value if grand mean scaling is not required. - nidm_json.('nidm_Data__nidm_grandMeanScaling') = true; + nidm_json('nidm_Data/nidm_grandMeanScaling') = true; if (iGloNorm==2) % Proportional scaling str = 'PropSca global mean to'; else @@ -386,10 +386,10 @@ function snpm_ui(varargin) case 'gmsca_no', GM = 50; end - nidm_json.('nidm_Data__nidm_targetIntensity') = GM; + nidm_json('nidm_Data/nidm_targetIntensity') = GM; elseif (iGMsca==1) % No grand mean scaling GM = 0; - nidm_json.('nidm_Data__nidm_grandMeanScaling') = false; + nidm_json('nidm_Data/nidm_grandMeanScaling') = false; end @@ -579,17 +579,16 @@ function snpm_ui(varargin) %-Construct full design matrix and name matrices for display %----------------------------------------------------------------------- [nHCBG,HCBGnames] = spm_DesMtx('Sca',H,Hnames,C,Cnames,B,Bnames,G,Gnames); -nidm_json.('nidm_DesignMatrix__prov_value') = nHCBG; +nidm_json('nidm_DesignMatrix/prov_value') = nHCBG; %-Setup is complete - save SnPMcfg Mat file %----------------------------------------------------------------------- s_SnPMcfg_save = ['s_SnPMcfg_save H C B G HCBGnames P PiCond ',... 'sPiCond bhPerms sHCform iGloNorm sGloNorm GM rg GX GMscale CONT ',... 'THRESH MASK ImMASK TH bVarSm vFWHM sVarSm bVolm bST sDesFile sDesign ',... - 'V pU_ST_Ut df1 nidm_json nPiCond_mx ', ... + 'V pU_ST_Ut df1 nidm_json nPiCond_mx nidm_json ', ... 'sDesSave ',sDesSave]; eval(['save SnPMcfg ',s_SnPMcfg_save]) -spm_jsonwrite('snpm_nidm.json',nidm_json) if ~spm_get_defaults('cmdline') From 8287c418732b4b1a04a0177bd027f186853d4cac Mon Sep 17 00:00:00 2001 From: cmaumet Date: Thu, 20 Apr 2017 14:22:01 +0100 Subject: [PATCH 11/29] Save regressor names to nidm_json --- snpm_ui.m | 1 + 1 file changed, 1 insertion(+) diff --git a/snpm_ui.m b/snpm_ui.m index 2ce66c8..f44cd40 100644 --- a/snpm_ui.m +++ b/snpm_ui.m @@ -580,6 +580,7 @@ function snpm_ui(varargin) %----------------------------------------------------------------------- [nHCBG,HCBGnames] = spm_DesMtx('Sca',H,Hnames,C,Cnames,B,Bnames,G,Gnames); nidm_json('nidm_DesignMatrix/prov_value') = nHCBG; +nidm_json('nidm_DesignMatrix/nidm_regressorNames') = HCBGnames; %-Setup is complete - save SnPMcfg Mat file %----------------------------------------------------------------------- From c2fb4ee4e09bd3c8bb4ef4ce11dc0ed8a60e85ac Mon Sep 17 00:00:00 2001 From: cmaumet Date: Thu, 20 Apr 2017 14:22:11 +0100 Subject: [PATCH 12/29] Use map rather than structure for nidm_json --- snpm_pp.m | 48 ++++++++++++++++++++++++------------------------ snpm_ui.m | 2 +- 2 files changed, 25 insertions(+), 25 deletions(-) diff --git a/snpm_pp.m b/snpm_pp.m index 4d991a9..b511107 100644 --- a/snpm_pp.m +++ b/snpm_pp.m @@ -490,8 +490,8 @@ function snpm_pp(CWD,varargin) alph_FWE = NaN; % FWE rate of a specified u threshold alph_FDR = NaN; % FDR rate of a specified alpha_ucp -nidm_inference = struct(); % NIDM structre storing information about inference -nidm_inference.('nidm_Inference__nidm_hasAlternativeHypothesis') = 'nidm_OneTailedTest'; +nidm_inference = containers.Map(); % Map storing information about inference +nidm_inference('nidm_Inference/nidm_hasAlternativeHypothesis') = 'nidm_OneTailedTest'; if BATCH bSpatEx = isfield(job.Thr,'Clus'); if ~bSpatEx @@ -501,18 +501,18 @@ function snpm_pp(CWD,varargin) case 'Pth' alpha_ucp = BoundCheck(job.Thr.Vox.VoxSig.Pth,[0 1],'Invalid Uncorrected P'); alph_FDR = snpm_P_FDR(alpha_ucp,[],'P',[],sSnPMucp'); - nidm_inference.('nidm_HeightThreshold__prov_type') = 'nidm_PValueUncorrected'; - nidm_inference.('nidm_HeightThreshold__prov_value') = alpha_ucp; + nidm_inference('nidm_HeightThreshold/prov:type') = 'nidm_PValueUncorrected'; + nidm_inference('nidm_HeightThreshold/prov:value') = alpha_ucp; case 'TFth' u = BoundCheck(job.Thr.Vox.VoxSig.TFth,[0 Inf],'Negative Threshold!'); alph_FWE = sum(MaxT > u -tol) / nPermReal; - nidm_inference.('nidm_HeightThreshold__prov_type') = 'obo_statistic'; - nidm_inference.('nidm_HeightThreshold__prov_value') = u; + nidm_inference('nidm_HeightThreshold/prov:type') = 'obo_statistic'; + nidm_inference('nidm_HeightThreshold/prov:value') = u; case 'FDRth' alph_FDR = BoundCheck(job.Thr.Vox.VoxSig.FDRth,[0 1],'Invalid FDR level'); alpha_ucp = snpm_uc_FDR(alph_FDR,[],'P',[],sSnPMucp'); - nidm_inference.('nidm_HeightThreshold__prov_type') = 'obo_qvalue'; - nidm_inference.('nidm_HeightThreshold__prov_value') = alph_FDR; + nidm_inference('nidm_HeightThreshold/prov:type') = 'obo_qvalue'; + nidm_inference('nidm_HeightThreshold/prov:value') = alph_FDR; case 'FWEth' alph_FWE = BoundCheck(job.Thr.Vox.VoxSig.FWEth,[0 1],'Invalid FWE level'); iFWE = ceil((1-alph_FWE)*nPermReal); @@ -522,12 +522,12 @@ function snpm_pp(CWD,varargin) C_MaxT = 0; end u = C_MaxT; - nidm_inference.('nidm_HeightThreshold__prov_type') = 'obo_FWERadjustedpvalue'; - nidm_inference.('nidm_HeightThreshold__prov_value') = alph_FWE; + nidm_inference('nidm_HeightThreshold/prov:type') = 'obo_FWERadjustedpvalue'; + nidm_inference('nidm_HeightThreshold/prov:value') = alph_FWE; end % No extent thresholding when voxelwise threshold is requested - nidm_inference.('nidm_ExtentThreshold__prov_type') = 'obo_statistic'; - nidm_inference.('nidm_ExtentThreshold__prov_value') = 0; + nidm_inference('nidm_ExtentThreshold/prov:type') = 'obo_statistic'; + nidm_inference('nidm_ExtentThreshold/prov:value') = 0; else % Cluster-wise inference if exist(fullfile(CWD,'SnPM_ST.mat'))~=2 & exist(fullfile(CWD,'STCS.mat'))~=2 @@ -543,11 +543,11 @@ function snpm_pp(CWD,varargin) ST_Ut_0 = ST_Ut; CFth=job.Thr.Clus.ClusSize.CFth; if (CFth<1) - nidm_inference.('nidm_ExtentThreshold__prov_type') = 'nidm_PValueUncorrected'; - nidm_inference.('nidm_ExtentThreshold__prov_value') = CFth; + nidm_inference('nidm_ExtentThreshold/prov:type') = 'nidm_PValueUncorrected'; + nidm_inference('nidm_ExtentThreshold/prov:value') = CFth; else - nidm_inference.('nidm_ExtentThreshold__prov_type') = 'obo_statistic'; - nidm_inference.('nidm_ExtentThreshold__prov_value') = CFth; + nidm_inference('nidm_ExtentThreshold/prov:type') = 'obo_statistic'; + nidm_inference('nidm_ExtentThreshold/prov:value') = CFth; end if (CFth<=0) @@ -612,17 +612,17 @@ function snpm_pp(CWD,varargin) switch tmp{1} case 'Cth' C_STCS = job.Thr.Clus.ClusSize.ClusSig.Cth; - nidm_inference.('nidm_ExtentThreshold__prov_type') = 'obo_statistic'; - nidm_inference.('nidm_ExtentThreshold__prov_value') = C_STCS; + nidm_inference('nidm_ExtentThreshold/prov:type') = 'obo_statistic'; + nidm_inference('nidm_ExtentThreshold/prov:value') = C_STCS; case 'PthC' alpha_ucp = BoundCheck(job.Thr.Clus.ClusSize.ClusSig.PthC,[0 1],'Invalid uncorrected P(k)'); - nidm_inference.('nidm_ExtentThreshold__prov_type') = 'nidm_PValueUncorrected'; - nidm_inference.('nidm_ExtentThreshold__prov_value') = alpha_ucp; + nidm_inference('nidm_ExtentThreshold/prov:type') = 'nidm_PValueUncorrected'; + nidm_inference('nidm_ExtentThreshold/prov:value') = alpha_ucp; case 'FWEthC' alph_FWE = BoundCheck(job.Thr.Clus.ClusSize.ClusSig.FWEthC,[0 1],'Invalid FWE level (cluster-level inference)'); iFWE = ceil((1-alph_FWE)*nPermReal); - nidm_inference.('nidm_ExtentThreshold__prov_type') = 'obo_FWERadjustedpvalue'; - nidm_inference.('nidm_ExtentThreshold__prov_value') = alph_FWE; + nidm_inference('nidm_ExtentThreshold/prov:type') = 'obo_FWERadjustedpvalue'; + nidm_inference('nidm_ExtentThreshold/prov:value') = alph_FWE; end end % END: Cluster-wise inference @@ -930,9 +930,9 @@ function snpm_pp(CWD,varargin) %-Save perm 1 stats for use later - [X;Y;Z;T;perm;STCno] tmp = spm_clusters(Locs_vox(1:3,:)); - nidm_inferences = nidm_json.('Inferences'); + nidm_inferences = nidm_json('Inferences'); nidm_inferences('nidm_ClusterDefinitionCriteria/nidm_hasConnectivityCriterion') = 'nidm_voxel18connected'; - nidm_json.('Inferences') = nidm_inferences; + nidm_json('Inferences') = nidm_inferences; if isPos==1 STCstats_Pos = [ SnPM_ST(:,tQ); tmp]; if bNeg==0 diff --git a/snpm_ui.m b/snpm_ui.m index f44cd40..5f1484a 100644 --- a/snpm_ui.m +++ b/snpm_ui.m @@ -579,7 +579,7 @@ function snpm_ui(varargin) %-Construct full design matrix and name matrices for display %----------------------------------------------------------------------- [nHCBG,HCBGnames] = spm_DesMtx('Sca',H,Hnames,C,Cnames,B,Bnames,G,Gnames); -nidm_json('nidm_DesignMatrix/prov_value') = nHCBG; +nidm_json('nidm_DesignMatrix/prov:value') = nHCBG; nidm_json('nidm_DesignMatrix/nidm_regressorNames') = HCBGnames; %-Setup is complete - save SnPMcfg Mat file From f5916cfcb0f0da54a9ac2a8244d490466650ac3b Mon Sep 17 00:00:00 2001 From: cmaumet Date: Thu, 20 Apr 2017 14:22:20 +0100 Subject: [PATCH 13/29] Save contrast images (and store in nidm_json) --- snpm_cp.m | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/snpm_cp.m b/snpm_cp.m index fd5cbe3..27757fe 100644 --- a/snpm_cp.m +++ b/snpm_cp.m @@ -459,17 +459,27 @@ function snpm_cp(CWD) if STAT=='T' VT_pos=snpm_clone_vol(Vt,'snpmT+.img',[str,' (+ve)']); VT_pos=spm_create_vol(VT_pos); + + VCON_pos=snpm_clone_vol(Vt,'con+.img',[str,' (+ve)']); + VCON_pos=spm_create_vol(VCON_pos); + nidm_contrasts = nidm_json('Contrasts'); nidm_contrast = nidm_contrasts(con_name); + + nidm_contrast('nidm_ContrastMap/prov:atLocation') = 'con+.img'; nidm_contrast('nidm_StatisticMap/prov:atLocation') = 'snpmT+.img'; nidm_contrasts(con_name) = nidm_contrast; % nidm_json('Contrasts').(con_name).('nidm_StatisticMap/prov:atLocation') = 'snpmT+.img'; VT_neg=snpm_clone_vol(Vt,'snpmT-.img',[str,' (-ve)']); VT_neg=spm_create_vol(VT_neg); + + VCON_neg=snpm_clone_vol(Vt,'con-.img',[str,' (+ve)']); + VCON_neg=spm_create_vol(VCON_neg); % nidm_json('Contrasts').(con_neg_name).('nidm_StatisticMap__prov_atLocation') = 'snpmT-.img'; nidm_neg_contrast = nidm_contrasts(con_neg_name); nidm_neg_contrast('nidm_StatisticMap/prov:atLocation') = 'snpmT-.img'; + nidm_neg_contrast('nidm_ContrastMap/prov:atLocation') = 'con-.img'; nidm_contrasts(con_neg_name) = nidm_neg_contrast; nidm_json('Contrasts') = nidm_contrasts; @@ -537,6 +547,8 @@ function snpm_cp(CWD) if STAT=='T' T_pos_image=repmat(NaN,1,WorkDim); T_neg_image=repmat(NaN,1,WorkDim); + CON_pos_image=repmat(NaN,1,WorkDim); + CON_neg_image=repmat(NaN,1,WorkDim); elseif STAT=='F' F_image=repmat(NaN,1,WorkDim); end @@ -656,9 +668,11 @@ function snpm_cp(CWD) %-Compute t-statistics for specified compounds of parameters %----------------------------------------------------------- T = zeros(1,size(BETA,2)); + CON = zeros(1,size(BETA,2)); Co = CONT; if STAT=='T' % t, as usual + CON(1,:) = Co*BETA; T(1,:) = Co*BETA./sqrt((ResSS*(Co*pinv([H C B G]'*[H C B G])*Co'))/df); else % F! @@ -699,6 +713,9 @@ function snpm_cp(CWD) if STAT=='T' T_pos_image(:,Q)=T; T_neg_image(:,Q)=-T; + + CON_pos_image(:,Q)=CON; + CON_neg_image(:,Q)=-CON; elseif STAT=='F' F_image(:,Q)=T; end @@ -739,6 +756,12 @@ function snpm_cp(CWD) T_neg_vol=reshape(T_neg_image,DIM(1),DIM(2),DIM(3)); spm_write_vol(VT_neg,T_neg_vol); + + CON_pos_vol=reshape(CON_pos_image,DIM(1),DIM(2),DIM(3)); + spm_write_vol(VCON_pos,CON_pos_vol); + + CON_neg_vol=reshape(CON_neg_image,DIM(1),DIM(2),DIM(3)); + spm_write_vol(VCON_neg,CON_neg_vol); elseif STAT=='F' F_vol=reshape(F_image,DIM(1),DIM(2),DIM(3)); From 56306a4ccc489d1e0aa5d5b8ea9b61ac814e665a Mon Sep 17 00:00:00 2001 From: cmaumet Date: Thu, 20 Apr 2017 14:31:02 +0100 Subject: [PATCH 14/29] Save contrast std error image --- snpm_cp.m | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/snpm_cp.m b/snpm_cp.m index 27757fe..96b04af 100644 --- a/snpm_cp.m +++ b/snpm_cp.m @@ -463,9 +463,14 @@ function snpm_cp(CWD) VCON_pos=snpm_clone_vol(Vt,'con+.img',[str,' (+ve)']); VCON_pos=spm_create_vol(VCON_pos); + VCONSE=snpm_clone_vol(Vt,'conse.img', str); + VCONSE=spm_create_vol(VCONSE); + nidm_contrasts = nidm_json('Contrasts'); nidm_contrast = nidm_contrasts(con_name); + nidm_contrast('nidm_ContrastStandardErrorMap/prov:atLocation') = 'conse.img'; + nidm_contrast('nidm_ContrastMap/prov:atLocation') = 'con+.img'; nidm_contrast('nidm_StatisticMap/prov:atLocation') = 'snpmT+.img'; nidm_contrasts(con_name) = nidm_contrast; @@ -480,6 +485,10 @@ function snpm_cp(CWD) nidm_neg_contrast = nidm_contrasts(con_neg_name); nidm_neg_contrast('nidm_StatisticMap/prov:atLocation') = 'snpmT-.img'; nidm_neg_contrast('nidm_ContrastMap/prov:atLocation') = 'con-.img'; + + % TODO: check it's fine to use same conse image for pos and neg contrasts + nidm_neg_contrast('nidm_ContrastStandardErrorMap/prov:atLocation') = 'conse.img'; + nidm_contrasts(con_neg_name) = nidm_neg_contrast; nidm_json('Contrasts') = nidm_contrasts; @@ -671,8 +680,9 @@ function snpm_cp(CWD) CON = zeros(1,size(BETA,2)); Co = CONT; if STAT=='T' - % t, as usual CON(1,:) = Co*BETA; + CONSE(1,:) = sqrt((ResSS*(Co*pinv([H C B G]'*[H C B G])*Co'))/df); + % t, as usual T(1,:) = Co*BETA./sqrt((ResSS*(Co*pinv([H C B G]'*[H C B G])*Co'))/df); else % F! @@ -716,6 +726,8 @@ function snpm_cp(CWD) CON_pos_image(:,Q)=CON; CON_neg_image(:,Q)=-CON; + + CONSE_image(:,Q)=CONSE; elseif STAT=='F' F_image(:,Q)=T; end @@ -763,6 +775,9 @@ function snpm_cp(CWD) CON_neg_vol=reshape(CON_neg_image,DIM(1),DIM(2),DIM(3)); spm_write_vol(VCON_neg,CON_neg_vol); + CONSE_vol=reshape(CONSE_image,DIM(1),DIM(2),DIM(3)); + spm_write_vol(VCONSE,CONSE_vol); + elseif STAT=='F' F_vol=reshape(F_image,DIM(1),DIM(2),DIM(3)); spm_write_vol(VF,F_vol); From 0d27f51496fee6839ec1ce09d756adce8c1904d0 Mon Sep 17 00:00:00 2001 From: cmaumet Date: Thu, 20 Apr 2017 15:17:47 +0100 Subject: [PATCH 15/29] Cluster forming threshold is height threshold --- snpm_pp.m | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/snpm_pp.m b/snpm_pp.m index b511107..440e6e9 100644 --- a/snpm_pp.m +++ b/snpm_pp.m @@ -524,6 +524,8 @@ function snpm_pp(CWD,varargin) u = C_MaxT; nidm_inference('nidm_HeightThreshold/prov:type') = 'obo_FWERadjustedpvalue'; nidm_inference('nidm_HeightThreshold/prov:value') = alph_FWE; + otherwise + error('Unknown threshold') end % No extent thresholding when voxelwise threshold is requested nidm_inference('nidm_ExtentThreshold/prov:type') = 'obo_statistic'; @@ -543,11 +545,11 @@ function snpm_pp(CWD,varargin) ST_Ut_0 = ST_Ut; CFth=job.Thr.Clus.ClusSize.CFth; if (CFth<1) - nidm_inference('nidm_ExtentThreshold/prov:type') = 'nidm_PValueUncorrected'; - nidm_inference('nidm_ExtentThreshold/prov:value') = CFth; + nidm_inference('nidm_HeightThreshold/prov:type') = 'nidm_PValueUncorrected'; + nidm_inference('nidm_HeightThreshold/prov:value') = CFth; else - nidm_inference('nidm_ExtentThreshold/prov:type') = 'obo_statistic'; - nidm_inference('nidm_ExtentThreshold/prov:value') = CFth; + nidm_inference('nidm_HeightThreshold/prov:type') = 'obo_statistic'; + nidm_inference('nidm_HeightThreshold/prov:value') = CFth; end if (CFth<=0) @@ -625,7 +627,7 @@ function snpm_pp(CWD,varargin) nidm_inference('nidm_ExtentThreshold/prov:value') = alph_FWE; end end % END: Cluster-wise inference - + else % GUI, interative inference specification str_img =[STAT,'|P']; From 038165dd1d3a7c566267fb0480bedc4d457891a3 Mon Sep 17 00:00:00 2001 From: cmaumet Date: Thu, 20 Apr 2017 16:15:29 +0100 Subject: [PATCH 16/29] Add missing NIDM info --- snpm_cp.m | 31 ++++++++++++++++++++----- snpm_pp.m | 69 +++++++++++++++++++++++++++++++++++++++++-------------- 2 files changed, 77 insertions(+), 23 deletions(-) diff --git a/snpm_cp.m b/snpm_cp.m index 96b04af..04c37d9 100644 --- a/snpm_cp.m +++ b/snpm_cp.m @@ -231,6 +231,8 @@ function snpm_cp(CWD) nidm_json('nidm_ErrorModel/nidm_dependenceMapWiseDependence') = 'nidm_IndependentParameter'; end +% TODO: check this is correct +nidm_json('nidm_ModelParameterEstimation/nidm_withEstimationMethod') = 'obo_ordinaryleastsquaresestimation'; if isempty([H C]) error('SnPM:NoModel', 'No model specified; [H C] empty'); @@ -341,12 +343,10 @@ function snpm_cp(CWD) df = q - r; %-Residual degrees of freedom nPerm = size(PiCond,1); %-# permutations -nidm_json('Inferences') = containers.Map(... - {'nidm_NonParametricNullDistribution/nidm_numberOfPermutations',... - 'nidm_NonParametricNullDistribution/nidm_hasResamplingScheme',... - 'nidm_NonParametricNullDistribution/nidm_hasApproximationMethod',... - 'nidm_NonParametricNullDistribution/nidm_maximumNumberOfPermutations'},... - { nPerm, 'nidm_Permutation', 'nidm_MonteCarlo', nPiCond_mx}); +nidm_json('nidm_NonParametricNullDistribution/nidm_numberOfPermutations') = nPerm; +nidm_json('nidm_NonParametricNullDistribution/nidm_hasResamplingScheme') = 'nidm_Permutation'; +nidm_json('nidm_NonParametricNullDistribution/nidm_hasApproximationMethod') = 'nidm_MonteCarlo'; +nidm_json('nidm_NonParametricNullDistribution/nidm_maximumNumberOfPermutations') = nPiCond_mx; con_names = nidm_json('Contrasts').keys; for i = 1:numel(con_names) @@ -456,6 +456,17 @@ function snpm_cp(CWD) [rank(CONT) df],vFWHM); end end + +Vmask=snpm_clone_vol(Vt,'mask.img',str); +Vmask=spm_create_vol(Vmask); +nidm_json('nidm_MaskMap/prov:atLocation') = 'mask.img'; + +% TODO: grand mean needs to be computed +warning('Grand mean is not computed') +Vgm=snpm_clone_vol(Vt,'gm.img',str); +Vgm=spm_create_vol(Vgm); +nidm_json('nidm_GrandMeanMap/prov:atLocation') = 'gm.img'; + if STAT=='T' VT_pos=snpm_clone_vol(Vt,'snpmT+.img',[str,' (+ve)']); VT_pos=spm_create_vol(VT_pos); @@ -614,6 +625,7 @@ function snpm_cp(CWD) else Q = find(all(X>TH) & any(diff(X)) & Wt); end + if length(Q) @@ -762,6 +774,13 @@ function snpm_cp(CWD) ResSS_vol=reshape(ResSS_image,DIM(1),DIM(2),DIM(3)); spm_write_vol(VResMS, ResSS_vol); + % Analysis mask + mask_vol=reshape(Q,DIM(1),DIM(2),DIM(3)); + spm_write_vol(Vmask,mask_vol); + + % TODO: replace mask volume by computed grand mean + spm_write_vol(Vgm,mask_vol); + if STAT=='T' T_pos_vol=reshape(T_pos_image,DIM(1),DIM(2),DIM(3)); spm_write_vol(VT_pos,T_pos_vol); diff --git a/snpm_pp.m b/snpm_pp.m index 440e6e9..24d57a1 100644 --- a/snpm_pp.m +++ b/snpm_pp.m @@ -529,7 +529,7 @@ function snpm_pp(CWD,varargin) end % No extent thresholding when voxelwise threshold is requested nidm_inference('nidm_ExtentThreshold/prov:type') = 'obo_statistic'; - nidm_inference('nidm_ExtentThreshold/prov:value') = 0; + nidm_inference('nidm_ExtentThreshold/nidm_clusterSizeInVoxels') = 0; else % Cluster-wise inference if exist(fullfile(CWD,'SnPM_ST.mat'))~=2 & exist(fullfile(CWD,'STCS.mat'))~=2 @@ -604,9 +604,19 @@ function snpm_pp(CWD,varargin) end ST_Ut = CFth; else % Threshold *was* set in snpm_ui. + +% TODO check why job.Thr.Clus.ClusSize.CFth is not stored in CFth if ~isnan(job.Thr.Clus.ClusSize.CFth) error('SnPM:InvalidClusterFormingThresh', sprintf('ERROR: Cluster-forming threshold of T=%0.2f was already set during analysis configuration; hence, in results, cluster-forming threshold must be left as "NaN".\nRe-run results with cluster-forming threshold set to NaN.\n',ST_Ut)) end + if (job.Thr.Clus.ClusSize.CFth<1) + nidm_inference('nidm_HeightThreshold/prov:type') = 'nidm_PValueUncorrected'; + nidm_inference('nidm_HeightThreshold/prov:value') = job.Thr.Clus.ClusSize.CFth; + else + nidm_inference('nidm_HeightThreshold/prov:type') = 'obo_statistic'; + nidm_inference('nidm_HeightThreshold/prov:value') = job.Thr.Clus.ClusSize.CFth; + end + end u=ST_Ut; % Flag use of a statistic-value threshold % Inference details... @@ -807,6 +817,10 @@ function snpm_pp(CWD,varargin) %-Calculate distribution of Maximum Suprathreshold Cluster size %-Calculate critical Suprathreshold Cluster Size %======================================================================= + +% TODO: check this is always true +nidm_json('nidm_ClusterDefinitionCriteria/nidm_hasConnectivityCriterion') = 'nidm_voxel18connected'; + if bSpatEx fprintf('Working on spatial extent...\n'); %-Compute suprathreshold voxels - check there are some @@ -932,9 +946,6 @@ function snpm_pp(CWD,varargin) %-Save perm 1 stats for use later - [X;Y;Z;T;perm;STCno] tmp = spm_clusters(Locs_vox(1:3,:)); - nidm_inferences = nidm_json('Inferences'); - nidm_inferences('nidm_ClusterDefinitionCriteria/nidm_hasConnectivityCriterion') = 'nidm_voxel18connected'; - nidm_json('Inferences') = nidm_inferences; if isPos==1 STCstats_Pos = [ SnPM_ST(:,tQ); tmp]; if bNeg==0 @@ -1393,17 +1404,17 @@ function snpm_pp(CWD,varargin) 'nidm_Peak/nidm_qValueFDR', ... 'nidm_Peak/prov:value', ... 'nidm_Peak/nidm_pValueUncorrected', ... + 'nidm_Peak/nidm_equivalentZStatistic',... 'nidm_Coordinate/nidm_coordinateVector', ... }, ... - {Pt(i), Pfdr(i), STC_SnPMt(i), Pu(i), STC_XYZ(1:3,i)}); + {Pt(i), Pfdr(i), STC_SnPMt(i), Pu(i), norminv(1-Pu(i)), STC_XYZ(1:3,i)}); y = y -1; %-Print up to 3 secondary maxima (>8mm apart) %------------------------------------------------------------------- - nidm_inferences = nidm_json('Inferences'); - nidm_inferences('nidm_PeakDefinitionCriteria/nidm_minDistanceBetweenPeaks') = 8; - nidm_json('Inferences') = nidm_inferences; + nidm_json('nidm_PeakDefinitionCriteria/nidm_minDistanceBetweenPeaks') = 8; + nidm_json('nidm_PeakDefinitionCriteria/nidm_maxNumberOfPeaksPerCluster') = 3; [null, k] = sort(-STC_SnPMt(j)); % Sort on t value D = i; @@ -1432,9 +1443,10 @@ function snpm_pp(CWD,varargin) 'nidm_Peak/nidm_qValueFDR', ... 'nidm_Peak/prov:value', ... 'nidm_Peak/nidm_pValueUncorrected', ... + 'nidm_Peak/nidm_equivalentZStatistic',... 'nidm_Coordinate/nidm_coordinateVector', ... }, ... - {Pt(d), Pfdr(d), STC_SnPMt(d), Pu(d), STC_XYZ(1:3,d)}); + {Pt(d), Pfdr(d), STC_SnPMt(d), Pu(d), norminv(1-Pu(d)), STC_XYZ(1:3,d)}); end @@ -1445,12 +1457,8 @@ function snpm_pp(CWD,varargin) r = r + 1; % Next region end - nidm_inferences = nidm_json('Inferences'); - nidm_inference = containers.Map(); nidm_inference('Clusters') = nidm_clusters; - nidm_inferences(contrast_id) = nidm_inference; - nidm_json('Inferences') = nidm_inferences; - + clear i j k D d r %-Footnote with SnPM parameters @@ -1496,7 +1504,14 @@ function snpm_pp(CWD,varargin) y=y-0.8; text(0,y,sprintf('Design: %s',sDesign),'FontSize',8); y=y-0.8; - text(0,y,sprintf('Search vol: %d cmm, %d voxels',S*abs(prod(VOX)),S), 'FontSize',8) + search_vol_cmm = S*abs(prod(VOX)); + search_vol_vox = S; + text(0,y,sprintf('Search vol: %d cmm, %d voxels',search_vol_cmm,search_vol_vox), 'FontSize',8) + nidm_inference('nidm_SearchSpaceMaskMap/nidm_searchVolumeInVoxels') = search_vol_vox; + +% TODO convert back to units + nidm_inference('nidm_SearchSpaceMaskMap/nidm_searchVolumeInUnits') = search_vol_cmm; + y=y-0.8; text(0.7,y,sprintf('Voxel size: [%5.2f, %5.2f, %5.2f] mm',abs(VOX)), ... 'FontSize', 8) @@ -1571,16 +1586,36 @@ function snpm_pp(CWD,varargin) Vs = sf_close_vol(Vs); clear t + % TODO SearchSpaceMaskMap can be different from analysis mask + nidm_inference('nidm_SearchSpaceMaskMap/prov:atLocation') = 'mask.img'; + + % TODO: always stationary?? + nidm_inference('nidm_SearchSpaceMaskMap/nidm_randomFieldStationarity') = true; - nidm_inferences = nidm_json('Inferences'); - nidm_inference = nidm_inferences(contrast_id); nidm_inference('nidm_ExcursionSetMap/prov:atLocation') = Fname; + + nidm_inferences = containers.Map(); nidm_inferences(contrast_id) = nidm_inference; nidm_json('Inferences') = nidm_inferences; + + nidm_json('nidm_NIDMResultsExporter/nidm_softwareVersion') = snpm('ver'); + nidm_json('nidm_NeuroimagingAnalysisSoftware/nidm_softwareVersion') = snpm('ver'); + + % TODO: are other units possible in SnPM?? + nidm_json('nidm_CoordinateSpace/units') = {'mm', 'mm', 'mm'}; + + % TODO: these should be filled in by the user --- + nidm_json('nidm_CoordinateSpace/nidm_inWorldCoordinateSystem') = 'nidm_StandardizedCoordinateSystem'; + nidm_json('nlx_Imaginginstrument/prov:type') = 'nlx_Magneticresonanceimagingscanner'; + % --- + + nidm_json('nidm_NeuroimagingAnalysisSoftware/prov:type') = 'src_SnPM'; + nidm_json('nidm_NeuroimagingAnalysisSoftware/prov:label') = 'SnPM'; % TODO: This temp file should only be produced if NIDM export is requested jsonwrite('snpm_nidm_thresh.json', nidm_json, ... struct('indent',' ', 'escape', false)); + spm_nidmresults(nidm_json, CWD) end %-Reset Interactive Window From 3352874cb295196d6b542d0a8efe71d9d77975f6 Mon Sep 17 00:00:00 2001 From: nicholst Date: Mon, 24 Apr 2017 14:34:51 +0100 Subject: [PATCH 17/29] Changes needed to write the grand mean --- snpm_cp.m | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/snpm_cp.m b/snpm_cp.m index e322346..d194316 100644 --- a/snpm_cp.m +++ b/snpm_cp.m @@ -260,7 +260,7 @@ function snpm_cp(CWD) %-Delete files from previous analyses, if they exist %----------------------------------------------------------------------- -files = { '^ResMS\..{3}$','^beta_.{4}\..{3}', '^lP_.{4}\..{3}',... +files = { '^ResMS\..{3}$','^beta_.{4}\..{3}','^GrandMean', '^lP_.{4}\..{3}',... '^lP.{1}\..{3}','^snpm.{2}\..{3}','^snpm.{1}\..{3}'}; for i=1:length(files) @@ -365,6 +365,7 @@ function snpm_cp(CWD) % %-Initialize image structures. % +Vgmean=snpm_clone_vol(Vt,'GrandMean','GrandMean'); for ii=1:p fname= sprintf('beta_%04d.img',ii); descrip=sprintf('beta_%04d hats',ii); @@ -442,6 +443,7 @@ function snpm_cp(CWD) %-Initialize the image data for this slice/volume %--------------------------------------------------------------------- + gmean_image=repmat(NaN,1,WorkDim); BETA_image=repmat(NaN,p,WorkDim); ResSS_image=repmat(NaN,1,WorkDim); if STAT=='T' @@ -531,6 +533,7 @@ function snpm_cp(CWD) % Use pseudo inverse rather than BETA=inv(D'*D)*D'*X for % D = DesMtx, to allow for non-unique designs. See matlab help. %----------------------------------------------------------------- + gmean = mean(X); BETA = pinv([H C B G])*X; ResSS = sum((X - [H C B G]*BETA).^2); @@ -604,6 +607,7 @@ function snpm_cp(CWD) % %- New! Write out data images. %- Input image data. + gmean_image(:,Q)=gmean; BETA_image(:,Q)=BETA; ResSS_image(:,Q)=ResSS; if STAT=='T' @@ -635,6 +639,9 @@ function snpm_cp(CWD) %-The image of the volume or the slice should be written out no matter length(Q)=1 %or 0. if bVolm + gmean_vol = reshape(gmean_image,DIM(1),DIM(2),DIM(3)); + spm_write_vol(gmean,gmean_vol); + for ii=1:p BETA_vol=reshape(BETA_image(ii,:),DIM(1),DIM(2),DIM(3)); spm_write_vol(Vbeta(ii),BETA_vol); @@ -661,6 +668,10 @@ function snpm_cp(CWD) end else + + gmean_plate=reshape(gmean_image, DIM(1), DIM(2)); + spm_write_plane(Vgmean,gmean_plate,i); + for ii=1:p BETA_plate=reshape(BETA_image(ii,:), DIM(1), DIM(2)); spm_write_plane(Vbeta(ii),BETA_plate,i); From cb5f43cb1baf984633c2aec496fdafc38c7205d1 Mon Sep 17 00:00:00 2001 From: cmaumet Date: Thu, 27 Apr 2017 17:12:00 +0100 Subject: [PATCH 18/29] Add all files to delete list --- snpm_cp.m | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/snpm_cp.m b/snpm_cp.m index 2e51ccd..5eabb43 100644 --- a/snpm_cp.m +++ b/snpm_cp.m @@ -311,7 +311,10 @@ function snpm_cp(CWD) %-Delete files from previous analyses, if they exist %----------------------------------------------------------------------- -files = { '^ResMS\..{3}$','^beta_.{4}\..{3}','^GrandMean', '^lP_.{4}\..{3}',... +files = { '^ResMS\..{3}$','^beta_.{4}\..{3}','^GrandMean\..{3}', ... + '^mask\..{3}', '^conse\..{3}',... + '^con.{2}\..{3}','^con.{1}\..{3}',... + '^lP_.{4}\..{3}',... '^lP.{1}\..{3}','^snpm.{2}\..{3}','^snpm.{1}\..{3}'}; for i=1:length(files) @@ -433,7 +436,7 @@ function snpm_cp(CWD) % %-Initialize image structures. % -Vgmean=snpm_clone_vol(Vt,'GrandMean','GrandMean'); +Vgmean=snpm_clone_vol(Vt,'GrandMean.img','GrandMean'); for ii=1:p fname{ii}= sprintf('beta_%04d.img',ii); descrip{ii}=sprintf('beta_%04d hats',ii); From 4d1033531db0ae97b3ee7d46f8a35af7add2ff0c Mon Sep 17 00:00:00 2001 From: cmaumet Date: Thu, 27 Apr 2017 17:17:02 +0100 Subject: [PATCH 19/29] Compute grand mean and all images slice-by-slice --- snpm_cp.m | 36 ++++++++++++++++++++++-------------- 1 file changed, 22 insertions(+), 14 deletions(-) diff --git a/snpm_cp.m b/snpm_cp.m index 5eabb43..868847e 100644 --- a/snpm_cp.m +++ b/snpm_cp.m @@ -436,7 +436,6 @@ function snpm_cp(CWD) % %-Initialize image structures. % -Vgmean=snpm_clone_vol(Vt,'GrandMean.img','GrandMean'); for ii=1:p fname{ii}= sprintf('beta_%04d.img',ii); descrip{ii}=sprintf('beta_%04d hats',ii); @@ -465,11 +464,9 @@ function snpm_cp(CWD) Vmask=spm_create_vol(Vmask); nidm_json('nidm_MaskMap/prov:atLocation') = 'mask.img'; -% TODO: grand mean needs to be computed -warning('Grand mean is not computed') -Vgm=snpm_clone_vol(Vt,'gm.img',str); -Vgm=spm_create_vol(Vgm); -nidm_json('nidm_GrandMeanMap/prov:atLocation') = 'gm.img'; +Vgmean=snpm_clone_vol(Vt,'GrandMean.img','GrandMean'); +Vgmean=spm_create_vol(Vgmean); +nidm_json('nidm_GrandMeanMap/prov:atLocation') = 'GrandMean.img'; if STAT=='T' VT_pos=snpm_clone_vol(Vt,'snpmT+.img',[str,' (+ve)']); @@ -773,9 +770,6 @@ function snpm_cp(CWD) %-The image of the volume or the slice should be written out no matter length(Q)=1 %or 0. if bVolm - gmean_vol = reshape(gmean_image,DIM(1),DIM(2),DIM(3)); - spm_write_vol(gmean,gmean_vol); - for ii=1:p BETA_vol=reshape(BETA_image(ii,:),DIM(1),DIM(2),DIM(3)); spm_write_vol(Vbeta(ii),BETA_vol); @@ -787,9 +781,10 @@ function snpm_cp(CWD) % Analysis mask mask_vol=reshape(Q,DIM(1),DIM(2),DIM(3)); spm_write_vol(Vmask,mask_vol); - - % TODO: replace mask volume by computed grand mean - spm_write_vol(Vgm,mask_vol); + + % Grand mean + gmean_vol = reshape(gmean_image,DIM(1),DIM(2),DIM(3)); + spm_write_vol(Vgmean,gmean_vol); if STAT=='T' T_pos_vol=reshape(T_pos_image,DIM(1),DIM(2),DIM(3)); @@ -817,8 +812,12 @@ function snpm_cp(CWD) spm_write_vol(VlwP, lwP_vol); end - else - + else + % Analysis mask + mask_plate=reshape(Q,DIM(1),DIM(2)); + spm_write_plane(Vmask,mask_plate); + + % Grand mean gmean_plate=reshape(gmean_image, DIM(1), DIM(2)); spm_write_plane(Vgmean,gmean_plate,i); @@ -836,6 +835,15 @@ function snpm_cp(CWD) T_neg_plate=reshape(T_neg_image, DIM(1), DIM(2)); spm_write_plane(VT_neg,T_neg_plate,i); + + CON_pos_plate=reshape(CON_pos_image,DIM(1),DIM(2)); + spm_write_plane(VCON_pos,CON_pos_plate); + + CON_neg_plate=reshape(CON_neg_image,DIM(1),DIM(2)); + spm_write_plane(VCON_neg,CON_neg_plate); + + CONSE_plate=reshape(CONSE_image,DIM(1),DIM(2)); + spm_write_plane(VCONSE,CONSE_plate); elseif STAT=='F' F_plate=reshape(F_image, DIM(1), DIM(2)); From b61ff14250805bac21c15723ea5ac9bb16f99750 Mon Sep 17 00:00:00 2001 From: cmaumet Date: Fri, 28 Apr 2017 12:55:05 +0100 Subject: [PATCH 20/29] Add nidm_ContrastExplainedMeanSquareMap --- snpm_cp.m | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/snpm_cp.m b/snpm_cp.m index 868847e..2b6fb5f 100644 --- a/snpm_cp.m +++ b/snpm_cp.m @@ -257,6 +257,7 @@ function snpm_cp(CWD) con_name = 'Positive'; nidm_json('Contrasts') = containers.Map({con_name}, {contrastMap}); + con_neg_name = ''; else contrast_pos = containers.Map(); contrast_neg = containers.Map(); % There is no deep copy of @@ -513,6 +514,10 @@ function snpm_cp(CWD) nidm_contrast('nidm_StatisticMap/prov:atLocation') = 'snpmF.img'; nidm_contrasts(con_name) = nidm_contrast; nidm_json('Contrasts') = nidm_contrasts; + + VFnum=snpm_clone_vol(Vt,'snpmFnum.img',str); + VFnum=spm_create_vol(VFnum); + nidm_contrast('nidm_ContrastExplainedMeanSquareMap/prov:atLocation') = 'snpmFnum.img'; end VlP_pos=snpm_clone_vol(Vt, 'lP+.img', '-log10(uncor. non-para. P, +ve)'); @@ -692,6 +697,7 @@ function snpm_cp(CWD) %-Compute t-statistics for specified compounds of parameters %----------------------------------------------------------- T = zeros(1,size(BETA,2)); + Fnum = zeros(1,size(BETA,2)); CON = zeros(1,size(BETA,2)); Co = CONT; if STAT=='T' @@ -704,6 +710,8 @@ function snpm_cp(CWD) pX = pinv([H C B G]); T(1,:) = (sum(((Co*BETA)'*inv(Co*pinv([H C B G]'*[H C B G])*Co'))' .* ... (Co*BETA),1)/rank(Co)) ./ (ResSS/df); + Fnum(1,:) = (sum(((Co*BETA)'*inv(Co*pinv([H C B G]'*[H C B G])*Co'))' .* ... + (Co*BETA),1)/rank(Co)); end %-Save Max T statistic @@ -746,6 +754,7 @@ function snpm_cp(CWD) CONSE_image(:,Q)=CONSE; elseif STAT=='F' F_image(:,Q)=T; + Fnum_image(:,Q)=Fnum; end if bVarAlph @@ -805,6 +814,9 @@ function snpm_cp(CWD) elseif STAT=='F' F_vol=reshape(F_image,DIM(1),DIM(2),DIM(3)); spm_write_vol(VF,F_vol); + + Fnum_vol=reshape(Fnum_image,DIM(1),DIM(2),DIM(3)); + spm_write_vol(VFnum,Fnum_vol); end if bVarAlph @@ -848,6 +860,9 @@ function snpm_cp(CWD) elseif STAT=='F' F_plate=reshape(F_image, DIM(1), DIM(2)); spm_write_plane(VF,F_plate,i); + + Fnum_plate=reshape(Fnum_image, DIM(1), DIM(2)); + spm_write_plane(VFnum,Fnum_plate,i); end if bVarAlph From 729a52356c0954e834054793b8caef778eb4f0c4 Mon Sep 17 00:00:00 2001 From: cmaumet Date: Fri, 28 Apr 2017 13:15:06 +0100 Subject: [PATCH 21/29] Save max number of permutations in all modules --- snpm_cp.m | 2 +- snpm_ui.m | 13 ++++++++++++- 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/snpm_cp.m b/snpm_cp.m index 2b6fb5f..becd723 100644 --- a/snpm_cp.m +++ b/snpm_cp.m @@ -350,7 +350,7 @@ function snpm_cp(CWD) nidm_json('nidm_NonParametricNullDistribution/nidm_numberOfPermutations') = nPerm; nidm_json('nidm_NonParametricNullDistribution/nidm_hasResamplingScheme') = 'nidm_Permutation'; nidm_json('nidm_NonParametricNullDistribution/nidm_hasApproximationMethod') = 'nidm_MonteCarlo'; -nidm_json('nidm_NonParametricNullDistribution/nidm_maximumNumberOfPermutations') = nPiCond_mx; +nidm_json('nidm_NonParametricNullDistribution/nidm_maximumNumberOfPermutations') = nPerm_max; con_names = nidm_json('Contrasts').keys; for i = 1:numel(con_names) diff --git a/snpm_ui.m b/snpm_ui.m index 5f1484a..21732f2 100644 --- a/snpm_ui.m +++ b/snpm_ui.m @@ -303,6 +303,17 @@ function snpm_ui(varargin) %-Total #observations nScan = size(P,1); +%-Max number of possible permutations +if exist('nPiCond_mx', 'var') + nPerm_max = nPiCond_mx; +elseif exist('nPiStud_mx', 'var') + nPerm_max = nPiStud_mx; +elseif exist('nPiSubj_mx', 'var') + nPerm_max = nPiSubj_mx; +else + error('snpm:MissingMaxPerm', 'Maximum number of permutations not found') +end + %-Get general analysis & data parameters %======================================================================= @@ -587,7 +598,7 @@ function snpm_ui(varargin) s_SnPMcfg_save = ['s_SnPMcfg_save H C B G HCBGnames P PiCond ',... 'sPiCond bhPerms sHCform iGloNorm sGloNorm GM rg GX GMscale CONT ',... 'THRESH MASK ImMASK TH bVarSm vFWHM sVarSm bVolm bST sDesFile sDesign ',... - 'V pU_ST_Ut df1 nidm_json nPiCond_mx nidm_json ', ... + 'V pU_ST_Ut df1 nidm_json nPerm_max nidm_json ', ... 'sDesSave ',sDesSave]; eval(['save SnPMcfg ',s_SnPMcfg_save]) From 78dbe416f60d04c63207a1eaa095af218fe94e22 Mon Sep 17 00:00:00 2001 From: cmaumet Date: Fri, 7 Jul 2017 09:02:01 +0100 Subject: [PATCH 22/29] containers to struct --- snpm_cp.m | 124 +++++++++++++++++++----------------------------- snpm_pp.m | 138 +++++++++++++++++++++++------------------------------- snpm_ui.m | 16 +++---- 3 files changed, 116 insertions(+), 162 deletions(-) diff --git a/snpm_cp.m b/snpm_cp.m index becd723..4e80a76 100644 --- a/snpm_cp.m +++ b/snpm_cp.m @@ -217,22 +217,22 @@ function snpm_cp(CWD) if strcmp(sDesFile, 'snpm_pi_OneSampT') || ... strcmp(sDesFile, 'snpm_pi_ANOVAwithinS') % Sign flipping - nidm_json('nidm_ErrorModel/nidm_hasErrorDistribution') = {'obo_nonparametricdistribution', 'obo_symmetricdistribution'}; - nidm_json('nidm_ErrorModel/nidm_errorVarianceHomogeneous') = false; - nidm_json('nidm_ErrorModel/nidm_varianceMapWiseDependence') = 'nidm_IndependentParameter'; - nidm_json('nidm_ErrorModel/nidm_hasErrorDependence') = 'nidm_IndependentError'; + nidm.ErrorModel_hasErrorDistribution = {'obo_nonparametricdistribution', 'obo_symmetricdistribution'}; + nidm.ErrorModel_errorVarianceHomogeneous = false; + nidm.ErrorModel_varianceMapWiseDependence = 'nidm_IndependentParameter'; + nidm.ErrorModel_hasErrorDependence = 'nidm_IndependentError'; else % Permutation - nidm_json('nidm_ErrorModel/nidm_hasErrorDistribution') = 'obo_nonparametricdistribution'; - nidm_json('nidm_ErrorModel/nidm_errorVarianceHomogeneous') = true; - nidm_json('nidm_ErrorModel/nidm_varianceMapWiseDependence') = 'nidm_IndependentParameter'; + nidm.ErrorModel_hasErrorDistribution = 'obo_nonparametricdistribution'; + nidm.ErrorModel_errorVarianceHomogeneous = true; + nidm.ErrorModel_varianceMapWiseDependence = 'nidm_IndependentParameter'; % TODO: the 'obo_exchangeable' term is not yet in STATO - nidm_json('nidm_ErrorModel/nidm_hasErrorDependence') = 'obo_exchangeable'; - nidm_json('nidm_ErrorModel/nidm_dependenceMapWiseDependence') = 'nidm_IndependentParameter'; + nidm.ErrorModel_hasErrorDependence = 'obo_exchangeable'; + nidm.ErrorModel_dependenceMapWiseDependence = 'nidm_IndependentParameter'; end % TODO: check this is correct -nidm_json('nidm_ModelParameterEstimation/nidm_withEstimationMethod') = 'obo_ordinaryleastsquaresestimation'; +nidm.ModelParameterEstimation_withEstimationMethod = 'obo_ordinaryleastsquaresestimation'; if isempty([H C]) error('SnPM:NoModel', 'No model specified; [H C] empty'); @@ -244,44 +244,38 @@ function snpm_cp(CWD) error('SnPM:InvalidContrast','Contrast problem; wrong number of columns'); end -if size(CONT,1) > 1 - contrastMap = containers.Map(); - contrastMap('obo_contrastweightmatrix/prov:value') = CONT; +if size(CONT,1) > 1 % F-contrast + nidm.Contrasts(1).contrastweightmatrix_value = CONT; warning('SnPM:FContrast', ... 'F contrast! F statistic images are being created.'); STAT = 'F'; - contrastMap('nidm_StatisticMap/nidm_statisticType') = 'obo_Fstatistic'; + nidm.Contrasts(1).StatisticMap_statisticType = 'obo_Fstatistic'; + if (CONT(1,:) == -CONT(2,:)) CONT = CONT(1,:); end con_name = 'Positive'; - - nidm_json('Contrasts') = containers.Map({con_name}, {contrastMap}); + nidm.Contrasts(1).StatisticMap_contrastName = con_name; con_neg_name = ''; else - contrast_pos = containers.Map(); - contrast_neg = containers.Map(); % There is no deep copy of - % containers.Map so we need a separate object from the start con_name = 'Positive'; + nidm.Contrasts(1).StatisticMap_contrastName = con_name; + con_neg_name = 'Negative'; + nidm.Contrasts(2).StatisticMap_contrastName = con_neg_name; + STAT = 'T'; if bVarSm - contrast_pos('nidm_StatisticMap/nidm_statisticType') = 'nidm_smoothedtstatistic'; - contrast_neg('nidm_StatisticMap/nidm_statisticType') = 'nidm_smoothedtstatistic'; + nidm.Contrasts(1).StatisticMap_statisticType = 'nidm_smoothedtstatistic'; + nidm.Contrasts(2).StatisticMap_statisticType = 'nidm_smoothedtstatistic'; else - contrast_pos('nidm_StatisticMap/nidm_statisticType') = 'obo_tstatistic'; - contrast_neg('nidm_StatisticMap/nidm_statisticType') = 'obo_tstatistic'; + nidm.Contrasts(1).StatisticMap_statisticType = 'obo_tstatistic'; + nidm.Contrasts(2).StatisticMap_statisticType = 'obo_tstatistic'; end - contrast_pos('nidm_ContrastMap/nidm_contrastName') = ['Positive T-Contrast: [' mat2str(CONT) ']']; - contrast_pos('obo_contrastweightmatrix/prov:value') = CONT; - contrast_neg('nidm_ContrastMap/nidm_contrastName') = ['Negative T-Contrast: [' mat2str(-CONT) ']']; - contrast_neg('obo_contrastweightmatrix/prov:value') = -CONT; - -% nidm_json('Contrasts').(con_name) = contrast_pos; -% nidm_json('Contrasts').(con_neg_name) = contrast_neg; - - nidm_json('Contrasts') = containers.Map({con_name, con_neg_name}, ... - {contrast_pos, contrast_neg}); + nidm.Contrasts(1).ContrastMap_contrastName = ['Positive T-Contrast: [' mat2str(CONT) ']']; + nidm.Contrasts(1).contrastweightmatrix_value = CONT; + nidm.Contrasts(2).ContrastMap_contrastName = ['Negative T-Contrast: [' mat2str(-CONT) ']']; + nidm.Contrasts(2).contrastweightmatrix_value = -CONT; end if rank(CONT) u -tol) / nPermReal; - nidm_inference('nidm_HeightThreshold/prov:type') = 'obo_statistic'; - nidm_inference('nidm_HeightThreshold/prov:value') = u; + nidm.Inferences(1).HeightThreshold_type = 'obo_statistic'; + nidm.Inferences(1).HeightThreshold_value = u; case 'FDRth' alph_FDR = BoundCheck(job.Thr.Vox.VoxSig.FDRth,[0 1],'Invalid FDR level'); alpha_ucp = snpm_uc_FDR(alph_FDR,[],'P',[],sSnPMucp'); - nidm_inference('nidm_HeightThreshold/prov:type') = 'obo_qvalue'; - nidm_inference('nidm_HeightThreshold/prov:value') = alph_FDR; + nidm.Inferences(1).HeightThreshold_type = 'obo_qvalue'; + nidm.Inferences(1).HeightThreshold_value = alph_FDR; case 'FWEth' alph_FWE = BoundCheck(job.Thr.Vox.VoxSig.FWEth,[0 1],'Invalid FWE level'); iFWE = ceil((1-alph_FWE)*nPermReal); @@ -522,14 +521,14 @@ function snpm_pp(CWD,varargin) C_MaxT = 0; end u = C_MaxT; - nidm_inference('nidm_HeightThreshold/prov:type') = 'obo_FWERadjustedpvalue'; - nidm_inference('nidm_HeightThreshold/prov:value') = alph_FWE; + nidm.Inferences(1).HeightThreshold_type = 'obo_FWERadjustedpvalue'; + nidm.Inferences(1).HeightThreshold_value = alph_FWE; otherwise error('Unknown threshold') end % No extent thresholding when voxelwise threshold is requested - nidm_inference('nidm_ExtentThreshold/prov:type') = 'obo_statistic'; - nidm_inference('nidm_ExtentThreshold/nidm_clusterSizeInVoxels') = 0; + nidm.Inferences(1).ExtentThreshold_type = 'obo_statistic'; + nidm.Inferences(1).ExtentThreshold_clusterSizeInVoxels = 0; else % Cluster-wise inference if exist(fullfile(CWD,'SnPM_ST.mat'))~=2 & exist(fullfile(CWD,'STCS.mat'))~=2 @@ -545,11 +544,11 @@ function snpm_pp(CWD,varargin) ST_Ut_0 = ST_Ut; CFth=job.Thr.Clus.ClusSize.CFth; if (CFth<1) - nidm_inference('nidm_HeightThreshold/prov:type') = 'nidm_PValueUncorrected'; - nidm_inference('nidm_HeightThreshold/prov:value') = CFth; + nidm.Inferences(1).HeightThreshold_type = 'nidm_PValueUncorrected'; + nidm.Inferences(1).HeightThreshold_value = CFth; else - nidm_inference('nidm_HeightThreshold/prov:type') = 'obo_statistic'; - nidm_inference('nidm_HeightThreshold/prov:value') = CFth; + nidm.Inferences(1).HeightThreshold_type = 'obo_statistic'; + nidm.Inferences(1).HeightThreshold_value = CFth; end if (CFth<=0) @@ -610,11 +609,11 @@ function snpm_pp(CWD,varargin) error('SnPM:InvalidClusterFormingThresh', sprintf('ERROR: Cluster-forming threshold of T=%0.2f was already set during analysis configuration; hence, in results, cluster-forming threshold must be left as "NaN".\nRe-run results with cluster-forming threshold set to NaN.\n',ST_Ut)) end if (job.Thr.Clus.ClusSize.CFth<1) - nidm_inference('nidm_HeightThreshold/prov:type') = 'nidm_PValueUncorrected'; - nidm_inference('nidm_HeightThreshold/prov:value') = job.Thr.Clus.ClusSize.CFth; + nidm.Inferences(1).HeightThreshold_type = 'nidm_PValueUncorrected'; + nidm.Inferences(1).HeightThreshold_value = job.Thr.Clus.ClusSize.CFth; else - nidm_inference('nidm_HeightThreshold/prov:type') = 'obo_statistic'; - nidm_inference('nidm_HeightThreshold/prov:value') = job.Thr.Clus.ClusSize.CFth; + nidm.Inferences(1).HeightThreshold_type = 'obo_statistic'; + nidm.Inferences(1).HeightThreshold_value = job.Thr.Clus.ClusSize.CFth; end end @@ -624,17 +623,17 @@ function snpm_pp(CWD,varargin) switch tmp{1} case 'Cth' C_STCS = job.Thr.Clus.ClusSize.ClusSig.Cth; - nidm_inference('nidm_ExtentThreshold/prov:type') = 'obo_statistic'; - nidm_inference('nidm_ExtentThreshold/prov:value') = C_STCS; + nidm.Inferences(1).ExtentThreshold_type = 'obo_statistic'; + nidm.Inferences(1).ExtentThreshold_value = C_STCS; case 'PthC' alpha_ucp = BoundCheck(job.Thr.Clus.ClusSize.ClusSig.PthC,[0 1],'Invalid uncorrected P(k)'); - nidm_inference('nidm_ExtentThreshold/prov:type') = 'nidm_PValueUncorrected'; - nidm_inference('nidm_ExtentThreshold/prov:value') = alpha_ucp; + nidm.Inferences(1).ExtentThreshold_type = 'nidm_PValueUncorrected'; + nidm.Inferences(1).ExtentThreshold_value = alpha_ucp; case 'FWEthC' alph_FWE = BoundCheck(job.Thr.Clus.ClusSize.ClusSig.FWEthC,[0 1],'Invalid FWE level (cluster-level inference)'); iFWE = ceil((1-alph_FWE)*nPermReal); - nidm_inference('nidm_ExtentThreshold/prov:type') = 'obo_FWERadjustedpvalue'; - nidm_inference('nidm_ExtentThreshold/prov:value') = alph_FWE; + nidm.Inferences(1).ExtentThreshold_type = 'obo_FWERadjustedpvalue'; + nidm.Inferences(1).ExtentThreshold_value = alph_FWE; end end % END: Cluster-wise inference @@ -819,7 +818,7 @@ function snpm_pp(CWD,varargin) %======================================================================= % TODO: check this is always true -nidm_json('nidm_ClusterDefinitionCriteria/nidm_hasConnectivityCriterion') = 'nidm_voxel18connected'; +nidm.ClusterDefinitionCriteria_hasConnectivityCriterion = 'nidm_voxel18connected'; if bSpatEx fprintf('Working on spatial extent...\n'); @@ -1192,8 +1191,8 @@ function snpm_pp(CWD,varargin) end % Display only if *not* in command line mode -nidm=true; -if ~spm_get_defaults('cmdline') || nidm +nidm_export=true; +if ~spm_get_defaults('cmdline') || nidm_export %======================================================================= %-D I S P L A Y : Max report @@ -1362,11 +1361,8 @@ function snpm_pp(CWD,varargin) %----------------------------------------------------------------------- r = 1; bUsed = zeros(size(STC_SnPMt)); - nidm_clusters = containers.Map(); while max(STC_SnPMt.*(~bUsed)) & (y > 3) - nidm_cluster = containers.Map(); - nidm_peaks = containers.Map(); [null, i] = max(STC_SnPMt.*(~bUsed)); % Largest t value j = find(STC_r == STC_r(i)); % Maxima in same region @@ -1397,24 +1393,21 @@ function snpm_pp(CWD,varargin) text(tCol(9),y,sprintf(Fmtst{9},STC_XYZ(2,i)),'UserData',STC_XYZ(:,i),StrAttrB{:}) text(tCol(10),y,sprintf(Fmtst{10},STC_XYZ(3,i)),'UserData',STC_XYZ(:,i),StrAttrB{:}) - nidm_cluster('nidm_SupraThresholdCluster/nidm_clusterSizeInVoxels') = STC_N(i); - nidm_peaks('Peak_1') = containers.Map(... - { ... - 'nidm_Peak/nidm_pValueFWER', ... - 'nidm_Peak/nidm_qValueFDR', ... - 'nidm_Peak/prov:value', ... - 'nidm_Peak/nidm_pValueUncorrected', ... - 'nidm_Peak/nidm_equivalentZStatistic',... - 'nidm_Coordinate/nidm_coordinateVector', ... - }, ... - {Pt(i), Pfdr(i), STC_SnPMt(i), Pu(i), norminv(1-Pu(i)), STC_XYZ(1:3,i)}); + % TODO: this is overwritten by 1st peak i below + nidm.Inferences(1).Clusters(r).SupraThresholdCluster_clusterSizeInVoxels = STC_N(i); + nidm.Inferences(1).Clusters(r).Peaks(1).Peak_pValueFWER = Pt(i); + nidm.Inferences(1).Clusters(r).Peaks(1).Peak_pValueFDR = Pfdr(i); + nidm.Inferences(1).Clusters(r).Peaks(1).Peak_value = STC_SnPMt(i); + nidm.Inferences(1).Clusters(r).Peaks(1).Peak_pValueUncorrected = Pu(i); + nidm.Inferences(1).Clusters(r).Peaks(1).Peak_equivalentZStatistic = norminv(1-Pu(i)); + nidm.Inferences(1).Clusters(r).Peaks(1).Coordinate_coordinateVector = STC_XYZ(1:3,i); y = y -1; %-Print up to 3 secondary maxima (>8mm apart) %------------------------------------------------------------------- - nidm_json('nidm_PeakDefinitionCriteria/nidm_minDistanceBetweenPeaks') = 8; - nidm_json('nidm_PeakDefinitionCriteria/nidm_maxNumberOfPeaksPerCluster') = 3; + nidm.PeakDefinitionCriteria_minDistanceBetweenPeaks = 8; + nidm.PeakDefinitionCriteria_maxNumberOfPeaksPerCluster = 3; [null, k] = sort(-STC_SnPMt(j)); % Sort on t value D = i; @@ -1436,28 +1429,19 @@ function snpm_pp(CWD,varargin) y = y -1; end end + - nidm_peaks(['Peak_' num2str(i)]) = containers.Map(... - { ... - 'nidm_Peak/nidm_pValueFWER', ... - 'nidm_Peak/nidm_qValueFDR', ... - 'nidm_Peak/prov:value', ... - 'nidm_Peak/nidm_pValueUncorrected', ... - 'nidm_Peak/nidm_equivalentZStatistic',... - 'nidm_Coordinate/nidm_coordinateVector', ... - }, ... - {Pt(d), Pfdr(d), STC_SnPMt(d), Pu(d), norminv(1-Pu(d)), STC_XYZ(1:3,d)}); - + nidm.Inferences(1).Clusters(r).Peaks(i).Peak_pValueFWER = Pt(d); + nidm.Inferences(1).Clusters(r).Peaks(i).Peak_pValueFDR = Pfdr(d); + nidm.Inferences(1).Clusters(r).Peaks(i).Peak_value = STC_SnPMt(d); + nidm.Inferences(1).Clusters(r).Peaks(i).Peak_pValueUncorrected = Pu(d); + nidm.Inferences(1).Clusters(r).Peaks(i).Peak_equivalentZStatistic = norminv(1-Pu(d)); + nidm.Inferences(1).Clusters(r).Peaks(i).Coordinate_coordinateVector = STC_XYZ(1:3,d); end - - nidm_cluster('Peaks') = nidm_peaks; - nidm_clusters(['Cluster_' num2str(r)]) = nidm_cluster; - + bUsed(j) = (bUsed(j) | 1 ); %-Mark maxima as "used" r = r + 1; % Next region end - - nidm_inference('Clusters') = nidm_clusters; clear i j k D d r @@ -1507,10 +1491,10 @@ function snpm_pp(CWD,varargin) search_vol_cmm = S*abs(prod(VOX)); search_vol_vox = S; text(0,y,sprintf('Search vol: %d cmm, %d voxels',search_vol_cmm,search_vol_vox), 'FontSize',8) - nidm_inference('nidm_SearchSpaceMaskMap/nidm_searchVolumeInVoxels') = search_vol_vox; + nidm.SearchSpaceMaskMap_searchVolumeInVoxels = search_vol_vox; % TODO convert back to units - nidm_inference('nidm_SearchSpaceMaskMap/nidm_searchVolumeInUnits') = search_vol_cmm; + nidm.SearchSpaceMaskMap_searchVolumeInUnits = search_vol_cmm; y=y-0.8; text(0.7,y,sprintf('Voxel size: [%5.2f, %5.2f, %5.2f] mm',abs(VOX)), ... @@ -1540,7 +1524,7 @@ function snpm_pp(CWD,varargin) %- Image output? %======================================================================= %-Write out filtered SnPMt? -if WrtFlt || nidm +if WrtFlt || nidm_export Fname = WrtFltFn; @@ -1587,30 +1571,26 @@ function snpm_pp(CWD,varargin) clear t % TODO SearchSpaceMaskMap can be different from analysis mask - nidm_inference('nidm_SearchSpaceMaskMap/prov:atLocation') = 'mask.img'; + nidm.SearchSpaceMaskMap_atLocation = 'mask.img'; % TODO: always stationary?? - nidm_inference('nidm_SearchSpaceMaskMap/nidm_randomFieldStationarity') = true; - - nidm_inference('nidm_ExcursionSetMap/prov:atLocation') = Fname; + nidm.SearchSpaceMaskMap_randomFieldStationarity = true; - nidm_inferences = containers.Map(); - nidm_inferences(contrast_id) = nidm_inference; - nidm_json('Inferences') = nidm_inferences; + nidm.Inferences(1).ExcursionSetMap_atLocation = Fname; - nidm_json('nidm_NIDMResultsExporter/nidm_softwareVersion') = snpm('ver'); - nidm_json('nidm_NeuroimagingAnalysisSoftware/nidm_softwareVersion') = snpm('ver'); + nidm.NIDMResultsExporter_softwareVersion = snpm('ver'); + nidm.NeuroimagingAnalysisSoftware_softwareVersion = snpm('ver'); % TODO: are other units possible in SnPM?? - nidm_json('nidm_CoordinateSpace/units') = {'mm', 'mm', 'mm'}; + nidm.CoordinateSpace_units = {'mm', 'mm', 'mm'}; % TODO: these should be filled in by the user --- - nidm_json('nidm_CoordinateSpace/nidm_inWorldCoordinateSystem') = 'nidm_StandardizedCoordinateSystem'; - nidm_json('nlx_Imaginginstrument/prov:type') = 'nlx_Magneticresonanceimagingscanner'; + nidm.CoordinateSpace_inWorldCoordinateSystem = 'nidm_StandardizedCoordinateSystem'; + nidm.Imaginginstrument_type = 'nlx_Magneticresonanceimagingscanner'; % --- - nidm_json('nidm_NeuroimagingAnalysisSoftware/prov:type') = 'src_SnPM'; - nidm_json('nidm_NeuroimagingAnalysisSoftware/prov:label') = 'SnPM'; + nidm.NeuroimagingAnalysisSoftware_type = 'src_SnPM'; + nidm.NeuroimagingAnalysisSoftware_label = 'SnPM'; % TODO: This temp file should only be produced if NIDM export is requested jsonwrite('snpm_nidm_thresh.json', nidm_json, ... diff --git a/snpm_ui.m b/snpm_ui.m index 21732f2..ce38701 100644 --- a/snpm_ui.m +++ b/snpm_ui.m @@ -123,7 +123,7 @@ function snpm_ui(varargin) % V Memory mapping handles % MASK Filename of explicit mask image % ImMASK Implicit masking; 0=none; 1=zeros are equivalent to NaN -% nidm_json json structure storing minimal information required for a +% nidm json structure storing minimal information required for a % NIDM-Results export % % df degrees of freedom due to error @@ -193,7 +193,7 @@ function snpm_ui(varargin) cd(job.dir{1}) end -nidm_json = containers.Map(); +nidm = struct(); %-Definitions & Design parameters %======================================================================= @@ -383,7 +383,7 @@ function snpm_ui(varargin) end if (iGMsca==2) % CHANGED from 1 to 2 as should not ask for a value if grand mean scaling is not required. - nidm_json('nidm_Data/nidm_grandMeanScaling') = true; + nidm.Data_grandMeanScaling = true; if (iGloNorm==2) % Proportional scaling str = 'PropSca global mean to'; else @@ -397,10 +397,10 @@ function snpm_ui(varargin) case 'gmsca_no', GM = 50; end - nidm_json('nidm_Data/nidm_targetIntensity') = GM; + nidm.Data_targetIntensity = GM; elseif (iGMsca==1) % No grand mean scaling GM = 0; - nidm_json('nidm_Data/nidm_grandMeanScaling') = false; + nidm.Data_grandMeanScaling = false; end @@ -590,15 +590,15 @@ function snpm_ui(varargin) %-Construct full design matrix and name matrices for display %----------------------------------------------------------------------- [nHCBG,HCBGnames] = spm_DesMtx('Sca',H,Hnames,C,Cnames,B,Bnames,G,Gnames); -nidm_json('nidm_DesignMatrix/prov:value') = nHCBG; -nidm_json('nidm_DesignMatrix/nidm_regressorNames') = HCBGnames; +nidm.DesignMatrix_value = nHCBG; +nidm.DesignMatrix_regressorNames = HCBGnames; %-Setup is complete - save SnPMcfg Mat file %----------------------------------------------------------------------- s_SnPMcfg_save = ['s_SnPMcfg_save H C B G HCBGnames P PiCond ',... 'sPiCond bhPerms sHCform iGloNorm sGloNorm GM rg GX GMscale CONT ',... 'THRESH MASK ImMASK TH bVarSm vFWHM sVarSm bVolm bST sDesFile sDesign ',... - 'V pU_ST_Ut df1 nidm_json nPerm_max nidm_json ', ... + 'V pU_ST_Ut df1 nidm nPerm_max nidm_json ', ... 'sDesSave ',sDesSave]; eval(['save SnPMcfg ',s_SnPMcfg_save]) From 76a20a96be56caaf5d3c3a507cd90dc28539cc54 Mon Sep 17 00:00:00 2001 From: cmaumet Date: Fri, 7 Jul 2017 10:11:22 +0100 Subject: [PATCH 23/29] fixes to make all tests pass --- snpm_cp.m | 15 +++++++++------ snpm_pp.m | 20 ++++++++++---------- snpm_ui.m | 4 +++- 3 files changed, 22 insertions(+), 17 deletions(-) diff --git a/snpm_cp.m b/snpm_cp.m index 4e80a76..5cd6904 100644 --- a/snpm_cp.m +++ b/snpm_cp.m @@ -545,6 +545,8 @@ function snpm_cp(CWD) gmean_image=repmat(NaN,1,WorkDim); BETA_image=repmat(NaN,p,WorkDim); ResSS_image=repmat(NaN,1,WorkDim); + mask_image=repmat(NaN,1,WorkDim); + CONSE_image=repmat(NaN,1,WorkDim); if STAT=='T' T_pos_image=repmat(NaN,1,WorkDim); T_neg_image=repmat(NaN,1,WorkDim); @@ -718,6 +720,7 @@ function snpm_cp(CWD) gmean_image(:,Q)=gmean; BETA_image(:,Q)=BETA; ResSS_image(:,Q)=ResSS; + mask_image(:,Q)=1; if STAT=='T' T_pos_image(:,Q)=T; T_neg_image(:,Q)=-T; @@ -762,7 +765,7 @@ function snpm_cp(CWD) spm_write_vol(VResMS, ResSS_vol); % Analysis mask - mask_vol=reshape(Q,DIM(1),DIM(2),DIM(3)); + mask_vol=reshape(mask_image,DIM(1),DIM(2),DIM(3)); spm_write_vol(Vmask,mask_vol); % Grand mean @@ -800,8 +803,8 @@ function snpm_cp(CWD) else % Analysis mask - mask_plate=reshape(Q,DIM(1),DIM(2)); - spm_write_plane(Vmask,mask_plate); + mask_plate=reshape(mask_image,DIM(1),DIM(2)); + spm_write_plane(Vmask,mask_plate,i); % Grand mean gmean_plate=reshape(gmean_image, DIM(1), DIM(2)); @@ -823,13 +826,13 @@ function snpm_cp(CWD) spm_write_plane(VT_neg,T_neg_plate,i); CON_pos_plate=reshape(CON_pos_image,DIM(1),DIM(2)); - spm_write_plane(VCON_pos,CON_pos_plate); + spm_write_plane(VCON_pos,CON_pos_plate,i); CON_neg_plate=reshape(CON_neg_image,DIM(1),DIM(2)); - spm_write_plane(VCON_neg,CON_neg_plate); + spm_write_plane(VCON_neg,CON_neg_plate,i); CONSE_plate=reshape(CONSE_image,DIM(1),DIM(2)); - spm_write_plane(VCONSE,CONSE_plate); + spm_write_plane(VCONSE,CONSE_plate,i); elseif STAT=='F' F_plate=reshape(F_image, DIM(1), DIM(2)); diff --git a/snpm_pp.m b/snpm_pp.m index b274d3e..ac71f2b 100644 --- a/snpm_pp.m +++ b/snpm_pp.m @@ -390,9 +390,9 @@ function snpm_pp(CWD,varargin) % Find corresponding contrast name in nidm json structure if bNeg == 0; - contrast_id = con_name; + nidm.Inferences(1).StatisticMap_contrastName = {con_name}; else - contrast_id = con_neg_name; + nidm.Inferences(1).StatisticMap_contrastName = {con_neg_name}; end %-Take MaxT for increases or decreases according to bNeg @@ -624,7 +624,7 @@ function snpm_pp(CWD,varargin) case 'Cth' C_STCS = job.Thr.Clus.ClusSize.ClusSig.Cth; nidm.Inferences(1).ExtentThreshold_type = 'obo_statistic'; - nidm.Inferences(1).ExtentThreshold_value = C_STCS; + nidm.Inferences(1).ExtentThreshold_clusterSizeInVoxels = C_STCS; case 'PthC' alpha_ucp = BoundCheck(job.Thr.Clus.ClusSize.ClusSig.PthC,[0 1],'Invalid uncorrected P(k)'); nidm.Inferences(1).ExtentThreshold_type = 'nidm_PValueUncorrected'; @@ -1491,10 +1491,10 @@ function snpm_pp(CWD,varargin) search_vol_cmm = S*abs(prod(VOX)); search_vol_vox = S; text(0,y,sprintf('Search vol: %d cmm, %d voxels',search_vol_cmm,search_vol_vox), 'FontSize',8) - nidm.SearchSpaceMaskMap_searchVolumeInVoxels = search_vol_vox; + nidm.Inferences(1).SearchSpaceMaskMap_searchVolumeInVoxels = search_vol_vox; % TODO convert back to units - nidm.SearchSpaceMaskMap_searchVolumeInUnits = search_vol_cmm; + nidm.Inferences(1).SearchSpaceMaskMap_searchVolumeInUnits = search_vol_cmm; y=y-0.8; text(0.7,y,sprintf('Voxel size: [%5.2f, %5.2f, %5.2f] mm',abs(VOX)), ... @@ -1571,10 +1571,10 @@ function snpm_pp(CWD,varargin) clear t % TODO SearchSpaceMaskMap can be different from analysis mask - nidm.SearchSpaceMaskMap_atLocation = 'mask.img'; + nidm.Inferences(1).SearchSpaceMaskMap_atLocation = 'mask.img'; % TODO: always stationary?? - nidm.SearchSpaceMaskMap_randomFieldStationarity = true; + nidm.Inferences(1).SearchSpaceMaskMap_randomFieldStationarity = true; nidm.Inferences(1).ExcursionSetMap_atLocation = Fname; @@ -1582,7 +1582,7 @@ function snpm_pp(CWD,varargin) nidm.NeuroimagingAnalysisSoftware_softwareVersion = snpm('ver'); % TODO: are other units possible in SnPM?? - nidm.CoordinateSpace_units = {'mm', 'mm', 'mm'}; + nidm.CoordinateSpace_voxelUnits = {'mm', 'mm', 'mm'}; % TODO: these should be filled in by the user --- nidm.CoordinateSpace_inWorldCoordinateSystem = 'nidm_StandardizedCoordinateSystem'; @@ -1593,9 +1593,9 @@ function snpm_pp(CWD,varargin) nidm.NeuroimagingAnalysisSoftware_label = 'SnPM'; % TODO: This temp file should only be produced if NIDM export is requested - jsonwrite('snpm_nidm_thresh.json', nidm_json, ... + jsonwrite('snpm_nidm_thresh.json', nidm, ... struct('indent',' ', 'escape', false)); - spm_nidmresults(nidm_json, CWD) + spm_nidmresults(nidm, CWD) end %-Reset Interactive Window diff --git a/snpm_ui.m b/snpm_ui.m index ce38701..a1498cf 100644 --- a/snpm_ui.m +++ b/snpm_ui.m @@ -310,6 +310,8 @@ function snpm_ui(varargin) nPerm_max = nPiStud_mx; elseif exist('nPiSubj_mx', 'var') nPerm_max = nPiSubj_mx; +elseif exist('nOfPerm', 'var') + nPerm_max = nOfPerm; else error('snpm:MissingMaxPerm', 'Maximum number of permutations not found') end @@ -598,7 +600,7 @@ function snpm_ui(varargin) s_SnPMcfg_save = ['s_SnPMcfg_save H C B G HCBGnames P PiCond ',... 'sPiCond bhPerms sHCform iGloNorm sGloNorm GM rg GX GMscale CONT ',... 'THRESH MASK ImMASK TH bVarSm vFWHM sVarSm bVolm bST sDesFile sDesign ',... - 'V pU_ST_Ut df1 nidm nPerm_max nidm_json ', ... + 'V pU_ST_Ut df1 nidm nPerm_max nidm ', ... 'sDesSave ',sDesSave]; eval(['save SnPMcfg ',s_SnPMcfg_save]) From b3656ff83572452dc410a02431706ea5b12d59ef Mon Sep 17 00:00:00 2001 From: cmaumet Date: Fri, 7 Jul 2017 11:09:59 +0100 Subject: [PATCH 24/29] Add user-specified info: space, scanner type, groups --- config/snpm_bch_pp.m | 103 +++++++++++++++++++++++++++++++- snpm_pp.m | 70 +++++++++++++++++++++- test/common/generic_test_snpm.m | 66 ++++++++++++++++++++ 3 files changed, 236 insertions(+), 3 deletions(-) diff --git a/config/snpm_bch_pp.m b/config/snpm_bch_pp.m index 5d803a4..b3aaf2a 100644 --- a/config/snpm_bch_pp.m +++ b/config/snpm_bch_pp.m @@ -234,10 +234,111 @@ '"FDR report" shows the histogram of voxel-wise uncorrected P-values, and the log-log plot of observed P-values versus (null hypothesis) expected P-values. The latter is the basis of the voxel-wise FDR-corrected P-values.' }; +% ---- Export to NIDM, adapted from SPM's spm_cfg_esults +%-------------------------------------------------------------------------- +% nsubj Number of subjects +%-------------------------------------------------------------------------- +nsubj = cfg_entry; +nsubj.tag = 'nsubj'; +nsubj.name = 'Number of subjects'; +nsubj.help = {'Number of subjects.'}; +nsubj.strtype = 'r'; +nsubj.num = [1 1]; + +%-------------------------------------------------------------------------- +% label Label +%-------------------------------------------------------------------------- +grplabel = cfg_entry; +grplabel.tag = 'label'; +grplabel.name = 'Label'; +grplabel.help = {'Group label.'}; +grplabel.strtype = 's'; +grplabel.num = [0 Inf]; + +%-------------------------------------------------------------------------- +% group +%-------------------------------------------------------------------------- +group = cfg_branch; +group.tag = 'group'; +group.name = 'Group'; +group.val = {nsubj grplabel}; +group.help = {['Number of subjects and labels per group. ', ... + 'For a single subject analysis, enter "1" and "single subject".']}; + +%-------------------------------------------------------------------------- +% groups +%-------------------------------------------------------------------------- +groups = cfg_repeat; +groups.tag = 'groups'; +groups.name = 'Groups'; +groups.help = {['Number of groups. ', ... + 'For a single subject analysis, specify one group.']}; +groups.values = {group}; +groups.num = [1 Inf]; + +%-------------------------------------------------------------------------- +% modality Modality +%-------------------------------------------------------------------------- +modality = cfg_menu; +modality.tag = 'modality'; +modality.name = 'Modality'; +modality.help = {'Modality.'}; +modality.labels = {'Anatomical MRI',... + 'Functional MRI',... + 'Diffusion MRI',... + 'PET',... + 'SPECT',... + 'EEG',... + 'MEG' +}'; +modality.values = {'AMRI','FMRI','DMRI','PET','SPECT','EEG','MEG'}; + +%-------------------------------------------------------------------------- +% refspace Reference space +%-------------------------------------------------------------------------- +refspace = cfg_menu; +refspace.tag = 'refspace'; +refspace.name = 'Reference space'; +refspace.help = {['Reference space. For an experiment completed only ',... + 'within SPM, choose one of the first four options.']}; +refspace.labels = {'Subject space (no normalisation)',... + 'Normalised space (using segment)',... + 'Normalised space (using old segment)',... + 'Customised space',... + 'Other normalised MNI space',... + 'Other normalised Talairach space',... +}'; +refspace.values = {'subject','ixi','icbm','custom','mni','talairach'}; + +%-------------------------------------------------------------------------- +% export Export results to NIDM +%-------------------------------------------------------------------------- +export_no = cfg_const; +export_no.tag = 'export_no'; +export_no.name = 'No'; +export_no.val = {0}; +export_no.help = {'Do not export to NIDM.'}; + + +export_nidm = cfg_branch; +export_nidm.tag = 'nidm'; +export_nidm.name = 'Export to NIDM'; +export_nidm.val = {modality refspace groups}; +export_nidm.help = {'NIDM (Neuroimaging Data Model)'}; + +export = cfg_choice; +export.tag = 'export'; +export.name = 'Export results to NIDM'; +export.help = {['Export your results to NIDM and share them easily with your collaborators.']}; +export.values = {export_nidm, export_no}; +export.val = {export_nidm}; + +% ---- + snpmpp = cfg_exbranch; snpmpp.name = 'Inference'; snpmpp.tag = 'inference'; -snpmpp.val = {snpmres ThrType posneg WrtFilt Report}; +snpmpp.val = {snpmres ThrType posneg WrtFilt Report export}; snpmpp.prog = @snpm_run_pp; %snpmpp.vout snpmpp.help = {'Examine the results of the SnPM computation.'}; diff --git a/snpm_pp.m b/snpm_pp.m index ac71f2b..2a81862 100644 --- a/snpm_pp.m +++ b/snpm_pp.m @@ -1190,10 +1190,76 @@ function snpm_pp(CWD,varargin) end end -% Display only if *not* in command line mode -nidm_export=true; +% NIDM export +nidm_export=isfield(job.export, 'nidm'); + +% Display only if *not* in command line mode or for NIDM export if ~spm_get_defaults('cmdline') || nidm_export + if nidm_export + % ---- Code adapted from SPM's spm_results_nidm ----- + %-NIDM Export + %---------------------------------------------------------------------- + %-Reference space + %-------------------------------------------------------------------------- + switch job.export.nidm.refspace + case 'subject' + coordsys = 'nidm_SubjectCoordinateSystem'; + case 'ixi' + coordsys = 'nidm_Ixi549CoordinateSystem'; + case 'icbm' + coordsys = 'nidm_IcbmMni152LinearCoordinateSystem'; + case 'custom' + coordsys = 'nidm_CustomCoordinateSystem'; + case 'mni' + coordsys = 'nidm_MNICoordinateSystem'; + case 'talairach' + coordsys = 'nidm_TalairachCoordinateSystem'; + otherwise + error('Unknown reference space.'); + end + nidm.CoordinateSpace_inWorldCoordinateSystem = coordsys; + + %-Data modality + %-------------------------------------------------------------------------- + MRIProtocol = ''; + switch job.export.nidm.modality + case 'AMRI' + ImagingInstrument = 'nlx_Magneticresonanceimagingscanner'; + MRIProtocol = 'nlx_AnatomicalMRIprotocol'; + case 'FMRI' + ImagingInstrument = 'nlx_Magneticresonanceimagingscanner'; + MRIProtocol = 'nlx_FunctionalMRIprotocol'; + case 'DMRI' + ImagingInstrument = 'nlx_Magneticresonanceimagingscanner'; + MRIProtocol = 'nlx_DiffusionMRIprotocol'; + case 'PET' + ImagingInstrument = 'nlx_Positronemissiontomographyscanner'; + case 'SPECT' + ImagingInstrument = 'nlx_Singlephotonemissioncomputedtomographyscanner'; + case 'EEG' + ImagingInstrument = 'nlx_Electroencephalographymachine'; + case 'MEG' + ImagingInstrument = 'nlx_Magnetoencephalographymachine'; + otherwise + error('Unknown modality.'); + end + nidm.Imaginginstrument_type = ImagingInstrument; + if ~isempty(MRIProtocol) + nidm.Data_hasMRIProtocol = MRIProtocol; + end + + %-Subject/Group(s) + %-------------------------------------------------------------------------- + groups = job.export.nidm.group; + if ~isequal(groups.N,1) + for i=1:numel(groups.N) + nidm.Groups(i).studygrouppopulation_groupName = groups.name{i}; + nidm.Groups(i).studygrouppopulation_numberOfSubjects = groups.N(i); + end + end + end + %======================================================================= %-D I S P L A Y : Max report %======================================================================= diff --git a/test/common/generic_test_snpm.m b/test/common/generic_test_snpm.m index 5900867..42bbcc7 100644 --- a/test/common/generic_test_snpm.m +++ b/test/common/generic_test_snpm.m @@ -102,6 +102,12 @@ function update_basis_matlabbatch(testCase) testCase.matlabbatch{3}.spm.tools.snpm.inference.Tsign = 1; testCase.matlabbatch{3}.spm.tools.snpm.inference.WriteFiltImg.name = 'SnPM_filtered_10none.nii'; + % NIDM export + testCase.matlabbatch{3}.spm.tools.snpm.inference.export.nidm.modality = 'FMRI'; + testCase.matlabbatch{3}.spm.tools.snpm.inference.export.nidm.refspace = 'ixi'; + testCase.matlabbatch{3}.spm.tools.snpm.inference.export.nidm.group.nsubj = 18; + testCase.matlabbatch{3}.spm.tools.snpm.inference.export.nidm.group.label = 'Control'; + if testCase.compaWithSpm % SPM batch testCase.spmBatch{2}.spm.stats.fmri_est.spmmat(1) = cfg_dep; @@ -312,6 +318,12 @@ function additional_results(testCase) testCase.matlabbatch{end}.spm.tools.snpm.inference.Thr.Vox.VoxSig.TFth = 1.6; testCase.matlabbatch{end}.spm.tools.snpm.inference.WriteFiltImg.name = 'SnPMt_filtered_vox_unc_t16.nii'; + % NIDM export + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.modality = 'FMRI'; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.refspace = 'ixi'; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.group.nsubj = 18; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.group.label = 'Control'; + % FWE voxel-wise p<0.5 testCase.matlabbatch{end+1}.spm.tools.snpm.inference.SnPMmat(1) = cfg_dep; testCase.matlabbatch{end}.spm.tools.snpm.inference.SnPMmat(1).tname = 'SnPM.mat results file'; @@ -322,6 +334,12 @@ function additional_results(testCase) testCase.matlabbatch{end}.spm.tools.snpm.inference.Thr.Vox.VoxSig.FWEth = 0.1; testCase.matlabbatch{end}.spm.tools.snpm.inference.WriteFiltImg.name = 'SnPMt_filtered_vox_fwe_p10.nii'; + % NIDM export + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.modality = 'FMRI'; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.refspace = 'ixi'; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.group.nsubj = 18; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.group.label = 'Control'; + % FDR voxel-wise p<0.5 testCase.matlabbatch{end+1}.spm.tools.snpm.inference.SnPMmat(1) = cfg_dep; testCase.matlabbatch{end}.spm.tools.snpm.inference.SnPMmat(1).tname = 'SnPM.mat results file'; @@ -331,6 +349,12 @@ function additional_results(testCase) testCase.matlabbatch{end}.spm.tools.snpm.inference.SnPMmat(1).src_output = substruct('.','SnPM'); testCase.matlabbatch{end}.spm.tools.snpm.inference.Thr.Vox.VoxSig.FDRth = 0.7; testCase.matlabbatch{end}.spm.tools.snpm.inference.WriteFiltImg.name = 'SnPMt_filtered_vox_fdr_p70.nii'; + + % NIDM export + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.modality = 'FMRI'; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.refspace = 'ixi'; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.group.nsubj = 18; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.group.label = 'Control'; end function additional_cluster_results(testCase) @@ -345,6 +369,12 @@ function additional_cluster_results(testCase) testCase.matlabbatch{end}.spm.tools.snpm.inference.Thr.Clus.ClusSize.ClusSig.PthC = 0.1; testCase.matlabbatch{end}.spm.tools.snpm.inference.WriteFiltImg.name = 'SnPMt_filtered_clus_4_unc_p10.nii'; + % NIDM export + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.modality = 'FMRI'; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.refspace = 'ixi'; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.group.nsubj = 18; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.group.label = 'Control'; + % Uncorrected (cluster forming u>4) cluster-wise k>6 testCase.matlabbatch{end+1}.spm.tools.snpm.inference.SnPMmat(1) = cfg_dep; testCase.matlabbatch{end}.spm.tools.snpm.inference.SnPMmat(1).tname = 'SnPM.mat results file'; @@ -356,6 +386,12 @@ function additional_cluster_results(testCase) testCase.matlabbatch{end}.spm.tools.snpm.inference.Thr.Clus.ClusSize.ClusSig.Cth = 6; testCase.matlabbatch{end}.spm.tools.snpm.inference.WriteFiltImg.name = 'SnPMt_filtered_clus_4_unc_k6.nii'; + % NIDM export + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.modality = 'FMRI'; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.refspace = 'ixi'; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.group.nsubj = 18; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.group.label = 'Control'; + % FWE (cluster forming u>4) cluster-wise p<0.5 testCase.matlabbatch{end+1}.spm.tools.snpm.inference.SnPMmat(1) = cfg_dep; testCase.matlabbatch{end}.spm.tools.snpm.inference.SnPMmat(1).tname = 'SnPM.mat results file'; @@ -367,6 +403,12 @@ function additional_cluster_results(testCase) testCase.matlabbatch{end}.spm.tools.snpm.inference.Thr.Clus.ClusSize.ClusSig.FWEthC = 0.5; testCase.matlabbatch{end}.spm.tools.snpm.inference.WriteFiltImg.name = 'SnPMt_filtered_clus_4_fwe_p50.nii'; + % NIDM export + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.modality = 'FMRI'; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.refspace = 'ixi'; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.group.nsubj = 18; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.group.label = 'Control'; + % FWE (cluster forming u>5) cluster-wise p<0.5 testCase.matlabbatch{end+1}.spm.tools.snpm.inference.SnPMmat(1) = cfg_dep; testCase.matlabbatch{end}.spm.tools.snpm.inference.SnPMmat(1).tname = 'SnPM.mat results file'; @@ -377,6 +419,12 @@ function additional_cluster_results(testCase) testCase.matlabbatch{end}.spm.tools.snpm.inference.Thr.Clus.ClusSize.CFth = 5; testCase.matlabbatch{end}.spm.tools.snpm.inference.Thr.Clus.ClusSize.ClusSig.FWEthC = 0.5; testCase.matlabbatch{end}.spm.tools.snpm.inference.WriteFiltImg.name = 'SnPMt_filtered_clus_5_fwe_p50.nii'; + + % NIDM export + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.modality = 'FMRI'; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.refspace = 'ixi'; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.group.nsubj = 18; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.group.label = 'Control'; end function additional_predefined_cluster_results(testCase) @@ -394,6 +442,12 @@ function additional_predefined_cluster_results(testCase) testCase.matlabbatch{end}.spm.tools.snpm.inference.Thr.Clus.ClusSize.ClusSig.PthC = 0.1; testCase.matlabbatch{end}.spm.tools.snpm.inference.WriteFiltImg.name = 'SnPMt_filtered_clus_p10_unc_p10.nii'; + % NIDM export + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.modality = 'FMRI'; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.refspace = 'ixi'; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.group.nsubj = 18; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.group.label = 'Control'; + % FWE (cluster forming p<0.1) cluster-wise p<0.5 testCase.matlabbatch{end+1}.spm.tools.snpm.inference.SnPMmat(1) = cfg_dep; testCase.matlabbatch{end}.spm.tools.snpm.inference.SnPMmat(1).tname = 'SnPM.mat results file'; @@ -404,6 +458,12 @@ function additional_predefined_cluster_results(testCase) testCase.matlabbatch{end}.spm.tools.snpm.inference.Thr.Clus.ClusSize.CFth = NaN; testCase.matlabbatch{end}.spm.tools.snpm.inference.Thr.Clus.ClusSize.ClusSig.FWEthC = 0.5; testCase.matlabbatch{end}.spm.tools.snpm.inference.WriteFiltImg.name = 'SnPMt_filtered_clus_p10_fwe_p50.nii'; + + % NIDM export + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.modality = 'FMRI'; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.refspace = 'ixi'; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.group.nsubj = 18; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.group.label = 'Control'; end function additional_cluster_mass_results(testCase) @@ -418,6 +478,12 @@ function additional_cluster_mass_results(testCase) testCase.matlabbatch{end}.spm.tools.snpm.inference.Thr.Clus.ClusMass.Theta = 0.5; testCase.matlabbatch{end}.spm.tools.snpm.inference.Tsign = 1; testCase.matlabbatch{end}.spm.tools.snpm.inference.WriteFiltImg.name = 'SnPMt_filtered_cluss_mass.nii'; + + % NIDM export + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.modality = 'FMRI'; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.refspace = 'ixi'; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.group.nsubj = 18; + testCase.matlabbatch{end}.spm.tools.snpm.inference.export.nidm.group.label = 'Control'; end function compare_batch_with_inter(testCase, zeroingNaNs) From 04e4ee026b51c376a726f1cbf3fdab4a365eb3e9 Mon Sep 17 00:00:00 2001 From: cmaumet Date: Fri, 7 Jul 2017 11:17:50 +0100 Subject: [PATCH 25/29] Rename tags in agreement with batch --- snpm_pp.m | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/snpm_pp.m b/snpm_pp.m index 2a81862..bf1832d 100644 --- a/snpm_pp.m +++ b/snpm_pp.m @@ -1252,10 +1252,10 @@ function snpm_pp(CWD,varargin) %-Subject/Group(s) %-------------------------------------------------------------------------- groups = job.export.nidm.group; - if ~isequal(groups.N,1) - for i=1:numel(groups.N) - nidm.Groups(i).studygrouppopulation_groupName = groups.name{i}; - nidm.Groups(i).studygrouppopulation_numberOfSubjects = groups.N(i); + if ~isequal(groups.nsubj,1) + for i=1:numel(groups) + nidm.Groups(i).studygrouppopulation_groupName = groups(i).label; + nidm.Groups(i).studygrouppopulation_numberOfSubjects = groups(i).nsubj; end end end From d87429c74f8f982b8a22d974ca3cff4b753cba1c Mon Sep 17 00:00:00 2001 From: cmaumet Date: Fri, 7 Jul 2017 11:37:03 +0100 Subject: [PATCH 26/29] Test results with no NIDM export --- test/test_snpm_pp.m | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/test_snpm_pp.m b/test/test_snpm_pp.m index ab9a3e0..c9b1422 100644 --- a/test/test_snpm_pp.m +++ b/test/test_snpm_pp.m @@ -81,6 +81,8 @@ function test_onesample_cluster_err_bigfile(testCase) testCase.matlabbatch{1}.spm.tools.snpm.inference.Thr.Clus.ClusSize.CFth = 4; testCase.matlabbatch{1}.spm.tools.snpm.inference.Thr.Clus.ClusSize.ClusSig.PthC = 0.1; testCase.matlabbatch{1}.spm.tools.snpm.inference.WriteFiltImg.name = 'SnPMt_filtered_clus_4_unc_p10.nii'; + testCase.matlabbatch{1}.spm.tools.snpm.inference.export.export_no = 0; + testCase.warningId = 'SnPM:SnPMSTFileNotLOaded'; end From cd7a689932924a74076034db5db2d8d8edd71112 Mon Sep 17 00:00:00 2001 From: cmaumet Date: Fri, 7 Jul 2017 14:33:07 +0100 Subject: [PATCH 27/29] Update camel case for OBO and NeuroLex terms --- snpm_cp.m | 14 +++++++------- snpm_pp.m | 36 ++++++++++++++++++------------------ 2 files changed, 25 insertions(+), 25 deletions(-) diff --git a/snpm_cp.m b/snpm_cp.m index 5cd6904..59ea13e 100644 --- a/snpm_cp.m +++ b/snpm_cp.m @@ -217,22 +217,22 @@ function snpm_cp(CWD) if strcmp(sDesFile, 'snpm_pi_OneSampT') || ... strcmp(sDesFile, 'snpm_pi_ANOVAwithinS') % Sign flipping - nidm.ErrorModel_hasErrorDistribution = {'obo_nonparametricdistribution', 'obo_symmetricdistribution'}; + nidm.ErrorModel_hasErrorDistribution = {'obo_NonParametricDistribution', 'obo_SymmetricDistribution'}; nidm.ErrorModel_errorVarianceHomogeneous = false; nidm.ErrorModel_varianceMapWiseDependence = 'nidm_IndependentParameter'; nidm.ErrorModel_hasErrorDependence = 'nidm_IndependentError'; else % Permutation - nidm.ErrorModel_hasErrorDistribution = 'obo_nonparametricdistribution'; + nidm.ErrorModel_hasErrorDistribution = 'obo_NonParametricDistribution'; nidm.ErrorModel_errorVarianceHomogeneous = true; nidm.ErrorModel_varianceMapWiseDependence = 'nidm_IndependentParameter'; % TODO: the 'obo_exchangeable' term is not yet in STATO - nidm.ErrorModel_hasErrorDependence = 'obo_exchangeable'; + nidm.ErrorModel_hasErrorDependence = 'obo_Exchangeable'; nidm.ErrorModel_dependenceMapWiseDependence = 'nidm_IndependentParameter'; end % TODO: check this is correct -nidm.ModelParameterEstimation_withEstimationMethod = 'obo_ordinaryleastsquaresestimation'; +nidm.ModelParameterEstimation_withEstimationMethod = 'obo_OrdinaryLeastSquaresEstimation'; if isempty([H C]) error('SnPM:NoModel', 'No model specified; [H C] empty'); @@ -249,7 +249,7 @@ function snpm_cp(CWD) warning('SnPM:FContrast', ... 'F contrast! F statistic images are being created.'); STAT = 'F'; - nidm.Contrasts(1).StatisticMap_statisticType = 'obo_Fstatistic'; + nidm.Contrasts(1).StatisticMap_statisticType = 'obo_FStatistic'; if (CONT(1,:) == -CONT(2,:)) CONT = CONT(1,:); @@ -269,8 +269,8 @@ function snpm_cp(CWD) nidm.Contrasts(1).StatisticMap_statisticType = 'nidm_smoothedtstatistic'; nidm.Contrasts(2).StatisticMap_statisticType = 'nidm_smoothedtstatistic'; else - nidm.Contrasts(1).StatisticMap_statisticType = 'obo_tstatistic'; - nidm.Contrasts(2).StatisticMap_statisticType = 'obo_tstatistic'; + nidm.Contrasts(1).StatisticMap_statisticType = 'obo_TStatistic'; + nidm.Contrasts(2).StatisticMap_statisticType = 'obo_TStatistic'; end nidm.Contrasts(1).ContrastMap_contrastName = ['Positive T-Contrast: [' mat2str(CONT) ']']; nidm.Contrasts(1).contrastweightmatrix_value = CONT; diff --git a/snpm_pp.m b/snpm_pp.m index bf1832d..fafced4 100644 --- a/snpm_pp.m +++ b/snpm_pp.m @@ -505,12 +505,12 @@ function snpm_pp(CWD,varargin) case 'TFth' u = BoundCheck(job.Thr.Vox.VoxSig.TFth,[0 Inf],'Negative Threshold!'); alph_FWE = sum(MaxT > u -tol) / nPermReal; - nidm.Inferences(1).HeightThreshold_type = 'obo_statistic'; + nidm.Inferences(1).HeightThreshold_type = 'obo_Statistic'; nidm.Inferences(1).HeightThreshold_value = u; case 'FDRth' alph_FDR = BoundCheck(job.Thr.Vox.VoxSig.FDRth,[0 1],'Invalid FDR level'); alpha_ucp = snpm_uc_FDR(alph_FDR,[],'P',[],sSnPMucp'); - nidm.Inferences(1).HeightThreshold_type = 'obo_qvalue'; + nidm.Inferences(1).HeightThreshold_type = 'obo_QValue'; nidm.Inferences(1).HeightThreshold_value = alph_FDR; case 'FWEth' alph_FWE = BoundCheck(job.Thr.Vox.VoxSig.FWEth,[0 1],'Invalid FWE level'); @@ -521,13 +521,13 @@ function snpm_pp(CWD,varargin) C_MaxT = 0; end u = C_MaxT; - nidm.Inferences(1).HeightThreshold_type = 'obo_FWERadjustedpvalue'; + nidm.Inferences(1).HeightThreshold_type = 'obo_FWERAdjustedPValue'; nidm.Inferences(1).HeightThreshold_value = alph_FWE; otherwise error('Unknown threshold') end % No extent thresholding when voxelwise threshold is requested - nidm.Inferences(1).ExtentThreshold_type = 'obo_statistic'; + nidm.Inferences(1).ExtentThreshold_type = 'obo_Statistic'; nidm.Inferences(1).ExtentThreshold_clusterSizeInVoxels = 0; else % Cluster-wise inference @@ -547,7 +547,7 @@ function snpm_pp(CWD,varargin) nidm.Inferences(1).HeightThreshold_type = 'nidm_PValueUncorrected'; nidm.Inferences(1).HeightThreshold_value = CFth; else - nidm.Inferences(1).HeightThreshold_type = 'obo_statistic'; + nidm.Inferences(1).HeightThreshold_type = 'obo_Statistic'; nidm.Inferences(1).HeightThreshold_value = CFth; end @@ -612,7 +612,7 @@ function snpm_pp(CWD,varargin) nidm.Inferences(1).HeightThreshold_type = 'nidm_PValueUncorrected'; nidm.Inferences(1).HeightThreshold_value = job.Thr.Clus.ClusSize.CFth; else - nidm.Inferences(1).HeightThreshold_type = 'obo_statistic'; + nidm.Inferences(1).HeightThreshold_type = 'obo_Statistic'; nidm.Inferences(1).HeightThreshold_value = job.Thr.Clus.ClusSize.CFth; end @@ -623,7 +623,7 @@ function snpm_pp(CWD,varargin) switch tmp{1} case 'Cth' C_STCS = job.Thr.Clus.ClusSize.ClusSig.Cth; - nidm.Inferences(1).ExtentThreshold_type = 'obo_statistic'; + nidm.Inferences(1).ExtentThreshold_type = 'obo_Statistic'; nidm.Inferences(1).ExtentThreshold_clusterSizeInVoxels = C_STCS; case 'PthC' alpha_ucp = BoundCheck(job.Thr.Clus.ClusSize.ClusSig.PthC,[0 1],'Invalid uncorrected P(k)'); @@ -632,7 +632,7 @@ function snpm_pp(CWD,varargin) case 'FWEthC' alph_FWE = BoundCheck(job.Thr.Clus.ClusSize.ClusSig.FWEthC,[0 1],'Invalid FWE level (cluster-level inference)'); iFWE = ceil((1-alph_FWE)*nPermReal); - nidm.Inferences(1).ExtentThreshold_type = 'obo_FWERadjustedpvalue'; + nidm.Inferences(1).ExtentThreshold_type = 'obo_FWERAdjustedPValue'; nidm.Inferences(1).ExtentThreshold_value = alph_FWE; end end % END: Cluster-wise inference @@ -1225,22 +1225,22 @@ function snpm_pp(CWD,varargin) MRIProtocol = ''; switch job.export.nidm.modality case 'AMRI' - ImagingInstrument = 'nlx_Magneticresonanceimagingscanner'; - MRIProtocol = 'nlx_AnatomicalMRIprotocol'; + ImagingInstrument = 'nlx_MagneticResonanceImagingScanner'; + MRIProtocol = 'nlx_AnatomicalMRIProtocol'; case 'FMRI' - ImagingInstrument = 'nlx_Magneticresonanceimagingscanner'; - MRIProtocol = 'nlx_FunctionalMRIprotocol'; + ImagingInstrument = 'nlx_MagneticResonanceImagingScanner'; + MRIProtocol = 'nlx_FunctionalMRIProtocol'; case 'DMRI' - ImagingInstrument = 'nlx_Magneticresonanceimagingscanner'; - MRIProtocol = 'nlx_DiffusionMRIprotocol'; + ImagingInstrument = 'nlx_MagneticResonanceImagingScanner'; + MRIProtocol = 'nlx_DiffusionWeightedImagingProtocol'; case 'PET' - ImagingInstrument = 'nlx_Positronemissiontomographyscanner'; + ImagingInstrument = 'nlx_PositronEmissionTomographyScanner'; case 'SPECT' - ImagingInstrument = 'nlx_Singlephotonemissioncomputedtomographyscanner'; + ImagingInstrument = 'nlx_SinglePhotonEmissionComputedTomographyScanner'; case 'EEG' - ImagingInstrument = 'nlx_Electroencephalographymachine'; + ImagingInstrument = 'nlx_ElectroencephalographyMachine'; case 'MEG' - ImagingInstrument = 'nlx_Magnetoencephalographymachine'; + ImagingInstrument = 'nlx_MagnetoencephalographyMachine'; otherwise error('Unknown modality.'); end From e33f18c7b3b0cf1b412c7c3eca33f46a2762045e Mon Sep 17 00:00:00 2001 From: cmaumet Date: Fri, 7 Jul 2017 14:33:29 +0100 Subject: [PATCH 28/29] Fix: info now provided by the user in the batch --- snpm_pp.m | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/snpm_pp.m b/snpm_pp.m index fafced4..4c3a46d 100644 --- a/snpm_pp.m +++ b/snpm_pp.m @@ -1649,12 +1649,7 @@ function snpm_pp(CWD,varargin) % TODO: are other units possible in SnPM?? nidm.CoordinateSpace_voxelUnits = {'mm', 'mm', 'mm'}; - - % TODO: these should be filled in by the user --- - nidm.CoordinateSpace_inWorldCoordinateSystem = 'nidm_StandardizedCoordinateSystem'; - nidm.Imaginginstrument_type = 'nlx_Magneticresonanceimagingscanner'; - % --- - + nidm.NeuroimagingAnalysisSoftware_type = 'src_SnPM'; nidm.NeuroimagingAnalysisSoftware_label = 'SnPM'; From 111b09495dae6e5e6a07913a7926465d42e2b942 Mon Sep 17 00:00:00 2001 From: cmaumet Date: Fri, 7 Jul 2017 15:25:47 +0100 Subject: [PATCH 29/29] More use of camel case for OBO and NeuroLex terms --- snpm_cp.m | 6 +++--- snpm_pp.m | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/snpm_cp.m b/snpm_cp.m index 59ea13e..983cf1c 100644 --- a/snpm_cp.m +++ b/snpm_cp.m @@ -245,7 +245,7 @@ function snpm_cp(CWD) end if size(CONT,1) > 1 % F-contrast - nidm.Contrasts(1).contrastweightmatrix_value = CONT; + nidm.Contrasts(1).ContrastWeightMatrix_value = CONT; warning('SnPM:FContrast', ... 'F contrast! F statistic images are being created.'); STAT = 'F'; @@ -273,9 +273,9 @@ function snpm_cp(CWD) nidm.Contrasts(2).StatisticMap_statisticType = 'obo_TStatistic'; end nidm.Contrasts(1).ContrastMap_contrastName = ['Positive T-Contrast: [' mat2str(CONT) ']']; - nidm.Contrasts(1).contrastweightmatrix_value = CONT; + nidm.Contrasts(1).ContrastWeightMatrix_value = CONT; nidm.Contrasts(2).ContrastMap_contrastName = ['Negative T-Contrast: [' mat2str(-CONT) ']']; - nidm.Contrasts(2).contrastweightmatrix_value = -CONT; + nidm.Contrasts(2).ContrastWeightMatrix_value = -CONT; end if rank(CONT)