function mapping=two_frame_mapping(P_0,P_1,reg_of_motion)
% Finds the mapping between particles detailed in P_0 and P_1 based on
% closest distance and relative size.
% P_0, P_1 - Matrices listing the labels, centers of mass position and
% sizes of particles.
% reg_of_motion - maximum distance between particles at different times to
% be considered to be identical
% mapping - two column matrix listing labels of particles in P_0 and P_1
% identified to be identical.
%
%For bugs, comments, etc. please contact Kapil Krishan :
%kkrishan@uci.edu
%www.physics.uci.edu/~foams
%version: I

%Do mapping going 0->1
map_01=[];
for ptcl_ct=P_0(:,1)'
    [mapped_ptcl,min_dist,n_sol]=two_frame_1ptcl_mapping(ptcl_ct,reg_of_motion,P_0,P_1);
    map_01=[map_01;ptcl_ct mapped_ptcl(1) min_dist n_sol];
end

%Do mapping going 1->0
map_10=[];
for ptcl_ct=P_1(:,1)'
    [mapped_ptcl,min_dist,n_sol]=two_frame_1ptcl_mapping(ptcl_ct,reg_of_motion,P_1,P_0);
    map_10=[map_10;ptcl_ct mapped_ptcl(1) min_dist n_sol];
end

mapping=intersect(map_01,map_10(:,[2 1 3 4]),'rows');

mapping(find(mapping(:,3)<0),:)=[];%Dont map vanishing particles
mapping(find(mapping(:,4)>1),:)=[];%Dont map multiple solutions

mapping=mapping(:,[1 2]);

function [mapped_ptcl,min_dist,n_sol]=two_frame_1ptcl_mapping(ptcl_ct,reg_of_motion,P_0,P_1)

orig_ptcl=P_0(ptcl_ct,:);

c1=find(abs(P_1(:,2)-orig_ptcl(2))<reg_of_motion(1)/2);
c2=find(abs(P_1(:,3)-orig_ptcl(3))<reg_of_motion(2)/2);

ptcl_rom=P_1(intersect(c1,c2),:);

if size(ptcl_rom,1)>0
    [min_dist,idx]=min(sqrt(sum(((ptcl_rom(:,[2 3])-ones(size(ptcl_rom,1),1)*orig_ptcl([2 3])).^2)')));
    clst_ptcls=ptcl_rom(idx,:);

    idx=1;
    if size(clst_ptcls,1)>1
        [blah,idx]=min(abs(clst_ptcls(:,4)-orig_ptcl(4)));
    end

    n_sol=length(idx);
    mapped_ptcl=clst_ptcls(idx(1),:);
        
else
    mapped_ptcl=zeros(size(orig_ptcl));
    n_sol=0;min_dist=-1;
end