TRACKPixx Demo 7 – Running a simple saccade to target task with data visualization

Note

This demo requires TRACKPixx Revision 18 or later. You can check for recent firmware updates at vpixx.com/whatsnew

This demo is a simple eight-trial saccade to target experiment. During a trial, a fixation cross is presented in the center of the display, and an eye-tracking schedule is started. We use the IsSubjectFixating flags and GetEyePosition to check whether the participant is fixating on the cross. When a central fixation is detected, the cross is replaced with a target in one of four locations. Participants must saccade to the target as quickly as possible. We again use the fixation flags and eye position to determine whether the target was successfully fixated. On a successful fixation, the time is logged and a 500 ms visual mask appears before the next trial. During this time, gaze data from the trial is imported from the I/O controller to the main computer and saved. After the experiment ends we plot saccade duration, xy position, velocity, acceleration, and pupil size for each trial.

At the beginning of the demo there is an optional calibration step. This step calls TPxTrackpixx3CalibrationTesting, which implements a standard MATLAB TRACKPixx3 calibration script. The tracker must be calibrated for every new participant. We also recommend calibrating after a participant moves away from the chinrest.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
function TPxSaccadeToTarget(initRequired)
%
% This demo flashes a sequence of 8 targets and tracks saccades to their
% position. 

% On a given trial, a fixation cross appears. When participants are
% fixating on the cross, a target appears in one of four locations (above,
% below, left, or right). Participants must saccade to the target as
% quickly as possible. When the saccade is complete, a 500ms mask is
% presented, and the next trial begins.

% Targets are circles with high or low luminance (high or low). Each luminance x
% position combination is repeated once, leading to 8 trials. After
% data is collected, data is plotted and output to a .csv file with headers
% for futher analysis 
%
% If initRequired is set to 1, the function first calls
% TPxTrackpixx3CalibrationTesting to connect to the TRACKPixx3 and
% calibrate the tracker. If you are calibrating your own setup, please
% check the calibration function to ensure your settings are correct.

% Made and tested with:
% -- TRACKPixx3 firmware revision 18 
% -- DATAPixx3 firmware revision 19 
% -- MATLAB version 9.6.0.1150989 (R2019a) 
% -- Psychtoolbox verison 3.0.15 
% -- Datapixx Toolbox version 3.7.5735
% -- Windows 10 version 1903, 64bit

% Sep 19, 2019  lef     Written
% Oct 17, 2019  lef     Revised
% Feb 24, 2020  lef     Revised, removes conversion NaN, now automated

%% Step 1 - Initialize (if needed)
if nargin==0
    initRequired=1;
end

% Get some user input
fileName= input('Enter participant name: ', 's');
fileID = [fileName '.mat'];

displayWidth = input('Enter width, in cm, of display monitor: ', 's');
displayWidth = str2double(displayWidth);

displayDistance = input('Enter distance, in cm, of participant to display monitor: ', 's');
displayDistance = str2double(displayDistance);

%If a calibration is needed, call the calibration script
if initRequired
    fprintf('\nInitialization required\n\nCalibrating the device...');
    TPxTrackpixx3CalibrationTesting;
end

%Connect to TRACKPixx3
Datapixx('Open');
Datapixx('SetTPxAwake');
Datapixx('SetupTPxSchedule');
Datapixx('RegWrRd');

%Set some flag settings. These values are the default settings but we
%specify them anyway.
fixationThreshold = 2500;               %pixels/s
framesForFixationFlag = 25;             %how many consecutive frames have to be below threshold for fixation flag to raise
saccadeThreshold = 10000;               %pixels/s
framesForSaccadeFlag = 10;              %how many consecutive frames have to be above threshold for saccade flag to raise

Datapixx('SetFixationThresholds' , fixationThreshold, framesForFixationFlag);
Datapixx('SetSaccadeThresholds' , saccadeThreshold, framesForSaccadeFlag);
Datapixx('RegWrRd');

%% Step 2 - Saccade to target task

%Set experiment parameters
maskPresentationTime = 0.5;
fixationPresentationTime = 0.5;
margin = 35;                                                % +/- margin of error, in pixels, within which the eyes can be considered in position

%open window
Screen('Preference', 'SkipSyncTests', 1 );
screenID = 2;                                              %change to switch display
[windowPtr, rect]=Screen('OpenWindow', screenID, [0,0,0]);
pixelSize = displayWidth/rect(4);
Screen('BlendFunction', windowPtr, 'GL_SRC_ALPHA', 'GL_ONE_MINUS_SRC_ALPHA');
Screen('Flip',windowPtr);

%describe targets and target locations
dotRadius = 20;
luminances = [1, 25,25,25;...                              %dark
              2, 250, 250, 250];                           %light
radius = 300;
center = rect(3:4)/2;
locations = [1, center(1), center(2)-radius;...            %top (in screen coordinates) 
             2, center(1)+radius, center(2);...            %right
             3, center(1), center(2)+radius;...            %bottom
             4, center(1)-radius, center(2)];              %left
         
%generate a stimuli list with target characteristics for each trial, and
%shuffle the order
numTargets = size(luminances(:,1),1);
numLocations = size(locations(:,1),1);
numReps = 1;
numTrials = numTargets*numLocations*numReps;
stimuli = [];

