function [minima,minimacell,prm] = getminima(varargin)
% GETMINIMA Find markers for cell segmentation
%
%   [MINIMA,MINIMACELL,PRM] = GETMINIMA(IM,MINIMACELL,PRM) 
%   Finds the markers for image IM, cell markers
%   MINIMACELL, and settings in struct PRM. Let PRM = [] for default
%   settings and let MINIMACELL = [] for no cell markers. If MINIMACELL is
%   non-empty and prm.method = 'automatic', the method only finds the
%   background markers.
%
%   [...] = GETMINIMA(...,IMNUCL) also read the nucleus image to find 
%   markers. For this option to be used, set prm.method = 'nucleus'. 
%
%   General fields for PRM
%
%   PRM.H (default = [0.5 0.5 1.5])            
%   Voxelsize 1x3 vector
% 
%   PRM.MINVOLFULL  (default = 5)
%   Lower cell volume in 1000 mcm3.
%
%   PRM.MAXVOLFULL (default = 50);
%   Upper cell volume in 1000 mcm3.
% 
%   PRM.METHOD (default = 'automated')
%   The markers are created automatically ('automated'), 
%                       by the nucleus channel ('nucleus') or the cell 
%                       markers are given manually ('manual'). This option
%                       requires a non-empty binary image for the cell
%                       markers, and only the background markers are found 
%                       by the program.
%
%   PRM.LEVEL (default = 0.3)
%   Normalized plane level 0 -> 1 to make markers. 
%                       If PRM.LEVEL = 'ALL' all planes 
%                       are used to create make markers. 
%
%   Specific settings for PRM.METHOD = 'automated':
%
%   PRM.ADAPTFILTTH (default = 0.01)
%   Threshold in adaptive filtering. Default 0.01
%
%   Specific fields in PRM for PRM.METHOD = 'nucleus'
%   PRM.NUCLEUS*      : Any option to SEGMCT for segmentation of the
%                       nuclei. See SEGMCT for options. Default settings
%                       are PRM.NUCLEUSMETHOD = 'itth', PRM.NUCLEUSSPLIT =
%                       0, PRM.NUCLEUSITTH = 1.
%
%
%   See also segmct, classifycells, segmsurfwat
% 	=======================================================================================
% 	CELLSEGM: Cell segmentation, Version 2012
% 	Copyright (c): Erlend Hodneland
%   Jonas Lies Vei 91,
%   5009 Bergen
%   Email: erlend.hodneland@biomed.uib.no
% 	=======================================================================================

msg = ['This is ' upper(mfilename) ' finding minima'];
disp(msg);


imsegm = varargin{1};
minimacell = varargin{2};
prmin = varargin{3};
imnucl = [];

if nargin == 4
    imnucl = varargin{4};
end;
dim = size(imsegm);
if numel(dim) == 2
    dim = [dim 1];
end;

%
% Minima parameters
%

% minimum and maximum volume
prm.minvolfull = 3;
prm.maxvolfull = 50;

% I change this 20120320 by adding 0.5 in front!
% prm.volfactor = 0.50;

% pixelsize
prm.h = [0.5 0.5 1.5];

% NB may not work if this is too high, flooding goes out
prm.rangemin = 0.4;

% how quickly to decrease the minima
prm.dimminball = 2;

% method
prm.method = 'automated';

% normalized level to make minima
prm.level = 0.30;

% adjust for 3D size
prm.just = 0.9;

% the threshold for the adaptive filter to find the minima(0.01-0.03)
prm.automatedadth = 0.01;

% settings for segmct in nucleus segmentation
prm.nucleusmethod = 'thrs';
prm.nucleusthrsth = 1.0;

if numel(dim) == 2
    dim = [dim 1];
end;
prm.dim = dim;

% assign the input variables
if ~isempty(prmin)
    prm = mergestruct(prm,prmin);
end;
prm.minvolfull = prm.minvolfull*1000;
prm.maxvolfull = prm.maxvolfull*1000;

% voxel volume
prm.voxelvol = prod(prm.h);

% adjust cell volume
[prm.minvol,prm.minvolvox,prm.maxvol,prm.maxvolvox] = ...
    cellsize(prm.minvolfull,prm.maxvolfull,prm.h,prm.just,dim(3));


