%------------------------------------------------------------------------------------------------
% Script will scan all images from choosen folder and make a position list of single cells 
% Script is addapted for 2 different microscope/software system 
%                                       1. ImageXpress / MetaXpress (Metamorph) 
%                                       2. Leica / LAS Matrix
% Any other microscope /software can be easily added 
% Before running - choose correct microscope/ software - line 38
%                - add correct parameter - lines 35
%
%     !!! Only if microscope N1 (ImageXpress is choosen)
%
%                - add correct parameters - lines 48-52
%                - check parameters - lines 55-58
%     !!! 
%
% When start running you will be asked to choose folder with 10x images (Nucleus staining) 
%
% 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: - choose option for manual thershold - line 40 
%                                 - add threshold value - line 41
%
% Position list 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 PosiyionList.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 position list.)


microscope = 1;  % Choose correct microscope (1 - ImageXpress / MetaXpress (Metamorph); 2 - Leica / LAS Matrix)

ThreshManual = 0; % Threshold definition ( 0 - automatically, 1 - manually)
ThreshValue = 100; % Threshold value (for manual only)

if (microscope ==1)


% Only if microscope N1 was choosen

obj=40;                     % Magnification that will be used for acquisition of single image. Put 60 for 60x and 40 for 40x
column = 17;                 % Number of columns when doing 10x acquisition (check microscope software) 
row = 10;                   % Number of rows when doing 10x acquisition (check microscope software)
ValX0 = 34344.18;           % The X position of left-up image when doing 10x acquisition (check microscope software)
ValY0 = 68384.14;           % The Y position of left-up image when doing 10x acquisition (check microscope software)


SizeX = 696;                %image width in px
SizeY = 520;                %image height in px
imgwUm = 897.84;            %image width in um
imghUm = 670.8;             %image height in um

%END adding parameters

SizeOnePixX = imgwUm/SizeX;
SizeOnePixY = imghUm/SizeY;

corX=imgwUm/2;                   
corY=imghUm/2;                   


%Coordinates map for images

yvalue = ValY0;
r=1;

for ck=1 : row
    xvalue = ValX0;
    for ckc=1 : column
                map(r,1)=xvalue;
        map(r,2)=yvalue;
        
       xvalue=xvalue+imgwUm;
       r=r+1;
    end
    yvalue=yvalue+imghUm;
end


%END of making coordinates map for images

% 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'];
% 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 ('numberOfNucl', 'nuclMeasurements','allNuclIntensities', 'keeperIndexes', 'keeperIndexes1', 'positions', 'mayBE','allNuclCentroids','allNuclAreas')


% Getting image number

Fname = tifFiles(fileNumber).name;
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);


disp (Fname);

% END getting image number

% Calculation of left-up position of the image 

LeftUpPositionX = map (numImage, 1)- corX;
LeftUpPositionY = map (numImage, 2)- corY;


% Read in the image.

originalImage = imread(fullFileName); 
subplot(2, 3, 1);
%Convert originalImage to 8 bit (just to be able to see better)
img8 = uint8(originalImage /2 );
%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;
axis image;


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

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

if (ThreshManual ~=0) % If manual thershhold value was choosen
thresholdValue = ThreshValue; % Alternatively, threshold value added manually (see line 41)
end

disp (thresholdValue);
disp (mean(min(originalImage)));
disp (mean(max(originalImage)));



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


% Do a "hole fill" to get rid of any background pixels inside the nucl.
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 nucl 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.\nObjects are numbered from top to bottom, then from left to right.');
title(caption);

% Get all the nucl properties.  Can only pass in originalImage in version R2008a and later.
nuclMeasurements = regionprops(labeledImage, originalImage, 'all');   
numberOfNucl = size(nuclMeasurements, 1);

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

