function [position0,position,messages]= delineate3D(w1,w2,w3,intreg,timenew,position0,position,intervalo,position1,position2,position3,indexes,samp,sig,messages)
%
%
% [position0,position,A,V,A2,V2]= delineate3D(w1,w2,w3,intreg,scale,timenew,position0,position,intervalo,position1,samp)
% 3D multilead delineation of T wave
%
% Input Parameters:
%   w1: matrix with WT scales 1 to 5 for lead 1
%   w2: matrix with WT scales 1 to 5 for lead 2
%   w3: matrix with WT scales 1 to 5 for lead 1
%  intreg: interval to be considered in fitting the optimal direction
%  scale to use in delineation
%   timenew:  QRS times in inedexes refering the interval included in the current excerpt (borders excluded)
%   position0: struct vector with the detected points for 3D multilead  step1
%   position: struct vector with the detected points for 3D multilead  step2
%   position1: struct vector with the detected points for lead1 - to be replaced for multilead marks
%   samp: samples included in the current excerpt (borders excluded)
%   intervalo: numeration of the beats processed in this segment
%
%Output Parameters:
%   actualized parameters: position, position0
% A - points of fitted lines in step 1
% A2 - points of fitted lines in step 2
%  V- vectors director to the fitted line
%  V2- vectors director to the fitted line
%
% Rute Almeida: trabalho de doutoramento 2.Dec.2004
% Last update: 05AGO2011
%
% MATLAB Version 6.5.0.180913a (R13)

global recursioncount OPT
A=[];A2=[];
V=[];V2=[];
if isempty(sig)
    
end
% Aon=[];Von=[];
% Vg=[];
% if nargin==15
% messages = new_empty_mesg;
% end
if ~isfield(messages.setup.wavedet,'rrant_mintol')
    messages.setup.wavedet.rrant_mintol  =  0.5;    % minimum time in sec admited for current RR to be used in T wave delineation or in Exponentially averaged RR
end
if ~isfield(messages.setup.wavedet,'rrant_maxtol')
    messages.setup.wavedet.rrant_maxtol  =  1.5;    % max time in sec admited for current RR to be used in Exponentially averaged RR
end
if ~isfield(messages.setup.wavedet,'rrant_average')
    messages.setup.wavedet.rrant_average  =  [0.8 0.2];    % Exponentially averaged RR weights
end
if ~isfield(messages.setup.wavedet,'scale')
    messages.setup.wavedet.scale=4;
end
if ~isfield(messages.setup.wavedet,'scale2')
    messages.setup.wavedet.scale2=5;
end
if ~isfield(messages.setup.wavedet,'T_CSE_tol')
    messages.setup.wavedet.T_CSE_tol=30.6*2/1000; % based on CSE tolerance
end
if ~isfield(messages.setup.wavedet,'Tconvergence_crit')
    messages.setup.wavedet.Tconvergence_crit=1; % 1 maks are acepted as the same if they difer by Tconvergence_crit
end

rrant_mintol=messages.setup.wavedet.rrant_mintol;
rrant_maxtol=messages.setup.wavedet.rrant_maxtol;
rrant_average=messages.setup.wavedet.rrant_average;
scale=messages.setup.wavedet.scale;
scale2=messages.setup.wavedet.scale2;
Tconvergence_crit=messages.setup.wavedet.Tconvergence_crit;

intervalo=intervalo(1):intervalo(end);

scalei=scale*ones(1,size(intreg,1));
if ~isempty(w3)
    %scalei(sum([position1.Tscale(intervalo);position2.Tscale(intervalo);position3.Tscale(intervalo)]==5)>1)=5;
    %Rute 17Mar06
    aux=[(indexes(1,intervalo));(indexes(2,intervalo));(indexes(3,intervalo))];
    aux(1,~isnan(aux(1,:)))=position1.Tscale(aux(1,~isnan(aux(1,:))));
    aux(2,~isnan(aux(2,:)))=position2.Tscale(aux(2,~isnan(aux(2,:))));
    aux(3,~isnan(aux(3,:)))=position3.Tscale(aux(3,~isnan(aux(3,:))));
    scalei(sum(aux==5)>1)=5;
    %%%%%%
else
    scalei(sum([position1.Tscale(intervalo);position2.Tscale(intervalo)]==scale2)>1)=scale2;
end

