Created
June 11, 2014 10:08
-
-
Save swederik/e23c3726ce311fab0bf7 to your computer and use it in GitHub Desktop.
CONN functional connectivity toolbox ROI-to-voxel stats exporting on Mac OSX (conn_vproject.m)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| 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. | |
| % keith.worsley@mcgill.ca , 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