function [cellbw,infocells] = classifycells(varargin)
% CLASSIFY Classify cells from user defined thresholds
%
%   [CELLBW, INFOCELLS] = CLASSIFY(WAT,IM,PRM) Classify the region in WAT from the image
%   IM. PRM is a struct array of settings used in the classification.
%   CELLBW is the cell image after classification and INFOCELLS contains
%   information about the classification.
%
%   PRM.H             : Voxel size. Default [0.5 0.5 1.5];
%   PRM.METHOD        : Defines the various modes of classification
%
%   1. PRM.METHOD = 'threshold':
%   PRM.MINVOLFULL    : Minimum volume of 3D cell. Default 3.
%   PRM.MAXVOLFULL    : Maximum volume of 3D cell. Default 100.   
%   PRM.INTINCELL     : Intensity in cell. Default 0.7.
%   PRM.INTBORDER     : Intensity of border. Default 1.20.
%   PRM.CONVAREA      : Convex area. Default 0.4. 
%   PRM.CONVPERIM     : Convex perimeter. Default 0.35.
%
%   INTINBORDER and INTINCELL are relative thersholds, compared to the
%   background. A reduced set of thresholds can be defined by a cell array
%   of strings in PRM.PROPNAME. By prm.propname = 'all', all available 
%   thresholds are used.
%
%   2. PRM.METHOD = 'minimacell':
%   The variable MINIMACELL must be defined as
%   [CELLBW, INFOCELLS] = CLASSIFY(WAT,IM,PRM,'MINIMACELL',MINIMACELL) 
%   Classify the region in WAT from the image im using information in
%   MINIMACELL for selecting cells. MINIMIACELL must be a binary image with
%   one white region inside each cell, and otherwise black.
% 
%
%   Ex:
%   cprm.minvolfull = 5;
%   cprm.maxvolfull = 50;
%   [cellbw,infocell] = classifycells(wat,im,cprm);
%
%   See also segmct, classifycells, findminima, selectclassifycells, combinetrainingset
% 	=======================================================================================
% 	CELLSEGM: Cell segmentation, Version 2012
% 	Copyright (c): Erlend Hodneland
%   Jonas Lies Vei 91,
%   5009 Bergen
%   Email: erlend.hodneland@biomed.uib.no
% 	=======================================================================================


msg = ['This is ' upper(mfilename) ' for classification of cells and background'];
disp(msg);

wat = varargin{1};
im = varargin{2};
prmin = varargin{3};
% minimacell = ones(size(im));
if nargin == 4
    % nucleus minima?
    minimacell = varargin{4}; 
end;
dim = size(im);
if numel(dim) == 2
    dim = [dim 1];
end;


%
% Classification parameters
%

% default method
prm.method = 'threshold';

% pixel sizes
prm.h = [0.5 0.5 1.5];

% adjust for 3D
prm.just = 0.9;

% to cut cells
prm.cut = 0;

% minimum and maximum volume
prm.minvolfull = 5;
prm.maxvolfull = 100;

% the intensity threshold relative to the mean of the background
% was at 1.15
prm.intincell = 1.35;

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

% the concavity measure, fully concave at 1, set to 0.4-0.5
prm.convexarea = 0.5;
prm.convexperim = 0.35;

% a cell in a marker image
prm.iscellmarker = 1;

% merge input
prm = mergestruct(prm,prmin);
prm.minvolfull = prm.minvolfull*1000;
prm.maxvolfull = prm.maxvolfull*1000;

% voxel volume
prm.voxelvol = prod(prm.h);


% adjust cell volume
[prm.minvol,prm.minvolvox,prm.maxvol,prm.maxvolvox] = cellsegm.cellsize(prm.minvolfull,prm.maxvolfull,prm.h,prm.just,dim(3));

msg = ['Using settings'];
disp(msg);
printstructscreen(prm);

% the integer values to loop over
valwat = unique(wat(wat > 0));

% if isequal(prm.method,'tree')
%     
%     msg = ['Using classification tree for classification'];
%     disp(msg);
% 
%     msg = ['Loading tree ' prm.treepathname];
%     disp(msg);
%     D = load(prm.treepathname);
%     prm.propname = D.propname;
%     tree = D.tree;    

