Applying Regionprops to Only the Region of Interest

I would like to use the "regionprops" function on a user-defined ROI from an image to generate a histogram that plots the frequencies of the diameters of the circular objects ONLY within the user-defined ROI. The image, labeled, "Image_No_ROI_Defined", is the output from my code which outlines circular objects from the entire image using "regionprops" and "viscircles" documentation. I had gotten to the point where my code allows the user to draw their ROI from this image and then marks the centroids of the circular objects within the user-defined ROI (See attached "Image_ROI_Defined" image file). Is there any possible way I would be able to use any functions to outline the circles within the user-defined ROI in a different color and then measure the diameters of the circles only within the yellow "ROI" region that have their centroids marked with the yellow crosses? I have attached the original image, the image with all circles outlined, and the image with the ROI I would only like regionprops to analyze to calculate diameter of the circles within the ROI.

Any help would be greatly appreciated. Thank you.

Answer 1

[Image Analyst]

Since it looks like you want even partial circles and the only criteria is that the circle centers lie within the polygon the user drew, I'd use inpolygon() to determine which measurements are inside. So if you have xCenters and yCenters arrays, and arrays of vertices, xv and yv, i'd do this

keeperIndexes = false(1, length(xCenters));
for k = 1 : length(xCenters)
    if inpolygon(xCenters(k), yCenters(k), xv(k), yv(k))
        keeperIndexes(k) = true;
    end
end
% Now extract only those that are inside the polygon
roix = xCenters(keeperIndexes);
roiy = yCenters(keeperIndexes);
roiRadii = radii(keeperIndexes);

[Michael DiStefano]

Thank you for your response. I am receiving an error: "Index exceeds the number of array elements (4)." Would you like me the portion of my code which incorporates your code?

[Image Analyst]

Of course.

[Michael DiStefano]

Below is the code. I have attached the "I3" image to this message just as reference. Thank you.

% Portion of code that marks all centroids and outlines all circular objects
cc = bwconncomp(I3);
stats = regionprops(cc,'Area');
idx = find([stats.Area] > 700 );
bw6 = ismember(labelmatrix(cc), idx);
stats2 = regionprops(bw6,'Centroid', 'MinorAxisLength', 'MajorAxisLength');
centroids = cat(1,stats2.Centroid);
Cx = centroids(:,1);
Cy = centroids(:,2);
diameters = mean([cell2mat({stats2.MajorAxisLength})' cell2mat({stats2.MinorAxisLength})'],2);
radii = diameters/2;
imshow(imread(file))
hold on
plot(Cx,Cy,'b*')
viscircles(centroids, radii);

% User-defined ROI Polygon
[ROI x y] = roipoly;

% Analysis of ROI
keeperIndexes = false(1, length(Cx));
for k = 1 : length(Cx)
    if inpolygon(Cx(k), Cy(k), x(k), y(k))
        keeperIndexes(k) = true;
    end
end
% Now extract only those that are inside the polygon
roix = Cx(keeperIndexes);
roiy = Cy(keeperIndexes);
roiRadii = radii(keeperIndexes);

[Michael DiStefano]

When I look at the vertex arrays, it appears as though the coordinates of the first vertex are repeated again as the final vertex since I am double clicking (or right clicking and selecting "Create Mask" which does the same thing) the first vertex to complete the polygon. Could this be the reason as to why there is an error?

[Image Analyst]

Michael, it doesn't work. I spent 5 minutes trying to fix it but failed. Please fix it and make sure it does what you expect with the image you posted.

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 15;
fprintf('Beginning to run %s.m ...\n', mfilename);

I3 = imread('I3.png');
if ndims(I3) > 1
	I3 = I3(:,:,1);
end
% Binarize.
I3 = I3 > 128;
% Portion of code that marks all centroids and outlines all circular objects
cc = bwconncomp(I3);
stats = regionprops(cc,'Area')
allAreas = sort([stats.Area], 'descend')
idx = find([stats.Area] > 700 );
bw6 = ismember(labelmatrix(cc), idx);
stats2 = regionprops(bw6,'Centroid', 'MinorAxisLength', 'MajorAxisLength');
centroids = cat(1,stats2.Centroid);
Cx = centroids(:,1);
Cy = centroids(:,2);
diameters = mean([cell2mat({stats2.MajorAxisLength})' cell2mat({stats2.MinorAxisLength})'],2);
radii = diameters/2;
% imshow(imread(file)) % Doesn't work
imshow(I3) % Attempted fix
hold on
plot(Cx,Cy,'b*')
viscircles(centroids, radii);

