% Notes and Recommendations:
% - This code calculates the total area of the object of interest across all frames of the selected videos.
% - The program allows selecting one or multiple videos to be processed in batches.
%   * Note: If more than one video is selected, it is recommended to save the variables of interest before the next video starts processing.
% - For the program to function properly, you must select the Excel file containing the applied stimuli.
%   * The file where the stimuli are saved must have the arrangement as seen in the supplementary guide for the LabVIEW application.
% - The variable "uhm" indicates a threshold where any value above it will be eliminated, in order to remove pixels with saturation or reflections.
% - The variable "arRef" is a threshold used to eliminate any object with an area smaller than this value.
%
% Code Structure:
% The code is divided into the following main stages:
%  1. Selection of videos and the applied stimuli.
%  2. Main loop for sequential processing (video by video):
%     2.1  Search for the applied stimulus corresponding to the video name.
%     2.2  Video file reading.
%     2.3  Filtering of all frames using morphological functions to isolate only the area of the object of interest.
%     2.4  Inversion of the obtained signal (heart area over time) so that minimum values represent relaxation and maximum values represent contraction.
%     2.5  Application of a band-pass filter (3-100 Hz) to suppress background noise and the motion artifact generated by the perfusion drop.
%     2.6  Signal normalization.
%          Note: It is recommended to save the unnormalized signal and later normalize all signals obtained from all videos against the global maximum and minimum of the experiment.
%     2.7  Calculation of the signal's frequency spectrum.
%     2.8  Smoothing of the heart area signal with respect to the freEs variable.
%     2.9  Peak detection in the heart area signal.
%     2.10 Calculation of Area-Area Intervals (AAI).
%     2.11 Calculation of heart rate.
%     2.12 Calculation of Poincaré Maps for all signals.
%          Note: For the FucPoincare(IAA, a3) function to be detected, it must be located in the same folder as this script.
%     2.13 Saving variables ('file3','CoraF','BPM','EstiVec','freEs','t','tE','FreVi')
%          Note: A .mat file is generated in the same folder as the video, using the following format [path,'Area',Na,'.mat'].

clear 
close all
clc
% 1. Prompt user for video(s) and stimulus Excel file
[file,path]=uigetfile('*.avi','Select the video to analyze','MultiSelect','on');
[file2,path2]=uigetfile('*.xlsx','Select the excel with the stimulus');

if iscell(file)
    h5=length(file);% number of videos to process
    a=1;
else
   h5=1;% number of videos to process
   a=0; 
end