for s = 1:numTargets
    for k = 1:numLocations
        for m = 1:numReps
            stimuli = [stimuli; luminances(s,:), locations(k, :)]; 
        end
    end
end

stimuli=stimuli(randperm(size(stimuli,1)),:);           %shuffle

%create a structure to hold our trial data
trials = struct('Trial', [],...
                'Target', [],...
                'TargetLocation', [],...
                'TrialStart', [],...
                'TargetOnset', [],...
                'SaccadeOnset', [],...
                'SaccadeEnd', [],...
                'SaccadeReactionTime', [],...
                'SaccadeDuration', [],...
                'EyeData', []);


%show instructions to participant
Screen('DrawLine', windowPtr, [255,255,255], center(1) + 8, center(2), center(1)-8, center(2), 2);
Screen('DrawLine', windowPtr, [255,255,255], center(1), center(2)+8, center(1), center(2)-8, 2);
text_to_draw = ['SACCADE TO TARGET EXPERIMENT:\n\nStare at the cross in the middle of the screen.'...
                '\nAs soon as you see a dot appear, move your eyes to look at it as quickly as possible!\n\nPress any key to start.'];
DrawFormattedText(windowPtr, text_to_draw, 'center', 700, 255);
Screen('Flip', windowPtr);

%wait for participant to continue
[~, ~, ~] = KbPressWait;
Screen('Flip', windowPtr);
WaitSecs(1);

start_time = Datapixx('GetTime');

for k=1:numTrials
    
    trials(k).Trial = k;
    trials(k).Target = stimuli(k,1:4);
    trials(k).TargetLocation = stimuli(k,5:7);

    %draw fixation cross
    Screen('DrawLine', windowPtr, [255,255,255], center(1) + 8, center(2), center(1)-8, center(2), 2);
    Screen('DrawLine', windowPtr, [255,255,255], center(1), center(2)+8, center(1), center(2)-8, 2);
    Screen('Flip', windowPtr);
    trials(k).TrialStart = Datapixx('GetTime');
    WaitSecs(fixationPresentationTime); 
    
    %start logging eye data
    Datapixx('StartTPxSchedule');
    Datapixx('RegWrRd');
    
    %wait until subject is fixating on fixation cross
    while 1
        
        %get new eye data
        Datapixx('RegWrRd');
 
        %check for a fixation
        [lFlag, rFlag] = Datapixx('IsSubjectFixating');
        if lFlag && rFlag
            %if fixating, get eye position
            [xScreenRight, yScreenRight, xScreenLeft, yScreenLeft, ~, ~, ~, ~, ~] = Datapixx('GetEyePosition');
            
            %convert screen coordinates in VPixx coordinates (origin in middle of
            %screen) to Psychtoolbox Screen coordinates (origin in top left)
            fixationLocs = Datapixx('ConvertCoordSysToCustom', [xScreenRight, yScreenRight; xScreenLeft, yScreenLeft]);
            
            %confirm fixation is on cross in Psychtoolbox Screen coordinates 
            if inpolygon(fixationLocs(1,1), fixationLocs(1,2), [center(1)-margin, center(1)+margin] , [center(2)-margin, center(2)+margin])...
               && inpolygon(fixationLocs(2,1), fixationLocs(2,2), [center(1)-margin, center(1)+margin] , [center(2)-margin, center(2)+margin])                
                break;
            end
        end
    end
    
    %draw target
    colour = trials(k).Target(2:4);
    x = trials(k).TargetLocation(2);
    y = trials(k).TargetLocation(3);
    Screen('FillOval', windowPtr, colour, [x-dotRadius, y-dotRadius, x+dotRadius, y+dotRadius]);
    
    %set a marker to get the exact time the screen flips
    Datapixx('SetMarker');
    Datapixx('RegWrVideoSync');
    Screen('Flip', windowPtr);
    
    %collect marker data
    Datapixx('RegWrRd');
    trials(k).TargetOnset = Datapixx('GetMarker');
    
    %Now that the target has been shown, start a while loop to repeatedly
    %check to see if participant is fixating on the target
    while 1
        
        %get new eye data
        Datapixx('RegWrRd');

        %check for a fixation
        [lFlag, rFlag] = Datapixx('IsSubjectFixating');
        if lFlag && rFlag
            [xScreenRight, yScreenRight, xScreenLeft, yScreenLeft, ~, ~, ~, ~, ~] = Datapixx('GetEyePosition');
            
            %convert screen coordinates in VPixx coordinates (origin in middle of
            %screen) to Psychtoolbox Screen coordinates (origin in top left)
            fixationLocs = Datapixx('ConvertCoordSysToCustom', [xScreenRight, yScreenRight; xScreenLeft, yScreenLeft]);
            
            %confirm fixation is on target in Psychtoolbox Screen coordinates 
            if inpolygon(fixationLocs(1,1), fixationLocs(1,2), [x-margin, x+margin] , [y-margin, y+margin])...
               && inpolygon(fixationLocs(2,1), fixationLocs(2,2),[x-margin, x+margin] , [y-margin, y+margin])
                
                %Stop recording
                Datapixx('StopTPxSchedule');
                Datapixx('RegWrRd');
                break;
            end
        end
    end
    
    %create and show visual mask - mask is a full-screen display of 5 pixel wide squares
    %with a random grayscale value.
    squaresize = 5; 
    numsquares=[rect(4)/squaresize, rect(3)/squaresize];
    pattern = rand([numsquares, 1])*128;                    
    mask = imresize(pattern, squaresize, 'nearest');
    textureIndex=Screen('MakeTexture', windowPtr, mask);
    Screen('DrawTexture', windowPtr, textureIndex); 
    Screen('Flip', windowPtr);
    WaitSecs(maskPresentationTime);
    
    %read in eye data
    Datapixx('RegWrRd');
    status = Datapixx('GetTPxStatus');
    toRead = status.newBufferFrames;
    [bufferData, ~, ~] = Datapixx('ReadTPxData', toRead);
    
    %bufferData is formatted as follows:
    %1      --- Timetag (in seconds)
    %2      --- Left Eye X (in pixels)
    %3      --- Left Eye Y (in pixels)
    %4      --- Left Pupil Diameter (in pixels)
    %5      --- Right Eye X (in pixels)
    %6      --- Right Eye Y (in pixels)
    %7      --- Right Pupil Diameter (in pixels)
    %8      --- Digital Input Values (24 bits)
    %9      --- Left Blink Detection (0=no, 1=yes)
    %10     --- Right Blink Detection (0=no, 1=yes) 
    %11     --- Digital Output Values (24 bits)
    %12     --- Left Eye Fixation Flag (0=no, 1=yes) 
    %13     --- Right Eye Fixation Flag (0=no, 1=yes)  
    %14     --- Left Eye Saccade Flag (0=no, 1=yes) 
    %15     --- Right Eye Saccade Flag (0=no, 1=yes)  
    %16     --- Message code (integer) 
    %17     --- Left Eye Raw X (in pixels) 
    %18     --- Left Eye Raw Y (in pixels)  
    %19     --- Right Eye Raw X (in pixels)  
    %20     --- Right Eye Raw Y (in pixels)
    
    %IMPORTANT: "RIGHT" and "LEFT" refer to the right and left eyes shown
    %in the console overlay. In tabletop and MEG setups, this view is
    %inverted. This means "RIGHT" in our labelling convention corresponds
    %to the participant's left eye. Similarly "LEFT" in our convention
    %refers to left on the screen, which corresponds to the participant's
    %right eye.

    %If you are using an MRI setup with an inverting mirror, "RIGHT" will
    %correspond to the participant's right eye.
    
    %save eye data from trial as a table in the trial structure
    trials(k).EyeData = array2table(bufferData, 'VariableNames', {'TimeTag', 'LeftEyeX', 'LeftEyeY', 'LeftPupilDiameter', 'RightEyeX', 'RightEyeY', 'RightPupilDiameter',...
                                    'DigitalIn', 'LeftBlink', 'RightBlink', 'DigitalOut', 'LeftEyeFixationFlag', 'RightEyeFixationFlag', 'LeftEyeSaccadeFlag', 'RightEyeSaccadeFlag',...
                                    'MessageCode', 'LeftEyeRawX', 'LeftEyeRawY', 'RightEyeRawX', 'RightEyeRawY'});
    %interim save
    save(fileID, 'trials');