% User-defined ROI Polygon
[ROI x y] = roipoly;

% Analysis of ROI
keeperIndexes = false(1, length(Cx));
for k = 1 : length(Cx)
    if inpolygon(Cx(k), Cy(k), x(k), y(k))
        keeperIndexes(k) = true;
    end
end
% Now extract only those that are inside the polygon
roix = Cx(keeperIndexes);
roiy = Cy(keeperIndexes);
roiRadii = radii(keeperIndexes);

fprintf('Done running %s.m ...\n', mfilename);
msgbox('Done!');

[Michael DiStefano]

I apologize for the confusion. I have posted my entire code below. Just to note, I have been analyzing large .tif images which I why I posted .PNG images. The "Original Image.PNG" file in my first post is the image I start with for my code. Thank you again for your help.

clear all; close all; clc;

[file, path] = uigetfile('*.tif','Select TEM Image of Interest');
figure, imshow(file), title('Original TEM Image');
TEM_Image = imadjust(imread(file));

Gauss = imgaussfilt(TEM_Image,6, 'padding','circular');
figure, imshow(Gauss)

Otsu = imbinarize(Gauss);
figure, imshow(Otsu)

Otsu_I = imcomplement(Otsu);
figure, imshow(Otsu_I)

L = watershed(Otsu_I);
Lrgb = label2rgb(L);
imshow(Lrgb)

imshow(imfuse(Otsu_I,Lrgb))
Otsu_I2 = ~bwareaopen(~Otsu_I, 100);
imshow(Otsu_I2)

D = -bwdist(~Otsu_I);
imshow(D,[])

Ld = watershed(D);
imshow(label2rgb(Ld))

Otsu_I2 = Otsu_I;
Otsu_I2(Ld == 0) = 0;
imshow(Otsu_I2)

mask = imextendedmin(D,8);
figure, imshowpair(Otsu_I,mask,'blend')

D2 = imimposemin(D,mask);
Ld2 = watershed(D2);
Otsu_I3 = Otsu_I;
Otsu_I3(Ld2 == 0) = 0;
imshow(Otsu_I3)


cc = bwconncomp(Otsu_I3);
stats = regionprops(cc,'Area');
idx = find([stats.Area] > 700 );
bw6 = ismember(labelmatrix(cc), idx);
stats2 = regionprops(bw6,'Centroid', 'MinorAxisLength', 'MajorAxisLength');
centroids = cat(1,stats2.Centroid);
Cx = centroids(:,1);
Cy = centroids(:,2);
diameters = mean([cell2mat({stats2.MajorAxisLength})' cell2mat({stats2.MinorAxisLength})'],2);
radii = diameters/2;
imshow(imread(file))
hold on
plot(Cx,Cy,'b*')
viscircles(centroids, radii);

[ROI x y] = roipoly;


keeperIndexes = false(1, length(Cx));
for k = 1 : length(Cx)
    if inpolygon(Cx(k), Cy(k), x(k), y(k))
        keeperIndexes(k) = true;
    end
end
% Now extract only those that are inside the polygon
roix = Cx(keeperIndexes);
roiy = Cy(keeperIndexes);
roiRadii = radii(keeperIndexes);

[Image Analyst]

Still needs fixing. I3 is color, not gray scale, and other things. Let me work on it.

[Michael DiStefano]

Thank you. Just a quick note, I3 is not as important as the "Original Image.PNG" which is the actual gray scale image I use at the beginning of my code. I just included the image I3 becuase that was the image used for the later part of my code.

[Image Analyst]

I fixed tons of errors and made it more robust and informative. I got this but I'm not sure its the same as what you got. Try running this much improved version and let me know if something is not correct. What is your algorithm supposed to find? It's certainly not finding the center of every circle because your circles are connected.

See Steve's blog:

