% Demo for plotting filtered and unfiltered fixation data 
infile='EggertStraube_FixationDriftData.mat';
load(infile,'results');
%% plot the trial shown in the Supporting Information S2 text Fig. A
plot_eyeposition(results);

%% plot all trials of a given task, block, and subject
SID=50;
task_i=1;
block_i=1;

% get the number of trials for this task, block, subject
SIDall=[results(task_i,block_i,:).SubjectIdentifier];
subj_i=find(SIDall==SID,1,'first');
if isempty(subj_i)
   error('SID=%d not found.');
end
NValidTrials=length(results(task_i,block_i,subj_i).VersVerg_deg);



for k=1:NValidTrials
   figure(k);
   plot_eyeposition(results,task_i,block_i,k,SID,k);
end


function plot_eyeposition(results,task_i,block_i,trial_number,SID,fignr)
% plots left and right eye position during valid fixation intervals
%
% Arguments:
% results     : structured array stored in EggertStraube_FixationDriftData.mat
%                   load data before calling : load('EggertStraube_FixationDriftData','results')
% task_i      : scalar integer specifying the task type:  1:stereo task, 2: acuity task (default: 1)
% block_i     : scalar integer specifying the measurement block: % 1: block A; 2: block B (default: 1)
% trial_number: scalar integer specifying the trial number (default: 6)
% SID         : integer vector containig the subject identifiers for which the plot is generated
%                  SID=[] data of all subjects are be plotted  (default: 50)
% fignr       : integer specifying the figure number: 0: plot in current figure; >0: plot in this figure(fignr)
%
% calling plot_eyeposition(results) creates the plot shown in Supporting Information S2 text (Fig. A)
if nargin<2 || isempty(task_i)
   task_i=1;   % 1:stereo, 2: acuity
end
if nargin<3 || isempty(block_i)
   block_i=1;  % 1: block A; 2: block B
end
if nargin<4 || isempty(trial_number)
   trial_number=6; % trial number
end
if nargin<5 || isempty(SID)
   SID=50;     % subject identifier
end

if nargin<5 || isempty(fignr)
   fignr=1;     % figure number
end

   dim=size(results);  % [NuberofTasks, NumberOfBlocks, NumberOfSubjects]

   SIDall=[results(1,1,:).SubjectIdentifier];
   if ~isempty(SID)
      % find the subject index from the subject identifier
      subj_i=NaN(size(SID));
      for j=1:numel(SID)
         ind=find(SIDall==SID(j),1,'first');
         if ~isempty(ind)
            subj_i(j)=ind;
         else
            fprintf('SID=%d not found.\n',SID(j));
         end
      end
      subj_i=subj_i(~isnan(subj_i));
      SIDall=SIDall(subj_i);
      results=results(:,:,subj_i);
   end   

   for subj_i=1:length(SIDall)
      NValidTrials=length(results(task_i,block_i,subj_i).VersVerg_deg);
      fprintf('number of valid trials: %d\n',NValidTrials);
      
      SamplingRate_Hz=results(task_i,block_i,subj_i).SamplingRate_Hz;
      
      
      
      
      
      
      if trial_number>NValidTrials
         fprintf('SID=%3d: specified trial number %d replaced by maximum trial number %d:\n',SIDall(subj_i),trial_number,NValidTrials);
         vtrialn=NValidTrials;
      else
         vtrialn=trial_number;
      end
      VersVerg=results(task_i,block_i,subj_i).VersVerg_deg{vtrialn};
      VersVergSlow=results(task_i,block_i,subj_i).VersVergSlow_deg{vtrialn};
      Segmentation=results(task_i,block_i,subj_i).Segmentation{vtrialn};
      PupilDiameter_mm=results(task_i,block_i,subj_i).PupilDiameter_mm_trace{vtrialn};

      SegmentClass=zeros(size(VersVergSlow,1),1);
      for fast_class=1:2
         sind=find([Segmentation(:).class]==fast_class);
         for sacc_i=1:length(sind)
            fast_i=floor((Segmentation(sind(sacc_i)).t_start_s+[0,Segmentation(sind(sacc_i)).dur_s])*SamplingRate_Hz)+1;
            fast_i=fast_i+[-1,1]*floor(0.001*SamplingRate_Hz);
            fast_i(2)=min(fast_i(2),length(SegmentClass));
            fast_i(1)=max(1,fast_i(1));
            SegmentClass(fast_i(1):fast_i(2))=fast_class;
         end
      end

      Elraw=VersVerg(:,1)+VersVerg(:,2)/2;  % raw version and vergence in left eye centred coordinates  (deg)
      Erraw=VersVerg(:,1)-VersVerg(:,2)/2;  % raw version and vergence in right eye centred coordinates  (deg)

      El=VersVergSlow(:,1)+VersVergSlow(:,2)/2;  % cleaned left eye centred coordinates
      Er=VersVergSlow(:,1)-VersVergSlow(:,2)/2;  % cleaned right eye centred coordinates
      
      t=(0:length(El)-1)'/SamplingRate_Hz;
      
      axh=zeros(3,1);
      ph=zeros(2,1);
      if fignr>0
         figure(fignr);
      end
      clf
      set(gcf,'Position',[53    64   580   794]);
      axh(1)=subplot(3,1,1);
      hold on
      set(gca,'Fontsize',12,'Linewidth',1.5);
      title(sprintf('SID=%d: trial #= %d',SIDall(subj_i),vtrialn));
      
      tmp=Elraw;
      tmp(SegmentClass==1)=NaN;
      ph(1)=plot(t,tmp,'-m','Linewidth',1.5);
      tmp=Elraw;
      tmp(SegmentClass==0)=NaN;
      plot(t,tmp,'--m','Linewidth',1.5);
      
      tmp=El;
      tmp(SegmentClass==1)=NaN;
      ylims=[1,2]+nanmean(tmp)-1.5;
      ph(2)=plot(t,tmp,'-b','Linewidth',1.5);
      tmp=El;
      tmp(SegmentClass==0)=NaN;
      plot(t,tmp,'--b','Linewidth',1.5);
      set(gca,'YLim',ylims);
      set(get(gca,'YLabel'),'string','left eye position (deg)');
      
      
      
      axh(1)=subplot(3,1,2);
      hold on
      set(gca,'Fontsize',12,'Linewidth',1.5);
      
      tmp=Erraw;
      tmp(SegmentClass==1)=NaN;
      ph(1)=plot(t,tmp,'-m','Linewidth',1.5);
      tmp=Erraw;
      tmp(SegmentClass==0)=NaN;
      plot(t,tmp,'--m','Linewidth',1.5);
      
      tmp=Er;
      tmp(SegmentClass==1)=NaN;
      ylims=[-2,-1]+nanmean(tmp)+1.5;
      ph(2)=plot(t,tmp,'-b','Linewidth',1.5);
      tmp=Er;
      tmp(SegmentClass==0)=NaN;
      plot(t,tmp,'--b','Linewidth',1.5);
      set(gca,'YLim',ylims);
      
      set(get(gca,'YLabel'),'string','right eye position (deg)');
      
      
      legend(ph,'non-filtered','spline filter');
      axh(2)=subplot(3,1,3);
      hold on
      set(gca,'Fontsize',12,'Linewidth',1.5);
      plot(t,PupilDiameter_mm,'-k','Linewidth',1.5);
      set(get(gca,'XLabel'),'string','time (s)');
      set(get(gca,'YLabel'),'string','pupil diameter (mm)');
      
      
      figure(gcf);
      if subj_i<length(SIDall)
         hstr=input('Press ENTER to continue!','s');
      end
   end
end