function [pos,picon,picoff,janelas,messages]=twave3D(w,timenew,i,freq,S,QRSoff,samp,rrmed,messages)
%
%[pos,picon,picoff,janelas]=twave3D(w,timenew,i,freq,S,QRSoff,samp,rrmed)
%
% multilead T wave delineation
%
%Input Parameters:
%   w: matrix with WT scales 1 to 5
%   timenew:  QRS times in inedexes refering the interval included in the current excerpt (borders excluded)
%   i: beat number
%   freq: sampling frequency (heasig.freq)
%   S: lattest S wave position found in any lead
%   QRSoff: lattest QRS end position found in any lead
%   samp: samples included in the current excerpt (borders excluded)
%   rrmed: Exponentially averaged RR
%
%Output Parameters:
%   pos: fiducial marks structure (position)
%   picon: position of the first relevant modulos maximum in the wavelet
%   picoff: position of the last relevant modulos maximum in the wavelet
%   janelas: T wave seach windows
%
% Rute Almeida
% Last update: Rute Almeida  05AGO2011
%
% Designed for MATLAB Version R12; tested with MATLAB Version R13
%

%%%%%% Constants and Thresholds  !!!!!!!!!!!!!!!!!!!!!!!!!
if ~isfield(messages.setup.wavedet,'umbraldetT')
    messages.setup.wavedet.umbraldetT = 0.25;  % We use umbraldet*sqrt(mean(w(time(i):time(i+1),4).^2))
end
if ~isfield(messages.setup.wavedet,'umbralsig')
    messages.setup.wavedet.umbralsig  =  1/8;
end
if ~isfield(messages.setup.wavedet,'Kton')
    messages.setup.wavedet.Kton = 4;      % 2 4!
end
if ~isfield(messages.setup.wavedet,'inivent_tol')
    messages.setup.wavedet.inivent_tol  =  0.1;
end
if ~isfield(messages.setup.wavedet,'inivent_tol_S')
    messages.setup.wavedet.inivent_tol_S=0.05; % sec
end
if ~isfield(messages.setup.wavedet,'finvent_tol')
    messages.setup.wavedet.finvent_tol=0.240;% sec
end
if ~isfield(messages.setup.wavedet,'finvent_max')
    messages.setup.wavedet.finvent_max=0.6;% sec
end
if ~isfield(messages.setup.wavedet,'min_vent')
    messages.setup.wavedet.min_vent=0.1;%sec
end
if ~isfield(messages.setup.wavedet,'Tmax_Tmin_time_min')
    messages.setup.wavedet.Tmax_Tmin_time_min=0.15;%sec
end
if ~isfield(messages.setup.wavedet,'Tmax_Tmin_bifasic')
    messages.setup.wavedet.Tmax_Tmin_bifasic=2.5;
end

Kton=messages.setup.wavedet.Kton;
Ktoff=messages.setup.wavedet.Kton;
messages.warnings=[messages.warnings {'Recall that in twave3D ktoff=kton.'}];
umbraldetT = messages.setup.wavedet.umbraldetT;
umbralsig = messages.setup.wavedet.umbralsig;
inivent_tol= messages.setup.wavedet.inivent_tol;
inivent_tol_S=messages.setup.wavedet.inivent_tol_S;
finvent_tol=messages.setup.wavedet.finvent_tol;
finvent_max= messages.setup.wavedet.finvent_max;
min_vent=messages.setup.wavedet.min_vent;
Tmax_Tmin_time_min= messages.setup.wavedet.Tmax_Tmin_time_min;%sec
Tmax_Tmin_bifasic=messages.setup.wavedet.Tmax_Tmin_bifasic;% extra criteria
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
janelas=[];%%31May05
% Initialization of auxiliary variables
T = []; Tprima=[]; picon=[]; picoff=[]; Ton=[]; Toff=[]; tipoT=[];
% minapos = []; minppos=[]; maxapos=[]; maxppos=[]; mina=[]; minp=[]; maxa=[];
% maxp=[];
picon_keep=[];%multileadchange
picoff_keep=[];%multileadchange
lead_keep=[];%multileadchange

inivent = round(inivent_tol*freq);   % Begining of window
if ~isempty(S),              % If there is an S wave
    inivent = max(inivent, S-timenew(i)+round(inivent_tol_S*freq));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% changed 13/06/02 Rute