Watershed transform question from tech support » Watershed transform question from tech support

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 15;
fprintf('Beginning to run %s.m ...\n', mfilename);

[baseFileName, folder] = uigetfile('*.png','Select TEM Image of Interest');
if baseFileName == 0
	return;
end
fullFileName = fullfile(folder, baseFileName)
grayImage = imread(fullFileName);
[rows, columns, numberOfColorChannels] = size(grayImage)
if numberOfColorChannels > 1
	grayImage = rgb2gray(grayImage);
end
subplot(2, 3, 1);
imshow(fullFileName);
title('Original TEM Image');
TEM_Image = imadjust(grayImage);

Gauss = imgaussfilt(TEM_Image,6, 'padding','circular');
subplot(2, 3, 2);
imshow(Gauss)
title('Gauss Image');

Otsu = imbinarize(Gauss);
subplot(2, 3, 3);
imshow(Otsu)
title('Otsu Image');

Otsu_I = imcomplement(Otsu);
subplot(2, 3, 4);
imshow(Otsu_I)
title('Otsu_I Image');

L = watershed(Otsu_I);
Lrgb = label2rgb(L);
subplot(2, 3, 5);
imshow(Lrgb)
title('Labeled Image');

subplot(2, 3, 6);
imshow(imfuse(Otsu_I,Lrgb))
Otsu_I2 = ~bwareaopen(~Otsu_I, 100);
figure;
subplot(2, 3, 1);
imshow(Otsu_I2)
title('Otsu_I2 Image');

D = -bwdist(~Otsu_I);
subplot(2, 3, 2);
imshow(D,[])
title('D Image');

Ld = watershed(D);
subplot(2, 3, 3);
imshow(label2rgb(Ld))
title('Ld Image');

Otsu_I2 = Otsu_I;
Otsu_I2(Ld == 0) = 0;
subplot(2, 3, 4);
imshow(Otsu_I2)
title('Otsu_I2 Image');

mask = imextendedmin(D,8);
subplot(2, 3, 5);
imshowpair(Otsu_I,mask,'blend')
title('Otsu_I and Mask Image');

D2 = imimposemin(D,mask);
Ld2 = watershed(D2);
Otsu_I3 = Otsu_I;
Otsu_I3(Ld2 == 0) = 0;
subplot(2, 3, 6);
imshow(Otsu_I3)
title('Otsu_I3 Image');

cc = bwconncomp(Otsu_I3);
stats = regionprops(cc,'Area');
idx = find([stats.Area] > 700 );
bw6 = ismember(labelmatrix(cc), idx);
stats2 = regionprops(bw6,'Centroid', 'EquivDiameter');
centroids = cat(1,stats2.Centroid);
Cx = centroids(:,1);
Cy = centroids(:,2);
diameters = [stats2.EquivDiameter];
radii = diameters/2;
figure;
imshow(grayImage)
hold on
hp = plot(Cx,Cy,'b*')
viscircles(centroids, radii);
g = gcf;
g.WindowState = 'maximized';

title('Draw a polygon.  Left click vertices, then right click to finish it.');
uiwait(helpdlg('Draw a polygon.  Left click vertices, then right click to finish it.'));
[mask, xPoly, yPoly] = roipolyold();

% Erase outside the polygon.
grayImage(~mask) = 0;
hold off;
cla;
imshow(grayImage)
hold on
hp = plot(Cx, Cy, 'b*', 'MarkerSize', 25, 'LineWidth', 2)

keeperIndexes = false(1, length(Cx));
for k = 1 : length(Cx)
    if inpolygon(Cx(k), Cy(k), xPoly, yPoly)
        keeperIndexes(k) = true;
    end
end
% Now extract only those that are inside the polygon
roix = Cx(keeperIndexes)
roiy = Cy(keeperIndexes)
roiRadii = radii(keeperIndexes)
hp = plot(roix, roiy,'c*', 'MarkerSize', 25, 'LineWidth', 2)
viscircles([roix, roiy], roiRadii);

fprintf('Done running %s.m ...\n', mfilename);
msgbox('Done!');

[Michael DiStefano]

