function res = spm_eeg_artefact_heartbeat(S)
% Detects eyeblinks in spm continuous data file
% S                     - input structure
% fields of S:
%    S.D                - M/EEG object
%    S.chanind          - vector of indices of channels that this plugin will look at.
%    S.threshold        - threshold parameter (in stdev)
%
%    Additional parameters can be defined specific for each plugin
% Output:
%  res -
%   If no input is provided the plugin returns a cfg branch for itself
%
%   If input is provided the plugin returns a matrix of size D.nchannels x D.ntrials
%   with zeros for clean channel/trials and ones for artefacts.
%______________________________________________________________________________________
% Copyright (C) 2008-2013 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% see http://fsl.fmrib.ox.ac.uk/eeglab/fmribplugin/
% $Id: spm_eeg_artefact_heartbeat.m 5592 2013-07-24 16:25:55Z vladimir $


%-This part if for creating a config branch that plugs into spm_cfg_eeg_artefact
% Any parameters can be specified and they are then passed to the plugin
% when it's called.
%--------------------------------------------------------------------------
if nargin == 0           
    excwin = cfg_entry;
    excwin.tag = 'excwin';
    excwin.name = 'Excision window';
    excwin.strtype = 'r';
    excwin.num = [1 1];
    excwin.val = {0};
    excwin.help = {'Window (in ms) to mark as bad around each heart beat, 0 to not mark data as bad'};
    
    heartbeat = cfg_branch;
    heartbeat.tag = 'heartbeat';
    heartbeat.name = 'Heart beats';
    heartbeat.val = {excwin};
    
    res = heartbeat;
    
    return
end

SVNrev = '$Rev: 5592 $';

%-Startup
%--------------------------------------------------------------------------
spm('sFnBanner', mfilename, SVNrev);
spm('FigName','M/EEG heartbeat detection');

if exist('fmrib_qrsdetect', 'file')~=2
    error('This tool requires FMRIB plugin, see http://fsl.fmrib.ox.ac.uk/eeglab/fmribplugin/');
end


if isequal(S.mode, 'reject')
    error('Only mark mode is supported by this plug-in, use event-based rejection to reject.');
end

D = spm_eeg_load(S.D);

chanind  =  S.chanind;

if length(chanind)~=1
    error('More than one channel - not currently supported')
end

% Detect QRS peaks using FMRIB plugin
%--------------------------------------------------------------------------
EEG = [];
EEG.data =  reshape(squeeze(D(chanind,:,:)), 1, []);
EEG.srate = D.fsample;
spikes = fmrib_qrsdetect(EEG,1);


% Update the event structure
%----------------------------------------------------------------------
if ~isempty(spikes)  
    for n = 1:D.ntrials
        cspikes   = spikes(spikes>(D.nsamples*(n-1)) & spikes<(D.nsamples*n));
        ctime  = D.trialonset(n)+(cspikes - D.nsamples*(n-1)-1)/D.fsample;
        ctime  = num2cell(ctime);
        
        ev = events(D, n);
        
        if iscell(ev)
            ev = ev{1};
        end
        
        
        if ~isempty(ev) && ~S.append
            ind1 = strmatch('artefact_heartbeat', {ev.type}, 'exact');
            if ~isempty(ind1)
                ind2 = strmatch(D.chanlabels(chanind), {ev(ind1).value}, 'exact');
                if ~isempty(ind2)
                    ev(ind1(ind2)) = [];
                end
            end
        end        
        
        Nevents = numel(ev);
        for i=1:numel(ctime)
            if ctime{i} == 0
                continue; %likely to be trial border falsely detected as heartbeat
            end
            ev(Nevents+i).type     = 'artefact_heartbeat';
            ev(Nevents+i).value    = char(D.chanlabels(chanind));
            if S.excwin == 0
                ev(Nevents+i).duration = [];
                ev(Nevents+i).time     = ctime{i};
            else
                ev(Nevents+i).time     = max(D.trialonset(n), ctime{i} - 5e-4*S.excwin);
                ev(Nevents+i).duration = min(1e-3*S.excwin, (D.time(end)-D.time(1))-(ev(Nevents+i).time-D.trialonset(n)))+...
                    min(ctime{i} - 5e-4*S.excwin, 0);               
            end                      
        end
        
        if ~isempty(ev)
            [dum, I] = sort([ev.time]);
            ev = ev(I);
            D = events(D, n, ev);
        end
    end    
else
    warning(['No heartbeat events detected in the selected channel']);
end

res = D;

spm('FigName','M/EEG heartbeat detection: done');