end

%Finish presentation
Screen('Closeall');
Datapixx('StopTPxSchedule');
Datapixx('RegWrRd');
finish_time = Datapixx('GetTime');


%% Step 3 - Read and process data
fprintf('\nRecording lasted %f seconds', finish_time-start_time);
fprintf('\nProcessing... this will take a while...');

%Data processing steps
% 1  - Adjust eye movement flags. The flags require calculating velocity
% over several frames and are therefore offset in a deterministic way, and
% need to be shifted to align with true event onset
% 2 - Get some saccade details - onset, reaction time, duration, end
% 3 - Convert x, y positions into degrees of visual angle
% 4 - Calculate Velocity (and filter) in degrees of visual angle
% 5 - Calculate Acceleration (and filter) in degrees of visual angle
% 6 - Create a "proportionOfTrial" variable for plotting purposes


% 1 - Loop throught and create adjusted eye movement flags that are
% correctly aligned with movement onset/end. This has three discrete steps:

% a) Both types of eye movement flags are based on the 4 previous frames,
% and the 4 subsequent frames. As a result, velocity for frame j is not
% calculated until frame j+5. All flags must be moved back 5 frames to
% account for this.

% b) The arguments framesForFixationFlag and framesForSaccadeFlag specify
% how many consecutive frames a velocity has to be above/below threshold
% for an eye movement flag to raise. This reduces noise but creates an
% offset of the flag onset, corresponding to the number of frames
% specified. We correct for this by adding flags to this number of frames
% prior to flag onset.

