Thread Subject: How can I find the largest inner circle of a arbitrary polygon?

This polygon is arbitrary. And the inner circle need not be tangent
with each edge of the polygon. So how can I get the center and radius
of this inner circle?

Thanks very much!

On Aug 24, 7:23=A0am, Eric <Yonli.Z...@gmail.com> wrote:
> This polygon is arbitrary. And the inner circle need not be tangent
> with each edge of the polygon. So how can I get the center and radius
> of this inner circle?
>
> Thanks very much!

-----------------------------------------------------------
Eric:
I think this might do it. Just average all the x and y coordinates
together to get the center. Then go through each point calculating
the distance of each polygon vertex from the center. Then take the 3
vertices with the shortest distances. Then fit a circle to those 3
points - I believe the way to do that has been discussed many times
before here and I think there's MATLAB file exchange code posted for
doing that.
Regards,
ImageAnalyst

ImageAnalyst <imageanalyst@mailinator.com> wrote in message
<a19898c7-2bae-42ec-a80d-
3390f2a448d7@e39g2000hsf.googlegroups.com>...
> On Aug 24, 7:23=A0am, Eric <Yonli.Z...@gmail.com> wrote:
> > This polygon is arbitrary. And the inner circle need not be tangent
> > with each edge of the polygon. So how can I get the center and radius
> > of this inner circle?

If it is not tangent to the edges, then how
will you define the circle?



> >
> > Thanks very much!
>
> ----------------------------------------------------------
-
> Eric:
> I think this might do it. Just average all the x and y coordinates
> together to get the center. Then go through each point calculating
> the distance of each polygon vertex from the center. Then take the 3
> vertices with the shortest distances. Then fit a circle to those 3
> points - I believe the way to do that has been discussed many times
> before here and I think there's MATLAB file exchange code posted for
> doing that.
> Regards,
> ImageAnalyst

I'll claim I can give a simple counterexample
where this will fail to compute the correct
incircle.

Even worse, I can give an example of a
non-convex polygon which will leave the
resulting circle with its center outside of
the polygon.

John

"John D'Errico" <woodchips@rochester.rr.com> wrote in message
<g8rn7t$o1n$1@fred.mathworks.com>...
> ImageAnalyst <imageanalyst@mailinator.com> wrote in message
> <a19898c7-2bae-42ec-a80d-

> > Eric:
> > I think this might do it. Just average all the x and y coordinates
> > together to get the center. Then go through each point calculating
> > the distance of each polygon vertex from the center. Then take the 3
> > vertices with the shortest distances. Then fit a circle to those 3
> > points - I believe the way to do that has been discussed many times
> > before here and I think there's MATLAB file exchange code posted for
> > doing that.
> > Regards,
> > ImageAnalyst
>
> I'll claim I can give a simple counterexample
> where this will fail to compute the correct
> incircle.

Consider the polygon defined by (x,y) pairs:

x = [0, 1, 1+eps, 1+eps, 1];
y = [0, -.1, -.05, .05 .1];

The centroid occurs at (0.8,0), so the closest
vertices are on the right side of the polygon.
But the circle passing though those points
has a nearly infinite radius, so it will clearly
contain the point (0,0) inside of it.

 
> Even worse, I can give an example of a
> non-convex polygon which will leave the
> resulting circle with its center outside of
> the polygon.

As easy to do.

x = [0 1 eps 0];
y = [0 0 eps 1];

Non-convex polygons will make the problem
quite difficult.

John.

On Aug 24, 9:21=A0pm, "John D'Errico" <woodch...@rochester.rr.com>
wrote:
> ImageAnalyst <imageanal...@mailinator.com> wrote in message
>
> <a19898c7-2bae-42ec-a80d-
> 3390f2a44...@e39g2000hsf.googlegroups.com>...
>
> > On Aug 24, 7:23=3DA0am, Eric <Yonli.Z...@gmail.com> wrote:
> > > This polygon is arbitrary. And the inner circle need not be tangent
> > > with each edge of the polygon. So how can I get the center and radius
> > > of this inner circle?
>
> If it is not tangent to the edges, then how
> will you define the circle?
>
>
>
>
>
> > > Thanks very much!
>
> > ----------------------------------------------------------
> -
> > Eric:
> > I think this might do it. =A0Just average all the x and y coordinates
> > together to get the center. =A0Then go through each point calculating
> > the distance of each polygon vertex from the center. =A0Then take the 3
> > vertices with the shortest distances. =A0Then fit a circle to those 3
> > points - I believe the way to do that has been discussed many times
> > before here and I think there's MATLAB file exchange code posted for
> > doing that.
> > Regards,
> > ImageAnalyst
>
> I'll claim I can give a simple counterexample
> where this will fail to compute the correct
> incircle.
>
> Even worse, I can give an example of a
> non-convex polygon which will leave the
> resulting circle with its center outside of
> the polygon.
>
> John

