-
Notifications
You must be signed in to change notification settings - Fork 8
/
Copy pathrest_NormalityTest.m
executable file
·173 lines (148 loc) · 7.23 KB
/
rest_NormalityTest.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
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
function rest_NormalityTest(ADataDir,Method,AResultFilename,AMaskFilename)
% NormalityTest
% ADataDir: The set of data (.img, .hdr) that rest_NormalityTest import
%
% Method: We suggest users to view MATLAB "help" for Lilliefors and
% Jarque-Bera goodness-of-fit tests. In rest_NormalityTest, Method = 1, Lilliefors test
% function will be applied; else Jarque-Bera test will functio will be applied.
%
% AResultFilename: AResultFilename will be a combination of the raw data
% name and the Method name.
%
% AMaskFilename: The brain mask (e.g. 61*73*61 mask file) that rest_NormalityTest applied before computing.
% Output: p value map (Normally distributed if p is larger than a threshold)
%_________________________________________________________________
% Zang Zhen-Xiang ([email protected]); Dong Zhang-Ye ([email protected]);
% 2011-12-18
theElapsedTime =cputime;
fprintf('\nComputing with:\t"%s"', ADataDir);
[AllVolume,VoxelSize,theImgFileList, Header,nVolumn] =rest_to4d(ADataDir);
% examin the dimensions of the functional images and set mask
nDim1 = size(AllVolume,1); nDim2 = size(AllVolume,2); nDim3 = size(AllVolume,3);
BrainSize = [nDim1 nDim2 nDim3];
sampleLength = nVolumn;
if Method ==1,
AResultFilename=[AResultFilename,'_','Lilliefors'];
else
AResultFilename=[AResultFilename,'_','JB'];
end
%20070512 Saving a big 3D+time Dataset to small pieces by its first dimension to make this process run at least
% put pieces of 4D dataset to the temp dir determined by the current time
theTempDatasetDirName =sprintf('Norm_%d_%s', fix((1e4) *rem(now, 1) ),rest_misc('GetCurrentUser'));
theTempDatasetDir =[tempdir theTempDatasetDirName] ;
ans=rmdir(theTempDatasetDir, 's');%suppress the error msg
mkdir(tempdir, theTempDatasetDirName); %Matlab 6.5 compatible
%Save 3D+time Dataset's pieces to disk after ROI time course retrieved
Save1stDimPieces(theTempDatasetDir, AllVolume, 'dim1_');
clear AllVolume;%Free large memory
%mask selection, added by Xiaowei Song, 20070421
fprintf('\n\t Load mask "%s".', AMaskFilename);
mask=rest_loadmask(nDim1, nDim2, nDim3, AMaskFilename);
fprintf('\n\t Build mask.\tWait...');
mask =logical(mask);%Revise the mask to ensure that it contain only 0 and 1
mask = repmat(mask, [1, 1, 1, sampleLength]);
%Save mask pieces to disk to make this program at least run
Save1stDimPieces(theTempDatasetDir, mask, 'mask_');
fprintf('\n\t Normality is computing.\tWait...');
NumPieces_Dim1 =4; %Constant number to divide the first dimension to "NumPieces_Dim1" pieces
NumComputingCount =floor(nDim1/NumPieces_Dim1);
if NumComputingCount< (nDim1/NumPieces_Dim1),
NumComputingCount =NumComputingCount +1;
else
end
for x=1:(NumComputingCount), %20071129
%for x=1:(floor(nDim1/NumPieces_Dim1) +1)
rest_waitbar(x/(floor(nDim1/NumPieces_Dim1) +1), ...
'Performing Normality Tests. Please wait...', ...
'REST working','Child','NeedCancelBtn');
%Load cached pieces of Datasets
theFilename =fullfile(theTempDatasetDir, sprintf('dim1_%.8d', x));
theDim1Volume4D =Load1stDimVolume(theFilename);
theDim1Volume4D =double(theDim1Volume4D);
%Load and Apply the pieces' mask
theFilename =fullfile(theTempDatasetDir, sprintf('mask_%.8d', x));
theDim1Mask4D =Load1stDimVolume(theFilename);
theDim1Volume4D(~theDim1Mask4D)=0;
ResultNormalityBrain=rest_Normality(theDim1Volume4D,Method);
%Save to file
theFilename =fullfile(theTempDatasetDir, sprintf('result%.2d_%.8d', x));
save(theFilename, 'ResultNormalityBrain');
end
fprintf('.');
clear theDim1Volume4D theDim1Mask4D ResultNormalityBrain;
%Construct the 3D+time Dataset from files again
fprintf('\n\t ReConstructing 3D Dataset.\tWait...');
%Construct the Normality map's filenames, 20070905
% [pathstr, name, ext, versn] = fileparts(AResultFilename);
% ResultMaps =[];
% ResultMaps =[ResultMaps;{[pathstr, filesep ,name, ext]}];
%Reconstruct the Result correlation map from pieces
theDataset3D=zeros(nDim1, nDim2, nDim3);
for x=1:(NumComputingCount)
rest_waitbar(x/(floor(nDim1/NumPieces_Dim1)+1), ...
'3D Brain reconstructing. Please wait...', ...
'REST working','Child','NeedCancelBtn');
theFilename =fullfile(theTempDatasetDir,sprintf('result%.2d_%.8d', x));
%fprintf('\t%d',x);% Just for debugging
if x~=(floor(nDim1/NumPieces_Dim1)+1),
theDataset3D(((x-1)*NumPieces_Dim1+1):(x*NumPieces_Dim1),:,:)=Load1stDimVolume(theFilename);
else
theDataset3D(((x-1)*NumPieces_Dim1+1):end,:,:)=Load1stDimVolume(theFilename);
end
fprintf('.');
fprintf('\n\t Saving Normality map.\tWait...');
rest_writefile(single(theDataset3D), ...
AResultFilename, ...
BrainSize,VoxelSize,Header, 'single'); %Revised by YAN Chao-Gan, 090321. Result data will be stored in 'single' format. %'double');
end%end for
theElapsedTime =cputime - theElapsedTime;
fprintf('\n\t Normality compution over, elapsed time: %g seconds.\n', theElapsedTime);
%After Band pass filter, remove the temporary files
ans=rmdir(theTempDatasetDir, 's');%suppress the error msg
%end
%Save the 1st dimension of the 4D dataset to files
function Save1stDimPieces(ATempDir, A4DVolume, AFilenamePrefix)
NumPieces_Dim1 =4; %Constant number to divide the first dimension to "NumPieces_Dim1" pieces
NumComputingCount =floor(size(A4DVolume,1)/NumPieces_Dim1);
if NumComputingCount< (size(A4DVolume,1)/NumPieces_Dim1),
NumComputingCount =NumComputingCount +1;
else
end
for x = 1:(NumComputingCount),
%for x = 1:(floor(size(A4DVolume,1)/NumPieces_Dim1)+1)
rest_waitbar(x/(floor(size(A4DVolume,1)/NumPieces_Dim1)+1), ...
'Cut one Big 3D+time Dataset into pieces of 3D+time Dataset. Please wait...', ...
'REST working','Child','NeedCancelBtn');
theFilename =fullfile(ATempDir, sprintf('%s%.8d',AFilenamePrefix, x));
if x~=(floor(size(A4DVolume,1)/NumPieces_Dim1)+1)
the1stDim = A4DVolume(((x-1)*NumPieces_Dim1+1):(x*NumPieces_Dim1), :,:,:);
else
the1stDim = A4DVolume(((x-1)*NumPieces_Dim1+1):end, :,:,:);
end
save(theFilename, 'the1stDim');
end
%Load the 1st dimension of the 4D dataset from files, return a Matrix not a struct
function Result=Load1stDimVolume(AFilename)
Result =load(AFilename);
theFieldnames=fieldnames(Result);
% Result =eval(sprintf('Result.%s',the1stField));%remove the struct variable to any named variable with a matrix
Result = Result.(theFieldnames{1});
%Normality Tests
function ResultNormalityBrain=rest_Normality(ABrain4D,Method);
[nDim1, nDim2, nDim3, nDim4]=size(ABrain4D);
ABrain4D=reshape(ABrain4D,nDim1*nDim2*nDim3,nDim4);
A=zeros(nDim1*nDim2*nDim3,1);
%Remove the mean
for i=1:nDim1*nDim2*nDim3,
T=ABrain4D(i,:);
if var(T)~=0,
if Method==1,
[h,p]=lillietest(T);
else
[h,p]=jbtest(T);
end
A(i,:)=p;
else A(i,:)=0;
end
end
ResultNormalityBrain=reshape(A,nDim1,nDim2,nDim3);