%% Skript pro extrahování příznaků z retinálních cév za účelem jejich klasifikace
% Načtené snímky a k nim příslušná data jsou zpracovávany do takové podoby,
% aby bylo možno extrahovat příznaky cév. Stěžejním výstupem skriptu jsou
% proměnné: priznak_mat_vess a label_mat_vess. Tyto proměnné jsou vstupy 
% do SVM klasifikátoru.

% Skript slouží pro klasifikaci cév sítnice z dat z video oftalmoskopu.
% Jelikož klasifikace předchází analýze pulzací řečiště, klasikace cév je
% realizována jen na úsecích cév, ve kterých jsou také analyzovány pulzace.

clc; clear;

%% Načtení cest k datům
% snímky sítnice
img_location = 'C:\Users\jana-\Documents\FEKT\DP\data\pulzace\klasifikace\train\images';    % cesta ke snímkům
addpath(img_location);
img_data = imageDatastore(img_location);
pocet = length(img_data.Files);                                                             % počet dat

% segmentace cév
seg_location = 'C:\Users\jana-\Documents\FEKT\DP\data\pulzace\klasifikace\train\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\train\classes'; % cesta k labelům cév
addpath(class_location);
class_data = imageDatastore(class_location);

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

%% For cycklus pro procházení jednotlivých snímků a jim příslušných dat (segmentace, labely, masky)
priznak_mat_vess = double([]);
label_mat_vess = double([]);

