%------------------------------------------------------------------------------------------------
% Script will scan all images from chosen folder and create data file (.xls) with mean and integrated
% intensities of nucleus and cytoplasm for every channel
% Before running - add correct parameters - lines 29-31
%                - check parameters - The title of DAPI image should contain 'w2DAPI'
%                                   
%                                   - Two other channels - 'w3Cy3' and 'w4FITC'
%                IF the file names are different then these, please correct
%                line 74 and lines 91-93.
%               
% When start running you will be asked to choose folder with your .tif
% files and destination folder where data file will be saved.  
%
% Software will automatically calculate the threshold of image using a logical operation to get a
% binary image, but there is a possibility to manually put threshold value
% if needed. For this you should: - coment line 138 
%                                 - activate line  140
%
% Data file (.xls) will be saved in destination folder you have chosen
%------------------------------------------------------------------------------------------------

tic; % Start timer.
clc; % Clear command window.
clear all; % Get rid of variables from prior run of this m-file.
disp('Running Intensity_Pattern_Not_labeled.m...'); % Message sent to command window.
workspace; % Show panel with all the variables.

%Add all parameters
exp='Esperiment-xx';  % Name/number of experiment. (Any name. Will be used only in the title of data file.)
Cy3='CyB';           % Name for channel I
FITC='CyA';          % Name for channel II


%END adding parameters


% Change the current folder to the folder of this m-file.

if(~isdeployed)
	cd(fileparts(which(mfilename)));
end


path =cd;
topLevelFolder = uigetdir(path);
if topLevelFolder == 0
return;
end

path =cd;
DestFolder = uigetdir(path);
if DestFolder == 0
return;
end


allSubFolders = genpath(topLevelFolder);
%disp('Here are the folders separated:');

% Scan through them separating them.
remain = allSubFolders;
listOfFolderNames = {};
while true % Code adapted from the help file.
[singleSubFolder, remain] = strtok(remain, ';');
if isempty(singleSubFolder), break; end
%disp(sprintf('%s', singleSubFolder))
listOfFolderNames = [listOfFolderNames singleSubFolder];
end
% Now we have a list of all folders and subfolders.

for kk = 1 : length(listOfFolderNames)
% Get the file pattern for dir().
%filePattern = [listOfFolderNames{kk}, '\*.tif'];
filePattern = [listOfFolderNames{kk}, '\*DAPI*.TIF'];
% Find out names of TIF files in this folder.
tifFiles = dir(filePattern); % Cell array
if ~isempty(tifFiles)
% If there are any TIF files, open them all
for fileNumber = 1 :length(tifFiles)
fullFileName = fullfile(listOfFolderNames{kk}, tifFiles(fileNumber).name);

% Get rid of variables from prior run of other image.
clear ('numberOfNucs', 'NucMeasurements','allNucIntensities', 'keeperIndexes', 'keeperIndexes1', 'positions', 'mayBE','allNucCentroids','allNucAreas')


% Geting image number

Fname = tifFiles(fileNumber).name;

disp (Fname);
fitcFile =regexprep (Fname, 'w1DAPI', 'w3FITC');
Cy3File =regexprep (Fname, 'w1DAPI', 'w2Cy3');


fullFileNameFITC = fullfile(listOfFolderNames{kk}, fitcFile);
fullFileNameCy3 = fullfile(listOfFolderNames{kk}, Cy3File);

disp (fullFileNameFITC);

