% SEGMCELL Cell segmentation
% SEGMCELL(FOLDER,VALSTACKS,PLANES,METHOD,DIAMLOW,DIAMHIGH,SEGMCH,PRMFILE) Segmentation of
% cells from stacks within VALSTACKS, the planes
% specified as an array in PLANES and the segmentation method in METHOD.
% DIAMLOW and DIAMHIGH are the limits for the diameter of a cell the 
% and SEGMCH is the channel to use for segmentation
%
% Choices for METHOD:
% METHOD = 1 : FINDEACHCELLSURFTH Using adaptive thresholding to find minima and 
% segmentation is done on surface stained cells filling of threshold image
% to find
% METHOD = 2 : SEGMCT Finding cells in CT image
% METHOD = 3: SEGMDEBREKMEANS Finding debris in image using k-means
% clustering
% METHOD = 4: SEGMDEBRETH Finding debris in image using adaptive
% thresholdind
%
% SEGMCELL(...'EVAL') adds a parameter for evaluation of the
% segmentation. 
%
% SEGMCELL(...'TIMELAPSE',TIMELAPSE) adds a parameter for a timelapse evaluation 
% where the results in the first frame are used as input for the next frames, 
% which repeats
%
% SEGMCELL(...'MANUALMINIMA') adds a parameter for using manual minima in
% the first frame. Is especially good in combination with TIMELAPSE option
%
% Ex: segmcell('.',1,1:10,1,30,100,1,'prmfilesegmsurfwat')
% Ex: segmcell('.',1,1:100,2,30,100,1,'prmfilesegmct');
%
function [] = segmcell(varargin)

prm.eval = [];
folder = varargin{1};
prm.valstacks = varargin{2};
prm.planes = varargin{3};
prm.method = varargin{4};
prm.diamlow = varargin{5};
prm.diamhigh = varargin{6};
prm.segmch = varargin{7};
prm.prmfile = varargin{8};

% prmin = varargin{8};

prm.timelapse = 0;
prm.manualminima = 0;
% prm.nucleusch = [];
% prm.startsegmplanes = 1;
prm.startplane = 1;
mprm.otsuth = 0.9;
cprm = [];
wprm = [];
smprm = [];
mprm.method = 'plane';

% the number of folders
if iscell(folder)
    numfolder = length(folder);
else
    numfolder = 1;
    a = folder;
    clear folder;
    folder{1} = a;
end;

%
% Load the parameter file if one is given
%
prmfile = prm.prmfile;
if ~isempty(prmfile)
    isfolder = ~isempty(findstr(prmfile,'/'));
    [a b c] = fileparts(prmfile);
    if isfolder        
        A = pwd;
        cd(a);        
    end;    
    msg = sprintf('Reading from parameter file %s',prmfile);
    disp(msg);
    cmd = [b c];
    [prmin,mprmin,cprmin,wprmin,smprmin] = eval(cmd);

    if isfolder
        cd(A);
    end;
    % merge the structs by overriding the existing values
    prm = mergestruct(prm,prmin);
    mprm = mergestruct(mprm,mprmin);
    cprm = mergestruct(cprm,cprmin);
    wprm = mergestruct(wprm,wprmin);
    smprm = mergestruct(smprm,smprmin);
end;


% % make drawings!!!
% if prm.manualminima == 1
%     for i = 1 : numfolder
%         % this folder
%         folderhere = deblank(folder(i,:));
%     
%         % get manual minima if we want them. Do this first otherwise it is
%         % unpratical for amny folders!!!
%         name = [folderhere '/' 'manualminima'];
% 
% %         manualminima = getmanualminima(folderhere,prm);
%         msg = ['Saving manual segmentation of first stack' name];
%         disp(msg);
%         save(name,'manualminima');
%     end; %for
% end;%if

for l = 1 : numfolder

    % this folder
    folderhere = folder{l};

    % load minima if it has been drawn
    name = [folderhere '/' 'manualminima'];    
    if prm.manualminima == 1 || prm.manualminima == 2
        msg = ['Loading manual segmentation ' name];
        disp(msg);
        load(name);
    end;

    
    % in this folder
    evalu = segmentfolder(folderhere,prm,mprm,cprm,wprm,smprm);

    % need this for returning the variable
    prm.eval = evalu;

    if evalu == 2
        disp('Terminating segmentation');
        break;
    end;
    
end;%for

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


function [eval] = segmentfolder(folderhere,prm,mprm,cprm,wprm,smprm)

eval = prm.eval;
% cf = pwd;
% cd(folderhere);

% cellbw = [];
numstacks = length(prm.valstacks);
% loop over all images
for i = 1 : numstacks
    % the stack here
    numStack = prm.valstacks(i);

    name = [folderhere '/'  'stack' int2str(numStack)];    
    
    % get the image
    [im,imsegm,imnucl] = getim(name,prm);

    if isempty(im)
        warning('Could not load image %s, continuing',name)
        continue;
    end;
    
    dim = size(im);

    if dim(1) == 0 || dim(2) == 0 || dim(3) == 0
        msg = ['Image empty one dimension, continuing'];
        disp(msg);
        continue;
    end;

    
    % information to return
    info = [];
    
    try
        
    % segmentation of surface stained cells
    if prm.method == 1       

            

        if isequal(mprm.method,'nucleus')  
            % subtract nucleus from WGA
            b1 = smooth3(imsegm);
            b2 = smooth3(imnucl);

            % for watershed we use subtracted image                
            D = load('ball3.mat');se = getball(D.ball,3,1);
            for j = 0.1:0.1:1
                imsegmtemp = b1 - j*b2;
                vol = imsegmtemp < 0;
                vol = imopen(vol,se);                    
                vol = sum(vol(:));
                % more than just noise?
                if vol > 0.05*numel(imsegm)
                    break;
                end;
                imsegm = imsegmtemp;
            end;
            clear imsegmtemp
            imsegm(imsegm < 0) = 0;

        end;       


        % manual minima given
        if isequal(mprm.method,'manual')
            minimaload = [folderhere '/' 'stack' int2str(numStack) '-mask.mat'];
            D = load(minimaload);
            range = prm.planes(1):min(dim(3),prm.planes(end));
            minima = D.minima(:,:,range,:);
            minimacell = D.minimacell(:,:,range,:);

        end;

        
        % dont take the lowest planes
        if prm.startplane > 1
            msg = ['Segmentation from plane ' int2str(prm.startplane)];
            disp(msg);

            if isequal(mprm.method,'manual')
                % check for manual labels in the lower planes; we use them
                % if they exist!
                lowim = minima(:,:,1:prm.startplane);
                if ~isempty(find(lowim,1))
                    prm.startplane = 1;
                end;
            end;
                
            if isequal(mprm.method,'manual')
                minima = minima(:,:,prm.startplane:end);
                minimacell = minimacell(:,:,prm.startplane:end);
            elseif isequal(mprm.method,'nucleus')
                imnucl = imnucl(:,:,prm.startplane:end);
            end;
            imsegm = imsegm(:,:,prm.startplane:end);
        end;

        
        % segmentation 
        if isequal(mprm.method,'manual')
            [cellbw,wat,minkappa,minima,info.cells] = segmsurfwat(imsegm,prm.diamlow,prm.diamhigh,'prm',prm,'mprm',mprm,'cprm',cprm,'wprm',wprm,'smprm',smprm,'minima',minima,'minimacell',minimacell);                        
            clear temp;
        elseif isequal(mprm.method,'nucleus')
            [cellbw,wat,minkappa,minima,info.cells] = segmsurfwat(imsegm,prm.diamlow,prm.diamhigh,'prm',prm,'mprm',mprm,'cprm',cprm,'wprm',wprm,'smprm',smprm,'imnucleus',imnucl);                                        
        else                  
            [cellbw,wat,minkappa,minima,info.cells] = segmsurfwat(imsegm,prm.diamlow,prm.diamhigh,'prm',prm,'mprm',mprm,'cprm',cprm,'wprm',wprm,'smprm',smprm);                        
        end;

    elseif prm.method == 2
        % CTG staining        
        [cellbw,wat] = segmct(imsegm,prm.diamlow,prm.diamhigh);
        minima = [];minkappa = [];
    elseif prm.method == 3
        % find debre in CT staining
        [cellbw,wat] = segmdebrekmeans(imsegm,prm.diamlow,prm.diamhigh);        
    elseif prm.method == 4
        % find debre in CT staining        
        [cellbw,wat] = segmdebreth(imsegm,prm.diamlow,prm.diamhigh);
        minkappa = [];minima = [];
    elseif prm.method == 5
        [cellbw,wat,minima,info] = segmctgrad(imsegm,prm.diamlow,prm.diamhigh,'prmfile',prm.prmfile);
        minkappa = [];
    else
        error('Wrong option');
    end;
        
    
    if prm.startplane > 1
        msg = ['Adding lower planes'];
        disp(msg);
        
        % fix the lower planes              
        cellbw = repeatplane(cellbw,prm.startplane);
        wat = repeatplane(wat,prm.startplane);                
        m = zeros(dim(1:3));
        m(:,:,prm.startplane:end) = minima;
        minima = m;        
        if isequal(mprm.method,'manual')
            m = zeros(dim(1:3));
            m(:,:,prm.startplane:end) = minimacell;
            minimacell = m;
        end;
    end;
    
    % for evaluating during segmentation
    if eval == 1
        mid = round(dim(3)/2);
        fh = figure(1);

        subplot(1,2,1);imshow(imsegm(:,:,mid),[0 255]);axis off;colormap(gray);title('Raw image')
        subplot(1,2,2);imshow(cellbw(:,:,mid),[0 1]);axis off;colormap(gray);title('Segmentation')
        subplot(1,2,2);imshow(cellbw(:,:,mid),[0 1]);axis off;colormap(gray);title('Segmentation')
        h1 = uicontrol('Style', 'pushbutton', 'String', 'Terminate','Units','Normalized','Position', ...
            [0.8 0.00 0.2 0.05], 'Callback', @terminate_callback);       
        h2 = uicontrol('Style', 'pushbutton', 'String', 'No evaluation','Units','Normalized','Position', ...
            [0.8 0.05 0.2 0.05], 'Callback', @noevaluate_callback);       
        h3 = uicontrol('Style', 'pushbutton', 'String', 'Evaluate more','Units','Normalized','Position', ...
            [0.8 0.10 0.2 0.05], 'Callback', @evaluate_callback);
        setappdata(fh,'eval',eval);


        % wait for input from user
        uiwait(fh)


        % close the window
        if eval == 0 || eval == 2
            close(fh);
        end;

    end

    catch
        wat = zeros(dim);
        cellbw = zeros(dim);
        minima = zeros(dim);
    end;
    
    % no more evaluation, terminate
    if eval == 2
        break;
    end;

    % for saving
    info.diamlow = prm.diamlow;
    info.diamhigh = prm.diamhigh;

    % planes start and stop
    info.plstart = prm.planes(1);
    info.plstop = prm.planes(end);
    info.segmch = prm.segmch;

    % save variables
    name = [folderhere '/' 'im' int2str(numStack) '-segmaut.mat' ];

    ver = version;   
    info.prm = prm;
    info.smprm = smprm;
    info.mprm = mprm;
    info.cprm = cprm;
    info.wprm = wprm;
    wat = double(wat);    
    if strcmp(ver(1:3),'6.5')
        save(name,'cellbw','wat','im','minima','info')    
        disp(sprintf('Saving %s',name))
    else
        save(name,'cellbw','wat','im','minima','info','-V6')    
        disp(sprintf('Saving %s\n',name))
    end;

