How to create a circles (smallest and biggest circle) based on points in an image given, then find the center and radius?

Hi,

I have extract some points from Image Scanning, I need to create 2 circle (smallest circle and largest circle) in order to calculate the error of circularity (Rmax - Rmin), in order to do that I need to find the center and radius. Did any have done this kind of case?

Below are some image as reference:

1. Points extract from image scanning

2. Example of circles need to be done

Any help will be appreciated.

Thanks in advance

Answer 1

[Image Analyst]

To get the minimal bounding outer circle, try this:

John D'Errico's File Exchange

I asked John for the largest interior circle, and he said that's a much more difficult problem and he doesn't have code for that.

[Image Analyst]

Jafar, are you still there? Hopefully you were able to get John's code to work. But it seems like you must be having trouble getting the FAQ code to work, and finding the minimum interior circle, so I did it for you. See the code below where I first generate some random points, then fit a circle through them using code from the FAQ, then use the Euclidean Distance Transform to find the smallest interior circle that touches 3 points. Both the fitted circle and interior circle are shown plotted over the image. The fitted circle is the outer one obviously. Adapt as needed. If this meets your needs, please click the "Accept this answer" link.

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;

xCenter = 12;
yCenter = 10;
numPoints = 12;
theta = linspace(0, 2 * pi - 2 * pi / numPoints, numPoints);
radius = 5 * ones(1, numPoints) + 3 * rand(1, numPoints);
x = radius .* cos(theta) + xCenter;
y = radius .* sin(theta) + yCenter;
subplot(2, 2, 1);
plot(x, y, 'b+', 'MarkerSize', 15, 'LineWidth', 2);
axis square;
xlim([0 20]);
ylim([0 20]);
grid on;

% First fit a circle
[xcFit, ycFit, RFit, a] = circfit(x,y)

% Now find the max and min fit
% Convert to a binary image
% Scale by 500 to get enough pixels to get 3 digits of precision.
scaleFactor = 500;
xScaled = scaleFactor*x;
yScaled = scaleFactor*y;
rows = ceil(scaleFactor * max(y));
columns = ceil(scaleFactor * max(x));
binaryImage = poly2mask(xScaled, yScaled, rows, columns);
subplot(2, 2, 3);
imshow(binaryImage);
title('binaryImage', 'fontSize', fontSize);
axis('on', 'image');
hold on;
plot(xScaled, yScaled, 'b+', 'MarkerSize', 15, 'LineWidth', 2);

% Get the Euclidean Distance Transform
edtImage = bwdist(~binaryImage);
% Display the image.
subplot(1, 2, 2);
imshow(edtImage, []);
title('Distance Transform', 'fontSize', fontSize);
axis('on', 'image');
% Find the max
minRadius = max(edtImage(:))
[rowMin, colMin] = find(edtImage == minRadius, 1, 'first');
hold on;
% Plot the center of the minimum interior circle.
plot(colMin, rowMin, 'r+', 'MarkerSize', 20);
% Plot the scaled x, y
plot(xScaled, yScaled, 'b+', 'MarkerSize', 15, 'LineWidth', 2);
% Plot the circles
viscircles([xcFit * scaleFactor, ycFit * scaleFactor; colMin, rowMin], [RFit * scaleFactor; minRadius]);
% Unscale to original dimensions
xCenterMin = colMin / scaleFactor
yCenterMin = rowMin / scaleFactor
minRadius = minRadius / scaleFactor

% Print out the average/fitted circle, and the min circle.
fprintf('Fitted circle : Radius = %f at (x, y) = (%f, %f).\n', RFit, xcFit, ycFit);
fprintf('Min circle    : Radius = %f at (x, y) = (%f, %f).\n', minRadius, xCenterMin, yCenterMin);

function [xc,yc,R,a] = circfit(x,y)
%CIRCFIT  Fits a circle in x,y plane
%
% [XC, YC, R, A] = CIRCFIT(X,Y)
% Result is center point (yc,xc) and radius R.  A is an optional
% output describing the circle's equation:
%
%   x^2+y^2+a(1)*x+a(2)*y+a(3)=0
         
% by Bucher izhak 25/oct/1991

n=length(x);  xx=x.*x; yy=y.*y; xy=x.*y;
A=[sum(x) sum(y) n;sum(xy) sum(yy) sum(y);sum(xx) sum(xy) sum(x)];
B=[-sum(xx+yy) ; -sum(xx.*y+yy.*y) ; -sum(xx.*x+xy.*y)];
a=A\B;
xc = -.5*a(1);
yc = -.5*a(2);
R  =  sqrt((a(1)^2+a(2)^2)/4-a(3));
end

[Jafar Nur Arafah]

Hi Image Analyst, thanks for your response and understanding.

