function [cellbw,wat,imsegm,minima,minimacell,info] = segmsurfwat(varargin)
% SEGMSURFWAT 3D segmentation of surface stained cells
%
%   [CELLBW,WAT,IMSEGM,MINIMA,MINIMACELL] = SEGMSURFWAT(IM,MINVOLFULL,MAXVOLFULL) 
%   Finding cells in a surface stained image IM. MINVOLFULL and MAXVOLFULL are the lower 
%   and upper expected volume of a cell in thousand cubic microns.
%   Returning the watershed image WAT, the classified cells in CELLBW, 
%   the ridge image IMRIDGE, the markers MINIMA, information of parameters
%   in INFO and the smoothed image IMSM.
%
%   [CELLBW,WAT,IMSEGM,MINIMA,MINIMACELL,INFO] = SEGMSURFWAT(IM,MINVOL,MAXVOL,'prmfile',PRMFILE)
%   specifies a parameter file that can be used for segmentation
%  
%   [CELLBW,WAT,IMSEGM,MINIMA,MINIMACELL,INFO] = SEGMSURFWAT(IM,MINVOL,MAXVOL,OPTIONNAME,OPTION)
%   specifies an OPTIONNAME for segmentation that can be of one of
%   these options:
%   PRM         : General parameters for SEGMSURFWAT
%   MINIMA      : All markers, must be same size as image
%   SMPRM       : Smoothing parameters
%   CPRM        : Classification parameters
%   MPRM        : Marker parameters
%   IMRIDGE     : Ridge enhanced image
%   MINIMACELL  : Markers for cells, for instance from manual drawings
%   IMNUCLEUS   : Nucleus image for defining markers
%
%   Example 
%   -------
%   
%   A segmentation with cells with volume of between 3000 and 50000 mcm3
%   [cellbw,wat,imsegmout,minima,minimacell,info] = segmsurfwat(im,3,50);
%
%
%   Example 
%   -------
%
%   Specifying a particular method for smoothing
%   smprm.method = 1;
%   [cellbw,wat,imsegmout,minima,minimacell,info] = segmsurfwat(im,5,50,'smprm',smprm);
%
%   Example 
%   -------
%
%   Segmentation in 3D. Using nucleus image for markers and classification.
%
%   load surfstain-and-nucleus-3D.mat
%   % Smoothing
%   smprm.method = 2;
%   % No ridge filtering
%   prm.filterridges = 0;
%   % Subtract the nucleus channel from the surface staining to reduce the
%   % cross talk effect. 
%   imsegm = imsegm - imnucl;
%   [cellbw,wat,imsegmout,minima,minimacell,info] = ...
%    segmsurfwat(imsegm,3,50,'imnucleus',imnucl,'smprm',smprm,'prm',prm);
%   % Show results: image, nucleus image, markers, watershed image, cells
%   showall(imsegm,imnucl,minima,wat,cellbw)
% 
%   Example SURFSTAIN_AND_NUCLEUS_2D
%   -------
%
%   The same as above but in 2D. Using nucleus image for markers and classification.
%
%   % load the data
%   load surfstain_and_nucleus_3D.mat
%   imsegm = imsegm(:,:,5);imnucl = imnucl(:,:,5);
% 
%   % Smoothing
%   smprm.method = 2;
% 
%   % No ridge filtering
%   prm.filterridges = 0;
% 
%   % Subtract the nucleus channel from the surface staining to reduce the
%   % cross talk effect. 
%   imsegm1 = imsegm;
%   imsegm = imsegm1 - imnucl;
%   [cellbw,wat,imsegmout,minima,minimacell,info] = ...
%    segmsurfwat(imsegm,3,80,'imnucleus',imnucl,'smprm',smprm,'prm',prm);
%
%   % Show results: image, nucleus image, markers, watershed image, cells
%   showall(imsegm,imnucl,minima,wat,cellbw)
%
%   See also segmct, classifycells, findminima
% 	=======================================================================================
% 	CELLSEGM: Cell segmentation, Version 2012
% 	Copyright (c): Erlend Hodneland
%   Jonas Lies Vei 91,
%   5009 Bergen
%   Email: erlend.hodneland@biomed.uib.no
% 	=======================================================================================
%

printmsg(['This is ' upper(mfilename) ' for segmentation of cells']);



% image for segmentation
imsegm = varargin{1};

% cell diameters, here after!
minvolfull = 1000*varargin{2};
maxvolfull = 1000*varargin{3};
imsegmini = imsegm;
prmfile = [];
minimacell = [];
imnucl = [];

% for visualilzation    
test = 0;
testfinal = 0;

