function [position0,position,messages]= multiqrsend(heasig,w1,w2,w3,auxqrspicoffset,QRSoff,time,timenew,position0,position,intervalo,samp,messages)
%[position0,position]= multiqrsend(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
% auxqrspicoffset: latest relevant modulus maximum location associated to QRS complex found in any lead
% QRSoff: latest QRS end location in any lead (initial QRS end 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
% Last update: 07FEB2012
global recursioncount OPT
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,'maxQRSoff')
messages.setup.wavedet.maxQRSoff =11.6*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;
maxQRSoff=messages.setup.wavedet.maxQRSoff;
auxqrsoff=nanmax(QRSoff-samp(1)+1)+round(maxQRSoff/1000*messages.setup.wavedet.freq);
for i=1:size(time,2)
position.QRSoffcriteria((i+intervalo(1)-1))=2;% default % 1 not considered
scale=2;
if size(w1,2)==3
pontos=[w1(auxqrspicoffset(i):min(auxqrsoff(i),end),scale)...
w2(auxqrspicoffset(i):min(auxqrsoff(i),end),scale)...
w3(auxqrspicoffset(i):min(auxqrsoff(i),end),scale)]; %4JUL2011 % protect with end
else
pontos=[w1(auxqrspicoffset(i):min(auxqrsoff(i),end),scale)...
w2(auxqrspicoffset(i):min(auxqrsoff(i),end),scale)];%4JUL2011 % protect with end
end
weight=ones(size(pontos,1),1);
[a,v] = optimline(pontos, weight,OPT); %#ok<ASGLU>
time0=round(nanmedian(time));% Rute 01.03.2010
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<ASGLU>
if isempty(qrspicoffset0)||isnan(qrspicoffset0) % RUTE 6ABRI09
% picoffall= NaN; %#ok<NASGU>
% amppicoffall= NaN; %#ok<NASGU>
% amppicoffaux= NaN; %#ok<NASGU>
recursioncount.QRSoff((i+intervalo(1)-1))=-2;
position.QRSoff(i+intervalo(1)-1) = NaN;
position.S(i+intervalo(1)-1) = NaN;
position_safe.R(i+intervalo(1)-1) = NaN; %13.AGO.07
position.RR_inQRSoff(i+intervalo(1)-1) = NaN;
position.Rprima(i+intervalo(1)-1)=NaN;
else
if size(w1,2)==3
pontos2=[w1(qrspicoffset0:min(end,(position0.QRSoff(i+intervalo(1)-1)-samp(1)+1+round(maxQRSoff/1000*messages.setup.wavedet.freq))),scale)...
w2(qrspicoffset0:min(end,(position0.QRSoff(i+intervalo(1)-1)-samp(1)+1+round(maxQRSoff/1000*messages.setup.wavedet.freq))),scale)...
w3(qrspicoffset0:min(end,(position0.QRSoff(i+intervalo(1)-1)-samp(1)+1+round(maxQRSoff/1000*messages.setup.wavedet.freq))),scale)];
else
pontos2=[w1(qrspicoffset0:min(end,(position0.QRSoff(i+intervalo(1)-1)-samp(1)+1+round(maxQRSoff/1000*messages.setup.wavedet.freq))),scale)...
w2(qrspicoffset0:min(end,(position0.QRSoff(i+intervalo(1)-1)-samp(1)+1+round(maxQRSoff/1000*messages.setup.wavedet.freq))),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=nan*ones(size(w1)); %Rute6NOV08
waux(intervalonewlead,1)=newleadbeat12;
waux(intervalonewlead,2)=newleadbeat22;
[positionaux,qrspiconset2,qrspicoffset2]= qrswavef(heasig,samp,time0(i),position,waux,[i+intervalo(1)-1 i+intervalo(1)-1],messages); %#ok<ASGLU>
position.QRSoff(i+intervalo(1)-1) = positionaux.QRSoff(i+intervalo(1)-1);
position.S(i+intervalo(1)-1) = positionaux.S(i+intervalo(1)-1);
position.R_inQRSoff(i+intervalo(1)-1) = positionaux.R(i+intervalo(1)-1) ;
position_safe.R(i+intervalo(1)-1) = positionaux.R(i+intervalo(1)-1) ;%13.AGO.07
position.Rprima(i+intervalo(1)-1) = positionaux.Rprima(i+intervalo(1)-1);
if isempty(qrspicoffset2)||isnan(qrspicoffset2) %16MAY2011
picoffall=[qrspicoffset0 NaN];
amppicoffall=[(newleadbeat2(qrspicoffset0-intervalonewlead(1)+1)) NaN];
amppicoffaux= NaN;
else
picoffall=[qrspicoffset0 qrspicoffset2];
amppicoffall=[(newleadbeat2(qrspicoffset0-intervalonewlead(1)+1)) (newleadbeat22(qrspicoffset2-intervalonewlead(1)+1))];
amppicoffaux= newleadbeat2(qrspicoffset2-intervalonewlead(1)+1);
end
newQRSoff=[NaN position0.QRSoff(i+intervalo(1)-1) position.QRSoff(i+intervalo(1)-1)];
newleadbeat22aux=newleadbeat22;
% stop criteria iii) | iv) ou novo
% if (positionaux.QRSoff(i+intervalo(1)-1)<=position.R(i+intervalo(1)-1)&...
%position0.QRSoff(i+intervalo(1)-1)>position.R(i+intervalo(1)-1))| ...
%((abs(amppicoffall(end))<abs(amppicoffall(end-1)))&((amppicoffaux(end)*amppicoffall(end))>0) &...
%~isnan(position0.S(i+intervalo(1)-1)))|( isnan(picoffall(end))) %| (abs(newQRSoff(end)-newQRSoff(end-1))<=1) % RUTE 9.07.2007
%13.AGO.07
% stop criteria iii) | iv) % opçao0 % parte da alternativa2
if ((abs(amppicoffall(end))<abs(amppicoffall(end-1)))&&((amppicoffaux(end)*amppicoffall(end))>0) &&...
~isnan(position0.S(i+intervalo(1)-1)))||( isnan(picoffall(end))) %| (abs(newQRSoff(end)-newQRSoff(end-1))<=1) % RUTE 9.07.2007
% stop criteria iii) | iv) | i) % alternativa1
%if ((abs(amppicoffall(end))<abs(amppicoffall(end-1)))&((amppicoffaux(end)*amppicoffall(end))>0) &...
%~isnan(position0.S(i+intervalo(1)-1)))|( isnan(picoffall(end))) | (abs(newQRSoff(end)-newQRSoff(end-1))<=1) % RUTE 9.07.2007
%alternative criteria stop criteria iii)without the S wave condition| iv) alternativa3
%if ((abs(amppicoffall(end))<abs(amppicoffall(end-1)))&((amppicoffaux(end)*amppicoffall(end))>0) )|( isnan(picoffall(end))) %| (abs(newQRSoff(end)-newQRSoff(end-1))<=1) % RUTE 9.07.2007
%if ( isnan(picoffall(end))) | (abs(newQRSoff(end)-newQRSoff(end-1))<=1) % RUTE 9.07.2007
if isnan(picoffall(end))
position.QRSoffcriteria(i+intervalo(1)-1)=4;
% alternativa2 - exclude criteria i)
elseif (abs(newQRSoff(end)-newQRSoff(end-1))<=1)
position.QRSoffcriteria(i+intervalo(1)-1)=1;
else
position.QRSoffcriteria(i+intervalo(1)-1)=3;
end
recursioncount.QRSoff((i+intervalo(1)-1))=-1;
position.QRSoff(i+intervalo(1)-1) = position0.QRSoff(i+intervalo(1)-1);
position.S(i+intervalo(1)-1) = position0.S(i+intervalo(1)-1);
position_safe.R(i+intervalo(1)-1) = position0.R(i+intervalo(1)-1) ; %13.AGO.07
position.Rprima(i+intervalo(1)-1) = position0.Rprima(i+intervalo(1)-1);
else
recursioncount.QRSoff((i+intervalo(1)-1))=0;
% stop criteria iv) | ii) | i) % opçao0 % parte da alternativa1
%while ~isnan(picoffall(end)) & nansum(newQRSoff(end)-newQRSoff(1:end-1)==0)<2 & (abs(newQRSoff(end)-newQRSoff(end-1))>1)%%%
%alternative criteria stop criteria iv) | ii)) % alternativa2e3
while ~isnan(picoffall(end)) && nansum(newQRSoff(end)-newQRSoff(1:end-1)==0)<2 %& (abs(newQRSoff(end)-newQRSoff(end-1))>1)%%%
recursioncount.QRSoff((i+intervalo(1)-1))=recursioncount.QRSoff((i+intervalo(1)-1))+1;
if size(w1,2)==3
pontosnew=[w1(qrspicoffset2:min(end,(position.QRSoff(i+intervalo(1)-1)-samp(1)+1+round(maxQRSoff/1000*messages.setup.wavedet.freq))),scale)...
w2(qrspicoffset2:min(end,(position.QRSoff(i+intervalo(1)-1)-samp(1)+1+round(maxQRSoff/1000*messages.setup.wavedet.freq))),scale)...
w3(qrspicoffset2:min(end,(position.QRSoff(i+intervalo(1)-1)-samp(1)+1+round(maxQRSoff/1000*messages.setup.wavedet.freq))),scale)];
else
pontosnew=[w1(qrspicoffset2:min(end,(position.QRSoff(i+intervalo(1)-1)-samp(1)+1+round(maxQRSoff/1000*messages.setup.wavedet.freq))),scale)...
w2(qrspicoffset2:min(end,(position.QRSoff(i+intervalo(1)-1)-samp(1)+1+round(maxQRSoff/1000*messages.setup.wavedet.freq))),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)
picoffall(end)=NaN;
recursioncount.QRSoff((i+intervalo(1)-1))=recursioncount.QRSoff((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<ASGLU>
picoffall=[picoffall qrspicoffset2]; %#ok<AGROW>
amppicoffall=[amppicoffall (newleadbeat22(qrspicoffset2-intervalonewlead(1)+1))]; %#ok<AGROW>
amppicoffaux= newleadbeat22aux(qrspicoffset2-intervalonewlead(1)+1);
newleadbeat22aux=newleadbeat22;
newQRSoff=[newQRSoff position.QRSoff(i+intervalo(1)-1)]; %#ok<AGROW>
% stop criteria iii) | iv) => iv) ou novo
%if (positionaux.QRSoff(i+intervalo(1)-1)<=position.R(i+intervalo(1)-1) &...
%position.QRSoff(i+intervalo(1)-1)>=position.R(i+intervalo(1)-1) ) |((abs(amppicoffall(end))<abs(amppicoffall(end-1)))&...
%((amppicoffaux(end)*amppicoffall(end))>0) & ~isnan(position.S(i+intervalo(1)-1)))|(isnan(picoffall(end)))
%13.AGO.07
% stop criteria iii) | iv) => iv) %opçao 0 parte da alternativa1e2
if ((abs(amppicoffall(end))<abs(amppicoffall(end-1)))&&((amppicoffaux(end)*amppicoffall(end))>0) &&...
~isnan(position.S(i+intervalo(1)-1)))||(isnan(picoffall(end)))
%alternative criteria stop criteria iii) without S condition | iv) => iv) alternativa3
%if ((abs(amppicoffall(end))<abs(amppicoffall(end-1)))&((amppicoffaux(end)*amppicoffall(end))>0) )|(isnan(picoffall(end)))
%alternative criteria stop criteria iv) => iv)
%if (isnan(picoffall(end)))
if isnan(picoffall(end))
position.QRSoffcriteria(i+intervalo(1)-1)=4;
else
position.QRSoffcriteria(i+intervalo(1)-1)=3;
end
picoffall(end)=NaN;
recursioncount.QRSoff((i+intervalo(1)-1))=recursioncount.QRSoff((i+intervalo(1)-1))-1;
else
%if ~(abs(newQRSoff(end)-newQRSoff(end-1))>1)
% position.QRSoffcriteria(i+intervalo(1)-1)=1;
%elseif nansum(newQRSoff(end)-newQRSoff(1:end-1)==0)>1
if nansum(newQRSoff(end)-newQRSoff(1:end-1)==0)>1 % alternativa2
position.QRSoffcriteria(i+intervalo(1)-1)=2;
end
position.QRSoff(i+intervalo(1)-1) = positionaux.QRSoff(i+intervalo(1)-1) ;
position.S(i+intervalo(1)-1) = positionaux.S(i+intervalo(1)-1);
%position.R_inQRSon(i+intervalo(1)-1) = positionaux.R(i+intervalo(1)-1) ;
position_safe.R(i+intervalo(1)-1) = positionaux.R(i+intervalo(1)-1) ;
position.Rprima(i+intervalo(1)-1) = positionaux.Rprima(i+intervalo(1)-1);
end
end
if (abs(newQRSoff(end)-newQRSoff(end-1))>1)
position.QRSoffcriteria(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.QRSoff(i+intervalo(1)-1)<=(1+position.qrs(i+intervalo(1)-1)-1)
position.QRSoff(i+intervalo(1)-1)=NaN;
position.QRSoffcriteria(i+intervalo(1)-1)=-1;
end
end