% FINDMINIMA Find the minima meanto prepare for segmemaxntation
% FINDMINIMA(MINKAPPA,MPRM) Finds the minima from 
% Fields to give:
% MPRM.ADAPTFILTTH   : Threshold in adaptive filterin
% MPRM.HX            : pixelsize in xy plane
% MPRM.HZ            : Pixelsize in z direction
% MPRM.DIAMLOW       : Lower limit for number of pixels in a minima
% MPRM.DIAMHIGH      : Upper limit for number of pixels in minima
% MPRM.MIDPLANE      : The midplane to make the minima around
% MPRM.RANGEMIN      : The range from the middle plane to dilate minima
% 
function [minima,minimacell,mprm] = findminima(varargin)
% IM is always IM and MINKAPPA can be either ridge enhanced or IM
disp('Finding minima');

global test

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

% default values
mprm.adaptfiltth = 0.01;
mprm.hx = 1;
mprm.hz = 4;
mprm.diamlow = 30;
mprm.diamhigh = 130;
% mprm.midplane = 0.4;
mprm.rangemin = 0.5;
mprm.dimminball = 1;
mprm.celltype = 'PC12';
mprm.level = 0.5;
mprm.method = 'plane';
mprm.thresholdlevel = 0.05;
mprm.level = 0.3;

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

% assign the input variables
if ~isempty(mprmin)
    mprm = mergestruct(mprm,mprmin);
end;


if isequal(mprm.level,'all')
    mprm.plane = 1 : dim(3);
else
    mprm.plane = round(mprm.level * dim(3));
%     mprm.plane = [max(mprm.plane-1,1) : min(mprm.plane+1,dim(3))];
end;


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


if isequal(mprm.method,'threshold')
    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,mprm.thresholdlevel);
    end;
    load ball2;se = getball(ball,2,2);
    mask = imclose(mask,se);
    mask = imopen(mask,se);
    
    minima = mask;
    just = 0.9;
    [numpixlow numpixhigh] = cellsize(mprm.diamlow,mprm.diamhigh,mprm.hx,mprm.hz,just,dim(3));
    minima = bwareaopen(minima,numpixlow);
    
elseif isequal(mprm.method,'nucleus')  
    
    
    % find background
    msg = ['Find background marker'];
    disp(msg);
    imhere = scale(imsegm);
    th = 0.01;
    totvol = numel(imsegm);
    dim = size(imsegm);
    mid = round(dim(3)/2);
    while 1
        
        minimabck = lt(imhere,th);
        minimabck(:,:,1:mid) = 0;
        vol = sum(minimabck(:));
        if vol/totvol > 0.05
            break;
        end;
        th = th + 0.01;
    end;
    
    % not the lower planes!!! Then you can get cells into it.
    minimabck(:,:,1:mid) = 0;
    
    minimabck = imfill(minimabck,'holes');
    load ball2;se = getball(ball,2,1);
    minimabck = imopen(minimabck,se);    
    minimabck = bwkeep(minimabck,3,6);
    minimabck = imclose(minimabck,se);    