% c) For a saccade flag to drop, the previous 2 frames must be
% below/above the threshold. We can remove these last 2 flags to get true
% movement end.
for k = 1:numel(trials) 
    trials(k).EyeData = addvars(trials(k).EyeData, nan(height(trials(k).EyeData), 1),...
                                            nan(height(trials(k).EyeData), 1),...
                                            nan(height(trials(k).EyeData), 1),...
                                            nan(height(trials(k).EyeData), 1),...
                                            'NewVariableNames', {'LeftEyeFixationAdjusted', 'RightEyeFixationAdjusted', 'LeftEyeSaccadeAdjusted', 'RightEyeSaccadeAdjusted'});
    
    % a) correct for velocity calculation offset                                    
    for s = 1: height(trials(k).EyeData)
        if height(trials(k).EyeData)-s >=5
            if trials(k).EyeData.LeftEyeFixationFlag(s+5) ==1
                trials(k).EyeData.LeftEyeFixationAdjusted(s)=trials(k).EyeData.LeftEyeFixationFlag(s+5);
            end
            if trials(k).EyeData.RightEyeFixationFlag(s+5) ==1
                trials(k).EyeData.RightEyeFixationAdjusted(s)=trials(k).EyeData.RightEyeFixationFlag(s+5);
            end
            if trials(k).EyeData.LeftEyeSaccadeFlag(s+5) ==1
                trials(k).EyeData.LeftEyeSaccadeAdjusted(s)=trials(k).EyeData.LeftEyeSaccadeFlag(s+5);
            end
            if trials(k).EyeData.RightEyeSaccadeFlag(s+5)
                trials(k).EyeData.RightEyeSaccadeAdjusted(s)=trials(k).EyeData.RightEyeSaccadeFlag(s+5);
            end
        end
    end
    
    %b) correct for flag onset offset
    for s = 1: height(trials(k).EyeData)
        if height(trials(k).EyeData)-s >= framesForFixationFlag
            %check for an onset a specficied number of frames ahead
            onset = s+framesForFixationFlag;
            if isnan(trials(k).EyeData.LeftEyeFixationAdjusted(onset-1)) && ~isnan(trials(k).EyeData.LeftEyeFixationAdjusted(onset))
                trials(k).EyeData.LeftEyeFixationAdjusted(s:onset)=1;
            end
            
            if isnan(trials(k).EyeData.RightEyeFixationAdjusted(onset-1)) && ~isnan(trials(k).EyeData.RightEyeFixationAdjusted(onset))
                trials(k).EyeData.RightEyeFixationAdjusted(s:onset)=1;
            end
        end
        
        if height(trials(k).EyeData)-s >= framesForSaccadeFlag
            %check for an onset a specficied number of frames ahead
            onset = s+framesForSaccadeFlag;
            if isnan(trials(k).EyeData.LeftEyeSaccadeAdjusted(onset-1)) && ~isnan(trials(k).EyeData.LeftEyeSaccadeAdjusted(onset))
                trials(k).EyeData.LeftEyeSaccadeAdjusted(s:onset)=1;
            end
            
            if isnan(trials(k).EyeData.RightEyeSaccadeAdjusted(onset-1)) && ~isnan(trials(k).EyeData.RightEyeSaccadeAdjusted(onset))
                trials(k).EyeData.RightEyeSaccadeAdjusted(s:onset)=1;
            end
        end
    end 
    
    % c) correct for flag end offset
    for s = 4: height(trials(k).EyeData)
        %check for a flag drop and clear last 2 flags
        if isnan(trials(k).EyeData.LeftEyeFixationAdjusted(s)) && ~isnan(trials(k).EyeData.LeftEyeFixationAdjusted(s-1))
            trials(k).EyeData.LeftEyeFixationAdjusted(s-2:s)=NaN;
        end
        
        if isnan(trials(k).EyeData.RightEyeFixationAdjusted(s)) && ~isnan(trials(k).EyeData.RightEyeFixationAdjusted(s-1))
            trials(k).EyeData.RightEyeFixationAdjusted(s-2:s)=NaN;
        end
        
        if isnan(trials(k).EyeData.LeftEyeSaccadeAdjusted(s)) && ~isnan(trials(k).EyeData.LeftEyeSaccadeAdjusted(s-1))
            trials(k).EyeData.LeftEyeSaccadeAdjusted(s-2:s)=NaN;
        end

        if isnan(trials(k).EyeData.RightEyeSaccadeAdjusted(s)) && ~isnan(trials(k).EyeData.RightEyeSaccadeAdjusted(s-1))
            trials(k).EyeData.RightEyeSaccadeAdjusted(s-2:s)=NaN;
        end
    end
end
    

% 2 - Get some saccade details - onset, reaction time, duration, end
for k = 1:numel(trials)
    saccadeStarted = NaN;
    saccadeEnded = NaN;
    
    for s = 1:height(trials(k).EyeData) 
       %If saccade flag is up for both eyes, log this as saccade start
       if isnan(saccadeStarted) && isnan(saccadeEnded) && ~isnan(trials(k).EyeData.LeftEyeSaccadeAdjusted(s)) && ~isnan(trials(k).EyeData.RightEyeSaccadeAdjusted(s))
           saccadeStarted = trials(k).EyeData.TimeTag(s);
       end

       %If saccade flag drops for both eyes, then saccade ended
       if ~isnan(saccadeStarted) && isnan(saccadeEnded) && isnan(trials(k).EyeData.LeftEyeSaccadeAdjusted(s)) && isnan(trials(k).EyeData.RightEyeSaccadeAdjusted(s))
           saccadeEnded = trials(k).EyeData.TimeTag(s);
       end
    end
    
    %get some saccade data 
    reactionTime = saccadeStarted - trials(k).TargetOnset;
    duration = saccadeEnded - saccadeStarted;
   
    trials(k).SaccadeOnset = saccadeStarted;
    trials(k).SaccadeEnd = saccadeEnded;
    trials(k).SaccadeReactionTime = reactionTime;
    trials(k).SaccadeDuration = duration;
