How to find the incircle and circumcircle for an irregular polygon?

Question

I was inspired by blog post where Gonzalo Ciruelos sorted all the countries by roundness, and looked at the Wikipedia page that gives a different definition of roundness to what Ciruelos used:

The ISO definition of roundness is based on the ratio between the inscribed and the circumscribed circles, i.e. the maximum and minimum sizes for circles that are just sufficient to fit inside and to enclose the shape

So then the question is, can we use Mathematica to find the incircle and circumcircle for a given polygon?

I can do something like this by hand just using Manipulate and a Locator, getting these results,

But I need some programmatic way to find the circles.

Answer 1

For the in-circle aspect of the question, we can use the fact that the center of the in-circle is the Voronoi node with the most negative signed distance to the polygon. Finding the Voronoi mesh and the negative signed distance of all of the Voronoi nodes are very fast:

incircle[poly_List?MatrixQ] := With[
    {nodes = VoronoiMesh[poly]["Coordinates"], sd = SignedRegionDistance[Polygon[poly]]},
    With[
        {ord = First @ Ordering[sd[nodes], 1]},

        Circle[
            nodes[[ord]],
            Abs @ sd[nodes[[ord]]]
        ]
    ]
]

Let's compare timings with the accepted answer:

poly = Entity["Country", "France"]["Polygon"][[1,1,1]];

inCircle[poly] //AbsoluteTiming
incircle[poly] //AbsoluteTiming

{2.41236, {{46.7076, 2.52144}, 3.47992}}

{0.116144, Circle[{46.7076, 2.52146}, 3.47992]}

Answer 2

