%% Trénování SVM klasifikátoru a predikce tříd - SNÍMKY Z VIDEO OFTALMOSKOPU
% Skript slouží pro natrénování klasifikátoru SVM a následnou predikci
% tříd cév.
% Skript je po natrénování klasifikátoru rozdělen do dvou částí, první část
% vyhodnocuje výsledek klasifikace na "raw" testovací množině. Druhá část
% skriptu slouží pro zpracování výsledku pro určitý testovací obraz, kde je
% zahrnuta i vizualizace výsledku.

% Tento skript slouží pro klasifikaci cév pouze okolo optického disku, je
% definována oblast zájmu klasifikace!

clc; clear;

%% Načtení trénovacích dat
load('vid_train_label_mat_vess.mat');
load('vid_train_priznak_mat_vess.mat');

%% Načtení již natrénovaného klasifikátoru SVM
load('SVMModel_rbf_vid.mat');
 
%% Trénování SVM
t = templateSVM('Standardize', 1, 'KernelFunction', 'rbf');
SVMModel_rbf = fitcecoc(vid_train_priznak_mat_vess, vid_train_label_mat_vess, 'Learners', t, 'FitPosterior', 1, 'ClassNames', [2, 3]);

%% Pro celou množinu testovacích dat
% Načtení testovacích dat
load('vid_test_label_mat_vess.mat');
load('vid_test_priznak_mat_vess.mat');

% Predikce tříd cév pomocí natrénovného SVM modelu
[label, ~] = predict(SVMModel_rbf_vid, vid_test_priznak_mat_vess);

% load('label_vid.mat');

vid_test_label_mat_vess = round(vid_test_label_mat_vess);
vid_test_label_mat_vess(vid_test_label_mat_vess == 1) = 2;

% Výpočet statistických metrik klasifikace
acc_all = double([]);
error_rate_all = double([]);
spravne = 0;
nespravne = 0;
for i = 1:length(label)
    if (label(i) == vid_test_label_mat_vess(i)) == 1
        spravne = spravne + 1;
    else
        nespravne = nespravne + 1;
    end
end

acc_pom_all = spravne / (length(label)) * 100;            % procento přesnosti klasifikace
acc_all = [acc_all; acc_pom_all];

error_rate_pom_all = nespravne / (length(label)) * 100;   % procento chybovosti klasifikace
error_rate_all = [error_rate_all; error_rate_pom_all];
   
% vytvoření matice záměn
C_all = confusionmat(vid_test_label_mat_vess, label, 'Order', [2, 3]);
CCH_all = confusionchart(C_all);

%% Pro jednotlivé testovací obrazy
% Načtení dat daného obrazu
load('vid_test_priznak_mat_vess1.mat');
load('vid_test_label_mat_vess1.mat');

% segmentace cév
seg_location = 'C:\Users\jana-\Documents\FEKT\DP\data\pulzace\klasifikace\test\vessels';   % cesta k maskám segmentovaných cév
addpath(seg_location);
seg_data = imageDatastore(seg_location);

% labely cév
class_location = 'C:\Users\jana-\Documents\FEKT\DP\data\pulzace\klasifikace\test\classes'; % cesta k labelům cév označneých experty
addpath(class_location);
class_data = imageDatastore(class_location);

% body OD
mask_location = 'C:\Users\jana-\Documents\FEKT\DP\data\pulzace\klasifikace\test\masks';    % cesta k ROI maskám
addpath(mask_location);
mask_data = fileDatastore(mask_location, 'ReadFcn', @load, 'FileExtensions', '.mat');

% načtení snímků a masek 
segmentace = im2double(imread('vid1_test_seg.tif'));
klas = imread('vid1_test_class.png');
OD = load('vid1_test_maskOD.mat');

%% Predikce tříd naučeného SVM modelu
[label_vess, ~] = predict(SVMModel_rbf_vid, vid_test_priznak_mat_vess1);

%% Vymezení oblasti zájmu jen na cévy
[~, g_seg, ~] = imsplit(segmentace);
g_seg(g_seg == 255) = 1; 
    
souradnice = OD.vyslednevrcholy2;
bw = poly2mask(souradnice(:, 1), souradnice(:, 2), size(g_seg, 1), size(g_seg, 2)); % definice polygonu procházejícího skrz body OD

