%% Příprava workspace
clear all;
clc;
close all;

%% Načtení datasetů
% load('/home/karela/data/StrokeRatBrain/batches/20190610_092401_190610_STROKE_RAT_BRAIN_AMCAO9_1_1__01/experiment1/conv/data.mat');
% data_reg = load('/home/karela/data/StrokeRatBrain/batches/20190610_092401_190610_STROKE_RAT_BRAIN_AMCAO9_1_1__01/experiment1/reg/data.mat');

% load('/home/karela/data/StrokeRatBrain/batches/20190611_093556_190611_STROKE_RAT_BRAIN_AMCAO10_1_1_01/experiment1/conv/data.mat');
% data_reg = load('/home/karela/data/StrokeRatBrain/batches/20190611_093556_190611_STROKE_RAT_BRAIN_AMCAO10_1_1_01/experiment1/reg/data.mat')

% load('/home/karela/data/StrokeRatBrain/batches/20190614_094953_20190614_Stroke_RAT_BRAIN_AMCAO12_1_1_01/experiment1/conv/data.mat');
% data_reg = load('/home/karela/data/StrokeRatBrain/batches/20190614_094953_20190614_Stroke_RAT_BRAIN_AMCAO12_1_1_01/experiment1/reg/data.mat');

% load('/home/karela/data/StrokeRatBrain/batches/20190821_092549_190821_STROKE_RAT_BRAIN_AMCAO14_1_1_01/experiment1/conv/data.mat')
% data_reg = load('/home/karela/data/StrokeRatBrain/batches/20190821_092549_190821_STROKE_RAT_BRAIN_AMCAO14_1_1_01/experiment1/reg/data.mat')

% load('/home/karela/data/StrokeRatBrain/batches/20200610_101545_AMCAO_AC05_1_1_20200610_101545_AMCAO_AC05_1_1_01/experiment1/conv/data.mat');
% data_reg = load('/home/karela/data/StrokeRatBrain/batches/20200610_101545_AMCAO_AC05_1_1_20200610_101545_AMCAO_AC05_1_1_01/experiment1/reg/data.mat');

%load('/home/karela/data/StrokeRatBrain/batches/20200611_114511_AMCAO_AC08_1_1_20200611_114511_AMCAO_AC08_1_1_01/experiment1/conv/data.mat');
%data_reg = load('/home/karela/data/StrokeRatBrain/batches/20200611_114511_AMCAO_AC08_1_1_20200611_114511_AMCAO_AC08_1_1_01/experiment1/reg/data.mat');

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

%% Simulovaná data
% load('SyntheticEchoes_cartesian_reco');
% load('SyntheticEchoes_cartesian_ReferenceMaps');

%% Angiografický dataset
%load('/nfs/bubla/data/dce-mri/upt/StrokeRatBrain/ResliceTOFonCS/20200610_101545_AMCAO_AC05_1_1/TOF_resliced/TOF_resliced.mat');
%load('/nfs/bubla/data/dce-mri/upt/StrokeRatBrain/ResliceTOFonCS/20190610_092401_190610_STROKE_RAT_BRAIN_AMCAO9_1_1/TOF_resliced/TOF_resliced.mat');

data_reg = data_reg.data2;
nSlices=size(data_reg,1);
nSamples=size(data_reg,2);
[nRows,nCols]=size(data_reg{1,1});
data4D_reg=zeros(nRows,nCols,nSlices,nSamples);
for sl=1:nSlices
    for t=1:nSamples
        data4D_reg(:,:,sl,t)=imrotate(data_reg{sl,t},90);
    end
end

%% 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

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

% data4D=permute(data4D,[2 1 3 4]);
% 
data_reg = data4D_reg;
data = data4D;