It looks like it's working very well! Thank you and I greatly appreciate your help! The goal of the code is to generate a histogram with the frequency of the diameters from the circles in the user-defined ROI. The diameters must be expressed in nanometers.

[Image Analyst]

Well that code doesn't find all the circles but maybe if you start with the actual original image instead of the binary image it works better. If it works fine for you as is, then click the "Accept this Answer" link. Otherwise, let's fix it.

[Michael DiStefano]

Is there any possible way that all "edge circles" within the ROI can be included or will that get too involved? This is the final image I am getting when I use the original gray scale image.

[Image Analyst]

Cool. Looks 3-D even though it's not (due to the longitudinal achromatic aberration in the human eye where blue has a different focal length than red).

You can get the equation of each line in your polygon. Then for each point that is inside the polygon, see if its "point-line distance" is less than the radius for the circle for that point. If it is, then it's too close and some part of the circle sticks outside the polygon, and then you can exclude that one.

See attached point-line distance m-file.

[Michael DiStefano]

Okay great, thanks again! I just wanted to run one thing by you regarding the conversion from pixels to nanometers. At the bottom left corner of each TEM image, there is a 200nm scale bar. I can use this scale bar in your spatial_calibration_demo.m code to determine the approximated conversion factor correct?

[Image Analyst]

Yes. If your original question was answered, could you please "Accept thsi answer"? Thanks in advance.

[Michael DiStefano]

Hello Image Analyst. Just out of curiosity, is there a way I can allow the user to select the centroids of the circles to exclude from further analysis? I've added an additional part of the code that gives the user the option to exclude a region using ~inpolygon, however, I believe that giving the user the option to just click on the centroids on circles to exclude them from further anaylysis would expedite this process. This would probably eliminate the need to draw the ROI for further analysis as well since the circles to be excluded will have had their centroids selected by the user.

I can place this question in a new discussion if you would like. Thanks again for your help so far.

[Image Analyst]

If you have the list of centers in (x,y), call ginput() to let them pick points, then I'd loop to find the closest circle to each clicked point.

message = sprintf('Click points and hit return when done entering all points.');
uiwait(helpdlg(message));
[xUser, yUser] = ginput(); % I think you hit return when done entering all points.
indexesToExclude = zeros(1, length(xCenters));
minDistance = zeros(1, length(xCenters));
for k = 1 : length(xUser)
    allDistances = sqrt((xUser(k) - xCenters).^2 + (yUser(k) - yCenters).^2);
    [minDistance(k), indexesToExclude(k)] = min(allDistances)
end

[Michael DiStefano]

Thank you again. I can include an "if" statement within the for loop to exclude the selected circles, based upon if the "minDistance" is less than the radius of each circle, from the plot and analysis or is this not necessary?

[Image Analyst]

Not necessary. Somehow you basically need to get a list of indexes either that

1. you want to exclude from analysis
2. you want to include in analysis

Of course those two sets of indexes are just inverses/complements of each other so you just need to do one.

So you can either

1. do your subsequent processing in that loop on the ones you want and skipping the ones you don't want with an if block, or
2. extract just the objects you want after the first loop, and then have another loop over only those remaining "good" ones.

[Michael DiStefano]

Okay. Just out of curiosity, I understand you said it wasn't necessary to include "minDistance(k) < radii(k)" in the "if" statement, but I am a bit confused as to why the code below wouldn't work properly. Although the number of circles I selected is correct, the code is excluding circles that I did not select. Below is the code and I've also attached the exported image.

message = sprintf('Click points and hit return when done entering all points.');
uiwait(helpdlg(message));
[xUser, yUser] = ginput(); 
indexesToExclude = zeros(1, length(Cx));
minDistance = zeros(1, length(Cx));
for k = 1 : length(xUser)
    allDistances = sqrt((xUser(k) - Cx).^2 + (yUser(k) - Cy).^2);
    [minDistance(k), indexesToExclude(k)] = min(allDistances); 
        if minDistance(k) < radii(k);
            indexesToExclude(k) = true;
        end    
end
roix = Cx(~indexesToExclude);
roiy = Cy(~indexesToExclude);
roiRadii = radii(~indexesToExclude);
hp = plot(roix, roiy,'y*', 'MarkerSize', 10, 'LineWidth', 2)
viscircles([roix, roiy], roiRadii);

