clear all;
close all;
clc;

fileIn = '/home/karela/data/StrokeRatBrain/batches/20200713_105447_AMCAO_AC10_1_1_20200713_105447_AMCAO_AC10_1_1_01/experiment1/conv/data.mat'
load(fileIn);

% convert cell data to 4D matrix (x,y,z,time)
nSlices=size(data2,1);
nSamples=size(data2,2);
[nRows,nCols]=size(data2{1,1});
data4D=zeros(nRows,nCols,nSlices,nSamples);
for sl=1:nSlices
    for t=1:nSamples
        data4D(:,:,sl,t)=imrotate(data2{sl,t},90);
    end
end

data = data4D;

outputFolder=[];
nSamples=666;
timeAxis=(0:1:nSamples-1)*info.acq.Ts/60;

%% Mean AIF model
aif = zeros(21,666);
aif05 = load('aif_05');
aif08 = load('aif_08');
aif09 = load('aif_09');
aif10_1 = load('aif_10_1');
aif10 = load('aif_10');
aif12 = load('aif_12');
aif14 = load('aif_14');

aif05 = aif05.aif_05;
aif08 = aif08.aif_08;
aif09 = aif09.aif_09;
aif10_1 = aif10_1.aif_10_1;
aif10 = aif10.aif_10;
aif12 = aif12.aif_12 ;
aif14 = aif14.aif_14;

aif(1,:) = aif05(1,1:666);
aif(2,:) = aif08(1,1:666);
aif(3,:) = aif09;
aif(4,:) = aif10(1,1:666);
aif(5,:) = aif10_1;
aif(6,:) = aif12;
aif(7,:) = aif14;

%% Mean AIF k-means
aif05 = load('aif05_k');
aif08 = load('aif08_k');
aif09 = load('aif09_k');
aif10_1 = load('aif10_1_k');
aif10 = load('aif10_k');
aif12 = load('aif12_k');
aif14 = load('aif14_k');

aif05 = aif05.aif05_k;
aif08 = aif08.aif08_k';
aif08 = aif08';
aif09 = aif09.aif09_k;
aif10_1 = aif10_1.aif10_1_k;
aif10 = aif10.aif10_k;
aif12 = aif12.aif12_k ;
aif14 = aif14.aif14_k;

aif(8,:) = aif05(1,1:666);
aif(9,:) = aif08(1,1:666);
aif(10,:) = aif09;
aif(11,:) = aif10(1,1:666);
aif(12,:) = aif10_1;
aif(13,:) = aif12;
aif(14,:) = aif14;

%% Mean AIF ref
aif05 = load('aif05_ref');
aif08 = load('aif08_ref');
aif09 = load('aif09_ref');
aif10_1 = load('aif10_1_ref');
aif10 = load('aif10_ref');
aif12 = load('aif12_ref');
aif14 = load('aif14_ref');

aif05 = aif05.aif05;
aif05 = aif05';
aif08 = aif08.aif08';
% aif08 = aif08';
aif09 = aif09.aif09;
aif09 = aif09';
aif10_1 = aif10_1.aif10_1;
aif10 = aif10.aif10;
aif10 = aif10';
aif12 = aif12.aif12 ;
aif12 = aif12';
aif14 = aif14.aif14;
aif14 = aif14';

aif(15,:) = aif05(1,1:666);
aif(16,:) = aif08(1,1:666);
aif(17,:) = aif09;
aif(18,:) = aif10(1,1:666);
aif(19,:) = aif10_1;
aif(20,:) = aif12;
aif(21,:) = aif14;

aif= aif./(sum(abs(aif),2)/size(aif,2));

% max_aif = max(aif(1:7,:),[],'all');
% aif(1:7,:) = aif(1:7,:)./max_aif;

% max_aif = max(aif(8:14,:),[],'all');
% aif(8:14,:) = aif(8:14,:)./max_aif;
% 
% max_aif = max(aif(15:21,:),[],'all');
% aif(15:21,:) = aif(15:21,:)./max_aif;

aif_mean_m = mean(aif(1:7,:),1);
aif_mean_k = mean(aif(8:14,:),1);

std_m = std(aif(1:7,:),1);
std_k = std(aif(8:14,:),1);

aif_mean_r = mean(aif(15:21,:),1);
std_r = std(aif(15:21,:),1);

%% Shade
% Model vs Ref
figure;
subplot 211
curve1 = aif_mean_m + std_m;
curve2 = aif_mean_m - std_m;
x2 = [timeAxis, fliplr(timeAxis)];
inBetween = [curve1, fliplr(curve2)];
fill(x2, inBetween, 'b','FaceAlpha',0.2, 'EdgeColor','none');
hold on;
plot(timeAxis, aif_mean_m, 'b', 'LineWidth', 1);

curve1 = aif_mean_r + std_r;
curve2 = aif_mean_r - std_r;
x2 = [timeAxis, fliplr(timeAxis)];
inBetween = [curve1, fliplr(curve2)];
fill(x2, inBetween, 'r','FaceAlpha',0.2, 'EdgeColor','none');
hold on;
plot(timeAxis, aif_mean_r, 'r', 'LineWidth', 1);
xlabel('Časová osa [min]')
ylabel('AIF [a.u.]')
legend('Bi-exponenciální model (std AIF)','Bi-exponenciální model (průměrná AIF)','Std MCA AIF','Průměrná MCA AIF')