measurements = regionprops(bw, 'Centroid');       % zisk středu masky
centroid = cat(1, measurements.Centroid);         % extrahování souřadnice středu OD
souradnice_x = floor(centroid(1));                % zaokrouhlení souřadnice x
souradnice_y = floor(centroid(2));                % zaokrouhlení souřadnice y

% vytvoření první masky
polomer_OD = souradnice_x - souradnice(1, 1);     % poloměr optického disku
prumer_OD = polomer_OD * 2;                       % průměr optického disku
[r, s] = size(g_seg);

% vytvoření kužnice 1.5 násobku OD
prumer_DD1 = prumer_OD * 1.5;
polomer_DD1 = floor(prumer_DD1 / 2); 

maska1 = zeros(r, s);                             % vytvoření masky1
[a, b] = meshgrid(1:s, 1:r); 
maska1((a - souradnice_x).^2 + (b - souradnice_y).^2 <= polomer_DD1.^2) = 1; % kružnice se středem v nalezených souřadnicích středu optického disku

% 
prumer_DD2 = prumer_OD * 1.8;
polomer_DD2 = floor(prumer_DD2 / 2); 

maska2 = zeros(r, s);                             % vytvoření masky2
[a, b] = meshgrid(1:s, 1:r); 
maska2((a - souradnice_x).^2 + (b - souradnice_y).^2 <= polomer_DD2.^2) = 1;

% rozdíl masek pro vymezení oblasti zájmu
maska_rozdil = maska2 - maska1;

ROI_cevy = immultiply(g_seg, maska_rozdil);  % násobení segmentační masky a ROI okolo OD

%% VESS predikce - převod vektoru do matice
ROI_cevy_pom = reshape(ROI_cevy, [], 1);
pom_lab_vess = zeros(size(ROI_cevy));
pom_lab_vess = reshape(pom_lab_vess, [], 1);
pom_ind1 = 1;
for m = 1:length(pom_lab_vess)
    if ROI_cevy_pom(m) == 1
        pom_lab_vess(m) = label_vess(pom_ind1);
        pom_ind1 = pom_ind1 + 1;
    end
end
pom_lab_vess = reshape(pom_lab_vess, size(ROI_cevy, 1), size(ROI_cevy, 2));
pom_final = pom_lab_vess;

%% Vizualizace výsledku klasifikace
vysledek_klas_R = size(pom_final);
vysledek_klas_G = size(pom_final);
vysledek_klas_B = size(pom_final);
for s = 1:size(pom_final, 1)
    for t = 1:size(pom_final, 2)
        if pom_final(s, t) == 0             % pixely pozadí -> černá
            vysledek_klas_R(s, t) = 0;
            vysledek_klas_G(s, t) = 0;
            vysledek_klas_B(s, t) = 0;
        elseif pom_final(s, t) == 2         % pixely tepen -> červená
            vysledek_klas_R(s, t) = 255;
            vysledek_klas_G(s, t) = 0;
            vysledek_klas_B(s, t) = 0;
        elseif pom_final(s, t) == 3         % pixely žil -> modrá
            vysledek_klas_R(s, t) = 0;
            vysledek_klas_G(s, t) = 0;
            vysledek_klas_B(s, t) = 255;
        end
    end
end      
   
vysledek_klas = cat(3, vysledek_klas_R, vysledek_klas_G, vysledek_klas_B); % vytvoření RGB obrazu  

klas_pom = im2double(klas);
[r_pom, g_pom, b_pom] = imsplit(klas_pom);
r_pom = immultiply(r_pom, maska_rozdil);
r_pom = immultiply(r_pom, g_seg);
g_pom = immultiply(g_pom, maska_rozdil);
g_pom = immultiply(g_pom, g_seg);
b_pom = immultiply(b_pom, maska_rozdil);
b_pom = immultiply(b_pom, g_seg);
klas_viz = cat(3, r_pom, g_pom, b_pom);

figure
subplot 121; imshow(klas_viz, []); title('Ground truth'); subplot 122; imshow(vysledek_klas, []); title('Klasifikace pomocí SVM');