[Image Analyst]

Please give me code that runs. This code below does not run because the Cx, Cy, and Radii vectors are missing. Perhaps attach a mat file with those in it.

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format short g;
format compact;
fprintf('Beginning to run %s.m ...\n', mfilename);

grayImage = imread('ExcludeExport.PNG');
imshow(grayImage);
message = sprintf('Click points and hit return when done entering all points.');
uiwait(helpdlg(message));
[xUser, yUser] = ginput();
indexesToExclude = zeros(1, length(Cx));
minDistance = zeros(1, length(Cx));
for k = 1 : length(xUser)
	allDistances = sqrt((xUser(k) - Cx).^2 + (yUser(k) - Cy).^2);
	[minDistance(k), indexesToExclude(k)] = min(allDistances);
	if minDistance(k) < radii(k);
		indexesToExclude(k) = true;
	end
end
roix = Cx(~indexesToExclude);
roiy = Cy(~indexesToExclude);
roiRadii = radii(~indexesToExclude);
hp = plot(roix, roiy,'y*', 'MarkerSize', 10, 'LineWidth', 2)
viscircles([roix, roiy], roiRadii);

fprintf('Done running %s.m.\n', mfilename);

[Michael DiStefano]

I apologize for that. Below is the entire code. I was also curious if there is a way that a small marker can be placed at each of the points that is selected using the ginput function so that the user does not lose track of which circles are to be excluded. Thank you again for your help.

clear all; close all; clc;
format long g;
format compact;
fontSize = 15;
fprintf('Beginning to run %s.m ...\n', mfilename);
[baseFileName, folder] = uigetfile('*.tif','Select TEM Image of Interest');
if baseFileName == 0
	return;
end
fullFileName = fullfile(folder, baseFileName)
grayImage = imread(fullFileName);
[rows, columns, numberOfColorChannels] = size(grayImage)
if numberOfColorChannels > 1
	grayImage = rgb2gray(grayImage);
end

subplot(2, 3, 1);
imshow(fullFileName);
title('Original TEM Image');
TEM_Image = imadjust(grayImage);
Gauss = imgaussfilt(TEM_Image,6, 'padding','circular');
subplot(2, 3, 2);
imshow(Gauss)
title('Gauss Image');
Otsu = imbinarize(Gauss);
subplot(2, 3, 3);
imshow(Otsu)
title('Otsu Image');
Otsu_I = imcomplement(Otsu);
subplot(2, 3, 4);
imshow(Otsu_I)
title('Otsu_I Image');
L = watershed(Otsu_I);
Lrgb = label2rgb(L);
subplot(2, 3, 5);
imshow(Lrgb)
title('Labeled Image');
subplot(2, 3, 6);
imshow(imfuse(Otsu_I,Lrgb))
Otsu_I2 = ~bwareaopen(~Otsu_I, 100);
figure;
subplot(2, 3, 1);
imshow(Otsu_I2)
title('Otsu_I2 Image');
D = -bwdist(~Otsu_I);
subplot(2, 3, 2);
imshow(D,[])
title('D Image');
Ld = watershed(D);
subplot(2, 3, 3);
imshow(label2rgb(Ld))
title('Ld Image');
Otsu_I2 = Otsu_I;
Otsu_I2(Ld == 0) = 0;
subplot(2, 3, 4);
imshow(Otsu_I2)
title('Otsu_I2 Image');
mask = imextendedmin(D,8);
subplot(2, 3, 5);
imshowpair(Otsu_I,mask,'blend')
title('Otsu_I and Mask Image');
D2 = imimposemin(D,mask);
Ld2 = watershed(D2);
Otsu_I3 = Otsu_I;
Otsu_I3(Ld2 == 0) = 0;
subplot(2, 3, 6);
imshow(Otsu_I3)
title('Otsu_I3 Image');

