function [M] = matchNeurons(masks1, masks2, threshold)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Match neurons between two sets of labels by first calculating an
%Intersection of the union based measurement of distance between 2 masks
%and then use the Hungarian algorithm to determine the best matches
% Inputs:
%     - masks1: MxNxT matrix of the masks by one grader 
%     - masks2: MxNxT matrix of the masks by a second grader on the same
%     video
%     - threshold: Threshold of the lowest intersection of the union that
%     will be tolerated for a match

% Output:
%     - M: Indexes of the matching neurons

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sm1_1 sm1_2 sm1_3] = size(masks1);
[sm2_1 sm2_2 sm2_3] = size(masks2);

for m1 = 1:sm1_3
    mask1 = masks1(:, :, m1);
    m1i{m1} = find(mask1);
end

for m2 = 1:sm2_3
    mask2 = masks2(:, :, m2);
    m2i{m2} = find(mask2);
end

dist = NaN(sm1_3, sm2_3);
for m1 = 1:sm1_3
    for m2 = 1:sm2_3
        m1ind = m1i{m1};
        m2ind = m2i{m2};
        union = length(unique(cat(1, m1ind, m2ind)));
        int1 = intersect(m1ind, m2ind);
        int2 = intersect(m2ind, m1ind);
        if length(m1ind) > 0 && length(m2ind) > 0
            if length(int1)==length(m1ind) || length(int2)==length(m2ind)
                dist(m1, m2) = 0;
            elseif length(int1)/union > threshold
                dist(m1, m2) = 1-length(int1)/union;
            elseif isempty(int1) == 1
                dist(m1, m2) = Inf;
            end
        end
    end
end
dist(isnan(dist)) = Inf;

M = matchpairs(dist,10);
matchind = sub2ind(size(dist), M(:, 1), M(:, 2));

end