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_cp.m b/snpm_cp.m index e322346..983cf1c 100644 --- a/snpm_cp.m +++ b/snpm_cp.m @@ -213,6 +213,27 @@ function snpm_cp(CWD) %-Load config file & catch all problem cases now %----------------------------------------------------------------------- load(CfgFile); + +if strcmp(sDesFile, 'snpm_pi_OneSampT') || ... + strcmp(sDesFile, 'snpm_pi_ANOVAwithinS') + % Sign flipping + 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_errorVarianceHomogeneous = true; + nidm.ErrorModel_varianceMapWiseDependence = 'nidm_IndependentParameter'; + % TODO: the 'obo_exchangeable' term is not yet in STATO + nidm.ErrorModel_hasErrorDependence = 'obo_Exchangeable'; + nidm.ErrorModel_dependenceMapWiseDependence = 'nidm_IndependentParameter'; +end + +% TODO: check this is correct +nidm.ModelParameterEstimation_withEstimationMethod = 'obo_OrdinaryLeastSquaresEstimation'; + if isempty([H C]) error('SnPM:NoModel', 'No model specified; [H C] empty'); end @@ -222,16 +243,41 @@ 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 + +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'; + nidm.Contrasts(1).StatisticMap_statisticType = 'obo_FStatistic'; + if (CONT(1,:) == -CONT(2,:)) CONT = CONT(1,:); end + con_name = 'Positive'; + nidm.Contrasts(1).StatisticMap_contrastName = con_name; + con_neg_name = ''; else + 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 + 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'; + end + 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)TH) & any(diff(X)) & Wt); end + if length(Q) @@ -531,6 +637,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); @@ -566,8 +673,12 @@ 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' + 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 @@ -575,6 +686,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 @@ -604,13 +717,21 @@ function snpm_cp(CWD) % %- New! Write out data images. %- Input image data. + 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; + + CON_pos_image(:,Q)=CON; + CON_neg_image(:,Q)=-CON; + + CONSE_image(:,Q)=CONSE; elseif STAT=='F' F_image(:,Q)=T; + Fnum_image(:,Q)=Fnum; end if bVarAlph @@ -643,16 +764,36 @@ 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(mask_image,DIM(1),DIM(2),DIM(3)); + spm_write_vol(Vmask,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)); spm_write_vol(VT_pos,T_pos_vol); 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); + + 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); + + Fnum_vol=reshape(Fnum_image,DIM(1),DIM(2),DIM(3)); + spm_write_vol(VFnum,Fnum_vol); end if bVarAlph @@ -660,7 +801,15 @@ function snpm_cp(CWD) spm_write_vol(VlwP, lwP_vol); end - else + else + % Analysis mask + 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)); + 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); @@ -675,10 +824,22 @@ 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,i); + + CON_neg_plate=reshape(CON_neg_image,DIM(1),DIM(2)); + spm_write_plane(VCON_neg,CON_neg_plate,i); + + CONSE_plate=reshape(CONSE_image,DIM(1),DIM(2)); + spm_write_plane(VCONSE,CONSE_plate,i); 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 diff --git a/snpm_pp.m b/snpm_pp.m index 2841c8e..fef1f82 100644 --- a/snpm_pp.m +++ b/snpm_pp.m @@ -370,7 +370,6 @@ function snpm_pp(CWD,varargin) load(fullfile(CWD,'SnPM')) load(fullfile(CWD,'SnPMucp')) - %-Ask whether positive or negative effects be analysed %----------------------------------------------------------------------- if BATCH @@ -389,6 +388,13 @@ function snpm_pp(CWD,varargin) end end +% Find corresponding contrast name in nidm json structure +if bNeg == 0; + nidm.Inferences(1).StatisticMap_contrastName = {con_name}; +else + nidm.Inferences(1).StatisticMap_contrastName = {con_neg_name}; +end + %-Take MaxT for increases or decreases according to bNeg MaxT = MaxT(:,bNeg+1); nPerm = size(PiCond,1); %nPerm is consistent with the one in snpm_cp @@ -484,6 +490,7 @@ 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.Inferences(1).Inference_hasAlternativeHypothesis = 'nidm_OneTailedTest'; if BATCH bSpatEx = isfield(job.Thr,'Clus'); if ~bSpatEx @@ -493,12 +500,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.Inferences(1).HeightThreshold_type = 'nidm_PValueUncorrected'; + nidm.Inferences(1).HeightThreshold_value = alpha_ucp; 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_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_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 +521,14 @@ function snpm_pp(CWD,varargin) C_MaxT = 0; end u = C_MaxT; + 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_clusterSizeInVoxels = 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.Inferences(1).HeightThreshold_type = 'nidm_PValueUncorrected'; + nidm.Inferences(1).HeightThreshold_value = CFth; + else + nidm.Inferences(1).HeightThreshold_type = 'obo_Statistic'; + nidm.Inferences(1).HeightThreshold_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 @@ -575,9 +603,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.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_value = job.Thr.Clus.ClusSize.CFth; + end + end u=ST_Ut; % Flag use of a statistic-value threshold % Inference details... @@ -585,14 +623,20 @@ 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_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'; + 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.Inferences(1).ExtentThreshold_type = 'obo_FWERAdjustedPValue'; + nidm.Inferences(1).ExtentThreshold_value = alph_FWE; end end % END: Cluster-wise inference - + else % GUI, interative inference specification str_img =[STAT,'|P']; @@ -772,6 +816,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.ClusterDefinitionCriteria_hasConnectivityCriterion = 'nidm_voxel18connected'; + if bSpatEx fprintf('Working on spatial extent...\n'); %-Compute suprathreshold voxels - check there are some @@ -896,6 +944,7 @@ 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,:)); + if isPos==1 STCstats_Pos = [ SnPM_ST(:,tQ); tmp]; if bNeg==0 @@ -1141,8 +1190,75 @@ function snpm_pp(CWD,varargin) end end -% Display only if *not* in command line mode -if ~spm_get_defaults('cmdline') +% 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_DiffusionWeightedImagingProtocol'; + 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.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 %======================================================================= %-D I S P L A Y : Max report @@ -1311,6 +1427,7 @@ function snpm_pp(CWD,varargin) %----------------------------------------------------------------------- r = 1; bUsed = zeros(size(STC_SnPMt)); + while max(STC_SnPMt.*(~bUsed)) & (y > 3) [null, i] = max(STC_SnPMt.*(~bUsed)); % Largest t value @@ -1341,10 +1458,23 @@ 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{:}) + + % 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.PeakDefinitionCriteria_minDistanceBetweenPeaks = 8; + nidm.PeakDefinitionCriteria_maxNumberOfPeaksPerCluster = 3; + [null, k] = sort(-STC_SnPMt(j)); % Sort on t value D = i; for i = 1:length(k) @@ -1365,14 +1495,22 @@ function snpm_pp(CWD,varargin) y = y -1; end end + + + 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 - + bUsed(j) = (bUsed(j) | 1 ); %-Mark maxima as "used" r = r + 1; % Next region end + clear i j k D d r - %-Footnote with SnPM parameters %======================================================================= line([0,1],[0.5,0.5],'LineWidth',1,'Color','r') @@ -1416,7 +1554,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.Inferences(1).SearchSpaceMaskMap_searchVolumeInVoxels = search_vol_vox; + +% TODO convert back to units + 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)), ... 'FontSize', 8) @@ -1445,7 +1590,7 @@ function snpm_pp(CWD,varargin) %- Image output? %======================================================================= %-Write out filtered SnPMt? -if WrtFlt +if WrtFlt || nidm_export Fname = WrtFltFn; @@ -1490,6 +1635,28 @@ function snpm_pp(CWD,varargin) end Vs = sf_close_vol(Vs); clear t + + % TODO SearchSpaceMaskMap can be different from analysis mask + nidm.Inferences(1).SearchSpaceMaskMap_atLocation = 'mask.img'; + + % TODO: always stationary?? + nidm.Inferences(1).SearchSpaceMaskMap_randomFieldStationarity = true; + + nidm.Inferences(1).ExcursionSetMap_atLocation = Fname; + + nidm.NIDMResultsExporter_softwareVersion = snpm('ver'); + nidm.NeuroimagingAnalysisSoftware_softwareVersion = snpm('ver'); + + % TODO: are other units possible in SnPM?? + nidm.CoordinateSpace_voxelUnits = {'mm', 'mm', 'mm'}; + + 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, ... + struct('indent',' ', 'escape', false)); + spm_nidmresults(nidm, CWD) end %-Reset Interactive Window @@ -1631,6 +1798,7 @@ function ShowDist(T,cT,aT,C,cC,aC,Typ) end + return diff --git a/snpm_ui.m b/snpm_ui.m index 7f146bf..a1498cf 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 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 = struct(); + %-Definitions & Design parameters %======================================================================= sDesigns=str2mat(... @@ -299,6 +303,19 @@ 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; +elseif exist('nOfPerm', 'var') + nPerm_max = nOfPerm; +else + error('snpm:MissingMaxPerm', 'Maximum number of permutations not found') +end + %-Get general analysis & data parameters %======================================================================= @@ -368,6 +385,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.Data_grandMeanScaling = true; if (iGloNorm==2) % Proportional scaling str = 'PropSca global mean to'; else @@ -381,8 +399,10 @@ function snpm_ui(varargin) case 'gmsca_no', GM = 50; end + nidm.Data_targetIntensity = GM; elseif (iGMsca==1) % No grand mean scaling GM = 0; + nidm.Data_grandMeanScaling = false; end @@ -572,18 +592,19 @@ 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.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 ', ... + 'V pU_ST_Ut df1 nidm nPerm_max nidm ', ... 'sDesSave ',sDesSave]; eval(['save SnPMcfg ',s_SnPMcfg_save]) - if ~spm_get_defaults('cmdline') %======================================================================= 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 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) 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