% CLASSIFY Classify cells from user defined thresholds
% CLASSIFY(WAT,IM,TH) Classify the region in WAT from the image
% IM. TH is an array of thresholds used in the classification, 
% the numbering as below:
% options for classification
% 1. NUMPIXLOW
% 2. NUMPIXHIGH  
% 3. INTINCELL   
% 4. INTBORDER   
% 5. CONVAREA    
% 6. CONVPERIM   
% INTINBORDER and INTINCELL are relative thersholds, compared to the
% background. Use NAN if the classification criteria should not be in use.
%
% Ex:
% th = [20 100 1.1 1.1 0.4 0.4];
% [cellbw,infocell] = classifycells(wat,im,th);
%
function [cellbw,infocells] = classifycells(varargin)

disp('Classification of cells and background')


wat = varargin{1};
im = varargin{2};
thin = varargin{3};
minimacell = ones(size(im));
cellmarker = 0;
if nargin == 4
    % nucleus minima?
    minimacell = varargin{4}; 
    cellmarker = 1;
end;
[M N O] = size(im);


%
% Default settings
%

% % names of options, logical decision to BE a cell!!
% th(1).name = 'numpixlow';  th(1).logic = 'gt';
% th(2).name = 'numpixhigh'; th(2).logic = 'lt';
% th(3).name = 'intincell';  th(3).logic = 'gt';
% th(4).name = 'intborder';  th(4).logic = 'gt';
% th(5).name = 'convarea';   th(5).logic = 'gt';
% th(6).name = 'convperim';  th(6).logic = 'gt';
% th(7).name = 'pyramide';   th(7).logic = 'ge';
% th(8).name = 'cellmarker'; th(8).logic = 'gt';

% number of pixels in cells
celldiamlow = 35;
celldiamhigh = 150;
[th.numpixlow th.numpixhigh] = cellsize(celldiamlow,celldiamhigh,1,1,0.9,O);


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

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

% the concavity measure, fully concave at 1, set to 0.4-0.5
th.convarea = 0.4;
th.convperim = 0.35;

% the fraction of the planes where a region has a pyrmaide like shape with
% fewer pixels as you go upward
th.pyramide = 0;

% is there overlap to nucleus image?
th.cellmarker = 0;

% name of thresholds
th.name = {'numpixlow','numpixhigh','intincell','intborder','convarea','convperim','pyramide','cellmarker'};
th.logic = {'gt','lt','gt','gt','gt','gt','ge','gt'};

% merge input thresholds
th = mergestruct(th,thin);

p = 15;
% number of thresholds
nth = length(th.name);
msg = ['Thresholds used:'];
disp(msg);
for i = 1 : nth
    name = th.name{i};
    logic = th.logic{i};
    val = th.(name);
    msg = [makestr(name,p) makestr(logic,p) makestr(num2str(val),p)];
    disp(msg);
end;

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


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

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

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

%
% Mean int background
%

% get mean intensity of background, therefore get the background as the
% darkest region
% NB Before I used the largest region as the backgroun, what is best?????
numpix = zeros(numwat,1);
meanint = zeros(numwat,1);
for j = 1 : numwat
    regHere = eq(wat,j);
    ind = find(regHere);
    numpix(j) = length(ind);
    meanint(j) = mean(im(ind));
end;
% minimum value 
% [minVal,ind] = max(numpix);
[minVal,ind] = min(meanint);

% mean int background
meanintbck = mean(im(wat == ind));

%
% Classify
%

if cellmarker 
    msg = ['Note that MINIMACELL is checked for volume in case of errors in labeling'];
    disp(msg);
end;