In your code, your creating your own randoms points, as for me, I already have the image with the points, so I change a little bit the code (by removing the coomand for plotting the random points and put the code to show the image and find the center, C (which I found that the coordoinate for C(x,y) = (217, 178)). Then I continue with the code you give me, but this error came out,

Can you check if my code is correct.

I = imread('test.GIF');
[ny,nx] = size(I) ;
imshow(I) ;
hold on
% Center 
C = round([nx ny]/2) ;
plot(C(1),C(2),'*r')   % mid point/ center of image 
hold on 
x = 217
y = 178
hold on
% First fit a circle
[xcFit, ycFit, RFit, a] = circfit(x,y)
% Now find the max and min fit
% Convert to a binary image
% Scale by 500 to get enough pixels to get 3 digits of precision.
scaleFactor = 500;
xScaled = scaleFactor*x;
yScaled = scaleFactor*y;
rows = ceil(scaleFactor * max(y));
columns = ceil(scaleFactor * max(x));
binaryImage = poly2mask(xScaled, yScaled, rows, columns);
subplot(2, 2, 3);
imshow(binaryImage);
title('binaryImage', 'fontSize', fontSize);
axis('on', 'image');
hold on;
plot(xScaled, yScaled, 'b+', 'MarkerSize', 15, 'LineWidth', 2);
% Get the Euclidean Distance Transform
edtImage = bwdist(~binaryImage);
% Display the image.
subplot(1, 2, 2);
imshow(edtImage, []);
title('Distance Transform', 'fontSize', fontSize);
axis('on', 'image');
% Find the max
minRadius = max(edtImage(:))
[rowMin, colMin] = find(edtImage == minRadius, 1, 'first');
hold on;
% Plot the center of the minimum interior circle.
plot(colMin, rowMin, 'r+', 'MarkerSize', 20);
% Plot the scaled x, y
plot(xScaled, yScaled, 'b+', 'MarkerSize', 15, 'LineWidth', 2);
% Plot the circles
viscircles([xcFit * scaleFactor, ycFit * scaleFactor; colMin, rowMin], [RFit * scaleFactor; minRadius]);
% Unscale to original dimensions
xCenterMin = colMin / scaleFactor
yCenterMin = rowMin / scaleFactor
minRadius = minRadius / scaleFactor
% Print out the average/fitted circle, and the min circle.
fprintf('Fitted circle : Radius = %f at (x, y) = (%f, %f).\n', RFit, xcFit, ycFit);
fprintf('Min circle    : Radius = %f at (x, y) = (%f, %f).\n', minRadius, xCenterMin, yCenterMin);
function [xc,yc,R,a] = circfit(x,y)
%CIRCFIT  Fits a circle in x,y plane
%
% [XC, YC, R, A] = CIRCFIT(X,Y)
% Result is center point (yc,xc) and radius R.  A is an optional
% output describing the circle's equation:
%
%   x^2+y^2+a(1)*x+a(2)*y+a(3)=0
         
% by Bucher izhak 25/oct/1991
n=length(x);  xx=x.*x; yy=y.*y; xy=x.*y;
A=[sum(x) sum(y) n;sum(xy) sum(yy) sum(y);sum(xx) sum(xy) sum(x)];
B=[-sum(xx+yy) ; -sum(xx.*y+yy.*y) ; -sum(xx.*x+xy.*y)];
a=A\B;
xc = -.5*a(1);
yc = -.5*a(2);
R  =  sqrt((a(1)^2+a(2)^2)/4-a(3));
end

I also have save a function for circfit.m as per shown by KSSV yesterday.

function [xc,yc,R,a] = circfit(x,y)
%CIRCFIT  Fits a circle in x,y plane
%
% [XC, YC, R, A] = CIRCFIT(X,Y)
% Result is center point (yc,xc) and radius R.  A is an optional
% output describing the circle's equation:
%
%   x^2+y^2+a(1)*x+a(2)*y+a(3)=0
         
% by Bucher izhak 25/oct/1991
n=length(x);  xx=x.*x; yy=y.*y; xy=x.*y;
A=[sum(x) sum(y) n;sum(xy) sum(yy) sum(y);sum(xx) sum(xy) sum(x)];
B=[-sum(xx+yy) ; -sum(xx.*y+yy.*y) ; -sum(xx.*x+xy.*y)];
a=A\B;
xc = -.5*a(1);
yc = -.5*a(2);
R  =  sqrt((a(1)^2+a(2)^2)/4-a(3));
end

I attached also the image I used 'test.GIF', for your reference.

Looking forward to your reply.

Thanks alot.

Nur Arafah Jafar

[Jafar Nur Arafah]

Hi Image Analyst,

As I run the code, figure below appear as the result.

Regards,

Nur Arafah

Answer 2

[KSSV]

YOu can try to fit a circle with the coordinates (x,y) you have. Check the below code:

function [xc,yc,R,a] = circfit(x,y)
%CIRCFIT  Fits a circle in x,y plane
%
% [XC, YC, R, A] = CIRCFIT(X,Y)
% Result is center point (yc,xc) and radius R.  A is an optional
% output describing the circle's equation:
%
%   x^2+y^2+a(1)*x+a(2)*y+a(3)=0
         
% by Bucher izhak 25/oct/1991

n=length(x);  xx=x.*x; yy=y.*y; xy=x.*y;
A=[sum(x) sum(y) n;sum(xy) sum(yy) sum(y);sum(xx) sum(xy) sum(x)];
B=[-sum(xx+yy) ; -sum(xx.*y+yy.*y) ; -sum(xx.*x+xy.*y)];
a=A\B;
xc = -.5*a(1);
yc = -.5*a(2);
R  =  sqrt((a(1)^2+a(2)^2)/4-a(3));

[Jafar Nur Arafah]

Hi KSSV,

Thanks for respond.

I tried to run the code, but error " Error: Function definitions are not permitted in this context." pop out. Then I tried to remove the "function [xc,yc,R,a] = circfit(x,y)", the circle did not come out (means that there is no circle produce overlap with the image). Do you have any idea on that?

Thanks.

Regards,

Nur Arafah

[KSSV]

Ou are running the code in a wrong way...you have to save the above code into a function and call the function.

[Image Analyst]

This fits a circle through the data, which could be part of the solution, but it does not find the max bounding circle (like my Answer).