Please , how can I find the center and the radius of the inscribed circle of a set of points

Please , how can I find the center and the radius of the inscribed circle of a set of points

[John D'Errico]

Do not add answers just to ask a question again!

Moved an answer into a comment:

"Are there any solution to get the inscribed circle inside an irregular polygone? thnak you"

[John D'Errico]

And my response to your question is that you never actually answered MY questions.

You need to define what it means for a circle to lie inside an irregular set of points. In some places, you state it as inside an irregular set of points. In others, it is now inside some general irregular polygon.

If you want an answer, you need to state CLEARLY what you need.

In any case, the computation you desire is NOT a trivial one.

Answer 1

[Image Analyst]

daeze, see if this is what you want. I solved it numerically by converting the data to a digital image, then using the Euclidean distance transform and finding the center of it. The max of the EDT is the center of the circle and its value there is the radius. If you need more than 3 digits of precision then you can increase the size of the image from 1000 to 10000 or larger.

% Initialization / clean-up code.
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 = 25;

x=[2.6381318;-18.66985584;-39.97784349;-69.8090262;-99.64020891;-92.131167993;-84.62315095;-57.4301004;-30.23704914;-3.65279788;22.93145337;34.09278024;45.2541071;23.94611945]
y=[75.28822303;63.72102973;52.15383644;43.02184173;33.88984702;19.48158871;5.07333039;5.07333039;5.07333039;5.07333039;5.07333039;25.77251893;46.4717064;60.87996471]
subplot(2, 2, 1);
plot(x, y, 'b.-', 'MarkerSize', 30);
grid on;
title('Original Points', 'FontSize', fontSize);
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]);

% Make data into a 1000x1000 image.
xMin = min(x)
xMax = max(x)
yMin = min(y)
yMax = max(y)
scalingFactor = 1000 / min([xMax-xMin, yMax-yMin])
x2 = (x - xMin) * scalingFactor + 1;
y2 = (y - yMin) * scalingFactor + 1;
mask = poly2mask(x2, y2, ceil(max(y2)), ceil(max(x2)));
% Display the image.
p2 = subplot(2, 2, 2);
imshow(mask);
axis(p2, 'on', 'xy');
title('Mask Image', 'FontSize', fontSize);

% Compute the Euclidean Distance Transform
edtImage = bwdist(~mask);
% Display the image.
p3 = subplot(2, 2, 3);
imshow(edtImage, []);
axis(p3, 'on', 'xy');
% Find the max
radius = max(edtImage(:))
% Find the center
[yCenter, xCenter] = find(edtImage == radius)
% Display circles over edt image.
viscircles(p3, [xCenter, yCenter], radius);
% Display polygon over image also.
hold on;
plot(x2, y2, 'r.-', 'MarkerSize', 30, 'LineWidth', 2);
title('Euclidean Distance Transform with Circle on it', 'FontSize', fontSize);

% Display the plot again.
subplot(2, 2, 4);
plot(x, y, 'b.-', 'MarkerSize', 30);
grid on;
% Show the circle on it.
hold on;
% Scale and shift the center back to the original coordinates.
xCenter = (xCenter - 1)/ scalingFactor + xMin
yCenter = (yCenter - 1)/ scalingFactor + yMin
radius = radius / scalingFactor
rectangle('Position',[xCenter-radius, yCenter-radius, 2*radius, 2*radius],'Curvature',[1,1]);
title('Original Points with Inscribed Circle', 'FontSize', fontSize);

[daeze zezezeze]

waw very nice work .dear Image Analyst , I thank you very much for your support.I will try to use this macro with others data points.It can give me the possiblity to get the inscribed circle for any shape ? thank you again:)

[daeze zezezeze]

Thank you very much 'Image Analyst" for your support.Your code works fine and this what I need to do.thanks again

Answer 2

[Jim Riggs]

I have worked this problem in the past, and the only solution I could find was a numerical searching method. There are a number of different ways to do it, but here is a suggested algorithm: It is assumed that the points are defined by a list of (x,y) coordinates.

