MATLAB Answers

finding the radii of circles knowing the center of those circles

38 views (last 30 days)
I have a 2D hologram which has 25 circles of diffrent radii. I have removed the noise and binarized the image and then I found the centre points of the circles.
Now, i need to finally find the radii of those circles without using advanced functions

Accepted Answer

Image Analyst
Image Analyst on 22 Jan 2021
Edited: Image Analyst on 22 Jan 2021
You can do this:
mask = grayImage > 50; % Or whatever
props = regionprops(mask, 'Centroid', 'EquivDiameter', 'Area');
allDiameters = [props.EquivDiameter]
allAreas = [props.Area]
xy = vertcat(props.Centroid)
% Plot centroid over the image.
for k = 1 : length(props)
txt = sprintf('#%d at (x, y)\n= (%.1f, %.1f)', k, xy(k, 1), xy(k, 2));
text(xy(k, 1), xy(k, 2), txt, 'Color', 'r', 'FontSize', 8, 'FontWeight', 'bold', 'HorizontalAlignment', 'center');
end
Below, and attached, is the full demo:
% Code to find the centroids of spots.
clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
imtool close all; % Close all imtool figures if you have the Image Processing Toolbox.
clear; % Erase all existing variables. Or clearvars if you want.
workspace; % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 16;
%===============================================================================
% Read in gray scale image.
folder = pwd;
baseFileName = 'Cropping_initial.jpg';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
if ~exist(fullFileName, 'file')
% Didn't find it there. Check the search path for it.
fullFileName = baseFileName; % No path this time.
if ~exist(fullFileName, 'file')
% Still didn't find it. Alert user.
errorMessage = sprintf('Error: %s does not exist.', fullFileName);
uiwait(warndlg(errorMessage));
return;
end
end
rgbImage = imread(fullFileName);
% Get the dimensions of the image. numberOfColorBands should be = 1.
[rows, columns, numberOfColorBands] = size(rgbImage);
% If it's RGB instead of grayscale, convert it to gray scale.
if numberOfColorBands > 1
grayImage = rgbImage(:, :, 1); % Take red channel.
else
grayImage = rgbImage;
end
% Display the original image.
subplot(2, 2, 1);
imshow(grayImage);
axis on;
impixelinfo; % Let user mouse around and see gray level.
caption = sprintf('Original Image : %s', baseFileName);
title(caption, 'FontSize', fontSize, 'Interpreter', 'none');
impixelinfo;
% Enlarge figure to full screen.
g = gcf;
g.WindowState = 'maximized';
g.NumberTitle = 'off';
g.Name = 'Demo by Image Analyst'
drawnow;
%===============================================================================
% Display the histogram
subplot(2, 2, 2);
imhist(grayImage);
grid on;
title('Original Gray Scale Image Histogram', 'FontSize', fontSize);
xlabel('Gray Level', 'FontSize', fontSize);
ylabel('Pixel Count', 'FontSize', fontSize);
%===============================================================================
% IMAGE SEGMENTATION
% Specify a threshold. The image is 16 bits in the range 0-65535 so the threshold will be high.
threshold = 50;
% Place line on histogram at the mean.
xline(threshold, 'Color', 'r', 'LineWidth', 2);
yl = ylim;
txt = sprintf(' Threshold = %d', threshold);
text(threshold, yl(1) + 0.76 * (yl(2)-yl(1)), txt, 'Color', 'r', 'FontSize', 18, 'FontWeight', 'bold', 'HorizontalAlignment', 'left');
% Create a binary image
mask = grayImage > threshold;
% Get rid of white surround or any blob touching the border (because it's not a complete blob).
mask = imclearborder(mask);
% Get rid of any blobs less than 1000 in area (there is one blob now that has an area of 1 and the rest are around 7000-9000).
mask = bwareaopen(mask, 1000);
subplot(2, 2, 3);
imshow(mask);
title('Final Mask', 'FontSize', fontSize);
drawnow;
%===============================================================================
% IMAGE ANALYSIS
props = regionprops(mask, 'Centroid', 'EquivDiameter', 'Area');
allDiameters = [props.EquivDiameter]
allAreas = [props.Area]
xy = vertcat(props.Centroid)
% Plot centroid over the image.
for k = 1 : length(props)
txt = sprintf('#%d at (x, y)\n= (%.1f, %.1f)', k, xy(k, 1), xy(k, 2));
text(xy(k, 1), xy(k, 2), txt, 'Color', 'r', 'FontSize', 8, 'FontWeight', 'bold', 'HorizontalAlignment', 'center');
end
msgbox('Done! Thank you Image Analyst!');

  6 Comments

Show 3 older comments
Image Analyst
Image Analyst on 24 Jan 2021
What I'd do is to use improfile() to get a profile through the centers. Then I'd use the FAQ
to fit the data to a circle. The data for one lens would go from the center outwards until the intensity leveled off. Or you could just go from the peak down to about 10% above the min value. That should be enough points to fit a circle to. The code from the FAQ will give you the radius. Do you think you can do it?
ARIJIT GANGULY
ARIJIT GANGULY on 24 Jan 2021
Nope. Can you edit it and add to the mat code(radii detection) I sent you.
Image Analyst
Image Analyst on 24 Jan 2021
Here's a start
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 = 20;
%--------------------------------------------------------------------------------------------------------
% READ IN IMAGE
folder = pwd;
baseFileName = 'cropping_initial.jpg';
grayImage = imread(baseFileName);
% Get the dimensions of the image.
% numberOfColorChannels should be = 1 for a gray scale image, and 3 for an RGB color image.
[rows, columns, numberOfColorChannels] = size(grayImage)
if numberOfColorChannels > 1
% It's not really gray scale like we expected - it's color.
% Use weighted sum of ALL channels to create a gray scale image.
grayImage = min(grayImage, [], 3);
end
%--------------------------------------------------------------------------------------------------------
% Display the image.
subplot(2, 1, 1);
imshow(grayImage, []);
axis('on', 'image');
title('Gray Scale Image', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;
hFig = gcf;
hFig.WindowState = 'maximized'; % May not work in earlier versions of MATLAB.
drawnow;
uiwait(msgbox('Click across a lens. Left click first point and right click second point.'));
[cx, cy, values] = improfile;
hold on;
plot(cx, cy, 'g-', 'LineWidth', 2);
% Convert to a line across the shape.
xAll = sqrt((cx - cx(1)) .^ 2 + (cy - cy(1)) .^ 2);
yAll = values;
subplot(2, 1, 2);
plot(xAll, yAll, 'b-', 'LineWidth', 2);
grid on;
% Fit a line between the endpoints.
linearCoefficients = polyfit([xAll(1), xAll(end)], [yAll(1), yAll(end)], 1);
yLineAcross = polyval(linearCoefficients, xAll);
hold on;
plot(xAll, yLineAcross, 'g-', 'LineWidth', 2);
% Get a first guess at where the circle is. Find out where the y value is above the line.
indexes = yAll > yLineAcross;
% Take the largest region only.
indexes = bwareafilt(indexes, 1);
xCircle = sqrt((cx(indexes) - cx(1)) .^ 2 + (cy(indexes) - cy(1)) .^ 2);
yCircle = values(indexes);
% Put lines up
plot(xCircle, yCircle, 'r-', 'LineWidth', 2);
grid on;
xline(xCircle(1), 'Color', 'r', 'LineWidth', 3);
xline(xCircle(end), 'Color', 'r', 'LineWidth', 3);
% That was a first guess. Now fit what's in between to a quadratic.
quadCoefficients = polyfit(xCircle, yCircle, 2);
yQuad = polyval(quadCoefficients, xAll);
hold on;
plot(xAll, yQuad, 'g-', 'LineWidth', 2);
% Find out where the y value is above the line.
indexes = yQuad > yLineAcross;
% Take the largest region only.
indexes = bwareafilt(indexes, 1);
xCircle = sqrt((cx(indexes) - cx(1)) .^ 2 + (cy(indexes) - cy(1)) .^ 2);
yCircle = values(indexes);
% Put lines up
hold on;
plot(xCircle, yCircle, 'r-', 'LineWidth', 2);
grid on;
xline(xCircle(1), 'Color', 'g', 'LineWidth', 3);
xline(xCircle(end), 'Color', 'g', 'LineWidth', 3);
% Get the radius
diameter = xCircle(end) - xCircle(1)
caption = sprintf('Diameter = %.1f, from %.1f to %.1f', diameter, xCircle(1), xCircle(end));
title(caption, 'FontSize', fontSize);

Sign in to comment.

More Answers (0)

Products


Release

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!