function [marksSLR,ind,messages]=SLR(marks,fa,QRSLIMIT,k,delta,messages)
% post processing rule for delineation
%
% [marksSLR,ind,messages]=SLR(marks,fa,QRSLIMIT,k,delta,messages)
%
% INPUT:
% marks: structure with a matrix of annotations in samples; in each field
%        each matrix of annotations has one line per beat and 10 colunms
%        corresponding respectively to P onset, P peak, P end, QRS onset, QRS
%        main peak, QRS end, T onset, T peak, T prima, T end
% fa: sampling frequency of the recording.
% QRSLIMIT: admited diference between marks of the same beat in samples (optional)
%           by default QRSLIMIT= 250 ms.
% k: number of neighbours to consider in SLR rule (optional)
%    by default k=3 
% delta: neighbourhood to consider in SLR rule in samples (optional)
%    by default delta = [4 nan 4 12 nan 10 12 nan nan 12] ms for 10 marks
%    or  delta = [4 nan nan 4 nan 12  nan nan nan nan 10 12 nan nan 12] ms for 15 marks
% messages:- errors and warnings (optional)
%
% OUTPUT:
% marksSLR: multilead based on marks obtained by SLR rule
% ind: indexes of the initial column vectors for each set, corresponding to the sincronized beats
%      each line of ind corresponds to a beat and each column to a set
%      if a beat is absent of a set
% messages.errors: errors during the processing
% messages.errors_desc:errors description
% messages.status: sucess of the procesing (1) ou fail (0)
% messages.status: parameters considered
% Last update: DEZ2011

marksSLR=[];
ind=[];
if nargin<2
    messages.errors= 'Fatal error in SLR.';
    messages.errors_desc= 'Inputs marks and fa required';
    messages.status=0;
    return
else
    if size(marks{1},2)~=10 && size(marks{1},2)~=15
        messages.errors= 'Fatal error in SLR.';
        messages.errors_desc= 'Input marks should include 10 or 15 marks.';
        messages.status=0;
        return
    end
    messages.status=1;
    if ~isfield(messages,'errors')
        messages.errors=[];
    end
    if ~isfield(messages,'errors_desc')
        messages.errors_desc=[];
    end
    if ~isfield(messages,'warnings')
        messages.warnings=[];
    end
    if ~isfield(messages,'status')
        messages.status=1;
    end
    if ~isfield(messages,'setup')
        messages.setup=[];
    end
    if ~isfield(messages.setup,'SLR')
        messages.setup.SLR=[];
    end
end

if nargin<5
    if size(marks{1},2)==10
        delta = [4 nan 4 12 nan 10 12 nan nan 12]*fa/1000;
    elseif size(marks{1},2)==15
        delta = [4 nan nan 4 nan 12  nan nan nan nan 10 12 nan nan 12]*fa/1000;
    end
    if nargin<4
        k=min(3,size(marks{1},1)-1);
        k=2;
        if nargin<3
            QRSLIMIT=250*fa/1000;
        end
    end
    messages.setup.SLR.QRSlimit=QRSLIMIT;
    messages.setup.SLR.delta=delta;
    messages.setup.SLR.k=k;
    if nargin<3
        if isfield(messages,'setup') & isfield(messages.setup,'SLR')
            if isfield(messages.setup.SLR,'QRSlimit')
                QRSLIMIT=messages.setup.SLR.QRSlimit;
            else
                QRSLIMIT=[];
            end
            if isfield(messages.setup.SLR,'k')
                k=messages.setup.SLR.k;
            else
                k=[];
            end
            if isfield(messages.setup.SLR,'delta')
                if length(delta)==size(marks{1},2)
                    delta= messages.setup.SLR.delta;
                else
                    delta=[];
                end
            end
        else
            QRSLIMIT=[];
            k=[];
            delta=[];
        end
    end
end

if isempty(QRSLIMIT) | ~isnumeric(QRSLIMIT) | isnan(QRSLIMIT)
    QRSLIMIT=250*fa/1000;
    messages.setup.SLR.QRSlimit=QRSLIMIT;
end
if isempty(k) | ~isnumeric(k) | isnan(k)
    k=min(3,size(marks{1},1)-1);
    messages.setup.SLR.k=k;
end
if isempty(delta) | size(delta)~=size(marks{1},2)
    if size(marks{1},2)==10
        delta = [4 nan 4 12 nan 10 12 nan nan 12]*fa/1000;
    elseif size(marks{1},2)==15
        delta = [4 nan nan 4 nan 12  nan nan nan nan 10 12 nan nan 12]*fa/1000;
    end
    messages.setup.SLR.delta=delta;
