-
Notifications
You must be signed in to change notification settings - Fork 1
/
spm_VOI.m
165 lines (149 loc) · 6.2 KB
/
spm_VOI.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
function TabDat = spm_VOI(SPM,VOL,Dis,Num,hReg)
% List of local maxima and adjusted p-values for a small Volume of Interest
% FORMAT TabDat = spm_VOI(SPM,VOL,Dis,Num,hReg)
%
% SPM - structure containing SPM, distribution & filtering details
% - required fields are:
% .swd - SPM working directory - directory containing current SPM.mat
% .Z - minimum of n Statistics {filtered on u and k}
% .n - number of conjoint tests
% .STAT - distribution {Z, T, X or F}
% .df - degrees of freedom [df{interest}, df{residual}]
% .u - height threshold
% .k - extent threshold {resels}
% .XYZ - location of voxels {voxel coords}
% .XYZmm - location of voxels {mm}
%
% VOL - structure containing details of volume analysed
% - required fields are:
% .S - search Volume {voxels}
% .R - search Volume {resels}
% .FWHM - smoothness {voxels}
% .M - voxels - > mm matrix
% .VOX - voxel dimensions {mm}
%
% Dis - Minimum distance between maxima
% (Passed to spm_list.m: Defaults on missing or empty)
% Num - Maxiumum number of local maxima tabulated per cluster
% (Passed to spm_list.m: Defaults on missing or empty)
% hReg - Handle of results section XYZ registry (see spm_results_ui.m)
%
% TabDat - Structure containing table data
% - see spm_list for definition
%
%_______________________________________________________________________
%
% spm_VOI is called by the SPM results section and takes variables in
% SPM and VOL to compute p-values corrected for a specified volume of
% interest.
%
% The volume of interest may be defined as a box or sphere centred on
% the current voxel, by the cluster in which it is embedded, or by an
% external mask image.
%
% If the VOI is defined by the cluster in which it is embedded, this
% cluster must have been defined independently of the SPM using a mask
% based on an orthogonal contrast and u = -Inf (i.e. p = 1)
%
% External mask images should be in the same orientation as the SPM
% (i.e. as the input used in stats estimation). The VOI is defined by
% voxels with values greater than 0.
%
% See also: spm_list
%_______________________________________________________________________
% @(#)spm_VOI.m 2.10 Karl Friston 99/11/17
%-Parse arguments
%-----------------------------------------------------------------------
if nargin<2, error('insufficient arguments'), end
if nargin<5, hReg = []; end
if nargin<4, Num = []; end
if isempty(Num), Num = 16; end
if nargin<3, Dis = []; end
if isempty(Dis), Dis = 04; end
%-Title
%-----------------------------------------------------------------------
spm('FigName',['SPM{',SPM.STAT,'}: Small Volume Correction']);
%-Get current location
%-----------------------------------------------------------------------
xyzmm = spm_results_ui('GetCoords');
%-Specify search volume
%-----------------------------------------------------------------------
str = sprintf(' at [%.0f,%.0f,%.0f]',xyzmm(1),xyzmm(2),xyzmm(3));
SPACE = spm_input('Search volume...',-1,'m',...
{['Sphere',str],['Box',str],'Nearest cluster',...
'Image'},['S','B','V','I']);
Q = ones(1,size(SPM.XYZmm,2));
vsc = [1 1 1]; %-Voxel size scaling for FWHM
switch SPACE, case 'S' % Sphere
%---------------------------------------------------------------
D = spm_input('radius of spherical VOI {mm}',-2);
str = sprintf('%0.1fmm sphere',D);
j = find(sum((SPM.XYZmm - xyzmm*Q).^2) <= D^2);
D = D./VOL.VOX;
S = (4/3)*pi*prod(D);
case 'B' % Box
%---------------------------------------------------------------
D = spm_input('box dimensions [k l m] {mm}',-2);
str = sprintf('%0.1f x %0.1f x %0.1f mm box',D(1),D(2),D(3));
j = find(all(abs(SPM.XYZmm - xyzmm*Q) <= D(:)*Q/2));
D = D(:)./VOL.VOX;
S = prod(D);
case 'V' %-Voxels
%---------------------------------------------------------------
if ~length(SPM.XYZ)
spm('alert!','No suprathreshold clusters!',mfilename,0);
spm('FigName',['SPM{',SPM.STAT,'}: Results']);
return
end
[xyzmm,i] = spm_XYZreg('NearestXYZ',xyzmm,SPM.XYZmm);
spm_results_ui('SetCoords',xyzmm);
A = spm_clusters(SPM.XYZ);
j = find(A == A(i));
str = sprintf('%0.0f voxel cluster',length(j));
D = SPM.XYZ(:,j);
S = length(j);
case 'I' % Image
%---------------------------------------------------------------
%-The VOI image must be in the same (real) space as the SPM:
%-Although the masking code below does take account of the
% .mat file for the VOI image, the orientation of the VOI image
% .mat file is *not* used in the resel calculation,
% (see spm_resels_vol.c) so any mat file should contain
% no rotations / shears relative to the SPM, for accuracy.
im = spm_get(1,'img','Image defining volume subset');
D = spm_vol(im);
tM = D.mat \ VOL.M;
if any(tM([2:5,7:10,12]))
spm('alert!',{ 'Mask image rotated/sheared!',...
'(relative to SPM image)',...
'Can''t use for SVC.'},mfilename,0);
spm('FigName',['SPM{',SPM.STAT,'}: Results']);
return
end
mXYZ = D.mat \ [SPM.XYZmm; ones(1, size(SPM.XYZmm, 2))];
j = spm_sample_vol(D, mXYZ(1,:),mXYZ(2,:),mXYZ(3,:),0) > 0;
str = sprintf('image mask: %s',spm_str_manip(im,'a30'));
%-Compute in-mask volume S:
% Correct for differences in mask and SPM voxel sizes
Y = spm_read_vols(D);
vsc = sqrt(sum(D.mat(1:3,1:3).^2)) ./ VOL.VOX';
S = sum(Y(:)>0) * prod(vsc);
end
spm('Pointer','Watch')
%-Select voxels within subspace
%-----------------------------------------------------------------------
SPM.Z = SPM.Z(j);
SPM.XYZ = SPM.XYZ(:,j);
SPM.XYZmm = SPM.XYZmm(:,j);
VOL.R = spm_resels(VOL.FWHM./vsc,D,SPACE);
VOL.S = S;
%-Tabulate p values
%-----------------------------------------------------------------------
str = sprintf('search volume: %s',str);
if any(strcmp(SPACE,{'S','B','V'}))
str = sprintf('%s at [%.0f,%.0f,%.0f]',str,xyzmm(1),xyzmm(2),xyzmm(3));
end
TabDat = spm_list('List',SPM,VOL,Dis,Num,str,hReg);
%-Reset title
%-----------------------------------------------------------------------
spm('FigName',['SPM{',SPM.STAT,'}: Results']);