% -------------------------------------------------------------------------------------------------
% BeatGroupingSC.m: Attempt to cluster ECG beats and average beat estimation.
%                   Rough version to be tested, tuned and refined.
%
% Copyright (C) 2004 Maurizio Varanini, Clinical Physiology Institute, CNR, Pisa, Italy
%
% This program is free software; you can redistribute it and/or modify it under the terms
% of the GNU General Public License as published by the Free Software Foundation; either
% version 2 of the License, or (at your option) any later version.
%
% This program is distributed "as is" and "as available" in the hope that it will be useful,
% but WITHOUT ANY WARRANTY of any kind; without even the implied warranty of MERCHANTABILITY
% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License along with this program;
% if not, write to the Free Software Foundation, Inc.,
% 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
%
% For any comment or bug report, please send e-mail to: maurizio.varanini@ifc.cnr.it
% -------------------------------------------------------------------------------------------------

function beatGroup= beatGroupingSC(x1, viq, npt)

vvqp=1:length(viq);
thrp=1;
thr=1.1;
%thr=0.5;
tha=0.2;
%tha=0.1;
ic=0;
while(~isempty(vvqp));
    ic=ic+1;
    viqcp=viq(vvqp);
    nvq=length(vvqp);
    % build a rough centroid of selected beats
    avgP1Q= roughAvgClass(x1, viqcp, npt);
    % evaluate the distance from the mean template
    [mid da]=dist2cAvg(x1, avgP1Q, viqcp);
    dam=meansc(mid,25,25);   % mean distance
    if(ic==1) tha=0.2*dam; end
    vqp=find(mid < tha +thrp*dam);  % array of index of similar beats
    viqp=viqcp(vqp);
    % build the centroid of selected beats
    avgC1Q(:,ic)= avgClass(x1, viqp, npt);
    % evaluate the distance from the centroid
%    [mid, mij, da]=dist2cAvgW(x1, x2, avgC1Q(:,ic), avgC2Q(:,ic), viqcp);
    [mid, mij, da, mih]=dist2cAvgWh(x1, avgC1Q(:,ic), viqcp);
    dam=meansc(mid,25,25);    % mean distance
    iq=find(mid < tha+thr*dam);  % array of index of similar beats
    vq{ic}=vvqp(iq);
    vij(vq{ic})=mij(iq); 
    vih(vq{ic})=mih(iq);
    vqcc(vq{ic})=ic;
    vvqp=vvqp(find(mid >= tha+thr*dam));
end
beatGroup.vq=vq;
beatGroup.vij=vij;
beatGroup.vih=vih;
beatGroup.vqcc=vqcc;
beatGroup.avgC1Q=avgC1Q;

% end function
%------------------------------------------------------------
% evaluate the distance from the mean template
function [mid, da]= dist2cAvg(x1, avg1Q, iw )
npt=length(avg1Q);
nq=length(iw);
for iq=1:nq
    mid(iq)=999999999999;
    j=iw(iq);
    d1= x1(j:j+npt-1)-avg1Q;
    da=sum(abs(d1));
    if(da<mid(iq)) mid(iq)=da;  end
end
% end function
%------------------------------------------------------------
% evaluate the distance from the mean template
function [mid, mij, da]=dist2cAvgW(x1, avg1Q, iw )
npt=length(avg1Q);
nq=length(iw);
for iq=1:nq
    mid(iq)=999999999999;
    ij=iw(iq)-20;
    fj=ij+20;
    for j=ij:fj
        d1= x1(j:j+npt-1)-avg1Q;
        da=sum(abs(d1));
        if(da<mid(iq)) mid(iq)=da; mij(iq)=j; end
    end
end
% end function
%------------------------------------------------------------
% evaluate the distance from the mean template
function [mid, mij, da, mih]=dist2cAvgWh(x1, avg1Q, iw )
npt=length(avg1Q);
nq=length(iw);
   mind = min(avg1Q);
   maxd = max(avg1Q);

for iq=1:nq
    mid(iq)=999999999999;
    ij=iw(iq)-20;
    fj=ij+20;
    for j=ij:fj
        for h=linspace((-maxd+mind)/10,(maxd-mind)/10,20)
            d1= x1(j:j+npt-1)-avg1Q+h;
            da=max(abs(d1));
            if(da<mid(iq)) mid(iq)=da; mij(iq)=j; mih(iq)=h; end
        end
    end
end
% end function
%------------------------------------------------------------
% build the centroid of selected beats
function avg1Q= avgClass(x1, viq, npt)
nnq=length(viq);
avg1Q=0;
for jq=1:nnq
    ij=viq(jq);
    fj=ij+npt-1;
    avg1Q= avg1Q+x1(ij:fj);
end
avg1Q=avg1Q/nnq;
% end function

%------------------------------------------------------------
% Build a rough centroid of selected beats
function avg1Q= roughAvgClass(x1, viq, npt)
% built matrix of QRS templates
nq=length(viq);
for iq=1:nq
    ji=viq(iq); jf=ji+npt-1;
    x1M(:,iq)= x1(ji:jf);
end
for is=1:npt
    avg1Q(is,1)=meansc(x1M(is,:),10,10);
end
% end function