% FINDEACHCELLSURFTH 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] = findEachCellSurtTh(im,30,100);
% Ex: [cellbw,wat] = findEachCellSurtTh(im,30,100);
%
function [cellbw,wat,minkappa,minima,info] = segmsurfwat(varargin)

printmsg('Segmentation of cells');

% parse inputs

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

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

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.smoothtype = 1;

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

% 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
mprm.diamlow = celldiamlow*0.7;
mprm.diamhigh = celldiamhigh*0.7;

% 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
mprm.adaptfiltth = 0.020;

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

%
% 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.0;

% 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
%

% the intensity threshold relative to the mean of the background
% was at 1.1
cprm.intincellth = 1.15;

% the intensity threshold on border relative to the mean of the image,
% higher more intensities are required
% was at 1.25
cprm.intborderth = 1.25;

% the concavity measure, fully concave at 1, set to 0.4-0.5
cprm.convareath = 0.4;
cprm.convperimth = 0.4;

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

% to cut cells
cprm.disconnPrm = 0;

% cell diameter, must be after loading or we overwrite
cprm.celldiamlow = celldiamlow;
cprm.celldiamhigh = celldiamhigh;

% classify or not
cprm.classifyprm = 1;

%
% Load the parameter file if one is given
%
if ~isempty(prmfile)
    msg = sprintf('Reading from parameter file %s',prmfile);
    disp(msg);
    cmd = [prmfile];
    [prmin,mprmin,cprmin,wprmin] = eval(cmd);
%     cmd = ['[prmin,mprmin,cprmin,wprmin] = ' prmfile];
%     eval(cmd);
    % merge the structs by overriding the existing values
    prm = mergestruct(prm,prmin);
    mprm = mergestruct(mprm,mprmin);
    cprm = mergestruct(cprm,cprmin);
    wprm = mergestruct(wprm,wprmin);
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
        case 'faserminima'
            faserminima = 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);     
        otherwise
            msg = sprintf('Wrong option: %s',namehere);
            error(msg)                
    end;
end;

global M N O test    

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

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


% to visualize during
test = 0;

% to visualize at the end
testFinal = 0;

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



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


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


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

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

% 
% Ridge enhancement
% 

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


%
% Finding minima
%


if ~exist('faserminima','var')
    % find the minima
    minima = findminima(im,mprm,minkappa);
else
    minima = gt(faserminima,0);
    disp('  Not creating minima inside function');
end;

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


%
% 3D watershed segmentation
%

% segmentation
disp('3D watershed segmentation')
% NB,har endret til minkappa!!!, var før im
Iimpose = imimposemin(minkappa,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,im,wprm.thintmerging,wprm.thconvmerging,optlog,optint);
    
    % re segment after fooling around with the image
    disp('  Resegment watershed')
    wat = resegmentwat(wat,im);
    
end;

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


%
% Classify cells
%

if cprm.classifyprm 
    options = [1 2 3 4 5 6]';
    
    cprm.th = [prm.numpixlow  ; 
           prm.numpixhigh ;
           cprm.intincellth ; 
           cprm.intborderth ; 
           cprm.convareath  ; 
           cprm.convperimth];
    % classify      
    [cellbw,info.classifycells] = classifycells(wat,im,cprm.th);
else
    warning('No cell classification')
    cellbw = zeros(size(im));
    info = [];
end;

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

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

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


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

        
% function [wat,diffbinwat] = testwatlines(wat,minkappa,wprm,cprm)
% disp('Testing significance of watershed lines');
% 
% % number of regions
% L = max(wat(:));
%     
% % the watershed lines
% binwat = eq(wat,0);
% binwatnew = binwat;
% 
% % the neighbors of each region
% reg = neighwat(wat);
%     
% checked = [];
% numrem = 0;
% for i = 1 : L
%     % the neighbors
%     neigh = reg(i).val;
%     numneigh = length(neigh);
% 
%     % this region dilated to find the middle watershed line
%     reghere = eq(wat,i);   
%     load ball3;se = getball(ball,3,3);
%     dilreghere = imdilatefast(reghere,se);
% 
%     % mean inside cell
%     meancell = mean(minkappa(reghere));
%     
%     for j = 1 : numneigh
%         
%         % neighbor region
%         neighreg = eq(wat,neigh(j));
%         
%         numpix = length(find(reghere));
%         numpixneigh = length(find(neighreg));
% %         [numpix numpixneigh prm.numpix*cprm.scnumpix]
%         
%         % dont open up to the background
%         if numpix > prm.numpix*cprm.scnumpix || numpixneigh  > prm.numpix*cprm.scnumpix
%             val = [i neigh(j)];
%             checked = [checked ; val];
%             continue;
%         end;
%         
%         ind = [];
%         if ~isempty(checked)
%             % check the earlier ones
% %             checked(:,1) == i 
% %             checked(:,2) == neigh(j)
%             ind = find(checked(:,1) == i & checked(:,2) == neigh(j) ...
%                      | checked(:,1) == neigh(j) & checked(:,2) == i);
%         end;
%         
%         % have checked this before, continue
%         if ~isempty(ind)
%             continue;
%         end;
%         middlewat = binwat .* dilreghere .* imdilatefast(neighreg,se);        
%         middlewat = logical(middlewat);
% 
%         % value inside this neighbor
%         meanneigh = mean(minkappa(neighreg));
%         
% %         % make snake around, must be 3D element!!!
% %         load ball3;se1 = getball(ball,3,6);
% %         load ball2;se2 = getball(ball,2,4);
% %         snake = imdilatefast(middlewat,se1) - imdilatefast(middlewat,se2);
% %         snake(binwat) = 0;
% %         snake = logical(snake);
%         
%         % relative intensities
%         meanwat = mean(minkappa(middlewat)); 
%         relval =  meanwat / mean([meancell  meanneigh]);
% 
% %         numpix
% %         numpixneigh
% %         [i neigh(j)]
% %         meancell
% %         meanneigh
% %         meanwat
% %         checked
% %         relval        
% %         show(maxprojimage(middlewat),3)
% %         show(maxprojimage(wat),4)
% %         pause
%         
%         % test the intensities
%         if relval < wprm.thRelIntWat
%             binwatnew(middlewat) = 0;
%             numrem = numrem + 1;
%         end;
%         val = [i neigh(j)];
%         checked = [checked ; val];
%     end;
%         
%     
% end;
% 
% load ball3;se1 = getball(ball,3,1);
% minima = imerode(imcomplement(binwatnew),se);
% Iimpose = imimposemin(minkappa,minima);
% disp('Repeating watershed segmentation')
% wat = watershed(Iimpose);
% binwatnew = eq(wat,0);
% 
% % the difference
% diffbinwat = binwat;
% diffbinwat(binwatnew == 1) = 0;
% 
% disp(sprintf('Removed %i watershed lines',numrem));
% 
