finding the radii of circles knowing the center of those circles

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

[Image Analyst]

What do you have in the digital (non-optical) domain? Do you have an digital image or something? MATLAB can't deal with a film - it needs digital data.

TUTORIAL: How to ask a question (on Answers) and get a fast answer

[ARIJIT GANGULY]

I have digital data.

Answer 1

[Image Analyst]

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

[ARIJIT GANGULY]

Hi, I showed it but I need to find the radius with direct use of hough transform or other things. I need to use loops not containing advanced functions. If you want I can send the code.

Can you correct the radii detection part in the code I have sent.

[Image Analyst]

Error using load

Unable to read file 'microlens_after_autofocusing.mat'. No such file or directory.

Error in radii_detection (line 1)

load microlens_after_autofocusing.mat

[ARIJIT GANGULY]

https://drive.google.com/file/d/1KqBVLnN3E8xUrzvgQQGOAtSru6xTch9e/view?usp=sharing

Link to the microlens image file. It's 75mb in size.

I have found the center points and saved them in the matrix and later placed the center points on the main image(As shown above with magenta center points). Now I need to find the radius of each circle and they are going to be 2-3 radius for each circle. I have average them and display radius, center points and the circle boundary from sobel (edge detected file) them in the main file. In sobel the circle are increased in size for better processing.

And finally display them as shown in the screenshot (104)

[Image Analyst]

What I'd do is to use improfile() to get a profile through the centers. Then I'd use the FAQ

https://matlab.fandom.com/wiki/FAQ#How_can_I_fit_a_circle_to_a_set_of_XY_data.3F

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]

Nope. Can you edit it and add to the mat code(radii detection) I sent you.

[Image Analyst]

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