for i = 1:pocet
    % načtení snímků a masek 
    img_RGB = im2double(imread(strcat('vid', num2str(i), '_train_img.png')));
    segmentace = im2double(imread(strcat('vid', num2str(i), '_train_seg.tif')));
    klas = imread(strcat('vid', num2str(i), '_train_class.png'));
    OD = load(strcat('vid', num2str(i), '_train_maskOD.mat'));

    %% RGB, HSV prostor
    [redChannel_nativ, greenChannel_nativ, blueChannel_nativ] = imsplit(img_RGB);
    img_HSV = rgb2hsv(img_RGB);
    [hueChannel_nativ, saturationChannel_nativ, valueChannel_nativ] = imsplit(img_HSV);
    
    % adaptivní ekvalizace histogramu
    redChannel = adapthisteq(redChannel_nativ);
    greenChannel = adapthisteq(greenChannel_nativ);
    blueChannel = adapthisteq(blueChannel_nativ);
    hueChannel = adapthisteq(hueChannel_nativ);
    saturationChannel = adapthisteq(saturationChannel_nativ);
    valueChannel = adapthisteq(valueChannel_nativ);
    
    %% Zisk středu OD
    souradnice = OD.vyslednevrcholy2;
    bw = poly2mask(souradnice(:, 1), souradnice(:, 2), size(img_RGB, 1), size(img_RGB, 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

    %% Vymezení ROI
    % 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(redChannel);

    % 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;
    
    %% Násobení masek
    [r, g_seg, b] = imsplit(segmentace);
    g_seg(g_seg == 255) = 1;                          % změna hodnot v segmentační masce na 1
    
    % vytvoření skeletonu
    BW = imbinarize(g_seg);                           % binarizace segmentovaného obrazu
    skeleton = bwskel(BW);                            % vytvoření skeletonu cévního stromu
    branch = bwmorph(skeleton, 'branchpoints');
    skeleton(branch) = 0;                             % odstranění bifurkací a bodů křížení
    
    ROI_cevy = immultiply(g_seg, maska_rozdil);
    ROI_skeleton = immultiply(ROI_cevy, skeleton);    % násobení ROI a skeletonu
    
    % vynásobení RGB obrazu a segmentační masky
    img_nasobeny_R_vess = immultiply(redChannel, ROI_cevy);
    img_nasobeny_G_vess = immultiply(greenChannel, ROI_cevy);
    img_nasobeny_B_vess = immultiply(blueChannel, ROI_cevy);
    
    % vynásobení RGB obrazu a skeletonu
    img_nasobeny_R_skeleton = immultiply(redChannel, ROI_skeleton);
    img_nasobeny_G_skeleton = immultiply(greenChannel, ROI_skeleton);
    img_nasobeny_B_skeleton = immultiply(blueChannel, ROI_skeleton);

    % vynásobení HSV obrazu a segmentační masky
    img_nasobeny_H_vess = immultiply(hueChannel, ROI_cevy);
    img_nasobeny_S_vess = immultiply(saturationChannel, ROI_cevy);
    img_nasobeny_V_vess = immultiply(valueChannel, ROI_cevy);
    
    % vynásobení HSV obrazu a skeletonu
    img_nasobeny_H_skeleton = immultiply(hueChannel, ROI_skeleton);
    img_nasobeny_S_skeleton = immultiply(saturationChannel, ROI_skeleton);
    img_nasobeny_V_skeleton = immultiply(valueChannel, ROI_skeleton);

    %% Intenzita centrálních pixelů v gaussian blurred RGB
    % sigma = 1
    gauss1_R_skeleton = imgaussfilt(img_nasobeny_R_skeleton, 1);
    gauss1_G_skeleton = imgaussfilt(img_nasobeny_G_skeleton, 1);
    gauss1_B_skeleton = imgaussfilt(img_nasobeny_B_skeleton, 1);
    gauss1_R_skeleton = immultiply(gauss1_R_skeleton, g_seg);
    gauss1_G_skeleton = immultiply(gauss1_G_skeleton, g_seg);
    gauss1_B_skeleton = immultiply(gauss1_B_skeleton, g_seg);
    norm_priznak_gauss1_R_skeleton = normalizace_priznaku(gauss1_R_skeleton, ROI_cevy);
    norm_priznak_gauss1_G_skeleton = normalizace_priznaku(gauss1_G_skeleton, ROI_cevy);
    norm_priznak_gauss1_B_skeleton = normalizace_priznaku(gauss1_B_skeleton, ROI_cevy);

    % sigma = 2
    gauss2_R_skeleton = imgaussfilt(img_nasobeny_R_skeleton, 2);
    gauss2_G_skeleton = imgaussfilt(img_nasobeny_G_skeleton, 2);
    gauss2_B_skeleton = imgaussfilt(img_nasobeny_B_skeleton, 2);
    gauss2_R_skeleton = immultiply(gauss2_R_skeleton, g_seg);
    gauss2_G_skeleton = immultiply(gauss2_G_skeleton, g_seg);
    gauss2_B_skeleton = immultiply(gauss2_B_skeleton, g_seg);
    norm_priznak_gauss2_R_skeleton = normalizace_priznaku(gauss2_R_skeleton, ROI_cevy);
    norm_priznak_gauss2_G_skeleton = normalizace_priznaku(gauss2_G_skeleton, ROI_cevy);
    norm_priznak_gauss2_B_skeleton = normalizace_priznaku(gauss2_B_skeleton, ROI_cevy);

    % sigma = 4
    gauss4_R_skeleton = imgaussfilt(img_nasobeny_R_skeleton, 4);
    gauss4_G_skeleton = imgaussfilt(img_nasobeny_G_skeleton, 4);
    gauss4_B_skeleton = imgaussfilt(img_nasobeny_B_skeleton, 4);
    gauss4_R_skeleton = immultiply(gauss4_R_skeleton, g_seg);
    gauss4_G_skeleton = immultiply(gauss4_G_skeleton, g_seg);
    gauss4_B_skeleton = immultiply(gauss4_B_skeleton, g_seg);
    norm_priznak_gauss4_R_skeleton = normalizace_priznaku(gauss4_R_skeleton, ROI_cevy);
    norm_priznak_gauss4_G_skeleton = normalizace_priznaku(gauss4_G_skeleton, ROI_cevy);
    norm_priznak_gauss4_B_skeleton = normalizace_priznaku(gauss4_B_skeleton, ROI_cevy);

    %% Průměr + směrodatná odchylka + rozptyl intenzity v červeném, zeleném a modrém kanálu cév
    okno = 3;

    mean_R_vess = conv2(img_nasobeny_R_vess, ones(okno)/okno^2, 'same');
    mean_G_vess = conv2(img_nasobeny_G_vess, ones(okno)/okno^2, 'same');
    mean_B_vess = conv2(img_nasobeny_B_vess, ones(okno)/okno^2, 'same');
    mean_R_vess = immultiply(mean_R_vess, g_seg);
    mean_G_vess = immultiply(mean_G_vess, g_seg);
    mean_B_vess = immultiply(mean_B_vess, g_seg);
    norm_priznak_mean_R_vess = normalizace_priznaku(mean_R_vess, ROI_cevy);
    norm_priznak_mean_G_vess = normalizace_priznaku(mean_G_vess, ROI_cevy);
    norm_priznak_mean_B_vess = normalizace_priznaku(mean_B_vess, ROI_cevy);

    std_R_vess = stdfilt(img_nasobeny_R_vess, ones(okno));
    std_G_vess = stdfilt(img_nasobeny_G_vess, ones(okno));
    std_B_vess = stdfilt(img_nasobeny_B_vess, ones(okno));
    std_R_vess = immultiply(std_R_vess, g_seg);
    std_G_vess = immultiply(std_G_vess, g_seg);
    std_B_vess = immultiply(std_B_vess, g_seg);
    norm_priznak_std_R_vess = normalizace_priznaku(std_R_vess, ROI_cevy);
    norm_priznak_std_G_vess = normalizace_priznaku(std_G_vess, ROI_cevy);
    norm_priznak_std_B_vess = normalizace_priznaku(std_B_vess, ROI_cevy);

    %% Průměr + směrodatná odchylka + rozptyl odstínu, saturace a intenzity cév
    mean_H_vess = conv2(img_nasobeny_H_vess, ones(okno)/okno^2, 'same');
    mean_S_vess = conv2(img_nasobeny_S_vess, ones(okno)/okno^2, 'same');
    mean_V_vess = conv2(img_nasobeny_V_vess, ones(okno)/okno^2, 'same');
    mean_H_vess = immultiply(mean_H_vess, g_seg);
    mean_S_vess = immultiply(mean_S_vess, g_seg);
    mean_V_vess = immultiply(mean_V_vess, g_seg);
    norm_priznak_mean_H_vess = normalizace_priznaku(mean_H_vess, ROI_cevy);
    norm_priznak_mean_S_vess = normalizace_priznaku(mean_S_vess, ROI_cevy);
    norm_priznak_mean_V_vess = normalizace_priznaku(mean_V_vess, ROI_cevy);

    std_H_vess = stdfilt(img_nasobeny_H_vess, ones(okno));
    std_S_vess = stdfilt(img_nasobeny_S_vess, ones(okno));
    std_V_vess = stdfilt(img_nasobeny_V_vess, ones(okno));
    std_H_vess = immultiply(std_H_vess, g_seg);
    std_S_vess = immultiply(std_S_vess, g_seg);
    std_V_vess = immultiply(std_V_vess, g_seg);
    norm_priznak_std_H_vess = normalizace_priznaku(std_H_vess, ROI_cevy);
    norm_priznak_std_S_vess = normalizace_priznaku(std_S_vess, ROI_cevy);
    norm_priznak_std_V_vess = normalizace_priznaku(std_V_vess, ROI_cevy);

    %% Vytvoření matice extrahovaných příznaků pro vstup do SVM
    priznak_pom_vess = double([norm_priznak_gauss1_R_skeleton, norm_priznak_gauss1_G_skeleton, norm_priznak_gauss1_B_skeleton, ...
        norm_priznak_gauss2_R_skeleton, norm_priznak_gauss2_G_skeleton, norm_priznak_gauss2_B_skeleton, ...
        norm_priznak_gauss4_R_skeleton, norm_priznak_gauss4_G_skeleton, norm_priznak_gauss4_B_skeleton, ...
        norm_priznak_mean_R_vess, norm_priznak_mean_G_vess, norm_priznak_mean_B_vess, ...
        norm_priznak_std_R_vess, norm_priznak_std_G_vess, norm_priznak_std_B_vess, ...
        norm_priznak_mean_H_vess, norm_priznak_mean_S_vess, norm_priznak_mean_V_vess, ...
        norm_priznak_std_H_vess, norm_priznak_std_V_vess]);
    priznak_mat_vess = [priznak_mat_vess; priznak_pom_vess];
    
    %% Vytvoření matice labelů pro vstup do SVM
    pom = labels_video(klas);
    label_pom = immultiply(pom, maska_rozdil);
    
    label_pom_vess = immultiply(label_pom, g_seg);
    label_pom_vess = reshape(label_pom_vess, [], 1);
    label_pom_vess(label_pom_vess == 0) = [];
    label_mat_vess = [label_mat_vess; label_pom_vess];
    
end