function [position0,position,messages]= multiqrson(heasig,w1,w2,w3,auxqrspiconset,QRSon,time,timenew,position0,position,intervalo,samp,messages)
% [position0,position]= multiqrson(heasig,w1,w2,w3,auxqrspiconset,QRSon,time,timenew,position0,position,intervalo,position1,position2,position3,samp)
%
% INPUT:
% heasig: header information
% w1,w2,w3: WT scales 1 to 3
% auxqrspiconset: earliest relevant modulus maximum location associated to QRS complex found in any lead
% QRSon:earliest QRS onset location in any lead (initial QRS onset location)
% time: beats processed in this segment in each lead
% timenew: beats processed in this segment (multilead)
% position: struct vector with the detected points at the first multilead
% step
% position: struct vector with the detected points
% intervalo: numeration of the beats processed in this segment
% position1,position2,position3: struct vectors with the detected points using each specified lead
% samp: samples included in the current excerpt (borders excluded)
%
%
% OUTPUT:
% actualized position0,position
%
% Rute Almeida 14.APR.2005 (based on multiqrs)
% Last update: Rute %05AGO2011
global recursioncount OPT
spf=heasig.spf_ecg; %24.03.09
if ~isfield(messages.setup.wavedet,'Kron')
messages.setup.wavedet.Kron = 20;% 3.5 5 8 10 15 20!
end
if ~isfield(messages.setup.wavedet,'Kroff')
messages.setup.wavedet.Kroff = 14;
end
if ~isfield(messages.setup.wavedet,'Kq')% Q onset, antiguo Kq==1; 2 2.5 3.5 5 10 15!
messages.setup.wavedet.Kq = 15;
end
if ~isfield(messages.setup.wavedet,'Ks')
messages.setup.wavedet.Ks = 8; % S offset (Proyeuto de Mapi) . 2.5 3 5 8;
end
if ~isfield(messages.setup.wavedet,'maxQRSon')
messages.setup.wavedet.maxQRSon =6.5*2; % based on CSE tolerance
end
% if ~isfield(messages.setup.wavedet,'tolsamepeak')
% messages.setup.wavedet.tolsamepeak=5; %%%% number in samples to consider that a different peak was chosen
% end
%tolsamepeak=messages.setup.wavedet.tolsamepeak;
% Kron=messages.setup.wavedet.Kron;
% Kroff=messages.setup.wavedet.Kroff;
% Kq=messages.setup.wavedet.Kq;
% Ks=messages.setup.wavedet.Ks;
maxQRSon=messages.setup.wavedet.maxQRSon;
auxqrson=nanmin(QRSon-samp(1)+1)-round(maxQRSon/1000*heasig.freq*spf);
%auxqrson=min(QRSon-samp(1)+1)-round(maxQRSon/1000*heasig.freq*spf); 28Abril06
for i=1:size(time,2)
position.QRSonsetcriteria((i+intervalo(1)-1))=2;% default % 1 not considered
scale=2;
if size(w1,2)==3
pontos=[w1(auxqrson(i):auxqrspiconset(i),scale) w2(auxqrson(i):auxqrspiconset(i),scale) w3(auxqrson(i):auxqrspiconset(i),scale)];
else
pontos=[w1(auxqrson(i):auxqrspiconset(i),scale) w2(auxqrson(i):auxqrspiconset(i),scale)];
end
%weight=ones(length(pontos),1); %RUTE 24MAR
weight=ones(size(pontos,1),1);
[a,v] = optimline(pontos, weight,OPT); %#ok<ASGLU>
%time0=max(time);% Rute 18.07.2007
%time0=round(nanmean(time));% Rute 18.07.2007
time0=round(nanmedian(time));% Rute 01.03.2010
%waux=w1; %Rute6NOV08
waux=nan*ones(size(w1)); %Rute6NOV08
if i==1
if size(timenew,2)>1 %4JUL08
intervalonewlead=1:timenew(i+1);
else
intervalonewlead=1:length(w1);
end
elseif i<length(timenew)
intervalonewlead=timenew(i-1):timenew(i+1);
else
intervalonewlead=timenew(i-1):length(w1);
end
if size(w1,2)==3
newleadbeat1=(w1(intervalonewlead,1).*v(1)+w2(intervalonewlead,1)*v(2)+w3(intervalonewlead,1)*v(3))./norm(v);
newleadbeat2=(w1(intervalonewlead,2).*v(1)+w2(intervalonewlead,2)*v(2)+w3(intervalonewlead,2)*v(3))./norm(v);
%signew=(sig(intervalonewlead,1).*v(1)+sig(intervalonewlead,2)*v(2)+sig(intervalonewlead,3)*v(3))./norm(v);
else
newleadbeat1=(w1(intervalonewlead,1).*v(1)+w2(intervalonewlead,1)*v(2))./norm(v);
newleadbeat2=(w1(intervalonewlead,2).*v(1)+w2(intervalonewlead,2)*v(2))./norm(v);
%signew=(sig(intervalonewlead,1).*v(1)+sig(intervalonewlead,2)*v(2))./norm(v);
end
waux(intervalonewlead,1)=newleadbeat1;
waux(intervalonewlead,2)=newleadbeat2;
[position0,qrspiconset0,qrspicoffset0]= qrswavef(heasig,samp,time0(i),position0,waux,[i+intervalo(1)-1 i+intervalo(1)-1],messages); %#ok<NASGU>
if isempty(qrspiconset0)||isnan(qrspiconset0) % RUTE 6ABRI09
% piconall= NaN;
% amppiconall= NaN;
% amppiconaux= NaN;
recursioncount.QRSon((i+intervalo(1)-1))=-2;
position.QRSon(i+intervalo(1)-1) = NaN;
position.Q(i+intervalo(1)-1) = NaN;
position_safe.R(i+intervalo(1)-1) = NaN; %13.AGO.07
position.R_inQRSon(i+intervalo(1)-1) = NaN;
else
if size(w1,2)==3
pontos2=[w1((position0.QRSon(i+intervalo(1)-1)-samp(1)+1-round(maxQRSon/1000*heasig.freq*spf)):qrspiconset0,scale) w2((position0.QRSon(i+intervalo(1)-1)-samp(1)+1-round(maxQRSon/1000*heasig.freq*spf)):qrspiconset0,scale) w3((position0.QRSon(i+intervalo(1)-1)-samp(1)+1)-round(maxQRSon/1000*heasig.freq*spf):qrspiconset0,scale)];
else
pontos2=[w1((position0.QRSon(i+intervalo(1)-1)-samp(1)+1-round(maxQRSon/1000*heasig.freq*spf)):qrspiconset0,scale) w2((position0.QRSon(i+intervalo(1)-1)-samp(1)+1-round(maxQRSon/1000*heasig.freq*spf)):qrspiconset0,scale)];
end
weight=ones(length(pontos2),1);
[a2,v2] = optimline(pontos2, weight,OPT); %#ok<ASGLU>
if i==1
if size(timenew,2)>1 %4JUL08
intervalonewlead=1:timenew(i+1);
else
intervalonewlead=1:length(w1);
end
elseif i<length(timenew)
intervalonewlead=timenew(i-1):timenew(i+1);
else
intervalonewlead=timenew(i-1):length(w1);
end
if size(w1,2)==3
newleadbeat12=(w1(intervalonewlead,1).*v2(1)+w2(intervalonewlead,1)*v2(2)+w3(intervalonewlead,1)*v2(3))./norm(v2);
newleadbeat22=(w1(intervalonewlead,2).*v2(1)+w2(intervalonewlead,2)*v2(2)+w3(intervalonewlead,2)*v2(3))./norm(v2);
%signew2=(sig(intervalonewlead,1).*v2(1)+sig(intervalonewlead,2)*v2(2)+sig(intervalonewlead,3)*v2(3))./norm(v2);
else
newleadbeat12=(w1(intervalonewlead,1).*v2(1)+w2(intervalonewlead,1)*v2(2))./norm(v2);
newleadbeat22=(w1(intervalonewlead,2).*v2(1)+w2(intervalonewlead,2)*v2(2))./norm(v2);
%signew2=(sig(intervalonewlead,1).*v2(1)+sig(intervalonewlead,2)*v2(2))./norm(v2);
end
%waux=w1; %Rute6NOV08
waux=nan*ones(size(w1)); %Rute6NOV08
waux(intervalonewlead,1)=newleadbeat12;
waux(intervalonewlead,2)=newleadbeat22;
% positionaux=position; %RUTE 11.07.2007
%[position,qrspiconset2,qrspicoffset2]= qrswavef(heasig,samp,time0(i),position,waux,[i+intervalo(1)-1 i+intervalo(1)-1],Kron,Kroff,Kq,Ks);
[positionaux,qrspiconset2,qrspicoffset2]= qrswavef(heasig,samp,time0(i),position,waux,[i+intervalo(1)-1 i+intervalo(1)-1],messages); %#ok<NASGU>
%%%%%%%%%%%%%%%%%%%%%%%RUTE 11.07.2007
position.QRSon(i+intervalo(1)-1) = positionaux.QRSon(i+intervalo(1)-1);
position.Q(i+intervalo(1)-1) = positionaux.Q(i+intervalo(1)-1);
position_safe.R(i+intervalo(1)-1) = positionaux.R(i+intervalo(1)-1) ;%13.AGO.07
position.R_inQRSon(i+intervalo(1)-1) = positionaux.R(i+intervalo(1)-1) ;
%%%%%%%%%%%%%%%%%%%%%%%
%if isempty(qrspiconset2) % RUTE 3ABRI09
if isempty(qrspiconset2)||isnan(qrspiconset2) % RUTE 3ABRI09
piconall=[qrspiconset0 NaN];
amppiconall=[(newleadbeat2(qrspiconset0-intervalonewlead(1)+1)) NaN];
amppiconaux= NaN;
else
piconall=[qrspiconset0 qrspiconset2];
amppiconall=[(newleadbeat2(qrspiconset0-intervalonewlead(1)+1)) (newleadbeat22(qrspiconset2-intervalonewlead(1)+1))];
amppiconaux= newleadbeat2(qrspiconset2-intervalonewlead(1)+1);
end
newQRSon=[NaN position0.QRSon(i+intervalo(1)-1) position.QRSon(i+intervalo(1)-1)];
newleadbeat22aux=newleadbeat22;
% stop criteria iii) | iv)) ou novo %13.AGO.07
%if position.QRSon(i+intervalo(1)-1)>=position.R(i+intervalo(1)-1)&((position0.QRSon(i+intervalo(1)-1)<position.R(i+intervalo(1)-1)))|((abs(amppiconall(end))<abs(amppiconall(end-1)))&((amppiconaux(end)*amppiconall(end))>0) & ~isnan(position0.Q(i+intervalo(1)-1)))|( isnan(piconall(end))) % RUTE 9.07.2007
% stop criteria iii) | iv)) % opçao0 % parte da alternativa2
if ((abs(amppiconall(end))<abs(amppiconall(end-1)))&&((amppiconaux(end)*amppiconall(end))>0) &&...
~isnan(position0.Q(i+intervalo(1)-1)))||( isnan(piconall(end))) % RUTE 9.07.2007
% stop criteria iii) | iv | i))) % alternativa1
%if ((abs(amppiconall(end))<abs(amppiconall(end-1)))&((amppiconaux(end)*amppiconall(end))>0) & ~isnan(position0.Q(i+intervalo(1)-1)))|( isnan(piconall(end))) | (abs(newQRSon(end)-newQRSon(end-1))<=1) % RUTE 9.07.2007
%alternative criteria stop criteria iii) without Q condition|iv)alternativa3
%if ((abs(amppiconall(end))<abs(amppiconall(end-1)))&((amppiconaux(end)*amppiconall(end))>0))|( isnan(piconall(end))) % | (abs(newQRSon(end)-newQRSon(end-1))<=1) % RUTE 9.07.2007
%alternative criteria stop criteria just iv)
%if ( isnan(piconall(end))) % | (abs(newQRSon(end)-newQRSon(end-1))<=1) % RUTE 9.07.2007
if isnan(piconall(end))
position.QRSonsetcriteria((i+intervalo(1)-1))=4;
elseif (abs(newQRSon(end)-newQRSon(end-1))<=1)% alternativa2 - exclude criteria i)
position.QRSonsetcriteria((i+intervalo(1)-1))=1;
else
position.QRSonsetcriteria((i+intervalo(1)-1))=3;
end
recursioncount.QRSon((i+intervalo(1)-1))=-1;
position.QRSon(i+intervalo(1)-1) = position0.QRSon(i+intervalo(1)-1);
position.Q(i+intervalo(1)-1) = position0.Q(i+intervalo(1)-1);
position_safe.R(i+intervalo(1)-1) = position0.R(i+intervalo(1)-1) ; %13.AGO.07
position.R_inQRSon(i+intervalo(1)-1) = position0.R(i+intervalo(1)-1);
else
recursioncount.QRSon((i+intervalo(1)-1))=0;
% stop criteria iv) | ii) | i) % opçao0 % parte da alternativa1
%while ~isnan(piconall(end)) & nansum(newQRSon(end)-newQRSon(1:end-1)==0)<2 & (abs(newQRSon(end)-newQRSon(end-1))>1 )%%%
%alternative criteria stop criteria iv) | ii)) % alternativa2e3
while ~isnan(piconall(end)) && nansum(newQRSon(end)-newQRSon(1:end-1)==0)<2 %& (abs(newQRSon(end)-newQRSon(end-1))>1 )%%%
recursioncount.QRSon((i+intervalo(1)-1))=recursioncount.QRSon((i+intervalo(1)-1))+1;
if size(w1,2)==3
pontosnew=[w1((position.QRSon(i+intervalo(1)-1)-samp(1)+1-round(maxQRSon/1000*heasig.freq*spf)):qrspiconset2,scale) w2((position.QRSon(i+intervalo(1)-1)-samp(1)+1-round(maxQRSon/1000*heasig.freq*spf)):qrspiconset2,scale) w3((position.QRSon(i+intervalo(1)-1)-samp(1)+1)-round(maxQRSon/1000*heasig.freq*spf):qrspiconset2,scale)];
else
pontosnew=[w1((position.QRSon(i+intervalo(1)-1)-samp(1)+1-round(maxQRSon/1000*heasig.freq*spf)):qrspiconset2,scale) w2((position.QRSon(i+intervalo(1)-1)-samp(1)+1-round(maxQRSon/1000*heasig.freq*spf)):qrspiconset2,scale)];
end
if size(pontosnew,1)>1
weight=ones(length(pontosnew),1);
if i==1
if size(timenew,2)>1 %4JUL08
intervalonewlead=1:timenew(i+1);
else
intervalonewlead=1:length(w1);
end
elseif i<length(timenew)
intervalonewlead=timenew(i-1):timenew(i+1);
else
intervalonewlead=timenew(i-1):length(w1);
end
[a2,v2] = optimline(pontosnew, weight,OPT); %#ok<ASGLU>
else
v2=NaN;
end
if isnan(v2)
piconall(end)=NaN;
recursioncount.QRSon((i+intervalo(1)-1))=recursioncount.QRSon((i+intervalo(1)-1))-1;
else
if size(w1,2)==3
newleadbeat12=(w1(intervalonewlead,1).*v2(1)+w2(intervalonewlead,1)*v2(2)+w3(intervalonewlead,1)*v2(3))./norm(v2);
newleadbeat22=(w1(intervalonewlead,2).*v2(1)+w2(intervalonewlead,2)*v2(2)+w3(intervalonewlead,2)*v2(3))./norm(v2);
%signew2=(sig(intervalonewlead,1).*v2(1)+sig(intervalonewlead,2)*v2(2)+sig(intervalonewlead,3)*v2(3))./norm(v2);
else
newleadbeat12=(w1(intervalonewlead,1).*v2(1)+w2(intervalonewlead,1)*v2(2))./norm(v2);
newleadbeat22=(w1(intervalonewlead,2).*v2(1)+w2(intervalonewlead,2)*v2(2))./norm(v2);
%signew2=(sig(intervalonewlead,1).*v2(1)+sig(intervalonewlead,2)*v2(2))./norm(v2);
end
waux(intervalonewlead,1)=newleadbeat12;
waux(intervalonewlead,2)=newleadbeat22;
positionaux=position;
positionaux.R(i+intervalo(1)-1)=position_safe.R(i+intervalo(1)-1); %13.AGO.07
[positionaux,qrspiconset2,qrspicoffset2]= qrswavef(heasig,samp,time0(i),positionaux,waux,[i+intervalo(1)-1 i+intervalo(1)-1],messages); %#ok<NASGU>
if isempty(qrspiconset2)||isnan(qrspiconset2) % RUTE 3ABRI09
piconall=[piconall NaN]; %#ok<AGROW>
amppiconall=[amppiconall NaN]; %#ok<AGROW>
amppiconaux= NaN;
else
piconall=[piconall qrspiconset2]; %#ok<AGROW>
amppiconall=[amppiconall (newleadbeat22(qrspiconset2-intervalonewlead(1)+1))]; %#ok<AGROW>
amppiconaux= newleadbeat22aux(qrspiconset2-intervalonewlead(1)+1);
end
newleadbeat22aux=newleadbeat22;
newQRSon=[newQRSon position.QRSon(i+intervalo(1)-1)]; %#ok<AGROW>
% stop criteria iii) | iv) => iv) ou novo
%if (positionaux.QRSon(i+intervalo(1)-1)>=position.R(i+intervalo(1)-1)& position.QRSon(i+intervalo(1)-1)<position.R(i+intervalo(1)-1))|((abs(amppiconall(end))<abs(amppiconall(end-1)))&((amppiconaux(end)*amppiconall(end))>0) & ~isnan(position.Q(i+intervalo(1)-1)))|( isnan(piconall(end)))
%13.AGO.07
% stop criteria iii) | iv) => iv) %opçao 0 % parte de alternativa1e2
if ((abs(amppiconall(end))<abs(amppiconall(end-1)))&&((amppiconaux(end)*amppiconall(end))>0) &&...
~isnan(position.Q(i+intervalo(1)-1)))||( isnan(piconall(end)))
%alternative criteria stop criteria iii) without Qcondition | iv) => iv) alternativa3
%if ((abs(amppiconall(end))<abs(amppiconall(end-1)))&((amppiconaux(end)*amppiconall(end))>0) )|( isnan(piconall(end)))
%alternative criteria stop criteria just iv) => iv)
%if ( isnan(piconall(end)))
if ( isnan(piconall(end)))
position.QRSonsetcriteria((i+intervalo(1)-1))=4;
else
position.QRSonsetcriteria((i+intervalo(1)-1))=3;
end
piconall(end)=NaN;
recursioncount.QRSon((i+intervalo(1)-1))=recursioncount.QRSon((i+intervalo(1)-1))-1;
else
% if (abs(newQRSon(end)-newQRSon(end-1))<2 )
% position.QRSonsetcriteria((i+intervalo(1)-1))=1;
%elseif nansum(newQRSon(end)-newQRSon(1:end-1)==0)>1
if nansum(newQRSon(end)-newQRSon(1:end-1)==0)>1 % so em alternativa2
position.QRSonsetcriteria((i+intervalo(1)-1))=2;
end
position.QRSon(i+intervalo(1)-1)=positionaux.QRSon(i+intervalo(1)-1) ;
position.Q(i+intervalo(1)-1) = positionaux.Q(i+intervalo(1)-1);
position_safe.R(i+intervalo(1)-1) = positionaux.R(i+intervalo(1)-1) ;
position.R_inQRSon(i+intervalo(1)-1) = positionaux.R(i+intervalo(1)-1) ;
end
end
if (abs(newQRSon(end)-newQRSon(end-1))>1 )
position.QRSoncriteria(i+intervalo(1)-1)=1;
end
end
end
end
%%% protection: if 2D is being used and QRSon was marked after R
%%% %05Sep07
if size(w1,2)==2 && position.QRSon(i+intervalo(1)-1)>=(-1+position.qrs(i+intervalo(1)-1)) %invalid nark
position.QRSon(i+intervalo(1)-1)=NaN;
position.QRSoncriteria(i+intervalo(1)-1)=-1;
end
position.R(i+intervalo(1)-1)=position_safe.R(i+intervalo(1)-1); %13.AGO.07
end