% SEGMEACHCELLSNAMES(IM,DIAMLOW,DIAMHIGH) Sedgmentation of image IM, finding
% cells with diameter between DIAMLOW and DIAMHIGH
%
% SEGMEACHCELLSNAMES(IM,DIAMLOW,DIAMHIGH,'MAXNITER',MAXNITER) Running at
% most to MAXNITEr iterations. Default is 500;
%
% SEGMEACHCELLSNAMES(IM,DIAMLOW,DIAMHIGH,'DT',DT) Running at timestep = DT.
% Default is 0.1
%
% SEGMEACHCELLSNAMES(IM,DIAMLOW,DIAMHIGH,'FASERMINIMA',FASERMINIMA) It is possible to
% define the minima used in the calculations. Must be a binary image of
% same size as IM.
%
% Journal of Microscopy, Vol. 201, Pt 3, March 2001, pp. 404�415.
% Received 17 January 2000; accepted 15 November 2000
% Segmentation of nuclei and cells using membrane related
% protein markers
% C. ORTIZ DE SOLORZANO, R. MALLADI, S. A. LELIEVRE & S. J. LOCKETT
%
% Ex: [faser,cellbw] = segmsurfsnakes(im,diamlow,diamhigh,'maxniter',100,'dt',0.1);
% Ex: [faser,cellbw] = segmsurfsnakes(im,diamlow,diamhigh,'maxniter',100,'dt',0.1);
%
function [faserminima,cellbw,infocells] = segmsurfsnakes(varargin)

im = varargin{1};
diamlow = varargin{2};
diamhigh = varargin{3};
maxniter = 100;
dt = 0.01;
% parameter before curvature term
% set to 0.2-0.4
eps = 10;
% parameter before reinforcement term
% set to 5-20
beta = 15;
% paramter in exponent in intensity term
% between [1,10],typically 2-4, stronger if larger, CHECKED!!!
alpha = 5;
% set to 1 for ridges filtering
filterridges = 1;
% set to 0 for no speedup
speedupfact = 0;

for i = 4 : 2 : nargin
    name = varargin{i};
    switch(name)
        case('maxniter')
            maxniter = varargin{i+1};
        case('dt')
            dt = varargin{i+1};
        case('faserminima')
            faserminima = varargin{i+1};
        case('speedupfact')
            speedupfact = varargin{i+1};
        case('filterridges')
            filterridges = varargin{i+1};
        case('eps')
            eps = varargin{i+1};
        case('beta')
            beta = varargin{i+1};
        case('alpha')
            alpha = varargin{i+1};
        otherwise
            error('Wrong option for input');                
    end;
end;
prm.dt = dt;
prm.eps = eps;
prm.beta = beta;
prm.alpha = alpha;
prm.filterridges = filterridges;

% step dimension in x, y and z direction
prm.hx = 1;
prm.hy = 1;
prm.hz = 4;

[M N O] = size(im);
if O == 1
    prm.conn = 8;
else
    prm.conn = 26;
end;

% to get same behaviour for different images
im = scale(im);

%
% Snake parameters
%

% Well working parameters:
% eps = 0
% beta = 5;
% alpha = 5;
% speedupfact = 0.5;

% or for 3D!!
% eps = 5:
% beta = 10;
% alpha = 5;
% speedupfact = 0.4


%
% Technical stuff
%

if sum(le(speedupfact,0) | gt(speedupfact,1)) > 0
    error('Speedupfact must be between [0 1]')
end;

% merging?
prm.merging = 0;

% when to stop
% diffpsith = 0.00001*maxim;
prm.diffbwth = 0;

% set to 0 if no reinitialization
prm.reini = 1;

% how ofter to reinitiliaze psi
prm.reiniit = 101;

% set to 1 if we want plotting of psi 
prm.plotpsi = 0;

% how often to show segmentation
prm.visit = 30;

%
% Merging parameters
%

% The signicance of watershed lines, stronger if larger, done on minkappa!!
% set to 1.05-1.10
% wprm.thint = 1.06;
wprm.thint = 2.6;

% The conexity measure for the merging, we should not allow a very much
% smaller convexity after merging
% ratioconv = 1 if no change in convexity after merging
% ratioconv > 1 if decrease in convexity after merging
% ratioconv < 1 if increase in convexity after merging
% merge for ratioconv < wprm.thconv
wprm.thconv = 1;


%
% Cell classificiation
%

% intensities inside
c.intincell = 0.9;

% intenisties on border
c.intborder = 1.3;

