% Gil Gelbert
% August 2019
% System for localizing bacteria in a picture of a petri dish

% Evgeniy Reznik
% January 2020

% main function.
warning('off','all'); close all; clc; clear;

%%get files
yourfolder = '';
segmentation_product_folder = '\Segmented';
mkdir(fullfile(yourfolder, segmentation_product_folder));
graphs_folder = '\Graphs';
mkdir(fullfile(yourfolder, graphs_folder));
files_BF = dir(fullfile(yourfolder, '*c2*.tif'));
files_GFP = dir(fullfile(yourfolder, '*c1*.tif'));
image_time_period = 6; %period of time between images
gfp_frequency = 4; %parameter that represnts frequency of gfp pictures
bf_frequency = 1; %parameter that represnts frequency of bf pictures
mycolormap_C = rand(10000,3);
linestyle = ["-";"--";"-.";];
microscope_resolution = 13.7; %1 micron is 13.7 pixel
bin_step = 100;
histBinLocations = 200:bin_step:800;

%% segment cells
for x=1:numel(files_BF)
    fprintf('\nImage number %d:\n',x);
    % count cells
    filename = fullfile(files_BF(x).folder, files_BF(x).name);
    [Num(x),Mean(x),CENT{x},ORIENT{x},label(:,:,x),BF_segmented(:,:,:,x), ~, ~] = Count_Cells(filename,0);
end

for x=1:numel(files_BF)
    filename = fullfile(yourfolder, segmentation_product_folder, sprintf('BF_segmented%d.jpg',x));
    imwrite(BF_segmented(:,:,:,x), filename, 'jpg');
end

save('BF_segmented.mat','BF_segmented');
save('Num.mat','Num');
save('label.mat', 'label');

%% get std and mean according to cell size
im_num_index = numel(files_BF);
countsCSH = zeros(numel(files_BF),size(histBinLocations,2));
centersCSH = zeros(numel(files_BF),size(histBinLocations,2));
for x=1:numel(files_BF)
    labelC = label(:,:,x);
    statsCSH = regionprops(labelC,'Area');
    statAreaCSH=[statsCSH.Area];
    statAreaCSH(statAreaCSH == 0) = NaN;
    [countsCSH(x,:),centersCSH(x,:)] = hist(statAreaCSH, histBinLocations);
end

fig = figure; hold on;  title('cell count'); xlabel('time [min]'); ylabel('cell count [a.u]');
% for x=1:size(histBinLocations,2)
%     plot((1:im_num_index).*image_time_period/bf_frequency, countsCSH(:,x),... 
%         linestyle(mod(x,3)+1), 'DisplayName', sprintf('Cell size: %d', floor(histBinLocations(x)./microscope_resolution)));
%     legend('-DynamicLegend');
%     legend('show');
%     drawnow;
% end
totalNumCells =sum(countsCSH(:,:), 2);
im_time_scale = (1:im_num_index).*image_time_period/bf_frequency;
plot((1:im_num_index).*image_time_period/bf_frequency, sum(countsCSH(:,:), 2), 'k', 'DisplayName', sprintf('Total'));
legend('-DynamicLegend');
legend('show');
filenameCNT = fullfile(yourfolder, graphs_folder, 'cell count.png');
saveas(fig,filenameCNT);
filenameCNT = fullfile(yourfolder, graphs_folder, 'cell count.fig');
saveas(fig,filenameCNT);

% [fitresult, gof] = cell_grouth_rate_fit(im_time_scale(1:20), totalNumCells(1:20));

%% build intensity mean and std depending on cell size
std_CSH = zeros(uint32(numel(files_GFP)/gfp_frequency),size(histBinLocations,2));
mean_CSH = zeros(uint32(numel(files_GFP)/gfp_frequency),size(histBinLocations,2));

im_num_index = uint32(numel(files_GFP)/gfp_frequency);
hist_values = {im_num_index};
hist_center = {im_num_index};
figure; hold on;
for x=1:im_num_index
    fprintf('\nGFP Image number %d:\n',x);
    filename = fullfile(files_GFP(gfp_frequency*x-(gfp_frequency-1)).folder, files_GFP(gfp_frequency*x-(gfp_frequency-1)).name);
    im_uint16 = imread(filename);
    im = im2double(im_uint16);
    labeli = label(:,:,gfp_frequency*x-(gfp_frequency-1));
    statsCSH = regionprops(labeli,'Area');
    statAreaCSH=[statsCSH.Area];
    statAreaCSH(statAreaCSH == 0) = nan;
    int_all = nan(1,numel(statAreaCSH));
    for xx=1:size(histBinLocations,2)
        fprintf('\nHist cell size %d:\n',histBinLocations(xx));
        if xx == 1
            cOB = ( (statAreaCSH <= (histBinLocations(xx) + bin_step/2)) );
        elseif xx == size(histBinLocations,2)
            cOB = ( (statAreaCSH > (histBinLocations(xx) + bin_step/2)) );
        else
            cOB = ((statAreaCSH > (histBinLocations(xx-1) + bin_step/2)) .* (statAreaCSH < (histBinLocations(xx) + bin_step/2)));
        end
        cellAVRG = zeros(1,numel(statAreaCSH));
        for cellIdx=1:numel(statAreaCSH)
            if cOB(cellIdx)==0
                cellAVRG(cellIdx) = nan;
                continue
            end
            mask = (labeli == cellIdx);
            masked_data = mask.*im;
            if isnan(statAreaCSH(cellIdx))
                cellAVRG(cellIdx) = nan;
            else
                cellAVRG(cellIdx) = sum(masked_data(:)); %./ statAreaCSH(cellIdx);
                int_all(cellIdx) = cellAVRG(cellIdx);
            end
        end
        std_CSH(x,xx) = nanstd(cellAVRG);
        mean_CSH(x,xx) = nanmean(cellAVRG);
    end
    h = histogram(int_all, 25);
    hist_values{x} = h.Values;
    shift_binEdge = [h.BinEdges]+h.BinWidth/2;
    hist_center{x} = shift_binEdge(1:end-1);