for i=1:size(intreg,1)

    if isnan(intreg(i,1))
        position.Tscale(intervalo(i))=NaN;
        position.Toff(intervalo(i))=NaN;
        position.T(intervalo(i))=NaN;
        position.Ttipooff(intervalo(i))=NaN;
        position.Ttipo(intervalo(i))=NaN;
        position.Tprima(intervalo(i))=NaN;
        position.Ton(intervalo(i))=NaN;
        position.Ttipoon(intervalo(i))=NaN;
    else

        scale=scalei(i);
        %marks taken as the last in the 3 leads

        aux2=find(~isnan(indexes(2,intervalo(i))));
        aux3=find(~isnan(indexes(3,intervalo(i))));
        S=max([position1.S(indexes(~isnan(indexes(1,intervalo(i))),intervalo(i))) position2.S(indexes(2*aux2,intervalo(i))) position3.S(indexes(3*aux3,intervalo(i))) ]-samp(1)+1);
        QRSoff=max([position1.QRSoff(indexes(~isnan(indexes(1,intervalo(i))),intervalo(i))) position2.QRSoff(indexes(2*aux2,intervalo(i))) position3.QRSoff(indexes(3*aux3,intervalo(i)))]-samp(1)+1);

        if ~isempty(w3)
            pontos=[w1(intreg(i,1):min(intreg(i,2),length(w1)),scale) w2(intreg(i,1):min(intreg(i,2),length(w2)),scale) w3(intreg(i,1):min(intreg(i,2),length(w3)),scale) ];
        else
            pontos=[w1(intreg(i,1):min(intreg(i,2),length(w1)),scale) w2(intreg(i,1):min(intreg(i,2),length(w2)),scale)];
        end
        weight=ones(size(pontos,1),1);
        [a,v] = optimline(pontos, weight,OPT);

        A=[A;a]; %#ok<AGROW>
        V=[V;v]; %#ok<AGROW>
    

        %4FEB2010
        if length(intervalo)>1
            rrant=position.qrs(intervalo(2))-position.qrs(intervalo(1));
        elseif intervalo(1)>1
            rrant=position.qrs(intervalo(1))-position.qrs(intervalo(1)-1);
        else
            rrant=messages.setup.wavedet.freq/2; %JUL2011
        end

        if (i>1),
            if (rrant_maxtol*rrant>(timenew(i)-timenew(i-1)))&&(timenew(i)-timenew(i-1)>rrant_mintol*rrant),
                rrmed = rrant_average(1)*rrant + rrant_average(2)*(timenew(i)-timenew(i-1));
            else
                rrmed = rrant;  % Exponentially averaged RR
            end
        else                  % Only for the first in each segment of the ECG
            rrmed = rrant;
        end
        rrant = rrmed;        %#ok<NASGU> % For next segment

        %WT projected: from the previous QRS complex to the following QRS complex
        if ~isempty(w3)
            if i<length(timenew),
                newleadbeat=(w1(timenew(i):timenew(i+1),scale).*v(1)+w2(timenew(i):timenew(i+1),scale)*v(2)+w3(timenew(i):timenew(i+1),scale)*v(3))./norm(v);
                %signew=(sig(timenew(i):timenew(i+1),1).*v(1)+sig(timenew(i):timenew(i+1),2)*v(2)+sig(timenew(i):timenew(i+1),3)*v(3))./norm(v);
            else                        % if last beat of the segment
                newleadbeat=(w1(timenew(i):end,scale).*v(1)+w2(timenew(i):end,scale)*v(2)+w3(timenew(i):end,scale)*v(3))./norm(v);
                %signew=(sig(timenew(i):end,1).*v(1)+sig(timenew(i):end,2)*v(2)+sig(timenew(i):end,3)*v(3))./norm(v);
            end
        else
            if i<length(timenew),
                newleadbeat=(w1(timenew(i):timenew(i+1),scale).*v(1)+w2(timenew(i):timenew(i+1),scale)*v(2))./norm(v);
                %signew=(sig(timenew(i):timenew(i+1),1).*v(1)+sig(timenew(i):timenew(i+1),2)*v(2))./norm(v);
            else                        % if last beat of the segment
                newleadbeat=(w1(timenew(i):end,scale).*v(1)+w2(timenew(i):end,scale)*v(2))./norm(v);
                %signew=(sig(timenew(i):end,1).*v(1)+sig(timenew(i):end,2)*v(2))./norm(v);
            end
        end
        [pos,picon,picoff,janelas1,messages]=twave3D( newleadbeat,timenew,i,messages.setup.wavedet.freq,S,QRSoff,samp,rrmed,messages); %#ok<ASGLU>