end
    
if ~isfield(messages.setup.SLR,'refrper')
    if isfield(messages.setup,'wavedet') & isfield(messages.setup.wavedet,'refrper')
        messages.setup.SLR.refrper=messages.setup.wavedet.refrper;
    else
        messages.setup.SLR.refrper=0.275;
    end
end


if ~isfield(messages.setup.SLR,'tolerance')
    messages.setup.SLR.tolerance=10*fa/1000; % 10 msec
end

if isstruct(marks)
    leads = fieldnames(marks);
    messages.setup.nleads=length(leads);
else
   messages.setup.nleads = size(marks,2);
end

if length(delta)==10
    messages.setup.marks_desc={'(','P',')','(','QRS',')','(','T','T''',')'};
%     nmarks=10;
    peaks=[2 5 8 9];
    onsets=[1 4 6];
    ends=[3 6 10];
    refmark=5; %JB 27JUL2011
elseif length(delta)==15
    messages.setup.marks_desc={'(','P','P''',')','QRS_first_peak','(','Q','R','S','R''',')','(','T','T''',')'};
%     nmarks=15;
    peaks=[2 3 8 13 14];
    onsets=[1 6 12];
    ends=[4 11 15];
   % refmark=8; %JB 27JUL2011
   refmark=[8 9 7 10];%RUTE %29JUL2011 porque la onda R puede ni existir en una derivacion especifica!
end
nmarks=length(delta);
beats=0;
% refmark=5;  %JB 27JUL2011
 refmark_n=refmark; %DEZ2011
for g=1:messages.setup.nleads
    if isstruct(marks)
        f = getfield(marks, char(leads(g))); %#ok<GFLD>
    else
        f = marks{g};
    end
    if length(f)> size(beats,1)
        beats((size(beats,1)+1):length(f),1:g)=zeros(length((size(beats,1)+1):length(f)),g);
    end
    % beats(1:size(f,1),g)=f(:,refmark(1));  %19ABRIL2010
    for n=1:size(f,1)%DEZ2011
        refmark(n)=refmark_n(find(~isnan(f(n,refmark_n)),1,'first'));%DEZ2011
        beats(n,g)=f(n,refmark(n)); %29JUL2011 %DEZ2011
    end
  %beats(1:size(f,1),g)=f(:,refmark); %29JUL2011     
end
beats(beats==0)=NaN;
[ind,messages1]=sincronizabeats(beats,QRSLIMIT); %27JUL2011
ind(sum(~isnan(ind),2)<k+1,:) = [];
% nbeats=size(ind,1);

%JB 27JUL2011
beats=NaN*ones(size(ind));
for g=1:size(ind,2)
    if isstruct(marks)
        f = getfield(marks, char(leads(g))); %#ok<GFLD>
    else
        f = marks{g};
    end
    beats(~isnan(ind(:,g)),g)=f(ind(~isnan(ind(:,g)),g),5);
end
%ind(diff(nanmedian(beats,2)) < 200,:) = [];%JB 27JUL2011

ind(find(diff(nanmedian(beats,2)) < messages.setup.SLR.refrper*fa)+1,:) = [];%RUTE 10AGO 2011


nbeats=size(ind,1);
marksSLR=NaN.*ones(nbeats,nmarks);
for mark=peaks
    beats=NaN*ones(size(ind));
    for g=1:size(ind,2)
        if isstruct(marks)
            f = getfield(marks, char(leads(g))); %#ok<GFLD>
        else
            f = marks{g};
        end
        beats(~isnan(ind(:,g)),g)=f(ind(~isnan(ind(:,g)),g),mark);
    end
    marksSLR(:,mark)=nanmedian(beats,2);
end