fontSize = 14;	% Used to control size of "nucl number" labels put atop the image.
labelShiftX = -7;	% Used to align the labels in the centers of the nucl.
nucleusECD = zeros(1, numberOfNucl);
% Print header line in the command window.
fprintf(1,'Nucl #      Mean Intensity  Area   Perimeter    Centroid       Diameter\n');
% Loop over all nucl printing their measurements to the command window.
for k = 1 : numberOfNucl           % Loop through all nucl.
	% Find the mean of each nucl.  
    thisNuclPixels = nuclMeasurements(k).PixelIdxList;  % Get list of pixels in current nucl.
    meanGL = mean(originalImage(thisNuclPixels)); % Find mean intensity (in original image!)
	meanGL2008a = nuclMeasurements(k).MeanIntensity; % Mean again, but only for version >= R2008a
	
	nuclArea = nuclMeasurements(k).Area;		% Get area.
	nuclPerimeter = nuclMeasurements(k).Perimeter;		% Get perimeter.
	nuclCentroid = nuclMeasurements(k).Centroid;		% Get centroid.
	nucleusECD(k) = sqrt(4 * nuclArea / pi);					% Compute ECD - Equivalent Circular Diameter.
    fprintf(1,'#%2d %17.1f %11.1f %8.1f %8.1f %8.1f % 8.1f\n', k, meanGL, nuclArea, nuclPerimeter, nuclCentroid, nucleusECD(k));
	% Put the "nucl number" labels on the "boundaries" grayscale image.
	text(nuclCentroid(1) + labelShiftX, nuclCentroid(2), num2str(k), 'FontSize', fontSize, 'FontWeight', 'Bold');
end


