% CLASSIFY Classify cells from user defined thresholds
% CLASSIFY(WAT,IM,OPTIONS,VAL) Classify the region in WAT from the image
% IM. OPTIONS is an array of which properties shall be used in
% classificatino, the numbering as below:
% options for classification
% NUMPIXLOW     : 1
% NUMPIXHIGH    : 2
% INTINCELL     : 3
% INTBORDER     : 4
% CONVAREA      : 5
% CONVPERIM     : 6
% VAL is an array of same size as OPTIONS with the thresholsd to be used.
% INTINBORDER and INTINCELL are relative thersholds, compared to the
% background
%
% option = [1 2 4];
% val = [100 1000 1.2];
% Ex [cellbw,infocell] = classifycells(wat,im,option,val);
% Ex [cellbw,infocell] = classifycells(wat,im,option,val);
%
function [cellbw,infocells] = classifycells(wat,im,options,val)

disp('Classification of cells and background')

[M N O] = size(im);

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

% the number of test options
numoptions = length(options);
for i = 1 : numoptions    
    
    % this option
    optionhere = options(i);
    
    % this value
    prm(optionhere).th = val(i);        
    
end;

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

% mean of background, find the largest and smallest region

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

% 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;
% [minVal,ind] = max(numpix);
[minVal,ind] = min(meanint);

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


%
% Classify
%


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

    % this value of WAT
    valwathere = valwat(j);
    
    % this region
    reghere = eq(wat,valwathere);
    
    % get parameters
    val = getprm(reghere,im,options,meanintbck);       
%     show(reghere,4)
%     numoptions
%     pause

    dec = zeros(numoptions,1);
    for i = 1 : numoptions
        
        % this option
        optionhere = options(i);
        
        % the logical decision
        logic = prm(optionhere).logic;
        
        % compare 1
        arg1 = val(i);
                
        % compare 2
        arg2 = prm(optionhere).th;
       

        % the command
        arg = [logic '(' num2str(arg1) ',' num2str(arg2) ')'];
        dec(i) = eval(arg);
        
        %
        % For saving
        %
        
        % the threshold
        infocells(j).th(i) = prm(optionhere).th;        
        
        % the name
        infocells(j).name{i} = prm(optionhere).name;
        
        % the value
        infocells(j).val(i) = val(i);
                
        % if this property was ok for being a cell or not
        infocells(j).iscellhere(i) = dec(i);
             
%         show(im,1)
%         show(reghere,2)
%         dec(i)
%         prm(optionhere).name
%         arg1
%         arg2
%         pause
        
    end;
    
    % if ok for all options (!) then a cell
    sumdec = sum(dec);
    iscell = sumdec == numoptions;
    if iscell
        cellbw(reghere) = 1;
        cellbw(reghere) = 1;
        infocells(j).iscellfinal = 1;
    else
        infocells(j).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)));
    
%     options
%     for i = 1 : numoptions
%         i
%         optionhere = options(i);
%         val(i)
%         prm(optionhere).th
%     end;
%     dec    
%     sumdec
%     show(im,1)
%     show(reghere,2)
%     show(cellbw,3)
%     show(wat,4)
%     pause
   

end;


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

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

[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;
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);

numoptions = length(options);
val = zeros(numoptions,1);
for k = 1 : numoptions
    % this option
    optionhere = options(k);

    % find the values of the options specified by the user
    switch(optionhere)
        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;
%             convperim
%             show(reghere,2)
%             show(regmaxconvperim,3)
%             pause
    end;
end;
% show(reghere,10)
% val
% pause

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