Asked by daeze zezezeze
on 18 Jan 2018

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

Answer by Image Analyst
on 1 Feb 2018

Accepted Answer

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
on 1 Feb 2018

daeze zezezeze
on 1 Feb 2018

Sign in to comment.

Answer by Jim Riggs
on 18 Jan 2018

Edited by Jim Riggs
on 18 Jan 2018

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
on 19 Jan 2018

Sign in to comment.

Answer by Image Analyst
on 18 Jan 2018

John D'Errico
on 18 Jan 2018

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
on 18 Jan 2018

John D'Errico
on 19 Jan 2018

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.

Sign in to comment.

Answer by John D'Errico
on 19 Jan 2018

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?

John D'Errico
on 1 Feb 2018

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
on 1 Feb 2018

Jan
on 1 Feb 2018

@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.

Sign in to comment.

Answer by daeze zezezeze
on 19 Jan 2018

Edited by Torsten
on 19 Jan 2018

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: woodchips@rochester.rr.com

% 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
on 1 Feb 2018

Sign in to comment.

Answer by Jan
on 1 Feb 2018

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.

Sign in to comment.

Answer by John D'Errico
on 1 Feb 2018

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.

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 2 Comments

## John D'Errico (view profile)

Direct link to this comment:https://www.mathworks.com/matlabcentral/answers/377838-please-how-can-i-find-the-center-and-the-radius-of-the-inscribed-circle-of-a-set-of-points#comment_530681

## John D'Errico (view profile)

Direct link to this comment:https://www.mathworks.com/matlabcentral/answers/377838-please-how-can-i-find-the-center-and-the-radius-of-the-inscribed-circle-of-a-set-of-points#comment_530682

Sign in to comment.