% adjust for reduced 3D
if ~isequal(prm.level,'all')    
    prm.minvol = round(prm.minvol/(dim(3)^prm.just));
    prm.minvolvox = round(prm.minvol/prm.voxelvol);
    prm.maxvol = round(prm.maxvol/(dim(3)^prm.just));    
    prm.maxvolvox = round(prm.maxvol/prm.voxelvol);    
end;

if isempty(prm.plane)
    if isequal(prm.level,'all')
        prm.plane = 1 : dim(3);
    else
        prm.plane = max(round(prm.level * dim(3)),1);
    end;
end;

msg = ['Using settings'];
disp(msg);
printstructscreen(prm);

msg = ['Minima plane ' int2str(prm.plane)];
disp(msg);

% if isequal(prm.method,'threshold')
%     msg = ['Using threshold method for markers'];
%     disp(msg);
% 
%     mask = zeros(dim);
%     % dont start at bottom! Too weak signals often
%     for i = min(dim(3),5) : dim(3)
%         imhere = imsegm(:,:,i);
%         imhere = scale(imhere);
%         mask(:,:,i) = lt(imhere,prm.thresholdlevel);
%     end;
%     load ball2;se = getball(ball,2,2);
%     mask = imclose(mask,se);
%     mask = imopen(mask,se);
%     
%     minima = mask;
% %     just = 0.9;
%     minima = bwareaopen(minima,prm.minvolvox);
%

if isequal(prm.method,'manual')
    msg = ['Using manually defined cell markers'];
    disp(msg);
    
    % get background
    minimabck = getbackground(imsegm,minimacell,prm);

    % add background and cell
    minima = gt(minimabck + minimacell,0);
    
elseif isequal(prm.method,'nucleus')  
    
    msg = ['Using nucleus method for markers'];
    disp(msg);
    
    % add nucleus information    
    minimacell = findminimanuecleus(imnucl,prm);
    
    % get background minima
    minimabck = getbackground(imsegm,minimacell,prm);
    
    % add background and cell
    minima = gt(minimabck + minimacell,0);

    if test
        'After remove small-final'
        showall(imhere,minimabck)
    end;

%         showall(minima,minimabck)
elseif isequal(prm.method,'automated')
    
    msg = ['Using automated method for markers'];
    disp(msg);

    if isempty(minimacell)
        [minima,prm] = findminimaplane(imsegm,prm);
    else
        msg = ['Only finding background since cells are known'];
        disp(msg);
        minimabck = getbackground(imsegm,minimacell,prm);
        minima = minimacell;
        minima(minimabck == 1) = 1;
    end;
        
else
    error(['Wrong option for prm.method in ' mfilename]);
end;

%--------------------------------------------------

function [im] = repeatplane(im,start,erorad)

% dim = size(im);
% add = start-1;
% dimnew = dim;
% dimnew(3) = dimnew(3) + add;
% imnew = zeros(dimnew);
% imnew(:,:,start:end) = im;
name = ['ball' int2str(erorad)];
load(name);se = getball(ball,erorad,1);
imhere = im(:,:,start);
for i = start-1:-1:1    
    imhere = imerode(imhere,se);
    im(:,:,i) = imhere;
end;

% showall(im)


%-----------------------

function [minimabck] = getbackground(imsegm,minimacell,prm)


test = 0;

% find background
imhere = scale(imsegm);    
th = 0.01;
totvol = numel(imsegm);
while 1        
    minimabck = lt(imhere,th);
    vol = nnz(minimabck);
    relvol = vol/totvol;
    if relvol > 0.10
        msg = ['Threshold for background: ' num2str(th)];
        disp(msg);
        break;
    end;
    th = th + 0.01;
end;

if test
    'After threshold'
    showall(imhere,minimabck)
end;

% close with one
load ball1;se = getball(ball,1,1);
minimabck = imclose(minimabck,se);

% open to disconnect to cells
load ball4;se = getball(ball,4,1);
minimabck = imopen(minimabck,se);

% remove small background minima
% was at 2
minimabck = bwareaopen(minimabck,round(prm.minvolvox/8));

