% "Olfactory behaviors assayed by computer tracking of Drosophila in a 
% four-quadrant olfactometer".
% 
% AUTHORS: 
% Chun-Chieh Lin, Olena Riabinina, Christopher J.Potter
% 
% The Solomon H. Snyder Department of Neuroscience, Center for Sensory 
% Biology, Johns Hopkins University School of Medicine, 
% Baltimore, MD 21205, USA. cpotter@jhmi.edu

% Remove all datapoints that move at a speed below 0.163 cm/s for more 
% than 3.3 sec 

% Set cut-off threshold for speed. For our spatial and temproal resolution
% of datam 0.163cm/s corresponds to 0.15 pixels/frame
v_threshold=0.15;

% Set the duration of the low speed window. For our data, 3.3s corresponds
% to 100 frames.
v_duration=100;


% load(data); % this should load a file with variables nx, ny and nt,
            % that represent x and y coordinates and time code of all tracked 
            % flies in the course of an experiment.  
            

            
% creat a copy of x,y coordinates and time
x=nx;
y=ny;
t=nt;

% calculate change of x,y, and t between successive datapoints
dx=diff(x);
dy=diff(y);
dt=diff(t);

% calculate vx and vy speeds as change of x and y coordinates,
% respectively, over time
vx=dx./dt;
vy=dy./dt;

% Calculate total speed as a vector sum of vx and vy. Apply 10-point
% running average to the values of speed
v=smooth(sqrt(vx.^2+vy.^2),10);

% v_filt will be equal to 0 for values of v that are above v_threshold, and
% equal to 1 for values below v_threshold
v_filt=(v<v_threshold);

% a_filt will be equal to 1 at time points when v_filt switched from 0 to
% 1. a_filt will be equal to -1 at time points when v_filt switched from 1 to
% 0. Other elements of a_filt will be equal to 0.
a_filt=diff(v_filt);

% ind1 will contain indexes of the elements of a_filt that are equal to 1
ind1=find(a_filt==1);

% ind_1 will contain indexes of the elements of a_filt that are equal to -1
ind_1=find(a_filt==-1);

% code below ensures that there is equal number of elements in ind1 and
% ind_1.
if ind1(1)>ind_1(1)
    ind1=[1; ind1];
end;

if ind1(end)>ind_1(end)
    ind_1=[ind_1; length(a_filt)];
end;

% calculate duration (in frames) of low speed bouts by subtracting values
% of indexes of ind_1 and ind1.
stop_dur=ind_1-ind1;

% find bouts that last longer than the cut-off duration v_duration
ind_stop=find(stop_dur>v_duration);


% create a vector ind_del, that contains indexes of x,y, and t that
% correspond to slow moving datapoints
ind_del=[];
for i=ind_stop
    ind_del=[ind_del ind1(i):ind_1(i)];
end;

% remove slow moving datapoints from the dataset
x(ind_del)=[];
y(ind_del)=[];
t(ind_del)=[];

% create new figure
figure

% plot remaining datapoints in grey, and removed datapoints - in red
plot(x,y,'.', 'MarkerSize', 6, 'MarkerEdgeColor', [0.7, 0.7, 0.7])
hold on
plot(nx(ind_del),ny(ind_del),'or', 'MarkerSize', 6, 'MarkerFaceColor', 'r')
axis equal
