% SEGMCT(IM,DIANLOW,DIAMHIGH) Segmentation of celltracker green image which has
% to be in the first dimension of the image function I which is loaded from
% [STACK NUM TYPE].  Takes the planes from PLSTART to PLSTOP. FOLDER is the
% folder where the segmentation is stored.
% Saves the segmentation and the images, and returns the
% segmentation CELLBW. RAD is estimated diameter of cells.
%
% SEGMCT(IM,DIAMLOW,DIAMHIGH,'OTSUTH',otsuth) defines the threshold in
% Otsus method. Its only used in the wahly approach.
%
function [cellBW,wat] = segmct(varargin)

vis = 0;

im = varargin{1};
diamlow = varargin{2};
diamhigh = varargin{3};
otsuth = 1.0;
otsuth = 2;
for i = 4 : 2 : nargin
    varhere = varargin{i};
    switch(varhere)
        case('otsuth')
            otsuth = varargin{i+1};
        otherwise
            error('')
    end;
end;

disp('Segmentation by SEGMCT.M');
[M N O] = size(im);

if O == 1 
    conn = 8;
else
    conn = 26;
end;

% parameters specific to the run
opDiam = 5;
stdVal = 3;

% expected number of pixels in a cell per plane
rad = diamlow / 2;
thpixlow = round(pi*rad*rad*O^0.9);
rad = diamhigh / 2;
thpixhigh = round(pi*rad*rad*O^0.9);

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


% filter image
filtopt = 0;


% for uint8!!
if max(im(:)) > 256
    im = im * 256/4095;
end;

imini = im;

% tophat image to remove uneven illumination
se = strel('disk',50);
im = imtophat(im,se);
% showall(imini,im)

disp('Filter image using Gaussian');
im = imfilter(im,fspecial('gaussian',7,stdVal),'replicate');


if filtopt == 1
    disp('Filter image using anisotrpic diffusion filtering');    
    for i = 1 : O
        disp(sprintf('Plane %i',i))
        im(:,:,i) = edgeenhdiff( im(:,:,i), 0.3, 20, 0.05);
    end;
end;

% showall(im,imini)
% pause


segmopt = 1;
if segmopt == 1
    % Almost 
    % High througput analysis of multispectral images of breast cancer
    % tissue, Umesh Adiga, IEEE, Vol8, No 8, 2006
    [cellBW,wat] = segmmethod1(im,thpixlow,thpixhigh,conn,opDiam,otsuth);
elseif segmopt == 2
    % Almost 
    % Combining intensity, edge and shape information for 2d and 3D
    % segmentaton fo cell nuclei in tissue sections, Wahlby, Journal of
    % microscopy, Vol 215, pp 67-76, 2004
    [cellBW,wat] = segmmethod2(im,thpixlow,thpixhigh,conn,opDiam,otsuth);
end;


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

function [cellBW,wat] = segmmethod2(im,thpixlow,thpixhigh,conn,opDiam,otsuth)

[M N O] = size(im);

% filt = fspecial('gaussian',13,7);
% im = imfilter(im,filt,'replicate');

test = 0;

% 
% Make the segmentaiton essemtially, foreground seeds
%

% whos
% pause
% must have this for the uint8 thing!!
im = 255*scale(im);
im = uint8(im);
% NB, lavere gir høyere!!!
cellBW = zeros(size(im));
level = otsuth*graythresh(im);
for i = 1 : O
    cellBW(:,:,i) = im2bw(im(:,:,i),level);
end;
im = double(im);

if test
    'graythresh'
    showall(im,cellBW)
end;

% remove small parts
name = ['ball' int2str(opDiam)];
load(name);se = getball(ball,opDiam,1);
cellBW = imopen(cellBW,se);

% remove small parts
disp('Remove small and big parts')
cellBW = bwareaopen(cellBW,thpixlow,conn);    
% cellBW = cellBW - bwareaopen(cellBW,thpixhigh,conn);    
% showall(im,cellBW)
if test
    'small'
    showall(im,cellBW)
end;

dist = bwdist(imcomplement(cellBW));   
filt = fspecial('gaussian',13,9);
dist = imfilter(dist,filt,'replicate');
dist = imcomplement(dist);
if test
    'dist'
    showall(im,cellBW,dist)
end;

% make minima here!!!
minima = imregionalmin(dist,conn);
load  ball4;se = getball(ball,4,O);
minima = imdilate(minima,se);
minima = imclose(minima,se);
minima(cellBW == 0) = 1;

% make the "ridge"
dist(cellBW == 0) = min(dist(:));
% for the VECT
dist = scale(dist);