% convex area
c.convarea = 0.5;

% convex border
c.convperim = 0.5;

%
% Filter ridges  
%
if prm.filterridges == 1
    disp('Ridge enhancement')
    minkappa = zeros(size(im));
    for i = 1 : O
        disp(sprintf('Plane %i',i))        
%         minkappa(:,:,i) = filterridgesHessian(im(:,:,i),1,1,2);    
        stdev = 2;
        rad = 5;
        minkappa(:,:,i) = ridgeenhcurv(im(:,:,i),1,1,3,stdev,rad);    
    end;
else
    warning('No ridge enhancement')
    minkappa = im;
end;

% make minima if not exists
if ~exist('faserminima','var')
    % Fields to give:
    % MPRM.ADAPTFILTTH   : Threshold in adaptive filterin
    % MPRM.ADAPTFILTDIM  : Dimension of filter in adaptive filtering
    % MPRM.HX            : pixelsize in xy plane
    % MPRM.HZ            : Pixelsize in z direction
    % MPRM.FACT          : The resize factor to make more robust
    % MPRM.DIAMLOW       : Lower limit for number of pixels in a minima
    % MPRM.DIAMHIGH      : Upper limit for number of pixels in minima
    % MPRM.RANGEMAXPROJ  : Range of max projection for making minima
    % MPRM.MIDPLANE      : The midplane to make the minima around    
    % MPRM.RANGEMIN      : The range from the middle plane to dilate minima
    mprm.adaptfiltth = 0.015;
    mprm.adaptfiltdim = 20;
    mprm.hx = 1;
    mprm.hz = 3;
    mprm.fact = 0.7;
    mprm.diamlow = diamlow;
    mprm.diamhigh = diamhigh;
    mprm.rangemaxproj = 1;
    mprm.midplane = 0.50;
    mprm.rangemin = 0.9;
    mprm.dimminball = 1;
    
    % find minima
    minima = findminima(minkappa,mprm);
    [faserminima,L] = bwlabeln(minima);
end;

% low and high thresholds for cell size
if O < diamlow/2
    disp('Reduced 3D for estimating size')
    minnumpix = round((O^0.95)*pi*(diamlow/2)^2);
    maxnumpix = round((O^0.95)*pi*(diamhigh/2)^2);
else    
    disp('True 3D for estimating size')
    minnumpix = round((4/3)*pi*(diamlow/2)^3);
    maxnumpix = round((4/3)*pi*(diamhigh/2)^3);
end;
prm.minnumpix = minnumpix;
prm.maxnumpix = maxnumpix;


% the number of scales
numscales = max(size(speedupfact));

% minimaini = minima;
% segmenation on a number of scales
for i = 1 : numscales

    % this speedup factor
    prm.speedupfacthere = speedupfact(i);
    
    disp(sprintf('Scale %6.2f\n',prm.speedupfacthere))

    % snake segmentation
    faserminima = snakesegm(faserminima,minkappa,maxniter(i),dt,prm,c,wprm,i);

    if i == numscales
        break;
    end;
    
    %
    % Make minima from wat
    %
    
    % this relation is empirically
    d = 1.0./prm.speedupfacthere;
    disp(sprintf('Erotion factor %6.2f',d))
    
    % minima from wat
    faserminima = wat2minima(faserminima,d,prm);
    
    
%     % remove small
%     minima = bwareaopen(minima,round(minnumpix/2));

end

if prm.merging
    optlog = 2;
    optint = 2;    
    
    % merging
    faserminima = mergefragments(faserminima,im,wprm.thint,wprm.thconv,optlog,optint);

    % relabel to get all numbers in FASER chronologically
    faserminima = relabel(faserminima);
end;

% classify
options = [1 2 3 4 5 6]';
val = [minnumpix maxnumpix c.intincell c.intborder c.convarea c.convperim];
disp('Classify regions')
[cellbw infocells] = classifycells(faserminima,im,options,val);
[cellbw infocells] = classifycells(faserminima,im,options,val);

% for printing
bw = gt(faserminima,0);

% the plane to show
mid = round(O/2);
prm.visim = mid;

printfinal(im,faserminima,bw,cellbw,prm)
printfinal(im,faserminima,bw,cellbw,prm)

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

function [faser] = snakesegm(faserminima,im,maxniter,dt,prm,c,wprm,scalenum)

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

% number of pixels
maxnumpix = prm.maxnumpix;
minnumpix = prm.minnumpix;