% %% Simul
% data2D=DataCells2Mat(data2);
% data2D_1 = imresize(data2D, [1024 1024]);
% 
% % 
% mean_baseline = mean(data2D(:,:,:,1:15),4);
% %% Převod pro simul
% for i = 1:size(data2D,4)
% % data(:,:,:,i) = (data2D(:,:,:,i)-data2D(:,:,:,1))./data2D(:,:,:,1);
% data(:,:,:,i) = (data2D(:,:,:,i)-mean_baseline)./mean_baseline;
% end
% 
% data4D = data;
% 
% mean_baseline = mean(SIs(:,:,1:2000),3);
% %% Převod pro simul
% for i = 1:size(SIs,3)
% % data(:,:,:,i) = (data2D(:,:,:,i)-data2D(:,:,:,1))./data2D(:,:,:,1);
% SIs_1(:,:,i) = (SIs(:,:,i)-mean_baseline)./mean_baseline;
% end

%% Normalizace dat
%% Maximem
% max_data = max(data,[],'all');
% data = data./max_data;

% %AUC
% data = data./(sum(abs(data),4)/size(data,4));

% % Simul
% data = imresize(data, [1024 1024]);

% Medián
AUC_data = sum(abs(data),4)/size(data,4);
AUC_data = AUC_data(:);
AUC_data(AUC_data==0)=[];
med_data = median(AUC_data,'all');
data = data./med_data;

% data(isinf(data))=0;
% data(isnan(data))=0;

% data = normalize(data);
% data = data./(sum(abs(data),4)/size(data,4));
% data(isinf(data))=0;
% data(isnan(data))=0;

% load('aif_test')
% aif_test = aif_test';
% [output] = fit_biexp(aif_test);
% figure;plot(squeeze(aif_test));
% hold on
% plot(output(1,1,1,8):933,bi_exp(squeeze(output),output(1,1,1,8):933))
% xlabel('Časová osa [vzorky]')
% ylabel('AIF [a.u.]')

%% Automatic detection

[output_common] = automatic_detection(data,Ts);

% rgbImage2 = cat(3, data(:,:,7,30), output_common(:,:,7), TOF_resliced(:,:,7));
% imshow(rgbImage2);

%% Visual vs angio

I1 = data_reg(:,:,8,61);
I1_1 = I1;
I1 = mat2gray(I1);
I1 = I1.*255;
I1 = uint8(I1);

a = output_common(:,:,8);
a = uint8(~a);

I2 = I1.*a;
I3 = I1;
I3(a==0)=255;


I1(:,:,1) = I3;
I1(:,:,2) = I2;
I1(:,:,3) = I2;

% figure;
% subplot 121
% imshow(I1);
% subplot 122
% TOF_resliced = imrotate(TOF_resliced,180);
% TOF_resliced = fliplr(TOF_resliced);
% imshow(TOF_resliced(:,:,8),[])

figure;
subplot 121
imshow(I1,[0 25]);
subplot 122
imshow(I1_1,[0 25])

% %% Simul visual
% 
% % I1 = data2D(:,:,:,61);
% I1 = data2D_1(:,:,:,61);
% I1 = mat2gray(I1);
% I1 = I1.*255;
% I1 = uint8(I1);
% 
% a = output_common(:,:,:);
% im_ref = output_common(:,:,:);
% a = uint8(~a);
% 
% I2 = I1.*a;
% I3 = I1;
% I3(a==0)=255;
% 
% 
% I1(:,:,1) = I3;
% I1(:,:,2) = I2;
% I1(:,:,3) = I2;
% 
% figure;
% subplot 121
% imshow(I1);
% subplot 122
% imshow(RoisMap,[])

% %% Kvant.
% 
% ref_map = RoisMap;
% ref_map = ref_map > 37 & ref_map < 54;
% %ref_map = imresize(ref_map, [128 128]);
% %im_ref = imresize(im_ref, [1024 1024]);
% 
% figure;
% subplot 121
% imshow(im_ref);
% subplot 122
% imshow(ref_map);
% 
% 
% % Dice
% dice_coeff = dice(im_ref,ref_map);
% disp(dice_coeff);


% figure;
% imshowpair(output_common(:,:,6),TOF_resliced(:,:,6));