if  nmarks==15;
    mark=7:10;
    QRS{1}=NaN*ones(size(ind));
    QRS{2}=QRS{1};
    QRS{3}=QRS{1};
    QRS{4}=QRS{1};

    for g=1:size(ind,2)
        if isstruct(marks)
            f = getfield(marks, char(leads(g))); %#ok<GFLD>
        else
            f = marks{g};
        end
        QRS{1}(~isnan(ind(:,g)),g)=f(ind(~isnan(ind(:,g)),g),mark(1)); 
        QRS{2}(~isnan(ind(:,g)),g)=f(ind(~isnan(ind(:,g)),g),mark(2)); 
        QRS{3}(~isnan(ind(:,g)),g)=f(ind(~isnan(ind(:,g)),g),mark(3)); 
        QRS{4}(~isnan(ind(:,g)),g)=f(ind(~isnan(ind(:,g)),g),mark(4)); 
    end
    

    for beat=1:nbeats
        AA=[QRS{1}(beat,:);QRS{2}(beat,:);QRS{3}(beat,:);QRS{4}(beat,:)];
        %look for small R waves that should go with Q waves
        %RS that is an inverted QR or RSR' that is an inverted QRS
        ii=find(AA(2,:)-nanmedian(AA(2,:))<-messages.setup.SLR.tolerance);
        QRS{1}(beat,ii)=QRS{2}(beat,ii); 
        QRS{2}(beat,ii)=QRS{3}(beat,ii); 
        QRS{3}(beat,ii)=QRS{4}(beat,ii); 
        QRS{4}(beat,ii)=NaN; 
        AA=[QRS{1}(beat,:);QRS{2}(beat,:);QRS{3}(beat,:);QRS{4}(beat,:)];

        %look for Q waves that should go with R waves
        %QS complexes or QR that is an inverted RS
        ii=find(AA(1,:)-nanmedian(AA(1,:))>messages.setup.SLR.tolerance);
        QRS{4}(beat,ii)=QRS{3}(beat,ii); 
        QRS{3}(beat,ii)=QRS{2}(beat,ii); 
        QRS{2}(beat,ii)=QRS{1}(beat,ii); 
        QRS{1}(beat,ii)=NaN; 
        AA=[QRS{1}(beat,:);QRS{2}(beat,:);QRS{3}(beat,:);QRS{4}(beat,:)];

        %look for R' waves that should go with S waves
        ii=find(AA(4,:)-nanmedian(AA(4,:))<-messages.setup.SLR.tolerance);
        QRS{1}(beat,ii)=QRS{2}(beat,ii); 
        QRS{2}(beat,ii)=QRS{3}(beat,ii); 
        QRS{3}(beat,ii)=QRS{4}(beat,ii); 
        QRS{4}(beat,ii)=NaN;  
        AA=[QRS{1}(beat,:);QRS{2}(beat,:);QRS{3}(beat,:);QRS{4}(beat,:)];
        marksSLR(beat,mark)=nanmedian(AA');
    end
end

for mark=onsets, %onsets
    beats=NaN.*ones(nbeats,nmarks);
    for g=1:size(ind,2)
        if isstruct(marks)
            f = getfield(marks, char(leads(g))); %#ok<GFLD>
        else
            f = marks{g};
        end
        beats(~isnan(ind(:,g)),g)=f(ind(~isnan(ind(:,g)),g),mark);
    end

    aux= sum(~isnan(beats),2);
    for beat=1:nbeats
        if aux(beat)>k, % if the mark was found in more than k of the leads
            kk = beats(beat,:);
            kk(isnan(kk))=[];
            kk =sort(kk); % sorted the mark found in the leads in beat
            while (length(kk)>k && kk(k+1)>(kk(1)+delta(mark))), %eliminate the first mark until kk satisfy the rule or length(kk)<k
                kk(1)=[];
            end
            %if length(kk)<k there is not a multilead mark
            if length(kk)==k
                marksSLR(beat,mark)=nan;
            else
                marksSLR(beat,mark)=kk(1);
            end
        else
            marksSLR(beat,mark)=nan;
        end
    end

end
for mark=ends, % ends
    beats=NaN.*ones(nbeats,nmarks);
    for g=1:size(ind,2)
        if isstruct(marks)
            f = getfield(marks, char(leads(g))); %#ok<GFLD>
        else
            f = marks{g};
        end
        beats(~isnan(ind(:,g)),g)=f(ind(~isnan(ind(:,g)),g),mark);
    end
    aux= sum(~isnan(beats),2);
    for beat=1:nbeats
        if aux(beat)>k, %
            kk = beats(beat,:);
            kk(isnan(kk))=[];
            kk =flipud(sort(kk)')'; %Rute 19AB2010
            while (length(kk)>k && kk(k+1)<(kk(1)-delta(mark))),
                kk(1)=[];
            end
            if length(kk)==k
                marksSLR(beat,mark)=nan;
            else
                marksSLR(beat,mark)=kk(1);
            end
        else
            marksSLR(beat,mark)=nan;
        end
    end
end