% 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(wat,im,thin)

disp('Classification of cells and background')

[M N O] = size(im);


%
% Default settings
%

celldiamlow = 30;
celldiamhigh = 100;

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

% number of pixels in cells
[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.1
th.intincell = 1.15;

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

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

% merge the structs to assign the user defined values
th = mergestruct(th,thin);

% thold = th;
% for i = 1 : numth
%     th(i).val = eval(['thold.' thi(i).name]);   
% end;

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

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


% 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,th,meanintbck);       
%     show(reghere,4)
%     numoptions
%     pause

    dec = zeros(numth,1);
    for i = 1 : length(th)
        %
        % For saving
        %
        
        % the name
        infocells(j,i).name{i,1} = th(i).name;
        
        % the value
        infocells(j,1).val(i,1) = val(i);          
        
        % the threshold
        infocells(j,1).th(i,1) = th(i);        
        
        if isnan(th(i).val)
            % the threshold
            infocells(j,i).th(i,1) = NaN;
            
            % if this property was ok for being a cell or not
            infocells(j,1).iscellhere(i,1) = 1;
            
            continue;
        end;
        
        % the logical decision
        logic = th(i).logic;
        
        % compare 1
        arg1 = val(i);
                
        % compare 2
        arg2 = th(i).val;
       

        % 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 == numth;
    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)

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

numth = length(th);
val = zeros(numth,1);
for k = 1 : numth
    
    % 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;
    end;
end;

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