The circle need not be tangent with all edges, but sure be tangent
with some edges. and the radius must be largest.

Eric <Yonli.Zhan@gmail.com> wrote in message
<8f147863-db3b-4b0c-b073-3f4c785d6196@i24g2000prf.googlegroups.com>...

>
> The circle need not be tangent with all edges, but sure be
tangent
> with some edges.

Really? If the polygone is a symmtric star with five sides
(such as morocco or vietnam flag), does the largest inner
circle tangent to any edge?

Bruno

Eric <Yonli.Zhan@gmail.com> wrote in message <8f147863-db3b-4b0c-
b073-3f4c785d6196@i24g2000prf.googlegroups.com>...
> On Aug 24, 9:21=A0pm, "John D'Errico" <woodch...@rochester.rr.com>
> wrote:
> > ImageAnalyst <imageanal...@mailinator.com> wrote in message
> >
> > <a19898c7-2bae-42ec-a80d-
> > 3390f2a44...@e39g2000hsf.googlegroups.com>...
> >
> > > On Aug 24, 7:23=3DA0am, Eric <Yonli.Z...@gmail.com> wrote:
> > > > This polygon is arbitrary. And the inner circle need not be tangent
> > > > with each edge of the polygon. So how can I get the center and
radius
> > > > of this inner circle?
> >
> > If it is not tangent to the edges, then how
> > will you define the circle?
> >
> >
> >
> >
> >
> > > > Thanks very much!
> >
> > > ------------------------------------------------------
----
> > -
> > > Eric:
> > > I think this might do it. =A0Just average all the x and y coordinates
> > > together to get the center. =A0Then go through each point calculating
> > > the distance of each polygon vertex from the center. =A0Then take
the 3
> > > vertices with the shortest distances. =A0Then fit a circle to those 3
> > > points - I believe the way to do that has been discussed many times
> > > before here and I think there's MATLAB file exchange code posted for
> > > doing that.
> > > Regards,
> > > ImageAnalyst
> >
> > I'll claim I can give a simple counterexample
> > where this will fail to compute the correct
> > incircle.
> >
> > Even worse, I can give an example of a
> > non-convex polygon which will leave the
> > resulting circle with its center outside of
> > the polygon.
> >
> > John
>
> The circle need not be tangent with all edges, but sure be tangent
> with some edges. and the radius must be largest.

As Bruno points out, if the polygon is
non-convex, then the maximal in-circle
need not be tangent to any edge.

Are you willing to restrict this to a convex
polygon?

John

On Aug 24, 9:37=A0am, "John D'Errico" <woodch...@rochester.rr.com>
wrote:
> "John D'Errico" <woodch...@rochester.rr.com> wrote in message
>
> <g8rn7t$o1...@fred.mathworks.com>...
>
>
>
>
>
> > ImageAnalyst <imageanal...@mailinator.com> wrote in message
> > <a19898c7-2bae-42ec-a80d-
> > > Eric:
> > > I think this might do it. =A0Just average all the x and y coordinates
> > > together to get the center. =A0Then go through each point calculating
> > > the distance of each polygon vertex from the center. =A0Then take the=
 3
> > > vertices with the shortest distances. =A0Then fit a circle to those 3
> > > points - I believe the way to do that has been discussed many times
> > > before here and I think there's MATLAB file exchange code posted for
> > > doing that.
> > > Regards,
> > > ImageAnalyst
>
> > I'll claim I can give a simple counterexample
> > where this will fail to compute the correct
> > incircle.
>
> Consider the polygon defined by (x,y) pairs:
>
> x =3D [0, 1, 1+eps, 1+eps, 1];
> y =3D [0, -.1, -.05, .05 .1];
>
> The centroid occurs at (0.8,0), so the closest
> vertices are on the right side of the polygon.
> But the circle passing though those points
> has a nearly infinite radius, so it will clearly
> contain the point (0,0) inside of it.
>
> > Even worse, I can give an example of a
> > non-convex polygon which will leave the
> > resulting circle with its center outside of
> > the polygon.
>
> As easy to do.
>
> x =3D [0 1 eps 0];
> y =3D [0 0 eps 1];
>
> Non-convex polygons will make the problem
> quite difficult.
>
> John.- Hide quoted text -
>
> - Show quoted text -

