function numCells = countCells(varargin),
%NUMCELLS count number of cells in the image.
%   numCells = countCells(cellImage,cellCenters,cellDimension) counts number
%   of cells in the image. The Algorithm learns about what do you call a
%   cell from the examples of cells in the image. The user provides
%   position of the center for several cells (5-7 is usually enough) and
%   a typical cell dimension (in number of pixels that fits across a cell).
%   The Algorithm identifies and counts the cells via iterative process.
%
%   cellImage   is an image consisting of cells to be counted. cellImage
%               can be image of any class or a file name of the image to be
%               read by imread. See imread.m for file type specifics.
%
%   cellCenters is a 2D numerical array of (x,y) location of the centers of
%               typical cells. The algorithm learns about what do you call
%               a cell based on these control cells. Make sure these cells
%               are visible clearly in the image and do not overlap with
%               other features. cellCenters has size
%               2-by-numberOfControlCells. A typical number of control
%               cells should be 5-7.
%
%   cellDimension  is an avarage number of pixels that fit across a typical
%               cell. 
%
%   Example
%   -------
%   cellImage = 'redBloodCells.jpg';
%   cellCenters = [81 27;63 23;46 80;91 78;94 72;46 80];
%   cellDimension = 7;
%   numCells = countCells(cellImage,cellCenters,cellDimension);
%   % with this image, the numCells should be 2610.
%   See also IMREAD.

%   Copyright 2004 Svetlana V. Panasyuk  

[aMatrix, aCell] = parse_inputs(varargin{:});

[len1 len2] = size(aMatrix);
siz = size(aCell);
num1 = -(siz(1)-3)/2:(siz(1)-3)/2;
num2 = -(siz(2)-3)/2:(siz(2)-3)/2;
cellCount = 0;
x2 = 1;nn = 0;X2 =~ aMatrix;
while sum(x2(:))>0,
    x1 = imfilter(aMatrix.*(aMatrix<-std(aMatrix(:))),aCell,'rep','same');
    x2 = bwmorph(x1>max([max(x1(:))-3*std(x1(:)) 0.25]),'shrink');
    cellCount = cellCount + sum(x2(:));
    X2 = X2 + x2;
    [i j] = find(x2);
    for k = 1:length(i),
        aMatrix(min(max(i(k)+num1,1),len1),min(max(j(k)+num2,1),len2)) = 0;
    end
    aMatrix = aMatrix - mean(aMatrix(:));
    nn(end+1) = sum(x2(:));
end
numCells = sum(sum(X2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Function: parse_inputs
function [aMatrix, aCell] = parse_inputs(varargin),
% Outputs:  aMatrix 2D numerical (double) array same as input image, but
%                       scaled between 0 and 1 
%           aCell   2D numerical (double) array representing an image of a
%                   typical cell and of the same size as a typical 
%                   cell.
%
nargchk(3,3,nargin);

% First, check the image or its file & assign the image array
if ischar(varargin{1}),
    aMatrix = imread(varargin{1});
elseif isnumeric(varargin{1}),
    aMatrix = varargin{1};
else error('first argument must be a filename string or an image array.');
end
aMatrix = im2double(aMatrix(:,:,1));
aMatrix = (aMatrix - min(aMatrix(:)))/(max(aMatrix(:)) - min(aMatrix(:)));
aMatrix = aMatrix - mean(aMatrix(:));

% Second, check & assign the typical cell info
str = {'cellImage','cellCenters','cellDimension'};
for k = 2:3,
    if ~isnumeric([varargin{k}]),
        error(['The ' str{k} ' argument must be numeric.']);
    end
end
siz = size(varargin{2});
if isempty(find(siz==2)),
    error('Second argument must be X-by-2 or 2-by-X array.');
end
cellCenters = permute(varargin{2},[find(siz==2),find(siz~=2)]);
cellRad = ceil((varargin{3}-1)/2);
num = -cellRad:cellRad;
for k = 1:length(cellCenters(1,:)),
    aCell(:,:,k) = aMatrix(cellCenters(1,k)+num,cellCenters(2,k)+num);
end
aCell = mean(aCell,3) - mean(aCell(:));