function [minima,minimacell,prm] = findminima(varargin)
% FINDMINIMA Find markers for cell segmentation
%
%   FINDMINIMA(IM,PRM) Finds the markers for image IM and settings in
%   struct PRM. Let MPRM = [] for default settings.
%
%   Allowed fields in PRM
%   PRM.ADAPTFILTTH   : Threshold in adaptive filterin. Default 0.01
%   PRM.H             : Voxelsize 1x3 vector. Default [0.5 0.5 1.5]
%   PRM.MINVOLFULL    : Lower cell volume in 1000 mcm3. Default 3.
%   PRM.MAXVOLFULL    : Upper cell volume in 1000 mcm3. Default 50.
%   PRM.LEVEL         : Normalized level 0 -> 1 to make markers. Default is
%                       0.3. If PRM.LEVEL = 'ALL' all planes are used to 
%                       make markers. 
%   PRM.METHOD        : Options
%                       PLANE   -> Uses segmentation image for markers. One 
%                       plane at level LEVEL (Default) or
%                       all planes if PRM.LEVEL = 'ALL'.
%                       NUCLEUS -> Uses nucleus image for markers. Uses all
%                       planes if PRM.LEVEL = 'ALL'.
%   PRM.OTSUTH        : Otsu threshold for segmentation, used in nucleus
%                       markers. Default 1.2.
% 	=======================================================================================
% 	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};
prmin = varargin{2};
imnucl = [];

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

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

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

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

% method
prm.method = 'plane';

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

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

% otsu threshold for segmentaiton
prm.otsuth = 1.2;

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

% % get the adjusted 3D cell volume in voxel size
% prm.minvol = round(prm.minvol*prm.volfactor);
% prm.maxvol = round(prm.maxvol*prm.volfactor);

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

if isequal(prm.level,'all')
    prm.plane = 1 : dim(3);
else
    prm.plane = max(round(prm.level * dim(3)),1);
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);
    
elseif isequal(prm.method,'nucleus')  
    
    msg = ['Using nucleus method for markers'];
    disp(msg);
    
    % add nucleus information    
    minimacell = findminimanuecleus(imnucl,prm);

    % find background
    imhere = scale(imsegm);
    th = 0.01;
    totvol = numel(imsegm);

    mid = round(dim(3)/2);
    while 1
        
        minimabck = lt(imhere,th);
        if dim(3) > 1
            minimabck(:,:,1:mid) = 0;
        end;
        vol = nnz(minimabck);

        if vol/totvol > 0.05
            msg = ['Aborting iterative threshold for background marker'];
            disp(msg);
            break;
        end;
        th = th + 0.02;
    end;
    
    % not the lower planes!!! Then you can get cells into it.
    if dim(3) > 1
        minimabck(:,:,1:mid) = 0;
    end;
    
    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)

        
    load ball2;se = getball(ball,2,2);
    minimabck(imdilate(minimacell,se) == 1) = 0;
    
    % add background and cell
    minima = gt(minimabck + minimacell,0);
    
%     showall(imnucl,minimacell,minimabck)
    
elseif isequal(prm.method,'plane')
%     if strcmp(prm.celltype,'PC12')    
    %     if prm.midplane == Inf
    %         minima = findminima3d(im,prm);
    %     else
    msg = ['Using automated method for markers'];
    disp(msg);

    [minima,prm] = findminimaPC12(imsegm,prm);
    %     end;
    
%     elseif strcmp(prm.celltype,'MDA231')
%         minima = findminimamda231(imsegm,prm);
%     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(imnucl,prm)

vis = 0;

dim = size(imnucl);

% find minima in the nucldus staining using cell segmentation of CT cells!
prm.otsuth = prm.otsuth;
% msg = ['In FINDMINIMANUCLEUS using plane ' int2str(prm.plane)];
% disp(msg);
imhere = imnucl(:,:,prm.plane);
p = 10;
[minimahere,wat] = segmct(imhere,prm.minvolfull/p,prm.maxvolfull);


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;

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


% remove very small ones
% 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] = findminimaPC12(im,prm)


test = 0;

dim = prm.dim;

minvol = prm.minvol;
maxvol = prm.maxvol;

% 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);
% [midmin O]

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

% connectivity
conn = 4;    


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

% remove very small
% was at 50
thim = bwareaopen(thim,round(0.1*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,3,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);
    
        
        % fill                
        filledRed = logical(zeros(size(regHereRed)));                
        for  k = 1 : o
            filledRed(:,:,k) = logical(imfill(regHereRed(:,:,k),'holes') - regHereRed(:,:,k));                    
        end;
            
%         show(filledRed,2)
%         pause
        
        filledRed = bwareaopen(filledRed,prm.minvolvox,conn);
%         show(filledRed,10)
        
        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;
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*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,minvol,conn);

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

%
% Make background
%
% 
% th = 0;

% background is probably lower intensities than MINIMA, use this here
meanminima = mean(imhere(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(imhere,th*meanminima);
    numpixbck = length(find(bck));
    % was at 4!!
    if numpixbck > 4*numpixminima || th >= 2
        break;
    end;
    th = th + 0.1;    
end;

% keep largest
bck = bwkeep(bck,1,conn);

% showall(imhere,background)

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

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

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

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

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

% dont want small background minima
bck = bwkeep(bck,1,conn);

% no overlap to minima in background!
% load ball5;se = getball(ball,5,2);
% [faser,L] = bwlabeln(imdilate(bck,se),conn);
% overlap = faser(minima == 1);
% overlap(overlap == 0) = [];
% overlap = unique(overlap);
% for i = 1 : numel(overlap)
%     bck(faser == overlap(i)) = 0;
% end;
% bck(imdilate(minima,se) == 1) = 0;

% add to minima
minima(bck == 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! 
minima = makeMinimaRange(minima,prm,midmin,dim(3));

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


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