end
       

% 3 - Convert x, y positions into degrees of visual angle
for k = 1:numel(trials)    
    trials(k).EyeData = addvars(trials(k).EyeData, nan(height(trials(k).EyeData), 1),...
                                        nan(height(trials(k).EyeData), 1),...
                                        nan(height(trials(k).EyeData), 1),...
                                        nan(height(trials(k).EyeData), 1),...
                                        'NewVariableNames', {'LeftEyeXDegrees', 'LeftEyeYDegrees', 'RightEyeXDegrees', 'RightEyeYDegrees'});

    for s = 1:height(trials(k).EyeData)
        xL = trials(k).EyeData.LeftEyeX(s) * pixelSize;
        trials(k).EyeData.LeftEyeXDegrees(s) = 2*atand(xL/(2*displayDistance));
        
        yL = trials(k).EyeData.LeftEyeY(s) * pixelSize;
        trials(k).EyeData.LeftEyeYDegrees(s) = 2*atand(yL/(2*displayDistance));

        xR = trials(k).EyeData.RightEyeX(s) * pixelSize;
        trials(k).EyeData.RightEyeXDegrees(s) = 2*atand(xR/(2*displayDistance));
        
        yR = trials(k).EyeData.RightEyeY(s) * pixelSize;
        trials(k).EyeData.RightEyeYDegrees(s) = 2*atand(yR/(2*displayDistance));
    end
end



% 4 - Calculate velocity
for k = 1:numel(trials)
    trials(k).EyeData = addvars(trials(k).EyeData, nan(height(trials(k).EyeData), 1),...
                                        nan(height(trials(k).EyeData), 1),...
                                        nan(height(trials(k).EyeData), 1),...
                                        nan(height(trials(k).EyeData), 1),...
                                        'NewVariableNames', {'LeftEyeVelocityRaw', 'RightEyeVelocityRaw', 'LeftEyeVelocity', 'RightEyeVelocity'});

    for s = 5:height(trials(k).EyeData)

        %Formula for velocity (as the tracker calulates it)
        %where i is the current frame 
        %deltaX = (sum(x(i+1):x(i+4)) - sum(x(i-4):x(i-1))) / 0.01;
        %deltaY = (sum(y(i+1):y(i+4)) - sum(y(i-4):y(i-1))) / 0.01;
        %velocity = sqrt(deltaX^2 + deltaY^2) 
        
        if height(trials(k).EyeData)-s >= 4
            
            leftX1 = sum(trials(k).EyeData.LeftEyeXDegrees(s-4:s-1));
            leftX2 = sum(trials(k).EyeData.LeftEyeXDegrees(s+1:s+4));
            
            leftY1 = sum(trials(k).EyeData.LeftEyeYDegrees(s-4:s-1));
            leftY2 = sum(trials(k).EyeData.LeftEyeYDegrees(s+1:s+4));
            
            trials(k).EyeData.LeftEyeVelocityRaw(s)= sqrt(((leftX2 - leftX1)/0.01)^2 + ((leftY2 - leftY1)/0.01)^2); 
            
            rightX1 = sum(trials(k).EyeData.RightEyeXDegrees(s-4:s-1));
            rightX2 = sum(trials(k).EyeData.RightEyeXDegrees(s+1:s+4));
            
            rightY1 = sum(trials(k).EyeData.RightEyeYDegrees(s-4:s-1));
            rightY2 = sum(trials(k).EyeData.RightEyeYDegrees(s+1:s+4));
            
            trials(k).EyeData.RightEyeVelocityRaw(s) = sqrt(((rightX2 - rightX1)/0.01)^2 + ((rightY2 - rightY1)/0.01)^2); 
        end

    end

    %Do some filtering on our data using a moving window average. Reverse
    %and forward filter to remove any phase-shift.
    windowSize = 25; 
    b = (1/windowSize)*ones(1,windowSize);
    a = 1;

    leftEye = trials(k).EyeData.LeftEyeVelocityRaw(:);
    rightEye = trials(k).EyeData.RightEyeVelocityRaw(:);

    %reverse filter
    leftEye = filter(b, a, flip(leftEye));
    rightEye = filter(b, a, flip(rightEye));

    %forward filter
    leftEye = filter(b, a, flip(leftEye));
    rightEye = filter(b, a, flip(rightEye));

    %assign to table
    trials(k).EyeData.LeftEyeVelocity = leftEye;
    trials(k).EyeData.RightEyeVelocity = rightEye;

end