%!!!
SplitOfFnamex=(regexp (Fname, '_', 'split'));
c=size(SplitOfFnamex');
PartOfFname=(SplitOfFnamex(c(1)));
PartOfFname=regexprep (PartOfFname, '.tif', '');
PartOfFname=regexprep (PartOfFname, '.TIF', '');
PartOfFname=regexprep (PartOfFname, 's', '');
numImage = str2double(PartOfFname);


%!!!

%END geting image number

% Read in the image.

originalImage = imread(fullFileName); 

originalImageFITC = imread(fullFileNameFITC);
originalImageCy3 = imread(fullFileNameCy3);

subplot(2, 3, 1);
%Convert originalImage to 8 bit (just to be able to see better)
img8 = uint8(originalImage /8 );
%img8 = uint8(double(originalImage)/(2^16-1))*255;
imshow(img8);
% Maximize the figure window.
set(gcf, 'Position', get(0, 'ScreenSize'));
% 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 = sprintf('Original image.');
title(caption); 
%axis square; % Make sure image is not artificially stretched because of screen's aspect ratio.
axis image;


% Threshold the image using a logical operation to get a binary image.

thresholdValue=round (mean(min(originalImage))+mean(max(originalImage))*0.6);

%thresholdValue = 100; % Alternatively, threshold value added manually 

disp (median(min(originalImage)));
disp (thresholdValue);

binaryImage = originalImage > thresholdValue; % % Select bright objects (nucl).

% Do a "hole fill" to get rid of any background pixels inside the Nucs.
binaryImage = imfill(binaryImage, 'holes');
%binaryImagePattern = imfill(binaryImagePattern, 'holes');


% Display the binary image.
subplot(2, 3, 2); imagesc(binaryImage); colormap(gray(256)); 
title('Binary Image, obtained by thresholding'); 
%axis square;
axis image;

labeledImage = bwlabel(binaryImage, 8);     % Label each Nuc so we can make measurements of it

coloredLabels = label2rgb (labeledImage, 'hsv', 'k', 'shuffle'); % pseudo random color labels



subplot(2, 3, 3); 
imagesc(coloredLabels);
%axis square;
axis image;
caption = sprintf('Pseudo colored labels, from label2rgb().\nObjects are numbered from top to bottom, then from left to right.');
title(caption);

% Get all the Nuc properties.  Can only pass in originalImage in version R2008a and later.
labeledImage1=labeledImage; %making logical matrix to use less memory
NucMeasurements = regionprops(labeledImage1, originalImage, 'all'); 

 


numberOfNucs = size(NucMeasurements, 1);


checkFor2big=0;
for cck = 1 : numberOfNucs
    if (NucMeasurements(cck).Area > 250)
        checkFor2big=checkFor2big+1;
    end
end


% bwboundaries() returns a cell array, where each cell contains the row/column coordinates for an object in the image.


fontSize = 14;	% Used to control size of "Nuc number" labels put atop the image.
labelShiftX = -7;	% Used to align the labels in the centers of the coins.
NucECD = zeros(1, numberOfNucs);
% Print header line in the command window.
fprintf(1,'Nuc #      Mean Intensity  Area   Perimeter    Centroid       Diameter\n');
% Loop over all Nucs printing their measurements to the command window.
for k = 1 : numberOfNucs           % Loop through all Nucs.
	% Find the mean of each Nuc.  (R2008a has a better way where you can pass the original image
	% directly into regionprops.  The way below works for all versions including earlier versions.)
    thisNucsPixels = NucMeasurements(k).PixelIdxList;  % Get list of pixels in current Nuc.
    meanGL = mean(originalImage(thisNucsPixels)); % Find mean intensity (in original image!)
	meanGL2008a = NucMeasurements(k).MeanIntensity; % Mean again, but only for version >= R2008a
	
	NucArea = NucMeasurements(k).Area;		% Get area.
	NucPerimeter = NucMeasurements(k).Perimeter;		% Get perimeter.
	NucCentroid = NucMeasurements(k).Centroid;		% Get centroid.
	NucECD(k) = sqrt(4 * NucArea / pi);					% Compute ECD - Equivalent Circular Diameter.
    fprintf(1,'#%2d %17.1f %11.1f %8.1f %8.1f %8.1f % 8.1f\n', k, meanGL, NucArea, NucPerimeter, NucCentroid, NucECD(k));
	% Put the "Nuc number" labels on the "boundaries" grayscale image.
	%text(NucCentroid(1) + labelShiftX, NucCentroid(2), num2str(k), 'FontSize', fontSize, 'FontWeight', 'Bold');





% Extract only those objects that meet our criteria, and 
% eliminate those objects that don't meet our criteria. 

if (NucMeasurements(k).Area > 250)&&(checkFor2big < 2)
    
  %Find max
%patternMes=regionprops(labeledImagePattern, 'all');

%patternCentroid=regionprops(labeledImagePattern, 'Centroid');

ImgHH = size(originalImage);
ImgWW = size(originalImage');
ImgH=(ImgHH(1,1));
ImgW=(ImgWW(1,1));



%labeledRegForKin = zeros(ImgH,ImgW);
labeledRegForPat = zeros(ImgH,ImgW);
%labeledRegForKinBG  = zeros(ImgH,ImgW);
labeledRegForPatBG  = zeros(ImgH,ImgW);
%labeledRegForKinBG1 = zeros(ImgH,ImgW);
labeledRegForPatBG1 = zeros(ImgH,ImgW);

%Protection for empty files and checking if nucleus

if (NucCentroid(1,1)>30)&& (NucCentroid(1,2)>30)&& (abs(ImgW-NucCentroid(1,1))>30)&& (abs(ImgH-NucCentroid(1,2))>30)
    
for ira =1 : ImgW
    for iras =1 : ImgH
        
              
        
        if (abs(NucCentroid(1,1)-ira)<84)&&(abs(NucCentroid(1,2)-iras)<25)
            labeledImagePatDr(iras,ira)=1;
            labeledRegForPat(iras,ira)=1;
                      
        end
        
         
         if (abs(NucCentroid(1,1)-ira)<87)&&(abs(NucCentroid(1,2)-iras)<28)
            labeledRegForPatBG(iras,ira)=1;
                                  
        end
        
         
          if (abs(NucCentroid(1,1)-ira)<90)&&(abs(NucCentroid(1,2)-iras)<31)
            labeledRegForPatBG1(iras,ira)=1;
          end
       
         
         
         
        
    end
    
end





labeledRegForPatBG1=labeledRegForPatBG1-labeledRegForPatBG;



labeledRegForCell=labeledRegForPat-labeledImage;
labeledRegForCell(labeledRegForCell<0)=0; % negative values to 0

CellCy3Measurement=regionprops(labeledRegForCell, originalImageCy3, 'All');
CellFITCMeasurement=regionprops(labeledRegForCell, originalImageFITC, 'All');
NucCy3Measurement=regionprops(labeledImage, originalImageCy3, 'All');
NucFITCMeasurement=regionprops(labeledImage, originalImageFITC, 'All');

Cy3MeasurementBGG=regionprops(labeledRegForPatBG1, originalImageCy3, 'All');
FITCMeasurementBGG=regionprops(labeledRegForPatBG1, originalImageFITC, 'All');
DAPIMeasurementBGG=regionprops(labeledRegForPatBG1, originalImage, 'All');
IntegrCy3=(CellCy3Measurement.MeanIntensity-Cy3MeasurementBGG.MeanIntensity)*CellCy3Measurement.Area;
IntegrFITC=(CellFITCMeasurement.MeanIntensity-FITCMeasurementBGG.MeanIntensity)*CellFITCMeasurement.Area;
IntegrNucCy3=(NucCy3Measurement(k).MeanIntensity-Cy3MeasurementBGG.MeanIntensity)*NucMeasurements(k).Area;
IntegrNucFITC=(NucFITCMeasurement(k).MeanIntensity-FITCMeasurementBGG.MeanIntensity)*NucMeasurements(k).Area;
IntegrDAPI=(NucMeasurements(k).MeanIntensity-DAPIMeasurementBGG.MeanIntensity)*NucMeasurements(k).Area;
% End (Correction for MEAN )


    
positions (1,1)=numImage;
positions (1,2)=NucMeasurements(k).Area; %DAPI area
positions (1,3)=CellCy3Measurement.Area; %cyt area
positions (1,4)=NucMeasurements(k).MeanIntensity; %DAPI mean
positions (1,5)=NucCy3Measurement(k).MeanIntensity; %Nuc Cy3 mean
positions (1,6)=NucFITCMeasurement(k).MeanIntensity; %Nuc FITC mean
positions (1,7)=CellCy3Measurement.MeanIntensity; %cyt Cy3 mean
positions (1,8)=CellFITCMeasurement.MeanIntensity; %cyt FITC mean
positions (1,9)=DAPIMeasurementBGG.MeanIntensity; %DAPI background
positions (1,10)=Cy3MeasurementBGG.MeanIntensity; %Cy3 background
positions (1,11)=FITCMeasurementBGG.MeanIntensity; %FITC background
positions (1,12)=IntegrDAPI; %Integr DAPI
positions (1,13)=IntegrNucCy3; %Integr Nuc Cy3
positions (1,14)=IntegrNucFITC; %Integr Nuc FITC
positions (1,15)=IntegrCy3; %Integr cyt Cy3
positions (1,16)=IntegrFITC; %Integr cyt FITC







% Plot the borders of the cell on the original grayscale image using the coordinates returned by bwboundaries.
subplot(2, 3, 4); imagesc(originalImageCy3);
title('Original image Cy3'); 
%axis square;
axis image;
hold on;
boundaries = bwboundaries(labeledRegForCell);
	
numberOfBoundaries = size(boundaries);
for k = 1 : numberOfBoundaries
	thisBoundary = boundaries{k};
	plot(thisBoundary(:,2), thisBoundary(:,1), 'g', 'LineWidth', 1);
end
hold off;

% Plot the borders of the cell on the original grayscale image using the coordinates returned by bwboundaries.
subplot(2, 3, 5); imagesc(originalImageFITC);
title('Original image FITC'); 
%axis square;
axis image;
hold on;
boundaries = bwboundaries(labeledRegForCell);	
numberOfBoundaries = size(boundaries);
for k = 1 : numberOfBoundaries
	thisBoundary = boundaries{k};
	plot(thisBoundary(:,2), thisBoundary(:,1), 'g', 'LineWidth', 1);
end
hold off;



if exist('positions', 'var')
    if exist('FinalPositions', 'var')
FinalPositions=[FinalPositions ; positions];
    else
        FinalPositions=positions;
    end
end


end %(protection for empty files)









end




end
end
end
end

if exist('FinalPositions', 'var')
    
    
    head{1,1}='ImageNum';
    head{1,2}='Area_N';
    head{1,3}='Area_Cyt';
    head{1,4}='DAPI_Mean';
    head{1,5}=sprintf('N_%s_Mean', Cy3);
    head{1,6}=sprintf('N_%s_Mean', FITC);
    head{1,7}=sprintf('Cyt_%s_Mean', Cy3);
    head{1,8}=sprintf('Cyt_%s_Mean', FITC);
    head{1,9}='Back_DAPI';
    head{1,10}=sprintf('Back_%s_Mean', Cy3);
    head{1,11}=sprintf('Back_%s_Mean', FITC);
    head{1,12}='Int_DAPI';
    head{1,13}=sprintf('Int_%s_N', Cy3);
    head{1,14}=sprintf('Int_%s_N', FITC);
    head{1,15}=sprintf('Int_%s_Cyt', Cy3);
    head{1,16}=sprintf('Int_%s_Cyt', FITC);
    
   
    
   
    
  

 FileNameOut=sprintf('%s/All_Data%s.xls', DestFolder, exp); 



out = [head; num2cell(FinalPositions)];
s = xlswrite(FileNameOut, out);



disp('GAME OVER! )) File Name:');
disp (FileNameOut);
else
  disp ('Sorry no any position is found');  
end

