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

22 次查看（过去 30 天）
daeze zezezeze
daeze zezezeze 2018-1-18
编辑： Image Analyst 2020-3-15
Please , how can I find the center and the radius of the inscribed circle of a set of points
  2 个评论
John D'Errico
John D'Errico 2018-2-1
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
John D'Errico 2018-2-1
编辑：John D'Errico 2018-2-1
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Image Analyst
Image Analyst 2018-2-1
编辑：Image Analyst 2020-3-15
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);
  2 个评论
daeze zezezeze
daeze zezezeze 2018-2-1
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
daeze zezezeze 2018-2-1
Thank you very much 'Image Analyst" for your support.Your code works fine and this what I need to do.thanks again

请先登录，再进行评论。

更多回答（6 个）

Jim Riggs
Jim Riggs 2018-1-18
编辑：Jim Riggs 2018-1-18
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.)
  1 个评论
John D'Errico
John D'Errico 2018-1-19
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...

请先登录，再进行评论。

Image Analyst
Image Analyst 2018-1-18
  3 个评论
显示 1更早的评论
Image Analyst
Image Analyst 2018-1-18
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
John D'Errico 2018-1-19
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.

请先登录，再进行评论。

John D'Errico
John D'Errico 2018-1-19
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?
  5 个评论
显示 3更早的评论
daeze zezezeze
daeze zezezeze 2018-2-1
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
Jan 2018-2-1
编辑：Jan 2018-2-1
@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.

请先登录，再进行评论。

daeze zezezeze
daeze zezezeze 2018-1-19
编辑：Torsten 2018-1-19
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.
Jan
Jan 2018-2-1
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.
John D'Errico
John D'Errico 2018-2-1
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.

类别

在 Help Center 和 File Exchange 中查找有关 Bounding Regions 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by