if isequal(prm.method,'threshold')

    msg = ['Using thresholds for classification'];
    disp(msg);
    
    prm.propname = {'volume',...
                    'volume',...
                    'intincell',...
                    'intborder',...
                    'convexarea',...
                    'convexperim'};
    prm.thname = {'minvol',...
                'maxvol',...
                'intincell',...
                'intborder',...
                'convexarea',...
                'convexperim'};
    prm.logic = {'gt','lt','gt','gt','gt','gt'};
    
    
elseif isequal(prm.method,'minimacell')
    msg = ['Using minimacell and volume for classification'];
    disp(msg);
    prm.propname = {'volume',...                    
                    'volume',...
                    'cellmarker'};
    prm.thname = {'minvol',...
                'maxvol',...
                'iscellmarker'};

    prm.logic = {'gt','lt','eq'};
    
%     nth = 3;
                
%     nth = numel(prm.propname);
                
%     for i = 1 : nth
%         name = prm.name{i};
%         logic = prm.logic{i};
%         val = prm.(name);
%         msg = [makestr(name,p) makestr(logic,p) makestr(num2str(val),p)];
%         disp(msg);
%         % the thresholds used
%         infocells.thused{i,1} = name;
%     end;
    
end;

if isequal(prm.method,'threshold') || isequal(prm.method,'minimacell')
    msg = ['Thresholds used for classification:'];
    disp(msg);
    nth = numel(prm.thname);
    p = 15;
    for i = 1 : nth
        name = prm.propname{i};
        logic = prm.logic{i};
        val = prm.(prm.thname{i});
        msg = [makestr(name,p) makestr(logic,p) makestr(num2str(val),p)];
        disp(msg);
    end;
end;

% number of properties to use
prm.npropname = numel(prm.propname);

%
% Mean int background to adjust the thresholds for intensites
%

% if the image is empty
if isempty(valwat)
    cellbw = zeros(size(wat));
    infocells = [];
    return;
end;

% the number of regions
nwat = length(valwat);

% get properties
prop = cellsgm.cellprop(im,wat,valwat,prm.propname,prm.h);

% classify
cellbw = zeros(dim);
for i = 1 : nwat

    % this value of WAT
    valwathere = valwat(i);
    
    % this region
    reghere = eq(wat,valwathere);
    
%     % get parameters
%     if isequal(prm.method,'tree')
%         val = NaN(1,prm.npropname);
%         for j = 1 : prm.npropname
%             val(j) = prop.(prm.propname{j})(i);
%         end;
%         yfit = treeval(tree,val);
%         iscell = yfit > 0.5;
%         if iscell
%             dec = 1;
%         else
%             dec = 0;
%         end;
%         infocells.iscellhere(i,1) = dec;
%         infocells.iscellfinal(i,1) = dec;
        
    if isequal(prm.method,'threshold') || isequal(prm.method,'minimacell')
            
        if isequal(prm.method,'automated')
            % find overlap to markers
            overlap = minimacell .* reghere;            
            prop.cellmarker(i,1) = ~isempty(find(overlap,1));  
%             prop.cellmarker(i,1)
%             showall(minimacell,reghere)
            
        end;
        
        nth = numel(prm.thname);
        dec = NaN(1,nth);
        for j = 1 : nth
            str1 = prop.(prm.propname{j})(i);      
            str2 = prm.(prm.thname{j});
            arg = [prm.logic{j} '(' num2str(str1) ',' num2str(str2) ')'];        
            dec(1,j) = eval(arg);
        end;
        infocells.iscellhere(i,:) = dec;

        % is it a cell?
        iscell = sum(dec) == nth;

        % store
        infocells.iscellfinal(i,1) = 0;
        if iscell
            infocells.iscellfinal(i,1) = 1;
        end;            
        
    end;
       
    if iscell
        cellbw(reghere) = 1;       
        infocells.iscellfinal(i,1) = 1;
    else
        infocells.iscellfinal(i,1) = 0;
    end;

    if iscell == 1
        res = 'cell';
    else
        res = 'background';
    end;
    disp(sprintf('  Classifying object %i (number %i) out of %i as %s, DEC = %s',i,valwathere,nwat,res,num2str(dec)));       

end;
% to save
infocells.prm = prm;
infocells.prop = prop;
infocells.propname = prm.propname;


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