if rrmed <= freq,
    %if rrmed >= freq,
    finvent = round(finvent_max*freq);     % End of window
else
    finvent = round(rrmed*finvent_max);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %% changed 13/06/02 Rute

% using timenew from diffferent leads ty ca dist less tahn 0.24 sec!!!
if i~=length(timenew),                 % For last beat in the segment
    finvent = min(finvent,timenew(i+1)-timenew(i)-round(finvent_tol*freq));
else
    finvent = min(finvent,timenew(i)-timenew(i-1)-round(finvent_tol*freq));
end

% We work at scale scale (in general).
if isempty(QRSoff),   % Should never happen, but...
    % begwin = inivent + timenew(i);
    begwin = inivent;
else
    %  begwin = max(inivent + timenew(i), QRSoff+1);
    begwin = max(inivent , QRSoff+1-timenew(i));
end
%endwin = min(finvent + timenew(i),length(w)); %% Rute 20/05/02
endwin = min(finvent ,length(w)); %% Rute 20/05/02
janelas=[janelas; i begwin endwin];%31May05
if ~isempty(begwin+1:endwin ) % 22.04.05
    
    % Positive maxima and negative minima in the window
    maxpos = begwin + modmax(w(begwin+1:endwin ),2,0,+1);
    minpos = begwin + modmax(w(begwin+1:endwin ),2,0,-1);
    [maxim ind] = max(w(maxpos ));   % The biggest of the positive
    maxpos = maxpos(ind);
    [minim ind] = min(w(minpos ));   % The biggest of the negative
    minpos = minpos(ind);
    
    if isempty(maxpos),               % If no local positive maximum
        % The maximum will be the first
        % or the last sample
        if (w(begwin )>=w(endwin )) && w(begwin )>0,
            maxpos = begwin; maxim = w(maxpos );
        elseif (w(endwin )>=w(begwin )) && w(endwin )>0,
            maxpos = endwin; maxim = w(maxpos );
        end
    end
    if isempty(minpos),               % if no local negative minimum
        % the minimum will be the first
        % or the last sample
        if (w(begwin )<=w(endwin )) && w(begwin )<0,
            minpos = begwin; minim = w(minpos );
        elseif (w(endwin )<=w(begwin )) && w(endwin )<0,
            minpos = endwin; minim = w(minpos );
        end
    end
    
    absmax = abs(maxim);
    absmin = abs(minim);
    
    if i<length(timenew),
        veficaz = sqrt(mean(w .^2)); % all w 03.Dec.04 interval was restricted before
    else                        % if last beat of the segment
        veficaz = sqrt(nanmean(w.^2)); % all w 03.Dec.04 interval was restricted before
    end
    
    hay_onda = ((absmax>umbraldetT*veficaz)|(absmin>umbraldetT*veficaz));
    
    % Rute 18/06/02
    if endwin-begwin<(min_vent*freq) % se a janela tem amplitude menor do que min_vent seg enato nao existe onda T
        hay_onda =0;
    end
    
    % Is there a wave?
    if hay_onda,
        if absmax >= absmin,    % the greatest modulus maximum is the maximum
            % Now we search the two minima nearest to maxpos, one before and one after
            minapos = max(modmax(w(begwin+1:maxpos-1 ),2,0,-1));
            minapos = begwin + minapos; % Position of the negative minimum before the maximum
            if isempty(minapos) && (maxpos ~= begwin) && (w(begwin )<0),
                minapos = begwin;% If no local minimum before the maximum, take the first sample
            end
            minppos = min(modmax(w(maxpos+1:endwin ),2,0,-1));
            minppos = maxpos + minppos; % Position of the positive maximum after the minimum
            if isempty(minppos) && (maxpos ~= endwin) && (w(endwin )<0),
                minppos = endwin;        % If no local minimum after the maximum, take the last sample
            end
            
            mina = abs(w(minapos ));     % Amplitude of minimum before maximum
            minp = abs(w(minppos ));     % Amplitude of minimum after maximum
            
            if (mina < umbralsig*absmax), % If mina is not big enough
                mina =[];                  % forget it
            elseif (maxpos-minapos>Tmax_Tmin_time_min*freq), %or if ther are more than 150 ms to maxpos
                mina = [];
            end
            if (minp < umbralsig*absmax), % If minp is not big enough
                minp =[];                  % forget it
            elseif (minppos-maxpos>Tmax_Tmin_time_min*freq), % and also if there are more than 150 ms to maxpos
                minp =[];
            end
            
            if ~isnan(mina)&~isnan(minp),      %#ok<AND2> %%% NUEVO JP
                if (mina >= minp)&&(minp < umbralsig*absmax*Tmax_Tmin_bifasic),
                    minp = [];
                elseif (minp> mina)&&(mina < umbralsig*absmax*Tmax_Tmin_bifasic),
                    mina = [];
                end
            end
            
            % Test which modulus maxima are significative and find zero crossings
            if isempty(mina),
                if isempty(minp),
                    tipoT = 2;      % only upwards T wave
                    if maxpos - minapos > 2,         % if not !!!!!!!!!!!
                        ind = zerocros(flipud(w(minapos:maxpos )));   %Scale 3 !!! % also in scale 4!!!!!!! 03.Dec.04
                        T = maxpos - ind +1;                           % Zero crossing = T wave position
                        picoff = maxpos;                               %wavelet  peak to detect offset
                    elseif isempty(minapos);  % If there were no minimum, there is no zero crossing
                        T = picant (w(begwin:maxpos ),maxpos);       % Take the minimum at scale 4
                        if ~isempty(T) %%%%%%%%%%%% 14/06/02 Rute
                            picoff = maxpos;
                        else
                            picoff = []; % if did not exist a peak in scale 4 there is no T
                        end %%%%%%%%%%%% 14/06/02 Rute
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave.'}];
                    end
                else      % minp exists (is significative) but mina not
                    tipoT = 0;  		%normal T wave
                    if minppos -maxpos >2,       % if not!!!!!!???
                        ind = zerocros(w(maxpos:minppos ));  %% 07/06/02 Rute % also in scale 4!!!!!!! 03.Dec.04
                        %                         if isempty(ind) %%%%%%%%%%%%%%Rute 03/09/02
                        %                             ind = zerocros(w(maxpos:minppos));
                        %                         end %%%%%%%%%%%%%%Rute 03/09/02
                        T = maxpos + ind -1;
                        picon = maxpos;		% For determining onset and offset
                        picoff = minppos;
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave.'}];
                    end
                end
            else
                if isempty(minp),   %  mina exists (is significative) but minp not
                    tipoT = 1;  	%inverted T wave
                    if maxpos -minapos >2,  % if not !!!!!!!!!
                        ind = zerocros(w(minapos:maxpos));  % wavelet zero crossing is T wave peak  %% 07/06/02 Rute % also in scale 4!!!!!!! 03.Dec.04
                        %                         if isempty(ind)%%%%%%%%%%%%%%Rute 03/09/02
                        %                             ind = zerocros(w(minapos:maxpos ));  % wavelet zero crossing is T wave peak  %% 03/09/02 Rute
                        %                         end %%%%%%%%%%%%%%Rute 03/09/02
                        T = minapos + ind -1;
                        picon = minapos;
                        picoff = maxpos;
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave.'}];
                    end
                else    % both mina and minp are significative.  Biphasic wave.
                    tipoT = 5;	% biphasic neg-pos T wave
                    if maxpos - minapos > 2,    %!!!!!!!!!!!!
                        ind = zerocros(flipud(w(minapos:maxpos ))); %% 07/06/02 Rute % also in scale 4!!!!!!! 03.Dec.04
                        %                         if isempty(ind) %%%%%%%%%%%%%%Rute 03/09/02
                        %                             ind = zerocros(flipud(w(minapos:maxpos ))); %% 03/09/02 Rute
                        %                         end%%%%%%%%%%%%%%Rute 03/09/02
                        T = maxpos - ind +1;
                        picon = minapos;
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave.'}];
                    end
                    if minppos - maxpos > 2,    %!!!!!!!!!!!!
                        ind = zerocros(flipud(w(maxpos:minppos ))); %% 07/06/02 Rute % also in scale 4!!!!!!! 03.Dec.04
                        %                         if isempty(ind)%%%%%%%%%%%%%%Rute 03/09/02
                        %                             ind = zerocros(flipud(w(maxpos:minppos ))); %% 03/09/02 Rute
                        %                         end %%%%%%%%%%%%%%Rute 03/09/02
                        %Tprima = maxpos + ind -1; %18.11.05
                        Tprima = minppos- ind +1;
                        picoff = minppos;
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave.'}];
                    end
                end
            end
        else        % If the greatest modulus maximum is the minimum
            % Search two maxima, one before and one after the minimum
            maxapos = max(modmax(w(begwin+1:minpos-1 ),2,0,1));
            maxapos = begwin + maxapos;
            if isempty(maxapos) && (minpos ~= begwin) && (w(begwin )>0),
                maxapos = begwin;
            end
            maxppos = min(modmax(w(minpos+1:endwin ),2,0,1));
            maxppos = minpos + maxppos;
            if isempty(maxppos) && (minpos ~= endwin) && (w(endwin )>0),
                maxppos = endwin;
            end
            maxa = abs(w(maxapos ));
            maxp = abs(w(maxppos )) ;                                   % See if they are significative
            if (maxa < umbralsig*absmin)
                maxa =[];
            elseif (minpos-maxapos>Tmax_Tmin_time_min*freq),
                maxa = [];
            end
            if (maxp < umbralsig*absmin),
                maxp =[];
            elseif (maxppos-minpos>Tmax_Tmin_time_min*freq),
                maxp = [];
            end
            if ~isnan(maxa)&~isnan(maxp),      %#ok<AND2> %%% NUEVO JP
                if (maxa >= maxp)&&(maxp < umbralsig*absmin*Tmax_Tmin_bifasic),
                    maxp = [];
                elseif (maxp> maxa)&&(maxa < umbralsig*absmin*Tmax_Tmin_bifasic),
                    maxa = [];
                end
            end
            % Test which modulus maxima are significative and find zero crossings
            if isempty(maxa),
                if isempty(maxp),
                    tipoT = 3;      % only downwards T wave
                    if minpos - maxapos > 2, %!!!!!!!!
                        ind = zerocros(flipud(w(maxapos:minpos )));   %Scale 3 % also in scale 4!!!!!!! 03.Dec.04
                        %                         if isempty(ind) %%%%%%%%%%%%%%Rute 03/09/02
                        %                             ind = zerocros(flipud(w(maxapos:minpos )));
                        %                         end %%%%%%%%%%%%%%Rute 03/09/02
                        T = minpos - ind +1;
                        picoff = minpos;
                    elseif isempty(maxapos);  % If there were no maximum, there is no zero crossing.
                        T = picant (w(begwin:minpos ),minpos);  % menimo en escala 4.
                        if ~isempty(T) %%%%%%%%%%%% 14/06/02 Rute
                            picoff = minpos;
                        else
                            picoff = []; % if did not exist a peak in scale 4 there is no T
                        end %%%%%%%%%%%% 14/06/02 Rute
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave.'}];
                    end
                else      % maxp is signficative, but not maxa
                    tipoT = 1;  %inverted T wave
                    if maxppos -minpos >2,  % !!!!!!
                        ind = zerocros(w(minpos:maxppos ));  %% 07/06/02 Rute % also in scale 4!!!!!!! 03.Dec.04
                        %                         if isempty(ind) %%%%%%%%%%%%%%Rute 03/09/02
                        %                             ind = zerocros(w(minpos:maxppos ));  %% 03/09/02 Rute
                        %                         end %%%%%%%%%%%%%%Rute 03/09/02
                        T = minpos + ind -1;
                        picon = minpos;		% For calculating onset and offset
                        picoff = maxppos;
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave.'}];
                    end
                end
            else
                if isempty(maxp),   % maxa is significative, but not maxp
                    tipoT = 0;       %normal T wave
                    if minpos -maxapos >2,
                        ind = zerocros(w(maxapos:minpos )); %% 07/06/02 Rute % also in scale 4!!!!!!! 03.Dec.04
                        %                         if isempty(ind) %%%%%%%%%%%%%%%%%%%% 03/09/02 Rute
                        %                             ind = zerocros(w(maxapos:minpos )); %% 03/09/02 Rute
                        %end %%%%%%%%%%%%%%%%%%%% 03/09/02 Rute
                        T = maxapos + ind -1;
                        picon = maxapos;
                        picoff = minpos;
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave.'}];
                    end
                else    % both maxa and maxp are significative.  Biphasic wave.
                    tipoT = 4;	% biphasic pos-neg T wave
                    if minpos - maxapos > 2,  %!!!!!!!!!!!
                        ind = zerocros(flipud(w(maxapos:minpos ))); %% 07/06/02 Rute % also in scale 4!!!!!!! 03.Dec.04
                        %                         if isempty(ind) %%%%%%%%%%%%%%%%%%%% 03/09/02 Rute
                        %                             ind = zerocros(flipud(w(maxapos:minpos ))); %% 03/09/02 Rute
                        %                         end %%%%%%%%%%%%%%%%%%%% 03/09/02 Rute
                        T = minpos - ind +1;
                        picon = maxapos;
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        %                         [i tipoT]
                        disp('caso nao previsto; onda T nao marcada!')
                    end
                    if maxppos - minpos > 2,    %!!!!!!!!!!!!
                        ind = zerocros(flipud(w(minpos:maxppos ))); %% 07/06/02 Rute % also in scale 4!!!!!!! 03.Dec.04
                        %                         if isempty(ind) %%%%%%%%%%%%%%%%%%%% 03/09/02 Rute
                        %                             ind = zerocros(flipud(w(minapos:maxpos ))); %% 03/09/02 Rute
                        %                         end %%%%%%%%%%%%%%%%%%%% 03/09/02 Rute
                        %Tprima = minpos + ind -1; %18.11.05 DUVIDA
                        Tprima = maxppos- ind +1;
                        picoff = maxppos;
                    else %%%%%%%%%%%%%%%%%%%%% nao faz nada!!!!!!! % 12/06/02  Rute
                        messages.warnings=[messages.warnings {'unknown case: unable to classify T wave.'}];
                    end
                end
            end
        end
        
        % T wave onset and offset detection
        %if isempty(T), picon=[]; picoff=[]; end
        
        picon_keep=[picon_keep picon]; %multileadchange
        picoff_keep=[picoff_keep picoff];%multileadchange
        lead_keep=[lead_keep 4];%multileadchange
        
        if ~isempty(picon),
            Ton = searchon (picon, w(max(begwin,picon-0.12*freq):picon ), Kton);
        end
        if ~isempty(picoff),
            Toff=searchoff(picoff, w(picoff:min([size(w,1) picoff+0.12*freq ])  ) , Ktoff);
            if (Toff > endwin),
                Toff = endwin;
            end
        end
    else
        tipoT=9;
        messages.warnings=[messages.warnings {'unknown case: unable to find T wave in multilead approach using scale 4.'}];
    end       %% utilizar a escala 5 07/06/02 Rute %%%%%%%%%%%%%%%%%%%%%%%%%%%
    
