%% Skript sloužící pro analýzu pulzací cévního řečiště
% Analýza sestává z načtení videa a definice oblasti zájmu ve snímcích, kde
% bude analýza pulzací provedena. Z cév v oblasti zájmu jsou pak získány
% průměrné intenzity. Tyto průměrné intenzity cév ze všech snímků videa
% jsou zobrazeny v čase.

clc; clear

%% Načtení cest k datům
% videa
vid_location = 'C:\Users\jana-\Documents\FEKT\DP\data\pulzace\analýza\videos';   % cesta ke videím
addpath(vid_location);
vid_data = fileDatastore(vid_location, 'ReadFcn', @load, 'FileExtensions', '.avi');
pocet = length(vid_data.Files);                                                  % počet dat

% průměr
av_location = 'C:\Users\jana-\Documents\FEKT\DP\data\pulzace\analýza\averages';  % cesta ke zprům. videím
addpath(av_location);
av_data = imageDatastore(av_location);

% segmentační masky
seg_location = 'C:\Users\jana-\Documents\FEKT\DP\data\pulzace\analýza\vessels';  % cesta k maskám segmentovaných cév
addpath(seg_location);
seg_data = imageDatastore(seg_location);

% body OD
OD_location = 'C:\Users\jana-\Documents\FEKT\DP\data\pulzace\analýza\masks';     % cesta k bodům OD
addpath(OD_location);
OD_data = fileDatastore(OD_location, 'ReadFcn', @load, 'FileExtensions', '.mat');

%% Načtení dat vybraného videa
video = VideoReader(vid_data.Files{3});
seg = imread(seg_data.Files{3});
OD = load(OD_data.Files{3});

%% Zpracování videa
FR = video.FrameRate;
nFrames = video.NumFrames;
prum_int = [];
for j = 1:nFrames
    frame = read(video, j);
    data_frame(:, :, j) = frame(:, :, 2); 

    snimek = data_frame(:, :, j);     % j-tý snímek videa

    %% Zisk středu OD
    souradnice = OD.vyslednevrcholy2;
    bw = poly2mask(souradnice(:, 1), souradnice(:, 2), size(snimek, 1), size(snimek, 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

    %% 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(snimek);

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

    %% Zobrazení dat
    % figure
    % subplot 121; imshow(snimek, []); hold on; plot(souradnice(:, 1), souradnice(:, 2), 'r.', 'LineWidth', 2, 'MarkerSize', 15); hold on; plot(souradnice_x, souradnice_y, 'r+', 'LineWidth', 1, 'MarkerSize', 7);
    % subplot 122; imshow(bw); hold on; plot(souradnice_x, souradnice_y, 'r+', 'LineWidth', 1, 'MarkerSize', 7);
    % 
    % figure
    % imshowpair(snimek, maska_rozdil, 'blend'); title('Oblast zájmu okolo OD');

    %% Násobení masek
    [r, g_seg, b] = imsplit(seg);
    g_seg(g_seg == 255) = 1;                       % změna hodnot v segmentační masce na 1
    snimek_cevy = immultiply(snimek, g_seg);       % násobení obrazu se segmentační maskou - zisk pouze cév
    snimek_cevy_ROI = immultiply(im2double(snimek_cevy), maska_rozdil);

    %% Zisk průměrné intenzity z jednotlivých úseků cév
    pom = immultiply(double(g_seg), maska_rozdil);
    [labeledImage, n] = bwlabel(pom);
    prum_int_pom = regionprops(labeledImage, snimek_cevy_ROI, 'MeanIntensity'); % průměrná intenzita každého úseku cévy
    prum_int = [prum_int, prum_int_pom];

    % zobrazení pomocí labeloverlay
    % C = labeloverlay(snimek_cevy_ROI, labeledImage);
    % figure; imshow(C, []); title('Zobrazení labelů pixelů');
    
end

%% Zpracování vybraných signálů
vislabels(labeledImage)                   % zobrazení výsledku fce bwlabel s přiřazeným číslem
% funkce vislabels citace: Steve Eddins (2021). Visualize output of BWLABEL (https://www.mathworks.com/matlabcentral/fileexchange/19665-visualize-output-of-bwlabel), MATLAB Central File Exchange. Retrieved May 17, 2021.

matrix = cell2mat(struct2cell(prum_int));
matrix = squeeze(matrix(1, :, :));        % uložení signálů do jedné matice

%% Zobrazení vybraných signálů

signal_tepna = matrix(17, :);             % výběr signálu tepny
signal_zila = matrix(14, :);              % výběr signálu žíly

c = corrcoef(signal_tepna, signal_zila)   % korelační koeficient

[envHigh_tepna, envLow_tepna] = envelope(signal_tepna, 1, 'peak');
[envHigh_zila, envLow_zila] = envelope(signal_zila, 1, 'peak');

[pk, loc] = findpeaks(envHigh_tepna, 'MinPeakDistance', 10);

figure
hold on
plot(signal_tepna, 'r');
plot(signal_zila, 'b');
hold off
xlabel('Počet videosekvencí'); ylabel('Průměrná intenzita');
grid on
grid minor
    