if test
    'After remove small'
    showall(imhere,minimabck,minimacell)
end;

% remove the overlap to cell markers
[faser,L] = bwlabeln(minimabck);
% load ball1;se = getball(ball,1,1);
% overlap = faser .* imdilate(minimacell,se);
overlap = faser .* minimacell;
overlap = overlap(:);
overlap(overlap == 0) = [];
overlap = unique(overlap);
for i = 1 : numel(overlap)
    minimabck(faser == overlap(i)) = 0;
end;

if test
    'After remove overlap'
    showall(imhere,minimabck,minimacell)
end;


% minimabck = imfill(minimabck,'holes');
minimabck = fillholes(minimabck,0,prm.minvolvox,4);
load ball2;se = getball(ball,2,1);
minimabck = imopen(minimabck,se);    
minimabck = imclose(minimabck,se);  

if test
    'After fill,open.close'
    showall(imhere,minimabck)
end;


% remove small background minima
minimabck = bwareaopen(minimabck,round(prm.minvolvox/2));

if test
    'Final'
    showall(imhere,minimabck)
end;

%---------------------------------------------------

function [minima] = findminimanuecleus(imnucl,prm)

vis = 0;

dim = size(imnucl);

% find minima in the nucldus staining using cell segmentation of CT cells!
% msg = ['In FINDMINIMANUCLEUS using plane ' int2str(prm.plane)];
% disp(msg);

imhere = imnucl(:,:,prm.plane);

% was at 10
% p1 = 1;

% segmct
prmin = selectfield(prm,'nucleus');
prmin.h = prm.h;
% [prmin.minvol,prmin.minvolvox,prmin.maxvol,prmin.maxvolvox] = ...
%     cellsize(prm.minvolfull,prm.maxvolfull,prm.h,prm.just,dim(3));
[minimahere,wat] = segmct(imhere,0.0001*prm.minvolfull,0.001*prm.maxvolfull,'prm',prmin);

if vis
    'after segmct'
    showall(imnucl(:,:,prm.plane),minimahere)
end;

% % can be too big sometimes and go outside the border!
% load ball1;se = getball(ball,1,1);
% minimahere = imerode(minimahere,se);
% dilate the weak pixels outside the nucleus to get away from the WGA
% around the nucleus

% imhere = imsegm(:,:,prm.plane);
% for i = 1 : 3
%     minimahere = dilateweak(minimahere,imhere,1);
% end;

% save test minimahere
% 
% load ball2;se = getball(ball,2,1);
% minimahere = imopen(minimahere,se);

if vis
    'after open'
    showall(imnucl(:,:,prm.plane),minimahere)
end;


minima = zeros(dim);
minima(:,:,prm.plane) = minimahere;

% repeat the plane downward to improve the segmentation
if numel(prm.plane) == 1
    minima = repeatplane(minima,prm.plane,1);
end;

% showall(minima)
% % cut the nucleus marker in z direction to not connect with the background
% msg = ['Cut lower and upper parts of minimacell since Hoechst extends the cell'];
% disp(msg);
% [faser,L] = bwlabeln(minima);
% minima = zeros(dim);
% for i = 1 : L
%     
%     reghere = faser == i;
%     [range,minz,maxz] = bwrange(reghere);
%     if range >= 5
%         reghere(:,:,min(minz:minz+2,dim(3))) = 0;
%         reghere(:,:,max(maxz-3:maxz,1)) = 0;
%     end;
%     
%     minima(reghere) = 1;
%     
% end;


% remove very small ones after erotion
% % minima = bwareaopen(minima,round(prm.minvolvox*0.01));
% minima = bwareaopen(minima,round(prm.minvolvox*0.10));

if vis
    'after remove small'
    showall(imnucl(:,:,prm.plane),minimahere)
end;

%-------------------------------------------------------

function [bw] = dilateweak(bw,im,p)

name = ['ball' int2str(p)];load(name);
se1 = getball(ball,p,1);

name = ['ball' int2str(p+1)];load(name);
se2 = getball(ball,p+1,1);