---------------------------------------------------------------------------=
--------------
John:
Yes, you're right. Another counterexample would be a blob shaped like
the letter "C." I thought of another way that would work with C-
shaped blobs (involving Euclidean DIstance Maps, which MATLAB calls
bwdist). It is to take the EDM, find the location of the pixel having
the greatest EDM value, and place a circle centered there with a
radius of the EDM value. I haven't thought it out completely or test
it but do you think it would work in general? Interesting problem -
interested in your solution(s).
Regards,
ImageAnalyst
---------------------------------------------------------------------------=
----------

ImageAnalyst <imageanalyst@mailinator.com> wrote in message <ad67cdfc-
cbde-451b-b154-bd5f009aa7e1@k7g2000hsd.googlegroups.com>...

> John:
> Yes, you're right. Another counterexample would be a blob shaped like
> the letter "C." I thought of another way that would work with C-
> shaped blobs (involving Euclidean DIstance Maps, which MATLAB calls
> bwdist). It is to take the EDM, find the location of the pixel having
> the greatest EDM value, and place a circle centered there with a
> radius of the EDM value. I haven't thought it out completely or test
> it but do you think it would work in general? Interesting problem -
> interested in your solution(s).
> Regards,
> ImageAnalyst

I don't know of a solution for a non-convex
polygon. I'd bet this is a non-trivial problem
for large, non-convex polygons. I'd probably
employ a Monte Carlo solution, much like the
idea that you described.

For a general polygon, use a random
sampling, generating lots of random points
inside the polygon. Just oversample the region,
then a test using inpolygon to pick out the
points inside the polygon. Now find the nearest
neighbor for each point inside the polygon to a
point on the boundary of the polygon.

The point with the largest distance to its
identified nearest neighbor is a good
approximation to the center of the incircle.
To make this work well, the polygon must
not have any long edges in it. I don't know
the IPT at all well, but I'll bet this is the same
as what you described.

John

Hi,

cause i liked this puzzle, I took some of the suggested
ideas together and wrote this snippet. Maybe I should
better put it on the FEX, but than I should do more
commenting and error checking...

For the parallel sections problem cited in the snippet, I
think one could deal with this by cutting those segments to
a length of their distance... I have to rethink that.

Bye Gwen

<SNIP>
function oldCircle = c_scLargestInnerCircle
(x,y,maxIter,fPlot)
% c_scLargestInnerCircle finds the largest circle in an
arbitrary polygon
%
% Input:
% x,y column vectors of the vertices of the polygon to
investigate
% maxIter: maximum number of iterations to investigate
% fPlot: boolean value indicating wether to plot
intermediate results
% or not
%
% Output:
% oldCircle returns the largest inner circle with
oldCircle(1:2) being
% the center of the circle and oldCircle(3) its
radius
%
% CAUTION: the function is far from failsave, especially it
fails if
% no voronoi point is inside the polygon (fix? maybe use
the midpoint
% between two outside points as voronoi point... or better
the
% intersection with the delaunay triangles) e.g.
% x = [0 1 0.1 0 0];
% y = [0 0 0.1 1 0];
% c_scLargestInnerCircle(x,y)
%
% if the polygon contains large parallel structures, the
iteration will
% take very long...
% x = [0 5 5 0.1 0.1 5 5 0 0]';
% y = [0 0 1 1 1.1 1.1 2 2 0]';
% c_scLargestInnerCircle(x,y,50)

% Examples:
% x = [0 -1 -1.5 -2 -1.5 -2 -0.5 0 1 1.5*sqrt
(0.5) 1.5 2 1 0]';
% y = [1.5 2 1 0.5 0 -1.5 -1.5 -2 -1.5 -1.5*sqrt
(0.5) -1 0.5 1 1.5]';
%
% x = [0 -1 -1.5 -1.5 -0.5 0 1 1.5*sqrt(0.5) 1.5
2 1 0]';
% y = [1.5 2 1 -1.5 -1.5 -2 -1.5 -1.5*sqrt(0.5) -1
0.5 1 1.5]';

% function inspired by a discussion on the cssm
%
http://www.mathworks.com/matlabcentral/newsreader/view_threa
d/235050#597170

