% SEGMSURFWAT Segmentation of surface stained cells
%
% [CELLBW,WAT] = FINDEACHCELLSURFTH(IM,CELLDIAMLOW, CELLDIAMHIGH) Finding cells in 
% a surface stained image. Robust against endocytosis, using adaptive
% filtering instead of filling to define markers. CELLDIAMLOW and
% CELLDIAMHIGH are the minimum and maximum allowed cell diameters
% Returning the watershed image WAT and the classified cells in CELLBW
%
% [CELLBW,WAT] = FINDEACHCELLSURFTH(IM,CELLDIAMLOW, CELLDIAMHIGH,MINIMA)
% MINIMA is here given by the user, as the initialization for the watershed
%
% Ex: [cellbw,wat,minkappa,minima,info] = segmsurfwat(im,30,100);
%
function [cellbw,wat,minkappa,minima,info] = segmsurfwat(varargin)

printmsg('Segmentation of cells');

% parse inputs

% image
imsegm = varargin{1};
% cell diameters
celldiamlow = varargin{2};
celldiamhigh = varargin{3};
prmfile = [];

% for visualilzation    
test = 0;
testfinal = 0;

% initizliations
imnucl = [];

if celldiamhigh < celldiamlow
    error('The upper cell diameter is lower than the lower, check the input')
end;

%
% General parameters, defualt settings
%

% set to 1 if we want ridges filtering
prm.filterridges = 1;

% extra dilation, used for sp5 due to wide borders. If not used, difficult
% to detect TNTs since they are so far from the cells
prm.dilation = 0;

%
% Smoothing parameters, defualt settings
%

% set to 1 if we want pre filtering
smprm.method = 2;

% 2D or 3D smoothing
smprm.smoothdim = 2;

% number of iterations of smoothing
smprm.numit = 3;

% how much in percentage [0 1] of the stack height where to select the middle
% plane for finding minima.
% pixel relations
prm.hx = 1;
prm.hy = 1;
prm.hz = 4;
prm.h = [1 1 4];

% adjust for 3D not complete shape of cells
prm.just = 0.9;

%
% Minima parameters, defualt settings
%

% % midplane for minima selection
% mprm.midplane = 0.4;

% cell diameter
% I change this 20120320 by adding 0.5 in front!
mprm.diamlow = 0.75*celldiamlow;
mprm.diamhigh = 0.75*celldiamhigh;

% pixelsize
mprm.hx = prm.hx;
mprm.hy = prm.hy;
mprm.hz = prm.hz;

% NB may not work if this is too high, flooding goes out
mprm.rangemin = 0.4;

% the threshold for the adaptive filter to find the minima(0.01-0.03)
% wat at 0.008 for the tissue segmentation
% was at 0.015 for the standard 3D PC12
% was at 0.020
mprm.adaptfiltth = 0.01;

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

% method
mprm.method = 'plane';

%
% Watershed parameters
%

% To test intensity of wat lines
wprm.merging = 0;

% The signicance of watershed lines, stronger if larger, done on minkappa!!
% if the raw image, then use much smaller, around 1
% set to 1.05-1.10
% wprm.thint = 1.06;
wprm.thintmerging = 1.3;

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

%
% Classify parameters, defualt settings
%


% pixel sizes
cprm.hx = prm.hx;
cprm.hy = prm.hy;
cprm.hz = prm.hz;

% to cut cells
cprm.disconnPrm = 0;

% classify or not
cprm.classifyprm = 1;

% read from parameter file
for i = 4 : 2 : nargin
    namehere = varargin{4};
    varhere = varargin{5};
    switch(namehere)
        case('prmfile')
            prmfile = varhere;            
        otherwise
            % nothing
    end;
end;

%
% Load the parameter file if one is given
%
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;

%
% The command line has the highest priority, merge the variables here
%
for i = 4 : 2 : nargin
    
    % this variable
    namehere = varargin{i};
    varhere = varargin{i+1};
    if strcmp(varhere,'default')
        continue;
    end;

    switch namehere
        % hvorfor FASERMINIMA her før???
        case 'minima'
            minima = varhere;
        case 'smoothtype'
            smprm.smoothtype = varhere;
        case('minkappa')
            minkappa = varhere;
        case('prmfile')
            % do nothing here!!!
        case('smprm')
            prmin = varhere;
            smprm = mergestruct(smprm,prmin);            
        case('mprm')
            prmin = varhere;
            mprm = mergestruct(mprm,prmin);            
        case('wprm')
            prmin = varhere;
            wprm = mergestruct(wprm,prmin);            
        case('prm')
            prmin = varhere;
            prm = mergestruct(prm,prmin);            
        case('cprm')
            prmin = varhere;
            cprm = mergestruct(cprm,prmin);   
        case('minimacell')
            minimacell = varhere;
        case('imnucleus')
            imnucl = varhere;           
        otherwise
            msg = sprintf('Wrong option: %s',namehere);
            error(msg)                
    end;
end;

global M N O test    

disp('Using settings:');
prm
mprm
wprm
cprm
smprm



% dimension of image
[M N O] = size(imsegm);

% get the 3D cell size
[prm.numpixlow prm.numpixhigh] = cellsize(celldiamlow,celldiamhigh,prm.hx,prm.hz,prm.just,O);

if ~isfield(cprm,'numpixlowth')
    cprm.numpixlow = prm.numpixlow;