dim = size(bw);
[faser,L] = bwlabeln(bw);
dilall = zeros(dim);
overlap = zeros(dim);
for i = 1 : L
    reghere = faser == i;
    
    dil = imdilatefast(reghere,se1);
    dil(reghere) = 0;
    meanint = mean(im(reghere));
    
    
    % remove the strong ones
    ind = im > 0.8*meanint;    
    dil(ind)  = 0;    
    dilall = dilall + dil;
    
    % check for overlap
    dil = imdilatefast(reghere,se2);
    overlap = overlap + dil;
    
%     showall(im,reghere,dil,dilall)
end;
overlap = gt(overlap,1);


% remove overlap
dilall(overlap == 1) = 0;
% showall(im,bw,dilall,overlap)
bw(dilall == 1) = 1;




%--------------------------------------------------



% FINDMINIMA3D Find the minima to prepare for segmemaxntation
% FINDMINIMA3D(MINKAPPA,prm) Finds the minima from 
% Fields to give:
% prm.ADAPTFILTTH   : Threshold in adaptive filterin
% prm.HX            : pixelsize in xy plane
% prm.HZ            : Pixelsize in z direction
% prm.DIAMLOW       : Lower limit for number of pixels in a minima
% prm.DIAMHIGH      : Upper limit for number of pixels in minima
% 
function [minima] = findminima3d(im,prm)

global test

test = 0;

[M N O] = size(im);

% adjust for the not full 3D
just = 0.9;

% dimension of adaptive filter
adaptfiltdim = 20;


% % get the true 3D size
% [minvol*prm.voxelvol maxvol*prm.voxelvol] = ...
%     cellsize(prm.diamlow,prm.diamhigh,prm.hx,prm.hz,just,O);

% connectivity
if O == 1
    conn = 4;    
else
    conn = 6;    
end;

% thresholding
thim = logical(adaptfilt(im,adaptfiltdim,prm.adaptfiltth,prm.hx,prm.hz));

test = 0;
if test == 1
    disp('Before remove small')
    showall(im,thim)
    pause
end;


% remove very small
% was at 50
se = strel('disk',2);
thim = imopen(thim,se);
thHere = round(70);
thim = bwareaopen(thim,thHere,conn);

test = 0;
if test == 1
    'after remove small and fix border'
    showall(im,thim)
    pause
end;


%
% find the small regions that can be filled, and keep those as thim
%
% thim = logical(zeros(size(im)));
disp('  Closing and filling');
filled = zeros(size(thim));
step = round(3);
high = round(12);

% make constant around to also include the cells at border. 
% May also create problems in the classification due to half cells 
for i = 1 : O
    thim(:,:,i) = constborder(thim(:,:,i),3,1);
end;
thimOld = thim;

% iteratice closing
for i = 0 : step : high
    
    thim = thimOld;
    if i > 0                
        name = ['ball' int2str(i)];
        load(name);se = getball(ball,i,1);
        thim = imclose(thim,se);
    end;
        
    % fill each region
    [faser,L] = bwlabeln(thim,conn);  
    
    % go back to the old version again
    for j = 1 : L
        regHere = eq(faser,j);    
        
        % reduce to save time
        [regHereRed,box] = boundreg(regHere,1,0);
        [m n o] = size(regHereRed);
    
        
        % fill                
        filledRed = logical(zeros(size(regHereRed)));                
        for  k = 1 : o
            filledRed(:,:,k) = logical(imfill(regHereRed(:,:,k),'holes') - regHereRed(:,:,k));                    
        end;
            
        
        % Remove small and large
        % We need a factor since the filled is inside the border, therefore
        % not the cellsize!!
        filledRed = bwareaopen(filledRed,prm.minvolvox,conn);
        filledRed = filledRed - bwareaopen(filledRed,2*maxvolvox,conn);              
        filledRed = logical(filledRed);
        
        
        % reduced old filled regions
        filledRedOld = filled(box(1,1):box(1,2),box(2,1):box(2,2),box(3,1):box(3,2));
        