%         janelas1=janelas1+timenew(i)-1;
        position0.Tscale(intervalo(i))=scale;
        position0.Ton(intervalo(i))=pos.Ton;
        position0.Toff(intervalo(i))=pos.Toff;
        position0.T(intervalo(i))=pos.T;
        position0.Ttipo(intervalo(i))=pos.Ttipo;
        position0.Tprima(intervalo(i))=pos.Tprima;
        
        %position0.Wavedet(intervalo(i))=newleadbeat; % derivada wavedet Maikel 30 marzo 2009
        %position0.direccion(c)=V; % vector direcci�n Maikel 30 marzo 2009
        %pontos2=[w1(picoff:intreg(i,2),scale) w2(picoff:intreg(i,2),scale) w3(picoff:intreg(i,2),scale) ];

        if ~isnan(picoff)
            if picoff>(intreg(i,2)-1)
                messages.warnings=[messages.warnings {'Tend out of the searching interval, in the fist step of multilead!'}];
                position.Tscale(intervalo(i))=NaN;
                position.Toff(intervalo(i))=NaN;
                position.T(intervalo(i))=NaN;
                position.Ttipooff(intervalo(i))=NaN;
                position.Ttipo(intervalo(i))=NaN;
                position.Tprima(intervalo(i))=NaN;
                recursioncount.Toff(intervalo(i))=NaN;
                position.contadorToff(intervalo(i))=NaN; %Rute
            else
                if (pos.Toff-samp(1)+1)==picoff
                    endaux=min([(position0.Toff(intervalo(i))-samp(1)+1+round(messages.setup.wavedet.T_CSE_tol*messages.setup.wavedet.freq)) length(w1)]); %3JUL

                else
                    endaux=min([intreg(i,2) (position0.Toff(intervalo(i))-samp(1)+1+messages.setup.wavedet.T_CSE_tol*messages.setup.wavedet.freq) ]);
                end
                if ~isempty(w3)
                    pontos2=[w1(picoff:endaux,scale) w2(picoff:endaux,scale) w3(picoff:endaux,scale)];
                    %pontos2=[w1(janelas1(2):janelas1(3),scale) w2(janelas1(2):janelas1(3),scale) w3(janelas1(2):janelas1(3),scale)];
                else
                    pontos2=[w1(picoff:endaux,scale) w2(picoff:endaux,scale)];
                end
                weight2=ones(size(pontos2,1),1);
                [a,v] = optimline(pontos2, weight2,OPT);
                A2=[A2;a]; %#ok<AGROW>
                V2=[V2;v]; %#ok<AGROW>
                if ~isempty(w3)
                    if i<length(timenew),
                        newleadbeat2=(w1(timenew(i):timenew(i+1),scale).*v(1)+w2(timenew(i):timenew(i+1),scale)*v(2)+w3(timenew(i):timenew(i+1),scale)*v(3))./norm(v);
                        %signew2=(sig(timenew(i):timenew(i+1),1).*v(1)+sig(timenew(i):timenew(i+1),2)*v(2)+sig(timenew(i):timenew(i+1),3)*v(3))./norm(v);
                    else                        % if last beat of the segment
                        newleadbeat2=(w1(timenew(i):end,scale).*v(1)+w2(timenew(i):end,scale)*v(2)+w3(timenew(i):end,scale)*v(3))./norm(v);
                        %signew2=(sig(timenew(i):end,1).*v(1)+sig(timenew(i):end,2)*v(2)+sig(timenew(i):end,3)*v(3))./norm(v);
                    end
                else
                    if i<length(timenew),
                        newleadbeat2=(w1(timenew(i):timenew(i+1),scale).*v(1)+w2(timenew(i):timenew(i+1),scale)*v(2))./norm(v);
                        %signew2=(sig(timenew(i):timenew(i+1),1).*v(1)+sig(timenew(i):timenew(i+1),2)*v(2))./norm(v);
                    else                        % if last beat of the segment
                        newleadbeat2=(w1(timenew(i):end,scale).*v(1)+w2(timenew(i):end,scale)*v(2))./norm(v);
                        %signew2=(sig(timenew(i):end,1).*v(1)+sig(timenew(i):end,2)*v(2))./norm(v);
                    end
                end
                [pos,picon2,picoff2,janelas2,messages]=twave3D(newleadbeat2,timenew,i,messages.setup.wavedet.freq,S,QRSoff,samp,rrmed,messages); %#ok<ASGLU>
