function [temp]=AttractionIndex(x,y, t, n, xc, yc, xa, ya) 

% "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

% This function assigns data to quadrants, and calculates the Attraction Index
% and the Percentage Index. The input data are:
% x,y, t - filtered and cleaned up x,y coordinates and time of all 
% datapoints;
% n - a number between 1 and 4, that indicates which quadrant contained 
% test odorant
% xc, yc- coordinates of the four circles, calculated in the fitting
% procedure in MaskSpatialFiltering.m 
% xa, ya - coordinates of the centre of the arena, calculated in the
% fitting procedure in MaskSpatialFiltering.m 


% select datapoints that correspond to the last 5 mins (300 s) of the 
% recording. Data recorded at 30 frames per second.

inds=find(t<= max(t) & t>max(t)-300*30);
x=x(inds);
y=y(inds);
t=t(inds);

% calculate coordinates of vectors from the centre on the arena to each 
% of the remaining datapoints
data_r=[(x-xa), (y-ya)];

% calculate coordinates of vectors from the centre of the arena to the
% centres of masking circles
circle_r=[(xc-xa), (yc-ya)];

% Create a matrix that contains signs (1 or -1) of vector multiplication 
% of the datapoint vectors by vectors to the centres of circles. Two of
% these values will always be +1, and two will be -1. If a datapoint lays
% in Quadrant 1, its vector will be between the first two vectors to
% centres of circles, and the matrix row that corresponds to this datapoint
% will be [1 1 -1 -1]. Datapoints in the second quadrant will corespond to
% [-1 1 1 -1], third - [-1 -1 1 1], and fourth [1 -1 -1 1].

matr=sign(data_r*circle_r');

% check if n is a valid number, if not - assign it to 1
if any(n==[1,2,3,4])
    n
else n=1
end;

% find indexes of matr that correspond to the test quadrant n
switch n
    case 1
        ind=find(matr(:,1)==1 & matr(:,2)==1);
    case 2
        ind=find(matr(:,2)==1 & matr(:,3)==1);
    case 3
        ind=find(matr(:,3)==1 & matr(:,4)==1);
    case 4
        ind=find(matr(:,4)==1 & matr(:,1)==1);
end;
  
% calculate attraction index and percentage index
ai=(numel(ind)-(numel(x)-numel(ind))/3)./(numel(ind)+(numel(x)-numel(ind))/3);
perc_i=numel(ind)./numel(x);

% display values for attraction index. See Protocol for explanations of the
% two alternative formulas for AI.
disp(sprintf('AI=%6.4f', ai))
disp(sprintf('PI=%6.4f', perc_i))


% to plot Attraction Index as a function of time, create vector time_bins,
% and calculate how many datapoints in total, and how many datapoints from
% the test quadrant occured during each time bin

time_step=300;% 10sec intervals
time_bins=min(t):time_step:max(t); 

ai=[];
perc_i=[];

for i=time_bins(1:end-1)
    % find indexes of all datapoints that occured during the time bin
    ind_tot=find(t>i & t<=i+time_step);
    % find indexes of datapoints from the test quadrant 
    ind_test=intersect(ind, ind_tot);
    
    % calculate attraction indexes for this time bin, andadd them to the
    % array of attraction indexes
    ai=[ai; (numel(ind_test)-(numel(ind_tot)-numel(ind_test))./3)./...
        (numel(ind_test)+(numel(ind_tot)-numel(ind_test))./3)];
    
    perc_i=[perc_i; numel(ind_test)/numel(ind_tot)];
end;

% plot out attraction indexes as a function of time
figure
plot(time_bins(1:end-1)/30, ai, '-or');
hold on
plot(time_bins(1:end-1)/30, perc_i, '-*b');
axis([min(t)/30 max(t)/30 -1 1])
axis square
xlabel('Time, s')
ylabel('Behavioral Indexes')
legend('AI', 'PI')
st{1}='Attraction Index and Percentage Index,';
st{2}='for the last 5 mins of experiment';
title(st,'Interpreter', 'none')
FormatAndSaveFig('AI_PI');
end;










