# Find biggest Circle having all points close to a List of Polygons

## Setting

Imagine a List of polygons polys (type: BoundaryMeshRegion, polygons may also overlap):

pnts = {{0.3, 0}, {4, 2}, {4, -2}};
rotpnts = Table[RotationTransform[2*Pi/3*i][pnts], {i, 0, 2}];
polys = BoundaryMeshRegion[#, Line[{1, 2, 3, 1}]] & /@ rotpnts;
Show[polys] ## Question

What is the radius of the biggest circle (center in {0,0}), whose inner points do not exceed a given (i.e. freely definable) distance dmaxto at least one of the polygons of polys?

Ideas

• The circle may continuously grow and then checked, if still every point on it is closer than dmax to one of the polygons
• "some" Voronoi approach ## Extra Question (solved)

Return the indices of the polys-list of the n closest polygons to a given point p.

Illustration For n=2 I would get a list of indices {3,1}, as polys[] is the nearest polygon to p and polys[] the second closest.

Illustrating Manipulate - Solution

pnts = {{0.3, 0}, {4, 2}, {4, -2}};
rotpnts = Table[RotationTransform[2*Pi/3*i][pnts], {i, 0, 2}];
polys = BoundaryMeshRegion[#, Line[{1, 2, 3, 1}]] & /@ rotpnts;

Manipulate[
testPt = p;
allPts = RegionNearest[#, p] & /@ polys;
mPts = allPts //
TakeSmallestBy[Norm[p - #] &, 3];
Show[polys,
Table[
Graphics[{
Blend[{Green, Red}, (i - 1)/(Length[mPts] - 1)]
, Line[{p, mPts[[i]]}]}]
, {i, 1, Length[mPts]}
]
]
, {{p, {0, 0}}, Locator}
]

Dynamic[Table[Norm[testPt - mPts[[i]]], {i, 1, Length[mPts]}]]
Dynamic[Ordering[Norm[testPt - #] & /@ allPts, 3]]

• Something like a circumcircle? Jul 17, 2016 at 11:42
• No, I want to define the maximum distance dmax of any inner points of the circle. The resulting circle may be smaller than the extents of the polygons...
– DPF
Jul 17, 2016 at 11:44
• And what have you tried so far? Jul 17, 2016 at 14:47
• @MarcoB I was thinking of a Voronoi-Approach, but I lack any experience with the application of it.
– DPF
Jul 18, 2016 at 6:34
• If you have not done so already I recommend that you ask this question on Mathematics for suggestions on the algorithm itself, which we can then help you implement optimally in Mathematica. Jul 18, 2016 at 12:37

# Final

pnts = {{0.3, 0}, {4, 2}, {4, -2}};
rotpnts =     Table[RotationTransform[2*Pi/3*i][pnts], {i, 0, 2}];
polys = BoundaryMeshRegion[#, Line[{1, 2, 3, 1}]] & /@ rotpnts;
reg = RegionUnion[polys]

regnp=ImplicitRegion[RegionDistance[reg,{x,y}]<2,{x,y}]

regn=NDSolveFEMToElementMesh[regnp,   MaxCellMeasure -> {"Length" -> .1}, "MaxBoundaryCellMeasure" -> .001]

-SignedRegionDistance[regn,{0,0}]


The following code can generate you a correct result, but at a low speed:

pnts = {{0.3, 0}, {4, 2}, {4, -2}};
rotpnts = Table[RotationTransform[2*Pi/3*i][pnts], {i, 0, 2}];
polys = BoundaryMeshRegion[#, Line[{1, 2, 3, 1}]] & /@ rotpnts;
reg = RegionUnion[polys]

c[r_?NumericQ] :=
NMaximize[
RegionDistance[reg, {x, y}], {x, y} \[Element] Circle[{0, 0}, r],
Method -> "NelderMead", PrecisionGoal -> 6, WorkingPrecision -> 10][]
NMinimize[{r, c[r] > 2, 2<=r<=5}, r,Method -> "SimulatedAnnealing", PrecisionGoal -> 4, WorkingPrecision -> 10]

(*3.5439*)


The constraint is written in c[r]>2, you can customize the function to make it wider applicable.

The main idea of the code is quite simple, for each curve with coefficient r, in this case a circle with radius r, we can find out the largest distance we can get by moving the point on the curve by RegionDistance and NMaximize.

I strongly suggest adding Method->"NelderMead" as sometimes other method like "DifferentialEvolution" will find the distance always 0 (the searching point keeps in the region and never get out.) and stop evaluation.

Then we can find out the shortest r we need by NMinimize.

There's still one problem left, If we use StepMonitor to track the process, we can find out that no matter which method we choose, the process will not be ideal, so I think maybe we can create a simple form of NMaximize or so specially designed for this problem. Also, we can find out that sometimes the value has already reached the most optimized value but it still won't generate an output.

# Edit1

The second problem can be solved by ConnectedMeshComponents and RegionNearest.

res=NMaximize[ RegionDistance[reg, {x, y}], {x, y} \[Element] Circle[{0, 0}, 3.5439(*the result*)], Method -> "NelderMead", PrecisionGoal -> 6, WorkingPrecision -> 10]
pt={x,y}/.res[];
dis=res[];
near=RegionNearest[#,pt]&/@ConnectedMeshComponents@DiscretizeRegion@RegionIntersection[Disk[pt,dis+.1],reg]
Show[RegionPlot[reg],Graphics[{Line[{pt,#}]&/@near,Red,Point@pt,Blue,Point/@near}]] # Edit2

I think you've got an answer by yourself already. This code can do the job you want:

pt = {3, 3};
n = 2;
near = RegionNearest[#, pt] & /@ ConnectedMeshComponents@reg
sel = TakeSmallestBy[near, EuclideanDistance[pt, #] &, n]
Show[RegionPlot[reg],
Graphics[{Line[{pt, #}] & /@ sel, Red, Point@pt, Blue,
Point /@ sel}]] # Edit3

Here gives a faster code optimized for this occasion:

pnts = {{0.3, 0}, {4, 2}, {4, -2}}; rotpnts = Table[RotationTransform[2*Pi/3*i][pnts], {i, 0, 2}]; polys = BoundaryMeshRegion[#, Line[{1, 2, 3, 1}]] & /@ rotpnts; reg = RegionUnion[polys]

c[r_?NumericQ] :=
NMaximize[
SignedRegionDistance[reg, {x, y}], {x, y} \[Element]
Circle[{0, 0}, r], Method -> "NelderMead", PrecisionGoal -> 6,
WorkingPrecision -> 10][]
find[tar_, range : {min_, max_}, prec_: .0001] :=
Round[NestWhile[
If[c@Mean@# > tar, Most, Rest][Subdivide[#[], #[], 2]] &,
N@range, #[] - #[] >= prec &][], prec]
find[2, {2, 5}]

(*3.5439*)


Use simple and violent method to find out the smallest r in a given range, but this method can only apply when c[r] is monotony in the given range.

In this case, this is correct, so it can generate a correct result in a short amout of time.

Will this help?

# Edit4

A totally new method: inflation and calculate the minimum distance.

regn=Quiet@DiscretizeRegion[ImplicitRegion[RegionDistance[reg,{x,y}]<2,{x,y}],MaxCellMeasure->0.001]

-SignedRegionDistance[regn,{0,0}]


This method is widely applicable as long as your precision goal is not that high. Firstly, we create a new region to include every point in a distance to the region, then we find out at which distance finally there's a point not satisfying that, then this distance is just what you want!

Or maybe this will be even better in performance?

regnp=ImplicitRegion[RegionDistance[reg,{x,y}]<2,{x,y}]

regn=NDSolveFEMToElementMesh[regnp,   MaxCellMeasure -> {"Length" -> .1}, "MaxBoundaryCellMeasure" -> .001]

• atm, I'm not at a PC, but I wonder how r is used in your c[r_]...
– DPF
Jul 17, 2016 at 19:33
• @DPF erer, sorry, I copied the test version, it should be Circle[{0,0},r]
– Wjx
Jul 17, 2016 at 22:53
• @DPF corrected~ now this code works. :) extra is not hard as well, I will try to update soon
– Wjx
Jul 17, 2016 at 22:54
• @DPF I'm currently without my computer, but I think maybe SignedRegionDistance will work better.
– Wjx
Jul 17, 2016 at 23:05
• @DPF then you may try RegionNearest and TakeSmallest.
– Wjx
Jul 18, 2016 at 7:18