function [positionSLR,marksSLR,ind,messages]=SLR2(marks,messages)
% post processing rule for delineation
%
% [positionSLR,marksSLR,ind,messages]=SLR2(marks,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
% messages.setup.SLR.freq: sampling frequency of the recording.
% messages.setup.SLR.QRSlimit: admited diference between marks of the same beat in ms (optional)
%           by default QRSlimit= 250 ms.
% messages.setup.SLR.k: number of neighbours to consider in SLR rule (optional)
%    by default k=3 
% messages.setup.SLR.delta: neighbourhood to consider in SLR rule in ms (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.setup.wavedet.refrper: minimum time between two admisible beats, by default 257ms
% messages.setup.SLR.tolerance = time in ms to define each individual QRS wave, by default 10ms
% 
% OUTPUT:
%positionSLR: struture with multilead based on marks obtained by SLR rule
% marksSLR: matrix 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: 27JAN2012
positionSLR=[];
marksSLR=[];
ind=[];
if nargin<2
    messages.errors= 'Fatal error in SLR.';
    messages.errors_desc= 'Inputs marks and messages.setup.SLR.freq 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
    if ~isfield(messages,'setup') || ((~isfield(messages.setup,'SLR') ||...
            (~isfield(messages.setup.SLR,'freq') && ~isfield(messages.setup.SLR,'fa')))  &&...
            (~isfield(messages.setup.wavedet,'wavedet') || ( isfield(messages.setup.wavedet,'wavedet') &&...
            ~isfield(messages.setup.wavedet,'freq')))) %
        messages.errors= 'Fatal error in SLR.';
        messages.errors_desc= 'Inputs messages.setup.SLR.freq required';
        messages.status=0;
        return
    else
        try
            fa=messages.setup.SLR.freq;
        catch me
            me.message = 'Not exist messages.setup.SLR.freq';
            try
                fa=messages.setup.SLR.fa;
            catch me
                me.message = 'Not exist messages.setup.SLR.fa';
                fa= messages.setup.wavedet.freq;
            end
        end
    end
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.SLR,'k')
    messages.setup.SLR.k=min(3,size(marks{1},1)-1);
end
if ~isfield(messages.setup.SLR,'QRSlimit')
    messages.setup.SLR.QRSlimit=250;
end
if ~isfield(messages.setup.SLR,'delta')
    if size(marks{1},2)==10
        messages.setup.SLR.delta = [4 nan 4 12 nan 10 12 nan nan 12];
    elseif size(marks{1},2)==15
        messages.setup.SLR.delta = [4 nan nan 4 nan 12  nan nan nan nan 10 12 nan nan 12];
    end
end

QRSlimit= messages.setup.SLR.QRSlimit*fa/1000;
delta= messages.setup.SLR.delta*fa/1000;
k= messages.setup.SLR.k;
if isempty(k) || ~isnumeric(k) || isnan(k)
    k=min(3,size(marks{1},1)-1);
    messages.setup.SLR.k=k;
end
if isempty(QRSlimit) || ~isnumeric(QRSlimit) || isnan(QRSlimit)
    messages.setup.SLR.QRSlimit=250;
    QRSlimit=messages.setup.SLR.QRSlimit*fa/1000;
end
if isempty(delta) | size(delta)~=size(marks{1},2) %#ok<OR2>
    if size(marks{1},2)==10
        messages.setup.SLR.delta = [4 nan 4 12 nan 10 12 nan nan 12];
    elseif size(marks{1},2)==15
        messages.setup.SLR.delta = [4 nan nan 4 nan 12  nan nan nan nan 10 12 nan nan 12];
    end
    delta= messages.setup.SLR.delta*fa/1000;
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*fa/1000;
    else
        messages.setup.SLR.refrper=275*fa/1000;
    end
end
if ~isfield(messages.setup.SLR,'tolerance')
    messages.setup.SLR.tolerance=10*fa/1000;
end

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

if length(delta)==10
    messages.setup.SLR.marks_desc={'(','P',')','(','QRS',')','(','T','T''',')'};
    peaks=[2 5 8 9];
    onsets=[1 4 6];
    ends=[3 6 10];
    refmark=5; 
elseif length(delta)==15
    messages.setup.SLR.marks_desc={'(','P','P''',')','QRS_first_peak','(','Q','R','S','R''',')','(','T','T''',')'};
    peaks=[2 3 8 13 14];
    onsets=[1 6 12];
    ends=[4 11 15];
    refmark=[8 9 7 10];%RUTE % porque la onda R puede ni existir en una derivacion especifica!
end
nmarks=length(delta);
beats=0;
refmark_n=refmark;
 
for g=1:messages.setup.SLR.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
    for n=1:size(f,1)
        refmark(n)=refmark_n(find(~isnan(f(n,refmark_n)),1,'first')); %#ok<AGROW>
        beats(n,g)=f(n,refmark(n));  %#ok<AGROW>
    end     
end

beats(beats==0)=NaN;
[ind,messages1]=sincronizabeats(beats,QRSlimit); 
if messages1.status == 0
    return;
end
ind(sum(~isnan(ind),2)<k+1,:) = [];

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(find(diff(nanmedian(beats,2)) < messages.setup.SLR.refrper)+1,:) = [];

nbeats=size(ind,1);
marksSLR=NaN.*ones(nbeats,nmarks);
marksSLR_lead=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));  %#ok<AGROW>
        QRS{2}(~isnan(ind(:,g)),g)=f(ind(~isnan(ind(:,g)),g),mark(2)); %#ok<AGROW>
        QRS{3}(~isnan(ind(:,g)),g)=f(ind(~isnan(ind(:,g)),g),mark(3)); %#ok<AGROW>
        QRS{4}(~isnan(ind(:,g)),g)=f(ind(~isnan(ind(:,g)),g),mark(4)); %#ok<AGROW>
    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); %#ok<AGROW>
        QRS{2}(beat,ii)=QRS{3}(beat,ii); %#ok<AGROW>
        QRS{3}(beat,ii)=QRS{4}(beat,ii); %#ok<AGROW>
        QRS{4}(beat,ii)=NaN; %#ok<AGROW>
        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); %#ok<AGROW>
        QRS{3}(beat,ii)=QRS{2}(beat,ii); %#ok<AGROW>
        QRS{2}(beat,ii)=QRS{1}(beat,ii); %#ok<AGROW>
        QRS{1}(beat,ii)=NaN; %#ok<AGROW>
        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); %#ok<AGROW>
        QRS{2}(beat,ii)=QRS{3}(beat,ii); %#ok<AGROW>
        QRS{3}(beat,ii)=QRS{4}(beat,ii); %#ok<AGROW>
        QRS{4}(beat,ii)=NaN;  %#ok<AGROW>
        AA=[QRS{1}(beat,:);QRS{2}(beat,:);QRS{3}(beat,:);QRS{4}(beat,:)];
        marksSLR(beat,mark)=nanmedian(AA,2)';
    end