%  test for empty minima
if isempty(find(faserminima,1))
    error('The minima image is empty')
end;


% minima = keeplargest(minima,4);
% minima = minima - keeplargest(minima,4);


[Mini Nini Oini] = size(im);
faserminimaold = faserminima;
imold = im;
if ne(prm.speedupfacthere,0)
    disp(sprintf('  Speedup of %6.2f',prm.speedupfacthere'))

    % readjust the pixel size in z to adjust for scaling in xy direction
    prm.hz = prm.hz * prm.speedupfacthere;
    
    minnumpix = round(minnumpix*prm.speedupfacthere^2);
    maxnumpix = round(maxnumpix*prm.speedupfacthere^2);
    [M N O] = size(im);

    % resize    
%     im = imresize3d(im,prm.speedupfact,'cubic');
%     faserminima = imresize3d(faserminima,prm.speedupfact,'nearest');
    for i = 1 : O
        imnew(:,:,i) = imresize(im(:,:,i),prm.speedupfacthere,'bicubic');
        faserminimanew(:,:,i) = imresize(faserminima(:,:,i),prm.speedupfacthere,'nearest');
    end;    
    im = imnew;    
    faserminima = faserminimanew;
    clear imhew faserminimanew

    % to avoid negative image values!
    minim = min(im(:));
    if minim < 0
        im = im - minim;
    end;
end;
[M N O] = size(im);
% whos
% showall(imold,faserminimaold)
% showall(im,faserminima)

% the plane to show
mid = round(O/2);
prm.visim = mid;

% % make sure that the number of minima does not increase after rescaling!!
% [faserold,Lold] = bwlabeln(minimaold,prm.conn);
% [faser,L] = bwlabeln(minima,prm.conn);
% faser = one2one(faserold,faser);
% minima = gt(faser,0);


% ???????????
% remove small artifacts from rescaling
val = unique(faserminima(faserminima > 0));
numpsi = length(val);
for i = 1 : numpsi
    reghere = faserminima == val(i);
    numpix = length(find(reghere));
    if numpix < 0.5*minnumpix
        faserminima(reghere) = 0;
    end;
end;

% minnumpix
% showall(faserminima)

% initiell psi around markers
val = unique(faserminima(faserminima > 0));
numpsi = length(val);
psi = zeros(M,N,O,numpsi);
for i = 1 : numpsi
    reghere = eq(faserminima,val(i));
    psi(:,:,:,i) = reghere;
end
psi = reinipsi(psi);

% initiell segmentation
faser = getsegm(psi);


% smooth image to make stop-function
% wat at 5,3
filt = fspecial('gaussian',3,1);
imfilt = imfilter(im,filt,'replicate');
% function to detect high intensities
power = -prm.alpha*abs(imfilt);
g = exp(power);
if O > 1
    [gy gx gz] = gradient(g,prm.hx,prm.hy,prm.hz);
else
    [gy gx] = gradient(g,prm.hx,prm.hy);
    gz = zeros(size(gx));
end;



% bw = gt(faser,0);
    
% initilaization
niter = 0;
numpixchange = zeros(maxniter,1);
while 1
    niter = niter + 1;        
   
    % update PDE    
    psiold = psi;
    psi = updatepsi(psi,im,gx,gy,gz,g,dt,prm,niter);
        
    % ajdust the psi
    psi = adjpsi(psiold,psi,niter);
    clear psiold;
    
    % reinitialize psi
    if prm.reini && round(niter/prm.reiniit) == niter/prm.reiniit
        psi = reinipsi(psi); 
    end;

    
    % the segmentation from psi
    faserold = faser;
    faser = getsegm(psi);
%     bw = gt(faser,0);
    

%     showall(minima,faser,faserminima,bw)

%     if prm.reini && round(niter/prm.reiniit) == niter/prm.reiniit
%         show(kold,6)
%         show(k,7)
%         show(psiold(:,:,:,1),8)
%         show(psi(:,:,:,1),9)
%         pause
%     end;

%     if niter > 160
%         niter
%         show(psi(:,:,:,1),5)
%         show(faser,6)
%         show(bw,7)
%         pause
%     end;

    % the change in psi to stop iteration
    numpixchange(niter) = getdiff(faserold,faser);       
    clear faserold;
    
    % draw images
    if round(niter/prm.visit) == niter/prm.visit && prm.visit < Inf
        drawim(im,psi,faser,faserminima,prm);        
    end;

    disp(sprintf('  Change of pixels: %i',numpixchange(niter)))
    disp(sprintf('  Number of psi: %i',numpsi));
    disp(sprintf('  Number of iterations: %i, scale number %i',niter,scalenum))

    % test for stop criteria
    p = 20;
    changepixlast = Inf(p,1);
    if niter > p
        changepixlast = numpixchange([niter-p:niter]);
    end;

    if sum(changepixlast) < 3
        disp(sprintf('Terminating iterations since snake does not move more'));
        break;
    elseif niter > maxniter
        disp(sprintf('Terminating iterations due to max number of iterations reached'));
        break;
    end;    
    
    
end;
% draw image again        
drawim(im,psi,faser,faserminima,prm);             

% % also get the rest of the regions as a region, important for classifying when 
% % we find the background value from the largest region!! 
% bck = imcomplement(gt(faser,0));
% maxfaser = max(faser(:));
% faser(bck) = maxfaser + 1;

% resize to normal size again
fasernew = zeros(Mini,Nini,Oini);
for i = 1 : Oini
    fasernew(:,:,i) = imresize(faser(:,:,i),[Mini Nini],'nearest');
end;    
faser = fasernew;
% faser = imresize3d(faser,[Mini Nini Oini],'nearest');

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

function [numpixchange] = getdiff(faserold,faser)

% the watershed values
val = unique(faser(faser > 0));
numval = length(val);

numpixchange = 0;
for i = 1 : numval
    reghere = eq(faser,i);        
    regold = eq(faserold,i);
    diffreg = abs(reghere - regold);

    % number of pixels that have changed
    numpixchange = numpixchange + length(find(diffreg));        

end;

% numpixchange
% show(faserold,4)
% show(faser,5)
% pause

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

function [faser] = relabel(faser)


val = unique(faser(faser > 0));
fasernew = zeros(size(faser));
for i = 1 :length(val)
    fasernew(faser == val(i)) = i;
end;
faser = fasernew;

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


function [] = printfinal(im,faser,bw,cellbw,prm)

im = im(:,:,prm.visim);
faser = faser(:,:,prm.visim);
bw = bw(:,:,prm.visim);
cellbw = cellbw(:,:,prm.visim);
cellbw = cellbw(:,:,prm.visim);

figure(2);
subplot(2,2,1);imagesc(im);colormap(gray);title('Image');
subplot(2,2,2);imagesc(faser);colormap(gray);title('Faser');
overlay(:,:,1) = imoverlay(im,bw,0.7);
overlay(:,:,2) = im;
overlay(:,:,3) = im;
subplot(2,2,3);imagesc(overlay);title('Overlay');
subplot(2,2,4);imagesc(cellbw);colormap(gray);title('Cells');
subplot(2,2,4);imagesc(cellbw);colormap(gray);title('Cells');

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


function [psi] = adjpsi(psiold,psi,niter)

[M N O numpsi] = size(psi);
psinew = zeros(M,N,O,numpsi);

if numpsi == 1
    psinew = psi;
    return;
end;

if niter/2 == round(niter/2)
    start = 1;
    step = 1;
    stop = numpsi;
else
    start = numpsi;
    step = -1;
    stop = 1;
end;

% avoid collisions
for i = start : step : stop
    
    % adjust this psi
    psihere = psi(:,:,:,i);
    
    psioldhere = -Inf(M,N,O);
    for j = start : step : stop
        
        % dont check against itself
        if i == j
            continue;
        end;

        % the psi to check against
        psioldhere = max(psioldhere,psiold(:,:,:,j));
        
    end;
    
    % adjust
    psinew(:,:,:,i) = min(psihere,-psioldhere);                

%     showall(psiold(:,:,:,i),psinew(:,:,:,i),psinew(:,:,:,i)-psiold(:,:,:,i))
end;



% to return
psi = psinew;

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

function [faser] = getsegm(psi)

[M N O numpsi] = size(psi);

faser = zeros(M,N,O);
for i = 1 : numpsi
    
    bwhere = gt(psi(:,:,:,i),0);
    faser(bwhere == 1) = i;
%     i
%     whos
%     showall(psi(:,:,:,i),faser,bwhere,faserminima)

    %     if niter > 150
%         i
%         numpsi
%         show(psi(:,:,:,i),5)
%         show(bwhere,6)
%         show(faser,7)
%         pause
%     end;
%     bw = gt(bw + bwhere,0)  ;
%     showall(psi(:,:,i),bwhere)
end;


% showall(faser)

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

function [] = drawim(im,psi,faser,faserminima,prm)
imini = im;
faserini = faser;
faserminimaini = faserminima;
psiini = psi;

% make one plane
im = im(:,:,prm.visim);
faser = faser(:,:,prm.visim);
faserminima = faserminima(:,:,prm.visim);
psi = psi(:,:,prm.visim,:);


overlay = imoverlayrgb(im,faser,1);
if prm.plotpsi
    figure(1);
    [M N O P] = size(psi);
    dimx = P+1;
    count = 0;
    for i = 1 : P
        count = count + 1;    
        msg = sprintf('Minima %i',i);
        subplot(dimx,3,count);imagesc(faserminima == i);colormap(gray);axis off;title(msg);    
        count = count + 1;    
        msg = sprintf('Levelset %i',i);
        subplot(dimx,3,count);imagesc(psi(:,:,:,i));colormap(gray);axis off;title(msg);    
        count = count + 1;
        msg = sprintf('Segmentation levelset %i',i);    
        subplot(dimx,3,count);imagesc(faser == i);colormap(gray);axis off;title(msg);        
    end;
    count = count + 1;
    subplot(dimx,3,count);imagesc(overlay);colormap(gray);axis off;title('Overlay segmentation and image')        
    count = count + 1;
    subplot(dimx,3,count);imagesc(faser);colormap(gray);axis off;title('Segmentation all')
end;
% cl
% showall(imini,faserini,faserminimaini)
figure(3);imagesc(overlay);colormap(gray);axis off;title('Overlay');drawnow
figure(4);imagesc(faser);colormap(gray);axis off;title('Faser');drawnow
figure(5);imagesc(faserminima);colormap(gray);axis off;title('Faserminima');drawnow
% pause

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

function [psi] = reinipsi(psi)


% filt = fspecial('gaussian',5,3);
disp('  Reinitialize');

[M N O P] = size(psi);
psiold = psi;
for i = 1 : P
    % this psi
    psihere = psi(:,:,:,i);
%     psihereold = psihere;
    
    % binary to compute the distance function
    bw =  gt(psihere,0);

    % the border
    perim = bwperim(bw);
    % need this, otherwise we get a border at image border, a wrong psi and
    % a strange H (Curvature) which spreads and creates instabilities
    for j = 1 : O
        perim(:,:,j) = constborder(perim(:,:,j),3,0);
    end;
    
    
    % distance functions
    psi1 = bwdist(bw);
    psi2 = bwdist(imcomplement(bw) + perim);
    
%     show(imcomplement(bw) + perim,10)
%     % get the border to fine adjust the distance function
%     for j = 1 : O
%         bwhere = bw(:,:,j);
%         perimhere = bwperim(bwhere);
%         perimhere = constborder(perimhere,3,0);
%         psi1(:,:,j) = bwdist(bwhere);
%         psi2(:,:,j) = bwdist(imcomplement(bwhere) + perimhere);
%     end;

    % add both
    psihere = -psi1 + psi2;
    
    % adjust for emtpy image
%     psihere(isinf(psihere)) = 0;

    

% %     smooth a little??
%     psi(:,:,:,i) = imfilter(psihere,filt,'replicate');
    psi(:,:,:,i) = psihere + 1;
    
%     showall(-psi1,psi2,psi(:,:,:,i),bw,gt(psi(:,:,:,i),0))
    
%     show(gt(psi(:,:,:,i),0),10)
%     show(gt(psihereold,0),11);
%     show(psi(:,:,:,i),12);
%     show(psihereold,13)
%     pause
end;

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


function [psinew] = updatepsi(psi,im,gx,gy,gz,g,dt,prm,niter)

[M N O P] = size(psi);
hx = prm.hx;
hy = prm.hy;
hz = prm.hz;


psitemp = zeros(size(psi));
for i = 1 : P
    % this function
    psihere = psi(:,:,:,i);
    
%     prm.eps = 5;
    % curvature term
    H = curv(psihere,hx,hy,hz);
    
    %
    % Fra Ian Mitchell: http://www.cs.ubc.ca/~mitchell/ToolboxLS/index.html
    %
    % lagt på desktop og path endret, men må skrive ved hver oppstart
    % 'addpath(genpath(a))' der 'a = pwd' i toolboxen
%     grid = processGrid([M N]', psihere);
%     [ H, gradMag ] = curvatureSecond(grid, psihere);
    
    if O > 1
        [psiherey psiherex psiherez] = gradient(psihere,hx,hy,hz);
    else
        [psiherey psiherex] = gradient(psihere,hx,hy);    
        psiherez = zeros(size(psiherey));
    end;
    
    % absolute gradient
    absgradpsi = sqrt(psiherex.^2 + psiherey.^2 + psiherez.^2);    

    % update the PDE
    update = g.*(-1 - prm.eps*H).*absgradpsi - prm.beta * (gx.*psiherex + gy.*psiherey + gz.*psiherez);   
%     update = g.*(-1 - prm.eps*H).*absgradpsi;
    psitemp(:,:,:,i) = psihere - dt*update;

%     showall(psitemp(:,:,:,i),H)    
end;

% return variable
psinew = psitemp;


% showall(im,g,gx,gy,gz)
% if niter > 1000
%     showall(im,psi,absgradpsi,g,gx,gy,gz,psiherex,psiherey,psiherez)
% end;

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

% function [minima] = one2wat(wat,minima,prm);
% 
% [faser,L] = bwlabeln(minima,prm.conn);
% 
% val = unique(wat(wat > 0));
% numval = length(val);
% for i = 1 : numval
%     
%     valhere = val(i);
%     reghere = wat == valhere;
%     
%     
%     
% end;

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

function [watnew] = wat2minima(wat,r,prm)
val = unique(wat(wat > 0));
numval = length(val);
r = round(r);
se = strel('disk',r);
watnew = zeros(size(wat));

for i = 1 : numval
    valhere = val(i);
    
    % this region
    reghere = eq(wat,valhere);
    
    % erode region
    reghere = imerode(reghere,se);
    
    % the new watershed
    watnew(reghere) = valhere;
    
end;

% dont want the number of regions to increase due to erotion!!!!!!
watnew = one2one(wat,watnew);

% cl
% showall(wat,watnew)

% 
% minima = gt(watnew,0);

% showall(wat,watnew)


return;

% [faserminima,L] = bwlabeln(minima,prm.conn);
% r = round(r);
% se = strel('disk',r);
% 
% for i = 1 : L
%     reghere = eq(faserminima,i);
%     
%     % the overlap values to the watershed
%     val = wat(reghere);
%     
%     % how many of each of the unique?
%     [numele,B] = numunique(val);
%     
%     [maxval,I] = max(numele);
%     
%     % the watershed region with the maximum overlap maximum 
%     valhere = B(I);
%     
%     % make the minima from the watershed region of most overlap to the
%     % minima, and erode it a little
%     regminima = wat == valhere;
% 
%     % want smaller region
%     regminima = imerode(regminima,se);
%     
%     % only largest since the erotion may create more regions
%     regminima = keeplargest(regminima,prm.conn);
%     minima(regminima) = 1;
%     
% end;

% showall(minima,wat)

% val = unique(wat(wat > 0));
% numval = length(val);
% minima = zeros(size(wat));
% r = round(r);
% se = strel('disk',r);
% for i = 1 : numval
%     reghere = eq(wat,val(i));    
% 
%     
%     
%     reghere = imerode(reghere,se);
%     reghere = keeplargest(reghere,prm.conn);
%     minima(reghere == 1) = 1;
% 
% end;
% % cl
% val
% showall(minima,wat)
%---------------------------------------------



% function [minima] = makeminima(im,minnumpix,prm)
% 
% test = 0;
% 
% [M N O] = size(im);
% 
% % fill image
% for i = 1 : O
%     filled(:,:,i) = imfill(im(:,:,i),8,'holes');
% end;
% 
% [fx fy] = gradient(filled);
% absgrad = sqrt(fx.^2 + fy.^2);
% 
% % the minima
% minima = eq(absgrad,0);
% if test == 1
%     'Gradient == 0'
%     showall(im,minima);
% end;
% 
% meanminima = mean(im(minima == 1));
% bck = lt(im,meanminima*0.6);
% minima = gt(minima + bck,0);
% if test == 1
%     'After adding background'
%     showall(im,minima);
% end;
% 
% % remove small
% minima = bwareaopen(minima,round(0.5*minnumpix),prm.conn);
% if test == 1
%     'After reemove small'
%     showall(im,minima);
% end;
% 
% 
% % close
% load ball2;se = getball(ball,2,1);
% minima = imclose(minima,se);
% % showall(minima)
% % % must be removed!!!
% % minima = keeplargest(minima,conn);
% % showall(im,filled,absgrad,minima)
% if test == 1
%     'Final'
%     showall(im,minima);
% end;