%                 janelasnew=janelas2+timenew(i)-1;
                picoffall=[picoff picoff2];
                if isnan(picoff2)
                    amppicoffall=[abs(newleadbeat(picoff-timenew(i)+1)) NaN];
                else
                    amppicoffall=[abs(newleadbeat(picoff-timenew(i)+1)) abs(newleadbeat2(picoff2-timenew(i)+1))];
                end
                newToff=[position0.Toff(intervalo(i)) pos.Toff];

                if (amppicoffall(end)<amppicoffall(end-1))
                    position.contadorToff(intervalo(i))=-1;
                    %Vg(i,:)=V(end,:);
                    position.Tscale(intervalo(i))=position0.Tscale(intervalo(i));
                    position.Toff(intervalo(i))=position0.Toff(intervalo(i));
                    position.T(intervalo(i))= position0.T(intervalo(i));
                    position.Ttipooff(intervalo(i))=position0.Ttipo(intervalo(i));
                    position.Ttipo(intervalo(i))=NaN;
                    position.Tprima(intervalo(i))=position0.Tprima(intervalo(i));
                    %position.contadorToff=recursioncount.Toff;
                    position.direccionToffopt(:,intervalo(i))=V(end,:);
                    %position.Wavedet(intervalo(i))=newleadbeat2; % derivada wavedet Maikel 30 marzo 2009
                    %position.direccion(:,intervalo(i))=V2; % vector direcci�n Maikel 30 marzo 2009
                else
                    position.contadorToff(intervalo(i))=0;
                    %Vg(i,:)=V2(end,:);
                    position.direccionToffopt(:,intervalo(i))=V2(end,:);
                    %%%%%%%%%%%%%%%%%%%%%%recursion
                    while (~isnan(picoffall(end)) &&  abs(newToff(end)-newToff(end-1))>Tconvergence_crit) && (amppicoffall(end)>amppicoffall(end-1))
                        position.contadorToff(intervalo(i))=position.contadorToff(intervalo(i))+1;
                        %maxTintervalend=maxTintervalend; % teste 14
                        if (pos.Toff-samp(1)+1)==picoff
                            endaux=(pos.Toff-samp(1)+1+round(messages.setup.wavedet.T_CSE_tol*messages.setup.wavedet.freq));
                        else
                            endaux=min([intreg(i,2) (pos.Toff-samp(1)+1+round(messages.setup.wavedet.T_CSE_tol*messages.setup.wavedet.freq))]);
                        end
                        if ~isempty(w3)
                            pontosnew=[w1(picoffall(end):endaux,scale) w2(picoffall(end):endaux,scale) w3(picoffall(end):endaux,scale)];
                            %pontosnew=[w1(picoff(end):endaux,scale) w2(picoff(end):endaux,scale) w3(picoff(end):endaux,scale)];
                            %pontosnew=[w1(janelasnew(2):janelasnew(3),scale) w2(janelasnew(2):janelasnew(3),scale) w3(janelasnew(2):janelasnew(3),scale)];
                        else
                            pontosnew=[w1(picoffall(end):endaux,scale) w2(picoffall(end):endaux,scale)];
                        end

                        if   size(pontosnew,1)>1
                            weightnew=ones(size(pontosnew,1),1);
                            [a,v] = optimline(pontosnew, weightnew,OPT);