end
if size(marksSLR,2) == 15
    positionSLR.P = marksSLR(:,2);
    positionSLR.Pprima = marksSLR(:,3);
    positionSLR.qrs = marksSLR(:,8);
    positionSLR.QRSon = marksSLR(:,6);
    positionSLR.Q = marksSLR(:,7);
    positionSLR.R = marksSLR(:,8);
    positionSLR.S = marksSLR(:,9);
    positionSLR.Rprima = marksSLR(:,10);
    positionSLR.T = marksSLR(:,13);
    positionSLR.Tprima = marksSLR(:,14);
else
    positionSLR.P = marksSLR(:,2);
    positionSLR.qrs = marksSLR(:,5);
    positionSLR.T = marksSLR(:,8);
    positionSLR.Tprima = marksSLR(:,9);
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 kki] =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)=[];
                kki(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);
                marksSLR_lead(beat,mark)=kki(1);
            end
        else
            marksSLR(beat,mark)=nan;
        end
    end
end
positionSLR.Pon = marksSLR(:,onsets(1));
positionSLR.QRSon = marksSLR(:,onsets(2));
positionSLR.Ton = marksSLR(:,onsets(3));
positionSLR.Pon_lead = marksSLR_lead(:,onsets(1));
positionSLR.QRSon_lead = marksSLR_lead(:,onsets(2));
positionSLR.Ton_lead = marksSLR_lead(:,onsets(3));

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 kki] =sort(kk,'descend');
           % kk =flipud(sort(kk)')'; 
            while (length(kk)>k && kk(k+1)<(kk(1)-delta(mark))),
                kk(1)=[];
                kki(1)=[];
            end
            if length(kk)==k
                marksSLR(beat,mark)=nan;
            else
                marksSLR(beat,mark)=kk(1);
                marksSLR_lead(beat,mark)=kki(1);
            end
        else
            marksSLR(beat,mark)=nan;
        end
    end
end
positionSLR.Poff = marksSLR(:,ends(1));
positionSLR.QRSoff = marksSLR(:,ends(2));
positionSLR.Toff = marksSLR(:,ends(3));
positionSLR.Poff_lead = marksSLR_lead(:,ends(1));
positionSLR.QRSoff_lead = marksSLR_lead(:,ends(2));
positionSLR.Toff_lead = marksSLR_lead(:,ends(3));