cellbw = logical(zeros(size(im)));
for j = 1 : numwat

    % this value of WAT
    valwathere = valwat(j);
    
    % this region
    reghere = eq(wat,valwathere);

    % we classify by manual labels
	if cellmarker == 1
        val = unique(minimacell(reghere == 1));
        
        % we have a cell        
        infocells(j,1).iscellfinal = 0;
        if sum(val == 1) > 0
            % check the size at least
            vol = nnz(reghere);
            if vol > th.numpixlow && vol < th.numpixhigh
                cellbw(reghere == 1) = 1;
                infocells(j,1).iscellfinal = 1;
            end;
        end;
        continue;
    end;


    % get parameters
    val = getprm(reghere,im,th,meanintbck,minimacell);       
    dec = zeros(nth,1);
    for i = 1 : nth
        thname = th.name{i};
        thval = th.(thname);
        %
        % For saving
        %
        
        % the name
        infocells(j,i).name{i,1} = thname;
        
        % the value
        infocells(j,1).val(i,1) = val(i);          
        
        % the threshold
        infocells(j,1).th(i,1) = thval; 
                        
        % the logical decision
        logic = th.logic{i};
        
        % compare 1
        arg1 = val(i);
                
        % compare 2
        arg2 = thval;
       
        % the command
        arg = [logic '(' num2str(arg1) ',' num2str(arg2) ')'];        
        dec(i) = eval(arg);

        %
        % For saving
        %
                                
        % if this property was ok for being a cell or not
        infocells(j,1).iscellhere(i,1) = dec(i);
                     
    end;
    
    % if ok for all options (!) then a cell
    sumdec = sum(dec);
    iscell = sumdec == nth;
    if iscell
        cellbw(reghere) = 1;       
        infocells(j,1).iscellfinal = 1;
    else
        infocells(j,1).iscellfinal = 0;
    end;

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

end;


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

function [val] = getprm(reghere,im,th,meanintbck,minimacell)

[M N O] = size(im);

% indices of this region
ind = find(reghere);

% size
numpix = length(ind);

% perimeter of region
perim = zeros(size(im));
for i = 1 : O
    perim(:,:,i) = bwperim(reghere(:,:,i));
end;
% I took this away 20111029 since in confocal images the cells get very
% dark
% load ball1;se = getball(ball,1,1);
% perim = imdilate(perim,se);   

% % maximum image for convexity considerations
% regmax = maxprojimage(reghere);

% the range
[range,minz,maxz] = bwrange(reghere);

nth = numel(th.name);
val = zeros(nth,1);
for k = 1 : nth
    
    % find the values of the options specified by the user
    switch(k)
        case(1)
            val(k) = numpix;
        case(2)
            val(k) = numpix;
        case(3)
            % intensities inside
            inside = reghere;
            inside(perim == 1) = 0;
            intincell = mean(im(inside));

            % output
            val(k) = intincell/meanintbck;
            
        case(4)
            % intensities on border       
            intborder = mean(im(perim == 1));

            % output
            val(k) = intborder/meanintbck;
        case(5)


            % convexity 
            convarea = 0;
            for i = minz:maxz
                if ~isempty(find(reghere(:,:,i)))
                    props = regionprops(double(reghere(:,:,i)),'ConvexArea');
                    convarea = convarea + props.ConvexArea;   
                end;
            end;

            % convex area
            convarea = numpix/convarea;

            % output
            val(k) = convarea;
        case(6)
            %
            % convexity of perimeter    
            %

            % middle plane of connection
            mid = round(mean([minz maxz]));
            regmid = reghere(:,:,mid);
            
            [regmax,box] = boundreg(regmid,0,0);
            props = regionprops(double(regmax),'ConvexImage');
            load ball1;se = getball(ball,1,1);
            regmaxconvperim = bwperim(imerode(props.ConvexImage,se));   

            % find perimeter which is convex
            regmaxperim = bwperim(regmax);

            % the convex parts of the boundary
            load ball1;se = getball(ball,1,1);
            regmaxconvperim = imdilate(regmaxconvperim,se) .* regmaxperim;

            % the ratio of convex boundary
            convperim =  length(find(regmaxconvperim)) / length(find(regmaxperim));

            % output
            val(k) = convperim;
        case(7)

            % the increasing cell volume as you go down in stack - or
            % decreasing as you go up                      

            npix = zeros(O,1);
            for i = O : -1 : 1
                reg = reghere(:,:,i);
                npix(i) = sum(reg(:));                
            end;
            dnpix = npix([2:end end]) - npix;
            dnpix = lt(dnpix,0);
            
            % the fraction of how many planes where the cells are
            % decreasing in size as you go up
            [rangereghere,a,b] = bwrange(reghere);
            val(k) = sum(dnpix) / rangereghere;
        case(8)
                        
            % overlap to nucleus?
            reg = reghere .* minimacell;
            val(k) = sum(reg(:));            
            
    end;
end;



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