%Fundamentalna funkcja algorytmu wykrywania zalamkow R.
%Zakladamy, ze najwyrazniejszym i tym samym najlatwiejszym do wykrycia
%elementem sygnalu EKG jest odcinek RS. Funkcja ta dobiera bruteforcem
%parametry charakterystyczne dla odcinka RS dla kazdego sygnalu
%indywidualnie oraz oddaje wykryte z ich uzyciem anotacje.

%Wiemy, ze amplituda zespolu QRS w zapisach jakie analizujemy jest 
%w przyblizeniu stala. Dlatego mozemy pozwolic sobie rowniez na okreslenie 
%w jakim zakresie wartosci powinien wystepowac zalamek R (czyli poczatek poszukiwanego przez nas spadku).
%Poszukiwanie wartosci R bylo przydatne kiedys, teraz jest to do
%usuniecia.

%Funkcja dobiera parametry: zakres dlugosci trwania (w probkach)
%odcinka RS oraz w jakim zakresie wartosci moze wystepowac zalamek R (czyli najwieksza wartosc poszukiwanego spadku).
%Funkcja korzysta z innej funkcji: findDeccelerations, aby wykryc w ogole
%wszystkie spadki o zadanych parametrach maksymalnej wartosci oraz dlugosci
%spadku.


function [ANNourREF,lowLengthBest,highLengthBest, fullSkelet,lowAmpBest,highAmpBest] = rsSlopeDetectionEKG(dobrySYGNAL,sfreq, GRANICE, GraniceAMP, ifecg)
%Funkcja pobiera:
%dobrySYGNAL - przygotowany wstepnie sygnal (po usunieciu ruchomej mediany), 
%sfreq - czestotliwosc probkowania, 
%dolnaLOW, dolnaHIGH - zakres poszukiwan optymalnej granicy dolnej czasu trwania odcinka RS 
%(tzn. w jakim zakresie bruteforce bedzie szukal optymalnej dolnej granicy czasu RS).

%gornaLOW, gornaHIGH - zakres poszukiwan optymalnej granicy gornej czasu trwania odcinka RS 
%(tzn. w jakim zakresie bruteforce bedzie szukal optymalnej gornej granicy czasu RS).

%stepLn - krok wyszukiwania czasu RS
%lowg, highg - gorna i dolna granica (w procentach, na podstawie kwantyli)
%wartosci R (jest to czesc do usuniecia, aktualnie nieuzywana). Jej
%usuniecie trzeba zgrac ze zmianami w funkcji findDecelerations.

%Otrzymujemy:
%ANNourREF - Nowe, optymalne wg naszego bruteforca anotacje (fullSkelet po
%removeOutliers).
%lowLengthBest - dla jakiej granicy dolnej dlugosci (w probkach) odcinka RS
%zostaly otrzymane anotacje
%highLengthBest - dla jakiej granicy gornej dlugosci (w probkach) odcinka RS
%zostaly otrzymane anotacje
%fullSkelet - pelny szkielet przed usunieciem smieci


    iterGranice=size(GRANICE,1);
    iterAmp=size(GraniceAMP,1);
    
    celsOfBest2 = cell(iterAmp*iterGranice,6); %inicjalizacja macierzy roboczej. Jest to macierz w ktorej beda zbierane wszystkie wyniki posrednie dzialania bruteforce'a.
    counterek=1; %iterator do chodzenia po macierzy roboczej

    for iiiAmp=1:iterAmp   %bruteforce poszukujacy optymalnej granicy dolnej dlugosci trwania odcinka RS (w probkach)
    
        for iiiTime = 1:iterGranice %bruteforce poszukujacy optymalnej granicy gornej dlugosci trwania odcinka RS (w probkach)
                    
                    Rvalue = QRSParameters(dobrySYGNAL,GraniceAMP(iiiAmp,1),GraniceAMP(iiiAmp,2),GRANICE(iiiTime,1),GRANICE(iiiTime,2),0); %pikifound - wszystkie spadki wykryte przy danym zestawie parametrow
                    
                    %findDecelerations - zewnetrzna funkcja, wyszukuje
                    %wszystkie spadki o zadanych parametrach. Ostatni
                    %parametr (tutaj 0.0) oznacza wymuszenie koniecznosci
                    %przejscia wyszukanego spadku przez ta wskazana wartosc.
                    
                    
                    celsOfBest2(counterek,1)={GRANICE(iiiTime,1)}; %dolna granica dlugosci RS (w probkach)
                    celsOfBest2(counterek,2)={GRANICE(iiiTime,2)}; %gorna granica dlugosci RS (w probkach)
                    celsOfBest2(counterek,3)={Rvalue}; %wszystkie wykryte deceleracje
                    celsOfBest2(counterek,4)={Rvalue}; %wszystkie wykryte deceleracje
                    celsOfBest2(counterek,5)={GraniceAMP(iiiAmp,1)}; %dolna granica dlugosci RS (w probkach)
                    celsOfBest2(counterek,6)={GraniceAMP(iiiAmp,2)}; %gorna granica dlugosci RS (w probkach)
                    counterek=counterek+1;
                    
        end
    
    end
    
    
    
    tabliceSzkieletow2=zeros(size(celsOfBest2(:,1),1),1); %inicjalizacja nowej macierzy roboczej. Beda w niej zbierane tylko dlugosci (!) gotowych szkieletow (po oczyszczeniu za pomoca removeOutliers).
    
    for iter=1:size(celsOfBest2(:,1)) %przejscie po wszystkich otrzymanych zbiorach deceleracji
        
        
        A = cell2mat(celsOfBest2(iter,3)); %zamiana formatu, aby umozliwic dalsze obliczenia
            
        
        if(isempty(A)) %gdy brak wykrytych anotacji nadajacych sie na szkielet, to dlugosc szkieletu = 0
            ANNour=[];
        else