% circle is a function on the FEX by Zhenhai Wang
%
http://www.mathworks.com/matlabcentral/fileexchange/loadFile
.do?objectId=
% 2876&objectType=file
 
 oldCircle = [];

 if nargin < 3
   maxIter = 10;
 end
 
 if nargin<4
   fPlot = true;
 end
 
 if isempty(x) || isempty(y)
   return
 end
 
 % is polygon closed?
 if (x(1)~=x(end)) || (y(1)~=y(end))
   x = [x; x(1)];
   y = [y; y(1)];
 end

jj=0;
while jj < maxIter
  
  TRI = delaunay(x,y);
  [vx, vy] = voronoi(x,y,TRI);
    
  if fPlot
    figure
    hold all
    plot(0,0,'*g')
    plot(x,y,'*-r')
    if ~isempty(oldCircle)
      circle(oldCircle(1:2),oldCircle(3),'b');
    end
    triplot(TRI,x,y,':c')
    plot(vx,vy,'g-')
  end
  
  ttt = [vx(:) vy(:)];
  ttt = unique(ttt,'rows');

  isinpoly = inpolygon(ttt(:,1),ttt(:,2),x,y);
  ttt(~isinpoly,:)=[];

  len1 = size(ttt,1);
  len2 = size(x,1);
  mindist = 0;
  for ii=1:len1
    tmp = repmat(ttt(ii,:),len2,1);
    dist = min((x(:)-tmp(:,1)).^2+(y(:)-tmp(:,2)).^2);
    if dist > mindist
      val = ii;
      mindist=dist;
    end
  end

  radCirc = sqrt(mindist);

  oldCircle = [ttt(val,1) ttt(val,2) radCirc];
  if fPlot
    circle(oldCircle(1:2),oldCircle(3),3600,'k');
    plot(oldCircle(1),oldCircle(2),'*k');
    axis equal
  end

  inPoint = [];
  for ii=1:len2-1
    [d,P] = distance([x(ii) y(ii)],[x(ii+1) y(ii+1)],...
                     [ttt(val,1),ttt(val,2)]);
    if d<radCirc
      if fPlot
        plot(P(1),P(2),'*k');
      end
      inPoint = [inPoint; ii, P(1:2)'];
    end
  end

  if ~isempty(inPoint)
    nInPoint = size(inPoint,1);
    newPoints = nan(len2+nInPoint,2);
    idx1 = [1; inPoint(:,1)+1];
    idx2 = [inPoint(:,1); len2];
    for ii=1:length(idx1)
      newPoints((idx1(ii):idx2(ii))+(ii-1),:) = [x(idx1
(ii):idx2(ii)) y(idx1(ii):idx2(ii))];
    end
    newPoints(inPoint(:,1)+(1:nInPoint)',:) = inPoint
(:,2:3);

    x = newPoints(:,1);
    y = newPoints(:,2);

    jj = jj + 1;
  else
    jj=maxIter;
  end

end % while jj < maxIter


function [d , P] = Distance( A , B , C)
% DISTANCE Calculate the distance d of the point C from the
line segment AB
%
% unfortunatly I've no idea whom to credit for this
function...

% if size(A,1)~=2
% A = A';
% end

A=A(:);
B=B(:);
C=C(:);

A(3,:) = 0;
B(3,:) = 0;
C(3,:) = 0;

u = B-A;
s = dot( (C-A),u )/norm(u)^2;

if s<0
    d = norm(A-C);
    P = A;
elseif s<1
    d = abs(cross((C-A),(B-A)))/norm(B-A);
    d = d(3);
    P = A + s*u;
else
    d = norm(B-C);
    P = B;
end

function H=circle(center,radius,NOP,style)
%-----------------------------------------------------------
----------------------------------
% H=CIRCLE(CENTER,RADIUS,NOP,STYLE)
% This routine draws a circle with center defined as
% a vector CENTER, radius as a scaler RADIS. NOP is
% the number of points on the circle. As to STYLE,
% use it the same way as you use the rountine PLOT.
% Since the handle of the object is returned, you
% use routine SET to get the best result.
%
% Usage Examples,
%
% circle([1,3],3,1000,':');
% circle([2,4],2,1000,'--');
%
% Zhenhai Wang <zhenhai@ieee.org>
% Version 1.00
% December, 2002
%-----------------------------------------------------------
----------------------------------

if (nargin <3),
 error('Please see help for INPUT DATA.');
elseif (nargin==3)
    style='b-';
end;
THETA=linspace(0,2*pi,NOP);
RHO=ones(1,NOP)*radius;
[X,Y] = pol2cart(THETA,RHO);
X=X+center(1);
Y=Y+center(2);
H=plot(X,Y,style);
axis square;
<SNAP>