%%

cfg = [];
cfg.method = 'channel';
data_clean = ft_rejectvisual(cfg, data);

disp(data_clean.cfg.artfctdef.channel.artifact)


%%

cfg = [];
cfg.channel = 'EOG';
data_eog = ft_selectdata(cfg, data);

cfg = [];
cfg.method = 'channel';
data_eog_clean = ft_rejectvisual(cfg, data_eog);

disp(data_eog_clean.cfg.artfctdef.channel.artifact)

% remove the marked segments from the complete MEG+EEG data structure
cfg = [];
cfg.artfctdef = data_eog_clean.cfg.artfctdef; % use the artifacts that we just identified
data_all_clean = ft_rejectartifact(cfg, data);

%%

cfg = [];
cfg.artfctdef.eog.channel = 'EOG';
cfg.artfctdef.eog.interactive = 'yes';
cfg = ft_artifact_eog(cfg, data);

disp(cfg.artfctdef.eog.artifact)

% remember it for later inspection
cfg_automatic_eog.artfctdef = cfg.artfctdef;

% remove the marked segments
cfg.artfctdef.reject = 'complete';
data_clean = ft_rejectartifact(cfg, data);

%%

cfg = [];
cfg.artfctdef.eog.channel = 'EEG';  % use the EEG channels to detect eye artifacts
cfg.artfctdef.eog.interactive = 'yes';
cfg.artfctdef.eog.cutoff = 20;
cfg = ft_artifact_eog(cfg, data);

%%

cfg = [];
cfg.artfctdef.threshold.channel = 'EOG';
cfg.artfctdef.threshold.demean = 'yes';
cfg.artfctdef.threshold.min = -100 * 1e-6;
cfg.artfctdef.threshold.max =  100 * 1e-6;
cfg.artfctdef.threshold.bpfilter = 'no';
cfg = ft_artifact_threshold(cfg, data);

disp(cfg.artfctdef.threshold.artifact)

% remember it for later inspection
cfg_automatic_threshold.artfctdef = cfg.artfctdef;

% remove the marked segments
cfg.artfctdef.reject = 'complete';
data_clean = ft_rejectartifact(cfg, data);


%%

cfg = [];
cfg.channel = 'EOG';
cfg = ft_databrowser(cfg, data);

disp(cfg.artfctdef.visual.artifact)

% remove the marked segments
data_clean = ft_rejectartifact(cfg, data);

%%

cfg = [];
cfg.channel = 'EOG';
cfg.continuous = 'yes';
cfg.blocksize = 30;
cfg.ylim = [-200 200]*1e-6;

% to get you started here are some that I marked 
cfg.artfctdef.visual.artifact = [
       52484       52546
       63142       63542
       67315       67479
       69299       69409
       70303       70388
       70877       70993
       71821       71961
       73430       73540
       75916       76014
       80548       80682
       ];

cfg = ft_databrowser(cfg, data);

% remove the marked segments
data_clean = ft_rejectartifact(cfg, data);

%%


cfg = [];
cfg.channel = 'EOG';
cfg.continuous = 'yes';
cfg.blocksize = 30;
cfg.ylim = [-200 200]*1e-6;

% reuse the automatic artifacts from above
cfg.artfctdef.eog       = cfg_automatic_eog.artfctdef.eog;
cfg.artfctdef.threshold = cfg_automatic_threshold.artfctdef.threshold; 

% inspect and adjust the marked segments
cfg = ft_databrowser(cfg, data);

% remove the marked segments
data_clean = ft_rejectartifact(cfg, data);

%%

cfg = [];
cfg.method = 'summary';
cfg.channel = 'EOG';
cfg.keepchannel = 'yes';
data_clean = ft_rejectvisual(cfg, data);

cfg = [];
cfg.method = 'summary';
cfg.channel = 'EEG';
cfg.keepchannel = 'yes';
data_clean = ft_rejectvisual(cfg, data);

%%

cfg = [];
cfg.channel = 'eog';
data_eog = ft_selectdata(cfg, data);

cfg = [];
cfg.method = 'summary';
cfg.keepchannel = 'yes';
data_eog_clean = ft_rejectvisual(cfg, data_eog);

cfg = [];
cfg.channel = 'EEG*';  % 'eeg' as such is not recognized
data_eeg = ft_selectdata(cfg, data);

cfg = [];
cfg.method = 'summary';
cfg.keepchannel = 'yes';
data_eeg_clean = ft_rejectvisual(cfg, data_eeg);

cfg = [];
cfg.channel = 'megmag';
data_megmag = ft_selectdata(cfg, data);

cfg = [];
cfg.method = 'summary';
cfg.keepchannel = 'yes';
data_megmag_clean = ft_rejectvisual(cfg, data_megmag);

cfg = [];
cfg.method = 'summary';
cfg.channel = 'megplanar';
data_megplanar = ft_selectdata(cfg, data);

cfg = [];
cfg.method = 'summary';
cfg.keepchannel = 'yes';
data_megplanar_clean = ft_rejectvisual(cfg, data_megplanar);

cfg = [];
cfg.artfctdef.eog.artifact       = data_eog_clean.cfg.artfctdef.summary.artifact;
cfg.artfctdef.eeg.artifact       = data_eeg_clean.cfg.artfctdef.summary.artifact;
cfg.artfctdef.megmag.artifact    = data_megmag_clean.cfg.artfctdef.summary.artifact;
cfg.artfctdef.megplanar.artifact = data_megplanar_clean.cfg.artfctdef.summary.artifact;

% inspect and adjust the marked segments
cfg.channel = 'megmag'; % look only at the 102 magnetometers
cfg = ft_databrowser(cfg, data);

% remove the marked segments
cfg = rmfield(cfg, 'channel'); % this is not an option for ft_rejectartifact
data_clean = ft_rejectartifact(cfg, data);

%%

% this is what we can see in the BIDS events.tsv file
Famous      = [5 6 7];
Unfamiliar  = [13 14 15];
Scrambled   = [17 18 19];

sum(ismember(data_clean.trialinfo, Famous))
sum(ismember(data_clean.trialinfo, Unfamiliar))
sum(ismember(data_clean.trialinfo, Scrambled))

%%

nsamples = sum(cellfun(@length, data.trial));
artifact = zeros(1, nsamples);
for i=1:size(cfg_automatic_threshold.artfctdef.threshold.artifact, 1)
  begsample = cfg_automatic_threshold.artfctdef.threshold.artifact(i,1);
  endsample = cfg_automatic_threshold.artfctdef.threshold.artifact(i,1);
  artifact(begsample:endsample) = 1;
end

% compute the time in the experiment, this does not account for the time between trials or runs
time = ((1:nsamples)-1)/data.fsample;

% plot where in time the artifacts happen
plot(time, artifact)

% make a smooth representation of where in time the artifacts happen
% this smoothes with 3000 samples, which is 10 seconds
plot(time, conv(artifact, ones(1,3000), 'same'))

%%

cfg = [];
cfg.channel = 'VEOG'; % a single channel is enough
data_eog = ft_selectdata(cfg, data);

cfg = [];
cfg.artfctdef = cfg_automatic_threshold.artfctdef;
cfg.artfctdef.reject = 'nan';
data_artifact = ft_rejectartifact(cfg, data_eog);
for i=1:numel(data_artifact.trial)
  % find where the nans are in every trial
  data_artifact.trial{i} = isnan(data_artifact.trial{i});
end

cfg = [];
avg_artifact = ft_timelockanalysis(cfg, data_artifact);

plot(avg_artifact.time, avg_artifact.avg);









