12
$\begingroup$

This question asks to find the largest circle within a convex polygon, specifically this one:

Graphics[
 {Line[{{0,1}, {0,6}, {4,10}, {8,10}, {11,7}, {11,4}, {7,0}, {1,0}, {0,1}}]}
]

enter image description here

One can optimize over the $(x,y)$ location of its center maximizing its radius $r$, but this seems quite a non-Mathematica approach. Presumably there is a method based on RegionFunctions, Centroids, and such to exploit Mathematica's region functionality.

Alas, the natural primitives such as Centroid and all the RegionMeasures and RegionMemberQ I could think of didn't help.

I have a feeling there is a one-line solution to this, but cannot find it. It would be especially great to get an analytic solution (based on the equations of the polygon edges, say), rather than just some numerical center location and radius.

$\endgroup$
6
  • 1
    $\begingroup$ I think the answers here might be helpful. I get a good result using inCircle from this answer. $\endgroup$ – Jason B. Oct 6 '20 at 1:19
  • $\begingroup$ Partly related too mathematica.stackexchange.com/q/226058/72682 though the method there was to test random seed points and ovoid other obstacles packing the interior. $\endgroup$ – flinty Oct 6 '20 at 1:21
  • $\begingroup$ I wouldn't be so certain it will be a one-liner. Looking at this site I'm struck by this sentence: "the solution will contact three edges". That seems true, but I can't say for sure. But if it is true - finding the incircle of a triangle is easy and analytic. Just find the incircle for all sets of three edges, and grab the largest one. $\endgroup$ – Jason B. Oct 6 '20 at 1:52
  • $\begingroup$ The solution disk must touch at least three edges of the convex polygon. If the two touched sides are not parallel, one can always displace the disk one of the two candidate directions to enable a larger radius. In the (rare) case the touched sides are parallel (as in a long rectangle), then one can shift the solution one direction until it touches one of the remaining sides (even though it does not change the solution area). $\endgroup$ – David G. Stork Oct 6 '20 at 1:55
  • $\begingroup$ @DavidG.Stork Very nice problem. Perhaps there exists a unique "minimal" solution with minimal distance between the centroids of disk and polygon? $\endgroup$ – Ulrich Neumann Oct 6 '20 at 8:37
14
$\begingroup$
poly = Polygon @ {{0, 1}, {0, 6}, {4, 10}, {8, 10}, {11, 7}, {11, 4}, 
   {7, 0}, {1,  0}, {0, 1}};

dsk = Disk[{x, y}, r];

We can use RegionWithin[poly, dsk] as the constraint in ArgMax:

sol = Quiet @ ArgMax[{r, RegionWithin[poly, dsk]}, {x, y, r}] 
{43/8, 39/8, 13/(2 Sqrt[2])}  
Graphics[{EdgeForm[Black], FaceForm[Yellow], poly, 
   Red, Circle[Most @ sol, Last @sol], PointSize[Large], Point[Most @ sol]}] 

enter image description here

$\endgroup$
9
  • $\begingroup$ Great, that's the one-liner ! I tried numerically NArgMax[{ r, RegionWithin[poly, dsk]}, {x, y, r}] but MMA doesn't evaluate. Any idea? Thanks. $\endgroup$ – Ulrich Neumann Oct 6 '20 at 8:00
  • $\begingroup$ @UlrichNeumann, I also tried NArgMax and NMaximize and had to abort because it was taking too long. I don't know why. $\endgroup$ – kglr Oct 6 '20 at 8:20
  • $\begingroup$ Mystery MMA... Thanks! $\endgroup$ – Ulrich Neumann Oct 6 '20 at 8:32
  • 3
    $\begingroup$ Yep... that counts as a one-liner! I certainly didn't know about RegionWithin. Once again the mysterious kglr solves what few others can. ($\checkmark$) $\endgroup$ – David G. Stork Oct 6 '20 at 15:58
  • $\begingroup$ @kglr I'm wondering , the code sol=... doesn't evaluate anymore in Mathematica v12.2, it only shows ArgMax[...] with several conditions like (Indeterminate | Indeterminate) \[Element] Reals && Indeterminate <= -1 . $\endgroup$ – Ulrich Neumann Mar 15 at 12:58