% median
% output_med5 = output_common(:,:,5);
% output_med5 = medfilt2(output_med5,[2 2]);
% 
% figure
% imshow(output_common(:,:,8))

aifs_mean = zeros(size(data,3),size(data,4));

for j = 1:11
[rows, columns]=find(output_common(:,:,j)==1);

%% Get AIF
aif = zeros(size(data,4)*size(data,3),length(rows));
aif = [];
aifs = [];
aif_s = [];

for i = 1:length(rows)
        %aifs = squeeze(data(rows(i,1),columns(i,1),:,:));  
        aifs= squeeze(data(rows(i,1),columns(i,1),j,:));
        aif_s(i,:) = aifs; 
        %aif_s(1,1:size(data,4))=aifs(1,:);
        %for k = 1:10
            %aif_s(1,k*size(data,4)+1:(k+1)*size(data,4))=aifs(k+1,:);
        %end
        %aif(:,i) = aif_s;
end

figure;
plot(aif_s');
xlabel('Časová osa [min]')
ylabel('AIF [a.u.]')

aif_mean = mean(aif_s,1);
figure;
plot(aif_mean);
xlabel('Časová osa [min]')
ylabel('AIF [a.u.]')

% hold on
% plot(aif_ref); 
% xlabel('Časová osa [vzorky]')
% ylabel('Koncentrace k.l. [a.u.]')
aif_e = isempty(aif_mean);

if aif_e==1;
aifs_mean(j,:) = zeros(1,size(data,4)); 
else
aifs_mean(j,:) = aif_mean; 
end

% area = trapz(aif);
% aif1 = aif/area;

% for z = 1:size(aif,2)
% aif1(:,z) = aif(:,z)/max(aif(:,z));
% end
% 
% aif1(isnan(aif1))=0;

% aif2 = zeros()
% for z = 1:size(aif,2)
% aif2(:,z) = aif1(~isnan(aif1(:,z)));
% end

% figure;
% %plot(aif1*Fs);
% plot(aifs(j,:));
% 
% 
% aif_mean = mean(aifs(j,:),1);
% figure;
% %plot(aif_mean*Fs);
% plot(aif_mean);

end

aifs_mean = mean(aifs_mean,1);
figure;
plot(timeAxis,aifs_mean);
xlabel('Časová osa [min]')
ylabel('AIF [a.u.]')

%aif_09 = aifs_mean;
%aif_10_1 = aifs_mean;
%aif_12 = aifs_mean;
%aif_14 = aifs_mean;
%aif_05 = aifs_mean;
%aif_08 = aifs_mean;
% aif_10 = aifs_mean;

%% Pseudoref AIF
% % AMCAO9
% figure;
% plot(timeAxis,squeeze(data4D(34,44,8,:)));
% hold on
% plot(timeAxis,aifs_mean);
% xlabel('Časová osa [s]')
% ylabel('AIF [a.u.]')

%output_over(:,:,:) = output1(output1(:,:,5)==output2(:,:,5) & output1(:,:,5)==output3(:,:,5))

% imshowpair(rgbImage1,data(:,:,5,end));

% full_image = rgbImage1 + data(:,:,5,end);
% full_image = imrotate(full_image,90);
% imshow(full_image);

%aif1 = squeeze(data(66,97,5,:));

% for i = 1:size(output,1)
%     for j = 1:size(output,2)
%         for k = 1:size(output,3)
%         if output(i,j,k,3)>1.5
%            output1(i,j,k) = output(i,j,k,3);
%         else
%            output1(i,j,k) = 0;
%         end
%         end
%     end
% end

% output1 = output(:,:,:,3);
% output2 = zeros(size(output1,1),size(output1,2));
% 
% for i = 1:size(output1,1)
%     for j = 1:size(output1,2)
%         if output1(i,j)< 1.5
%            output2(i,j) = 0;
%         else
%            output2(i,j) = output1(i,j);
%         end
%     end
% end
% for i=1:11
%         figure;
%         imagesc(output(:,:,i,3)-data(:,:,i,end),[0,300]);
% end