https://github.com/Chaogan-Yan/REST
Tip revision: 5f9de23b90b10fd5eca9e9bc1016bbde75d43ab6 authored by Chaogan-Yan on 16 June 2013, 16:06:36 UTC
Fixed a bug in temporal correlation of two groups of images in Image Calculator;The midline of VMHC results were set to zero.
Fixed a bug in temporal correlation of two groups of images in Image Calculator;The midline of VMHC results were set to zero.
Tip revision: 5f9de23
rest_AlphaSim.m
function rest_AlphaSim(maskfile,outdir,outname,rmm,s,pthr,iter)
% rest_AlphaSim(maskfile,outdir,outname,rmm,s,pthr,iter)
% Monte Carlo simulation program similar to the AlphaSim in AFNI.
% The mechanism is based on AFNI's AlphaSim, please see more details from http://afni.nimh.nih.gov/pub/dist/doc/manual/AlphaSim.pdf
%------------------------------------------------------------------------------------------------------------------------------
% Input:
% maskfile - The image file indicate which voxels to analyze.
% outdir - The path to save the result file.
% outname - The filename of the result file.
% rmm - Cluster connection radius (mm).
% s - Gaussian filter width (FWHM, in mm).
% pthr - Individual voxel threshold probability
% iter - Number of Monte Carlo simulations.
% Output:
% outdir/outname.txt - You can find the cluster-size thresholds and voxel-wise p-value resulting a corrected P value in this file.
% By YAN Chao-Gan, Dong Zhang-Ye and ZHU Wei-Xuan 091108.
% State Key Laboratory of Cognitive Neuroscience and Learning, Beijing Normal University, China, 100875
% http://www.restfmri.net
% Mail to Authors: <a href="ycg.yan@gmail.com">YAN Chao-Gan</a>; <a href="dongzy08@gmail.com">DONG Zhang-Ye</a> ; <a href="zhuweixuan@gmail.com">ZHU Wei-Xuan</a>
% Version=1.0;
% Release=20091215;
% Modified by YAN Chao-Gan 0901215: Fixed the error when rmm=6.
%------------------------------------------------------------------------------------------------------------------------------
if ~(exist('spm_conv_vol.m'))
uiwait(msgbox('This function is based on SPM, please install SPM5 or later version at first.','REST AlphaSim'));
return
end
VariableLine=10000; %YAN Chao-Gan 091215. %VariableLine=3000;
[maskpath,maskname,masketc]=fileparts(maskfile);
[mask,voxdim,header]=rest_readfile(maskfile);
[nx,ny,nz] = size(mask);
mask = logical(mask);
nxyz = numel(find(mask));
dvoxel=max(voxdim);
outfilename=outname;
outfile=strcat(outdir,filesep,outfilename,'.txt');
if s == 4
s = 4.55; %afni 4 : spm 4.55
end
if rmm <= dvoxel*sqrt(2)
connect= 6;
else if rmm <= dvoxel*sqrt(3)
connect =18;
else
connect =26; %Revised by YAN Chao-Gan 091215. Fixed the error when rmm=6. %connect =27 ;
end
end
ft=zeros(1,VariableLine);
mt=zeros(1,VariableLine);
count=0;
suma=0;
sumsq=0;
for nt=1:iter
foneimt=zeros(1,VariableLine);
fim=normrnd(0,1,nx,ny,nz);
fim = fim.*mask;
if s ~= 0
fim = gauss_filter(s,fim,voxdim);
end
fimca=reshape(fim,1,[]);
count=count+nxyz;
suma=sum(fimca)+suma;
sumsq=sum(fimca.*fimca)+sumsq;
mean=suma/count;
sd = sqrt((sumsq - (suma * suma)/count) / (count-1));
zthr =norminv(1 - pthr);
xthr=sd*zthr+mean;
fim(fim<=xthr)=0;
fim(fim>xthr)=1;
a=numel(find(fim==1))/nxyz;
[theCluster, theCount] =bwlabeln(fim, connect);
for i=1:theCount
foneimt(numel(find(theCluster==i)))=foneimt(numel(find(theCluster==i)))+1;
end
for i=1:theCount
ft(numel(find(theCluster==i)))=ft(numel(find(theCluster==i)))+1;
end
mt(max(find(foneimt)))=mt(max(find(foneimt)))+1;
fprintf('iter=%d pvoxel=%f zthr=%f mc=%d mean=%f\n',nt,a,xthr,max(find(foneimt)),mean);
end
g_max_cluster_size = max(find(mt));
total_num_clusters = sum(ft);
divisor=iter*nxyz;
prob_table=zeros(1,g_max_cluster_size);
alpha_table=zeros(1,g_max_cluster_size);
cum_prop_table=zeros(1,g_max_cluster_size);
for i = 1:g_max_cluster_size
prob_table(i) = i * ft(i) / divisor;
alpha_table(i) = mt(i) / iter;
cum_prop_table(i) = ft(i) / total_num_clusters;
end
for i = 1:g_max_cluster_size-1
j = g_max_cluster_size - i +1;
prob_table(j-1) = prob_table(j)+prob_table(j-1);
alpha_table(j-1) = alpha_table(j)+alpha_table(j-1);
cum_prop_table(i+1) = cum_prop_table(i)+cum_prop_table(i+1);
end
fid=fopen(sprintf('%s',outfile),'w');
if(fid)
if s == 4.55
s=4;
end
fprintf(fid,'REST AlphaSim\nMonte Carlo simulation program similar to the AlphaSim in AFNI\nBy YAN Chao-Gan (ycg.yan@gmail.com), DONG Zhang-Ye (dongzy08@gmail.com) and ZHU Wei-Xuan (zhuweixuan@gmail.com).\nThe mechanism is based on AFNI''s AlphaSim, please see more details from http://afni.nimh.nih.gov/pub/dist/doc/manual/AlphaSim.pdf\nState Key Laboratory of Cognitive Neuroscience and Learning, Beijing Normal University, China, 100875\nhttp://www.restfmri.net\nVersion=1.0;\nRelease=20091105;\n\n\n');
fprintf(fid,'Mask filename = %s\n',maskname);
fprintf(fid,'Voxels in mask = %d\n',nxyz);
fprintf(fid,'Gaussian filter width (FWHM, in mm) = %.3f\n',s);
fprintf(fid,'Cluster connection radius: rmm = %.2f\n',rmm);
fprintf(fid,'Individual voxel threshold probability = %.3f\n',pthr);
fprintf(fid,'Number of Monte Carlo simulations = %d\n',iter);
fprintf(fid,'Output filename = %s\n\n\n',outfilename);
fprintf(fid,'Cl Size\tFrequency\tCum Prop\tp/Voxel\tMax Freq\tAlpha\n');
for i=1:g_max_cluster_size
fprintf(fid,'%d\t\t%d\t\t%f\t%f\t%d\t\t%f\n',i,ft(i),cum_prop_table(i),prob_table(i),mt(i),alpha_table(i));
end
fclose(fid);
end
function Q=gauss_filter(s,P,VOX)
if length(s) == 1; s = [s s s];end
s = s./VOX; % voxel anisotropy
s = max(s,ones(size(s))); % lower bound on FWHM
s = s/sqrt(8*log(2)); % FWHM -> Gaussian parameter
x = round(6*s(1)); x = [-x:x];
y = round(6*s(2)); y = [-y:y];
z = round(6*s(3)); z = [-z:z];
x = exp(-(x).^2/(2*(s(1)).^2));
y = exp(-(y).^2/(2*(s(2)).^2));
z = exp(-(z).^2/(2*(s(3)).^2));
x = x/sum(x);
y = y/sum(y);
z = z/sum(z);
i = (length(x) - 1)/2;
j = (length(y) - 1)/2;
k = (length(z) - 1)/2;
Q=P;
spm_conv_vol(P,Q,x,y,z,-[i,j,k]);