10
$\begingroup$

As an area maximization problem:

reg = Polygon[{{0, 1}, {0, 6}, {4, 10}, {8, 10}, {11, 7}, {11, 4}, {7, 0}, {1, 0}, {0, 1}}];
rnf = RegionNearest[RegionBoundary[reg]];
gendisk[{x_, y_}] := Disk[{x, y}, EuclideanDistance[{x, y}, rnf[{x, y}]]]
cost[{x_?NumericQ, y_?NumericQ}] := Area[gendisk[{x, y}]]
{err, sol} = NMaximize[cost[{x, y}], {x, y} \[Element] reg,
  Method -> "RandomSearch"];
Graphics[{FaceForm[None], EdgeForm[Black], reg, Yellow, 
  FaceForm[Yellow], gendisk[Values@sol], Red, Point[Values[sol]]}]

disk in convex

There's a ridge at the centre of the distance transform, so I don't think there is a unique solution but a family of disks with maximal area.

ImageAdjust@DistanceTransform@ColorNegate@Graphics[reg]

distance transform

$\endgroup$
1
  • $\begingroup$ This is nice, but ultimately equivalent to the general form I considered: Searching through $(x,y)$ while keep $r$ as large as possible. (+1)... but I have a feeling a more elegant solution remains to be shown. $\endgroup$ – David G. Stork Oct 6 '20 at 1:32
3
$\begingroup$

Conjecture

We believe that the center of the max circle should be lie in the line segment.

pts = {{0, 1}, {0, 6}, {4, 10}, {8, 10}, {11, 7}, {11, 4}, {7, 0}, {1,
     0}};