end
% Filling the structure with positions
if isempty(Ton), Ton = NaN; end;
if isempty(Toff), Toff = NaN;  end;
if isempty(T), T=NaN; end;
if isnan(T),
    picon_keep=[picon_keep NaN]; %multileadchange
    picoff_keep=[picoff_keep NaN];%multileadchange
    lead_keep=[lead_keep NaN];%#ok<NASGU> %multileadchange
end
if isempty(Tprima), Tprima=NaN; end;
if isempty(tipoT), tipoT=NaN; end;
pos.Ton= Ton+samp(1)-1+timenew(i)-1;
pos.Toff= Toff+samp(1)-1+timenew(i)-1;
pos.T= T+samp(1)-1+timenew(i)-1;
pos.Tprima= Tprima+samp(1)-1+timenew(i)-1;
pos.Ttipo= tipoT;
% T = []; Tprima=[]; picon=[]; picoff=[]; Ton=[]; Toff=[]; tipoT=[];
% minapos = []; minppos=[]; maxapos=[]; maxppos=[]; mina=[]; minp=[]; maxa=[];
% maxp=[];
%end
% position.Ton(intervalo(1):intervalo(2)) = pos.Ton+samp(1)-1;
% position.Toff(intervalo(1):intervalo(2))= pos.Toff+samp(1)-1;
% position.T(intervalo(1):intervalo(2)) = pos.T+samp(1)-1;
% position.Tprima(intervalo(1):intervalo(2)) = pos.Tprima+samp(1)-1;
% position.Ttipo(intervalo(1):intervalo(2)) = pos.Ttipo;
% pos.Ton
% Ton = pos.Ton(i)
% Toff= pos.Toff(i)
% T= pos.T(i)
% Tprima = pos.Tprima(i)
% Ttipo= pos.Ttipo(i);
%timenew(i)
%picon_keep
picoff=picoff_keep+timenew(i)-1;
picon=picon_keep+timenew(i)-1;