-
Notifications
You must be signed in to change notification settings - Fork 2
/
GraphCutOnStarplusData.m
67 lines (64 loc) · 2.96 KB
/
GraphCutOnStarplusData.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
% ÷èòàåì äàííûå
data = Read4DArrayFromStarplus(4847, 5);
mdata = mean(data,4);
imagesc(mdata(:,:,1));
roi_names = ["CALC", "LIPL_LT_LDLPFC", "LTRIA", "LOPER", "LIPS" ];
%roi_names = ["CALC_LFEF_RFEF_LSGA_RSGA_LSPL_RSPL_LIT_RIT", "LDLPFC_RDLPFC_LOPER_ROPER", "LIFG_LPPREC_RPPREC_LTRIA_RTRIA_LIPL_RIPL", "LT_RT", "SMA_LIPS_RIPS"];
%roi_names = ["CALC_LFEF_RFEF_LSGA_RSGA_LSPL_RSPL_LIT_RIT", "LDLPFC_RDLPFC_LOPER_ROPER"];
k = 6;
% ðàñêëàäûâàåì äàííûå ïî ðÿäàì è íîðìàëèçóåì
timeseries_by_rows = reshape(data, [size(data, 1) * size(data, 2) * size(data, 3), size(data, 4)]);
timeseries_by_rows = double(timeseries_by_rows);
% íàõîäèì íåïóñòûå âîêñåëè, íóæíî, ÷òîáû óáèðàòü èç ðàññìîòðåíèÿ ïóñòûå
non_empty_rows = sqrt(sum(timeseries_by_rows .^ 2, 2)) ~= 0;
%timeseries_by_rows = timeseries_by_rows - mean(timeseries_by_rows(:));
timeseries_by_rows = NormalizeToUnitLength(timeseries_by_rows);
% èíèöèàëèçèðóåì ñòðóêòóðû äëÿ ïàðàìåòðîâ ðàñïðåäåëåíèÿ
vmf = vmffactory(size(data, 4));
mus = zeros(k, size(data, 4));
kappas = zeros(1, k);
% êîìáèíèðîâàííàÿ ìàñêà âñåõ ROI, íóæíà, ÷òîáû íàéòè íå ROI îáëàñòü
sum_roi_mask = zeros([size(timeseries_by_rows,1), 1]);
% íàñòîÿùåå ðàçáèåíèå íà ROI
groundTruth = zeros([size(data, 1), size(data, 2), size(data, 3)]);
for i=1:(k-1)
roi_name = roi_names(i);
% äîñòàåì äàííûå ñîîòâåòñòâóþùèå ðàçìå÷åííîé îáëàñòè
roi_mask = ReadROIFromStarplus(4847, roi_name);
groundTruth = groundTruth + roi_mask * i;
roi_mask = roi_mask(:);
sum_roi_mask = sum_roi_mask | roi_mask;
roi = timeseries_by_rows(non_empty_rows & roi_mask==1, :);
roi = roi(sqrt(sum(roi .^ 2, 2)) ~= 0, :);
[theta] = vmf.estimatedefault(roi');
mus(i, :) = squeeze(theta.mu);
kappas(i) = theta.kappa;
end
%
% ðàñïðåäåëåíèå íå ROI îáëàñòè
non_roi_mask = ~sum_roi_mask;
groundTruth = groundTruth + reshape(k .* non_roi_mask .* (sqrt(sum(timeseries_by_rows .^ 2, 2)) ~= 0), size(groundTruth));
roi = timeseries_by_rows(non_roi_mask==1, :);
roi = roi(sqrt(sum(roi .^ 2, 2)) ~= 0, :);
[theta] = vmf.estimatedefault(roi');
mus(k, :) = squeeze(theta.mu);
kappas(k) = theta.kappa;
% kappas = ones(1,3);
% ñðàâíåíèå âåêòîðîâ mu
% dist = squareform(pdist(mus, 'cosine'));
% ñ÷èòàåì âåðîÿòíîñòè P è -logP
[probs, logprobs] = CalculateLikelihoodProbabilities(timeseries_by_rows, k, kappas, mus);
% íàõîäèì îöåíêó ìàêñèìàëüíîãî ïðàâäîïîäîáèÿ
segment_init = MLE(timeseries_by_rows, probs);
% ïóñòûì âîêñåëÿì ïðèñâàèâàåì ìåòêó 0
segment_init(~non_empty_rows) = 0;
% òðàíñôîðìèðóåì ñåãìåíòàöèþ â òðåõìåðíóþ ìàòðèöó
segment_init = reshape(segment_init, [size(data,1), size(data,2), size(data,3)]);
% èùåì MAP-îöåíêó
[map, energy] = MRF_MAP_GraphCutABSwap(segment_init, logprobs, 1, k, 10);
[beta, mu, kappa] = EstimateParametersOneSample(reshape(timeseries_by_rows, size(data)), map, energy, k, 13, 1, mus, kappas);
mdata = mean(data, 4);
z = 3;
imagesc([segment_init(:,:,z),map(:,:,z),groundTruth(:,:,z)]);
figure();
imagesc(mdata(:,:,z));