With the help of J.M. and Chip Hurst, I was able to get this up and running with a reasonably fast (meaning it's still slow) and clean code. The problem of the circumcircle is solved by the new-in-version-10.4 function BoundingRegion, and the incircle is found by using NMinimize along with SignedRegionDistance.

In the code below, I overload the definitions to account for the fact that when you call CountryData[<country>,"Polygon"] you get your coordinates wrapped in GeoPosition, where the x and y coordinates are backwards. But if you call CountryData[<country>,{"Polygon","Mercator"}] then it's a normal polygon.

I also wrote the code so that it will find the total circumcircle for a collection of polygons, so it takes account for all the land masses in a disjointed geographic region.

ClearAll[inCircle, outCircle];
inCircle[Polygon[GeoPosition[a_]]] := 
  inCircle[Polygon[a /. {x_?NumericQ, y_?NumericQ} :> {y, x}]];
inCircle[Polygon[a_ /; Depth[a] == 4]] := 
  Last@(SortBy[inCircle /@ a, Last]);
inCircle[a_List] := inCircle[Polygon@a];
inCircle[Polygon[a_ /; Depth[a] == 3]] := inCircle[Polygon[a]] =
  Module[{interior},
    interior = BoundaryDiscretizeGraphics@Polygon[a];
    {{x, y} /. #2, Abs@#1} & @@ 
     NMinimize[
      SignedRegionDistance[interior, {x, y}], {x, y} ∈ 
       interior]
    ];

outCircle[p_List] := List @@ BoundingRegion[p, "MinDisk"];
outCircle[Polygon[GeoPosition[a_]]] := 
  outCircle[Polygon[a /. {x_?NumericQ, y_?NumericQ} :> {y, x}]];
outCircle[Polygon[a_]] := outCircle[Flatten[a, Depth[a] - 3]];

Just for educational value, I also made a version of this that runs on Mathematica 9, posted as a separate answer. Here is the result for Japan,

Show[Graphics[{Blue, Disk @@ outCircle[#]}], Graphics@#, 
   Graphics[{Red, Disk @@ inCircle[#]}]] &@
 CountryData["Japan", {"Polygon", "Mercator"}]

And here I sort the countries by their roundness:

countries = CountryData[#, "StandardName"] & /@ CountryData[];
resultsList = Table[
      With[{pol = CountryData[country, {"Polygon", "Mercator"}]},
       {CountryData[country, "Name"], Last[inCircle[pol]]^2/
        Last[outCircle[pol]]^2, 
        Show[Graphics[{Blue, Disk @@ outCircle[pol]}], Graphics@pol, 
         Graphics[{Red, Disk @@ inCircle[pol]}]]}],
      {country, countries}] // SortBy[#[[2]] &] // 
    Reverse;~Monitor~country
resultsList[[;;20]] // TableForm

The 20 roundest countries are below, and the full list is posted on gist, and can be downloaded via << "https://git.io/v6mt2"

Answer 3

Let me present another approach. The code should be self-explanatory.

country = "Mauritius";
imgwidth = 500; (* increase this for better fit of inner circle *)
allpols = Polygon /@ First @ CountryData[country, {"Polygon", "Mercator"}];
areas = Area /@ allpols;
pol = First @ Pick[allpols, areas, Max[areas]]; (* just the biggest connected land *)
polpoints = First @ pol;
img = Composition[
  ImageCrop,
  Image[#, ImageSize -> imgwidth] &,
  Graphics
] @ pol
skeleton = Composition[
  ImageData,
  ImageAdjust,
  SkeletonTransform,
  ColorNegate
] @ img;
Image @ skeleton (* just to show the skeleton *)
centre = FirstPosition[skeleton, 1.] (* in 'List coordinates' *)
dim = Dimensions[skeleton]
newcentre = {Last @ centre, First @ dim - First @ centre} (* in Cartesian coordinates *)
{xmin, xmax} = Through[{Min, Max}[First /@ polpoints]]
ymin = Min[Last /@ polpoints]
scalefactor = First @ ImageDimensions[img] / First @ Differences[{xmin, xmax}]
scaler = ScalingTransform[ConstantArray[scalefactor, 2], {0, 0}]
translationvector = - {xmin, ymin}
translator = TranslationTransform[translationvector]
newpolpoints = (scaler @* translator) /@ polpoints;
boundary = MeshRegion[
  newpolpoints,
  Line /@ Partition[Append[Range[Length @ newpolpoints], 1], 2, 1]
];
innerradius = RegionDistance[boundary, newcentre]
outercircle = BoundingRegion[newpolpoints, "MinDisk"]
outerradius = Last @ outercircle
Graphics[{
  Blue,
  outercircle,
  White,
  Polygon[newpolpoints],
  Red,
  Disk[newcentre, innerradius]
}]
roundness = (innerradius/outerradius)^2;
Print["Roundness of ", country, ": ", roundness]

Notes

1. pol is taken to be the biggest connected land of a country, in the selected projection. For Japan in the Mercator projection, for example, pol is Hokkaido (which is actually smaller than Honshu). You may want to choose a representative area differently.

2. For some countries like Fiji, you would really need to crank up imgwidth.

Answer 4

At Rahul's suggestion I looked into the "smallest circle problem", which covers the circumcircle. And for the incircle I needed to read about the "largest empty circle".

I wrote some code that will find the necessary circles for version 9, which doesn't have the nifty built-in functions that I use above. I borrowed a function from cormullion here, and adapted some code from DrBob here which uses an algorithm from here.

The circumcircle is found via a recursive algorithm that looks at each point and decides whether it is inside the circle so far, if not then it will be part of the boundary points.

For the incircle, it is trickier. I use the fact that for a convex polygon, the center of the incircle will be on one of the vertices of the Voronoi diagram of the polygon vertices (this may not be generally true, I'm awful at proofs, but it seems to work here). Generating these vertices takes a long time using the ComputationalGeometry package, but it is doable.

Needs["ComputationalGeometry`"];
outCircle[Polygon[a_]] := outCircle[Flatten[a, Depth[a] - 3]];
outCircle[points_List] := 
  outCircle[points] = 
   Module[{calcDisk, moveToFront, circleThrough3Points},
    circleThrough3Points[{p1_, p2_, p3_}] := 
     Module[{ax, ay, bx, by, cx, cy, a, b, c, d, e, f, g, centerx, 
       centery, r}, {ax, ay} = p1;
      {bx, by} = p2;
      {cx, cy} = p3;
      a = bx - ax;
      b = by - ay;
      c = cx - ax;
      d = cy - ay;
      e = a (ax + bx) + b (ay + by);
      f = c (ax + cx) + d (ay + cy);
      g = 2 (a (cy - by) - b (cx - bx));
      If[g == 0, 
       False, {centerx = (d e - b f)/g, centery = (a f - c e)/g, 
        r = Sqrt[(ax - centerx)^2 + (ay - centery)^2]}];
      {{centerx, centery}, r}];
    calcDisk[{}] := {0, 0};
    calcDisk[{{x_, y_}}] := {0, 0};
    calcDisk[{{x_, y_}, {x2_, y2_}}] := {Mean[{{x, y}, {x2, y2}}], 
      Norm[{x, y} - {x2, y2}]/2};
    calcDisk[{{x_, y_}, {x2_, y2_}, {x3_, y3_}}] := 
     circleThrough3Points[{{x, y}, {x2, y2}, {x3, y3}}];
    moveToFront[pts_, b_] := Module[{copy = pts, c, rs},
      {c, rs} = calcDisk[b];
      If[Length@b != 3, 
       Do[If[Norm[copy[[n]] - c] > rs, {c, rs} = 
          moveToFront[copy[[;; n - 1]], Union[b, copy[[{n}]]]];
         copy = Prepend[Delete[copy, n], copy[[n]]]], {n, 
         Length@pts}]];
      {c, rs}];
    moveToFront[points, {}]];
inCircle[Polygon[a_ /; Depth[a] == 4]] := 
  Last@SortBy[inCircle /@ a, Last];
inCircle[Polygon[a_ /; Depth[a] == 3]] := inCircle[a];
inCircle[p_List] := inCircle[p] = Quiet@Module[{vdpoints, dist},

    vdpoints = 
     p // VoronoiDiagram // First // Cases[#, _List] & // 
      Select[#, (Graphics`Mesh`InPolygonQ[p, #] &)] &;
    dist[pt_] :=
     Min[EuclideanDistance[pt, #] & /@ p];
    {#, dist@#} &@Last[SortBy[vdpoints, dist]]
    ]

Here it is applied to a couple of examples (be patient, it takes a while). You can see the results match up reasonably well with those above.

pol = CountryData["Japan", {"Polygon", "Mercator"}];
Graphics[{Blue,
    Disk @@ outCircle[#], Black, #,
    Red, Disk @@ inCircle[#]}] &@pol
Last[inCircle[pol]]^2/Last[outCircle[pol]]^2

pol = CountryData["SierraLeone", {"Polygon", "Mercator"}];
Graphics[{Blue,
    Disk @@ outCircle[#], Black, #,
    Red, Disk @@ inCircle[#]}] &@pol
Last[inCircle[pol]]^2/Last[outCircle[pol]]^2