% Notes and Recommendations:
% - This code calculates the area of each atrium over time and obtains an average of the atrial response.
% - 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 "mask" is a flag. If "mask" == 0, the two ROIs (Regions of Interest) are manually selected during the first second of video 1. If "mask" > 0, two previously saved ROIs are selected.
% - The function filterPeaks(Ven, locs, pks, inCua2, FreVi, uhmAm, TiemMin) filters all atrial events and estimates the amplitude of each one.
%   * 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.2.1 If "mask" == 0: The first selected video is used to choose two ROIs (one for each atrium); these will be applied to all other videos.
%       2.2.2 If "mask" > 0: Two previously saved ROIs are selected to be used across all selected videos.
%     2.3 Filtering of each ROI individually across all frames using morphological functions to isolate only the area of the object of interest.
%     2.4 Calculation of the object of interest's area, considering any pixel with a value greater than 0 as part of the area.
%     2.5 Calculation of the frequency spectrum of the applied stimulus.
%     2.6 Application of a band-pass filter (3-100 Hz) to suppress background noise and the motion artifact generated by the perfusion drop.
%     2.7 Inversion of the obtained signal (atrial area) so that minimum values represent relaxation and maximum values represent contraction.
%     2.8 Averaging of the two atrial signals: "AFS" = ("AFS1" + "AFS2") / 2;
%     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 atrial 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 Atrium-Atrium Intervals ("IAT").
%     2.14 Calculation of heart rate.
%     2.15 Saving variables ('file3', 'NorA', 'AF1', 'AF2', 'BPM', 'Estivec', 'Fe', 'inter', 't', 'tE', 'FreVi', 'Fm', 'b1', 'b2').
%       * Note: A .mat file is generated in the same folder as the video, using the following format: [path, 'Atria', 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');
file3=file;

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

% 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
BPMT=[];
MinAmpli=[];
ampliDer=[];
ampli=[];
uhm=175;% Maximum reference value to remove artifacts with great intensity, everything above this value will be eliminated
se = strel('disk', 1);
mask=0; % Flag indicating whether new ROIs are generated or one from another video is used; if "mask" > 1, the mask from another video is used as a reference.
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
    readerobj = VideoReader([path file3]);  
    disp(['Video in... ',file3,' process']);
    N=round(readerobj.NumberOfFrames);
    FreVi=round(readerobj.FrameRate);
    if v==1 % Use the first selected video as a reference to choose the ROIs.
        im=read(readerobj, FreVi+1);
        im=im2gray(im);
        h1=size(im);
        if mask==0 % 2.2.1. Selecting the ROIs
            % Reference image to place the roi
            figure(1)
            imshow(im)
            title(file3)
            % Select the ROI
            a1=drawellipse;
            pause
            disp('Atrium 1')
            a2=drawellipse;
            pause
            disp('Atrium 2')
            b1=createMask(a1);% Roi mask 1
            b2=createMask(a2);% Roi mask 2
        else % 2.2.2. Select ROIs from another video
            [file4,path4]=uigetfile('*.mat','Select the .mat file containing the masks.');
            load([path4,file4],'b1','b2');
        end
        b3=b1+b2; % Mask together
        we=find(b3==0);
        imref=im;
        imref(we)=0;
        % Reference image of ROI
        figure(2)
        imshow(imref)
        title(file3)
        pause(1)
        we2=find(b1==0);
        we3=find(b2==0);
    end
    
    A1=[];% atrium vector 1
    A2=[];% atrium vector 2
    for i=1:N % 2.3 Filtering of all frames using morphological functions
        im=read(readerobj,i);
        im=im2gray(im);
        im=double(im);
        imA1=im;
        imA2=im;
        imA1(we2)=0;
        imA2(we3)=0;

        % ROI Filtering 1
        BA1 = imsharpen(imA1);
        B3A1=AumCont(BA1,1); % Contrast enhancement function
        B2A1=edge(B3A1,"canny");
        we4=find(B2A1>0);
        A1=[A1,length(we4)]; % 2.4. Calculation of the area
    
        % ROI Filtering 2
        BA2 = imsharpen(imA2);
        B3A2=AumCont(BA2,1);
        B2A2=edge(B3A2,"canny");
        we5=find(B2A2>0);
        A2=[A2,length(we5)];% 2.4. Calculation of the area

    end
    
    % Filtering and calculation of parameters
    t=(1/FreVi)*(1:1:N);
    tE=(1/Fm)*(0:1:(h-2));
    % 2.5. 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);
    f1=0;
    if isempty(we10)
        Fe=6; % Baseline heart rate
        f1=1;
    else
        Fe=fE(we10(1)); %stimulus frequency 
    end
    inter=1/round(Fe);
    inCua=round(inter*FreVi); % Window adjustable to the stimulation frequency
    
    if rem(inCua,2)==0
        inCua2=(inCua)+1;
    else
        inCua2=(inCua);
    end
    
    % 2.6. 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,A1);
    iref31 = fliplr(iref3);
    iref32= filter(b,a2,iref31);
    AF1 = fliplr(iref32); 
    AFS1=sgolayfilt(AF1,5,inCua2); % Signal smoothed with respect to the Fe variable

    iref3=filter(b,a2,A2);
    iref31 = fliplr(iref3);
    iref32= filter(b,a2,iref31);
    AF2 = fliplr(iref32);
    AFS2=sgolayfilt(AF2,5,inCua2); % Signal smoothed with respect to the Fe variable
    
    % 2.7. Invert the signal of the event
    maxA1=max(AFS1);
    AFS1=-(AFS1-maxA1);
    maxA2=max(AFS2);
    AFS2=-(AFS2-maxA2);
    
    AFS=(AFS1+AFS2)/2; % 2.8. Prom Atria
    
    
    % 2.9. Normalizes the stimulus and area signal, this example only normalizes against the maximum and minimum value of the same video
    minA=min(AFS(FreVi+1:end));
    maxA=max(AFS(FreVi+1:end));
    NorA=(AFS-minA)/(maxA-minA);
    
    % 2.10. Calculation of the frequency spectrum atria
    hA=length(NorA);
    xfA=fft(NorA);
    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.11. Peak detection for stimulus and atria signals
    if ~(f1==1)
        [Epks,Elocs,Ew,Ep]=findpeaks(EstF,tE,'MinPeakDistance',inter*.6);
    end
    [Atpks,Atlocs,Atw,Atp]=findpeaks(NorA,t,'MinPeakDistance',inter*.6);

    % 2.12. Calculation of the amplitude of each peak and the atria interval
    [locAt,vecAmpAt]=filterPeaks(NorA,Atlocs,Atpks,inCua2,FreVi,uhmAm,TiemMin);
    Atlocven3=Atlocs(locAt);
    Atpks3=Atpks(locAt);
    promA=mean(vecAmpAt(:,1));
    promAD=mean(vecAmpAt(:,2));
    minA=min(vecAmpAt(:,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<1);
        Elocs(we3)=[];
        Epks(we3)=[];
        IE=diff(Elocs);
    end
    IAT=diff(Atlocven3); % 2.13. Intervals between atria and atria event
    BPM=60/mean(IAT); % 2.14. Frecuency of atria
    BPMT=[BPMT,BPM];
    Fe=Fe*60;
    % Graph of the stimulus, atria and detected peaks
    figure
    plot(t,NorA)
    hold on
    plot(Atlocven3,Atpks3,'o')
    if ~(f1==1)
        hold on
        plot(tE,EstF)
        hold on
        plot(Elocs,Epks,'o')
        legend("Atria","Peaks atria","Stimulus","Peaks Stimilus")
        Na=num2str(Fe);
    else
        Na='Basal';
    end
    ID='Atria'; % Code Identifier (Atria,Ventricules,Area)
    % 2.15. Save results and key variables for this video
    save([path,'Atria',Na,'.mat'],'ID','file3','NorA','AFS1',"AFS2",'BPM','EstiVec','Fe','inter','t','tE','FreVi','Fm','b1','b2')
end