poly = Polygon[pts];
fig1 = Graphics[{{LightGreen, poly}, {Red, Point[pts]}, Blue, 
    Text[#, RegionCentroid[RegionDifference[Disk[#, 1.3], poly]]] & /@
      pts}];
p1 = {x, y} /. (Reduce[
       RegionDistance[InfiniteLine[{{0, 6}, {4, 10}}], {x, y}] == 
         RegionDistance[InfiniteLine[{{7, 0}, {11, 4}}], {x, y}] == 
         RegionDistance[InfiniteLine[{{1, 0}, {7, 0}}], {x, y}] && {x,
           y} ∈ poly, Reals] // ToRules) // Simplify;
p2 = {x, y} /. (Reduce[
       RegionDistance[InfiniteLine[{{0, 6}, {4, 10}}], {x, y}] == 
         RegionDistance[InfiniteLine[{{7, 0}, {11, 4}}], {x, y}] == 
         RegionDistance[
          InfiniteLine[{{4, 10}, {8, 10}}], {x, y}] && {x, 
          y} ∈ poly, Reals] // ToRules) // Simplify;
Show[fig1, 
 Graphics[{Text[p1, p1, {1, 1}], 
   Text[p2, p2, {-1, -1}], {Red, Point[{p1, p2}], Line[{p1, p2}]}}]]

enter image description here

pts = {{0, 1}, {0, 6}, {4, 10}, {8, 10}, {11, 7}, {11, 4}, {7, 0}, {1,
     0}};
poly = Polygon[pts];
bds = InfiniteLine /@ Partition[pts, 2, 1, 1];
Maximize[{Min[RegionDistance[#, {x, y}] & /@ bds], {x, y} ∈ 
     poly}, {x, y}] // Simplify;
ContourPlot[
 Min[RegionDistance[#, {x, y}] & /@ bds], {x, y} ∈ poly, 
 Contours -> {1, 1.5, 2.5, 3.5, 4, 4.5}, ContourShading -> Automatic, 
 PlotPoints -> 50, MaxRecursion -> 2]

enter image description here

Edit II

pts = {{0, 1}, {0, 6}, {4, 10}, {8, 10}, {11, 7}, {11, 4}, {7, 0}, {1,
     0}};
poly = Polygon[pts];
bds = InfiniteLine /@ Partition[pts, 2, 1, 1];
Maximize[{Min[RegionDistance[#, {x, y}] & /@ bds], {x, y} ∈ 
    poly}, {x, y}] // Simplify

$$\left\{\frac{13}{2 \sqrt{2}},\left\{x\to \frac{1}{4} \left(42-13 \sqrt{2}\right),y\to 10-\frac{13}{2 \sqrt{2}}\right\}\right\}$$

pts = {{0, 1}, {0, 6}, {4, 10}, {8, 10}, {11, 7}, {11, 4}, {7, 0}, {1,
     0}};
poly = Polygon[pts];
bds = InfiniteLine /@ Partition[pts, 2, 1, 1];
sol = Maximize[{r, 
    Sequence @@ 
     Table[EuclideanDistance[RegionNearest[bd, {x, y}], {x, y}] >= 
       r, {bd, bds}], {x, y} ∈ poly}, {r, x, y}] // Simplify
Graphics[{{Opacity[0.1], poly}, Point[{x, y}], Circle[{x, y}, r]} /. 
   Last[sol]] // Timing

enter image description here $$\left\{\frac{13}{2 \sqrt{2}},\left\{r\to \frac{13}{2 \sqrt{2}},x\to \frac{11}{2},y\to 5\right\}\right\}$$

Edit I

pts = {{0, 1}, {0, 6}, {4, 10}, {8, 10}, {11, 7}, {11, 4}, {7, 0}, {1,
     0}};
poly = Polygon[pts];
bds = InfiniteLine /@ Partition[pts, 2, 1, 1];
sol = Maximize[{r, 
    Table[RegionDistance[bd, {x, y}] >= r, {bd, bds}], {x, 
      y} ∈ poly}, {r, x, y}] // Simplify
Graphics[{{Opacity[0.2], poly}, Point[{x, y}], Circle[{x, y}, r]} /. 
  Last[sol]]

$$\left\{\frac{13}{2 \sqrt{2}},\left\{r\to \frac{13}{2 \sqrt{2}},x\to \frac{43}{8},y\to \frac{39}{8}\right\}\right\}$$

If we add condition to y such as y>=5,the result is

$$\left\{\frac{13}{2 \sqrt{2}},\left\{r\to \frac{13}{2 \sqrt{2}},x\to \frac{355}{64},y\to \frac{323}{64}\right\}\right\}$$

So It must be exist a line attain the maximum.

$\endgroup$
1
  • $\begingroup$ This is great (+1)... far more than was required by my request for a "one-liner." I knew the solution admitted the circle center being on a line (as can be found by Skeletonization) but your approach is analytic, which is almost always superior. Thanks. $\endgroup$ – David G. Stork Jan 15 at 19:01
3
$\begingroup$

I played around a bit with this problem and noticed that all numeric functionalities become significantly faster when you use a polygon with machine-precision points. So if you want to use NMaximize or the like, I highly recommend this.

Furthermore, here's a numerical implementation that uses SignedRegionDistance:

Clear[x];
poly = Polygon @ N[{{0, 1}, {0, 6}, {4, 10}, {8, 10}, {11, 7}, {11, 4}, {7, 0}, {1, 0}, {0, 1}}];
sol = With[{
   sgnDist = SignedRegionDistance[poly]
  },
  NMinimize[sgnDist[x], x \[Element] poly]
]
Graphics[
 {
  {Red, Disk[x /. Last[sol], Abs[First @ sol]]},
  {FaceForm[None], EdgeForm[Black], poly}
 }
]

enter image description here

$\endgroup$
3
  • $\begingroup$ I also test SignedRegionDistance, your approach is best than my original method. $\endgroup$ – cvgmt Jan 15 at 9:34
  • $\begingroup$ Thanks (+1). A helpful improvement for larger problems. $\endgroup$ – David G. Stork Jan 15 at 16:38
  • $\begingroup$ @DavidG.Stork There are some more optimizations you can make, potentially. For ConvexPolygonQ polygons, you can use FindMinimum using RegionCentroid as an initial guess. That should be faster for more complicated polygons. $\endgroup$ – Sjoerd Smit Jan 15 at 17:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.