function [event] = read_fcdc_event(filename)

% READ_FCDC_EVENT reads all events from an EEG/MEG dataset and returns them
% in a well defined structure. It is a wrapper around different EEG/MEG file
% importers, directly supported formats are CTF, Neuromag, EEP, BrainVision, 
% Neuroscan and Neuralynx.
%
% Use as
%   [event] = read_fcdc_event(filename)
%
% This function returns a structure with the following fields
%   event.type     = string
%   event.sample   = expressed in samples, first sample of file is 1
%   event.value    = number or string
%   event.offset   = expressed in samples
%   event.duration = expressed in samples
%
% Some of these fields can be empty, depending on the type of event file.
% After reading the event structure, you can use following tricks to
% extract information about those events in which you are interested.
%
% Determine the different event types
%   unique({event.type})
%
% Get the index of all trial events
%   find(strcmp('trial', {event.type}))
%
% Make a vector with all triggers that occurred on the backpanel
%   [event(find(strcmp('backpanel trigger', {event.type}))).value]
%
% Find the events that ocurred in trial t=26
%   t = 26; find([event.sample]>trl(t,1) & [event.sample]<trl(t,2))
%
% See also READ_FCDC_DATA, READ_FCDC_HEADER

% Copyright (C) 2004, Robert Oostenveld
%
% $Log: read_fcdc_event.m,v $
% Revision 1.29  2005/10/05 06:30:25  roboos
% added support for MPI datasets and DAP files
%
% Revision 1.28  2005/09/15 07:38:38  roboos
% added support for ctf_res4 and ctf_meg4 by renaming filename into dataset
%
% Revision 1.27  2005/09/08 09:25:53  roboos
% added support for the Yokogawa MEG data format
%
% Revision 1.26  2005/09/06 12:45:00  roboos
% renamed plextor into plexon (incorrect company name)
%
% Revision 1.25  2005/08/17 19:34:58  roboos
% added a check for the presence of the file
% added a TTLValue event for each trigger in Neuralynx dataset
%
% Revision 1.24  2005/06/17 10:56:59  roboos
% fixed event.type for CTF stimulus channels (should be string instead of cell)
%
% Revision 1.23  2005/05/24 07:37:04  roboos
% implemented reading of other trigger channels (like UPPT001, UTRG001) for CTF event detection
%
% Revision 1.22  2005/05/19 07:11:32  roboos
% added support for neuralynx_nev
%
% Revision 1.21  2005/05/18 15:50:31  roboos
% repeated last bugfix, now hopefully correct (thanks to Juliet)
%
% Revision 1.20  2005/05/18 10:26:43  roboos
% fixed "Operands to the || and && operators must be convertible to logical scalar values" problem by first checking whether variable is a string or numeric
%
% Revision 1.19  2005/05/17 17:50:38  roboos
% changed all "if" occurences of & and | into && and ||
% this makes the code more compatible with Octave and also seems to be in closer correspondence with Matlab documentation on shortcircuited evaluation of sequential boolean constructs
%
% Revision 1.18  2005/05/11 08:27:50  roboos
% implemeted the correct offset (in samples) for CTF markerfiles, sofar the offset was always set to zero
%
% Revision 1.17  2005/03/31 07:05:58  roboos
% implemented support for besa_avr using my own read_besa_avr
%
% Revision 1.16  2005/02/16 07:48:59  roboos
% added support for Plexon nex files, using read_nex_event
%
% Revision 1.15  2005/02/02 14:33:16  roboos
% renamed event.type of eep_avg to 'average', implemented triggers for eep_cnt (using *.trg file)
%
% Revision 1.14  2004/12/07 18:13:04  roboos
% added initial support for Neuroscan cnt files
%
% Revision 1.13  2004/11/17 09:01:50  roboos
% added check for the presence of CTF MarkerFile, to avoid unneccesary warning in readmarkerfile
%
% Revision 1.12  2004/11/15 09:14:09  roboos
% implemented support for Neuroscan *.eeg (epoched) files
% applied auto-indentation, which changed some whitespace in other parts of the file
%
% Revision 1.11  2004/10/01 09:53:25  roboos
% added Neuromag to the help as one of the supported formats
%
% Revision 1.10  2004/09/27 15:37:11  roboos
% renamed the function readClassFile to read_ctf_cls (part of file_io)
%
% Revision 1.9  2004/09/27 13:55:24  roboos
% fixed bugs for CTF marker
%
% Revision 1.8  2004/09/24 07:13:57  roboos
% added support for CTF Markerfiles, using the readmarkerfile function I got from Tom Holroyd
%
% Revision 1.7  2004/09/21 13:09:46  roboos
% added 2 digital trigger channels for fif, analog triggers in fif file are only scanned if no digital triggers are present
%
% Revision 1.6  2004/09/20 11:48:31  roboos
% improved support for Neuromag, threshold all STI-channels at =>5 and treat STI014 as digital channel
%
% Revision 1.5  2004/09/03 09:02:29  roboos
% added raw>5 for flank detection for neuromag trigger channels
%
% Revision 1.4  2004/08/20 09:18:14  roboos
% changed event(i).sample for CTF trials
% implemented Neuromag format: trials and "STI xxx" channels
%
% Revision 1.3  2004/06/30 16:05:24  roboos
% added support for BrainVision marker file on a more general way than implemented in the read_brainvision_vmrk function
% added subfunction tokenize
% renamed dataset into filename
%
% Revision 1.2  2004/06/21 19:21:10  roberto
% renamed eventfile into dataset, added check for empty CTF classfile
%
% Revision 1.1  2004/06/03 09:15:22  roberto
% initial version, supported are CTF and eep_avr
%

