% SEGMDEBRE Segmentation of debree. 
% FUNCTION [CELLBW,WAT] = SEGMDEBRE(IM,DIAM) finds the debree in the
% image IM, with diameter DIAM using Kmeans. Returns the segmented image in CELLBW 
% (Like the cell segmentation) and the label image WAT.
%
function [cellBW,wat] = segmDebreKmeans(im,diam)

% dimension of image
[M N O] = size(im);

% set to 1 if we want to remove satureated pixels
sat = 1;

% set to 1 if speedup
spUp = 1;

% number of pixels
thNumPix = (4/3)*pi*(diam/2)^3;

% how small and big around mean
sc = 7;

% standard deviation in gaussian filter
stdVal = 1;

% scale to 0-1
im = scale(im);

% remove saturated pixels
% disp('Removing saturated pixels')
% if sat == 1
%     im(find(im == 1)) = 0;
% end;

% filter image
disp('Filter image');
im = imfilter(im,fspecial('gaussian',3,stdVal),'replicate');


% mid plane
midZ = round(O/2);

% clustering
disp('Cluster');
if spUp == 1
    % finding the threshold for one plane, and use for the rest of the
    % planes
    [kIm,th] = clustKmeans(im(:,:,midZ),2);
    kIm = ge(im,th(1));
    thCT = th(1);
else
    % working on the whole range
    [kIm,th] = clustKmeans(im,2);
    p = 1;
    kIm = ge(kIm,p);
    thCT = p;
end;

% show(kIm,2)

% remove small and big
thLow = round(thNumPix/sc);
kIm = bwareaopen(kIm,thLow);
thHigh = round(thNumPix*sc);
kIm = kIm - bwareaopen(kIm,thHigh);

% close small holes
load ball1;se = getball(ball,1,1);
kIm = imclose(kIm,se);


% fill small holes
for j = 1 : O
    regHere = kIm(:,:,j);
    diff = imfill(regHere,'holes') - regHere;
    diff = diff - bwareaopen(diff,thLow);
    regHere(find(diff)) = 1;    
    kIm(:,:,j) = regHere;
end;


% out variable
[wat,L] = bwlabeln(kIm);
cellBW = kIm;

%show(im,1)
%show(kIm,3)
%pause