%         filledredbef = filledRed;
        filledRed = removeoverlap(filledRed,filledRedOld,conn);
        
        % put back
        filled(box(1,1):box(1,2),box(2,1):box(2,2),box(3,1):box(3,2)) = ...
            gt(filled(box(1,1):box(1,2),box(2,1):box(2,2),box(3,1):box(3,2)) + filledRed,0);       
        
    end;

end;

% the minima becomes the filled regions
minima = filled;

% showall(im,minima)

% and erode!
se = strel('disk',1);
minima = imerode(minima,se);

% MUST open
dim = round(8);
name = ['ball' int2str(dim)];
load(name);se  = getball(ball,dim,1);
minima = imopen(minima,se);

% must remove small after open
minima = bwareaopen(minima,minvol/prm.voxelvol);

% clear border
for i = 1 : O
    minima(:,:,i) = constborder(minima(:,:,i),3,0);
end;
minima = logical(minima);

%
% Close each region separately
%

[faser,L] = bwlabeln(minima);
for j = 1 : L
    reghere = eq(faser,j);
    load ball3;se = getball(ball,3,1);
    reghere = imclose(reghere,se);
    minima(reghere) = 1;
end;


% fill small holes 
minima = fillholes(minima,0,minvol/prm.voxelvol,conn);


% % only keep those that have low intensities,  this i relevant in  sp5
% % where the bordesr are so wide
% [faser,L] = bwlabeln(minima);
% meanim = mean(im(:));
% for j = 1 : L
%     reghere = faser == j;
%     meanint = mean(im(reghere));
%     
%     % if too weak then remove
%     if meanint > 3*meanim
%         minima(reghere) = 0;
%     end;
% end;

if test == 1
    'After closing and filling'
    showall(im,minima)
end

%
% Make background
%

% background is probably lower intensities than MINIMA, use this here
meanminima = mean(im(minima));
numpixminima = length(find(minima));
th = 0.1;
% need to do iterativt to get a proper value, sometimes we have im,
% sometimes minkappa

while 1
    bck = lt(im,th*meanminima);
    numpixbck = length(find(bck));

    if numpixbck > 3*numpixminima || th >= 1
        break;
    end;
    th = th + 0.2;    
end;
% remove the minima from the background not to mix them, also 3D dilation
% here!!
load ball5;se = getball(ball,5,5);
dilminima = logical(imdilate(minima,se));
bck(dilminima == 1) = 0;

% showall(minima,bck)

% must close to get away all small holes
dim = round(5);
se = strel('disk',dim);
bck = imclose(bck,se);

% to remove small parts
bck = imopen(bck,se);

% fill the holes
bck = fillholes(bck,0,prm.minvolvox,conn);

% erode
dim = round(8);
se = strel('disk',dim);
bck = imerode(bck,se);

% dont want small background minima
bck = bwareaopen(bck,round(1.5*maxvolvox),conn);

% add to minima
minima(bck == 1) = 1; 

% 3D again!!
test = 0;
if test == 1  
    'After fill holes in background and clear border'
    showall(im,minima)
end



%---------------------------------------------------------

function [minima] = findminimamda231(im,minkappa,prm)

% % find standard minima
% minima = findminima2d(minkappa,prm);


% and now add the bright nucleus!
[M N O] = size(im);

thim = adaptfilt(im,50,0.1);

% % remove small parts
% vollow = round(pi*prm.diamlow^2);
% th = round(0.1*vollow);
% thim = bwareaopen(thim,th);

% open
load ball1;se = getball(ball,1,1);
thim = imopen(thim,se);

% assign changes
% minima(thim == 1) = 1;
minima = thim;


%----------------------------------------------------

function [minima,prm] = findminimaplane(im,prm)


test = 0;

dim = prm.dim;

% % always!!!
% prm.rangemaxproj = 1;

% dimension of adaptive filter
adaptfiltdim = 20;

