function Align = FindFiducials(ecg, q, s, Align)
% FindFiducials.m
% Author: Alexander Khaustov; alexander dot khaustov at gmail dot com 
% Copyright (C) 2008 St.-Petersburg Institute of Cardiological Technics (Incart), www.incart.ru
% This software is released under the terms of the GNU General
% Public License (http://www.gnu.org/copyleft/gpl.html).
% 
% finds the lead with maximum QRS amplitude, takes a beat as QRS template and finds fiducial points such that
% aligned beats have maximum cross-correlation (over Q to S or S to T interval) in the chosen
% lead against the template
%
% no attempt has been made to make this code computationally effective

global Param
q = q(1:min(length(q), length(s)));
s = s(1:min(length(q), length(s)));

leadSorted = SortLeads(ecg, q, s, Align.st);

for ilead = 1:length(leadSorted)
    Align.lead = leadSorted(ilead);

    Align.orientation = MaxFrontOrientation(ecg(:, Align.lead), q, s);

    Align.fidBase = FindFidBase(ecg(:, Align.lead), q, s, Align.orientation);


    Align.f2s = floor(median(s - Align.fidBase) + Param.stAdjIntv);   % reserve to avoid including part of QRS in ST interval computations
    Align.q2f = floor(median(Align.fidBase - q));
    Align.fid = [];
    NFailed = 0;

    Align.template = 0;
    while isempty(Align.fid)
        %    template = floor(rand * length(s)) + 1;
        Align.template = Align.template + 1;
        %disp(['template beat - ' num2str(Align.template)]);

        Align = Adjust(ecg, Align);
        if isempty(Align.fid)

            %       disp('doesn''t do...');
            NFailed = NFailed + 1;
            if (NFailed > 0.1 * length(s))
                %          disp('what a mess..');
                break;
            end;
        end;
    end;
    if (~isempty(Align.fid))
        for i = 1:length(Align.fid)
            for lead=1:size(ecg, 2)
                Align.amp(i, lead) = mean(ecg(Align.fidQRS(i) - Align.q2f : Align.fidQRS(i) + Align.f2s, lead));  % align to zero mean on QS segment
            end;
        end;
        Align = CheckLeadCorrelations(ecg, Align);
        disp(['FindFiducials: Alignment in lead ' num2str(Align.lead)]);
        break;
    end;
    disp(['FindFiducials: lead ' num2str(Align.lead) ' is too noisy or ectopic..']);
end;
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function leads = SortLeads(ecg, q, s, st)

NLead = size(ecg, 2);

global Param
%   QRS or T amplitudes in each lead
ptp = zeros(min(length(q), 8),NLead);
for i = 1:min(length(q), 8)
    for l = 1:NLead
        if strcmp(Param.Alignment, 'st')
            ptp(i, l) = max(ecg(s(i):s(i) + st, l)) - min(ecg(s(i):s(i) + st, l));
        else
            ptp(i, l) = max(ecg(q(i):s(i), l)) - min(ecg(q(i):s(i), l));
        end;
    end;
end;
ptpmed = zeros(1,NLead);
for i = 1:NLead
    ptpmed(i) = median(ptp(:, i));
end;

[sorted indices] = sort(ptpmed);
leads = indices(length(indices):-1:1);
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Align = Adjust(ecg, Align)

global Param

qs_templ = ecg((Align.fidBase(Align.template) - Align.q2f):(Align.fidBase(Align.template) + Align.f2s), Align.lead);
st_templ = ecg((Align.fidBase(Align.template) + Align.f2s):(Align.fidBase(Align.template) + Align.f2s + Align.st), Align.lead);

[a b c d] = AdjustFiducials(ecg, Align.fidBase, Align.q2f, Align.f2s, qs_templ, Align.st, st_templ);

if ~isempty(c)
    Align.fidQRS = a;
    Align.QRScorr(:, Align.lead) = b;
    Align.fid = c;
    Align.Tcorr(:, Align.lead) = d;
    amp = mean(qs_templ);
    Align.qs_templ(:, Align.lead) = qs_templ - amp;
    Align.st_templ(:, Align.lead) = st_templ - amp;
    
    for i = 1:length(Align.fid)
        Align.valid(i, Align.lead) = (Align.QRScorr(i, Align.lead) >= Param.corrQRS) && (Align.Tcorr(i, Align.lead) >= Param.corrT);
    end;
    Align.validleads(Align.lead) = true;
end;

return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function orientation = MaxFrontOrientation(ecg, q, s)

NPos = 0;

for i = 1:min(20, length(s))
    [ma maI] = max(ecg(q(i):s(i)));
    [mi miI] = min(ecg(q(i):s(i)));

    NPos = NPos + (maI > miI);
end;

orientation = (NPos > 10);

return;