cc = bwconncomp(Otsu_I3);
stats = regionprops(cc,'Area');
idx = find([stats.Area] > 600 );
bw6 = ismember(labelmatrix(cc), idx);
stats2 = regionprops(bw6,'Centroid', 'MinorAxisLength','MajorAxisLength');
centroids = cat(1,stats2.Centroid);
Cx = centroids(:,1);
Cy = centroids(:,2);
diameters = mean([cell2mat({stats2.MajorAxisLength})' cell2mat({stats2.MinorAxisLength})'],2);
radii = diameters/2;
figure;
imshow(grayImage)
hold on
hp = plot(Cx,Cy,'b*')
viscircles(centroids, radii);

message = sprintf('Click points and hit return when done entering all points.');
uiwait(helpdlg(message));
[xUser, yUser] = ginput(); 
indexesToExclude = false(1, length(Cx));
minDistance = false(1, length(Cx));
for k = 1 : length(xUser)
    allDistances = sqrt((xUser(k) - Cx).^2 + (yUser(k) - Cy).^2);
    [minDistance(k), indexesToExclude(k)] = min(allDistances); 
        if minDistance(k) < radii(k);
            indexesToExclude(k) = true;
        end    
end
roix = Cx(~indexesToExclude);
roiy = Cy(~indexesToExclude);
roiRadii = radii(~indexesToExclude);
hp = plot(roix, roiy,'y*', 'MarkerSize', 10, 'LineWidth', 2)
viscircles([roix, roiy], roiRadii);

[Image Analyst]

Michael, please let me help you. The image you attached original_image.png, does not seem to work:

It's not finding the circles or centers. Please attach an image that will work.

[Michael DiStefano]

I've been using .tif extension images and have been sending you my screenshots of these images in .png format. For some reason, the code is not working for .png images but I have attached what it looks like for a .tif image I have been using.

[Michael DiStefano]

It seems as though the code is finishing much faster for the .png image than for the .tif image. Based on your expertise, what do you think this could be attributed to and is this the reason why the code is having difficulty in detecting the centroids and outlining the circles from the .png image as opposed to the .tif image? Thank you for your help so far.

[Image Analyst]

I don't know. There should be no difference except for maybe the time to read it in from disk and decompress it. But you forgot to attach a tiff file so I can't compare them.

[Michael DiStefano]

I am having trouble submitting the zipped form of the .tif file becuase the zip exceeds 5 MB. I am confused becuase in the Help section, it says the image should not exceed 15 MB but it is less than 15 MB. I apologize for all of this confusion but I cannot seem to be able to upload the image. Do you have any suggestions?

[Image Analyst]

Crop out the upper left quadrant of it.

[Michael DiStefano]

Due to the size of the file, I had to crop out a majority of the original image but this should do the trick. I've run the code on the image I have attached and it works decently. The only problem is that the circles on the border of the image are not outlined as effectively as they are in the actual (non-cropped) image but that is not my major concern at the moment. I am re-attaching the code below. Just to reiterate the problem I need some assistance with, although the number of circles I selected is correct, the code is excluding circles that I did not select with the "ginput" function. Thank you for your help and I apologize for the confusion with the images not working properly.

clear all; close all; clc;
format long g;
format compact;
fontSize = 15;
fprintf('Beginning to run %s.m ...\n', mfilename);

[baseFileName, folder] = uigetfile('*.tif','Select TEM Image of Interest');
if baseFileName == 0
	return;
end
fullFileName = fullfile(folder, baseFileName)
grayImage = imread(fullFileName);
[rows, columns, numberOfColorChannels] = size(grayImage)
if numberOfColorChannels > 1
	grayImage = rgb2gray(grayImage);
end