end;
if ~isfield(cprm,'numpixhighth')
    cprm.numpixhigh = prm.numpixhigh;
end;

%--------------------------------------------%
%-----------Cell segmentation start----------%
%--------------------------------------------%

%
% Smooth image by coherent anisotropic diffusion for tnt detection
%

imsegmini = imsegm;
% smprm.method = 10;
% types of smoothing
if smprm.method == 1 
    disp(sprintf('Coherence enhancing diffusion'))    
    if smprm.smoothdim == 2
        for i = 1 : O
            imhere = imsegm(:,:,i);
            imsegm(:,:,i) = cohenhdiff(imhere,0.1,smprm.numit,0.0001);       
        end;
    elseif smprm.smoothdim == 3
        % true 3D36
        imsegm = cohenhdiff(imsegm,0.1,smprm.numit,0.0001,'opt','num');       
    else
        error('Wrong option to SEGMSURFWAT')        
    end;
elseif smprm.method == 2
    disp(sprintf('Directional coherence enhancement filtering'))    
    if smprm.smoothdim == 2
        for i = 1 : O
            imhere = imsegm(:,:,i);
            imsegm(:,:,i) = dircohenh(imhere,13,1,1,3);
        end;
    elseif smprm.smoothdim == 3
        % true 3D
         imsegm = dircohenh(imsegm,13,1,1,3);
    else
        error('Wrong option to SEGMSURFWAT')
    end;

else
    
    msg = 'No filtering in SEGMSURFWAT, FILTIM is assigned as IM';
    disp(msg);    
end;


% 
% Ridge enhancement
% 
if prm.filterridges == 1
    disp('Ridge enhancement')
    minkappa = ridgeenhhessian(imsegm,prm.h);    
    
%     minkappa = zeros(size(imsegm));
%     for i = 1 : O
%         disp(sprintf('Plane %i',i))        
%         minkappa(:,:,i) = ridgeenhhessian(imsegm(:,:,i),1,1,2);    
%         stdev = 1;
%         rad = 3;
% %         minkappa(:,:,i) = ridgeenhcurv(imsegm(:,:,i),1,1,3,stdev,rad);    
%     end;
else
    if ~exist('minkappa','var')
        minkappa = imsegm;        
    end;
    msg = 'No ridge enhancement';
    disp(msg);
end;

%
% Finding minima
%


if ~exist('minima','var')
    % find the minima
    if isempty(imnucl)
        [minima,minimacell,mprm] = findminima(imsegm,mprm);
    else
        [minima,minimacell,mprm] = findminima(imsegm,mprm,imnucl);
    end;
else
    disp('  Not creating minima inside function, using manually minima');
end;

% showall(imsegm,minima)
if test == 1
    showall(minkappa,minima)
end;

%
% 3D watershed segmentation
%


% segmentation
disp('3D watershed segmentation')
Iimpose = imimposemin(imsegm,minima);
wat = watershed(Iimpose);clear Iimpose


if wprm.merging == 1
    % Test the significance of watershed lines, the new version
    % optlog = 1 : merge the strong lines (for CT)
    % optlog = 2 : merge the weak lines (for WGA)
    % optint = 1 : use the region based merging
    % optint = 2 : use the snake based merging
    optlog = 2;
    optint = 2;
    disp('  Merging watershed regions')
    wat = mergefragments(wat,imsegm,wprm.thintmerging,wprm.thconvmerging,optlog,optint);
    
    % re segment after fooling around with the image
    disp('  Resegment watershed')
    wat = resegmentwat(wat,imsegm);
    
end;


if test == 1
    show(minima,10)
    show(wat,11)
end;


%
% Classify cells
%

if cprm.classifyprm 
    % classify
    if isfield(cprm,'intincell')
        thin.intincell = cprm.intincell;
    end;
    if isfield(cprm,'intborder')
        thin.intborder = cprm.intborder;
    end;
    if isfield(cprm,'convperim')
        thin.convperim = cprm.convperim;
    end;
    if isfield(cprm,'convarea')
        thin.convarea = cprm.convarea;
    end;
    if isfield(cprm,'numpixlow')
        thin.numpixlow = cprm.numpixlow;
    end;
    if isfield(cprm,'numpixhigh')
        thin.numpixhigh = cprm.numpixhigh;
    end;
    if isfield(cprm,'pyramide')
        thin.pyramide = cprm.pyramide;
    end;
    
    if ~isempty(find(minimacell,1))
        [cellbw,info.classifycells] = classifycells(wat,imsegmini,thin,minimacell);
    else
        [cellbw,info.classifycells] = classifycells(wat,imsegmini,thin);
    end;
else
    warning('No cell classification')
    cellbw = zeros(size(imsegm));
    info = [];
end;

% dilation of strong regions
if prm.dilation
    disp('Dilation of strong pixels')
    cellbw = dilatestrong(cellbw,imsegm);
else
    disp('No dilation of strong pixels')
end;

if test == 1 || testfinal == 1
    if O == 1
        show(imsegm,1);title('Image');
        show(wat,10);title('Wat');
        show(minima,11);title('Minima')
        show(cellbw,12);title('Final segmentation');
        pause
    else
       
        showall(imsegm,cellbw,wat,minima)
    end;
    
end;

% return variables
info.prm = prm;
info.mprm = mprm;
info.wprm = wprm;
info.cprm = cprm;


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

        