%5 - Derive acceleration from velocity data. Acceleration is the change in
%velocity across a frame, which has a fixed duration of 1/2000th of a
%second
for k = 1:numel(trials)
    trials(k).EyeData = addvars(trials(k).EyeData, nan(height(trials(k).EyeData), 1),...
                                        nan(height(trials(k).EyeData), 1),...
                                        nan(height(trials(k).EyeData), 1),...
                                        nan(height(trials(k).EyeData), 1),...
                                        'NewVariableNames', {'LeftEyeAccelerationRaw', 'RightEyeAccelerationRaw', 'LeftEyeAcceleration', 'RightEyeAcceleration'});

    for s = 2:height(trials(k).EyeData)
        trials(k).EyeData.LeftEyeAccelerationRaw(s) = (trials(k).EyeData.LeftEyeVelocityRaw(s) - trials(k).EyeData.LeftEyeVelocityRaw(s-1))/(1/2000);
        trials(k).EyeData.RightEyeAccelerationRaw(s) = (trials(k).EyeData.RightEyeVelocityRaw(s) - trials(k).EyeData.RightEyeVelocityRaw(s-1))/(1/2000);
    end

    %Do some filtering on our data using a moving window average. Reverse
    %and forward filter to remove any phase-shift.
    windowSize = 25; 
    b = (1/windowSize)*ones(1,windowSize);
    a = 1;

    leftEye = trials(k).EyeData.LeftEyeAccelerationRaw(:);
    rightEye = trials(k).EyeData.RightEyeAccelerationRaw(:);

    %reverse filter
    leftEye = filter(b, a, flip(leftEye));
    rightEye = filter(b, a, flip(rightEye));

    %forward filter
    leftEye = filter(b, a, flip(leftEye));
    rightEye = filter(b, a, flip(rightEye));

    %assign to table
    trials(k).EyeData.LeftEyeAcceleration = leftEye;
    trials(k).EyeData.RightEyeAcceleration = rightEye;

end

% 6 Finally, for plotting purposes, create a column with proportion of trial time
for k = 1:numel(trials)
    trials(k).EyeData = addvars(trials(k).EyeData, nan(height(trials(k).EyeData), 1),...
                                        'NewVariableNames', 'ProportionTrial', 'After', 'TimeTag');

    trialStart = trials(k).EyeData.TimeTag(1);
    trialDuration = trials(k).EyeData.TimeTag(end)- trialStart;
    for s = 1:height(trials(k).EyeData)
        trials(k).EyeData.ProportionTrial(s) = (trials(k).EyeData.TimeTag(s)- trialStart)/trialDuration;
    end
end

save(fileID, 'trials');


%% Step 4 - Plot results

%Figure 1 - Saccade Reaction Time
%Figure 2 - Saccade Duration
%Figure 3 - X,Y screen position over time for target 1
%Figure 4 - X,Y screen position over time for target 2
%Figure 5 - Filtered velocity (degrees/second) for target 1
%Figure 6 - Filtered velocity (degrees/second) for target 2
%Figure 7 - Filtered acceleration (degrees/second^2) for target 1
%Figure 8 - Filtered acceleration (degrees/second^2) for target 2
%Figure 9 - Pupil dilation for target 1
%Figure 10 - Pupil dilation for target 2

lineColours = rand(numLocations, 3);

%Figure 1 - plot saccade reaction time
figure();
for k = 1:numTargets
    for s = 1:numLocations
        temp = [];
        for h= 1:numel(trials)
            if trials(h).Target(1) == k &&  trials(h).TargetLocation(1) == s
                temp = [temp; trials(h).SaccadeReactionTime];
            end
        end
        
    average = mean(temp);
    sterr = std(temp)/sqrt(numel(temp));
    subplot(1,2,k)
    hold on
    errorbar(s, average, sterr, 'Color', lineColours(s, :), 'MarkerFaceColor', lineColours(s, :), 'Marker', 'o');
    end
    
    xticks(1:numLocations);
    xticklabels({'Up', 'Right', 'Down', 'Left'});
    titleText=['Saccade Reaction times for target ' int2str(k)];
    title(titleText);
    ylabel('Reaction time (s)');
    xlabel('Saccade direction');
    xlim([0,numLocations+1]);
    ylim([0,1]);
end


%Figure 2 - plot saccade duration
figure();
for k = 1:numTargets
    for s = 1:numLocations
        temp = [];
        for h = 1:numel(trials)
            if trials(h).Target(1) == k &&  trials(h).TargetLocation(1) == s
                temp = [temp; trials(h).SaccadeDuration];
            end
        end
        
    average = mean(temp);
    sterr = std(temp)/sqrt(numel(temp));
    subplot(1,2,k)
    hold on
    errorbar(s, average, sterr, 'Color', lineColours(s, :), 'MarkerFaceColor', lineColours(s, :), 'Marker', 'o');
    end
    
    xticks(1:numLocations);
    xticklabels({'Up', 'Right', 'Down', 'Left'});
    titleText=['Saccade duration for target ' int2str(k)];
    title(titleText);
    ylabel('Duration (s)');
    xlabel('Saccade direction');
    xlim([0,numLocations+1]);
    ylim([0,0.05]);
end


%Figure 3 & 4 - Plot x and y position on each trial, by location, for each
%target type. If you have a lot of reps this graph will get very messy. Consider
%creating multiple figures to handle large datasets.
numRows = 2;
numCols = 2;