% % make a "fake" 3D minima!!!
% mid = round(O/2);
% if round(O/2) == O/2
%     vect = [mid:-1:1 1:1:mid];
% else
%     vect = [mid:-1:2 1 2:1:mid];
% end;
% vect = vect.^2;
% vect = 0.02*scale(vect);
% % dist = scale(dist);
% % vect = vect * 0.05;
% disp(sprintf('Weighting the planes to create a 3D minima using vect = %s',num2str(vect)));
% % add the weights
% for i = 1 : O
%     dist(:,:,i) = dist(:,:,i) + vect(i);
% end;

% %
% % To combine minima we make markers!!
% %
% minima = imregionalmin(dist,conn);
% % minima = eroLarMin(minima
% load  ball5;se = getball(ball,5,O);
% minima = imdilate(minima,se);
% minima = imclose(minima,se);

% showall(im,dist,minima)

% this watershed captures the boundaries found in graythresh, but also
% divides the connected cells!!!
Iimpose = imimposemin(dist,minima);
wat = watershed(Iimpose);
% save test dist wat im
% showall(im,dist,minima,wat)
if test
    'wat'
    showall(im,minima,wat)
end;


cellBWold = cellBW;
watold = wat;
%
% Merge
%
% if we make it smaller we merge weaker boundaries, i.e. if too little
% merging make it smaller, if too much merging (we loose tru regions to
% background) make it bigger
thint = 0.86;
% smaller than zero means that the region must become more convex after
% merging
thconv = 0.87;
optlog = 1;
optint = 1;
wat = mergefragments(wat,im,thint,thconv,optlog,optint);

% re segment after fooling around with the image
wat = resegmentwat(wat,dist);

% make cells from watershed
disp('Classify cells')
cellBW = classify(wat,im,thpixlow,thpixhigh);

% showall(im,cellBW,dist,minima,wat,cellBWold,watold)
% pause

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

function cellBW = classify(wat,im,thpixlow,thpixhigh)

meanall = mean(im(:));
L = max(wat(:));

cellBW = zeros(size(im));
for i = 1 : L
    reghere = eq(wat,i);
   
    % mean intensities in this region
    meanint = mean(im(reghere));
    
    % the reatio to the rest of the image
    ratio = meanint/meanall;
    
    % number of pixels
    numpix = length(find(reghere));
    
%     show(im,1)
%     show(reghere,2)
%     ratio

%     pause
    % was at 1.2
    th = 1.2;
    if ratio > th && numpix > thpixlow && numpix < thpixhigh
        cellBW(reghere) = 1;
    else        
        disp('Not cell,')
        disp(sprintf('Ratio intensity: %6.2f [ ''>'' %6.2f ]',ratio,th))
        disp(sprintf('Number pixels: %6.2f [%6.2f  ''>'' %6.2f ]',numpix,thpixhigh,thpixlow))
        
    end;
end;



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

function [cellBW,wat] = segmmethod1(im,thpixlow,thpixhigh,conn,opDiam,otsuth)

[M N O] = size(im);
im = 255*scale(im);
im = uint8(im);
% NB, lavere gir høyere!!!
cellBW = zeros(size(im));
% th = 1.0;
level = otsuth*graythresh(im);

for i = 1 : O
    cellBW(:,:,i) = im2bw(im(:,:,i),level);
end;
im = double(im);
% showall(im,cellBW)



% remove small parts
name = ['ball' int2str(opDiam)];
load(name);se = getball(ball,opDiam,1);
cellBW = imopen(cellBW,se);

% % remove small parts
% disp('Remove small and big parts')
% cellBW = bwareaopen(cellBW,thpixlow,conn);    
% % cellBW = cellBW - bwareaopen(cellBW,thpixhigh,conn);    
% showall(cellBW)

%
% Close objects
%
disp('Close objects')
[faser,L] = bwlabeln(cellBW,6);
load ball4;se = getball(ball,4,1);
for i = 1 : L
    regHere = eq(faser,i);
    regHere = imclose(regHere,se);
    cellBW(regHere) = 1;
end;

%
% Fill holes
%
disp('Fill holes')
[faser,L] = bwlabeln(cellBW,6);
load ball4;se = getball(ball,4,1);
for i = 1 : L
    regHere = eq(faser,i);
    fill = imfill(regHere,'holes') - regHere;
    fill = fill - bwareaopenrange(fill,thpixlow,conn);
    cellBW(fill == 1) = 1;
end;

vis = 0;
if vis == 1
    disp('AFter filling of holes and before distance function')
    showall(im,cellBW);
end;

%
% Make distance function to divide the cells
%
disp('Make distance function to cut')
dist = bwdist(imcomplement(cellBW));
dist = imfilter(dist,fspecial('gaussian',20,5),'replicate');    


% find maximum points
maxBW = imregionalmax(dist,conn);

% close them to avoid over segmentaiton
load ball3;se = getball(ball,3,1);
maxBW = imdilate(maxBW,se);
% % NB must be in 3D!!!!????? (only if 3D)
load ball7;se = getball(ball,7,O);
maxBW = imclose(maxBW,se);

if vis == 1
    disp('The maximum image MAXBW')
    showall(im,maxBW);
end;

% dilate strong pixels without merging wit neighbours
% larger: less dilation
intth = 4;
% intth = 1.5;

% intth = 0.4;
disp('Dilate the strong pixels without merging with the neighbours')
% this method is dilating one pixel at a time, dont think it makes a big
% differnece to dilatestrong
% cellBW = dilatestrong2(maxBW,im,intth,conn);
% use this if other code
cellBW = dilatestrong(cellBW,maxBW,im,intth,conn);

if vis == 1
    disp('After dilation')
    showall(im,cellBW);
end;

% remove small parts
disp('Remove small and big parts')
cellBW = bwareaopen(cellBW,thpixlow,conn);    
cellBW = cellBW - bwareaopenrange(cellBW,thpixhigh,conn);    


test = 0;
if test == 1    
    showall(im,cellBW)
end;
 
% label for output
[wat,L] = bwlabeln(cellBW);

if vis == 1
    disp('Final output')
    showall(im,cellBW);
end;
% ----------------------------------------------
% this method is dilating one pixel at a time, dont think it makes a big
% differnece to dilatestrong
function [cellbw] = dilatestrong2(bw,im,intth,conn)

imini = im;
load ball2;se = getball(ball,2,1);
niter = 0;
while 1
    niter = niter + 1;
    bwold = bw;
    
    % label the binary regions
    [faser,L] = bwlabeln(bw,conn);
%     showall(im,bwold,bw,faser)
        
    % dilated border
%     border = zeros(size(im));
    faserth = zeros(L,1);
    for i = 1 : L
        reghere = faser == i;
        
        % the dilation
        borderhere = imdilate(reghere,se) - reghere;
        
%         % update BORDER to use it later 
%         border(borderhere == 1) = 1;
        
        % update the FASERBORDER to know which FASER the dilated pixel
        % belongs to
        faser(borderhere == 1) = i;
        
        % the threshold for this FASER
        faserth(i) = mean(im(reghere));
    end;
    
%     showall(im,bw,faser)
    % the border to sort the values
    border = gt(faser,0) - bw;

    % sort the pixels
    ind = find(border == 1);
    val = im(ind);
    [sorted,indsorted] = sort(val,1,'descend');
    
    % loop over the sorted pixels to start wiht the smallest one!!
    for i = 1 : length(val)
        indarray = indsorted(i);
        indimage = ind(indarray);
        
        % this image value
        valhere = val(indarray);
        
        % this faser value to find the connection to FASER to check for
        % threshold value
        faservalhere = faser(indimage);

        % the threshold for this FASER
        thhere = faserth(faservalhere)*intth;
        
        % neighbour values to look for crash!!
        point = zeros(size(im));
        point(indimage) = 1;
        point = imdilatefast(point,ones(5,5));
        faservalaround = faser(point == 1);
        
        % take away itself and 0, that is no crash
        faservalaround(faservalaround == faservalhere) = [];
        faservalaround(faservalaround == 0) = [];

%         if indimage == 38279
%             faserhere = faser;
%             im(indimage) = max(im(:));
%             faserhere(indimage) = max(faser(:)) + 2;
%             valhere
%             thhere
%             indimage
%             showall(im,faserhere)
%             im = imini;
%         end;
        
        % too close neighbour region
        if length(faservalaround)> 0
%             disp('Too close neighbour, continue')
            continue;
        % too weak, we stop since the array is sorted
        elseif valhere < thhere
            disp('Too weak, break');
            
            break;            
        end;
        
        % assign the value
        bw(indimage) = 1;
        
    end;
    
%     faserth
%     showall(im,bwold,bw,faser)
    
    % stop criterion
    diffbw = bw-bwold;
    numchange = sum(diffbw(:));
    disp(sprintf('Number of dilations: %i, changing pixels: %i',niter,numchange))
    if numchange == 0
        break;
    end;    
end;
% output
cellbw = bw;
cellbw = bw;

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

function [cellBW] = dilatestrong(cellBW,BW,im,intth,conn)

[M N O] = size(im);

[faser,numcells] = bwlabeln(BW);
cellBWini = cellBW;
[fasercellBW,LcellBW] = bwlabeln(cellBW);

%
% Get the means to have a stop criteria
%
borderth = zeros(M,N,O);
% calculate distance function to divide the image into zones
for i = 1 : O
    cellBWred(:,:,i) = imresize(cellBW(:,:,i),0.3,'nearest');
end;
boundimred = distim(cellBWred,conn);
boundim = zeros(M,N,O);
for i = 1 : O
    boundim(:,:,i) = imresize(boundimred(:,:,i),[M N],'nearest');
end;
% showall(boundim)


borderth = ones(size(im))*mean(im(:));

% to be similar as below
cellBW = gt(faser,0);   
% radius of dilation between each step
rad = 2;
name = ['ball' int2str(rad)];
load(name);se = getball(ball,rad,O);
niter = 0;
while 1
    niter = niter + 1;
    
    % dilate  to find the border of the cell
    % NOTE: BORDERFASER is the border of the cells, not the image!!!!
    borderfaser = zeros(M,N,O);
    for i = 1 : numcells 
        reghere = eq(faser,i);
        borderfaser((imdilatefast(reghere,se) - reghere) == 1) = i;
        % update the thresholds
        borderth(reghere) = mean(im(reghere));
    end;
    % all border after dilation
    border = gt(borderfaser,0);
%     show(im,4)
%     show(borderth,3)

    % remove the pixels with too low intensities
    border = border.*ge(im,intth*borderth);
    % update borderfaser also after removing pixels
    borderfaser = borderfaser.*border;

    
    %
    % Now we want to find out if the regions are close to other regions
    %
    
    % do filtering to check neighborship for another number of FASERBORDER
    d = rad*3;
    for i = 1 : O
        filtmax(:,:,i) = ordfilt2(borderfaser(:,:,i),d*d,ones(d,d));
    end;
    filtmax = filtmax .* border;
    borderfaserold = borderfaser;

    showall(borderfaser,filtmax)
    
    % to fix the zero problem gets the min
    borderfaser(borderfaser == 0) = Inf;
    for i = 1 : O
        filtmin(:,:,i) = ordfilt2(borderfaser(:,:,i),1,ones(d,d));
    end;
    % must remove the Inf again since the ordfilt will make Inf in
    % homogenoius regions
    filtmin(filtmin == Inf) = 0;
    filtmin = filtmin .* border;
    borderfaser = borderfaserold;
    
%     showall(filtmax,filtmin,faser,border)
    
    % the ones to remove from wrong neighborship
    borderfaser(filtmax ~= borderfaser) = 0;
    borderfaser(filtmin ~= borderfaser) = 0;
    
    % add the ok pixels to the faser image
    faser = faser + borderfaser;
    
    % the cells as binary
    cellBWold = cellBW;
    cellBW = gt(faser,0);
    
    diffBW = cellBW-cellBWold;
    numchange = sum(diffBW(:));
    disp(sprintf('Number of dilations: %i, changing pixels: %i',niter,numchange))
    if numchange == 0
        break;
    end;
%     show(cellBW,6)
%     pause

end


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

function [cellBW] = convexmerge(cellBW)

[M N O] = size(cellBW);

[faser,L] = bwlabeln(cellBW);

for i = 1 : L
    reghere = eq(faser,i);
    
    % get the neighbours
    load ball3;se = getball(ball,3,1);
    dilreghere = imdilate(reghere,se);
    
    neigh = unique(faser(dilreghere == 1));
    neigh(neigh == i) = [];
    neigh(neigh == 0) = [];
    numneigh = length(neigh);
    
    
    % loop over the neighbrous to test convexity
    ratioboth = zeros(numneigh,1);
    for j = 1 : numneigh
        j
        
        % this neighbour
        neighhere = neigh(j);                
        regneigh = eq(faser,neighhere);

        % the combination together
        load ball3;se = getball(ball,3,1);
        regboth = imclose(regneigh + reghere,se);
        
        % the convex are of the combination of this region and the
        % neighbour
        convareaboth = calcconvarea(regboth);        
        areaboth = length(find(regboth));
        
        % the ratio
        ratioboth(j) = areaboth / convareaboth; 
        
    end;
                
    % convex area of thie region
    convareareghere = calcconvarea(reghere);
    areahere = length(find(reghere));
    ratiohere = areahere / convareareghere;
    
    
   numneigh
    show(reghere,2)
    show(faser,3)
    neigh
    ratioboth
    ratiohere
    pause
    
    
end;

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

function [convarea] = calcconvarea(BW)

[M N O] = size(BW);
convarea = 0;
for j = 1 : O
    props = regionprops(double(BW(:,:,j)),'ConvexArea');
    convarea = convarea + props.ConvexArea;        
end;