subplot(2, 3, 1);
imshow(fullFileName);
title('Original TEM Image');
TEM_Image = imadjust(grayImage);
Gauss = imgaussfilt(TEM_Image,6, 'padding','circular');
subplot(2, 3, 2);
imshow(Gauss)
title('Gauss Image');
Otsu = imbinarize(Gauss);
subplot(2, 3, 3);
imshow(Otsu)
title('Otsu Image');
Otsu_I = imcomplement(Otsu);
subplot(2, 3, 4);
imshow(Otsu_I)
title('Otsu_I Image');
L = watershed(Otsu_I);
Lrgb = label2rgb(L);
subplot(2, 3, 5);
imshow(Lrgb)
title('Labeled Image');
subplot(2, 3, 6);
imshow(imfuse(Otsu_I,Lrgb))
Otsu_I2 = ~bwareaopen(~Otsu_I, 100);
figure;
subplot(2, 3, 1);
imshow(Otsu_I2)
title('Otsu_I2 Image');
D = -bwdist(~Otsu_I);
subplot(2, 3, 2);
imshow(D,[])
title('D Image');
Ld = watershed(D);
subplot(2, 3, 3);
imshow(label2rgb(Ld))
title('Ld Image');
Otsu_I2 = Otsu_I;
Otsu_I2(Ld == 0) = 0;
subplot(2, 3, 4);
imshow(Otsu_I2)
title('Otsu_I2 Image');
mask = imextendedmin(D,8);
subplot(2, 3, 5);
imshowpair(Otsu_I,mask,'blend')
title('Otsu_I and Mask Image');
D2 = imimposemin(D,mask);
Ld2 = watershed(D2);
Otsu_I3 = Otsu_I;
Otsu_I3(Ld2 == 0) = 0;
subplot(2, 3, 6);
imshow(Otsu_I3)
title('Otsu_I3 Image');

cc = bwconncomp(Otsu_I3);
stats = regionprops(cc,'Area');
idx = find([stats.Area] > 600 );
bw6 = ismember(labelmatrix(cc), idx);
stats2 = regionprops(bw6,'Centroid', 'MinorAxisLength','MajorAxisLength');
centroids = cat(1,stats2.Centroid);
Cx = centroids(:,1);
Cy = centroids(:,2);
diameters = mean([cell2mat({stats2.MajorAxisLength})' cell2mat({stats2.MinorAxisLength})'],2);
radii = diameters/2;
figure;
imshow(grayImage)
hold on
hp = plot(Cx,Cy,'b*')
viscircles(centroids, radii);

message = sprintf('Click points and hit return when done entering all points.');
uiwait(helpdlg(message));
[xUser, yUser] = ginput(); 
indexesToExclude = false(1, length(Cx));
minDistance = false(1, length(Cx));
for k = 1 : length(xUser)
    allDistances = sqrt((xUser(k) - Cx).^2 + (yUser(k) - Cy).^2);
    [minDistance(k), indexesToExclude(k)] = min(allDistances); 
        if minDistance(k) < radii(k);
            indexesToExclude(k) = true;
        end    
end
roix = Cx(~indexesToExclude);
roiy = Cy(~indexesToExclude);
roiRadii = radii(~indexesToExclude);
hp = plot(roix, roiy,'y*', 'MarkerSize', 10, 'LineWidth', 2)
viscircles([roix, roiy], roiRadii);

[Image Analyst]

Do you want circles if their centroid is within the user-drawn polygon?

Or only if the entire circle is within the polygon?

If the first, then you can use inpolygon().

If the second, you have to use binary image operations.

[Michael DiStefano]

I have similar code that uses inpolygon() to either outline an ROI for further analysis or exclude an ROI to prevent it from being included in further analysis. To make it even easier for the user and to also allow the user to select circles that should not be included in further analysis, I was wondering if there was a way to allow the user to select the circular objects to exclude from further analysis using ginput(), as you had previously mentioned to me through this code below:

message = sprintf('Click points and hit return when done entering all points.');
uiwait(helpdlg(message));
[xUser, yUser] = ginput(); % I think you hit return when done entering all points.
indexesToExclude = zeros(1, length(xCenters));
minDistance = zeros(1, length(xCenters));
for k = 1 : length(xUser)
    allDistances = sqrt((xUser(k) - xCenters).^2 + (yUser(k) - yCenters).^2);
    [minDistance(k), indexesToExclude(k)] = min(allDistances)
end

However, the code in my previous message, which incorporates this portion of code you had sent, attempts to do that but there is an issue with the last for loop in that it does not exclude the circles that the user selected. I was also considering writing a script that allows the user to outline multiple ROIs to either include or exclude but I thought that selecting circles would be more efficient. If there is any way you can help recommend how I can fix the mistake in my for loop in my previous message so that the user can just select circles to exclude, I would greatly appreciate it. Thank you again.

[Image Analyst]

That will only compare distances to where you clicked, not the entire outline. I think you should use imfreehand() instead of ginput(). See attached demos.