for m = 1:numTargets
    figure();
    for k = 1:numel(trials)
        if trials(k).Target(1) == m
            subplot(numRows, numCols,  trials(k).TargetLocation(1));
            plot3(trials(k).EyeData.ProportionTrial, trials(k).EyeData.LeftEyeXDegrees, trials(k).EyeData.LeftEyeYDegrees, 'b');
            hold on
            plot3(trials(k).EyeData.ProportionTrial, trials(k).EyeData.RightEyeXDegrees, trials(k).EyeData.RightEyeYDegrees, 'r');
            
            %plot target position at onset
            trialDuration = trials(k).EyeData.TimeTag(end)- trials(k).EyeData.TimeTag(1);
            targetx = (trials(k).TargetOnset- trials(k).EyeData.TimeTag(1))/trialDuration;
            targety = (trials(k).TargetLocation(2)-rect(3)/2) * pixelSize;
            targety =  2*atand(targety/(2*displayDistance));
            targetz = -(trials(k).TargetLocation(3)-rect(4)/2) * pixelSize;
            targetz =  2*atand(targetz/(2*displayDistance));
            plot3(targetx, targety, targetz, 'ok'); 
            
            zlabel('Y (degrees)');
            ylabel('X (degrees)');
            xlabel('Proportion of Trial');
            zlim([-15, 15]);
            ylim([-15, 15]);

            titleText=['Position of target ' int2str(trials(k).Target(1)) ' location ' int2str(trials(k).TargetLocation(1))];
            title(titleText);
            legend({'Left eye', 'Right eye', 'Target onset'});
        end
    end
end


%Figure 6 & 7- Velocity on each trial, by location, for each target type.
%If you have a lot of reps this graph will get very messy. Consider
%creating multiple figures to handle large datasets.
numCols = numLocations;
numRows = numReps;

for m = 1:numTargets
    figure();
    repCount= zeros(1,numLocations);
    for k = 1:numel(trials)
        if trials(k).Target(1) == m
            
            %get which iteration this is
            location = trials(k).TargetLocation(1);
            repCount(location)= repCount(location)+1;
            iteration = repCount(location);
            plotPosition = location+(numLocations*(iteration-1));
            
            subplot(numRows, numCols, plotPosition);
            plot(trials(k).EyeData.ProportionTrial, trials(k).EyeData.LeftEyeVelocity, 'b');
            hold on     
            plot(trials(k).EyeData.ProportionTrial, trials(k).EyeData.RightEyeVelocity, 'r');

            %Put some markers to delineate saccade when the participant was
            %saccading or fixating. Lift them up over data.
            statusL = 470; 
            statusR = 490;

            plot(trials(k).EyeData.ProportionTrial, (trials(k).EyeData.LeftEyeFixationAdjusted*statusL), ':b');
            plot(trials(k).EyeData.ProportionTrial, (trials(k).EyeData.RightEyeFixationAdjusted*statusR), ':r');
            plot(trials(k).EyeData.ProportionTrial, (trials(k).EyeData.LeftEyeSaccadeAdjusted*statusL), 'b', 'LineWidth', 5);
            plot(trials(k).EyeData.ProportionTrial, (trials(k).EyeData.RightEyeSaccadeAdjusted*statusR), 'r',  'LineWidth', 5);
            
            %plot target onset as a vertical line
            trialDuration = trials(k).EyeData.TimeTag(end)- trials(k).EyeData.TimeTag(1);
            targetx = (trials(k).TargetOnset- trials(k).EyeData.TimeTag(1))/trialDuration;
            line([targetx, targetx], [0, 500]); 

            %add peak velocity values to plot
            [mL,iL]=max(trials(k).EyeData.LeftEyeVelocity);
            xL = trials(k).EyeData.ProportionTrial(iL);
            plot(xL, mL, 'ob');
            txtL = ['Peak left eye: ' num2str(mL) ' degrees/s'];
            text(0.1, 100, {txtL}, 'FontSize', 8, 'Color', 'b');
            
            [mR,iR]=max(trials(k).EyeData.RightEyeVelocity);
            xR = trials(k).EyeData.ProportionTrial(iR);
            plot(xR, mR, 'or');
            txtR = ['Peak right eye: ' num2str(mR) ' degrees/s'];
            text(0.1, 80, {txtR}, 'FontSize', 8, 'Color', 'r');
            
            ylabel('Velocity (degrees/s)');
            xlabel('Proportion Trial');
            xlim([0,1]);
            ylim([0, 500]);
            if iteration == 1
                legend({'Left eye velocity', 'Right eye velocity', 'Left eye fixation', 'Right eye fixation', 'Left eye saccade', 'Right eye saccade', 'Target onset'}, 'Location', 'southoutside' );
            end

            titleText=['Trial velocity for target ' int2str(trials(k).Target(1)) ' location ' int2str(trials(k).TargetLocation(1))];
            title(titleText);
            
            txt = {'**Data forward and reverse filtered with';'a 10-frame moving window average**'};
            text(.1, 50,txt, 'FontSize', 8, 'Color', [.5, .5, .5])
        end
    end
end


%Figure 7 & 8- Acceleration on each trial, by location, for each target
numCols = numLocations;
numRows = numReps;