% detect the external toolboxes that optionally can be used
hasbiosig = (exist('sopen') & exist('sread'));                 % see http://biosig.sourceforge.net/
hasmegpd  = (exist('rawdata') & exist('channames'));           % see http://www.kolumbus.fi/kuutela/programs/meg-pd/
haseegsf  = (exist('ctf_read_meg4') & exist('ctf_read_res4')); % see http://eeg.sourceforge.net/

% test whether the file exists
if ~exist(filename)
  error(sprintf('file ''%s'' does not exist', filename));
end

% start with an empty event structure
event = [];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if filetype(filename, 'ctf_ds') || filetype(filename, 'ctf_meg4') || filetype(filename, 'ctf_res4')

  % obtain the dataset name
  if filetype(filename, 'ctf_meg4') ||  filetype(filename, 'ctf_res4')
    filename = fileparts(filename); 
  end

  [path, name, ext] = fileparts(filename);
  headerfile = fullfile(path, [name ext], [name '.res4']);
  datafile   = fullfile(path, [name ext], [name '.meg4']);
  classfile  = fullfile(path, [name ext], 'ClassFile.cls');

  hdr = read_ctf_res4(headerfile);

  % read the trigger codes from the STIM channel, usefull for (pseudo) continuous data
  % this splits the trigger channel into the lowers and highest 16 bits, 
  % corresponding with the front and back panel of the electronics cabinet at the Donders Centre
  [backpanel, frontpanel] = read_ctf_trigger(filename);
  for i=find(backpanel(:)')
    event(end+1).type   = 'backpanel trigger';
    event(end  ).sample = i;
    event(end  ).value  = backpanel(i);
  end
  for i=find(frontpanel(:)')
    event(end+1).type   = 'frontpanel trigger';
    event(end  ).sample = i;
    event(end  ).value  = frontpanel(i);
  end

  % determine the trigger channels from the header
  if isfield(hdr, 'sensType')
    for i=find(hdr.sensType(:)'==11)
      % read the trigger channel as raw data
      trig = read_ctf_meg4(datafile, hdr, 1, hdr.nTrials*hdr.nSamples, i);
      % correct for reading it as signed integer, whereas it should be an unsigned int
      trig(find(trig<0)) = trig(find(trig<0)) + 2^32;
      % convert the trigger into an event with a value at a specific sample
      for j=find(diff([0 trig(:)'])>1)
        event(end+1).type   = hdr.label{i};
        event(end  ).sample = j;
        event(end  ).value  = trig(j);
      end
    end
  end

  % make an event for each trial as defined in the header
  for i=1:hdr.nTrials
    event(end+1).type     = 'trial';
    event(end  ).sample   = (i-1)*hdr.nSamples + 1;
    event(end  ).offset   = -hdr.nSamplesPre;
    event(end  ).duration =  hdr.nSamples;
  end

  % read the classification file and make an event for each classified trial
  [condNumbers,condLabels] = read_ctf_cls(classfile);
  if ~isempty(condNumbers)
    Ncond = length(condLabels);
    for i=1:Ncond
      for j=1:length(condNumbers{i})
        event(end+1).type     = 'classification';
        event(end  ).value    = condLabels{i};
        event(end  ).sample   = (condNumbers{i}{j}-1)*hdr.nSamples + 1;
        event(end  ).offset   = -hdr.nSamplesPre;
        event(end  ).duration =  hdr.nSamples;
      end
    end
  end

  if exist(fullfile(filename,'MarkerFile.mrk'))
    % read the marker file and make an event for each marker
    % this depends on the readmarkerfile function that I got from Tom Holroyd
    % I have not tested this myself extensively, since at the FCDC we
    % don't use the marker files
    mrk = readmarkerfile(filename);
    for i=1:mrk.number_markers
      for j=1:mrk.number_samples(i)
        % determine the location of the marker, expressed in samples
        trialnum = mrk.trial_times{i}(j,1);
        synctime = mrk.trial_times{i}(j,2);
        begsample = (trialnum-1)*hdr.nSamples + 1;    % of the trial, relative to the start of the datafile
        endsample = (trialnum  )*hdr.nSamples;        % of the trial, relative to the start of the datafile
        offset    = round(synctime*hdr.Fs);           % this is the offset (in samples) relative to time t=0 for this trial
        offset    = offset + hdr.nSamplesPre;         % and time t=0 corrsponds with the nSamplesPre'th sample
        % store this marker as an event
        event(end+1).type    = mrk.marker_names{i};
        event(end ).value    = [];
        event(end ).sample   = begsample + offset;
        event(end ).duration = 0;
        event(end ).offset   = offset;
      end
    end
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif filetype(filename, 'eep_avr')
  % the headerfile and datafile are the same
  hdr = read_fcdc_header(filename);
  event(end+1).type     = 'average';
  event(end  ).sample   = 1;
  event(end  ).duration = hdr.nSamples;
  event(end  ).offset   = -hdr.nSamplesPre;
  event(end  ).value    = [];

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif filetype(filename, 'eep_cnt')
  % try to read external trigger file in EEP format
  trgfile = [filename(1:(end-3)), 'trg'];
  if exist(trgfile, 'file')
    hdr = read_fcdc_header(filename);
    tmp = read_eep_trg(trgfile);
    % translate the EEProbe trigger codes to events
    for i=1:length(tmp)
      event(i).type     = 'trigger';
      event(i).sample   = round((tmp(i).time/1000) * hdr.Fs) + 1;    % convert from ms to samples
      event(i).value    = tmp(i).code;
      event(i).offset   = 0;
      event(i).duration = 0;
    end
  else
    warning('no triggerfile was found');
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif filetype(filename, 'brainvision_vhdr') || filetype(filename, 'brainvision_vmrk')

  if filetype(filename, 'brainvision_vhdr')
    % read the headerfile belonging to the dataset and try to determine the corresponding markerfile
    hdr = read_brainvision_vhdr(filename);

    % replace the filename with the filename of the markerfile
    if ~isfield(hdr, 'MarkerFile') || isempty(hdr.MarkerFile)
      filename = [];
    else
      [p, f, e] = fileparts(filename);
      filename = fullfile(p, hdr.MarkerFile);
    end
  end

  fid=fopen(filename,'rt');
  if fid==-1,
    error('cannot open BrainVision marker file')
  end

  line = [];
  while ischar(line) || isempty(line)
    line = fgetl(fid);
    if ~isempty(line) && ~(isnumeric(line) && line==-1)
      if strncmpi(line, 'Mk', 2)
        % this line contains a marker
        tok = tokenize(line, '=');
        if length(tok)~=2
          warning('skipping unexpected formatted line in BrainVision marker file');
        else
          % the line looks like "MkXXX=YYY", which is ok
          % the interesting part now is in the YYY, i.e. the second token
          tok = tokenize(tok{2}, ',');
          if isempty(tok{1})
            tok{1}  = [];
          end
          if isempty(tok{2})
            tok{2}  = [];
          end
          event(end+1).type     = tok{1};
          event(end  ).value    = tok{2};
          event(end  ).sample   = str2num(tok{3});
          event(end  ).duration = str2num(tok{4});
          event(end  ).offset=[];   % added by H. F. Kwok, University of Birmingham Dec 2005 for compatibility with definetrial.m
        end
      end
    end
  end
  fclose(fid);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif filetype(filename, 'yokogawa_sqd') || filetype(filename, 'yokogawa_con') || filetype(filename, 'yokogawa_raw')
  event = read_yokogawa_event(filename);
  

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif filetype(filename, 'mpi_ds') || filetype(filename, 'mpi_dap')
  hdr = read_fcdc_header(filename);
  % determine the DAP files that compromise this dataset
  if isdir(filename)
    ls = dir(filename);
    dapfile = {};
    for i=1:length(ls)
      if ~isempty(regexp(ls(i).name, '.dap$'))
        dapfile{end+1} = fullfile(filename, ls(i).name);
      end
    end
    dapfile = sort(dapfile);
  elseif iscell(filename)
    dapfile = filename;
  else
    dapfile = {filename};
  end
  % assume that each DAP file is accompanied by a dat file
  % read the trigger values from the separate dat files
  trg = [];
  for i=1:length(dapfile)
    datfile = [dapfile{i}(1:(end-4)) '.dat'];
    trg = cat(1, trg, textread(datfile, '', 'headerlines', 1));
  end
  % construct a event structure, one 'trialcode' event per trial
  for i=1:length(trg)
    event(i).type     = 'trialcode';            % string
    event(i).sample   = (i-1)*hdr.nSamples + 1; % expressed in samples, first sample of file is 1
    event(i).value    = trg(i);                 % number or string
    event(i).offset   = 0;                      % expressed in samples
    event(i).duration = hdr.nSamples;           % expressed in samples
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif filetype(filename, 'ns_eeg')
  hdr = read_fcdc_header(filename);
  for i=1:hdr.nTrials
    event(end+1).type     = 'trial';
    event(end  ).sample   = (i-1)*hdr.nSamples + 1;
    event(end  ).value    = [];
    event(end  ).offset   = -hdr.nSamplesPre;
    event(end  ).duration =  hdr.nSamples;
    % read the data to determine manually accepted/rejected trials
    tmp = read_ns_eeg(filename, i);
    if tmp.sweep.accept
      event(end).value = 'accept';
    else
      event(end).value = 'reject';
    end
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif filetype(filename, 'ns_cnt')
  warning('handling of event structure has not been fully tested for Neuroscan *.cnt files');
  % this is a very inefficient implementation
  % since it reads not only the event table, but also all data from the file
  tmp = read_ns_cnt(filename);
  % translate the cnt event table into known FieldTrip event types
  for i=1:tmp.nevent
    event(i).type     = 'trigger';
    event(i).sample   = tmp.event.frame(i);
    event(i).value    = tmp.event.stimtype(i);
    event(i).offset   = 0;
    event(i).duration = 0;
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif filetype(filename, 'besa_avr')
  hdr = read_fcdc_header(filename);
  event(end+1).type     = 'average';
  event(end  ).sample   = 1;
  event(end  ).duration = hdr.nSamples;
  event(end  ).offset   = -hdr.nSamplesPre;
  event(end  ).value    = [];

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif filetype(filename, 'neuromag_fif')
  hdr = read_fcdc_header(filename);
  % add the trials to the event structure
  for i=1:hdr.nTrials
    event(end+1).type     = 'trial';
    event(end  ).sample   = (i-1)*hdr.nSamples + 1;
    event(end  ).value    = [];
    event(end  ).offset   = -hdr.nSamplesPre;
    event(end  ).duration =  hdr.nSamples;
  end

  % add triggers based on the binary trigger channel, this is based the following email
  %
  % On 14 Sep 2004, at 16:33, Lauri Parkkonen wrote:
  % > Anke's file is probably quite recent one. It should have normal binary
  % > coding on STI014 and its "analog" counterparts on STI1 ... STI6,
  % > swinging between 0 ... +5 volts. These latter channels are virtual,
  % > i.e., they are generated from STI014 for backwards compatibility. STI1
  % > is the LSB of STI014 etc. For all the new stuff, I would recommend
  % > using STI014 as that's the way the triggers will be stored in the future
  % > (maybe the channel name is going to change to something more reasonable
  % > like TRIG). Unfortunately, some older files from the 306-channel system
  % > have the STI014 coded in a different way - the two 8-bit halves code
  % > the input and output triggers separately.

  triglab  = {'STI 014', 'STI 015', 'STI 016'};
  trigindx = match_str(hdr.label, triglab);
  if length(trigindx)>0
    % pad with the last sample of the previous block, start with zeros
    pad = zeros(length(trigindx),1);
    for i=1:hdr.nTrials
      begsample = (i-1)*hdr.nSamples + 1;
      endsample = (i  )*hdr.nSamples;
      raw = read_fcdc_data(filename, hdr, begsample, endsample, trigindx);
      % detect the flank of the trigger
      flank = ((diff([raw pad], [], 2) .* raw)~=0);
      % for each flank in the TRIG channels, add an event
      for j=1:length(trigindx)
        for k=find(flank(j,:))
          event(end+1).type     = 'trigger';
          event(end  ).sample   = begsample + k;
          event(end  ).value    = raw(j,k);
          event(end  ).offset   = [];
          event(end  ).duration = [];
        end
      end
      % remember the last sample to ensure continuity with the next block
      pad = raw(:,end);
    end
  end % if length(trigindx)

  % look for the analog triggers only if no binary trigger channel is present
  if isempty(trigindx)
    % add the triggers to the event structure based on trigger channels with the name "STI xxx"
    % this is not a very efficient implementation, since it has to read all data
    % furthermore, there are some issues with noise on these analog trigger channels
    stimindx = [];
    stimlab  = {};
    for i=1:length(hdr.label)
      if all(hdr.label{i}(1:3) == 'STI')
        stimindx(end+1) = i;
        stimlab{end+1}  = hdr.label{i};
      end
    end
    if length(stimindx)>0
      % pad with the last sample of the previous block, start with zeros
      pad = zeros(length(stimindx),1);
      for i=1:hdr.nTrials
        begsample = (i-1)*hdr.nSamples + 1;
        endsample = (i  )*hdr.nSamples;
        raw = read_fcdc_data(filename, hdr, begsample, endsample, stimindx);
        % detect the flank of the trigger
        flank = ((diff([raw pad], [], 2) .* raw)~=0) & (raw>=5);
        % Note: the original code read
        %   flank = (diff([raw pad], [], 2) .* raw)~=0;
        % but according to Joachim, real events always have triggers > 5
        % for each flank in the STIM channels, add an event
        for j=1:length(stimindx)
          for k=find(flank(j,:))
            event(end+1).type     = stimlab{j};
            event(end  ).sample   = begsample + k;
            event(end  ).value    = raw(j,k);
            event(end  ).offset   = [];
            event(end  ).duration = [];
          end
        end
        % remember the last sample to ensure continuity with the next block
        pad = raw(:,end);
      end
    end % if length(stimindx)
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif filetype(filename, 'plexon_nex')
  event = read_nex_event(filename);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif filetype(filename, 'neuralynx_nev') || filetype(filename, 'neuralynx_ds')
  if ~isdir(filename);
    filename = fileparts(filename); % only the directory should be used here
  end
  % read the header and the events
  hdr = read_neuralynx_header(filename);
  nev = read_neuralynx_event(filename);
  AcqStartTimeStamp = max(hdr.FirstTimeStamp); % the first sample that is present in all files is called "1"
  for i=1:length(nev)
    % add an event containing the string as value
    event(end+1).type     = 'EventString';
    event(end  ).sample   = round(hdr.SamplingFrequency * (nev(i).TimeStamp - AcqStartTimeStamp)/1e6) + 1;
    event(end  ).value    = nev(i).EventString;
    event(end  ).offset   = [];
    event(end  ).duration = [];
    % add an event containing the TTL as value
    event(end+1).type     = 'TTLValue';
    event(end  ).sample   = round(hdr.SamplingFrequency * (nev(i).TimeStamp - AcqStartTimeStamp)/1e6) + 1;
    event(end  ).value    = nev(i).TTLValue;
    event(end  ).offset   = [];
    event(end  ).duration = [];
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % added by H. F. Kwok for University of Birmingham data
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif filetype(filename,'bham_bdf')
    event=read_biosemi_event(filename);
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
else
  error('unknown filetype for events');
end % different filetypes

if ~isempty(event)
  % sort the events on the sample on which they occur
  [dum, indx] = sort([event.sample]);
  event = event(indx);
else
  warning(sprintf('no events found in %s', filename));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SUBFUNCTION that cuts a string into pieces on a user-defined separator character
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [tok] = tokenize(str, sep)
tok = {};
f = find(str==sep);
f = [0, f, length(str)+1];
for i=1:(length(f)-1)
  tok{i} = str((f(i)+1):(f(i+1)-1));
end
