function [gangmask, idx] = makeGanglia_kmeans(matname)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Goes through all of the videos for which masks have been selected and
%assembles the masks into ganglia by identifying nearby neighbors
% Inputs:
%     - matname: the name of the .mat file containing the masks
% Output:
%     - gangmask: cell array where each cell contains the masks of the
%     neurons that belong to that cell's ganglia
%     - idx: cluster indices corresponding to unique ganglia
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


gangmask = {};
alpha = 0.01;
count = 1;
load(matname)

[X, Y] = findCOM(double(masks));
COMpos = cat(2, X, Y);
inertia = [];
metric = [];
for clust = 1:min([30 size(masks, 3)])
    [idx, ~, ~, D] = kmeans(COMpos, clust, 'Distance', 'cityblock');
    inertia(clust) = trace(D(:, idx));
    metric(clust) = (inertia(clust)/inertia(1))+alpha*clust;
end

[~, metind] = min(metric);
[idx, ~, ~, ~] = kmeans(COMpos, metind, 'Distance', 'cityblock');

for clust = 1:max(idx)
    cind = find(idx == clust);
    gangmask{count} = masks(:, :, cind);        
    count = count+1;
end

amask = [];
gcount = [];
amask = double(logical(sum(masks, 3)));
gcount = 2;
for cind = 1:length(gangmask)
    smask = logical(sum(gangmask{cind}, 3));
    amask(find(smask)) = gcount;
    gcount = gcount + 1;
end