-
Notifications
You must be signed in to change notification settings - Fork 0
/
rp_Smoothest.m
139 lines (123 loc) · 6.22 KB
/
rp_Smoothest.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
function [dLh,resels,FWHM, nVoxels]=rp_Smoothest(StatsImgFile, MaskFile, DOF, VoxelSize)
% function [dLh,resels,FWHM, nVoxels]=y_Smoothest(StatsImgFile, MaskFile, DOF, VoxelSize)
% Function to estimate smoothness like smoothest in FSL.
% Reference: Flitney, D.E., & Jenkinson, M. 2000. Cluster Analysis Revisited. Tech. rept. Oxford Centre for Functional Magnetic Resonance Imaging of the Brain, Department of Clinical Neurology, Oxford University, Oxford, UK. TR00DF1.
% Input:
% StatsImgFile - The Z statistcal image file name.
% MaskFile - The mask file name. If empty (i.e., ''), then all voxels are included.
% DOF - Degree of freedom if residule images are used. Only effective when the time points in StatsImgFile is bigger than 1.
% VoxelSize - if StatsImgFile is not given as the file name but as the data matrix, then VoxelSize need to be specified also.
% Output:
% dLh - Smoothness estimated as sqrt(det(Lambda)), can be used in inference.
% resels - The size of one Resel: volume of space with dimensions FWHMx, FWHMy and FWHMz (in voxels).
% FWHM - FWHM in x, y, z. Note: in mm
% nVoxels - Number of voxels in the mask
%___________________________________________________________________________
% Written by YAN Chao-Gan 120120.
% The Nathan Kline Institute for Psychiatric Research, 140 Old Orangeburg Road, Orangeburg, NY 10962, USA
% Child Mind Institute, 445 Park Avenue, New York, NY 10022, USA
% The Phyllis Green and Randolph Cowen Institute for Pediatric Neuroscience, New York University Child Study Center, New York, NY 10016, USA
% ycg.yan@gmail.com
if nargin<3
DOF=100; %Used when nDimTimePoints > 1
end
if ischar(StatsImgFile)
[AllVolume, VoxelSize, Header] =rp_readfile(StatsImgFile);
else
AllVolume = StatsImgFile;
end
[nDim1 nDim2 nDim3 nDimTimePoints]=size(AllVolume);
if ~isempty(MaskFile)
[MaskData,MaskVox,MaskHead]=rp_readfile(MaskFile);
else
MaskData=ones(nDim1,nDim2,nDim3);
end
AllVolume=single(AllVolume);
if nDimTimePoints>2
AllVolume = (AllVolume-repmat(mean(AllVolume,4),[1,1,1, nDimTimePoints]))./repmat(std(AllVolume,0,4),[1,1,1, nDimTimePoints]); %Zero mean and one std
AllVolume(isnan(AllVolume))=0;
end
SSminus=[0 0 0];
S2=[0 0 0];
N=0;
for x=2:nDim1
for y=2:nDim2
for z=2:nDim3
if MaskData(x, y, z) && MaskData(x-1, y, z) && MaskData(x, y-1, z) && MaskData(x, y, z-1)
N=N+1;
for t=1:nDimTimePoints
SSminus(1) = SSminus(1) + AllVolume(x, y, z, t) * AllVolume(x-1, y, z, t);
SSminus(2) = SSminus(2) + AllVolume(x, y, z, t) * AllVolume(x, y-1, z, t);
SSminus(3) = SSminus(3) + AllVolume(x, y, z, t) * AllVolume(x, y, z-1, t);
S2(1) = S2(1) + 0.5 * ((AllVolume(x, y, z, t)^2) + (AllVolume(x-1, y, z, t)^2));
S2(2) = S2(2) + 0.5 * ((AllVolume(x, y, z, t)^2) + (AllVolume(x, y-1, z, t)^2));
S2(3) = S2(3) + 0.5 * ((AllVolume(x, y, z, t)^2) + (AllVolume(x, y, z-1, t)^2));
end
end
end
end
fprintf('.');
end
if SSminus(1)>0.99999999*S2(1)
SSminus(1)=0.99999999*S2(1);
warning('possibly biased smootheness in X');
end
if SSminus(2)>0.99999999*S2(2)
SSminus(2)=0.99999999*S2(2);
warning('possibly biased smootheness in Y');
end
if SSminus(3)>0.99999999*S2(3)
SSminus(3)=0.99999999*S2(3);
warning('possibly biased smootheness in Z');
end
sigmasq(1) = -1 / (4 * log(abs(SSminus(1)/S2(1))));
sigmasq(2) = -1 / (4 * log(abs(SSminus(2)/S2(2))));
sigmasq(3) = -1 / (4 * log(abs(SSminus(3)/S2(3))));
dLh=((sigmasq(1)*sigmasq(2)*sigmasq(3))^-0.5)*(8^-0.5);
if nDimTimePoints > 1
fprintf('DLH %f voxels^-3 before correcting for temporal DOF\n',dLh);
% lut(5) = 1.5423138; lut(6) = 1.3757105; lut(7) = 1.2842680;
% lut(8) = 1.2272151; lut(9) = 1.1885232; lut(10) = 1.1606988;
% lut(11) = 1.1398000; lut(12) = 1.1235677; lut(13) = 1.1106196;
% lut(14) = 1.1000651; lut(15) = 1.0913060; lut(16) = 1.0839261;
% lut(17) = 1.0776276; lut(18) = 1.0721920; lut(19) = 1.0674553;
% lut(20) = 1.0632924; lut(25) = 1.0483053; lut(30) = 1.0390117;
% lut(40) = 1.0281339; lut(50) = 1.0219834; lut(60) = 1.0180339;
% lut(70) = 1.0152850; lut(80) = 1.0132621; lut(90) = 1.0117115;
% lut(100) = 1.0104851; lut(150) = 1.0068808; lut(200) = 1.0051200;
% lut(300) = 1.0033865; lut(500) = 1.0020191;
lut(6) = 1.5423138; lut(7) = 1.3757105; lut(8) = 1.2842680;
lut(9) = 1.2272151; lut(10) = 1.1885232; lut(11) = 1.1606988;
lut(12) = 1.1398000; lut(13) = 1.1235677; lut(14) = 1.1106196;
lut(15) = 1.1000651; lut(16) = 1.0913060; lut(17) = 1.0839261;
lut(18) = 1.0776276; lut(19) = 1.0721920; lut(20) = 1.0674553;
lut(21) = 1.0632924; lut(26) = 1.0483053; lut(31) = 1.0390117;
lut(41) = 1.0281339; lut(51) = 1.0219834; lut(61) = 1.0180339;
lut(71) = 1.0152850; lut(81) = 1.0132621; lut(91) = 1.0117115;
lut(101) = 1.0104851; lut(151) = 1.0068808; lut(201) = 1.0051200;
lut(301) = 1.0033865; lut(501) = 1.0020191;
y = lut(lut~=0);
x = find(lut~=0);
% xi=[1:500];
xi=[1:501];
lut_interpolated=interp1(x,y,xi,'linear');
if (DOF < 6)
dLh=dLh * 1.1;
elseif (DOF>500)
dLh=dLh * (1.0321/DOF +1)^0.5;
else
% dLh=dLh * lut_interpolated(DOF)^0.5;
retval=(lut_interpolated(floor(DOF)+1)-lut_interpolated(floor(DOF)))*(DOF-floor(DOF)+1) + ...
lut_interpolated(floor(DOF)+1);
dLh=dLh * retval^0.5;
end
end
FWHM(1) = sqrt(8 * log(2) * sigmasq(1));
FWHM(2) = sqrt(8 * log(2) * sigmasq(2));
FWHM(3) = sqrt(8 * log(2) * sigmasq(3));
resels = FWHM(1)*FWHM(2)*FWHM(3);
fprintf('\nFWHMx = %f voxels\nFWHMy = %f voxels\nFWHMz = %f voxels\n',FWHM(1),FWHM(2),FWHM(3));
FWHM=FWHM.*VoxelSize;
fprintf('FWHMx = %f mm\nFWHMy = %f mm\nFWHMz = %f mm\n',FWHM(1),FWHM(2),FWHM(3));
nVoxels=length(find(MaskData));
fprintf('DLH = %f\nVOLUME = %d\nRESELS = %f\n',dLh,nVoxels,resels);