%                             if i==6
%                                 title([ecgnr '(Tend-step' num2str(2+position.contadorToff(intervalo(i)))')'])
%                             else
%                                 close
%                             end
                            A2=[A2;a]; %#ok<AGROW>
                            V2=[V2;v]; %#ok<AGROW>
                            if ~isempty(w3)
                                if i<length(timenew),
                                    newleadbeatnew=(w1(timenew(i):timenew(i+1),scale).*v(1)+w2(timenew(i):timenew(i+1),scale)*v(2)+w3(timenew(i):timenew(i+1),scale)*v(3))./norm(v);
                                    %signewnew=(sig(timenew(i):timenew(i+1),1).*v(1)+sig(timenew(i):timenew(i+1),2)*v(2)+sig(timenew(i):timenew(i+1),3)*v(3))./norm(v);
                                else                        % if last beat of the segment
                                    newleadbeatnew=(w1(timenew(i):end,scale).*v(1)+w2(timenew(i):end,scale)*v(2)+w3(timenew(i):end,scale)*v(3))./norm(v);
                                    %signewnew=(sig(timenew(i):end,1).*v(1)+sig(timenew(i):end,2)*v(2)+sig(timenew(i):end,3)*v(3))./norm(v);
                                end
                            else
                                if i<length(timenew),
                                    newleadbeatnew=(w1(timenew(i):timenew(i+1),scale).*v(1)+w2(timenew(i):timenew(i+1),scale)*v(2))./norm(v);
                                    %signewnew=(sig(timenew(i):timenew(i+1),1).*v(1)+sig(timenew(i):timenew(i+1),2)*v(2))./norm(v);
                                else                        % if last beat of the segment
                                    newleadbeatnew=(w1(timenew(i):end,scale).*v(1)+w2(timenew(i):end,scale)*v(2))./norm(v);
                                    %signewnew=(sig(timenew(i):end,1).*v(1)+sig(timenew(i):end,2)*v(2))./norm(v);
                                end
                            end

                            [pos,picon2,picoffnew,janelasnew,messages]=twave3D(newleadbeatnew,timenew,i,messages.setup.wavedet.freq,S,QRSoff,samp,rrmed,messages); %#ok<ASGLU>
%                             janelasnew=janelas2+timenew(i)-1;
                        else
                            picoffnew=NaN;
                        end
                        picoffall=[picoffall picoffnew]; %#ok<AGROW>
                        if isnan(picoffnew)
                            amppicoffall=[amppicoffall NaN]; %#ok<AGROW>
                        else
                            amppicoffall=[amppicoffall abs(newleadbeatnew(picoffnew-timenew(i)+1))]; %#ok<AGROW>
                        end
                        newToff=[newToff pos.Toff]; %#ok<AGROW>

                        if (amppicoffall(end)<amppicoffall(end-1))
                            pos.Toff=NaN; % because the lead got worse
                        end
                        if isnan(pos.Toff); %01Jun05
                            pos.Toff=newToff(end-1);
                            position.contadorToff(intervalo(i))=position.contadorToff(intervalo(i))-1;
                        else
                            %Vg(i,:)=V2(end,:);
                            position.direccionToffopt(:,intervalo(i))=V2(end,:);
                            %contador=[contador; recursioncount.Toff];
                        end
                        if  prod(picoffall(1:end-1)-picoffall(end))==0  % ciclo infinito no teste 14
                            picoffall(end)=NaN;
                        end
                    end
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    position.Tscale(intervalo(i))=scale;
                    position.Toff(intervalo(i))=pos.Toff;
                    position.T(intervalo(i))=pos.T;
                    position.Ttipooff(intervalo(i))=pos.Ttipo;
                    position.Ttipo(intervalo(i))=NaN;
                    position.Tprima(intervalo(i))=pos.Tprima;
                    %position.contadorToff=recursioncount.Toff;
                    %position.Wavedet(intervalo(i))=newleadbeatnew; % derivada wavedet Maikel 30 marzo 2009
                    %position.direccion(:,intervalo(i))=V2; % vector direcci�n Maikel 30 marzo 2009
                
                end
                %
            end
        else
            position.Tscale(intervalo(i))=NaN;
            position.Toff(intervalo(i))=NaN;
            position.T(intervalo(i))=NaN;
            position.Ttipooff(intervalo(i))=NaN;
            position.Ttipo(intervalo(i))=NaN;
            position.Tprima(intervalo(i))=NaN;
            recursioncount.Toff(intervalo(i))=NaN;
            position.contadorToff(intervalo(i))=NaN; %Rute
            %position.Wavedet(intervalo(i))=NaN;  % derivada wavedet Maikel 30 marzo 2009
            %position.direccion(:,intervalo(i))=NaN; % vector direcci�n Maikel 30 marzo 2009   
        end
        
        Tbeginauxiliar
    end
end
% position0.A(intervalo(i))=A;
% position0.A2(intervalo(i))=A2;
% position0.V(intervalo(i))=V;
% position0.V2(intervalo(i))=V2;