% Notes and Recommendations:
% - This code estimates an ellipse that fits as closely as possible to the heart's perimeter, and calculates ventricular shortening by taking the maximum and minimum points of the ellipse along the Y-axis.
% - 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.
% - The variable "a1" indicates whether a representative video of the original video is generated with the calculated ellipse over time; if a1 = 1, this video will be generated.
%   * Note: If the video is generated, the variable "ReFPS" defines the frame rate (FPS) at which the video will be recorded.
% - The variable "canE" indicates how many maximum and minimum values along the axis will be taken; if this value is greater than 1, they will be averaged.
% - The function filterPeaks(Ven, locs, pks, inCua2, FreVi, uhmAm, TiemMin) filters all ventricular events and estimates the amplitude of each event.
%   * Requirement: For it to be detected, it must be located in the same folder as this script.
%   * The variable "uhmAm" is an amplitude threshold that eliminates any event that does not exceed the specified amplitude percentage.
%   * The variable "TiemMin" is a time threshold that eliminates any detected event below 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  Calculation of the ellipse for each frame.
%     2.5  Calculation of the ellipse's maximum and minimum values along the Y-axis.
%     2.6  Calculation of ventricular shortening.
%     2.7  Inversion of the obtained signal (ventricular shortening over time) so that minimum values represent relaxation and maximum values represent contraction.
%     2.8  Application of a band-pass filter (3-100 Hz) to suppress background noise and the motion artifact generated by the perfusion drop.
%     2.9  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.10 Calculation of the signal's frequency spectrum.
%     2.11 Detection of peaks in the signals.
%     2.12 Filtering of peaks using the "filterPeaks" function.
%     2.13 Calculation of Ventricle-Ventricle Intervals (VVI).
%     2.14 Calculation of heart rate.
%     2.15 Saving variables ('file3', 'VenF', 'BPM', 'EstiVec', 'freEs', 'inter', 't', 'tE', 'FreVi', 'Fm').
%          * Note: A .mat file is generated in the same folder as the video, using the following format: [path, 'Ventricles', Na, '.mat'].

clear
close all
clc