%     minimabck = bwareaopen(minimabck,
%     showall(imhere,minimabck)

    % add nucleus information    
    minimacell = findminimanuecleus(imsegm,imnucl,mprm);
        
    load ball2;se = getball(ball,2,2);
    minimabck(imdilate(minimacell,se) == 1) = 0;
    
    % add background and cell
    minima = gt(minimabck + minimacell,0);
    
    
    
elseif isequal(mprm.method,'plane')
    if strcmp(mprm.celltype,'PC12')    
    %     if mprm.midplane == Inf
    %         minima = findminima3d(im,mprm);
    %     else
            [minima,mprm] = findminimaPC12(imsegm,mprm);
            minimacell = zeros(dim);
    %     end;
    
    elseif strcmp(mprm.celltype,'MDA231')
        minima = findminimamda231(imsegm,mprm);
    end;
end;

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

function [im] = repeatplane(im,start)

% dim = size(im);
% add = start-1;
% dimnew = dim;
% dimnew(3) = dimnew(3) + add;
% imnew = zeros(dimnew);
% imnew(:,:,start:end) = im;
for i = 1 : start - 1
    im(:,:,i) = im(:,:,start);
end;


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

function [minima] = findminimanuecleus(imsegm,imnucl,mprm)


dim = size(imnucl);

just = 0.9;
[numpixlow numpixhigh] = ...
    cellsize(mprm.diamlow,mprm.diamhigh,mprm.hx,mprm.hz,just,dim(3));

% find minima in the nucldus staining using cell segmentation of CT cells!
prm.otsuth = mprm.otsuth;
% msg = ['In FINDMINIMANUCLEUS using plane ' int2str(mprm.plane)];
% disp(msg);
imhere = imnucl(:,:,mprm.plane);
[minimahere,wat] = segmct(imhere,0.5*mprm.diamlow,mprm.diamhigh,'prm',prm);

% 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(:,:,mprm.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);

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

% showall(minima)
% repeat the plane downward to improve the segmentation
if numel(mprm.plane) == 1
    minima = repeatplane(minima,mprm.plane);
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;



% load ball1;se = getball(ball,1,1);
% minima = imerode(minima,se);
minima = bwareaopen(minima,round(numpixlow/(dim(3)^just)));


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

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;




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

function [ minima ] = findminimafill(im,mprm)


dim = size(im);

% find cell size
just = 0.9;
[numpixlow numpixhigh] = cellsize(mprm.diamlow,mprm.diamhigh,mprm.hx,mprm.hz,just,dim(3));
numpixlowplane = round(numpixlow/(dim(3))^just)
numpixhighplane = round(numpixhigh/(dim(3))^just)

im = scale(im);
meanim = mean(im(:));
maxim = max(im(:));

load ball5;se5 = getball(ball,5,1);
load ball3;se3 = getball(ball,3,1);
load ball2;se2 = getball(ball,2,1);
for i = dim(3) : -1 : 1
    imhere = im(:,:,i);
    thim = gt(imhere,0.3);
    thim = bwareaopen(thim,numpixlowplane,4);
    thim = constborder(thim,2,1);
    filled = imfill(thim,'holes') - thim;
    filled = imerode(filled,se2);
    filled = imopen(filled,se3);
    filled = bwareaopen(filled,numpixlowplane,4);
%     filled = filled - bwareaopen(filled,numpixhighplane,4);
    
    % remove high intensity regions
    a = filled;
    [faser,L] = bwlabeln(filled);
    for j = 1 : L
        reghere = faser == j;
        val = mean(imhere(reghere == 1));
        if val > 2*meanim
            filled(reghere) = 0;
        end;        
    end;
    
            
    % check overlap
    wat = bwlabeln(filled);
    if i < dim(3)
        n = max(wat(:));
        removed = zeros(dim);
        for j = 1 : n
            reghere = wat == j;
%             vol = sum(reghere(:));
            val = unique(watold(reghere));
            val(val == 0) = [];
            nval = numel(val);
            volold = zeros(numel(val));
            for k = 1 : nval
                regold = watold == val(k);
                volold(k) = sum(regold(:));
            end;
            meanvolold = mean(volold);
            
            % two regions go into one for a small region --> a cell, that
            % is not ok
            if numel(val) > 1 & meanvolold < numpixhighplane
                filled(reghere) = 0;
                removed(reghere) = 1;
                'removed'
            end;
%             show(imhere,1);title('imhere');
%             show(watold,2);title('watold')
%             show(wat,3);title('wat')
%             show(reghere,4);title('reghere')
%             pause
            
        end;

        show(imhere,1);title('Image');
        show(filled,2);title('Filled');
        show(removed,3);title('Removed');
        
    showall(imhere,filled,removed)    
    end;
    watold = wat;
    minima(:,:,i) = filled;
                
end;
showall(im,minima)

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


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

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
[numpixlow numpixhigh] = ...
    cellsize(mprm.diamlow,mprm.diamhigh,mprm.hx,mprm.hz,just,O);

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

% thresholding
thim = logical(adaptfilt(im,adaptfiltdim,mprm.adaptfiltth,mprm.hx,mprm.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,numpixlow,conn);
        filledRed = filledRed - bwareaopen(filledRed,2*numpixhigh,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,numpixlow);

% 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,numpixlow,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,numpixlow,conn);

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

% dont want small background minima
bck = bwareaopen(bck,round(1.5*numpixhigh),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,mprm)

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


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

thim = adaptfilt(im,50,0.1);

% % remove small parts
% vollow = round(pi*mprm.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,mprm] = findminimaPC12(im,mprm)


% always!!!
mprm.rangemaxproj = 1;


[M N O] = size(im);

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

% dimension of adaptive filter
adaptfiltdim = 20;

% number of pixels in plane
numpixlow = round(pi*(mprm.diamlow/2)^2);
numpixhigh = round(pi*(mprm.diamhigh/2)^2);

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

% % midmin = max(1,round(O*mprm.midplane));
% % midmin = min(midmin,O);
% I took this away, do we need it?
% low = max(round(midmin)-mprm.rangemaxproj,1);
% high = min(round(midmin)+mprm.rangemaxproj,O);
% 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,O));
% nplane = size(imhere,3);

