%------------------------------------------------------------------------------------------------
% Script will scan all images from chosen folder and create data file (.xls) with mean and integrated
% intensities of nucleus, kinetochores and centrosomes
% Before running - add correct parameters - lines 29
%                - check parameters - The title of DAPI image should contain 'DAPI.w1DAPI'
%                                   - Channels for kinetochore and centrosome markers - 'Kin.w2FITC', 'Cen.w4Cy5'
%                                   - Channels for Ab on kinetochore and centrosome  - 'Plk1.w3Cy3' and 'Plk1.w5Cy3'
%                IF the file names are different then these, please correct
%                line 65 and lines 82-85.
%               
% When start running you will be asked to choose folder with your .tif
% files.  
%
% 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 lines 145, 146 
%                                 - activate lines  149, 150
%
% Data file (.xls) will be saved in MATLAB folder (root folder)
%------------------------------------------------------------------------------------------------

tic; % Start timer.
clc; % Clear command window.
clear all; % Get rid of variables from prior run of this m-file.
disp('Running IntensityKin_Cen_Zst.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.)


%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

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, 'DAPI.w1DAPI', 'Kin.w2FITC');
Plk1KinFile =regexprep (Fname, 'DAPI.w1DAPI', 'Plk1.w3Cy3');
CenFile =regexprep (Fname, 'DAPI.w1DAPI', 'Cen.w4Cy5');
Plk1CenFile =regexprep (Fname, 'DAPI.w1DAPI', 'Plk1.w5Cy3');

fullFileNameFITC = fullfile(listOfFolderNames{kk}, fitcFile);
fullFileNamePlk1Kin = fullfile(listOfFolderNames{kk}, Plk1KinFile);
fullFileNameCen = fullfile(listOfFolderNames{kk}, CenFile);
fullFileNamePlk1Cen = fullfile(listOfFolderNames{kk}, Plk1CenFile);
disp (fullFileNameFITC);

