% SEGMDEBRE Segmentation of debree. 
% [cellbw,WAT] = SEGMDEBRETH(IM,DIAMLOW,DIAMHIGH) finds the debree 
% in thE image IM, with diameter between DIAMLOW and DIAMHIGH using 
% adaptive thresholding. 
% Returns the segmented image in cellbw and the label image WAT. 
% 
% [cellbw,WAT] = SEGMDEBRETH(IM,DIAMLOW,DIAMHIGH,'RAD',RAD) uses RAD as the
% radius in the adaptive filtering
% 
% [cellbw,WAT] = SEGMDEBRETH(IM,DIAMLOW,DIAMHIGH,'ADTH',ADTH) uses ADTH as 
% the threshold in the adaptive filterin.
% 
% [cellbw,WAT] = SEGMDEBRETH(IM,DIAMLOW,DIAMHIGH,'FILTEROPT',FILTEROPT) uses
% FILTEROPT as the filter option for filtering.
% FILTEROPT = 'GAUSSIAN' or
% FILTEROPT = 'MEDIAN'
%
%
function [cellbw,wat] = segmdebreth(varargin)


test = 0;

% parse inputs
im = varargin{1};
diamlow = varargin{2};
diamhigh = varargin{3};
% default values
prm.rad = 20;
prm.adth = 0.15;
filteropt = 'gaussian';
conn = 4;
for i = 4:2:nargin
    name = varargin{i};
    switch(name)
        case('prm')
            prmin = varargin{i+1};
%         case('prm.rad')
%             rad = varargin{i+1};
%         case('adth')
%             adth = varargin{i+1};
%         case('filteropt')
%             filteropt = varargin{i+1};
        otherwise
            error('Wrong option to SEGMDEBRETH')
    end;
end;


% from input    
prm = mergestruct(prm,prmin);

disp(sprintf('Radius filter :  %6.3f',prm.rad));
disp(sprintf('Threshold filter :  %6.3f',prm.adth));            

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

% number of pixels
thnumpixlow = round((4/3)*pi*(diamlow/2)^3);
thnumpixhigh = round((4/3)*pi*(diamhigh/2)^3);

% standard deviation in gaussian filter
stdVal = 1;

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

imini = im;
% filter image
disp('Filter image');
if strcmp(filteropt,'gaussian')
    disp('Using gaussian filter')
    im = imfilter(im,fspecial('gaussian',3,stdVal),'replicate');
elseif strcmp(filteropt,'median')
    disp('Using median filter');
    im = medfilt2(im,[3 3]);
else
    error('Wrong filter option to SEGMDEBRETH')
end;


% % make the spots sharper
% se = strel('disk',3);
% imtop = imtophat(im,se);

% show(im,1)
% show(imtop,2)
% pause

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

% adaptive thresholding
% kim1 = adaptfilt(im,rad,th);
kim2 = adaptfilt(im,prm.rad,prm.adth);
kIm = kim2;
% kIm = gt(kim1 + kim2,0);
% 
% show(imini,1)
% show(im,2)
% show(kim2,4)
% show(kIm,5)
% pause



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

% open to disconnect small bridges

load(['ball' int2str(prm.opendiam)]);
se = getball(ball,prm.opendiam,1);
kIm = imopen(kIm,se);

% remove small and big
% use same CONN as in the tracking, else you get SMALL particles!!!!
kIm = bwareaopen(kIm,thnumpixlow,conn);
kIm = kIm - bwareaopen(kIm,thnumpixhigh,conn);
disp(sprintf('Threshold number pixels low : %6.2f',thnumpixlow))
disp(sprintf('Threshold number pixels high : %6.2f',thnumpixhigh))

% out variable
[wat,L] = bwlabeln(kIm);
cellbw = kIm;
if test == 1
    show(imini,1)
    show(im,2)
    show(kim2,3);title('After adaptive')
    show(cellbw,6);title('final segmentation')
    pause
end;