% connectivity
conn = 4;    

% test = 1;
test = 0;

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


% remove very small
% was at 50
thim = bwareaopen(thim,70,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 O > 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,3,1);
thimOld = thim;
% show(thim,2)
% pause

% 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);  
    
%     show(thim,3)
%     pause
    
    % 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);
    
        
        % fill                
        filledRed = logical(zeros(size(regHereRed)));                
        for  k = 1 : o
            filledRed(:,:,k) = logical(imfill(regHereRed(:,:,k),'holes') - regHereRed(:,:,k));                    
        end;
            
%         show(filledRed,9)
        % Remove small and large
        % We need a factor since the filled is inside the border, therefore
        % not the cellsize!!

        filledRed = bwareaopen(filledRed,numpixlow,conn);
%         show(filledRed,10)
        
        filledRed = filledRed - bwareaopen(filledRed,2*numpixhigh,conn);              
        filledRed = logical(filledRed);
%         show(filledRed,11)
%         pause

        % 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;

% if isnucleus == 1
%     msg = ['Using nucleus staining for finding minima'];
%     disp(msg);
%     
%     % find minima in the nucldus staining using cell segmentation of CT cells!
%     [miniman,watn] = segmct(imnucleus(:,:,midmin),mprm.diamlow,mprm.diamhigh);
%     
%     load ball3;se = getball(ball,3,1);
%     miniman = imerode(miniman,se);
%         
%     % add both
%     minima = minima == 1 | miniman == 1;
%     
%     % make 3D for later use!!
%     a = zeros(size(im));
%     a(:,:,midmin) = miniman;
%     miniman = a;    
% else
%     miniman = ones(size(im));
% end;
% minima = thim;

% show(filled,10)
% pause
if test == 1
    'After closing and filling'
    showall(imhere,minima)
end

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

% MUST open large regions!
dim = 12;
name = ['ball' int2str(dim)];
load(name);se  = getball(ball,dim,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*numpixlow
        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,numpixlow);

if test == 1
    'After erode and open'
    show(imhere,1)
    show(minima,2)
    pause
end
 
% must remove small after open
minima = bwareaopenrange(minima,100,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,numpixlow,conn);

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

%
% Make background
%

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


while 1
    bck = lt(imhere,th*meanminima);
    numpixbck = length(find(bck));
    % was at 4!!
    if numpixbck > 4*numpixminima || th >= 2
        break;
    end;
    th = th + 0.1;    
end;

% trengs denne her????
load ball2;se = getball(ball,2,1);
bck(imdilate(minima,se)) = 0;

if test == 1
    'Background'
    show(im,1)
    show(bck,2)
    show(minima,3)
    pause
end;
% save test
% must close to get away all small holes
dim = 5;
se = strel('disk',dim);
bck = imclose(bck,se);

% to remove small parts
se = strel('disk',dim+2);
bck = imopen(bck,se);

% fill the holes
minima = fillholes(minima,0,numpixlow,conn);
% filled = imfill(bck,'holes');
% holes = filled - bck;
% holes = holes - bwareaopen(holes,numpixlowplane);
% bck(holes == 1) = 1;

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

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

% showall(bck,imhere,minima)

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

% % make a littel smaller
% load ball1;se = getball(ball,1,1);
% minima = imerode(minima,se);

% whos minima bck
% pause

% showall(minima,imhere)


% 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! 
minima = makeMinimaRange(minima,mprm,midmin,O);

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

% % full size again
% minimafull = zeros(M,N,O);
% for i = 1 : size(minima,3)
%     minimafull(:,:,i) = imresize(minima(:,:,i),[M N],'nearest');
% end;
% minima = minimafull;


% remove small as final step
minima = bwareaopen(minima,numpixlow);


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,mprm,midmin,O)

[M N o] = size(minimamid);

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

% make 3D minima from 2D 
name  = ['ball' int2str(mprm.dimminball)];
load(name);se = getball(ball,mprm.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);


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