% % find midplane automatically
% if ~isfield(prm,'midplane')
%     msg = 'Finding midplane for minima automatically';
%     disp(msg);
%     
%     meanint = zeros(dim(3),1);
%     for i = 1 : dim(3)
%         imhere = im(:,:,i);
%         meanint(i) = mean(imhere(:));
%     end;
% %     [maxint,ind] = max(meanint);
%     sumall = sum(meanint);
%     for i = 1 : 1 : dim(3)
% 
%         % we have reached the "middle" of the cells
%         if sum(meanint(1:i)) > prm.level*sumall;
%             break;
%         end;
%     end;
%     midmin = i;
%     % for saving and info
%     prm.midplane = midmin/dim(3);
% else
%     % the plane where to make makem minima around
%     midmin = max(round(dim(3)*prm.midplane),1);
%     midmin = min(midmin,dim(3));
% end;
% msg = ['Minima plane ' int2str(midmin) ' (out of ' int2str(dim(3)) ')'];
% disp(msg);

% disp(sprintf('  Using planes from %i to %i in max proj for creating minima',low,high));
% imhere = maxprojimage(im(:,:,low:high));
% imhere = im(:,:,max(midmin-1,1):min(midmin+1,dim(3)));
% nplane = size(imhere,3);
% imhere = im(:,:,midmin);
imhere = im(:,:,prm.plane);

% connectivity
conn = 4;    


% thresholding
im = scale(im);
thim = logical(adaptfiltim(imhere,adaptfiltdim,prm.automatedadth));
if test == 1
    disp('Before remove small')
    showall(imhere,thim)
end;

% remove at border
thim = constborder(thim,5,0);

% remove very small
% was at 50
thim = bwareaopen(thim,round(0.2*prm.minvolvox),conn);
% load ball1;se = getball(ball,1,1);
% thim = imopen(thim,se);
if test == 1
    'after remove small and fix border'
    showall(imhere,thim)
end;

% % collect from all three planes
% thim = maxprojimage(thim);
% if dim(3) > 1
%     imhere = imhere(:,:,2);
% end;

%
% find the small regions that can be filled, and keep those as thim
%
% thim = logical(zeros(size(im)));
disp('  Closing and filling');
filled = zeros(size(thim));
step = 3;
high = 10;

% make constant around to also include the cells at border. 
% May also create problems in the classification due to half cells 
% NB this may also create some large minima in small images that make out the whole image!!!!

thim = constborder(thim,5,1);
thimOld = thim;


% iteratice closing
for i = 0 : step : high
    
    thim = thimOld;
    if i > 0                
        name = ['ball' int2str(i)];
        load(name);se = getball(ball,i,1);
        thim = imclose(thim,se);
    end;

%     show(thim,20)

    % fill each region
    [faser,L] = bwlabeln(thim,conn);  
    
    
    % go back to the old version again
    for j = 1 : L
        regHere = eq(faser,j);    
%         [j L]
        % reduce to save time
        [regHereRed,box] = boundreg(regHere,1,0);
        [m n o] = size(regHereRed);
    
%         save test
        
        % fill                
        filledRed = logical(zeros(size(regHereRed)));                
        for  k = 1 : o
            filledRed(:,:,k) = logical(imfill(regHereRed(:,:,k),'holes') - regHereRed(:,:,k));                    
        end;
%         nnz(regHereRed)
%         show(regHereRed,10)
%         show(filledRed,11)
%         pause

        filledRed = bwareaopen(filledRed,prm.minvolvox,conn);
%         show(filledRed,12)
%         prm
        
        filledRed = filledRed - bwareaopen(filledRed,2*prm.maxvolvox,conn);              
        filledRed = logical(filledRed);

        % reduced old filled regions
        filledRedOld = filled(box(1,1):box(1,2),box(2,1):box(2,2),box(3,1):box(3,2));
        
%         filledredbef = filledRed;
        filledRed = removeoverlap(filledRed,filledRedOld,conn);
        
        % put back
        filled(box(1,1):box(1,2),box(2,1):box(2,2),box(3,1):box(3,2)) = ...
            gt(filled(box(1,1):box(1,2),box(2,1):box(2,2),box(3,1):box(3,2)) + filledRed,0);       
        
    end;
%     showall(filled)
end;

% the minima becomes the filled regions
minima = filled;

if test == 1
    'After closing and filling'
    showall(imhere,minima)
end

% and erode!
% se = strel('disk',1);
% minima = erolarreg(minima,3,maxvol);
% minima = imerode(minima,se);