end
filenameMI_CSH = fullfile(yourfolder, graphs_folder, 'mean_CSH.mat');
save(filenameMI_CSH,'mean_CSH');
filenameSTD_CSH = fullfile(yourfolder, graphs_folder, 'std_CSH.mat');
save(filenameSTD_CSH,'std_CSH');
filenameHIST_CSH = fullfile(yourfolder, graphs_folder, 'hist_CSH.mat');
save(filenameHIST_CSH,'hist_values', 'hist_center');


%% print CSH stats intensity
im_num_index = uint32(numel(files_GFP)/gfp_frequency);

fig = figure; hold on; title('mean cell intensity'); xlabel('time [min]'); ylabel('intensity [a.u]');
for xx=1:size(histBinLocations,2)
    plot((1:im_num_index).*image_time_period*gfp_frequency, mean_CSH(:,xx),... 
        linestyle(mod(xx,3)+1), 'DisplayName', sprintf('Cell size: %d', floor(histBinLocations(xx)./microscope_resolution)));
    legend('-DynamicLegend');
    legend('show');
    drawnow;
end
filenameMI = fullfile(yourfolder, graphs_folder, 'mean cell intensity.png');
saveas(fig,filenameMI);
filenameMI = fullfile(yourfolder, graphs_folder, 'mean cell intensity.fig');
saveas(fig,filenameMI);

fig = figure; hold on; title('STD cell intensity'); xlabel('time [min]'); ylabel('intensity [a.u]');
for xx=1:size(histBinLocations,2)
    plot((1:im_num_index).*image_time_period*gfp_frequency, std_CSH(:,xx),... 
        linestyle(mod(xx,3)+1), 'DisplayName', sprintf('Cell size: %d', floor(histBinLocations(xx)./microscope_resolution)));
    legend('-DynamicLegend');
    legend('show');
    drawnow;
end
filenameSTD = fullfile(yourfolder, graphs_folder, 'std cell intensity.png');
saveas(fig,filenameSTD);
filenameSTD = fullfile(yourfolder, graphs_folder, 'std cell intensity.fig');
saveas(fig,filenameSTD);

fig = figure; hold on; title('SNR'); xlabel('time [min]'); ylabel('SNR');
for xx=1:size(histBinLocations,2)
    plot((1:im_num_index).*image_time_period*gfp_frequency, mean_CSH(:,xx)./std_CSH(:,xx),... 
        linestyle(mod(xx,3)+1), 'DisplayName', sprintf('Cell size: %d', floor(histBinLocations(xx)./microscope_resolution)));
    legend('-DynamicLegend');
    legend('show');
    drawnow;
end
filenameSNR = fullfile(yourfolder, graphs_folder, 'snr cell intensity.png');
saveas(fig,filenameSNR);
filenameSNR = fullfile(yourfolder, graphs_folder, 'snr cell intensity.fig');
saveas(fig,filenameSNR);

%% print intencity histagram for each image
im_num_index = uint32(numel(files_GFP)/gfp_frequency);
gauss_means = zeros(im_num_index,1);
gauss_vars = zeros(im_num_index,1);
figure; hold on; title('Gamma distribution of  microscopy'); xlabel('intensity [a.u]'); ylabel('cell amount [num]');
for x=10:10
    centersINT = double(hist_center{x});
    countINT = double(hist_values{x});
    scatter(centersINT, countINT);
    xLin = min(centersINT):0.1:max(centersINT);
    % initial paramter guesses (based on plot)
    initGuess(1) = 0.9*max(countINT);  % amplitude
    initGuess(2) = 400;       % mean centred on 0 degrees
    initGuess(3) = 100;      % SD in degrees
    % equation for Gamma distribution
    norm_func = @(p,x) p(1) .* exp(-((x - p(2))/p(3)).^2);
    % use nlinfit to fit Gamma using Least Squares
    [bestfit,resid]=fminsearch(@(p) norm(norm_func(p,centersINT)-countINT), initGuess);
    plot(xLin, norm_func(bestfit,xLin));
    if bestfit(2) < 1000
        %filter non realistic values
        gauss_means(x) = bestfit(2);
        gauss_vars(x) = bestfit(3);
    end
end
legend('-DynamicLegend');
legend('show');

figure; hold on;
plot(1:im_num_index, gauss_means, 'DisplayName', 'gauss mean');
plot(1:im_num_index, gauss_vars, 'DisplayName', 'gauss variance');
legend('-DynamicLegend');
legend('show');