Skip to content

Instantly share code, notes, and snippets.

@swederik

swederik/conn_vproject.m

Created June 11, 2014 10:08
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save swederik/e23c3726ce311fab0bf7 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save swederik/e23c3726ce311fab0bf7 to your computer and use it in GitHub Desktop.
Download ZIP
CONN functional connectivity toolbox ROI-to-voxel stats exporting on Mac OSX (conn_vproject.m)
Raw
conn_vproject.m
function [dataplot,infoplot,data1plot]=conn_vproject(a,b,c,d,views,projection,thres,res,box,select,mat,threshold,data1plot,spmfile,voxeltovoxel)
%%%%%
% Slight modification to CONN functional connectivity toolbox to fix stats exporting on Mac OSX
% uicontextmenu appears broken on OSX. May be due to resolution issues on retina macbooks?
%%%%%
%CONN_VPROJECT Volume display
%
if numel(a)==1 && ishandle(a), % callbacks from UI objects
init=0;
GCF=gcf;
if isstruct(b), ARG=b; end;
OPTION=c;
if nargin>3, OPTION2=d; else OPTION2=''; end
if isempty(OPTION), return; end
DATA=get(GCF,'userdata');
a=DATA.a;
b=DATA.b;
c=DATA.c;
d=DATA.d;
projection=DATA.projection;
thres=DATA.thres;
threshold=DATA.threshold;
res=DATA.res;
box=DATA.box;
views=DATA.views;
select=DATA.select;
mat=DATA.mat;
spmfile=DATA.spmfile;
if isfield(DATA,'voxeltovoxel'), voxeltovoxel=DATA.voxeltovoxel; else voxeltovoxel=0; end
data1plot=[];
if isstruct(OPTION),
else,
switch(OPTION),
case 'buttondown',
xlim=get(DATA.axes,'xlim');ylim=get(DATA.axes,'ylim');xpos=get(DATA.axes,'currentpoint');
if xpos(1,1)>=xlim(1)&&xpos(1,1)<=xlim(2)&&xpos(1,2)>=ylim(1)&&xpos(1,2)<=ylim(2),
DATA.buttondown=get(0,'pointerlocation');
DATA.selectiontype=get(GCF,'selectiontype');
set(GCF,'windowbuttonmotionfcn',{@conn_vproject,'buttonmotion'},'userdata',DATA);
res=DATA.res*3/2;
box=1;
else, DATA.buttondown=[]; set(GCF,'userdata',DATA); return; end
case 'buttonup',
p1=DATA.buttondown;
if ~isempty(p1),
p2=get(0,'pointerlocation')-p1;
if strcmp(DATA.selectiontype,'extend'),
ang=-.005*(p2(2)+p2(1));projection([1,2],:)=[cos(ang),sin(ang);-sin(ang),cos(ang)]*projection([1,2],:);
else,
%ang=.01*p2(2);projection([2,3],:)=[cos(ang),sin(ang);-sin(ang),cos(ang)]*projection([2,3],:);
ang=-.01*p2(1);projection([1,3],:)=[cos(ang),sin(ang);-sin(ang),cos(ang)]*projection([1,3],:);
end
DATA.projection=projection;
set(GCF,'windowbuttonmotionfcn',[],'userdata',DATA);
else, return; end
case 'buttonmotion',
p1=DATA.buttondown;
p2=get(0,'pointerlocation')-p1;
if strcmp(DATA.selectiontype,'extend'),
ang=-.005*(p2(2)+p2(1));projection([1,2],:)=[cos(ang),sin(ang);-sin(ang),cos(ang)]*projection([1,2],:);
else,
%ang=.01*p2(2);projection([2,3],:)=[cos(ang),sin(ang);-sin(ang),cos(ang)]*projection([2,3],:);
ang=-.01*p2(1);projection([1,3],:)=[cos(ang),sin(ang);-sin(ang),cos(ang)]*projection([1,3],:);
end
res=DATA.res*3/2;
box=1;
case 'resolution',
res=str2num(get(DATA.handles(2),'string'));
if isempty(res), res=DATA.res; end
set(DATA.handles(2),'string',num2str(res));
DATA.res=res;
set(GCF,'userdata',DATA);
init=-1;
case 'threshold',
threshold=str2num(get(DATA.handles(4),'string'));
if isempty(threshold), threshold=DATA.threshold; end
set(DATA.handles(4),'string',num2str(threshold));
DATA.threshold=threshold;
set(GCF,'userdata',DATA);
init=-1;
case 'keypress',
switch(lower(ARG.Character)),
case 'y',
DATA.rotation=DATA.rotation+[10,0,0]*(2*(length(ARG.Modifier)>0)-1);
case 'x',
DATA.rotation=DATA.rotation+[0,10,0]*(2*(length(ARG.Modifier)>0)-1);
case 'z',
DATA.rotation=DATA.rotation+[0,0,10]*(2*(length(ARG.Modifier)>0)-1);
otherwise,
return;
end
ang=-.01*DATA.rotation(1);projection([1,3],:)=[cos(ang),sin(ang);-sin(ang),cos(ang)]*projection([1,3],:);
ang=.01*DATA.rotation(2);projection([2,3],:)=[cos(ang),sin(ang);-sin(ang),cos(ang)]*projection([2,3],:);
ang=-.01*DATA.rotation(3);projection([1,2],:)=[cos(ang),sin(ang);-sin(ang),cos(ang)]*projection([1,2],:);
DATA.rotation=[0,0,0];
DATA.projection=projection;
set(GCF,'windowbuttonmotionfcn',[],'userdata',DATA);
case 'selectroi'
selectcluster=get(DATA.handles(8),'value');
if selectcluster>length(DATA.clusters), select=[]; else, select=DATA.clusters{selectcluster}; end
DATA.select=select; DATA.selectcluster=selectcluster; set(GCF,'userdata',DATA);
init=-2;
case 'vox-thr',
temp=str2num(get(DATA.handles(2),'string'));
if ~isempty(temp), thres{1}=temp; else, set(DATA.handles(2),'string',num2str(thres{1})); end
thres{2}=get(DATA.handles(3),'value');
DATA.selectcluster=[];
DATA.thres=thres;set(GCF,'userdata',DATA);
init=-1;
case 'clu-thr',
temp=str2num(get(DATA.handles(5),'string'));
if ~isempty(temp), thres{3}=temp; else, set(DATA.handles(5),'string',num2str(thres{3})); end
thres{4}=get(DATA.handles(6),'value');
DATA.selectcluster=[];
if thres{4}>4, set(DATA.handles(4),'string','peak threshold: peak p =');
elseif thres{4}==1, set(DATA.handles(4),'string','cluster threshold: cluster k =');
else set(DATA.handles(4),'string','cluster threshold: cluster p <');
end
DATA.thres=thres;set(GCF,'userdata',DATA);
init=-1;
case 'export_mask',
if isempty(OPTION2)
[filename,filepath]=uiputfile('*.img;*.nii','Save mask as',fileparts(spmfile));
else
[filepath,filename_name,filename_ext]=fileparts(OPTION2);
filename=[filename_name,filename_ext];
end
if ischar(filename),
[nill,filename_name,filename_ext]=fileparts(filename);
if isempty(filename_ext), filename=[filename,'.img']; end
V=struct('mat',mat{1},'dim',size(b(:,:,:,1)),'dt',[spm_type('float32') spm_platform('bigend')],'fname',fullfile(filepath,filename));
V=spm_write_vol(V,abs(d).*double(b(:,:,:,end)>0));
fprintf('Mask file saved as %s\n',V.fname);
V=struct('mat',mat{1},'dim',size(b(:,:,:,1)),'pinfo',[1;0;0],'dt',[spm_type('uint32') spm_platform('bigend')],'fname',fullfile(filepath,[filename_name,'.ROIs.img']));
btemp=zeros(size(b(:,:,:,1)));
for nbtemp=1:numel(DATA.clusters), btemp(DATA.clusters{nbtemp})=nbtemp; end
V=spm_write_vol(V,btemp);
fprintf('ROI file saved as %s\n',V.fname);
fh=fopen(fullfile(filepath,[filename_name,'.ROIs.txt']),'wt');
for nbtemp=1:numel(DATA.clusters),
fprintf(fh,'%+d ',[DATA.xyzpeak{nbtemp},1]*DATA.mat{1}(1:3,:)');
fprintf(fh,'\n');
end
fclose(fh);
[peaks,peaks_idx,C]=conn_peaks({d,mat{1}});
peaks_suprathreshold=b(peaks_idx+size(b,1)*size(b,2)*size(b,3))>0;
%peaks_F=d(peaks_idx);
%peaks_p=exp(-b(peaks_idx(peaks_suprathreshold)));
%dof=mat{3};
save(fullfile(filepath,[filename_name,'.PEAKS.MAT']),'peaks','peaks_suprathreshold');%,'peaks_F','peaks_p','dof');
%fprintf('Peaks file saved as %s\n',fullfile(filepath,[filename_name,'.PEAKS.MAT']));
if ~isempty(C)
V=struct('mat',mat{1},'dim',size(C),'pinfo',[1;0;0],'dt',[spm_type('uint32') spm_platform('bigend')],'fname',fullfile(filepath,[filename_name,'.SEGs.img']));
V=spm_write_vol(V,C);
fprintf('Segmentation file saved as %s\n',V.fname);
fh=fopen(fullfile(filepath,[filename_name,'.SEGs.txt']),'wt');
for nbtemp=1:numel(peaks_idx),
fprintf(fh,'%+d ',peaks(nbtemp,:));
fprintf(fh,'\n');
end
fclose(fh);
end
end
return;
case 'export_list'
% Added function for export list button
% uicontextmenu on MacOSX seems not to fire
cbo = DATA.handles(8);
header = get(DATA.handles(7),'string');
conn_exportlist(cbo,'',header);
return;
case 'surface_view'
filename=fullfile(fileparts(spmfile),'results.img');
conn_vproject(gcf,[],'export_mask',filename);
conn_mesh_display(filename,'');
return;
case 'volume_view'
filename=fullfile(fileparts(spmfile),'results.img');
conn_vproject(gcf,[],'export_mask',filename);
conn_mesh_display('',filename,[],[],[],.5);
return;
case 'spm_view'
[spmfile_path,spmfile_name]=fileparts(spmfile);
cwd=pwd;cd(spmfile_path);
load(spmfile,'SPM');
spm defaults fmri;
[hReg,xSPM,SPM] = spm_results_ui('setup',SPM);
assignin('base','hReg',hReg);
assignin('base','xSPM',xSPM);
assignin('base','SPM',SPM);
figure(spm_figure('GetWin','Interactive'));
cd(cwd);
return;
case 'cluster_view'
filename=fullfile(fileparts(spmfile),'results.img');
conn_vproject(gcf,[],'export_mask',filename);
[spmfile_path,spmfile_name]=fileparts(spmfile);
cwd=pwd;cd(spmfile_path);
filename=fullfile(fileparts(spmfile),'results.ROIs.img');
conn_rex(spmfile,filename,'output_type','saverex','level','clusters','select_clusters',0,'steps',{'extract','results'},'s',[],'gui',1);
%roi_extract(spmfile,fullfile(filepath,filename),'output_type','text','select_clusters','each','steps',{'extract'});
cd(cwd);
return;
case 'cluster_import',
filename=fullfile(fileparts(spmfile),'results.img');
conn_vproject(gcf,[],'export_mask',filename);
[spmfile_path,spmfile_name]=fileparts(spmfile);
cwd=pwd;cd(spmfile_path);
filename=fullfile(fileparts(spmfile),'results.ROIs.img');
htfig=msgbox('Loading connectivity values. Please wait','');
[y,name]=conn_rex(spmfile,filename,'level','clusters');
if ishandle(htfig), delete(htfig); end
name=regexprep(name,'^results\.ROIs\.','cluster ');
temp=load(DATA.spmfile,'SPM');
if isfield(temp.SPM.xX,'SelectedSubjects')&&~rem(size(y,1),nnz(temp.SPM.xX.SelectedSubjects))
ty=zeros(size(y,1)/nnz(temp.SPM.xX.SelectedSubjects)*numel(temp.SPM.xX.SelectedSubjects),size(y,2));
ty(repmat(logical(temp.SPM.xX.SelectedSubjects),size(y,1)/nnz(temp.SPM.xX.SelectedSubjects),1),:)=y;
y=ty;
end
conn_importl2covariate(name,num2cell(y,1));
cd(cwd);
return;
case 'connectome',
optionDistancePeaks=12; % minimum distance between extracted peaks (mm)
switch(DATA.side)
case 1, [peaks,peaks_idx]=conn_peaks({d,mat{1}},optionDistancePeaks);
case 2, [peaks,peaks_idx]=conn_peaks({-d,mat{1}},optionDistancePeaks);
case 3, [peaks,peaks_idx]=conn_peaks({abs(d),mat{1}},optionDistancePeaks);
end
if 0 % one peak per cluster
peaks_C=zeros(size(d));
for n1=1:numel(DATA.clusters), peaks_C(DATA.clusters{n1})=n1; end
withinClusterPeak=zeros(numel(DATA.clusters)+1,1);
withinClusterPeak(1+peaks_C(flipud(peaks_idx)))=flipud(peaks_idx);
otherwithinClusterPeak=find(peaks_C(peaks_idx)>0&~ismember(peaks_idx,withinClusterPeak(2:end)));
peaks_idx(otherwithinClusterPeak)=[];
peaks(otherwithinClusterPeak,:)=[];
end
peaks_suprathreshold=b(peaks_idx+size(b,1)*size(b,2)*size(b,3))>0;
[spmfile_path,spmfile_name]=fileparts(spmfile);
save(fullfile(spmfile_path,'PEAKS.mat'),'peaks','peaks_suprathreshold');
conn_process('extract_connectome',peaks(peaks_suprathreshold,:),[peaks(peaks_suprathreshold,:);peaks(~peaks_suprathreshold,:)],-1);
return
% ROI=conn_process('results_connectome',spmfile_path,-1);
% save(fullfile(spmfile_path,'ROI.mat'),'ROI');
% conn_displayroi('initfile','results_roi',0,fullfile(spmfile_path,'ROI.mat'));
%conn_displayroi('initfile','results_connectome',0,spmfile_path,peaks_suprathreshold,fullfile(spmfile_path,'ROI.mat'));
%conn_displayroi('init','results_connectome',0,spmfile_path,peaks_suprathreshold);
case 'switch',
side=get(DATA.handles(11),'value');
if side~=DATA.side,
switch(side),
case 1, b=c;
case 2, b=1-c;
case 3, b=2*min(c,1-c);
end
%b(b==0)=nan;
b=-log(max(eps,b));
DATA.side=side;
%b=-log(1-exp(-DATA.b(:,:,:,1)));
%DATA.side=-DATA.side;
init=-1;
DATA.selectcluster=[];
end
otherwise,
return
end
end
else, %called not from UI
init=1;
if nargin<11 || isempty(mat), mat={eye(4),numel(a)}; end
if nargin<10 || isempty(select), select=[]; end
if nargin<9 || isempty(box), box=0; end
if nargin<8 || isempty(res), res=1; end
if nargin<7 || isempty(thres), thres={.05,1,1,1}; end
if nargin<6 || isempty(projection), projection=[-1,0,0;0,0,-1;0,-1,0]; end;
%if nargin<3 || isempty(projection), projection=[-0.3529,0.9347,-0.0429;0.0247,-0.0365,-0.9990;-0.9353,-0.3537,-0.0103]; end
if nargin<5 || isempty(views), views='full'; end
if nargin<4, d=[]; end
if nargin<3, c=[]; end
if nargin<2, b=[]; end
if nargin<1 || isempty(a), a=spm_read_vols(spm_vol(fullfile(fileparts(which('spm')),'apriori','grey.nii'))); end
if nargin<12 || isempty(threshold), if isempty(b), threshold=[nanmean(a(:))]; else, threshold=[nanmean(a(:)),2]; end; end
if nargin<13 || isempty(data1plot), data1plot=[]; end
if nargin<14 || isempty(spmfile), spmfile=[]; end
if nargin<15 || isempty(voxeltovoxel), voxeltovoxel=0; end
GCF=gcf;
if ~strcmp(views,'none'),
a(:,:,[1,end],1)=0;a(:,[1,end],:,1)=0;a([1,end],:,:,1)=0;
clf;
DATA.a=a;
DATA.b=b;
DATA.c=c;
DATA.d=d;
DATA.projection=projection;
DATA.thres=thres;
DATA.threshold=threshold;
DATA.res=res;
DATA.box=box;
DATA.views=views;
DATA.rotation=[0,0,0];
DATA.selectcluster=[];
DATA.select=select;
DATA.mat=mat;
DATA.clusters={[]};
DATA.side=1;
DATA.spmfile=spmfile;
DATA.voxeltovoxel=voxeltovoxel;
DATA.peakFDR=0;
end
end
if init>0, % initialize
%map=[gray(64).^2;.5*repmat([1,1,0],[64,1])+.5*((gray(64).^2)*diag([1,1,0]));.5*repmat([1,0,0],[64,1])+.5*((gray(64).^2)*diag([1,0,0]))];
%map=[.5+.5*(gray(64).^2);.5*repmat([1,1,0],[64,1])+.5*((gray(64).^.5)*diag([1,1,0]));.5*repmat([0,0,1],[64,1])+.5*((gray(64).^.5)*diag([0,0,1]))];map(1,:)=1;
map=[.25+.75*(gray(64).^2)*diag([1,1,1]);0*repmat([1,1,1],[64,1])+1*((gray(64).^1)*diag([1,0,0]));.25*repmat([1,1,0],[64,1])+.75*((gray(64).^.5)*diag([1,1,0]))];map(1,:)=.25;
switch(views),
case 'full',
clf;
if isempty(b),
DATA.handles=[uicontrol('style','text','units','norm','position',[.6,.95,.1,.05],'string','resolution','backgroundcolor','w','foregroundcolor','b'),...
uicontrol('style','edit','units','norm','position',[.7,.95,.1,.05],'string',num2str(DATA.res),'callback',{@conn_vproject,'resolution'},'backgroundcolor','w','foregroundcolor','k'),...
uicontrol('style','text','units','norm','position',[.8,.95,.1,.05],'string','threshold','backgroundcolor','w','foregroundcolor','b'),...
uicontrol('style','edit','units','norm','position',[.9,.95,.1,.05],'string',num2str(DATA.threshold),'callback',{@conn_vproject,'threshold'},'backgroundcolor','w','foregroundcolor','b')];
DATA.axes=subplot(122); set(DATA.axes,'units','norm','position',[.55,.1,.4,.8],'visible','off');
set(GCF,'name','Volume display','numbertitle','off','color','w','units','norm','position',[.1,.5,.8,.4],'interruptible','off','busyaction','cancel','windowbuttondownfcn',{@conn_vproject,'buttondown'},'windowbuttonupfcn',{@conn_vproject,'buttonup'},'windowbuttonmotionfcn',[],'keypressfcn',{@conn_vproject,'keypress'},'colormap',map,'userdata',DATA);
else
uicontrol('style','frame','units','norm','position',[.39,.87,.45,.11],'backgroundcolor',.9*[1,1,1]);
DATA.handles=[...
uicontrol('style','text','units','norm','position',[.40,.92,.15,.04],'string','height thr: voxel p <','fontname','arial','fontsize',8,'foregroundcolor','k','backgroundcolor',.9*[1,1,1]),...
uicontrol('style','edit','units','norm','position',[.55,.93,.05,.03],'string',num2str(thres{1}),'fontname','arial','fontsize',8,'callback',{@conn_vproject,'vox-thr'},'tooltipstring','Height-threshold value: False-positive threshold value for each voxel'),...
uicontrol('style','popupmenu','units','norm','position',[.625,.935,.1,.03],'string',{'p-uncorrected','p-FDR corrected'},'fontname','arial','fontsize',8,'callback',{@conn_vproject,'vox-thr'},'value',thres{2},'tooltipstring','False-positive control type for individual voxels'),...
uicontrol('style','text','units','norm','position',[.4,.88,.15,.04],'string','cluster thr: cluster p <','fontname','arial','fontsize',8,'foregroundcolor','k','backgroundcolor',.9*[1,1,1]),...
uicontrol('style','edit','units','norm','position',[.55,.89,.05,.03],'string',num2str(thres{3}),'fontname','arial','fontsize',8,'callback',{@conn_vproject,'clu-thr'},'tooltipstring','Extent threshold value: False-positive threshold value for individual clusters (based on cluster size)'),...
uicontrol('style','popupmenu','units','norm','position',[.625,.895,.1,.03],'string',{'cluster k-value','cluster p-FWE corrected','cluster p-FDR corrected','cluster p-uncorrected','peak p-FWE corrected','peak p-FDR corrected','peak p-uncorrected'},'fontname','arial','fontsize',8,'callback',{@conn_vproject,'clu-thr'},'value',thres{4},'tooltipstring','False-positive control type for individual clusters'),...
uicontrol('style','text','units','norm','position',[.05,.40,.925,.025],'string',sprintf('%18s%15s%15s%15s%15s%15s%15s%15s','Clusters (x,y,z)','k','cluster p-FWE','cluster p-FDR','cluster p-unc','peak p-FWE','peak p-FDR','peak p-unc'),'fontname','arial','fontsize',8,'backgroundcolor','w','foregroundcolor','b','horizontalalignment','left','fontname','monospaced','fontsize',8),...
uicontrol('style','listbox','units','norm','position',[.05,.06,.925,.34],'string','','backgroundcolor','w','foregroundcolor','k','horizontalalignment','left','fontname','monospaced','fontsize',8,'callback',{@conn_vproject,'selectroi'}),...
uicontrol('style','pushbutton','units','norm','position',[.875,.885,.1,.04],'string','Surface display','fontname','arial','fontsize',8,'callback',{@conn_vproject,'surface_view'},'tooltipstring','Displays results projected to cortical surface'),...
0,...
uicontrol('style','popupmenu','units','norm','position',[.875,.945,.1,.03],'string',{'positive contrast','negative contrast','two-sided'},'fontname','arial','fontsize',8,'callback',{@conn_vproject,'switch'},'tooltipstring','Analysis results directionality'),...
uicontrol('style','listbox','units','norm','position',[.62,.45,.355,.28],'string','','fontname','arial','fontsize',9,'visible','off','backgroundcolor','w','foregroundcolor','k','horizontalalignment','left','max',2),...
uicontrol('style','text','units','norm','position',[.725,.92,.1,.04],'string','','fontname','arial','fontsize',8,'horizontalalignment','right','foregroundcolor','k','backgroundcolor',.9*[1,1,1]),...
uicontrol('style','text','units','norm','position',[.725,.88,.1,.04],'string','','fontname','arial','fontsize',8,'horizontalalignment','right','foregroundcolor','k','backgroundcolor',.9*[1,1,1]),...
uicontrol('style','pushbutton','units','norm','position',[.875,.02,.1,.04],'string','Explore clusters','fontname','arial','fontsize',8,'callback',{@conn_vproject,'cluster_view'},'tooltipstring','Explores results within each significant cluster (using REX)'),...
uicontrol('style','pushbutton','units','norm','position',[.875,.765,.1,.04],'string','Export mask','fontname','arial','fontsize',8,'callback',{@conn_vproject,'export_mask'},'tooltipstring','Exports mask of supra-threshold voxels'),...
uicontrol('style','pushbutton','units','norm','position',[.875,.845,.1,.04],'string','Volume display','fontname','arial','fontsize',8,'callback',{@conn_vproject,'volume_view'},'tooltipstring','Displays results on 3d brain'),...
uicontrol('style','pushbutton','units','norm','position',[.775,.02,.1,.04],'string','Import values','fontname','arial','fontsize',8,'callback',{@conn_vproject,'cluster_import'},'tooltipstring','Imports average cluster connectivity values for each subject into CONN toolbox as second-level covariates')];
uicontrol('style','pushbutton','units','norm','position',[.875,.805,.1,.04],'string','SPM display','fontname','arial','fontsize',8,'callback',{@conn_vproject,'spm_view'},'tooltipstring','Displays results in SPM'),...
%uicontrol('style','pushbutton','units','norm','position',[.875,.845,.1,.04],'string','export stats','fontname','arial','fontsize',12,'callback',{@conn_vproject,'export_stats'}),...
% This sets up the export button
uicontrol('style','pushbutton','units','norm','position',[.875,.705,.1,.04],'string','Export list','fontname','arial','fontsize',12,'callback',{@conn_vproject,'export_list'}),... %%% LIZETTE ADDED THIS LINE!!
%uicontrol('style','text','units','norm','position',[.2,.30,.1,.025],'string','k','backgroundcolor','k','foregroundcolor','y','horizontalalignment','left'),...
%uicontrol('style','text','units','norm','position',[.3,.30,.1,.025],'string','voxel p-unc','backgroundcolor','k','foregroundcolor','y','horizontalalignment','left'),...
%uicontrol('style','text','units','norm','position',[.4,.30,.1,.025],'string','voxel p-cor','backgroundcolor','k','foregroundcolor','y','horizontalalignment','left'),...
%uicontrol('style','listbox','units','norm','position',[.2,.05,.1,.25],'string','','backgroundcolor','k','foregroundcolor','w','horizontalalignment','left'),...
%uicontrol('style','listbox','units','norm','position',[.3,.05,.1,.25],'string','','backgroundcolor','k','foregroundcolor','w','horizontalalignment','left'),...
%uicontrol('style','listbox','units','norm','position',[.4,.05,.1,.25],'string','','backgroundcolor','k','foregroundcolor','w','horizontalalignment','left')];
if DATA.mat{6}~='T', set(DATA.handles(11),'value',3,'enable','off'); end;
if ~isfield(DATA,'peakFDR')||~DATA.peakFDR,
set(DATA.handles(6),'string',{'cluster k-value','cluster p-FWE corrected','cluster p-FDR corrected','cluster p-uncorrected','peak p-FWE corrected','peak p-uncorrected'});
set(DATA.handles(7),'string',sprintf('%18s%15s%15s%15s%15s%15s%15s','Clusters (x,y,z)','k','cluster p-FWE','cluster p-FDR','cluster p-unc','peak p-FWE','peak p-unc'));
end
hc1=uicontextmenu;
uimenu(hc1,'Label','Export table','callback',@(varargin)conn_exportlist(DATA.handles(8),'',get(DATA.handles(7),'string')));
set(DATA.handles(8),'uicontextmenu',hc1);
hc1=uicontextmenu;
uimenu(hc1,'Label','Export table','callback',@(varargin)conn_exportlist(DATA.handles(12)));
set(DATA.handles(12),'uicontextmenu',hc1);
%uicontrol('style','pushbutton','units','norm','position',[.775,.02,.1,.04],'string','Export table','fontname','arial','fontsize',8,'callback',@(varargin)conn_exportlist(DATA.handles(12)),'tooltipstring','Exports table'),...
%hc1=uicontextmenu;
%uimenu(hc1,'Label','Export stats','callback',{@conn_vproject,'export_stats'});
%set(DATA.handles(8),'uicontextmenu',hc1);
if voxeltovoxel
DATA.handles(10)=uicontrol('style','pushbutton','units','norm','position',[.875,.845,.1,.04],'string','post-hoc peak analyses','fontname','arial','fontsize',8,'callback',{@conn_vproject,'connectome'});
end
%if DATA.mat{6}~='T'||isempty(DATA.mat{2}), %enable Fcluster%
% delete(DATA.handles(6)); DATA.handles(6)=uicontrol('style','popupmenu','units','norm','position',[.625,.895,.1,.03],'string',{'peak p-FDR corrected','extent k-value'},'fontname','arial','fontsize',8,'callback',{@conn_vproject,'clu-thr'},'value',min(2,thres{4}));
%end
%if isempty(DATA.mat{2}), set(DATA.handles(6),'visible','off','value',4); DATA.thres{4}=4; DATA.thres{3}=10; thres=DATA.thres; set(DATA.handles(5),'string',num2str(thres{3})); set(DATA.handles(4),'string','extent threshold: cluster k ='); end
DATA.axes=subplot(222);set(DATA.axes,'units','norm','position',[.40,.45,.2,.38],'visible','off');%[.55,.35,.4,.6],'visible','off');
h0=get(0,'screensize');
set(GCF,'name',['Volume display ',DATA.spmfile],'numbertitle','off','color','w','units','pixels','position',[2,h0(4)-min(540,.7*h0(3))-48,h0(3)-2,min(540,.7*h0(3))],'interruptible','off','busyaction','cancel','windowbuttondownfcn',{@conn_vproject,'buttondown'},'windowbuttonupfcn',{@conn_vproject,'buttonup'},'windowbuttonmotionfcn',[],'keypressfcn',{@conn_vproject,'keypress'},'colormap',map,'userdata',DATA);
end
case '3d',
clf;
DATA.handles=[uicontrol('style','text','units','norm','position',[0,0,.1,.05],'string','resolution','backgroundcolor','k','foregroundcolor','w'),...
uicontrol('style','edit','units','norm','position',[.1,0,.1,.05],'string',num2str(DATA.res),'callback',{@conn_vproject,'resolution'},'backgroundcolor','k','foregroundcolor','w'),...
uicontrol('style','text','units','norm','position',[.2,0,.1,.05],'string','threshold','backgroundcolor','k','foregroundcolor','w'),...
uicontrol('style','edit','units','norm','position',[.3,0,.1,.05],'string',num2str(DATA.threshold),'callback',{@conn_vproject,'threshold'},'backgroundcolor','k','foregroundcolor','w')];
DATA.axes=gca;set(DATA.axes,'visible','off');
set(GCF,'name','Volume display','numbertitle','off','color','w','interruptible','off','busyaction','cancel','windowbuttondownfcn',{@conn_vproject,'buttondown'},'windowbuttonupfcn',{@conn_vproject,'buttonup'},'windowbuttonmotionfcn',[],'keypressfcn',{@conn_vproject,'keypress'},'colormap',map,'userdata',DATA);
case 'orth',
clf;
DATA.handles=[uicontrol('style','text','units','norm','position',[0,0,.1,.05],'string','resolution','backgroundcolor','k','foregroundcolor','w'),...
uicontrol('style','edit','units','norm','position',[.1,0,.1,.05],'string',num2str(DATA.res),'callback',{@conn_vproject,'resolution'},'backgroundcolor','k','foregroundcolor','w'),...
uicontrol('style','text','units','norm','position',[.2,0,.1,.05],'string','threshold','backgroundcolor','k','foregroundcolor','w'),...
uicontrol('style','edit','units','norm','position',[.3,0,.1,.05],'string',num2str(DATA.threshold),'callback',{@conn_vproject,'threshold'},'backgroundcolor','k','foregroundcolor','w')];
DATA.axes=gca;set(DATA.axes,'visible','off');
%DATA.axes=gca;
set(GCF,'color','w','windowbuttondownfcn',[],'windowbuttonupfcn',[],'windowbuttonmotionfcn',[],'keypressfcn',[],'colormap',map,'userdata',DATA);
case 'none',
if ~nargout, DATA.axes=gca;set(DATA.axes,'visible','off'); end
%DATA.axes=gca;
%set(GCF,'color','w','userdata',DATA);
end
end
% display
if strcmp(views,'full'),
if init~=0, % redraws orth views
if ~isempty(b),%length(threshold)>1,
if isempty(DATA.clusters)||init~=-2, % compute/display stats
bnew=b(:,:,:,1);bnew(~a)=nan;%bnew(a<=threshold(1))=nan;
[bnew,txt,xyzpeak,clusters,clusternames,ithres]=vproject_display(bnew,threshold,thres,mat,DATA.side,DATA.d,DATA.peakFDR);
b(:,:,:,2)=bnew; threshold(3)=0;%.5;
%for n1=1:4, set(DATA.handles(8+n1),'string',strvcat(txt{n1})); end
set(DATA.handles(8),'string',strvcat(txt,' '),'value',max(1,min(size(txt,1)+1,get(DATA.handles(8),'value'))));
DATA.stdprojections=cell(1,4);
DATA.txt=txt;
DATA.clusternames=clusternames;
DATA.ithres=ithres;
else
clusters=DATA.clusters;
xyzpeak=DATA.xyzpeak;
clusternames=DATA.clusternames;
end
if isfield(DATA,'selectcluster'), selectcluster=DATA.selectcluster;else, selectcluster=get(DATA.handles(8),'value'); end
if selectcluster<=length(clusters), select=clusters{selectcluster}; clusternames=clusternames{selectcluster}.txt;
elseif (isempty(selectcluster)||selectcluster==length(clusters)+1)&&~isempty(clusternames), select=[]; clusternames=clusternames{length(clusters)+1}.txt;
else, select=[]; clusternames=[]; end
end
%M={[-1,0,0;0,0,-1;0,-1,0],[0,-1,0;0,0,-1;-1,0,0],[-1,0,0;0,1,0;0,0,1]};
M={[1,0,0;0,0,-1;0,1,0],[0,-1,0;0,0,-1;1,0,0],[1,0,0;0,-1,0;0,0,-1]};
set(GCF,'pointer','watch'); drawnow
if ~isfield(DATA,'stdprojections') || init==-1, DATA.stdprojections=cell(1,4); end
tplot={};infoplot={};for n1=1:3, [tplot{n1},infoplot{n1},DATA.stdprojections{n1}]=conn_vproject(a,b,c,d,'none',M{n1},thres,res,0,select,mat,threshold,DATA.stdprojections{n1}); end;
%tplot=[tplot{1},tplot{2};tplot{3},ones(size(tplot{3},1),size(tplot{2},2))];
tplot=cat(1,cat(2,tplot{1},tplot{2}),cat(2,tplot{3},tplot{1}(1)*ones([size(tplot{3},1),size(tplot{2},2),size(tplot{2},3)])));
if ~isempty(b), h=subplot(3,4,1); set(h,'units','norm','position',[.05,.45,.3,.5]); else, h=subplot(3,4,1); set(h,'units','norm','position',[.05,.1,.4,.8]);end
%pres=.5;image(1:pres:size(tplot,2),1:pres:size(tplot,1),round(max((ceil(tplot(1:pres:end,1:pres:end)/64)-1)*64+1,convn(convn(tplot(1:pres:end,1:pres:end),conn_hanning(7)/4,'same'),conn_hanning(7)'/4,'same'))));axis equal; axis off;
pres=.25;
conf=conn_hanning(2/pres+1);conf=conf/sum(conf);
temp=convn(convn(tplot(round(1:pres:end),round(1:pres:end),:),conf,'same'),conf','same');
image(1:pres:size(tplot,2),1:pres:size(tplot,1),temp);axis equal; axis off;
%image(round(convn(convn(tplot(round(1:.5:end),round(1:.5:end),:),conn_hanning(3)/2,'same'),conn_hanning(3)'/2,'same')));axis equal; axis off;
%image(tplot);axis equal; axis off;
data1plot=DATA.stdprojections{4};
if ~isempty(clusternames),set(DATA.handles(12),'string',clusternames,'value',1,'visible','on');else, set(DATA.handles(12),'string','','value',1,'visible','off');end
if numel(DATA.mat{3})>1, set(DATA.handles(13),'string',[DATA.mat{6},'(',num2str(DATA.mat{3}(1)),',',num2str(DATA.mat{3}(2)),')','_min = ',num2str(DATA.ithres(1),'%0.2f')]);
else set(DATA.handles(13),'string',[DATA.mat{6},'(',num2str(DATA.mat{3}),')','_min = ',num2str(DATA.ithres(1),'%0.2f')]);
end
set(DATA.handles(14),'string',['k_min =',num2str(DATA.ithres(2))]);
if ~isempty(b),
DATA.b=b;DATA.threshold=threshold;DATA.xyzpeak=xyzpeak;DATA.infoplot=infoplot;DATA.select=select;DATA.clusters=clusters;
%set(GCF,'userdata',DATA);
end
end
%subplot(122);conn_vproject(a,'none',projection,threshold,res,box);
%if ~isempty(b), DATA.axes=subplot(222); else, DATA.axes=subplot(122); end
set(DATA.axes,'visible','off');
%views='none';
set(GCF,'pointer','arrow');%,'userdata',DATA);
end
if strcmp(views,'orth'),
%M={[-1,0,0;0,0,-1;0,-1,0],[0,-1,0;0,0,-1;-1,0,0],[-1,0,0;0,1,0;0,0,1]};
M={[1,0,0;0,0,-1;0,1,0],[0,1,0;0,0,-1;1,0,0],[1,0,0;0,1,0;0,0,-1]};
set(GCF,'pointer','watch'); drawnow
for n1=1:3, subplot(2,2,n1);conn_vproject(a,b,c,d,'none',M{n1},threshold,res,box); end;
set(GCF,'pointer','arrow');
else,
if ~nargout, set(GCF,'pointer','watch'); drawnow; end
% volume plot
[dataplot,idx,mask,B,b2,res,data1plot]=vproject_view(a,b,c,projection,threshold,res,select,data1plot);
%if size(dataplot,3)>1, dataplot=dataplot(:,:,1)+dataplot(:,:,2); end
infoplot=struct('boundingbox',b2,'res',res,'size',[size(dataplot,1),size(dataplot,2)],'projection',projection);
if ~nargout,
axes(DATA.axes);
%image(dataplot);axis equal;axis off;
pres=.25;conf=conn_hanning(2/pres+1);conf=conf/sum(conf); image(1:pres:size(dataplot,2),1:pres:size(dataplot,1),convn(convn(dataplot(round(1:pres:end),round(1:pres:end),:),conf,'same'),conf','same'));axis equal; axis off;
%pres=.5;image(1:pres:size(dataplot,2),1:pres:size(dataplot,1),round(max((ceil(dataplot(1:pres:end,1:pres:end)/64)-1)*64+1,convn(convn(dataplot(1:pres:end,1:pres:end),conn_hanning(5)/3,'same'),conn_hanning(5)'/3,'same'))));axis equal; axis off;
DATA.axes=gca;
DATA.infoplot{4}=infoplot;DATA.stdprojections{4}=data1plot; set(GCF,'userdata',DATA);
% bounding box
if box,
B3=B'*pinv(projection);
B3=[1+(B3(:,1)-b2(1,1))/res, 1+(B3(:,2)-b2(1,2))/res, 1+(B3(:,3)-b2(1,3))/res];
border={[1,3,7,5,1],[2,6,4,8,2],[1,3,2,8,1],[7,5,4,6,7]};
[minBt3]=max(B3(:,3));
for n0=1:length(border),
for n1=1:length(border{n0})-1,
color=(B3(border{n0}(n1),3)==minBt3|B3(border{n0}(n1+1),3)==minBt3);
Bt1=linspace(B3(border{n0}(n1),1),B3(border{n0}(n1+1),1),100);
Bt2=linspace(B3(border{n0}(n1),2),B3(border{n0}(n1+1),2),100);
Bt3=linspace(B3(border{n0}(n1),3),B3(border{n0}(n1+1),3),100);
if color,
Bt1(Bt3>idx(max(1,min(prod(size(mask)), round(Bt2)+size(idx,1)*round(Bt1-1)))) & ...
mask(max(1,min(prod(size(mask)), round(Bt2)+size(idx,1)*round(Bt1-1)))) )=nan;
end
hold on; h=plot(Bt1,Bt2,'-'); set(h,'color',.25*[0,0,1],'linewidth',4-3*color);hold off;
end
end
end
set(GCF,'pointer','arrow');
end
end
end
function [dataplot,idx,mask,B,b2,res,data1plot]=vproject_view(a,b,c,projection,threshold,res,select,data1plot);
%interpolate&project
if ~isempty(data1plot),
mask=data1plot.mask;
idx=data1plot.idx;
trans=data1plot.trans;
B=data1plot.B;
b2=data1plot.b2;
else,
if ~isempty(b), a=a+j*b(:,:,:,end); end
[x1,y1,z1]=meshgrid((1:size(a,2))-(size(a,2)+1)/2,(1:size(a,1))-(size(a,1)+1)/2,(1:size(a,3))-(size(a,3)+1)/2);
bb={[1,size(a,2)]-(size(a,2)+1)/2,[1,size(a,1)]-(size(a,2)+1)/2,[1,size(a,3)]-(size(a,3)+1)/2};
B=[];for n1=1:8,B=[B,[bb{1}(1+rem(n1,2));bb{2}(1+rem(floor(n1/2),2));bb{3}(1+rem(floor(n1/4),2))]]; end
B2=B'*pinv(projection);
b2=[min(B2,[],1);max(B2,[],1)];
%res=res*((prod(b2(2,:)-b2(1,:))./prod(size(a))).^(1/3));
if 1,%faster&avoids memory problems
nz2=b2(1,3):res:b2(2,3);
[x2,y2]=meshgrid(b2(1,1):res:b2(2,1),b2(1,2):res:b2(2,2));
N=numel(x2);
d2=[x2(:),y2(:)]*projection(1:2,:);
d2=d2-repmat([bb{2}(1),bb{1}(1),bb{3}(1)],[N,1]);
d2=reshape(d2,[size(x2),3]);
trans=zeros(size(x2));
idx=zeros([size(x2),1+(length(threshold)>1)]);
idx1=1:N;idx2=[];
for n1=1:length(nz2),
x3=round(d2(idx1)+nz2(n1)*projection(3,1));
y3=round(d2(idx1+N)+nz2(n1)*projection(3,2));
z3=round(d2(idx1+N*2)+nz2(n1)*projection(3,3));
%a3=interp3(x1,y1,z1,a,x3,y3,z3,'nearest');
%idxs=max(1,min(size(a,1), 1+round(y3-bb{2}(1)) )) + size(a,1)*max(0,min(size(a,2)-1, round(x3-bb{1}(1)) )) + size(a,1)*size(a,2)*max(0,min(size(a,3)-1, round(z3-bb{3}(1)) ));
idxs=max(1,min(size(a,1), 1+y3 )) + size(a,1)*max(0,min(size(a,2)-1, x3 )) + size(a,1)*size(a,2)*max(0,min(size(a,3)-1, z3 ));
a3=a(idxs);
idx0=real(a3)>threshold(1)|imag(a3)>threshold(end);
idx(idx1(idx0))=(length(nz2)+1-n1)/length(nz2);
idx2=[idx2,idx1(idx0)];
idx1=idx1(~idx0);
if length(threshold)>1,
x3=d2(idx2)+nz2(n1)*projection(3,1);
y3=d2(idx2+N)+nz2(n1)*projection(3,2);
z3=d2(idx2+N*2)+nz2(n1)*projection(3,3);
%a3=interp3(x1,y1,z1,a,x3,y3,z3,'nearest');
%idxt=max(1,min(size(a,1), 1+round(y3-bb{2}(1)) )) + size(a,1)*max(0,min(size(a,2)-1, round(x3-bb{1}(1)) )) + size(a,1)*size(a,2)*max(0,min(size(a,3)-1, round(z3-bb{3}(1)) ));
idxt=max(1,min(size(a,1), 1+round(y3) )) + size(a,1)*max(0,min(size(a,2)-1, round(x3) )) + size(a,1)*size(a,2)*max(0,min(size(a,3)-1, round(z3) ));
a3=a(idxt);
if 0,
idx00=find(imag(a3)>threshold(end));
idx(N+idx2(idx00))=max(idx(N+idx2(idx00)),imag(a3(idx00))); %n1;
idxtrans=find(trans(idx2(idx00))==0);
trans(idx2(idx00(idxtrans)))=idxt(idx00(idxtrans));
else,
idx00=find(imag(a3)>threshold(end));
idxtrans=find(~idx(N+idx2(idx00)));%idxtrans=find(imag(a3(idx00))>=idx(N+idx2(idx00)));
idx(N+idx2(idx00(idxtrans)))=(length(nz2)+1-n1)/length(nz2);%imag(a3(idx00(idxtrans)));
trans(idx2(idx00(idxtrans)))=idxt(idx00(idxtrans));
end
%idx2=idx2(~idx00);
end
end
mask=(idx>0);
%if any(any(mask(:,:,1)>0)), idx(~mask(:,:,1))=5*max(idx(mask(:,:,1)>0)); else, idx(~mask(:,:,1))=length(nz2); end
else,
[x2,y2,z2]=meshgrid(b2(1,1):res:b2(2,1),b2(1,2):res:b2(2,2),b2(1,3):res:b2(2,3));
d=[x2(:),y2(:),z2(:)]*projection;
x2=reshape(d(:,1),size(x2));
y2=reshape(d(:,2),size(y2));
z2=reshape(d(:,3),size(z2));
%a2=interp3(x1,y1,z1,a,x2,y2,z2,'nearest');
a2=a(max(1,min(size(a,1), 1+round(y2-bb{2}(1)) )) + size(a,1)*max(0,min(size(a,2)-1, round(x2-bb{1}(1)) )) + size(a,1)*size(a,2)*max(0,min(size(a,3)-1, round(z2-bb{3}(1)) )));
[mask,idx]=max(a2>threshold,[],3);idx(~mask)=max(idx(mask>0));
end
data1plot.mask=mask;
data1plot.idx=idx;
data1plot.trans=trans;
data1plot.B=B;
data1plot.b2=b2;
end
conf1=conn_hanning(3)/2;conf1=conf1/sum(conf1);
conf2=conn_hanning(3)/2;conf2=conf2/sum(conf2);
if size(idx,3)==1, dataplot=-idx; dataplot=1+63*(dataplot-min(dataplot(:)))/(max(dataplot(:))-min(dataplot(:)));
else,
%dataplot1=-idx(:,:,1); dataplot1=1+63*(dataplot1-min(dataplot1(:)))/(max(dataplot1(:))-min(dataplot1(:)));
%conjmask=mask(:,:,1)&mask(:,:,2);
%dataplot2=-idx(:,:,2); dataplot2=1+63*(dataplot2-min(dataplot2(:)))/(max(dataplot2(:))-min(dataplot2(:)));
%dataplot=dataplot1; dataplot(conjmask)=64+dataplot2(conjmask);
if 1,
dataplot1=idx(:,:,1);
dataplot2=idx(:,:,2);
%dataplot1=convn(convn(dataplot1,conf1,'same'),conf1','same');
%dataplot2=convn(convn(dataplot2,conf2,'same'),conf2','same');
k=.75;%.2;
dataplotA=(dataplot2==0).*(k+(1-k)*dataplot1.^4) + (dataplot2>0).*(k+(1-k)*dataplot1.^4);
dataplotA(dataplot1==0)=1;
dataplotB=(dataplot2==0).*dataplotA + (dataplot2>0).*(dataplotA.*max(0,min(.75,tanh(1.5*(1*dataplot1-dataplot2)))));
%dataplotB=(dataplot2==0).*dataplotA + (dataplot2>0).*(.25*dataplotA.*(dataplot2<=.9*dataplot1));
%dataplotB=(dataplot2==0).*dataplotA + (dataplot2>0).*(0*(dataplot2>.95*dataplot1) + .5*dataplotA.*(dataplot2<=.95*dataplot1));
dataplot=cat(3,dataplotA,dataplotB,dataplotB);
else,
dataplot1=idx(:,:,1);
dataplot1=0+1*(dataplot1-min(dataplot1(~isinf(dataplot1))))/(max(dataplot1(~isinf(dataplot1)))-min(dataplot1(~isinf(dataplot1))));
dataplot2=idx(:,:,2); if any(dataplot2(:)>0),mindataplot=.95*min(dataplot2(dataplot2>0));else,mindataplot=0;end;
%dataplot2=dataplot2>mindataplot;
dataplot2=max(0,dataplot2-mindataplot);
dataplot2=convn(convn(dataplot2,conf1,'same'),conf1','same');
dataplot2=0+1*dataplot2/max(eps,max(dataplot2(:)));
%dataplot1=max(prctile(dataplot1(dataplot1<max(dataplot1(:))),5),dataplot1);dataplot2=max(0,dataplot1-dataplot2);
%dataplot1=.75+.25*(dataplot1).^2;dataplot2=max(0,dataplot1-dataplot2);
dataplot1=convn(convn(max(0,min(1,dataplot1)),conf2,'same'),conf2','same');
dataplot1=0+1*(dataplot1).^4;
%dataplotA=max(0,min(1, dataplot1+.25*(dataplot2>0)));dataplotB=max(0,min(1, dataplot1-dataplot2+.0*(dataplot2>0)));
dataplotA=max(0,min(1, dataplot1.*(dataplot2==0)+(.25+dataplot1).*(dataplot2>0)));
dataplotB=max(0,min(1, dataplot1.*(dataplot2==0)+(.25+dataplot1).*(.5-.5*dataplot2).*(dataplot2>0)));
dataplot=cat(3,dataplotA,dataplotB,dataplotB);
%idxback=find(all(dataplot<.1,3)); dataplot(idxback)=1;dataplot(idxback+size(dataplot,1)*size(dataplot,2))=1;dataplot(idxback+2*size(dataplot,1)*size(dataplot,2))=1;
%dataplot=cat(3,dataplot1,dataplot2,dataplot2);
end
end
if ~isempty(select),
[transa,transb,transc]=unique(trans(:));%trans=a(c);
[nill,idxp]=intersect(transa(:),select(:));
d=logical(zeros(size(transa)));d(idxp)=1;
e=find(d(transc));
%%dataplot(e)=dataplot(e)+64;
dataplot3=zeros([size(dataplot,1),size(dataplot,2)]);dataplot3(e)=1;dataplot3=convn(convn(dataplot3,conf1,'same'),conf1,'same');e=find(dataplot3>0);
n=size(dataplot,1)*size(dataplot,2);
dataplot(e+1*n)=dataplot(e+1*n).*(1-dataplot3(e))+0*dataplot(e+0*n).*(dataplot3(e));
dataplot(e+2*n)=dataplot(e+2*n).*(1-dataplot3(e))+0*dataplot(e+0*n).*(dataplot3(e));
dataplot(e+0*n)=dataplot(e+0*n).*(1-dataplot3(e))+0*dataplot(e+0*n).*(dataplot3(e));
end
if size(idx,3)~=1,
n=size(dataplot,1)*size(dataplot,2);
idxtemp1=find(dataplot2>0);
idxtemp2=find(c(trans(idxtemp1))>.5);
idxtemp=idxtemp1(idxtemp2);
temp3=dataplot(idxtemp+2*n);
dataplot(idxtemp+2*n)=dataplot(idxtemp+0*n);
dataplot(idxtemp+0*n)=temp3;
end
end
function [anew,txt,xyzpeak,clusters,clusternames,ithres]=vproject_display(a,threshold,thres,mat,side,Z,peakfdr)
global CONN_gui;
persistent refsrois;
if isempty(refsrois),
if isfield(CONN_gui,'refs')&&isfield(CONN_gui.refs,'rois')&&isfield(CONN_gui.refs.rois,'filename')&&~isempty(CONN_gui.refs.rois.filename),
refsrois=CONN_gui.refs.rois;
else
%filename=fullfile(fileparts(which('conn')),'utils','otherrois','Td.img');
%filename=fullfile(fileparts(which('conn')),'rois','aal.nii');
filename=fullfile(fileparts(which('conn')),'rois','BA.img');
[filename_path,filename_name,filename_ext]=fileparts(filename);
V=spm_vol(filename);
refsrois=struct('filename',filename,'filenameshort',filename_name,'V',V,'data',spm_read_vols(V),'labels',{textread(fullfile(filename_path,[filename_name,'.txt']),'%s','delimiter','\n')});
end
[x,y,z]=ndgrid(1:size(a,1),1:size(a,2),1:size(a,3));xyz=[x(:),y(:),z(:)];
refsrois.data=spm_get_data(refsrois.V,pinv(refsrois.V.mat)*mat{1}*[xyz,ones(size(xyz,1),1)]');
maxrefsrois=max(refsrois.data(:)); refsrois.count=hist(refsrois.data(:),0:maxrefsrois);
refsrois.labels={'not-labeled',refsrois.labels{:}};
end
if nargin<5, side=1; end
if nargin<4, mat={eye(4),numel(a)}; end
idx0=find(~isnan(a));
p=exp(-a(idx0));
p0=a;p0(idx0)=p;
P=a;P(idx0)=conn_fdr(p(:));
ithres=[nan,nan];
switch(thres{2}),
case 1,%'vox-unc',
% voxels above height threshold
idxvoxels=find(a>-log(thres{1}));
anew=zeros(size(a));anew(idxvoxels)=a(idxvoxels)+log(thres{1});%1;
if ~isempty(idxvoxels),
[xt,yt,zt]=ind2sub(size(a),idxvoxels);C=spm_clusters([xt,yt,zt]');c=hist(C,1:max(C));
end
if mat{6}=='T'&&side<3, %T stat
u=spm_invTcdf(1-thres{1},mat{3}(2));
ithres(1)=u;
elseif mat{6}=='T'&&side==3, %two-sided T stat
u=spm_invFcdf(1-thres{1},mat{3}(1),mat{3}(2));
ithres(1)=sqrt(u);
else %F stat
u=spm_invFcdf(1-thres{1},mat{3}(1),mat{3}(2));
ithres(1)=u;
end
case 2,%'vox-FDR',
idxvoxels=find(P<thres{1}) ;
anew=zeros(size(a));anew(idxvoxels)=-log(P(idxvoxels))+log(thres{1});%1;
if ~isempty(idxvoxels),
[xt,yt,zt]=ind2sub(size(a),idxvoxels);C=spm_clusters([xt,yt,zt]');c=hist(C,1:max(C)); % C: cluster per voxel; c: #voxels per clusters
end
if isempty(idxvoxels), u=1; ithres(1)=inf;else,
if mat{6}=='T'&&side<3, %T stat
u=spm_invTcdf(1-exp(-min(a(idxvoxels))),mat{3}(2)); ithres(1)=u;
elseif mat{6}=='T'&&side==3, %two-sided T stat
u=spm_invFcdf(1-exp(-min(a(idxvoxels))),mat{3}(1),mat{3}(2)); ithres(1)=sqrt(u);
else %F stat
u=spm_invFcdf(1-exp(-min(a(idxvoxels))),mat{3}(1),mat{3}(2)); ithres(1)=u;
end
end
end
% if mat{6}=='T'&&side>2, %T stat two-sided
% ithres(1)=spm_invTcdf(1-(1-spm_Tcdf(ithres(1),mat{3}(2)))/2,mat{3}(2));
% end
txt=[];xyzpeak={};clusters={};clusternames={};
if ~isempty(idxvoxels),
xyz=[xt,yt,zt];
idx1=find(c>0);
idxn=[];idx2={};for n1=1:length(idx1),idx2{n1}=find(C==(idx1(n1))); idxn(n1)=length(idx2{n1}); end
[nill,sidxn]=sort(idxn(:));sidxn=flipud(sidxn);
xyzpeak=cell(1,length(idx1));
clusters=cell(1,length(idx1));
k=zeros(1,length(idx1));cp=k;cP=k;Ez=k;cPFDR=k; minP=k;minp=k;maxZ=k;pPFWE=k;
for n1=1:length(idx1),
k(n1)=idxn(sidxn(n1));
clusters{n1}=idxvoxels(idx2{sidxn(n1)});
[minp(n1),idxminp]=min(p0(clusters{n1}));
minP(n1)=P(clusters{n1}(idxminp));
maxZ(n1)=Z(clusters{n1}(idxminp));
if side==2, maxZ(n1)=-maxZ(n1); end
%[minP(n1),idxminp]=min(P(clusters{n1}));
%minp(n1)=exp(-max(a(clusters{n1})));
xyzpeak{n1}=xyz(idx2{sidxn(n1)}(idxminp),:);
if ~isempty(mat{2}),
if mat{6}=='T'&&side<3, %T stat
if isempty(mat{5}), [cP(n1),cp(n1)]=spm_P(1,k(n1)*mat{2}(end)/mat{4},u,mat{3},'T',mat{2},1,mat{4});
else, [cP(n1),cp(n1)]=spm_P(1,k(n1)*mat{5},u,mat{3},'T',mat{2},1,mat{4});
if isnan(cP(n1))&&~cp(n1),cP(n1)=0; end
end
pPFWE(n1)=spm_P(1,0,maxZ(n1),mat{3},'T',mat{2},1,mat{4});
elseif mat{6}=='T'&&side==3, %two-sided T stat
if isempty(mat{5}), [nill,cP(n1),nill,cp(n1)] =stat_thres(mat{2},mat{4},mat{7},[mat{3};inf,inf], .05, u, k(n1)*mat{2}(end)/mat{4});
else, [nill,cP(n1),nill,cp(n1)] =stat_thres(mat{2},mat{4},mat{7},[mat{3};inf,inf], .05, u, k(n1)*mat{5});
if isnan(cP(n1))&&~cp(n1),cP(n1)=0; end
end
pPFWE(n1)=spm_P(1,0,maxZ(n1).^2,mat{3},'F',mat{2},1,mat{4});
else %F stat
if 0,%enable Fcluster% %note:change this to 0 for NS toolbox implementation of cluster-level stats for F maps (assuming homogeneity)
cp(n1)=nan; cP(n1)=nan;
else
if isempty(mat{5}), [nill,cP(n1),nill,cp(n1)] =stat_thres(mat{2},mat{4},mat{7},[mat{3};inf,inf], .05, u, k(n1)*mat{2}(end)/mat{4});
else, [nill,cP(n1),nill,cp(n1)] =stat_thres(mat{2},mat{4},mat{7},[mat{3};inf,inf], .05, u, k(n1)*mat{5});
if isnan(cP(n1))&&~cp(n1),cP(n1)=0; end
end
end
pPFWE(n1)=spm_P(1,0,maxZ(n1),mat{3},'F',mat{2},1,mat{4});
end
else, cP(n1)=nan;cp(n1)=nan; pPFWE(n1)=nan; end
end
if peakfdr>0,
try
if peakfdr==1, % consider all peak voxels for peak-FDR correction
if side==2, tZ=-Z(idxvoxels); else tZ=Z(idxvoxels); end
mtZ=min(tZ);
[maxN,maxZ,maxXYZ,maxA]=spm_max(1-mtZ+tZ,xyz');
maxZ=maxZ+mtZ-1;
[maxOK,maxI]=min(sum(abs(conn_bsxfun(@minus,permute(maxXYZ',[1,3,2]),permute(cell2mat(xyzpeak'),[3,1,2]))).^2,3),[],1);
else % consider only onepeak voxel within each cluster for peak-FDR correction
maxI=1:numel(maxZ);
end
Ez=nan(size(maxZ));
if mat{6}=='T'&&side<3,
[nill,nill,Eu] = spm_P_RF(1,0,u,mat{3},'T',mat{2},1);
for n1=1:numel(maxZ)
[nill,nill,Ez(n1)] = spm_P_RF(1,0,maxZ(n1),mat{3},'T',mat{2},1);
end
elseif mat{6}=='T'&&side==3,
[nill,nill,Eu] = spm_P_RF(1,0,u,mat{3},'F',mat{2},1);
for n1=1:numel(maxZ)
[nill,nill,Ez(n1)] = spm_P_RF(1,0,maxZ(n1).^2,mat{3},'F',mat{2},1);
end
else
[nill,nill,Eu] = spm_P_RF(1,0,u,mat{3},'F',mat{2},1);
for n1=1:numel(maxZ)
[nill,nill,Ez(n1)] = spm_P_RF(1,0,maxZ(n1),mat{3},'F',mat{2},1);
end
end
pPFDR=conn_fdr(Ez/Eu);
pPFDR=pPFDR(maxI);
maxZ=maxZ(maxI);
catch
pPFDR=nan(size(k));
end
end
cPFDR=conn_fdr(cp);
for n1=1:length(idx1),
temp=['( ',sprintf('%+03.0f ',(mat{1}(1:3,:)*[xyzpeak{n1}';1])'),') '];
temp=[temp,repmat(' ',[1,18-length(temp)])];
if peakfdr>0,
txt=strvcat(txt,[...
temp,...
[sprintf('%15d',k(n1))],...
[sprintf('%15f',cP(n1))],...
[sprintf('%15f',cPFDR(n1))],...
[sprintf('%15f',cp(n1))],...
[sprintf('%15f',pPFWE(n1))],...
[sprintf('%15f',pPFDR(n1))],...
[sprintf('%15f',minp(n1))]]);
%[sprintf('%15f',minP(n1))]]);
else
txt=strvcat(txt,[...
temp,...
[sprintf('%15d',k(n1))],...
[sprintf('%15f',cP(n1))],...
[sprintf('%15f',cPFDR(n1))],...
[sprintf('%15f',cp(n1))],...
[sprintf('%15f',pPFWE(n1))],...
[sprintf('%15f',minp(n1))]]);
end
end
if ~peakfdr, thres4=thres{4}+(thres{4}>=6);
else thres4=thres{4};
end
switch(thres4),
case 1,%'k-value',
idxclu=find(k>=thres{3});
idxrem=find(~(k>=thres{3}));
case 2,%'clu-FWE',
idxclu=find(cP<thres{3});
idxrem=find(~(cP<thres{3}));
case 3,%'clu-FDR',
idxclu=find(cPFDR<thres{3});
idxrem=find(~(cPFDR<thres{3}));
case 4,%'clu-unc',
idxclu=find(cp<thres{3});
idxrem=find(~(cp<thres{3}));
case 5,%'peak-FWE',
idxclu=find(pPFWE<thres{3});
idxrem=find(~(pPFWE<thres{3}));
case 6,%'peak-FDR',
idxclu=find(pPFDR<thres{3});
idxrem=find(~(pPFDR<thres{3}));
case 7,%'peak-unc',
idxclu=find(minp<thres{3});
idxrem=find(~(minp<thres{3}));
end
if ~isempty(idxclu), ithres(2)=min(k(idxclu)); end
for n1=1:length(idxrem), anew(clusters{idxrem(n1)})=0; end
txt=txt(idxclu,:);
xyzpeak={xyzpeak{idxclu}};
clusters={clusters{idxclu}};
if ~isempty(idxclu),
for n1=1:length(idxclu)+1,
if n1<=length(idxclu),xyztemp=xyz(idx2{sidxn(idxclu(n1))},:);
else, xyztemp=xyz(cat(2,idx2{sidxn(idxclu)}),:); end;
v=spm_get_data(refsrois.V,pinv(refsrois.V.mat)*mat{1}*[xyztemp,ones(size(xyztemp,1),1)]');
uv=unique(v);
clusternames{n1}.uvb=zeros(1,length(uv));for n2=1:length(uv),clusternames{n1}.uvb(n2)=sum(v==uv(n2));end
clusternames{n1}.uvc=refsrois.count(1+uv);
clusternames{n1}.uvd={refsrois.labels{1+uv}};
[nill,uvidx]=sort(-clusternames{n1}.uvb+1e10*(uv==0));
clusternames{n1}.uvb=clusternames{n1}.uvb(uvidx);clusternames{n1}.uvc=clusternames{n1}.uvc(uvidx);clusternames{n1}.uvd={clusternames{n1}.uvd{uvidx}};
clusternames{n1}.txt=cell(1,length(uv)); for n2=1:length(uv),clusternames{n1}.txt{n2}=[num2str(clusternames{n1}.uvb(n2)),' voxels covering ',num2str(clusternames{n1}.uvb(n2)/clusternames{n1}.uvc(n2)*100,'%0.0f'),'% of ',refsrois.filenameshort,'.',clusternames{n1}.uvd{n2}]; end
clusternames{n1}.txt=strvcat(clusternames{n1}.txt{:});
end
end
% clusternames={clusternames{idxclu,1},clusternames{idxclu(idxclulast),2}};
%for n1=1:length(clusternames),disp(' ');disp(clusternames{n1}.txt);end
end
%if params.plotlist, figure('color','w'); end
%h=[];xlim=[];ylim=[];
%disp(txt);
% if params.permutation && ~isempty(params.L),
% txt{end}=[txt{end},' p=',num2str(mean(params.L(:,1)>=k),'%0.4f'),' (',num2str(length(unique(params.L(params.L(:,1)>=k,2)))/params.Ln,'%0.4f'),')'];
% end
% if params.plotlist,
% %subplot(length(idx1),1,n1);
% h(n1)=axes('units','norm','position',...
% [.9,.5*(1-n1/max(length(idx1)+1,5)+.05/max(length(idx1)+1,5)),.09,.5*.9/max(length(idx1)+1,5)]);
% mx1=mean(X1(:,idx(idx2{sidxn(n1)})),2);
% mx2=mean(X2(:,idx(idx2{sidxn(n1)})),2);
% plot(mx1,mx2,'.'); h0=ylabel(txt{end});set(h0,'fontsize',8,'rotation',0,'horizontalalignment','right')
% params.output.clusters{n1}=[mx1,mx2];
% xlim=cat(1,xlim,[min(mx1),max(mx1)]);ylim=cat(1,ylim,[min(mx2),max(mx2)]);
% end
% end
%if params.plotlist,
% for n1=1:length(idx1),
% set(h(n1),'xlim',[min(xlim(:,1))-.01,max(xlim(:,2))+.01],'ylim',[min(ylim(:,1))-.01,max(ylim(:,2))+.01]);
% if n1~=length(idx1), set(h(n1),'xticklabel',[],'yticklabel',[],'fontsize',6); end
% end
%end
end
function [peak_threshold, extent_threshold, peak_threshold_1, extent_threshold_1] = ...
stat_thres(search_volume, num_voxels, fwhm, df, p_val_peak, ...
cluster_threshold, p_val_extent, nconj, nvar, EC_file, expr)
%
% stat_thresh.m
% Modified version of STAT_THRESHOLD function by Keith Worsley. The original
% STAT_THRESHOLD function is part of the FMRISTAT packgage. The main use of the
% original version is to calculate peak and cluster thresholds, both corrected
% and uncorrected. Details on input and output arguments are found in the
% original STAT_THRESHOLD function available from Keith Worsley's web site
% at McGill University's Math & Stat department.
%
% This stat_thresh.m function is a customized version to be called by a function
% spm_list_nS for non-stationarity correction of RFT-based cluster size test.
% The input and output of this function is therefore modified for producing cluster
% p-values (corrected) under non-stationarity. The modification includes:
% -supressing the output from being displayed.
% -the number of cluster p-values it can calculate has been increased to 500 clusters
% (the default in the original version was 5).
% -the p_val_extent is treated as extent, no matter how small it is.
%
% stat_thresh is called by spm_list_nS in the following format:
% [PEAK_P CLUSTER_P PEAK_P_1 CLUSTER_P_1] =
% stat_thresh(V_RESEL,NUM_VOX,1,[DF_ER;DF_RPV],ALPHA,CL_DEF_TH,CL_RESEL);
% PARAMETERS:
% V_RESEL: The seach volume in terms of resels. It is a 1x4 vector
% describing the topological characteristic of the search
% volume.
% NUM_VOX: The number of voxels in the search volume.
% DF_ER: Degrees of freedom of error
% DF_RPV: Degrees of freedom of RPV image estimation. Usually the same
% as the error df.
% ALPHA: The significance level of the peak (arbitrarily set to 0.05)
% CL_DEF_TH: The cluster defining threshold. Can be entered in terms of
% a p-value (uncorrected) or a t-value.
% CL_RESEL: The cluster size in terms of resel
%
% PEAK_P: Peak p-value (FWE corrected). Not used for our purpose.
% PEAK_P_1: Peak p-value (uncorrected). Not used for our purpose.
% CLUSTER_P: Cluster p-value (FWE corrected)
% CLUSTER_P_1: Cluster p-value (uncorrected)
%
% ----------------
%
% More etails on non-stationary cluster size test can be found in
%
% Worsley K J, Andermann M, Koulis T, MacDonald D and Evans A C
% Detecting Changes in Nonisotropic Images
% Human Brain Mapping 8: 98-101 (1999)
%
% Hayasaka S, Phan K L, Liberzon I, Worsley K J, and Nichols T E
% Nonstationary cluster size inference with random-field and permutation methods
% NeuroImage 22: 676-687 (2004)
%
%
%-----------------------------------------------------------------------------------
% Version 0.76b Feb 19, 2007 by Satoru Hayasaka
%
%############################################################################
% COPYRIGHT: Copyright 2003 K.J. Worsley
% Department of Mathematics and Statistics,
% McConnell Brain Imaging Center,
% Montreal Neurological Institute,
% McGill University, Montreal, Quebec, Canada.
% [email protected] , www.math.mcgill.ca/keith
%
% Permission to use, copy, modify, and distribute this
% software and its documentation for any purpose and without
% fee is hereby granted, provided that this copyright
% notice appears in all copies. The author and McGill University
% make no representations about the suitability of this
% software for any purpose. It is provided "as is" without
% express or implied warranty.
%############################################################################
%############################################################################
% UPDATES:
%
% Variable nvar is rounded so that it is recognized as an integer.
% Feb 19, 2007 by Satoru Hayasaka
%############################################################################
% Defaults:
if nargin<1; search_volume=[]; end
if nargin<2; num_voxels=[]; end
if nargin<3; fwhm=[]; end
if nargin<4; df=[]; end
if nargin<5; p_val_peak=[]; end
if nargin<6; cluster_threshold=[]; end
if nargin<7; p_val_extent=[]; end
if nargin<8; nconj=[]; end
if nargin<9; nvar=[]; end
if nargin<10; EC_file=[]; end
if nargin<11; expr=[]; end
if isempty(search_volume); search_volume=1000000; end
if isempty(num_voxels); num_voxels=1000000; end
if isempty(fwhm); fwhm=0.0; end
if isempty(df); df=Inf; end
if isempty(p_val_peak); p_val_peak=0.05; end
if isempty(cluster_threshold); cluster_threshold=0.001; end
if isempty(p_val_extent); p_val_extent=0.05; end
if isempty(nconj); nconj=1; end
if isempty(nvar); nvar=1; end
if size(fwhm,1)==1; fwhm(2,:)=fwhm; end
if size(fwhm,2)==1; scale=1; else scale=fwhm(1,2)/fwhm(1,1); fwhm=fwhm(:,1); end;
isscale=(scale>1);
if length(num_voxels)==1; num_voxels(2,1)=1; end
if size(search_volume,2)==1
radius=(search_volume/(4/3*pi)).^(1/3);
search_volume=[ones(length(radius),1) 4*radius 2*pi*radius.^2 search_volume];
end
if size(search_volume,1)==1
search_volume=[search_volume; [1 zeros(1,size(search_volume,2)-1)]];
end
lsv=size(search_volume,2);
fwhm_inv=all(fwhm>0)./(fwhm+any(fwhm<=0));
resels=search_volume.*repmat(fwhm_inv,1,lsv).^repmat(0:lsv-1,2,1);
invol=resels.*(4*log(2)).^(repmat(0:lsv-1,2,1)/2);
for k=1:2
D(k,1)=max(find(invol(k,:)))-1;
end
% determines which method was used to estimate fwhm (see fmrilm or multistat):
df_limit=4;
% max number of pvalues or thresholds to print:
% it can print out a ton of stuff! (the original default was 5)
nprint=500;
if length(df)==1; df=[df 0]; end
if size(df,1)==1; df=[df; Inf Inf]; end
if size(df,2)==1; df=[df [0; df(2,1)]]; end
% is_tstat=1 if it is a t statistic
is_tstat=(df(1,2)==0);
if is_tstat
df1=1;
df2=df(1,1);
else
df1=df(1,1);
df2=df(1,2);
end
if df2 >= 1000; df2=Inf; end
df0=df1+df2;
dfw1=df(2,1);
dfw2=df(2,2);
if dfw1 >= 1000; dfw1=Inf; end
if dfw2 >= 1000; dfw2=Inf; end
if length(nvar)==1; nvar(2,1)=df1; end
nvar = round(nvar); %-to make sure that nvar is integer!
if isscale & (D(2)>1 | nvar(1,1)>1 | df2<Inf)
D
nvar
df2
fprintf('Cannot do scale space.');
return
end
Dlim=D+[scale>1; 0];
DD=Dlim+nvar-1;
% Values of the F statistic:
t=((1000:-1:1)'/100).^4;
% Find the upper tail probs cumulating the F density using Simpson's rule:
if df2==Inf
u=df1*t;
b=exp(-u/2-log(2*pi)/2+log(u)/4)*df1^(1/4)*4/100;
else
u=df1*t/df2;
b=exp(-df0/2*log(1+u)+log(u)/4-betaln(1/2,(df0-1)/2))*(df1/df2)^(1/4)*4/100;
end
t=[t; 0];
b=[b; 0];
n=length(t);
sb=cumsum(b);
sb1=cumsum(b.*(-1).^(1:n)');
pt1=sb+sb1/3-b/3;
pt2=sb-sb1/3-b/3;
tau=zeros(n,DD(1)+1,DD(2)+1);
tau(1:2:n,1,1)=pt1(1:2:n);
tau(2:2:n,1,1)=pt2(2:2:n);
tau(n,1,1)=1;
tau(:,1,1)=min(tau(:,1,1),1);
% Find the EC densities:
u=df1*t;
kk=(max(DD)-1+min(DD))/2;
uu=conn_bsxfun(@power,u,0:.5:kk);
[ii,jj]=ndgrid(0:kk,0:kk);
gammalnii=gammalni(0:max(DD)+kk);
for d=1:max(DD)
for e=0:min(min(DD),d)
s1=0;
cons=-((d+e)/2+1)*log(pi)+gammaln(d)+gammaln(e+1);
for k=0:(d-1+e)/2
%[i,j]=ndgrid(0:k,0:k);
i=ii(1:k+1,1:k+1); j=jj(1:k+1,1:k+1);
if df2==Inf
q1=log(pi)/2-((d+e-1)/2+i+j)*log(2);
else
q1=(df0-1-d-e)*log(2)+gammaln((df0-d)/2+i)+gammaln((df0-e)/2+j) ...
-gammalni(df0-d-e+i+j+k)-((d+e-1)/2-k)*log(df2);
end
%q2=cons-gammalni(i+1)-gammalni(j+1)-gammalni(k-i-j+1) ...
% -gammalni(d-k-i+j)-gammalni(e-k-j+i+1);
q2=cons-gammalnii(1+max(0,i+1))-gammalnii(1+max(0,j+1))-gammalnii(1+max(0,k-i-j+1)) ...
-gammalnii(1+max(0,d-k-i+j))-gammalnii(1+max(0,e-k-j+i+1));
s2=sum(sum(exp(q1+q2)));
if s2>0
%s1=s1+(-1)^k*u.^((d+e-1)/2-k)*s2;
s1=s1+(-1)^k*uu(:,1+2*((d+e-1)/2-k))*s2;
end
end
if df2==Inf
s1=s1.*exp(-u/2);
else
s1=s1.*exp(-(df0-2)/2*log(1+u/df2));
end
if DD(1)>=DD(2)
tau(:,d+1,e+1)=s1;
if d<=min(DD)
tau(:,e+1,d+1)=s1;
end
else
tau(:,e+1,d+1)=s1;
if d<=min(DD)
tau(:,d+1,e+1)=s1;
end
end
end
end
% For multivariate statistics, add a sphere to the search region:
a=zeros(2,max(nvar));
for k=1:2
j=(nvar(k)-1):-2:0;
a(k,j+1)=exp(j*log(2)+j/2*log(pi) ...
+gammaln((nvar(k)+1)/2)-gammaln((nvar(k)+1-j)/2)-gammaln(j+1));
end
rho=zeros(n,Dlim(1)+1,Dlim(2)+1);
for k=1:nvar(1)
for l=1:nvar(2)
rho=rho+a(1,k)*a(2,l)*tau(:,(0:Dlim(1))+k,(0:Dlim(2))+l);
end
end
if is_tstat
t=[sqrt(t(1:(n-1))); -flipdim(sqrt(t),1)];
rho=[rho(1:(n-1),:,:); flipdim(rho,1)]/2;
for i=0:D(1)
for j=0:D(2)
rho(n-1+(1:n),i+1,j+1)=-(-1)^(i+j)*rho(n-1+(1:n),i+1,j+1);
end
end
rho(n-1+(1:n),1,1)=rho(n-1+(1:n),1,1)+1;
n=2*n-1;
end
% For scale space:
if scale>1
kappa=D(1)/2;
tau=zeros(n,D(1)+1);
for d=0:D(1)
s1=0;
for k=0:d/2
s1=s1+(-1)^k/(1-2*k)*exp(gammaln(d+1)-gammaln(k+1)-gammaln(d-2*k+1) ...
+(1/2-k)*log(kappa)-k*log(4*pi))*rho(:,d+2-2*k,1);
end
if d==0
cons=log(scale);
else
cons=(1-1/scale^d)/d;
end
tau(:,d+1)=rho(:,d+1,1)*(1+1/scale^d)/2+s1*cons;
end
rho(:,1:(D(1)+1),1)=tau;
end
if D(2)==0
d=D(1);
if nconj>1
% Conjunctions:
b=gamma(((0:d)+1)/2)/gamma(1/2);
for i=1:d+1
rho(:,i,1)=rho(:,i,1)/b(i);
end
m1=zeros(n,d+1,d+1);
for i=1:d+1
j=i:d+1;
m1(:,i,j)=rho(:,j-i+1,1);
end
for k=2:nconj
for i=1:d+1
for j=1:d+1
m2(:,i,j)=sum(rho(:,1:d+2-i,1).*m1(:,i:d+1,j),2);
end
end
m1=m2;
end
for i=1:d+1
rho(:,i,1)=m1(:,1,i)*b(i);
end
end
if ~isempty(EC_file)
if d<3
rho(:,(d+2):4,1)=zeros(n,4-d-2+1);
end
fid=fopen(EC_file,'w');
% first 3 are dimension sizes as 4-byte integers:
fwrite(fid,[n max(d+2,5) 1],'int');
% next 6 are bounding box as 4-byte floats:
fwrite(fid,[0 0 0; 1 1 1],'float');
% rest are the data as 4-byte floats:
if ~isempty(expr)
eval(expr);
end
fwrite(fid,t,'float');
fwrite(fid,rho,'float');
fclose(fid);
end
end
if all(fwhm>0)
pval_rf=zeros(n,1);
for i=1:D(1)+1
for j=1:D(2)+1
pval_rf=pval_rf+invol(1,i)*invol(2,j)*rho(:,i,j);
end
end
else
pval_rf=Inf;
end
% Bonferroni
pt=rho(:,1,1);
pval_bon=abs(prod(num_voxels))*pt;
% Minimum of the two:
pval=min(pval_rf,pval_bon);
tlim=1;
if p_val_peak(1) <= tlim
peak_threshold=minterp1(pval,t,p_val_peak);
if length(p_val_peak)<=nprint
peak_threshold;
end
else
% p_val_peak is treated as a peak value:
P_val_peak=interp1(t,pval,p_val_peak);
peak_threshold=P_val_peak;
if length(p_val_peak)<=nprint
P_val_peak;
end
end
if fwhm<=0 | any(num_voxels<0)
peak_threshold_1=p_val_peak+NaN;
extent_threshold=p_val_extent+NaN;
extent_threshold_1=extent_threshold;
return
end
% Cluster_threshold:
%###-Changed so that cluster_threshold is considered as cluster extent no matter what.
if cluster_threshold > eps
tt=cluster_threshold;
else
% cluster_threshold is treated as a probability:
tt=minterp1(pt,t,cluster_threshold);
Cluster_threshold=tt;
end
d=D(1);
rhoD=interp1(t,rho(:,d+1,1),tt);
p=interp1(t,pt,tt);
% Pre-selected peak:
pval=rho(:,d+1,1)./rhoD;
if p_val_peak(1) <= tlim
peak_threshold_1=minterp1(pval,t, p_val_peak);
if length(p_val_peak)<=nprint
peak_threshold_1;
end
else
% p_val_peak is treated as a peak value:
P_val_peak_1=interp1(t,pval,p_val_peak);
peak_threshold_1=P_val_peak_1;
if length(p_val_peak)<=nprint
P_val_peak_1;
end
end
if D(1)==0 | nconj>1 | nvar(1)>1 | D(2)>0 | scale>1
extent_threshold=p_val_extent+NaN;
extent_threshold_1=extent_threshold;
if length(p_val_extent)<=nprint
extent_threshold;
extent_threshold_1;
end
return
end
% Expected number of clusters:
%###-Change tlim to a small number so that p_val_extent is considered as cluster extent
tlim = eps;
EL=invol(1,d+1)*rhoD;
cons=gamma(d/2+1)*(4*log(2))^(d/2)/fwhm(1)^d*rhoD/p;
if df2==Inf & dfw1==Inf
if p_val_extent(1) <= tlim
pS=-log(1-p_val_extent)/EL;
extent_threshold=(-log(pS)).^(d/2)/cons;
pS=-log(1-p_val_extent);
extent_threshold_1=(-log(pS)).^(d/2)/cons;
if length(p_val_extent)<=nprint
extent_threshold;
extent_threshold_1;
end
else
% p_val_extent is now treated as a spatial extent:
pS=exp(-(p_val_extent*cons).^(2/d));
P_val_extent=1-exp(-pS*EL);
extent_threshold=P_val_extent;
P_val_extent_1=1-exp(-pS);
extent_threshold_1=P_val_extent_1;
if length(p_val_extent)<=nprint
P_val_extent;
P_val_extent_1;
end
end
else
% Find dbn of S by taking logs then using fft for convolution:
ny=2^12;
a=d/2;
b2=a*10*max(sqrt(2/(min(df1+df2,dfw1))),1);
if df2<Inf
b1=a*log((1-(1-0.000001)^(2/(df2-d)))*df2/2);
else
b1=a*log(-log(1-0.000001));
end
dy=(b2-b1)/ny;
b1=round(b1/dy)*dy;
y=((1:ny)'-1)*dy+b1;
numrv=1+(d+1)*(df2<Inf)+d*(dfw1<Inf)+(dfw2<Inf);
f=zeros(ny,numrv);
mu=zeros(1,numrv);
if df2<Inf
% Density of log(Beta(1,(df2-d)/2)^(d/2)):
yy=exp(y./a)/df2*2;
yy=yy.*(yy<1);
f(:,1)=(1-yy).^((df2-d)/2-1).*((df2-d)/2).*yy/a;
mu(1)=exp(gammaln(a+1)+gammaln((df2-d+2)/2)-gammaln((df2+2)/2)+a*log(df2/2));
else
% Density of log(exp(1)^(d/2)):
yy=exp(y./a);
f(:,1)=exp(-yy).*yy/a;
mu(1)=exp(gammaln(a+1));
end
nuv=[];
aav=[];
if df2<Inf
nuv=[df1+df2-d df2+2-(1:d)];
aav=[a repmat(-1/2,1,d)];
end
if dfw1<Inf
if dfw1>df_limit
nuv=[nuv dfw1-dfw1/dfw2-(0:(d-1))];
else
nuv=[nuv repmat(dfw1-dfw1/dfw2,1,d)];
end
aav=[aav repmat(1/2,1,d)];
end
if dfw2<Inf
nuv=[nuv dfw2];
aav=[aav -a];
end
for i=1:(numrv-1)
nu=nuv(i);
aa=aav(i);
yy=y/aa+log(nu);
% Density of log((chi^2_nu/nu)^aa):
f(:,i+1)=exp(nu/2*yy-exp(yy)/2-(nu/2)*log(2)-gammaln(nu/2))/abs(aa);
mu(i+1)=exp(gammaln(nu/2+aa)-gammaln(nu/2)-aa*log(nu/2));
end
% Check: plot(y,f); sum(f*dy,1) should be 1
omega=2*pi*((1:ny)'-1)/ny/dy;
shift=complex(cos(-b1*omega),sin(-b1*omega))*dy;
prodfft=prod(fft(f),2).*shift.^(numrv-1);
% Density of Y=log(B^(d/2)*U^(d/2)/sqrt(det(Q))):
ff=real(ifft(prodfft));
% Check: plot(y,ff); sum(ff*dy) should be 1
mu0=prod(mu);
% Check: plot(y,ff.*exp(y)); sum(ff.*exp(y)*dy.*(y<10)) should equal mu0
alpha=p/rhoD/mu0*fwhm(1)^d/(4*log(2))^(d/2);
% Integrate the density to get the p-value for one cluster:
pS=cumsum(ff(ny:-1:1))*dy;
pS=pS(ny:-1:1);
% The number of clusters is Poisson with mean EL:
pSmax=1-exp(-pS*EL);
if p_val_extent(1) <= tlim
yval=minterp1(-pSmax,y,-p_val_extent);
% Spatial extent is alpha*exp(Y) -dy/2 correction for mid-point rule:
extent_threshold=alpha*exp(yval-dy/2);
% For a single cluster:
yval=minterp1(-pS,y,-p_val_extent);
extent_threshold_1=alpha*exp(yval-dy/2);
if length(p_val_extent)<=nprint
extent_threshold;
extent_threshold_1;
end
else
% p_val_extent is now treated as a spatial extent:
P_val_extent=interp1(y,pSmax,log(p_val_extent/alpha)+dy/2);
extent_threshold=P_val_extent;
% For a single cluster:
P_val_extent_1=interp1(y,pS,log(p_val_extent/alpha)+dy/2);
extent_threshold_1=P_val_extent_1;
if length(p_val_extent)<=nprint
P_val_extent;
P_val_extent_1;
end
end
end
return
end
function x=gammalni(n);
i=find(n>=0);
x=Inf+n;
if ~isempty(i)
x(i)=gammaln(n(i));
end
return
end
function iy=minterp1(x,y,ix);
% interpolates only the monotonically increasing values of x at ix
n=length(x);
mx=x(1);
my=y(1);
xx=x(1);
for i=2:n
if x(i)>xx
xx=x(i);
mx=[mx xx];
my=[my y(i)];
end
end
iy=interp1(mx,my,ix);
return
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment