%SwarMer Matlab graphing code
%Saumya Saurabh, Ph.D.
%saumya.saurabh@gmail.com

clc; % Clear command window.
close all;
clearvars; % Get rid of variables from prior run of this m-file.
fprintf('Running dropletAnalysis.m...\n'); % Message sent to command window.
workspace; % Make sure the workspace panel with all the variables is showing.
%imtool close all;  % Close all imtool figures.
format long g;
format compact;
captionFontSize = 14;
%uiwait(msgbox('Pahse image Segmetation; Three main numbers: thresholdValue, AreaHigh and AreaLow','Step1','modal'));
% Read in the phase image
%PathName = uigetdir(matlabroot,'Select folder with images');
analyzed_files = dir('*.tif');
%analyzed_files(2)=[]
%background_file = imread('DarkImageBackground-1.tif');


numSamples = length(analyzed_files);
figure(221)
for filenum = 1:length(analyzed_files)
    
FileName = char(analyzed_files(filenum,1).name); 
if ~isequal(FileName,0)
      originalImage1 = imread(FileName);
      %originalImage = imread('Project002_NoSat_451_noDP001_ch01.tif');
end

originalImage=originalImage1-5; %subtracting background constant; best to have an image
originalImage = rgb2gray(originalImage);
imagesize = size(originalImage);

% Display the grayscale image.

thresholdValue = 26;
%level = graythresh(originalImage);
%BW = imbinarize(originalImage,level);
BW = originalImage > thresholdValue;


subplot(3, numSamples, filenum);
imagesc(originalImage), colormap(hot);
% Maximize the figure window.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]);
% Force it to display RIGHT NOW (otherwise it might not display until it's all done, unless you've stopped at a breakpoint.)
drawnow;
caption = strcat('Replicate#',num2str(filenum));
title(caption, 'FontSize', captionFontSize);
axis image; % Make sure image is not artificially stretched because of screen's aspect ratio.

subplot(3, numSamples, filenum+numSamples);
imshow(BW)

BWnobord = uint16(imclearborder(BW, 4));


labeledImage = bwlabel(BWnobord, 4); 
blobMeasurements = regionprops(labeledImage, originalImage, 'all');
blobMeasurements_size = size(blobMeasurements,1);

allBlobAreas = [blobMeasurements.Area];
allowableAreaIndexes = allBlobAreas > 5000; % Take the small objects out.
keeperIndexes = find(allowableAreaIndexes);
keepercellsImage = ismember(labeledImage, keeperIndexes);
C1 = uint8(keepercellsImage).*(originalImage);
C = imfill(C1)
C_imadjust = imadjust(C);
subplot(3, numSamples, filenum+numSamples+numSamples);
imshow(C_imadjust);
labelednewImage = bwlabel(keepercellsImage, 8);  
Col_Area = regionprops(labelednewImage, C,'Area');
SwMotDiff(filenum,:) = [Col_Area.Area]
end

SwMotDiff(:,3) = SwMotDiff(:,1)-SwMotDiff(:,2);
%SwMotDiff(:,4) = SwMotDiff(:,2)./SwMotDiff(:,1);

% SwMotDiff(filenum+1,:) = [mean(SwMotDiff(:,1)) mean(SwMotDiff(:,2)) mean(SwMotDiff(:,3))]; %mean
% SwMotDiff(filenum+2,:) = [std(SwMotDiff(:,1)) std(SwMotDiff(:,2)) std(SwMotDiff(:,3))]; %Stdev
% SwMotDiff(filenum+3,:) = SwMotDiff(filenum+1,:)./SwMotDiff(filenum+1,1) %Normalized mean
% SwMotDiff(filenum+4,:) = SwMotDiff(filenum+2,:)./SwMotDiff(filenum+1,1) %Normalized Stdev
% %plot and p-value
% x = 1:2;
% data = [SwMotDiff(filenum+3,1) SwMotDiff(filenum+3,2)]
% errhigh = [(SwMotDiff(filenum+4,1)/2) (SwMotDiff(filenum+4,2)/2)];
% errlow = [(SwMotDiff(filenum+4,1)/2) (SwMotDiff(filenum+4,2)/2)];
% [h,p,ci,stats] = ttest2(SwMotDiff(1:filenum,1),SwMotDiff(1:filenum,2))
% figure(291)
% bar(x,data,'LineWidth',3,'FaceColor','w')
% hold on
% er = errorbar(x,data,errlow,errhigh,'LineWidth',3);  
% er.Color = [0 0 0];                            
% er.LineStyle = 'none';  
% datanames={'WT'; '\DeltaIDR'};
% set(gca,'xticklabel',datanames)
% set(gca,'LineWidth',3)
% ylabel('Normalized Swarming Motility')
% txt = strcat('P-value = ', num2str(p));
% text(2.5,1,txt)
% ylim([0 1.1])
% hold off

%save file
fileout = strcat(FileName(1:end-7),'_analysis.mat');
save(fileout);