end%for

% cd(cf);

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

end%segmentfolder

function [imnew] = 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
    imnew(:,:,i) = im(:,:,start);
end;

end





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

function [im,imsegm,imnucl] = getim(name,prm)

try
    load(name);
    disp(sprintf('Segmentation of %s using channel %i',name,prm.segmch));
catch ME1
    disp(sprintf('Unable to load %s',name));
    ME1
    im = [];
    imsegm = [];
    imnucl = [];
    return;
end;

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

% im = I(:,:,prm.planes(1):min(O,prm.planes(end)),:);
im = im(:,:,prm.planes(1):min(O,prm.planes(end)),:);
% % the image to segment
imsegm = im(:,:,:,prm.segmch);

imnucl = [];
if ~isempty(prm.nucleusch)
    imnucl = im(:,:,:,prm.nucleusch);
end;
end%getim

%----------------------------------------
function [manualminima] = getmanualminima(folder,prm)

numfolder = size(folder,1);
for i = 1 : numfolder
    folderhere = deblank(folder(i,:));

    % only the first stack
    stackhere = prm.valstacks(1);
    % get the image to segment
    name = [ folderhere '/' 'stack' int2str(stackhere)];               
    [im,imsegm] = getim(name,prm);

    [M N O] = size(im);
    mid = round(O/2);
    % draw minima on image
    bw = drawimage(imsegm(:,:,mid),1);
    
    %
    % Make background
    %
    bck = bckminima(imsegm,bw);

    % assign to minima
    bw(bck == 1) = 1;

    % assign
    manualminima(:,:,mid,i) = bw;                                                

end %for

end% manualminima

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


function [bck] = bckminima(imsegm,bw)


[M N O] = size(imsegm);
% mid = round(O/2);

if isempty(find(bw,1))
    msg = 'Using default settings in BCKMINIMA in SEGMCELL';
    warning(msg);
    maxminima = 100;
    numpixminima = 1000;
    meanbw = mean(imsegm(:));
else
    [dim faser] = bwsize(bw);
    maxminima = max(dim);
    % background is probably lower intensities than MINIMA, use this here
    meanbw = mean(imsegm(bw == 1));
    numpixminima = length(find(bw == 1));
end;


th = 0.3;
% need to do iterativt to get a proper value, sometimes we have im,
% sometimes minkappa


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

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

% must have no overlap to the given minima
load ball3;se = getball(ball,3,1);
bck(imdilate(bw,se) == 1) = 0;

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


%     % fill the holes
%     minima = fillholes(minima,0,numpixlow,conn);

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

% dont want small background minima
bck = bwareaopen(bck,maxminima,4);
        
end% bckminima

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

end% main function