%             ANNour=A(:,1);
                        ANNour = removeOutliers(A(:,1),sfreq); %usuniecie smieci ze zbioru deceleracji w celu uzyskania szkieletu (tzn. "pewnych" deceleracji)
        end
        
        
      
        if(isempty(ANNour)) %gdy brak wykrytych anotacji nadajacych sie na szkielet, to dlugosc szkieletu = 0
            tabliceSzkieletow2(iter)=0;
        else
            tabliceSzkieletow2(iter)=size(ANNour,1); %wpisujemy dlugosc powstalego szkieletu
        end
        celsOfBest2(iter,3)={ANNour};      %wpisanie poprawnych szkieletow do macierzy roboczej
    end    
   

    [~,I] = max(tabliceSzkieletow2); %poszukiwanie najdluzszego (czyt. najlepszego) szkieletu
    lowLengthBest = cell2mat(celsOfBest2(I(1,1),1)); %odczytanie granicy dolnej dlugosci odcinka RS, dla ktorej otrzymalismy ostateczny, najlepszy (najdluzszy) szkielet
    highLengthBest = cell2mat(celsOfBest2(I(1,1),2)); %odczytanie granicy gornej dlugosci odcinka RS, dla ktorej otrzymalismy ostateczny, najlepszy (najdluzszy) szkielet
    lowAmpBest = cell2mat(celsOfBest2(I(1,1),5)); %odczytanie granicy dolnej dlugosci odcinka RS, dla ktorej otrzymalismy ostateczny, najlepszy (najdluzszy) szkielet
    highAmpBest = cell2mat(celsOfBest2(I(1,1),6)); %odczytanie granicy gornej dlugosci odcinka RS, dla ktorej otrzymalismy ostateczny, najlepszy (najdluzszy) szkielet
     %ANNourREF = [];
   % ANNourREF = findDeccelerations(dobrySYGNAL,lowg,highg,lowLengthBest,highLengthBest,0.0);
%     if(ifecg==1)
%         [ANNourREF, Svalue, Qvalue] = QRSParameters(dobrySYGNAL,lowg,highg,lowLengthBest,highLengthBest,0);
%     else
%         [ANNourREF, Svalue, Qvalue] = QRSParameters(dobrySYGNAL,lowg,highg,lowLengthBest,highLengthBest,0);
%     end
    ANNourREF = cell2mat(celsOfBest2(I(1,1),3)); %ostateczna anotacja
        if(isempty(ANNourREF)) %gdy brak wykrytych anotacji nadajacych sie na szkielet, to dlugosc szkieletu = 0
            ANNourREF=[];
        else
            ANNourREF = removeOutliers(ANNourREF(:,1),sfreq); %usuniecie smieci ze zbioru deceleracji w celu uzyskania szkieletu (tzn. "pewnych" deceleracji)
            ANNourREF = ANNourREF(:,1);
        end
    fullSkelet = cell2mat(celsOfBest2(I(1,1),4)); %pelna lista deceleracji (przed removeOutliers) z ktorej powstal ostateczny szkielet
    if(~isempty(fullSkelet))
        fullSkelet = fullSkelet(:,1);
    end
	if(isempty(ANNourREF)) %Kolo ratunkowe na wypadek braku jakichkolwiek wykryc we wczesniejszych krokach.
		ANNourREF =[]; %QRSParameters(dobrySYGNAL,0.93,1.0,dolnaLOW,gornaHIGH,0.0); %Jako ostateczna anotacje wstawiamy wszystkie deceleracje wykryte przy najbardziej liberalnych granicach
    end
    if size(ANNourREF,1) < 3 %Kolo ratunkowe na wypadek super-krotkich wykrytych anotacji ostatecznych.
		ANNourREF =[]; % QRSParameters(dobrySYGNAL,0.93,1.0,dolnaLOW,gornaHIGH,0.0); %Jako ostateczna anotacje wstawiamy wszystkie deceleracje wykryte przy najbardziej liberalnych granicach.
    end
    
end