% MUST open large regions!
% 12
dimopen = 10;
name = ['ball' int2str(dimopen)];
load(name);se  = getball(ball,dimopen,1);
[faser,L] = bwlabeln(minima);
minima = zeros(size(minima));
for i = 1 : L
    reghere = faser == i;
%     [range,minZ,maxZ] = bwrange(reghere);
    vol = sum(reghere(:));
    if vol > 2*prm.minvol
        reghere = imopen(reghere,se);
    end;
    minima(reghere == 1) = 1;
end;
% must erode again after opeining to take apart from each other minima we
% have opened
minima = erolarreg(minima,2,prm.minvolvox);

if test == 1
    'After erode and open'
    show(imhere,1)
    show(minima,2)
    pause
end
 
% must remove small after open
minima = bwareaopenrange(minima,round(0.1*prm.minvolvox),conn);

% clear border
minima = logical(constborder(minima,5,0));

%
% Close each region separately
%

[faser,L] = bwlabeln(minima);
for j = 1 : L
    reghere = eq(faser,j);
    load ball3;se = getball(ball,3,1);
    reghere = imclose(reghere,se);
    minima(reghere) = 1;
end;


% fill small holes 
minima = fillholes(minima,0,prm.minvol,conn);

if test == 1
    'After closing each region separately'
    show(imhere,1)
    show(minima,2)
    pause
end


% get background minima
minimabck = getbackground(imhere,minima,prm);

% add to minima
minima(minimabck == 1) = 1; 

% make the minima over many planes if we use only one plane for minima
% NB this must come after the making of background, otherwise the
% background will flood the cells! 
% took this away 20120727
% minima = makeMinimaRange(minima,prm,prm.plane,dim(3));

% make 3D again
m = minima;
minima = zeros(dim);
minima(:,:,prm.plane) = m;

% repeat the plane downward to improve the segmentation
% NB only downward since this is often where the staining is missing
if numel(prm.plane) == 1
    minima = repeatplane(minima,prm.plane,1);
end;

% 3D again!!
if test == 1  
    'After fill holes in background and clear border'
    showall(im,minima)
end

% remove small as final step
minima = bwareaopen(minima,prm.minvolvox);

if test == 1
    'Final'
    showall(im,minima)
end

%--------------------------------------------

function [filledred] = removeoverlap(filledred,filledredold,conn)


[m n o] = size(filledred);

% loop over the reduced filled objects to find overlap, dont want
% overlap
for i = 1 : o
    % new region
    reg = filledred(:,:,i);
    
    % old region
    regold = filledredold(:,:,i);
    
    % count the objects in this plane of new region
    [faser,L] = bwlabeln(reg,conn);
    
    % check the overlap in this plane
    overlap = logical(reg.*regold);
    valoverlap = unique(faser(overlap));
    
    % if no overlap than continue to next object, we do nothing
    if isempty(valoverlap)
        continue;
    end;
        
    % remove the overlap regions
    for j = 1 : length(valoverlap)
        reg(faser == valoverlap(j)) = 0;
    end;
    
    % put back the info
    filledred(:,:,i) = reg;
    
end;

%----------------------------------------------------        
function [minima] =  makeMinimaRange(minimamid,prm,midmin,O)

[M N o] = size(minimamid);

% the true number of planes that are used on each side to make minima
rangeMinHere = round((prm.rangemin*O/2));
disp(sprintf('  Using %i planes on each side to extend minima',rangeMinHere))

% make 3D minima from 2D 
name  = ['ball' int2str(prm.dimminball)];
load(name);se = getball(ball,prm.dimminball,1);
% load ball1; se = getball(ball,1,1);
minima = zeros(M,N,O);
minima(:,:,midmin) = minimamid;

for i = 1 : rangeMinHere
    minimamid = imerode(minimamid,se); 
    try    
        minima(:,:,midmin+i) = minimamid; 
    catch
        disp(sprintf('Unable to make minima in plane %i',midmin + i))        
    end;
    
    try
        minima(:,:,midmin-i) = minimamid; 
    catch
        disp(sprintf('Unable to make minima in plane %i',midmin - i))        
    end;        
end;
minima = minima(:,:,1:O);


%----------------------------------------------------
