forked from Chaogan-Yan/REST
-
Notifications
You must be signed in to change notification settings - Fork 0
/
rest_AlphaSim.m
executable file
·150 lines (135 loc) · 5.85 KB
/
rest_AlphaSim.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
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="[email protected]">YAN Chao-Gan</a>; <a href="[email protected]">DONG Zhang-Ye</a> ; <a href="[email protected]">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 ([email protected]), DONG Zhang-Ye ([email protected]) and ZHU Wei-Xuan ([email protected]).\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]);