% 1 Ask user to pick one or more videos and stimulus Excel
[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

% Read stimulus sheet and prepare parameters
x1=readcell([path2,file2],"Sheet","STIMULUS");% Excel table name
hx1=size(x1);
Fm=round(1/x1{2,end}); %stimulus acquisition frequency
ampli=[]; 
MinAmpli=[];
ampliDer=[];
BPMT=[];

canE=20;% Number of lines to average
uhm=175;% Maximum reference value to remove artifacts with great intensity, everything above this value will be eliminated
se = strel('disk',5);
a1=0; %% Run Video elipcse a1=0 do not record, a1=1 record.
theta = linspace(0,2*pi,200); % Number of points on the entire ellipse
ReFPS=75; % Frequency of the representative video
arRef=500;% Threshold for small objects to be removed
uhmAm=0.1;% Threshold to filter out small peaks
TiemMin=1;% Time threshold to exclude values ​​below this time.
for v=1:h5 % 2 Loop to process each video
    if a>0
       file3=char(string(file(v)));
    else
        file3=file; 
    end
   
    ban=0;
    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

    % 2.2 Main cycle: open video and prepare output
    readerobj = VideoReader([path file3]);  
    disp(['Video in... ',file3,' process']);
    N=readerobj.NumberOfFrames; % number of frames
    FreVi=round(readerobj.FrameRate);% frames per second
    im=read(readerobj, FreVi+1);
    im=im2gray(im);
    h1=size(im);
    VideoObj=VideoWriter([path,filesep,'ViRef',file3(1:end-4),'.avi'],'Motion JPEG AVI');
    VideoObj.Quality=100;
    VideoObj.FrameRate=round(ReFPS); %resampling value 
    if a1==1
        open(VideoObj);
    end 
    Ven=[];% Ventricular shortening over time
    for i=1:N % 2.3 Filtering of all frames using morphological functions
        im=read(readerobj, i);
        im=im2gray(im);
        im2=double(im);
        we2=find(im2>=uhm);
        im2(we2)=0;
        imf=imfill(im2); % Gap filling function 
        h1=size(im2);
        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);
        % Detection of objects leaving only the one with the most area
        s = regionprops(im4,{'Centroid','Orientation','MajorAxisLength','MinorAxisLength','Area'});
        if isempty(s)
            if a1==1
                imshow(uint8(zeros(h1)))
            end
            Ven=[Ven,0];
        else
            % 2.4 Ellipse calculation
            ar=[s(:).Area];
            we3=find(ar==max(ar));
            col = ((s(we3).MajorAxisLength)/2)*cos(theta);
            row = ((s(we3).MinorAxisLength)/2)*sin(theta);
            M = makehgtform('translate',[s(we3).Centroid, 0],'zrotate',deg2rad(-1*s(we3).Orientation));
            D = M*[col;row;zeros(1,numel(row));ones(1,numel(row))];
            if a1==1
                imshow(uint8(im))
                hold on
                plot(D(1,:),D(2,:),'r','LineWidth',2)
                hold off
                FVE=getframe(gcf);
                writeVideo(VideoObj,FVE);
            end
            % 2.5 Calculation of maximum and minimum points
            hd=size(D);
            dvec=[];
            for i2=1:hd(2)/2
                if i2==1
                    cony=1;
                else
                    cony=hd(2)-(i2-1);
                end
                conx=(hd(2)/2)+i2;
                dvec=[dvec,sqrt(((D(2,cony)-D(2,conx))^2)+((D(1,cony)-D(1,cony))^2))];
            end
            VenRef=[];
            D3=D;
            for i3=1:canE
                % 2.6 Calculation of ventricular shortening
                xmin=min(D3(2,:));
                we2=find(D3(2,:)==xmin);
                xymin=D3(1,we2(1));
                xmax=max(D3(2,:));
                we2=find(D3(2,:)==xmax);
                xymax=D3(1,we2(1));
                D3(:,we2(1))=[];
                VenRef=[VenRef,sqrt(((xmax-xmin)^2)+((xymax-xymin)^2))];
            end
            Ven=[Ven,mean(VenRef)];
        end
    end
    if a1==1
        close(VideoObj);
    end
    t=(1/FreVi)*(1:1:N); % Video time vector
    tE=(1/Fm)*(0:1:(h-2));% Stimulus time vector
    % Filtering and calculation of parameters
    % 2.7 Invert area signal
    Ven2=Ven-max(Ven);
    Ven2=-Ven2;
    
    % 2.8 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,Ven2);
    iref31 = fliplr(iref3);
    iref32= filter(b,a2,iref31);
    VenF = fliplr(iref32);
    
    % 2.9 Normalizes the stimulus and area signal, this example only normalizes against the maximum and minimum value of the same video
    minV=min(VenF(FreVi+1:end));
    maxV=max(VenF(FreVi+1:end));
    VenF2=(VenF-minV)/(maxV-minV);
    
    % 2.10 Calculation of the frequency spectrum ventricules
    hV=length(VenF2);
    xfV=fft(VenF2);
    fV=(FreVi*(0:((hV-1)/2)))/(hV-1);
    p2V=abs(xfV/(hV-1));
    p1V=p2V(1:round(((hV-1)/2)+1));
    p1V(2:end-1)=2*p1V(2:end-1);
    
    mxP1V=max(p1V(2:end));
    we11=find(p1V==mxP1V);
    freVs=fV(we11(1));
    
    % Calculation of the frequency spectrum stimulus
    Fm=round(1/x1{2,end});
    EstiVec=EstiVec-min(EstiVec);
    minE=min(EstiVec);
    maxE=max(EstiVec);
    EstF=(EstiVec-minE)/(maxE-minE);
    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);
    na="";
    f1=0;
    if isempty(we10)
        freEs=6;% Baseline heart rate
        f1=1;
    else
        freEs=fE(we10(1)); %stimulus frequency
    end
    inter=1/round(freEs);
    
    % 2.11 Calculate peaks in stimulus and ventricular signals
    if ~(f1==1)  
        [Epks,Elocs,Ew,Ep]=findpeaks(EstF,tE,'MinPeakDistance',inter*.6);
    end
    [Vpks,Vlocs,Vw,Vp]=findpeaks(VenF2,t,'MinPeakDistance',inter*.6);
    inCua=round(inter*FreVi); % Window adjustable to the stimulation frequency
    if rem(inCua,2)==0
        inCua2=(inCua)+1;
    else
        inCua2=(inCua);
    end
    
    % 2.12 Calculation of the amplitude of each peak and the ventricular interval
    [locV,vecAmpV]=filterPeaks(VenF2,Vlocs,Vpks,inCua2,FreVi,uhmAm,TiemMin);
    Vlocven3=Vlocs(locV);
    Vpks3=Vpks(locV);
    promA=mean(vecAmpV(:,1));
    promAD=mean(vecAmpV(:,2));
    minA=min(vecAmpV(:,1));
    ampli=[ampli,promA]; % Left width of each average video
    ampliDer=[ampliDer,promAD]; % Right width of each average video
    MinAmpli=[MinAmpli,minA]; % Minimum width of each average video
    
    if ~(f1==1)
        we3=find(Elocs<TiemMin);
        Elocs(we3)=[];
        Epks(we3)=[];
        IEE=diff(Elocs);
    end
    IVV=diff(Vlocven3); % 2.13 Intervals between ventricles and ventricles event
    BPM=60/mean(IVV); % 2.14 Frequency of ventricles
    BPMT=[BPMT,BPM];
    
    % Graph of the stimulus, ventricles and detected peaks
    figure
    plot(t,VenF2)
    hold on
    plot(Vlocven3,Vpks3,'o')
    freEs=freEs*60;
    if ~(f1==1)
        hold on
        plot(tE,EstF)
        hold on
        plot(Elocs,Epks,'o')
        legend("Ventricles","Stimulus","Peaks ventricles","Peaks Stimilus")
        Na=num2str(freEs);
    else
        Na='Basal';
    end
    ID='Ventricles'; % Code Identifier (Atria,Ventricules,Area)
    % 2.15 Saving variables
    save([path,'Ventricles',Na,'.mat'],'ID','file3','VenF','BPM','EstiVec','freEs','inter','t','tE','FreVi','Fm')
end