1) Find the min and max of the x and the y coordinates. This defines a rectangle from (Xmin, Ymin) to (Xmax, Ymax).

2) Use the center of this rectangle as the initial point for the center of the circle. C = {(Xmax-Xmin)/2 , (Ymax-Ymin)/2 } You can compute the circle radius from C to any vertex of the box (e.g. (Xmin, Ymin).

3) Compute the distance from C to each data point.

4) Sort the distances to find the largest value. (Call this point P.) The circle radius may be reduced to this value, and it will still contain all of the data points, with one point, P, lying on the circle.

5) In order to further reduce the radius, you must move the center position. Move the center of the circle a small distance toward point P (call this new center point C2) and compute the distance from C2 to every other point. Continue adjusting the center point along a line from C to P until you have a second point that lies on the circle. (Call this second point Q)

6) Now define the line that goes through the center C2 and for which all points are equidistant to P and Q. Now move the circle center along this new line (reducing the radius value) until you have a third point which lies on the circle.

Since three points define a circle, you cannot reduce it any further. You have the center, and radius of a circle that inscribes all of the points and contains three of them.

Is it the smallest circle that can be found? I don't know. But this algorithm always converges to a solution.

(Watch out for the case where P and Q define the diameter of the circle. In this case, there will be no third point.)