%
% Standard parameters
%

% voxel size
prm.h = [0.5 0.5 1.5];

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

% filter ridges
prm.filterridges = 1;

% % dilation of strong pixels
% prm.dilation = 0;

% classify cells
prm.classify = 1;


%
% Minima parameters
%
mprm.level = 0.30;

% method
mprm.method = 'automated';

%
% Smoothing parameters
%

% 0 : no filtering
% 1 : coherence enhancing diffusion CED
% 2 : directional coherence enhancing diffusion DCED
smprm.method = 'dirced';

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

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

% GPU in smoothing
prm.gpu = 0;

% dimension image
dim = size(imsegm);
if numel(dim) == 2
    dim = [dim 1];
end;

% 
% Classification parameters
%
cprm.method = 'threshold';


%
% Watershed parameters
%

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

% The signicance of watershed lines, stronger if larger, done on ridgeim!!
% 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;


if maxvolfull < minvolfull
    error('The upper cell volume is lower than the lower, check the input')
end;

% 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};
    
    switch namehere
        case 'minima'
            minima = varhere;
        case('imridge')
            imridge = varhere;
        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;

prm.dim = size(imsegm);

% if you give in a nucleus image you want to use it!
if ~isempty(imnucl)
    mprm.method = 'nucleus';
end;

% test for credibility!
if isequal(mprm.method,'manual')
    if ~exist('minima','var') && ~exist('minimacell','var')
        error('You must specify either minima eller minimacell')
    end;
end;

if isequal(mprm.method,'nucleus')
    if ~exist('imnucl','var')
        error('You must specify a nucleu image');
    end;
end;
%
% Computed settings
%

% get the ADJUSTED 3D cell size
[prm.minvol,prm.minvolvox,prm.maxvol,prm.maxvolvox] = ...
    cellsize(minvolfull,maxvolfull,prm.h,prm.just,dim(3));

cprm.minvolfull = minvolfull/1000;
cprm.maxvolfull = maxvolfull/1000;
mprm.minvolfull = minvolfull/1000;
mprm.maxvolfull = maxvolfull/1000;
wprm.minvolfull = minvolfull/1000;
wprm.maxvolfull = maxvolfull/1000;
mprm.minvolfull = minvolfull/1000;
mprm.maxvolfull = maxvolfull/1000;
prm.minvolfull = minvolfull;
prm.maxvolfull = maxvolfull;

% pixel size
mprm.h = prm.h;
cprm.h = prm.h;
smprm.h = prm.h;
wprm.h = prm.h;


disp('Using settings:');
printstructscreen(prm);

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

%
% Correcting for uneven illumination
%
msg = ['Correcting uneven illumination'];
disp(msg);
se = strel('disk',12);
imsegm = imtophat(imsegm,se);
if isequal(mprm.method,'nucleus')    
    imnucl = imtophat(imnucl,se);
end;

%
% Smoothing
%
smprm.gpu = prm.gpu;
imsegm = smoothim(imsegm,smprm.method,'prm',smprm);

% 
% Ridge enhancement
% 
if prm.filterridges == 1
    disp('Ridge enhancement')
    imsegm = ridgeenhhessian(imsegm,prm.h);           
else    
    msg = 'No ridge enhancement';
    disp(msg);
end;

%
% Finding minima
%

% dont find markers if:
if isequal(mprm.method,'manual') && exist('minima','var')
    disp('Not creating any additional minima, only using manual minima');
else
    % find the minima
    if isequal(mprm.method,'nucleus')
        [minima,minimacell,mprm] = findminima(imsegm,minimacell,mprm,imnucl);
    else        
        [minima,minimacell,mprm] = findminima(imsegm,minimacell,mprm);
    end;
end;

% % add cell information
% if ~isempty(minimacell)
%     minima(minimacell == 1) = 1;
% end;


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

%
% 3D watershed segmentation
%

% segmentation
disp('3D watershed segmentation')
% save test imsegm minima minimacell
% pause
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
%
% save test
if prm.classify
    if ~isempty(minimacell)
        [cellbw,info.classifycells] = classifycells(wat,imsegmini,cprm,minimacell);
    else
        [cellbw,info.classifycells] = classifycells(wat,imsegmini,cprm);
    end;
else
    warning('No cell classification')
    cellbw = zeros(dim);
end;

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

if test == 1 || testfinal == 1
    if dim(3) == 1
        show(imsegm,1);title('Image');
        show(wat,10);title('Wat');
        show(minima,11);title('Minima')
        show(minimacell,12);title('Minimacell')
        show(cellbw,13);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;


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

        