% Select certain objects based using the ismember function.
% with an intensity between 80 and 4000 and an area less than 800 pixels.
allNuclIntensities = [nuclMeasurements.MeanIntensity];
allNuclAreas = [nuclMeasurements.Area];
% Get a list of coordinates for all objects
allNuclCentroids = [nuclMeasurements.Centroid];
% Get a list of minor and major AxisLength for all objects
allNuclMaxAxis = [nuclMeasurements.MajorAxisLength];
allNuclMinAxis = [nuclMeasurements.MinorAxisLength];
% Get a list of the nucl that meet our criteria and we need to keep. 
allowableIntensityIndexes = (allNuclIntensities > 80) & (allNuclIntensities < 4000);
allowableAreaIndexes = (allNuclAreas > 60) & (allNuclAreas < 800); % Take the small objects.
allowableMaxAxisIndexes = (allNuclMaxAxis > 5) & (allNuclMaxAxis < 28);
allowableMinAxisIndexes = (allNuclMinAxis > 3) & (allNuclMinAxis < 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. 
keeperNuclImage = ismember(labeledImage, keeperIndexes); 
% Re-label with only the keeper nucl kept.
labeledDimeImage = bwlabel(keeperNuclImage, 8);     % Label each object 

subplot(2, 3, 5); 
imshow(labeledDimeImage);
%axis square;
axis image;
title('"Keeper" objects (after "size" filter)');


keeperIndexes1=[];
keeperIndCol=keeperIndexes';
col=size(keeperIndCol);

for ch=1: col
    
  mayBE(ch,1:2) = nuclMeasurements(keeperIndexes(1,ch)).Centroid;
  mayBE(ch,3) = keeperIndexes(1,ch);
  
end


%Counting neighbors 


dist=50; %Distance to other cell

distB=20; %Distance to Image Border

%Correction for 40x center
corrX=20;
corrY=40;

%Correction for 60x center
if (obj == 60)
corrX=40;
corrY=20;
end

% Filter for cell too close to each other

i=1;
for nn=1 : col
raz=0;
fak=0;
   for mm=1: numberOfNucl
       znachenieX=abs(mayBE(nn,1)- allNuclCentroids(1,mm+raz));
       znachenieY=abs(mayBE(nn,2)- allNuclCentroids(1,mm+raz+1));
       
       
      
       if (znachenieX>0)||(znachenieY>0)
           
           if (znachenieX <dist)&&(znachenieY<dist)
           fak=1;
           
           break
          
           end
       end
        raz=raz+1;
     
   end
      
  % Filter for cell too close to image border
  
 if (fak==0)
     
     distFromImageBorderX=SizeX-mayBE(nn,1);
     distFromImageBorderY=SizeY-mayBE(nn,2);
     
     if (distFromImageBorderX>distB) && (distFromImageBorderX<SizeX-distB)&&(distFromImageBorderY>distB) && (distFromImageBorderY<SizeY-distB)
   
     keeperIndexes1(1,i)=mayBE(nn,3);
     positions(i,1:2)= mayBE(nn,1:2);
     
  % Keep positions of cells that passed filters
  
     interX= LeftUpPositionX + positions(i,1)*SizeOnePixX - corrX;
     interY= LeftUpPositionY + positions(i,2)*SizeOnePixY + corrY;
     positions(i,3) = interX;
     positions(i,4) = interY;
     
         i=i+1;
         
     end   
 end
end


% If position list already exist add new positions, 
% otherwise create new position list

if exist('positions', 'var')
    if exist('FinalPositions', 'var')
FinalPositions=[FinalPositions ; positions];
    else
        FinalPositions=positions;
    end
end
keeperNuclImage= ismember(labeledImage, keeperIndexes1);
% Re-label with only the keeper nucl.
labeledDimeImage= bwlabel(keeperNuclImage,8);     % Label each nucl
subplot(2, 3, 6); 
imshow(labeledDimeImage);
%axis square;
axis image;
title('Objects after "neighbor" filter');




end
end
end


% Create position list in special format (for microscope)

if exist('FinalPositions', 'var')
    
    head='"Stage Memory List", Version 5.0';
    head2='0, 0, 0, 0, 0, 0, 0, "um", "um"';
    footer='0';
    body='';
    PositionZ='9290';
    %PositionZ = num2str(PositionZ);
    for iii=1 : size (FinalPositions)
      
      PositionX=FinalPositions(iii,3);
      PositionX=round(PositionX*100)/100;
      PositionX=num2str(PositionX);
      
      PositionY=FinalPositions(iii,4);
      PositionY=round(PositionY*100)/100;
      PositionY=num2str(PositionY);
      
      %PositionY=num2str(FinalPositions(iii,4));
      %PositionX=num2str(FinalPositions(iii,3));
      body1=sprintf('"%s%d", %s, %s, %s, 0, 0, FALSE, -9999, TRUE, TRUE, 1\r\n', exp, iii, PositionX,PositionY,PositionZ);  
       body=sprintf('%s%s', body, body1); 
    end
    
 content=sprintf('%s\r\n%s\r\n%s\r\n%d\r\n%s', head , head2 , footer, iii, body );    

 % Save position list on the hard drive
 
 FileNameXmlOut=sprintf('Positions%s-%d.STG', exp, iii); 
fid=fopen(FileNameXmlOut,'w' ); 
fprintf(fid, content); 
fclose=(fid);

disp('GAME OVER! )) File Name:');
disp (FileNameXmlOut);
else
  disp ('Sorry no any single position was found');  
end

end


% Socond microscope


 if (microscope == 2)
        
   % 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'];
% 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 ('numberOfNucl', 'nuclMeasurements','allNuclIntensities', 'keeperIndexes', 'keeperIndexes1', 'positions', 'mayBE','allNuclCentroids','allNuclAreas','C')

% Get metadata (xml) file
addPart='\metadata';

pat1='--C00.ome.tif'; % Difference in file name.
xmlFilename=regexprep (tifFiles(fileNumber).name, pat1, '.ome.xml');

xmlFile = fullfile(listOfFolderNames{kk}, addPart, xmlFilename);
disp (xmlFile);
C=xml2struct(xmlFile);

% Get information from metadata file
PositionX = str2double(C.OME.Image.Pixels.Plane.StagePosition.Attributes.PositionX);

PositionY = str2double(C.OME.Image.Pixels.Plane.StagePosition.Attributes.PositionY);

SizeOnePixX = str2double(C.OME.Image.Pixels.Attributes.PhysicalSizeX);

SizeOnePixY = str2double(C.OME.Image.Pixels.Attributes.PhysicalSizeY);

SizeX = str2double(C.OME.Image.Pixels.Attributes.SizeX);
SizeY = str2double(C.OME.Image.Pixels.Attributes.SizeY);

LeftUpPositionX=PositionX*1000000 +SizeOnePixX * SizeX / 2; % Calculate left-up X coordinate of image
LeftUpPositionY=PositionY*1000000 -SizeOnePixY * SizeY / 2; % Calculate left-up Y coordinate of image


%END xmlFILE


% Read in the image.

originalImage = imread(fullFileName); 
subplot(2, 3, 1);
%Convert originalImage to 8 bit (just to be able to see better)
img8 = uint8(originalImage /64 );
%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 using a logical operation to get a binary image.

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

if (ThreshManual ~=0) % If manual thershhold value was choosen
thresholdValue = ThreshValue; % Alternatively, threshold value added manually (see line 41)
end 

disp (thresholdValue);
disp (mean(min(originalImage)));
disp (mean(max(originalImage)));



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


% Do a "hole fill" to get rid of any background pixels inside the nucl.
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;

labeledImage = bwlabel(binaryImage, 8);     % Label each nucl 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 nucl properties.  Can only pass in originalImage in version R2008a and later.
nuclMeasurements = regionprops(labeledImage, originalImage, 'all');   
numberOfNucl = size(nuclMeasurements, 1);

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

fontSize = 14;	% Used to control size of "nucl number" labels put atop the image.
labelShiftX = -7;	% Used to align the labels in the centers of the nucl.
nuclECD = zeros(1, numberOfNucl);
% Print header line in the command window.
fprintf(1,'Nucl #      Mean Intensity  Area   Perimeter    Centroid       Diameter\n');
% Loop over all blobs printing their measurements to the command window.
for k = 1 : numberOfNucl           % Loop through all nucl.
	% Find the mean of each nucl. 
    thisBlobsPixels = nuclMeasurements(k).PixelIdxList;  % Get list of pixels in current nucl.
    meanGL = mean(originalImage(thisBlobsPixels)); % Find mean intensity (in original image!)
	meanGL2008a = nuclMeasurements(k).MeanIntensity; % Mean again, but only for version >= R2008a
	
	nuclArea = nuclMeasurements(k).Area;		% Get area.
	nuclPerimeter = nuclMeasurements(k).Perimeter;		% Get perimeter.
	nuclCentroid = nuclMeasurements(k).Centroid;		% Get centroid.
	nuclECD(k) = sqrt(4 * nuclArea / pi);					% Compute ECD - Equivalent Circular Diameter.
    fprintf(1,'#%2d %17.1f %11.1f %8.1f %8.1f %8.1f % 8.1f\n', k, meanGL, nuclArea, nuclPerimeter, nuclCentroid, nuclECD(k));
	% Put the "nucl number" labels on the "boundaries" grayscale image.
	text(nuclCentroid(1) + labelShiftX, nuclCentroid(2), num2str(k), 'FontSize', fontSize, 'FontWeight', 'Bold');
end


% Select certain objects based using the ismember function.
% with an intensity between 100 and 65000 and an area less than 400 pixels.
allNuclIntensities = [nuclMeasurements.MeanIntensity];
allNuclAreas = [nuclMeasurements.Area];
% Get a list of coordinates for all objects
allNuclCentroids = [nuclMeasurements.Centroid];
% Get a list of minor and major AxisLength for all objects
allNuclMaxAxis = [nuclMeasurements.MajorAxisLength];
allNuclMinAxis = [nuclMeasurements.MinorAxisLength];
% Get a list of the blobs that meet our criteria and we need to keep. 
allowableIntensityIndexes = (allNuclIntensities > 100) & (allNuclIntensities < 65000);
allowableAreaIndexes = (allNuclAreas > 20) & (allNuclAreas < 400); % Take the small objects.
allowableMaxAxisIndexes = (allNuclMaxAxis > 6) & (allNuclMaxAxis < 16);
allowableMinAxisIndexes = (allNuclMinAxis > 3) & (allNuclMinAxis < 10);
keeperIndexes = find(allowableIntensityIndexes & allowableAreaIndexes & allowableMaxAxisIndexes & allowableMinAxisIndexes);

% Extract only those objects that meet our criteria, and 
% eliminate those objects that don't meet our criteria. 
keeperNuclImage = ismember(labeledImage, keeperIndexes); 
% Re-label with only the keeper nucl kept.
labeledDimeImage = bwlabel(keeperNuclImage, 8);     % Label each nucl 

subplot(2, 3, 5); 
imshow(labeledDimeImage);
axis square;
title('"Keeper" objects (after "size" filter)');


keeperIndexes1=[];
keeperIndCol=keeperIndexes';
col=size(keeperIndCol);

for ch=1: col
    
  mayBE(ch,1:2) = nuclMeasurements(keeperIndexes(1,ch)).Centroid;
  mayBE(ch,3) = keeperIndexes(1,ch);
  
end


%Counting neighbors 


dist=30; %Distance to other cell

distB=10; %Distance to Image Border

%Correction for 40x center
corrX=13;
corrY=33;

% Filter for cell too close to each other
i=1;
for nn=1 : col
raz=0;
fak=0;
   for mm=1: numberOfNucl
       znachenieX=abs(mayBE(nn,1)- allNuclCentroids(1,mm+raz));
       znachenieY=abs(mayBE(nn,2)- allNuclCentroids(1,mm+raz+1));
       
       
      
       if (znachenieX>0)||(znachenieY>0)
           
           if (znachenieX <dist)&&(znachenieY<dist)
           fak=1;
           
           break
          
           end
       end
        raz=raz+1;
     
   end
  
% Filter for cell too close to image border
   
 if (fak==0)
     
     distFromImageBorderX=SizeX-mayBE(nn,1);
     distFromImageBorderY=SizeY-mayBE(nn,2);
     
     if (distFromImageBorderX>distB) && (distFromImageBorderX<SizeX-distB)&&(distFromImageBorderY>distB) && (distFromImageBorderX<SizeY-distB)
    
     keeperIndexes1(1,i)=mayBE(nn,3);
     positions(i,1:2)= mayBE(nn,1:2);
     
% Keep positions of cells that passed filters

     interX= LeftUpPositionX - positions(i,1)*SizeOnePixX - corrX; 
     interY= positions(i,2)*SizeOnePixY + LeftUpPositionY - corrY;
     positions(i,3) = interX/1000000;
     positions(i,4) = interY/1000000;
     
         i=i+1;
     end   
 end
end

% If position list already exist add new positions, 
% otherwise create new position list

if exist('positions', 'var')
    if exist('FinalPositions', 'var')
FinalPositions=[FinalPositions ; positions];
    else
        FinalPositions=positions;
    end
end
keeperNuclImage= ismember(labeledImage, keeperIndexes1);
% Re-label with only the keeper nucl.
labeledDimeImage= bwlabel(keeperNuclImage,8);     % Label each nucl 
subplot(2, 3, 6); 
imshow(labeledDimeImage);
axis square;
title('Objects after "neighbor" filter');




end
end
end

% Create position list in special format (for microscope)

if exist('FinalPositions', 'var')
    
    head='<?xml version="1.0" encoding="UTF-8"?>';
    head2='<MarkAndFindPointList xmlns:xsd="http://www.w3.org/2001/XMLSchema" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">';
    footer='</MarkAndFindPointList>';
    body='';
    PositionZ=0.0024502658; % needs to be corrected for every microscope
    PositionZ = num2str(PositionZ);
    for iii=1 : size (FinalPositions)
      PositionY=num2str(FinalPositions(iii,4));
      PositionX=num2str(FinalPositions(iii,3));
      body1=sprintf('<MarkAndFindPoint FieldRotation="0" FieldZCoordinate="%s" FieldYCoordinate="%s" FieldXCoordinate="%s"/>\r\n', PositionZ,PositionY,PositionX);  
       body=sprintf('%s%s', body, body1); 
    end
    
 content=sprintf('%s\r\n%s\r\n%s%s', head , head2 , body , footer);    

 % Save position list on the hard drive
 
 FileNameXmlOut=sprintf('Positions-%s-%d.xml', exp, iii); 
fid=fopen(FileNameXmlOut,'w' ); 
fprintf(fid, content); 
fclose=(fid);

disp('File Name:');
disp (FileNameXmlOut);
else
  disp ('Sorry no any single position was found');  
end

    
        
        
end