function [n1,n2,vol1,vol2] = colocalize(bw1,bw2,relminvol,voxelvol)
% Function to find colocalization of vesicles
%

dim = size(bw1);
ndim = numel(dim);
if ndim == 2
    conn = 4;
else
    conn = 6;
end;

% label and count size
[vol1,faser1] = bwsize(bw1,conn);
L1 = max(faser1(:));
[vol2,faser2] = bwsize(bw2,conn);
L2 = max(faser2(:));

%
% Find overlap
%
val = faser1(bw2 == 1);
val(val == 0) = [];
valu = unique(val);
n = numel(valu);
N1 = NaN(n,1);
for i = 1 : n
    N1(i) = sum(val == valu(i));
end;
n1 = zeros(L1,1);
n1(valu) = N1;

val = faser2(bw1 == 1);
val(val == 0) = [];
valu = unique(val);
n = numel(valu);
N2 = NaN(n,1);
for i = 1 : n
    N2(i) = sum(val == valu(i));
end;
n2 = zeros(L2,1);
n2(valu) = N2;

% remove small colocalization particles
ind = n1./vol1 < relminvol;
n1(ind) = 0;
ind = n2./vol2 < relminvol;
n2(ind) = 0;

n1 = n1*voxelvol;
n2 = n2*voxelvol;
vol1 = vol1*voxelvol;
vol2 = vol2*voxelvol;