for m = 1:numTargets
    figure();
    repCount= zeros(1,numLocations);
    for k = 1:numel(trials)
        if trials(k).Target(1) == m
            
            %get which iteration this is
            location = trials(k).TargetLocation(1);
            repCount(location)= repCount(location)+1;
            iteration = repCount(location);
            plotPosition = location+(numLocations*(iteration-1));
            
            subplot(numRows, numCols, plotPosition);
            plot(trials(k).EyeData.ProportionTrial, trials(k).EyeData.LeftEyeAcceleration, 'b');   
            hold on     
            plot(trials(k).EyeData.ProportionTrial, trials(k).EyeData.RightEyeAcceleration, 'r');
         
            %Put some markers to delineate saccade when the participant was
            %saccading or fixating. Lift them up over data.
            statusL = 32000; 
            statusR = 30000;

            plot(trials(k).EyeData.ProportionTrial, (trials(k).EyeData.LeftEyeFixationAdjusted*statusL), ':b');
            plot(trials(k).EyeData.ProportionTrial, (trials(k).EyeData.RightEyeFixationAdjusted*statusR), ':r');

            plot(trials(k).EyeData.ProportionTrial, (trials(k).EyeData.LeftEyeSaccadeAdjusted*statusL), 'b', 'LineWidth', 5);
            plot(trials(k).EyeData.ProportionTrial, (trials(k).EyeData.RightEyeSaccadeAdjusted*statusR), 'r',  'LineWidth', 5);
            
            %plot target onset as a vertical line
            trialDuration = trials(k).EyeData.TimeTag(end)- trials(k).EyeData.TimeTag(1);
            targetx = (trials(k).TargetOnset- trials(k).EyeData.TimeTag(1))/trialDuration;
            line([targetx, targetx], [-35000,35000]); 

            ylabel('Acceleration (degrees/s^2)');
            xlabel('Proportion Trial');
            xlim([0,1]);
            ylim([-35000,35000]);
            if iteration ==1
                legend({'Left eye acceleration', 'Right eye acceleration', 'Left eye fixation', 'Right eye fixation', 'Left eye saccade', 'Right eye saccade', 'Target onset'}, 'Location', 'southoutside' );
            end
            titleText=['Trial Acceleration for target ' int2str(trials(k).Target(1)) ' location ' int2str(trials(k).TargetLocation(1))];
            title(titleText);
            
            txt = {'**Data forward and reverse filtered with';'a 10-frame moving window average**'};
            text(.1,-20000,txt, 'FontSize', 8, 'Color', [.5, .5, .5])
        end
    end
end



%Figure 9 & 10 - Dilation on each trial, by location, for each target type
numCols = numLocations;
numRows = numReps;

for m = 1:numTargets
    figure();
    repCount= zeros(1,numLocations);
    for k = 1:numel(trials)
        if trials(k).Target(1) == m
            
            %get which iteration this is
            location = trials(k).TargetLocation(1);
            repCount(location)= repCount(location)+1;
            iteration = repCount(location);
            plotPosition = location+(numLocations*(iteration-1));
            
            subplot(numRows, numCols, plotPosition);
            plot(trials(k).EyeData.ProportionTrial, trials(k).EyeData.LeftPupilDiameter, 'b');
            hold on     
            plot(trials(k).EyeData.ProportionTrial, trials(k).EyeData.RightPupilDiameter, 'r');
            
            %Put some markers to delineate saccade when the participant was
            %saccading or fixating. Lift them up over data.
            statusL = 70; 
            statusR = 75;

            plot(trials(k).EyeData.ProportionTrial, (trials(k).EyeData.LeftEyeFixationAdjusted*statusL), ':b');
            plot(trials(k).EyeData.ProportionTrial, (trials(k).EyeData.RightEyeFixationAdjusted*statusR), ':r');

            plot(trials(k).EyeData.ProportionTrial, (trials(k).EyeData.LeftEyeSaccadeAdjusted*statusL), 'b', 'LineWidth', 5);
            plot(trials(k).EyeData.ProportionTrial, (trials(k).EyeData.RightEyeSaccadeAdjusted*statusR), 'r',  'LineWidth', 5);

            ylabel('Dilation (pixels)');
            xlabel('Proportion Trial');
            xlim([0,1]);
            ylim([0,80]);
            if iteration ==1
                legend({'Left pupil diameter', 'Right pupil diameter', 'Left eye fixation', 'Right eye fixation', 'Left eye saccade', 'Right eye saccade'}, 'Location', 'southoutside' );
            end
            titleText=['Pupil dilation for target ' int2str(trials(k).Target(1)) ' location ' int2str(trials(k).TargetLocation(1))];
            title(titleText);
        end
    end
end




%% Step 5 - Write data to csv

rawResults = table();
for k = 1:numel(trials)
    newTrial = trials(k).EyeData;
    newTrial.TrialNumber = repmat(trials(k).Trial, [height(newTrial), 1]);
    newTrial.Target = repmat(trials(k).Target(1), [height(newTrial), 1]);
    newTrial.TargetLocation = repmat(trials(k).TargetLocation(1), [height(newTrial), 1]);
    newTrial.TrialStart = repmat(trials(k).TrialStart, [height(newTrial), 1]);
    newTrial.TargetOnset = repmat(trials(k).TargetOnset, [height(newTrial), 1]);
    newTrial.SaccadeOnset = repmat(trials(k).SaccadeOnset, [height(newTrial), 1]);
    newTrial.SaccadeEnd = repmat(trials(k).SaccadeEnd, [height(newTrial), 1]);
    newTrial.SaccadeReactionTime = repmat(trials(k).SaccadeReactionTime, [height(newTrial), 1]);
    newTrial.SaccadeDuration = repmat(trials(k).SaccadeDuration, [height(newTrial), 1]);
    newTrial = newTrial(:, [40:end, 1:37]);
    rawResults = [rawResults; newTrial];
end

csvFileID = [fileName '_Results.csv'];
writetable(rawResults, csvFileID);

%% Step 6 - Finish up

Datapixx('SetTPxSleep');
Datapixx('RegWrRd');
Datapixx('Close');

end