-
Notifications
You must be signed in to change notification settings - Fork 1
/
spmrt_corrsub.m
289 lines (243 loc) · 11 KB
/
spmrt_corrsub.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
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
function [rP,CIP,rC,CIC] = spmrt_corrsub(image1, image2, masks, metric)
% Computes the Pearson (vector-wise) and concordance correlations with CI
% between image1 and image2 for data included in the mask
%
% FORMAT [rP,CIP,rC,CIC]=spmrt_corrsub(image1,image1,masks)
%
% INPUT if no input the user is prompted
% image1 is the filename of an image (see spm_select)
% image2 is the filename of an image (see spm_select)
% masks is a cell array of filename of masks in same space as image1 and image2
% the assumption for masks that the 1st filename is the full
% brain mask, while the 2nd, 3rd and 4th are for gray matter,
% white matter, and CSF (input at least one filename)
% metric is 'Pearson', 'Concordance', or 'both'
%
% OUTPUT rP is the Pearson correlation coefficient
% CIP is the 95% Confidence interval of the Pearson correlation coefficient
% rC is the concordance correlation coefficient
% CIC is the 95% Confidence interval of the Concordance correlation coefficient
%
% Conditional input/output:
% - if masks is a single filename, it is assumed to be the full brain mask and
% rP, CIP, rC, CIC are returned as individual variables
% - if masks has multiple filenames, rP, CIP, rC, CIC are returned in a
% matrix format in which each row corresponds to the mask file
% - if masks 2,3,4 are not binary, it is assumed to be probability maps and
% correlation curves are computed
%
% Cyril Pernet
% --------------------------------------------------------------------------
% Copyright (C) spmrt
if nargin == 0
[image1,image2,masks]=spmrt_pickupfiles;
end
if ~exist('metric','var')
metric = 'both';
end
%% whole brain analysis
if iscell(masks)
M = masks{1};
N =size(masks,2);
else
M = masks(1,:);
N =size(masks,1);
end
disp('whole brain analysis')
[rP,CIP,rC,CIC]=spmrt_corr(image1,image2,M,metric); % call spmrt_corr
disp('--------------------')
%% analysis per tissue class
if N > 1 % size(masks,2)
% check all the masks are of the same data type
% ---------------------------------------------
binary_test = ones(1,N-1); % assume all is binary
for tc = 2:N
if iscell(masks)
M = spm_read_vols(spm_vol(masks{tc}));
else
M = spm_read_vols(spm_vol(masks(tc,:)));
end
if length(unique(M))>2
binary_test(tc-1) = 0;
end
end
if sum(binary_test)>0 && sum(binary_test)<N-1
error('tissue class masks are not of the same type')
end
% compute correlation per tissue class
% -------------------------------------
if sum(binary_test) == N-1 %% if all binary compute as for full brain
if strcmpi(metric,'Pearson') || strcmpi(metric,'Both')
tmp = NaN(N,1); tmp(1,1) = rP; rP = tmp; clear tmp
tmp = NaN(N,2); tmp(1,:) = CIP; CIP = tmp; clear tmp
end
if strcmpi(metric,'Concordance') || strcmpi(metric,'Both')
tmp = NaN(N,1); tmp(1,1) = rC; rC = tmp; clear tmp
tmp = NaN(N,2); tmp(1,:) = CIC; CIC = tmp; clear tmp
end
for tc = 2:N
if iscell(masks)
M = masks{tc};
else
M = masks(tc,:);
end
fprintf('computing correlations for tissue class %g\n',tc-1)
[tmpP,tmpcip,tmpC,tmpcic]=spmrt_corr(image1,image2,M,metric);
if strcmpi(metric,'Pearson') || strcmpi(metric,'Both')
rP(tc,1) = tmpP; CIP(tc,:) = tmpcip;
end
if strcmpi(metric,'Concordance') || strcmpi(metric,'Both')
rC(tc,1) = tmpC; CIC(tc,:) = tmpcic;
end
end
% if probability masks then iterate to create curves
% ---------------------------------------------------
elseif sum(binary_test) == 0
if strcmpi(metric,'Pearson') || strcmpi(metric,'Both')
tmp = NaN(N,10); tmp(1,1) = rP; rP = tmp; clear tmp
tmp = NaN(N,10,2); tmp(1,1,:) = CIP; CIP = tmp; clear tmp
end
if strcmpi(metric,'Concordance') || strcmpi(metric,'Both')
tmp = NaN(N,10); tmp(1,1) = rC; rC = tmp; clear tmp
tmp = NaN(N,10,2); tmp(1,1,:) = CIC; CIC = tmp; clear tmp
end
for tc = 2:N % for each tissue class
if iscell(masks)
M = masks{tc};
else
M = masks(tc,:);
end
for th=1:10 % for each threshold of the mask
fprintf('computing correlations for tissue class %g threshold %g \n',tc-1,th/10-0.1)
[tmpP,tmpcip,tmpC,tmpcic]=spmrt_corr(image1,image2,M,metric,0,th/10-0.1,0.005); % <-- note the adjustment for multiple comparisons 0.005
if strcmpi(metric,'Pearson') || strcmpi(metric,'Both')
rP(tc,th) = tmpP; CIP(tc,th,1) = tmpcip(1); CIP(tc,th,2) = tmpcip(2);
end
if strcmpi(metric,'Concordance') || strcmpi(metric,'Both')
rC(tc,th) = tmpC; CIC(tc,th,1) = tmpcic(1); CIC(tc,th,2) = tmpcic(2);
end
end % close the iterartion to crerate curves
end % close tissue type analysis
end
end
%% data viz
if N==1
figure('Name','Pearson Correlations per tissue class')
set(gcf,'Color','w','InvertHardCopy','off', 'units','normalized','outerposition',[0 0 1 1])
if iscell(masks);
M = masks;
else
M = masks;
end
X = spmrt_getdata(image1,image2,M);
scatter(X(:,1),X(:,2),50,[0 0 0]); hold on
h=lsline; set(h,'Color','k','LineWidth',4); % add the least square line
xlabel('img1','FontSize',14); ylabel('img2','FontSize',14); % label
box on;set(gca,'Fontsize',14); axis square; hold on
v = axis; plot([v(1):[(v(2)-v(1))/10]:v(2)],[v(3):[(v(4)-v(3))/10]:v(4)],'--k','LineWidth',2); % add diagonal
if strcmpi(metric,'Pearson')
title(['Whole image Pearsons'' corr =' num2str(rP)],'FontSize',12);
elseif strcmpi(metric,'Concordance')
title(['Whole image Concordance corr =' num2str(rC)],'FontSize',12);
elseif strcmpi(metric,'Both')
title(['Whole image Pearsons'' corr =' num2str(rP) ' Concordance corr =' num2str(rC)],'FontSize',12);
end
elseif sum(binary_test) == N-1
% for binary masks, simply plot data for all brain and tissue types
% ----------------------------------------------------------------
tricolors = [0 0 1; 1 0 0; 0 1 0];
if strcmpi(metric,'Pearson') || strcmpi(metric,'Both')
figure('Name','Pearson Correlations per tissue class')
set(gcf,'Color','w','InvertHardCopy','off', 'units','normalized','outerposition',[0 0 1 1])
for tc = N:-1:1
if iscell(masks);
M = masks{tc};
else
M = masks(tc,:);
end
subplot(ceil(N/2),2,tc);
if tc ~= 1
X{tc} = spmrt_getdata(image1,image2,M);
scatter(X{tc}(:,1),X{tc}(:,2),50,tricolors(tc-1,:)); grid on % plot pairs of observations
h=lsline; set(h,'Color','r','LineWidth',4); % add the least square line
else % whole brain mask
X{tc} = spmrt_getdata(image1,image2,M);
scatter(X{tc}(:,1),X{tc}(:,2),50,[1 1 1]); hold on
h=lsline; set(h,'Color','r','LineWidth',4); % add the least square line
for n=N:-1:2 % for each other tissue classes
scatter(X{n}(:,1),X{n}(:,2),50,tricolors(n-1,:)); grid on
end
end
xlabel('img1','FontSize',14); ylabel('img2','FontSize',14); % label
box on;set(gca,'Fontsize',14); axis square; hold on
v = axis; plot([v(1):[(v(2)-v(1))/10]:v(2)],[v(3):[(v(4)-v(3))/10]:v(4)],'--k','LineWidth',2); % add diagonal
title(['Mask ' num2str(tc) ': Pearsons'' corr =' num2str(rP(tc))],'FontSize',12);
end
end
if strcmpi(metric,'Concordance') || strcmpi(metric,'Both')
figure('Name','Concordance Correlations per tissue class')
set(gcf,'Color','w','InvertHardCopy','off', 'units','normalized','outerposition',[0 0 1 1])
for tc = N:-1:1
if iscell(masks);
M = masks{tc};
else
M = masks(tc,:);
end
subplot(ceil(N/2),2,tc);
if tc ~= 1
X{tc} = spmrt_getdata(image1,image2,M);
scatter(X{tc}(:,1),X{tc}(:,2),50,tricolors(tc-1,:)); grid on % plot pairs of observations
h=lsline; set(h,'Color','r','LineWidth',4); % add the least square line
else % whole brain mask
X{tc} = spmrt_getdata(image1,image2,M);
scatter(X{tc}(:,1),X{tc}(:,2),50,[1 1 1]); hold on
h=lsline; set(h,'Color','r','LineWidth',4); % add the least square line
for n=N:-1:2 % for each other tissue classes
scatter(X{n}(:,1),X{n}(:,2),50,tricolors(n-1,:)); grid on
end
end
xlabel('img1','FontSize',14); ylabel('img2','FontSize',14); % label
box on;set(gca,'Fontsize',14); axis square; hold on
v = axis; plot([v(1):[(v(2)-v(1))/10]:v(2)],[v(3):[(v(4)-v(3))/10]:v(4)],'--k','LineWidth',2); % add diagonal
title(['Mask ' num2str(tc) ': Concordance corr =' num2str(rC(tc))],'FontSize',12);
end
end
else
% for probability masks, plot data interactively depending on maks threshold
% --------------------------------------------------------------------------
if strcmpi(metric,'Both')
answer=questdlg('Plot concordance corr or Pearson''s corr?','Plotting option','Concordance','Pearson','Concordance');
else
answer = metric;
end
% create a large matrix with all data to plot
if strcmpi(answer,'concordance')
% Data dim
% N-1 tissue classes * 3 for CI low, r, CI high
% size(rC,2) for the number of thresholds
Data = NaN((N-1)*3,size(rC,2));
row_index = 0;
for n=1:N-1
Data(1+row_index,:) = CIC(n+1,:,1);
Data(2+row_index,:) = rC(n+1,:);
Data(3+row_index,:) = CIC(n+1,:,2);
row_index = row_index +3;
end
elseif strcmpi(answer,'pearson')
Data = NaN((N-1)*3,size(rP,2));
row_index = 0;
for n=1:N-1
Data(1+row_index,:) = CIP(n+1,:,1);
Data(2+row_index,:) = rP(n+1,:);
Data(3+row_index,:) = CIP(n+1,:,2);
row_index = row_index +3;
end
end
% get the data + tissue values and coordinates
for n=2:N
[tmp,x,y,z,M] = spmrt_getdata(image1,image2,masks(2,:));
tissue_values{n-1} = [tmp x y z M];
end
% now pass this info to the figure GUI
spmrt_tissue_reliability(Data,tissue_values)
end