BPMT=[]; % vector of frequency values
% Reading the stimulus with respect to the name of the video
x1=readcell([path2,file2],"Sheet","STIMULUS");% Excel table name
hx1=size(x1);
Fm=round(1/x1{2,end});% stimulus acquisition frequency
se = strel('disk',5);
uhm=175;% Maximum reference value to remove artifacts with great intensity, everything above this value will be eliminated
SD1HAT=[]; % Variability values ​​in short times
SD2HAT=[]; % Variability values ​​in long times
arRef=500; % Threshold for small objects to be removed
for i2=1:h5 % 2. Loop to process each video
    ban=0;
    if a>0
       file3=char(string(file(i2))); % Video's name
    else
        file3=file; 
    end

    for j=1:hx1(2)-1 % 2.1 Look for the stimulus related to the video's name.
        rNom=char(x1(1,j));
        if strcmp(file3,rNom)
            EstiVec=cell2mat(x1(2:end,j));
            ban=1;
            h=hx1(1);
        end
    end
    
    % Main cycle
    % 2.2 Read video and compute area-time series by morphological processing
    readerobj = VideoReader([path file3]);  
    disp(['Video in... ',file3,' process']);
    N=readerobj.NumberOfFrames;% number of frames
    FreVi=round(readerobj.FrameRate); % frames per second
    Acora=[];% heart area as a function of time
     for i=1:N % 2.3 Filtering of all frames using morphological functions
        im=read(readerobj, i);
        im=im2gray(im);
        im=double(im);
        we2=find(im>=uhm);
        im(we2)=0;
        imf=imfill(im); % Gap filling function 
        h1=size(im);
        im2=zeros(h1);
        we=find(imf>0);
        im2(we)=255;
        im3 = imbinarize(im2);
        im4 = bwareafilt(im3,arRef); % Eliminates objects with areas less than 500
        im4=imerode(im4,se);
        im4=imdilate(im4,se);
        we3=find(im4>0);
        Acora=[Acora,length(we3)];
    end
    % Filtering and calculation of parameters
    t=(1/FreVi)*(1:1:N); % Video time vector
    tE=(1/Fm)*(0:1:(h-2)); % Stimulus time vector
    % 2.4 Invert area signal
    mxC=max(Acora);
    nCora=-(Acora-mxC);
    % 2.5 Apply a bandpass filter at a frequency from 3 to 100 Hz
    Fc1=3;
    Fc2=100;
    [b,a2] = butter(5,[Fc1 Fc2]/(FreVi/2),'bandpass');
    iref3=filter(b,a2,nCora);
    iref31 = fliplr(iref3);
    iref32= filter(b,a2,iref31);
    CoraF = fliplr(iref32);
    
    % 2.6 Normalizes the stimulus and area signal, this example only normalizes against the maximum and minimum value of the same video
    minC=min(CoraF(FreVi+1:end));
    maxC=max(CoraF(FreVi+1:end));
    AreaF=(CoraF-minC)/(maxC-minC);
    minE=min(EstiVec);
    maxE=max(EstiVec);
    EstF=(EstiVec-minE)/(maxE-minE);
    
    % 2.7 Calculation of the frequency spectrum stimulus
    he=length(EstF);
    xfE=fft(EstF);
    fE=(Fm*(0:((he-1)/2)))/(he-1);
    p2E=abs(xfE/(he-1));
    p1E=p2E(1:round(((he-1)/2)+1));
    p1E(2:end-1)=2*p1E(2:end-1);
    mxP1E=max(p1E(2:end));
    we10=find(p1E==mxP1E);
    f1=0;
    if isempty(we10)
        freEs=6; % Baseline heart rate
        f1=1;
    else
        freEs=fE(we10(1)); % stimulus frequency
    end
    inter=1/round(freEs);
    inCua=round(inter*FreVi); % Window adjustable to the stimulation frequency

    if rem(inCua,2)==0
        inCua2=(inCua)+1;
    else
        inCua2=(inCua);
    end
    AreaF2=sgolayfilt(AreaF,5,inCua2); % 2.8 Signal smoothed with respect to the freEs variable
    % Calculation of the frequency spectrum area
    hA=length(AreaF2);
    xfA=fft(AreaF2);
    fA=(FreVi*(0:((hA-1)/2)))/(hA-1);
    p2A=abs(xfA/(hA-1));
    p1A=p2A(1:round(((hA-1)/2)+1));
    p1A(2:end-1)=2*p1A(2:end-1);
    mxP1A=max(p1A(2:end));
    we11=find(p1A==mxP1A);
    freAs=fA(we11(1));
    
    % 2.9 Calculation of the peaks of each cardiac event
    if ~(f1==1) 
        [Epks,Elocs,Ew,Ep]=findpeaks(EstF,tE,'MinPeakDistance',inter*.6);
    end
    [Apks,Alocs,Aw,Ap]=findpeaks(AreaF2,t,'MinPeakDistance',inter*.6);

    we=find(Alocs<1); %The first second of the signals is eliminated
    Alocs(we)=[];
    Apks(we)=[];
    Alocs(end)=[];
    Apks(end)=[];
    IAA=diff(Alocs); % 2.10 Intervals between area and area
    BPM=60/mean(IAA);% 2.11 Frequency of heart area
    BPMT=[BPMT,BPM];
    
    % Graph of the stimulus, area and detected peaks
    figure
    plot(t,AreaF2)
    hold on
    plot(tE,EstF)
    hold on
    plot(Alocs,Apks,'o')
    if ~(f1==1) 
        hold on
        plot(Elocs,Epks,'o')
    end
    title(string(file(i2)))
       
    % 2.12 Calculation of poicare maps
    a3=1; % If a3 > 0, a plot of the poincare map is generated
    if a3==1
        figure 
        hold on
    end
    [SDA1,SDA2,SAAREA]=FucPoincare(IAA,a3);
    SD1HAT=[SD1HAT,SDA1];
    SD2HAT=[SD2HAT,SDA2];

    if ~(f1==1) 
        we2=find(Elocs<1);
        Elocs(we2)=[];
        Elocs(end)=[];
        IEE=diff(Elocs);
        [SDE1,SDE2,SEAREA]=FucPoincare(IEE,a3);
        Na=num2str(freEs);
    else
        Na='Basal';
    end
    % Poincare map graph of stimulus and area
    if a3==1
        axis([0.06 0.2 0.06 0.2])
        title(file3)
        xlabel('IAA(s)')
        ylabel('IAAt+1(s)')
        hold off
    end
    freEs=freEs*60;
    ID='Area'; % Code Identifier (Atria,Ventricules,Area)
    % 2.13 Saving variables
    save([path,'Area',Na,'.mat'],'file3','ID','CoraF','BPM','EstiVec','freEs','t','tE','FreVi','Fm')
end
