mpib_memoreye/gaze_shift/detect_gaze_shifts_jld_size.m

function [num_shifts, shift_directions, sizeGa] = detect_gaze_shifts_jld_size(data)
% Input:
%   - data: a 2 x N matrix containing the eye movements in the vertical
%           and horizontal axes, respectively, over time.
% Output:
%   - num_shifts: the total number of gaze shifts detected in the data.
%   - shift_directions: a 1 x num_shifts vector containing the direction of
%                       each gaze shift, where 1 indicates a rightward shift,
%                       -1 indicates a leftward shift, 2 indicates an upward
%                       shift, and -2 indicates a downward shift. 
%                       NaN indicates that no direction was detected.
sizeGa=[];

% Set parameters
sampling_rate = 500; % in Hz
vel_threshold = 5 * median(sqrt(sum(diff(data,1,2).^2,1)))/ sampling_rate; %5 is the original values
min_saccade_interval = 100 / 1000 * sampling_rate; % in samples
smooth_window_size = 7 / 1000 * sampling_rate; % in samples

% Euclidean distance between succesives time points
edP=sqrt(sum(diff(data,1,2).^2,1));
edP=edP/sampling_rate;

% Smooth the data
data_smooth = smoothdata(edP', 'gaussian', smooth_window_size)'; %Euclidean D smooth
ori_data_smooth = smoothdata(data', 'gaussian', smooth_window_size)'; %Original data smooth

% data_vel = diff(data_smooth, 1, 2);
data_vel = data_smooth;


% Find saccade onsets
% vel_norm = sqrt(sum(data_vel.^2, 1));
saccade_onsets = find(data_vel > vel_threshold);
if isempty(saccade_onsets)
    num_shifts = 0;
    shift_directions = [];
    return;
end

% Remove saccades that are too close in time
saccade_diff = diff(saccade_onsets);
toremove=find(saccade_diff < min_saccade_interval)+1; %index to remove
saccade_onsets(toremove)=[];

% Find saccade offset (defined as velocity smaller than threshold)
all_saccade_offsets = find(data_vel < vel_threshold);

% Finding the next offset to each onset
saccade_offsets=[];
for i=1:length(saccade_onsets)
    found=find(all_saccade_offsets>saccade_onsets(i)); %Find the next closet offset
    if isempty(found) %only when the epoch finish
        saccade_offsets(i)=length(data_vel);
    else
        saccade_offsets(i)=all_saccade_offsets(found(1));        
    end
end

% Determine the direction of each gaze shift
num_shifts = length(saccade_onsets);
shift_directions = nan(1, num_shifts);
sizeGa = nan(1, num_shifts);

for i = 1:num_shifts
    
    saccade_start = saccade_onsets(i);
    saccade_end = saccade_offsets(i);


    posOnset=ori_data_smooth(:,saccade_onsets);
    posOffset=ori_data_smooth(:,saccade_end);
    
    dx = posOffset(1) - posOnset(1);
    dy = posOffset(2) - posOnset(2);
    angle = atan2(dy, dx); %Calculating the distance in radians

    shift_directions(i) = angle;
    sizeGa(i) = norm(posOnset - posOffset);
    
end
end