% SEGMCELL Cell segmentation
% SEGMCELL(FOLDER,VALSTACKS,PLANES,METHOD,DIAMLOW,DIAMHIGH,SEGMCH) 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',EVAL) adds a parameter for evaluation of the
% segmentation. 
%
% SEGMCELL(...'OTSUTH',OTSUTH) specifies the threshold for binarization
% in METHOD = 3.
% 
% SEGMCELL(...'ADTH',ADTH) specifies the threshold for adaptive
% thresholding for METHOD = 4.
% 
% SEGMCELL(...'RAD',RAD) specifies the radius of the adaptive filter for 
% METHOD = 4.

% Ex : segmcell('.',1,'',3,4,20,200,1)
%
function [] = segmcell(varargin)

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

% default values
prm.eval = 0;
prmsegmct.otsuth = 0.7;
prm.adth = 0.1;
prm.rad = 20;
for i = 8:2:nargin
    name = varargin{i};
    switch(name)
        case('eval')
            prm.eval = varargin{i+1};
        case('otsuth')
            prmsegmct.otsuth = varargin{i+1};
        case('adth')
            prm.adth = varargin{i+1};
        case('rad')
            prm.rad = varargin{i+1};
        otherwise
            error('Wront option to SEGMCELL')
    end;
end;

%
% Settings for SEGMCT
%

% default value is 3
prmsegmct.segmopt = 3;

% evalini = eval;
for i = 1 : size(folder,1)
%     % initial value of eval
%     eval = evalini;

    % this folder
    folderhere = deblank(folder(i,:));

    % in this folder
    eval = segmentFolder(folderhere,prm,prmsegmct);

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

    if eval == 2
        disp('Terminating segmentation');
        break;
    end;
    
end;

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

function [eval] = segmentFolder(folderhere,prm,prmsegmct)

%     global eval
    eval = prm.eval;

    % loop over all images
    for i = 1 : size(prm.valstacks,2)
        % the stack here
        numStack = prm.valstacks(i);

        name = [ folderhere '/' 'stack' int2str(numStack)];    

        try
            load(name);
            disp(sprintf('Segmentation of %s using channel %i',name,prm.segmch));
        catch
            disp(sprintf('Unable to load %s',name));
            continue
        end;

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

        if isequal(prm.planes,'all')
            im = I;
        elseif prm.planes(end) <= O 
            im = I(:,:,prm.planes(1):prm.planes(end),:);        
        elseif prm.planes(end) > O
            warning('Given planes extended outside');
            im = I(:,:,prm.planes(1):O,:);                    
        end;
        clear I;
        
        % the image to segment
        imsegm = im(:,:,:,prm.segmch);

        if M == 0 || N == 0 || O == 0
            msg = ['Image empty one dimension, continuing'];
            disp(msg);
            continue;
        end;
        
        % find cells
        info = [];
        if prm.method == 1
            % finding each ccell using adaptive thresholding and filling!
            [cellbw,wat,minkappa,minima,info.cells] = ...
                segmsurfwat(imsegm,prm.diamlow,prm.diamhigh,'intincellth',0.7,'intborderth',0.85);
        elseif prm.method == 2
            % CTG staining
            [cellbw,wat] = segmct(imsegm,prm.diamlow,prm.diamhigh,'otsuth',prmsegmct.otsuth,'segmopt',prmsegmct.segmopt);
            minima = [];minkappa = [];
        elseif prm.method == 3
            % find debre in CT staining
            [cellbw,wat] = segmdebrekmeans(imsegm,prm.diamlow,prm.diamhigh);        
        elseif method == 4
            % find debre in CT staining        
            [cellbw,wat] = segmdebreth(imsegm,prm.diamlow,prm.diamhigh,'adth',prm.adth,'rad',prm.rad,'filteropt','gaussian');
            minkappa = [];minima = [];
        else
            error('Wrong option');
        end;

        % for evaluating during segmentation
        if eval == 1
            mid = round(O/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);
        
%               h1 = uicontrol('Style', 'pushbutton', 'String', 'Terminate','Units','Normalized','Position', ...
%                 [0.8 0.00 0.2 0.05], 'Callback', 'global eval;eval = 2;uiresume;');       
%             h2 = uicontrol('Style', 'pushbutton', 'String', 'No evaluation','Units','Normalized','Position', ...
%                 [0.8 0.05 0.2 0.05], 'Callback', 'global eval;eval = 0;uiresume;');       
%             h3 = uicontrol('Style', 'pushbutton', 'String', 'Evaluate more','Units','Normalized','Position', ...
%                 [0.8 0.1 0.2 0.05], 'Callback', 'global eval;eval = 1;uiresume;');

            % wait for input from user
            uiwait(fh)


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

        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' ];

        ver = version;   
        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;



function noevaluate_callback(hObject,eventdata)
    eval = 0;
    setappdata(fh,'eval',eval);
    uiresume;
end
function evaluate_callback(hObject,eventdata)
    eval = 1;
    setappdata(fh,'eval',eval);
    uiresume;
end
function terminate_callback(hObject,eventdata)
    eval = 2;
    setappdata(fh,'eval',eval);
    uiresume;
end

end% segmentfolder

end %SEGMCELL
