function [NumCells,meanArea,cellCentroids,cellOrinetations,label_2,output_CLR,BW2,im_adj] = Count_Cells(Image_name,prints)
% Gil Gelbert
% august 2019

% Image_name should be the name of the uint16 grayscale TIF file.
% Image_name should be in the current directory.
% prints is a flag (1 for prints, 0 for no prints. (IMSHOWs)

% System for localizing bacteria in a picture of a petri dish

% Evgeniy Reznik
% august 2019


% Get Image
warning('off','all'); clck = tic;
disp('Loading Image ...');
im_uint32 = imread(Image_name);
im = im2double(im_uint32);
ORG_SIZE = size(im);

conv_kernel = ones(5,5)*(-1);
conv_kernel(3,3) = 24;

microscope_resolution = 13.7; %1 micron is 13.7 pixel
sclae_line_coordinates = [[2285 2559], [1963  1970]];

max_size_param = 4;
cell_MajorAxisLength_min = microscope_resolution*1; % 1 micron long
cell_MinorAxisLength_min = microscope_resolution*0.5; % 1/2 micron wide

cell_MajorAxisLength_max = cell_MajorAxisLength_min*max_size_param;
cell_MinorAxisLength_max = cell_MinorAxisLength_min*max_size_param;

max_cell_area = cell_MinorAxisLength_max*cell_MajorAxisLength_max;
min_cell_area = cell_MinorAxisLength_min*cell_MajorAxisLength_min;
% max_cell_area = 500;
% min_cell_area = 200;

%% counting cells in image
disp('Object segmentation ...');
%adjust contrast
im_adj = imadjust(im);
if (prints==1)
    figure('Name', 'Contrast stretch'); imshow(im_adj);hold on
    title('Contrast stretch')
end

im_adj_gauss = imgaussfilt(im_adj,4);
Csame = conv2(im_adj_gauss,conv_kernel,'same');
% level = adaptthresh(Csame,0.9);
level = graythresh(Csame);
im_bin = imbinarize(Csame, level*0.3); %0.3 1.4
% im_bin = Csame > 0.1;
im_bin_close = imclose(im_bin,strel('disk', 6));
im_rmsmall = bwareaopen(im_bin_close, floor(min_cell_area));
im_remove_ocular =im_rmsmall;
im_remove_ocular(1:2,:) = 0;
im_remove_ocular(end-2:end,:) = 0;
im_remove_ocular(:,1:2) = 0;
im_remove_ocular(:,end-2:end) = 0;
im_filled = imfill(im_remove_ocular, 'holes');
im_filled_med = medfilt2(im_filled);
im_cell_bnd = im_filled_med .* imcomplement(imbinarize(im_adj));

MASK = im_cell_bnd;

% segment separate cells
D = -bwdist(~MASK);
seed_mask = imextendedmin(D,2);
if (prints==1)
    figure; imshowpair(MASK,seed_mask,'blend'); title('Filter out tiny local minima using imextendedmin')
end
D2 = imimposemin(D,seed_mask);
Labc = watershed(D2);
if (prints==1)
    figure; imshow((Labc~=0).*MASK); title('Segment blobs by watershed')
end

% clolorize segmented cells
Lrgb = label2rgb(bwlabel((Labc~=0).*MASK),'jet','w','shuffle');
if (prints==1)
    figure; imshow(Lrgb); title('Colored Watershed Label Matrix')
    figure; imshow(im_adj);
    hold on
    himage = imshow(Lrgb);
    himage.AlphaData = 0.3; title('Colored Labels Superimposed Transparently on Original Image')
end

label_1 = bwlabel((Labc~=0).*MASK);

%% Morphological Filtering [DOES NOTHING NOW]
disp('Morphological filtering ...');

stats = regionprops(label_1,'Area', 'Eccentricity','FilledArea');
area = [stats.Area];
ecc = [stats.Eccentricity];
farea = [stats.FilledArea];
good_regions = find((area > min_cell_area) .* (area < max_cell_area) .* (ecc > 0.6) .* ((farea - area) < 1));
BW2 = ismember(label_1, good_regions);  

if (prints==1)
    figure; imshow(BW2,[]); title('Watershed Label Matrix Morphological filtering'); impixelinfo;
end

%% Get binary image of found objects, no hoels, can be unseparted
disp('Finalizing ...');
label_2 = bwlabel(BW2);
s = rng;
rng(1);
mycolormap = rand(100000,3);
rng(s);
RGB = label2rgb(label_2,mycolormap);
RR=RGB(:,:,1); GG=RGB(:,:,2); BB=RGB(:,:,3);

R = im2uint8(im_adj); R(BW2) = RR(BW2);
G = im2uint8(im_adj); G(BW2) = GG(BW2);
B = im2uint8(im_adj); B(BW2) = BB(BW2);

output_CLR(:,:,1)=R;
output_CLR(:,:,2)=G;
output_CLR(:,:,3)=B;
output_CLR(sclae_line_coordinates(3):sclae_line_coordinates(4), sclae_line_coordinates(1):sclae_line_coordinates(2),:) = 0;
if (prints==1)
    figure; imshow(output_CLR,[]); hold on;
    impixelinfo; title('Colored Labels Superimposed Transparently on Original Image'); hold off;
end

%% summary
disp('Statistics calc ...');

stats = regionprops(label_2,'Area','Centroid','Orientation');
statArea=[stats.Area];
statArea=statArea(statArea>0);
meanArea=round(mean(statArea));
NumCells=size(statArea,2);
allCentroids = [stats.Centroid];
centroidsX = allCentroids(1:2:end-1);
centroidsY = allCentroids(2:2:end); 
cellCentroids = transpose([centroidsX ; centroidsY]);
cellOrinetations = [stats.Orientation];
runtime = round(toc(clck));

fprintf('%d Cells were found.\n',NumCells);
fprintf('The average area of a cell is %d micron^2.\n',round(meanArea/microscope_resolution));
fprintf('Calculation took %d seconds.\n',runtime);

