function [num_shifts, shift_directions] = detect_gaze_shifts(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.

% Set parameters
sampling_rate = 500; % in Hz
vel_threshold = 5 * median(sqrt(sum(diff(data,1,2).^2,1)))/ sampling_rate;
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)=[];


% Determine the direction of each gaze shift
num_shifts = length(saccade_onsets);
shift_directions = nan(1, num_shifts);
for i = 1:num_shifts
    saccade_start = saccade_onsets(i);
    if i < num_shifts
        saccade_end = saccade_onsets(i+1) - 1;
    else
        saccade_end = saccade_onsets(i)+1;
    end

    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;
    
end
end