% K-means vs Ref

subplot 212
curve1 = aif_mean_k + std_k;
curve2 = aif_mean_k - std_k;
x2 = [timeAxis, fliplr(timeAxis)];
inBetween = [curve1, fliplr(curve2)];
fill(x2, inBetween, 'b','FaceAlpha',0.2, 'EdgeColor','none');
hold on;
plot(timeAxis, aif_mean_k, 'b', 'LineWidth', 1);

curve1 = aif_mean_r + std_r;
curve2 = aif_mean_r - std_r;
x2 = [timeAxis, fliplr(timeAxis)];
inBetween = [curve1, fliplr(curve2)];
fill(x2, inBetween, 'r','FaceAlpha',0.2, 'EdgeColor','none');
hold on;
plot(timeAxis, aif_mean_r, 'r', 'LineWidth', 1);
xlabel('Časová osa [min]')
ylabel('AIF [a.u.]')
legend('Shluková analýza k-means (std AIF)','Shluková analýza k-means (průměrná AIF)','Std MCA AIF','Průměrná MCA AIF')


%% Vyhodnocení
% 09
figure;
subplot 221
plot(timeAxis,aif(3,:));
hold on
plot(timeAxis,aif(10,:));
hold on
plot(timeAxis,aif(17,:));
title('AMCAO9')
xlabel('Časová osa [min]')
ylabel('AIF [a.u.]')
legend('Bi-exponenciální model (AIF)','Shluková analýza k-means (AIF)','MCA AIF')

% 10_1
subplot 222
plot(timeAxis,aif(5,:));
hold on
plot(timeAxis,aif(12,:));
hold on
plot(timeAxis,aif(19,:));
title('AMCAO10')
xlabel('Časová osa [min]')
ylabel('AIF [a.u.]')
legend('Bi-exponenciální model (AIF)','Shluková analýza k-means (AIF)','MCA AIF')

% 12
subplot 223
plot(timeAxis,aif(6,:));
hold on
plot(timeAxis,aif(13,:));
hold on
plot(timeAxis,aif(20,:));
title('AMCAO12')
xlabel('Časová osa [min]')
ylabel('AIF [a.u.]')
legend('Bi-exponenciální model (AIF)','Shluková analýza k-means (AIF)','MCA AIF')

% 14
subplot 224
plot(timeAxis,aif(7,:));
hold on
plot(timeAxis,aif(14,:));
hold on
plot(timeAxis,aif(21,:));
title('AMCAO14')
xlabel('Časová osa [min]')
ylabel('AIF [a.u.]')
legend('Bi-exponenciální model (AIF)','Shluková analýza k-means (AIF)','MCA AIF')

% 05
figure;
subplot 221
plot(timeAxis,aif(1,:));
hold on
plot(timeAxis,aif(8,:));
hold on
plot(timeAxis,aif(15,:));
title('AMCAO-AC05')
xlabel('Časová osa [min]')
ylabel('AIF [a.u.]')
legend('Bi-exponenciální model (AIF)','Shluková analýza k-means (AIF)','MCA AIF')

% 08
subplot 222
plot(timeAxis,aif(2,:));
hold on
plot(timeAxis,aif(9,:));
hold on
plot(timeAxis,aif(16,:));
title('AMCAO-AC08')
xlabel('Časová osa [min]')
ylabel('AIF [a.u.]')
legend('Bi-exponenciální model (AIF)','Shluková analýza k-means (AIF)','MCA AIF')

% 10
subplot 223
plot(timeAxis,aif(4,:));
hold on
plot(timeAxis,aif(11,:));
hold on
plot(timeAxis,aif(18,:));
title('AMCAO-AC10')
xlabel('Časová osa [min]')
ylabel('AIF [a.u.]')
legend('Bi-exponenciální model (AIF)','Shluková analýza k-means (AIF)','MCA AIF')

%% Simul
load('aif_simul_biexp');
load('aif_simul_kmeans');
load('aif_simul_ref');
load('aif_simul_MCA')

nSamples=500;
timeAxis=(0:1:nSamples-1)*info.acq.Ts/60;

aif = zeros(3,500);
aif(1,:) = aif_mean;
aif(2,:) = Lowest1StMomentSig;
aif(3,:) = aif_simul_ref;
aif(4,:) = aif_simul_MCA;

aif= aif./(sum(abs(aif),2)/size(aif,2));

% MCA AIF vs Ref
figure;
plot(timeAxis,aif(4,:));
hold on
plot(timeAxis,aif(3,:));
xlabel('Časová osa [min]')
ylabel('AIF [a.u.]')
legend('MCA AIF','Referenční AIF')


% Model vs. Ref
figure;
subplot 211
plot(timeAxis,aif(1,:));
hold on
plot(timeAxis,aif(3,:));
xlabel('Časová osa [min]')
ylabel('AIF [a.u.]')
legend('Bi-exponenciální model','Referenční AIF')

% K-means vs Ref
subplot 212
plot(timeAxis,aif(2,:));
hold on
plot(timeAxis,aif(3,:));
xlabel('Časová osa [min]')
ylabel('AIF [a.u.]')
legend('Shluková analýza k-means','Referenční AIF')






