Merge pull request fieldtrip#661 from robertoostenveld/thresholdpeak

ENH - also determine the offset (i.e. 3rd column) corresponding to th…

ft_artifact_threshold.m

@@ -177,7 +177,7 @@
 numtrl      = size(cfg.trl,1);
 channel     = ft_channelselection(artfctdef.channel, hdr.label);
 channelindx = match_str(hdr.label,channel);
-artifact    = zeros(0,2);
+artifact    = zeros(0,3);
 
 if ~isempty(artfctdef.onset) || ~isempty(artfctdef.offset)
   if artfctdef.onset>0 && artfctdef.offset>0
@@ -248,9 +248,30 @@
   if any(artval)
     begsample = find(diff([false artval])>0) + cfg.trl(trlop,1) - 1;
     endsample = find(diff([artval false])<0) + cfg.trl(trlop,1) - 1;
-    artifact  = cat(1, artifact, [begsample(:) endsample(:)]);
+    offset    = nan(size(begsample));
+
+    if size(dat,1)==1
+      % determine the offset of the peak value, this only works in case of a single channel
+      for i=1:numel(begsample)
+        seg = dat(begsample(i):endsample(i)); % get the segment of data
+        if all(seg>=artfctdef.max) || strcmp(direction, 'up')
+          [~, indx] = max(seg);
+          offset(i) = 1 - indx; % relative to the start of the segment, 0 is the first sample, -1 is the 2nd, etc.
+        elseif all(seg<=artfctdef.min) || strcmp(direction, 'down')
+          [~, indx] = min(seg);
+          offset(i) = 1 - indx; % relative to the start of the segment, 0 is the first sample, -1 is the 2nd, etc.
+        end % if up or down
+      end % for each artifact in this trial
+    end % if single channel
+
+    artifact  = cat(1, artifact, [begsample(:) endsample(:) offset(:)]);
   end
 
+end % for trllop
+
+if any(isnan(artifact(:,3)))
+  % don't keep the offset if it cannot be determined consistently
+  artifact = artifact(:,[1 2]);
 end
 
 ft_info('detected %d artifacts\n', size(artifact,1));
@@ -264,4 +285,3 @@
 % do the general cleanup and bookkeeping at the end of the function
 ft_postamble provenance
 ft_postamble previous data
-

test/test_ft_artifact_threshold.m

@@ -3,6 +3,8 @@
 % WALLTIME 00:10:00
 % MEM 1gb
 
+% this script serves to test the asymetric onset/offset and the peak-detection
+
 %%
 
 data = [];
@@ -23,6 +25,20 @@
 cfg.artfctdef.threshold.max = 400;
 [cfg1, artifact1] = ft_artifact_threshold(cfg, data);
 
+cfg = [];
+cfg.trl = artifact1;
+data1 = ft_redefinetrial(cfg, data);
+
+cfg = [];
+timelock1 = ft_timelockanalysis(cfg, data1);
+
+figure
+plot(timelock1.time, timelock1.avg, '.-');
+
+% the peak should be at t=0
+[val, indx] = max(timelock1.avg);
+assert(isalmostequal(timelock1.time(indx), 0));
+
 %%
 
 cfg = [];
@@ -43,6 +59,21 @@
 
 assert(~isequal(artifact3, artifact1));
 
+cfg = [];
+cfg.trl = artifact3;
+data3 = ft_redefinetrial(cfg, data);
+
+cfg = [];
+timelock3 = ft_timelockanalysis(cfg, data3);
+
+figure
+plot(timelock3.time, timelock3.avg, '.-'); % it should be asymetric
+
+% the peak should be at t=0
+[val, indx] = max(timelock3.avg);
+assert(isalmostequal(timelock3.time(indx), 0, 'abstol', 100*eps));
+
+
 %%
 
 cfg = [];
@@ -51,7 +82,7 @@
 cfg.artfctdef.threshold.offset = 1;
 [cfg4, artifact4] = ft_artifact_threshold(cfg, data);
 
-assert(isequal(artifact4, [1 1000]));
+assert(isequal(artifact4, [1 1000 -499]));
 
 %%
 % invert the data and invert the thresholds
@@ -66,4 +97,3 @@
 
 assert(isequal(artifact5, artifact3));
 
-