[John D'Errico]

I have a funny feeling that this scheme (I'd call it a greedy algorithm) usually works, but that it can be broken by careful choice of the point set, designed to push the iterations into the wrong circle. I've used methods like it myself. In fact, this scheme probably always works on a nicely convex point set. But throw an interior point into the mix...

Answer 3

[Image Analyst]

Tell me if John's program does what you want: https://www.mathworks.com/matlabcentral/fileexchange/34767-a-suite-of-minimal-bounding-objects

[John D'Errico]

Not to be a downer here, but incircle will probably not be what is wanted, for at least two reasons.

Problem 1 is it works on the convex hull. So if the set is not convex, then there will be points that lie inside the convex hull. So incircle might be insufficient in that regard. As well, it will fail for another reason.

Problem 2: Consider a set of 4 points, on the corners of a unit square, thus [0 1]X[0 1]. incircle will find a circle that touches each edge of the square. So the center will be [0.5 0.5], and the radius will also be 0.5.

But the incircle that has been requested is the one that lies inside a point set. So in this case, the desired circle will be centered at [0.5 0.5] again, but the desired radius will be sqrt(2)/2=0.7071...

So I predict incircle will be inadequate for the task. Good try though.

[Image Analyst]

I wasn't sure, and am still not, whether the circle was one that was outside and contained a roughly circular cluster of points, or whether there was a "hole" in the points (like a dotted donut) and he wanted the circle to be inside the points but containing the edge of the hole (essentially outlining the hole of the donut). A diagram of what he wants would certainly clarify this.

[John D'Errico]

It has been a while since I wrote incircle, so I had to check myself to be sure.

[C,R] = incircle([0 0 1 1],[0 1 0 1])
C =
          0.5          0.5
R =
          0.5

So incircle is definitely looking at the convex hull polygon.

Answer 4

[John D'Errico]

The largest incircle of a scattered set may be difficult to find. A simple heuristic, like that suggested by Jim will work if the set is a convex one. I'd predict it can also do strange things if you choose a nasty set designed to break it. In fact, you can probably make life interesting for any algorithm with a set of points like this:

x = [0.1 -0.1 0 0 -1 -1 1 1];
y = [0 0 .1 -.1 -1 1 -1 1];
plot(x,y,'o')
axis equal
grid on
axis square

Be careful, since the largest circle that passes through three of those points and does not contain any other point has its center outside of the convex hull.

So careful definition of the desired inscribed circle may be needed to resolve the question. Must the center lie inside the convex hull?

[daeze zezezeze]

I used to solidworks to check if the resukts of your code are ok or not and I found that the results were not bad but there are some error to be adjusted for the radius and the centers.For my case study I found that the radius error =0,26mm and the center also there are some errors !

[daeze zezezeze]

In my case the set of points are as follow :

x=[2.6381318;-18.66985584;-39.97784349;-69.8090262;-99.64020891;-92.131167993;-84.62315095;-57.4301004;-30.23704914;-3.65279788;22.93145337;34.09278024;45.2541071;23.94611945]

Y=[75.28822303;63.72102973;52.15383644;43.02184173;33.88984702;19.48158871;5.07333039;5.07333039;5.07333039;5.07333039;5.07333039;25.77251893;46.4717064;60.87996471]
In attachement u find the resukts that I got.it's not very good I think

function [C,R] = incircle(x,y)
% incircle: compute the maximal in-circle of the polygonal convex hull of a set of points in the plane
%
% [C,R] = incircle(x,y)
%    Called with a pair of vectors of the same length, 
%    incircle computes the convex hull in 2-dimensions,
%    then finds the maximal radius in-circle that lies
%    entirely within the polygon of the convex hull.
%
% [C,R] = incircle(xy)
%    Called with a single nx2 array, incircle treats
%    this as a list of points in 2 dimensions. Each row
%    of xy is treated as a single point. The convex hull
%    is computed and the maximal in-circle then computed
%    for the convex hull.
%
% [C,R] = incircle(xy,edges)
%    Called with two arguments, the first of which is an
%    nx2 array of data points, and the second a list of
%    edges of the convex hull as would be generated by
%    convhull (or convhulln), the maximal in-circle is
%    generated for the hull as provided. No test is made
%    to assure the hull is truly convex.
%
% All input forms return a row vector of length 2 (C)
% which defines the center of the in-circle, and the
% radius (R) of the in-circle.
%
% Example:
% 
% [C,R] = incircle(rand(1000,1),rand(1000,1))
%
% C =
%      0.50191      0.49999
%
% R =
%      0.49788
%
% Example:
%
% x = randn(100,1) + 1;
% y = randn(100,1) - 2.5;
% H = convhull(x,y);
% [C,R] = incircle([x,y],H)
%
% C =
%      1.0212      -2.3458
%
% R =
%      2.1329
%
% t = linspace(0,2*pi,200);
% xc = R*cos(t) + C(1);
% yc = R*sin(t) + C(2);
% plot(x(H),y(H),'-r',x,y,'ko',xc,yc,'-b')
% axis equal
% xlabel X
% ylabel Y
% title 'Convex hull, with polygon in-circle'
% grid on
%
%
% See also: minboundcicle, convhull
%
% Author: John D'Errico
% e-mail: [e-mail]
% Release: 1.0
% Release date: 6/14/09

if (nargin < 1) || (nargin > 2)
  error('INCIRCLE:improperarguments', ...
    'incircle requires exactly 1 or 2 arguments')
elseif (nargin == 2) && isvector(x) && isvector(y) && (numel(x) == numel(y))
  % a pair of vectors
  xy = [x(:),y(:)];
  % compute the hull. I prefer convhulln.
  edges = convhulln(xy);
elseif (nargin == 1) && (size(x,2) == 2)
  % a single list of points as rows of x
  xy = x;
  edges = convhulln(xy);
elseif (nargin == 2) && (size(x,2) == 2) && (size(y,2) == 2)
  % y must be a list of edges from convhulln
  xy = x;
  edges = y;
elseif (nargin == 2) && (size(x,2) == 2) && isvector(y) && (y(1) == y(end))
  % y must be a list of edges from convhull
  xy = x;
  I = y(:);
  % in case the first point was not wrapped in the polygon
  if I(1) ~= I(end)
    I = [I;I(1)];
  end
  edges = [I(1:(end-1)),I(2:end)];
else
  % none of the forms I allow seem to fit.
  % give up and throw an error.
  error('INCIRCLE:invaliddata', ...
    'x and y do not seem to fit into any of the allowed forms for incircle input')
end
ne = size(edges,1);

% the end points of each edge are...
A = xy(edges(:,1),:);
B = xy(edges(:,2),:);

% the normal vector to each edge
N = (B - A)*[0 1;-1 0];

% normalize to unit length
L = sqrt(sum(N.^2,2));
N = N./[L,L];

% a central point inside the hull itself
C0 = mean(A,1);

% ensure the normals point inwards. While
% I can use the signed area for a hull as
% generated by convhull (suggestion by Bruno)
% this test is also efficient, and it deals
% with the case of a hull provided from some
% unknown and less trusted source.
k = sum(N.*bsxfun(@minus,C0,A),2) < 0;
N(k,:) = -N(k,:);

% formulate the linear programming problem.
% given a point C inside the hull, the distance
% from C to the edge that contains A(i,:) is
% given by (dot(N(i,:),C - A(i,:))). If this
% number is negative for any edge of the hull,
% then C is outside of the hull.
%
% Now think of it as a set of slack variables,
% one for each edge of the hull. Given the
% unknown point C that defines the center of
% the in-circle,
%
%  dot(N(i,:),C-A(i,:)) - S(i) == 0
%
% Thus the vector S is of length ne, where ne is
% the number of edges in the convex hull. Every
% element of S must be positive.
%
% Create one more slack variable, a scalar R.
% R is the minimum value that S attains for a
% given point C. 
%
%   R >= 0
%   S(i) >= R
%
% The objective for our linear programming problem
% is simple. It is just -R. When we find the
% solution to the LP problem that satisfies the
% equality constraints above between C and S, the
% bound constraint on T, and the inequality
% constraints between S and R, the solution yields
% the maximal radius in-circle that fits inside
% the convex hull polygon.
%
% The unknowns for the LP are a list of ne+3
% variables. The first two variables represent
% the center coordinates of the circle. X(3) = R
% is the circle radius. The remainder of the
% variables are the slack variabls S(i).

% equality constraints defined by the dot products
% with the normal vectors.
Aeq = [N,zeros(ne,1),-eye(ne)];
beq = sum(N.*A,2);

% lower bounds only for the slack variables
LB = [-inf; -inf; 0; zeros(ne,1)];
% there are no upper bounds
UB = [];

% inequalities defined by the slack variable
% constraints
A = [zeros(ne,2),ones(ne,1),-eye(ne)];
b = zeros(ne,1);

% the objective is just -R
f = zeros(ne+3,1);
f(3) = -1;

% just turn off the output message
options = optimset('linprog');
options.Display = 'off';

% linprog does all of the hard work.
result = linprog(f,A,b,Aeq,beq,LB,UB,[],options);

% unpack the circle parameters
C = result(1:2)';
R = result(3)-0.26;

[John D'Errico]

I explicitly stated that incircle does not solve the problem you want to solve, and I have said that several times. Trying to use the code, and claiming there is a bug in the code is insanity, when you apparently have no clue as to what the code actually does.

E = convhull(x,Y)
E =
     1
     5
     7
     8
     9
    10
    11
    13
     1

[C,R] = incircle(x,Y)
C =
      -1.2818       38.087
R =
       33.013

t = linspace(0,2*pi,100);
plot(x,Y,'bo',x(E),Y(E),'r-',C(1) + R*cos(t),C(2) + R*sin(t),'g-')
axis equal

INCIRCLE works PERFECTLY WELL. It does EXACTLY what it is designed to do. You can see that the circle fits neatly into the red polygon created by the convex hull.

That incircle does not solve your problem does not mean there is a bug in the code, as you have stated.

INCIRCLE is designed to fit a circle of maximum radius inside the convex hull of your data. As you see, it did EXACTLY that.

It is not able to solve a problem for which it was not designed. Computers do not read your mind, changing their own programming to coincide with what you want them to do.

[daeze zezezeze]

Hi dear D'Erico , I understanf now what you want to say.It's ok nd sorrry because I have not experience with matlab ! just I use visual basic to solve my problems.Thank you any way

[Jan]

@John: "claiming there is a bug in the code is insanity" - no, it is not insane. It is a misunderstanding. Daeze tries to solve a problem and did not understand exactly, what the function incircle does.

I'd suggest to remove the code on incircle from this thread. The link to the FEX submission is better.

@daeze: Your questions are welcome in the forum. It is useful that you show your own effort, although the approach with incircle is not working for your problem. You try to solve your problem and the forum tries to help you. I appreciate your polite reaction.

Answer 5

[daeze zezezeze]

Sincerly I tried to apply this Matlab code :

function [C,R] = incircle(x,y)
% incircle: compute the maximal in-circle of the polygonal convex hull of a set of points in the plane
%
% [C,R] = incircle(x,y)
%    Called with a pair of vectors of the same length, 
%    incircle computes the convex hull in 2-dimensions,
%    then finds the maximal radius in-circle that lies
%    entirely within the polygon of the convex hull.
%
% [C,R] = incircle(xy)
%    Called with a single nx2 array, incircle treats
%    this as a list of points in 2 dimensions. Each row
%    of xy is treated as a single point. The convex hull
%    is computed and the maximal in-circle then computed
%    for the convex hull.
%
% [C,R] = incircle(xy,edges)
%    Called with two arguments, the first of which is an
%    nx2 array of data points, and the second a list of
%    edges of the convex hull as would be generated by
%    convhull (or convhulln), the maximal in-circle is
%    generated for the hull as provided. No test is made
%    to assure the hull is truly convex.
%
% All input forms return a row vector of length 2 (C)
% which defines the center of the in-circle, and the
% radius (R) of the in-circle.
%
% Example:
% 
% [C,R] = incircle(rand(1000,1),rand(1000,1))
%
% C =
%      0.50191      0.49999
%
% R =
%      0.49788
%
% Example:
%
% x = randn(100,1) + 1;
% y = randn(100,1) - 2.5;
% H = convhull(x,y);
% [C,R] = incircle([x,y],H)
%
% C =
%      1.0212      -2.3458
%
% R =
%      2.1329
%
% t = linspace(0,2*pi,200);
% xc = R*cos(t) + C(1);
% yc = R*sin(t) + C(2);
% plot(x(H),y(H),'-r',x,y,'ko',xc,yc,'-b')
% axis equal
% xlabel X
% ylabel Y
% title 'Convex hull, with polygon in-circle'
% grid on
%
%
% See also: minboundcicle, convhull
%
% Author: John D'Errico
% e-mail: [e-mail]
% Release: 1.0
% Release date: 6/14/09

if (nargin < 1) || (nargin > 2)
  error('INCIRCLE:improperarguments', ...
    'incircle requires exactly 1 or 2 arguments')
elseif (nargin == 2) && isvector(x) && isvector(y) && (numel(x) == numel(y))
  % a pair of vectors
  xy = [x(:),y(:)];
  % compute the hull. I prefer convhulln.
  edges = convhulln(xy);
elseif (nargin == 1) && (size(x,2) == 2)
  % a single list of points as rows of x
  xy = x;
  edges = convhulln(xy);
elseif (nargin == 2) && (size(x,2) == 2) && (size(y,2) == 2)
  % y must be a list of edges from convhulln
  xy = x;
  edges = y;
elseif (nargin == 2) && (size(x,2) == 2) && isvector(y) && (y(1) == y(end))
  % y must be a list of edges from convhull
  xy = x;
  I = y(:);
  % in case the first point was not wrapped in the polygon
  if I(1) ~= I(end)
    I = [I;I(1)];
  end
  edges = [I(1:(end-1)),I(2:end)];
else
  % none of the forms I allow seem to fit.
  % give up and throw an error.
  error('INCIRCLE:invaliddata', ...
    'x and y do not seem to fit into any of the allowed forms for incircle input')
end
ne = size(edges,1);

% the end points of each edge are...
A = xy(edges(:,1),:);
B = xy(edges(:,2),:);

% the normal vector to each edge
N = (B - A)*[0 1;-1 0];

% normalize to unit length
L = sqrt(sum(N.^2,2));
N = N./[L,L];

% a central point inside the hull itself
C0 = mean(A,1);

% ensure the normals point inwards. While
% I can use the signed area for a hull as
% generated by convhull (suggestion by Bruno)
% this test is also efficient, and it deals
% with the case of a hull provided from some
% unknown and less trusted source.
k = sum(N.*bsxfun(@minus,C0,A),2) < 0;
N(k,:) = -N(k,:);

% formulate the linear programming problem.
% given a point C inside the hull, the distance
% from C to the edge that contains A(i,:) is
% given by (dot(N(i,:),C - A(i,:))). If this
% number is negative for any edge of the hull,
% then C is outside of the hull.
%
% Now think of it as a set of slack variables,
% one for each edge of the hull. Given the
% unknown point C that defines the center of
% the in-circle,
%
%  dot(N(i,:),C-A(i,:)) - S(i) == 0
%
% Thus the vector S is of length ne, where ne is
% the number of edges in the convex hull. Every
% element of S must be positive.
%
% Create one more slack variable, a scalar R.
% R is the minimum value that S attains for a
% given point C. 
%
%   R >= 0
%   S(i) >= R
%
% The objective for our linear programming problem
% is simple. It is just -R. When we find the
% solution to the LP problem that satisfies the
% equality constraints above between C and S, the
% bound constraint on T, and the inequality
% constraints between S and R, the solution yields
% the maximal radius in-circle that fits inside
% the convex hull polygon.
%
% The unknowns for the LP are a list of ne+3
% variables. The first two variables represent
% the center coordinates of the circle. X(3) = R
% is the circle radius. The remainder of the
% variables are the slack variabls S(i).

% equality constraints defined by the dot products
% with the normal vectors.
Aeq = [N,zeros(ne,1),-eye(ne)];
beq = sum(N.*A,2);

% lower bounds only for the slack variables
LB = [-inf; -inf; 0; zeros(ne,1)];
% there are no upper bounds
UB = [];

% inequalities defined by the slack variable
% constraints
A = [zeros(ne,2),ones(ne,1),-eye(ne)];
b = zeros(ne,1);

% the objective is just -R
f = zeros(ne+3,1);
f(3) = -1;

% just turn off the output message
options = optimset('linprog');
options.Display = 'off';

% linprog does all of the hard work.
result = linprog(f,A,b,Aeq,beq,LB,UB,[],options);

% unpack the circle parameters
C = result(1:2)';
R = result(3);

So after applying this code and defined a set of points I found this result like in the picture.I think there are something wrong in this code because the circle is outside the contour.

[John D'Errico]

NO. It is not the code. ARGH. Listen to what I said. incircle does not solve the problem you want to solve.

Answer 6

[Jan]

Your example contains a small set of points only. Then a brute force approach can solve the problem in a bearable time: Create all possible sets of 3 points, calculate the corresponding circles and choose the one, which matches all criteria.

But what exactly are the criteria? John has asked already: Must the center be inside the convex hull? Which circle is wanted for John's wonderful example above? As long as it is not defined mathematically, what an "inscribed circle of a set of points" is, posting solutions is based on guessing. So please define your problem exactly.

This is a global optimization problem: Starting at one specific point and trying to let a circle grow as far as possible can stuck in a local maximum. For a reliable search, you must start from each triangle of the triangulated area. And for this triangulation it matters, if you mean an alpha shape (which could contain wholes) or the convex hull.

Answer 7

[John D'Errico]

If I were going to solve the problem of finding an incircle for a general non-degenerate closed polygon, I MIGHT try the following algorithm:...

1. Verify the polygon is closed. If not, then form the closure.

2. Verify the polygon does not cross itself, or is non-degenerate in some other way.

3. Verify the polygon is not already convex. If is is convex, then incircle does the job efficiently, so you are done.

4. Triangulate the polygon. Ear clipping is entirely adequate here. Fast, efficient and stable as long as the polygon is not degenerate.

5. Generate a large set of random points inside the polygon. 10000 points should suffice, and is efficiently done given a triangulation.

6. Find the internal point (among all the random set) which has a maximum minimal distance to any edge of the polygon.

7. Use fminsearch to optimize the center of a circle which maximizes the minimum distance to any edge of the polygon. fminsearch will be entirely adequate here, since it is only a 2-dimensional problem. The start point will be the best point chosen in step 6.