%!!!
SplitOfFnamex=(regexp (Fname, '--', 'split'));
%c=size(SplitOfFnamex');
cX=(SplitOfFnamex(1));
cY=(SplitOfFnamex(2));


 
cX=regexprep (cX, 'X', '');
cY=regexprep (cY, 'Y', '');



numImage1 = str2double(cX);
numImage2 = str2double(cY);
numImage = sprintf ('%d%d',numImage1,numImage2);
numImage =str2double(numImage);


%!!!

%END geting image number

% Read in the image.

originalImage = imread(fullFileName); 

originalImageFITC = imread(fullFileNameFITC);
originalImagePlk1Kin = imread(fullFileNamePlk1Kin);
originalImageCen = imread(fullFileNameCen);
originalImagePlk1Cen = imread(fullFileNamePlk1Cen);
subplot(2, 3, 1);
%Convert originalImage to 8 bit (just to be able to see better)
img8 = uint8(originalImage /128 );
%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.
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 to get a binary image (only 0's and 1's) of class "logical."

% Method #1: using im2bw()

doubleImage=double(originalImage)/65000;

normalizedThresholdValue = graythresh(doubleImage);
binaryImage = im2bw(originalImage, normalizedThresholdValue);   % One way to threshold to binary
   

%thresholdValue = 100;
%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');



% 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.
NucMeasurements = regionprops(labeledImage, originalImage, 'all'); 

 


numberOfNucs = size(NucMeasurements, 1);


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');
end


% Now select certain objects based using the ismember function.
% with an intensity between 100 and 4000 and an area less than 800 pixels.
allNucIntensities = [NucMeasurements.MeanIntensity];
allNucAreas = [NucMeasurements.Area];
% Get a list of coordinates for all objects
allNucCentroids = [NucMeasurements.Centroid];
% Get a list of minor and major AxisLength for all objects
allNucMaxAxis = [NucMeasurements.MajorAxisLength];
allNucMinAxis = [NucMeasurements.MinorAxisLength];
% Get a list of the Nucs that meet our criteria and we need to keep. 
allowableIntensityIndexes = (allNucIntensities > 100) & (allNucIntensities < 4000);
allowableAreaIndexes = (allNucAreas > 80) & (allNucAreas < 800); % Take the small objects.
allowableMaxAxisIndexes = (allNucMaxAxis > 5) & (allNucMaxAxis < 28);
allowableMinAxisIndexes = (allNucMinAxis > 3) & (allNucMinAxis < 14);
keeperIndexes = find(allowableIntensityIndexes & allowableAreaIndexes & allowableMaxAxisIndexes & allowableMinAxisIndexes);

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


if (NucMeasurements(k).Area > 400)
    
  %Find max
binaryImageFITC=originalImageFITC >= max(max(originalImageFITC));
labeledImageFITC = bwlabel(binaryImageFITC, 8);
kinCentroid=regionprops(labeledImageFITC, 'Centroid');

binaryImageCen=originalImageCen >= max(max(originalImageCen));
labeledImageCen = bwlabel(binaryImageCen, 8);
cenCentroid=regionprops(labeledImageCen, 'Centroid');

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

backROI=~labeledImageFITC;

labeledRegForKin = zeros(ImgH,ImgW);
labeledRegForCen = zeros(ImgH,ImgW);
labeledRegForKinBG  = zeros(ImgH,ImgW);
labeledRegForCenBG  = zeros(ImgH,ImgW);
labeledRegForKinBG1 = zeros(ImgH,ImgW);
labeledRegForCenBG1 = zeros(ImgH,ImgW);

%Protection for empty files

if (cenCentroid(1,1).Centroid(1,1)>10)&&(cenCentroid(1,1).Centroid(1,2)>10)&& (kinCentroid(1,1).Centroid(1,1)>10)&&(kinCentroid(1,1).Centroid(1,2)>10)

for ira =1 : ImgW
    for iras =1 : ImgH
        
              
        
        if (abs(cenCentroid(1,1).Centroid(1,1)-ira)<2)&&(abs(cenCentroid(1,1).Centroid(1,2)-iras)<2)
            labeledImageCen(iras,ira)=1;
            labeledRegForCen(iras-1,ira)=1;
                      
        end
         if (abs(kinCentroid(1,1).Centroid(1,1)-ira)<2)&&(abs(kinCentroid(1,1).Centroid(1,2)-iras)<2)
            labeledImageFITC(iras,ira)=1;
            labeledRegForKin(iras,ira)=1;
                        
         end
         
         if (abs(cenCentroid(1,1).Centroid(1,1)-ira)<3)&&(abs(cenCentroid(1,1).Centroid(1,2)-iras-1)<3)
            labeledRegForCenBG(iras,ira)=1;
                                  
        end
         if (abs(kinCentroid(1,1).Centroid(1,1)-ira)<3)&&(abs(kinCentroid(1,1).Centroid(1,2)-iras)<3)
            labeledRegForKinBG(iras,ira)=1;
                                    
         end
         
          if (abs(cenCentroid(1,1).Centroid(1,1)-ira)<4)&&(abs(cenCentroid(1,1).Centroid(1,2)-iras-1)<4)
            labeledRegForCenBG1(iras,ira)=1;
          end
       
         if (abs(kinCentroid(1,1).Centroid(1,1)-ira)<4)&&(abs(kinCentroid(1,1).Centroid(1,2)-iras)<4)
            labeledRegForKinBG1(iras,ira)=1;
                                    
         end
         
         
        %if (20<ira)&&(20<iras)&&(ira<ImgW-20)&&(iras<ImgH-20)
         if (10<ira)&&(10<iras)&&(ira<ImgW-10)&&(iras<ImgH-10)
           
            backROI(iras,ira)=0;
          
        end
    end
    
end




labeledRegForKinBG1=labeledRegForKinBG1-labeledRegForKinBG;
labeledRegForCenBG1=labeledRegForCenBG1-labeledRegForCenBG;
labeledRegForKinBG=labeledRegForKinBG-labeledRegForKin;
labeledRegForCenBG=labeledRegForCenBG-labeledRegForCen;



kinMeasurement=regionprops(labeledRegForKin, originalImagePlk1Kin, 'All');
cenMeasurement=regionprops(labeledRegForCen, originalImagePlk1Cen, 'All');
kinMeasurementBG=regionprops(labeledRegForKinBG, originalImagePlk1Kin, 'All');
cenMeasurementBG=regionprops(labeledRegForCenBG, originalImagePlk1Cen, 'All');
kinMeasurementBGG=regionprops(labeledRegForKinBG1, originalImagePlk1Kin, 'All');
cenMeasurementBGG=regionprops(labeledRegForCenBG1, originalImagePlk1Cen, 'All');
IntegrKin=(kinMeasurement.MeanIntensity-kinMeasurementBG.MeanIntensity)*kinMeasurement.Area;
IntegrCen=(cenMeasurement.MeanIntensity-cenMeasurementBG.MeanIntensity)*cenMeasurement.Area;

% End (Correction for MEAN )


nucPlk1= regionprops(labeledImage, originalImagePlk1Kin, 'all');
backPlk1= regionprops(backROI, originalImagePlk1Kin, 'all');
backDAPI = regionprops(backROI, originalImage, 'all'); 
    
positions (1,1)=numImage;
positions (1,2)=NucMeasurements(k).Area; %DAPI area
positions (1,3)=NucMeasurements(k).MeanIntensity; %DAPI mean
positions (1,4)=nucPlk1.MeanIntensity; %Plk1 mean
positions (1,5)=kinMeasurement.MaxIntensity;
positions (1,6)=kinMeasurement.MeanIntensity;
positions (1,7)=kinMeasurement.MinIntensity;
positions (1,8)=cenMeasurement.MaxIntensity;
positions (1,9)=cenMeasurement.MeanIntensity;
positions (1,10)=cenMeasurement.MinIntensity;
positions (1,11)=backPlk1.MeanIntensity;
positions (1,12)=backDAPI.MeanIntensity;
positions (1,13)=kinMeasurementBG.MeanIntensity;
positions (1,14)=cenMeasurementBG.MeanIntensity;
positions (1,15)=kinMeasurementBGG.MeanIntensity;
positions (1,16)=cenMeasurementBGG.MeanIntensity;
positions (1,17)=IntegrKin;
positions (1,18)=IntegrCen;




% Plot the borders of all the coins on the original grayscale image using the coordinates returned by bwboundaries.
subplot(2, 3, 4); imagesc(originalImageFITC);
title('Outlines, from bwboundaries()'); 
%axis square;
axis image;
hold on;
boundaries = bwboundaries(labeledImageFITC);	
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 all the coins on the original grayscale image using the coordinates returned by bwboundaries.
subplot(2, 3, 5); imagesc(originalImageCen);
title('Outlines, from bwboundaries()'); 
%axis square;
axis image;
hold on;
boundaries = bwboundaries(labeledImageCen);	
numberOfBoundaries = size(boundaries);
for k = 1 : numberOfBoundaries
	thisBoundary = boundaries{k};
	plot(thisBoundary(:,2), thisBoundary(:,1), 'g', 'LineWidth', 1);
end
hold off;


end %(protection for empty files)


end


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

% We have a labeled image of Nucs that meet our specified criteria.
subplot(2, 3, 6); 
imshow(originalImageCen);
%axis square;
axis image;
title('Objects after "neighbor" filter');




end
end
end

if exist('FinalPositions', 'var')
    
    
    head{1,1}='ImageNum';
    head{1,2}='Area_N';
    head{1,3}='DAPI_Mean';
    head{1,4}='Plk1_Mean';
    head{1,5}='Kin_Plk1_Max';
    head{1,6}='Kin_Plk1_Mean';
    head{1,7}='Kin_Plk1_Min';
    head{1,8}='Cen_Plk1_Max';
    head{1,9}='Cen_Plk1_Mean';
    head{1,10}='Cen_Plk1_Min';
    head{1,11}='Back_Plk1_Mean';
    head{1,12}='Back_DAPI_Mean';
    head{1,13}='Back_FOR_Kin_Mean';
    head{1,14}='Back_FOR_Cen_Mean';
    head{1,15}='Back1_FOR_Kin_Mean';
    head{1,16}='Back1_FOR_Cen_Mean';
    head{1,17}='Integr Kin';
    head{1,18}='Integr Cen';
    
    body='';
    
    
    for iii=1 : size (FinalPositions)
      
       body1=sprintf('%d, %s, %s, %s, %s, %s, %s, %s, %s, %s\r\n', FinalPositions(iii,1), num2str(FinalPositions(iii,2)), num2str(FinalPositions(iii,3)), num2str(FinalPositions(iii,4)), num2str(FinalPositions(iii,5)), num2str(FinalPositions(iii,6)), num2str(FinalPositions(iii,7)), num2str(FinalPositions(iii,8)), num2str(FinalPositions(iii,9)), num2str(FinalPositions(iii,10)));  
       body=sprintf('%s%s', body, body1); 
    end
    
 

 FileNameOut=sprintf('All_Data%s-%d.xls', exp, fileNumber); 



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



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

