Generate A Point Inside An Arbitrary Concave Polygon

Question

Imagine I have an arbitrary concave polygon. I need to find/generate a point inside of it, efficiently.

RandomPoint@Polygon@points is relatively fast only for simple such polygons.

n=1;(* toy polygon complexity *)

points=Join[Transpose[{Range[2*n+1]*5,Drop[Flatten@ConstantArray[{0,45},n+1],-1]}],Transpose[{Reverse@Range[2*n+1]*5-5,Drop[Flatten@ConstantArray[{55,10},n+1],-1]}]];

insidepoint=RandomPoint@Polygon@points;//RepeatedTiming

Show[
Graphics@Thread@Line@Transpose[{points,RotateLeft@points}],
Graphics@{Red,PointSize[Medium],Point@insidepoint}
]

n=1

> {0.00244676, Null}

n=100

{0.233274, Null}

n=1000

{2.84268, Null}

So, for complex polygons,RandomPoint@Polygon is very, very slow.

Is there a method to find a point inside a convex polygon efficiently, no matter the complexity?

UPD.

I’ve added timings here:

timings

Answer 1

We can use mesh regions instead. For $n = 1000$ here's timings on my machine:

SeedRandom[1];
RandomPoint@Polygon[points]; // AbsoluteTiming

{0.930068, Null}

SeedRandom[1];
RandomPoint@MeshRegion[points, 
  Polygon[Range[Length[points]]]]; // AbsoluteTiming

{0.018052, Null}

SeedRandom[1];
RandomPoint@BoundaryMeshRegion[points, 
  Line[Append[Range[Length[points]], 1]]]; // AbsoluteTiming

{0.004879, Null}

Answer 2

• decomposition the polygon into several triangles (or convex set) before select points.

triangles = PolygonDecomposition[Polygon@points, "Triangle"];
areas = Area /@ triangles;

(pt = RandomPoint@RandomChoice[areas -> triangles]) // RepeatedTiming
RegionMember[Polygon@points, pt]

True.

Answer 3

I posted a hasty answer which without having Mathematica. I have deleted that one. Here is a better version if someone else wants a go.

Take two points A, B on a polygon side. If A is B and the side has no length, go on to the next side.

Get the mid-point M = (A+B)/2;

Get a small perpendicular vector V = (A-B)*[0,1;-1,0]*0.0001:

Here are two points just off the mid-point...

P0 = M+V;

P1 = M-V;

If the winding number for P0 is non-zero then P0 is inside the polygon. If it is zero then P1 is inside the polygon unless you are really unlucky.

Answer 4

Here's a fast algo: 1) find convex vertices, 2) from those find one with the maximum Cross between its segments, 3) then pick a point using bisector and a small distance. This approach remotely resembles @Richard Kirk ‘s idea.

Testing code:

(* Generate a convex polygon *)
    n = 1;(* 1 | 100 | 1000 | 5000 | 10000 *)
points = 
  Join[Transpose[{Range[2*n + 1]*5, 
     Drop[Flatten@ConstantArray[{0, 45}, n + 1], -1]}], 
   Transpose[{Reverse@Range[2*n + 1]*5 - 5, 
     Drop[Flatten@ConstantArray[{55, 10}, n + 1], -1]}]];

(* OP slow solution *)
{t0, pt0} = RandomPoint@Polygon[points] // RepeatedTiming;

(* Greg  Hurst 1 - RandomPoint@MeshRegion *)
{t1, pt1} = 
  RandomPoint@MeshRegion[points, Polygon[Range[Length[points]]]] // 
   RepeatedTiming;
(* Greg  Hurst 2 - RandomPoint@
    BoundaryMeshRegion *)
{t2, pt2} = 
  RandomPoint@
    BoundaryMeshRegion[points, 
     Line[Append[Range[Length[points]], 1]]] // RepeatedTiming;

(* cvgmt *)
{t3, nothing} = RepeatedTiming[
   triangles = PolygonDecomposition[Polygon@points, "Triangle"];
   areas = Area /@ triangles;
   pt3 = RandomPoint@RandomChoice[areas -> triangles];
   ];

(* Anton *)
distance = 1*^-9;

generateInsidePoint[points_, distance_] := 
 Module[{segVectors, svNorm, svNormRL, cross, rotations, ccws, 
   ccwcrosses, vertexPos, vertex}, 
  segVectors = RotateLeft@points - points;
  svNorm = 
   segVectors/(segVectors[[All, 1]]^2 + segVectors[[All, 2]]^2)^.5;
  svNormRL = RotateLeft@svNorm;
  cross = 
   Abs[-svNorm[[All, 2]]*svNormRL[[All, 1]] + 
     svNorm[[All, 1]]*svNormRL[[All, 2]]];
  rotations = isitclockwise@points;
  ccws = Pick[Range@Length@points, rotations, -1];
  ccwcrosses = cross[[ccws]];
  vertexPos = ccws[[Ordering[ccwcrosses][[-1]]]];
  vertex = points[[vertexPos]];
  vertex + 
   Normalize[
     Mean@{vertex - 
         Normalize[vertex - RotateRight[points][[vertexPos]]], 
        vertex + 
         Normalize[RotateLeft[points][[vertexPos]] - vertex]} - 
      vertex]*distance]

generateInsidePoint2[points_, distance_] := 
 Module[{segVectors, svNorm, svNormRL, cross, rotations, ccws, 
   ccwcrosses, vertexPos, vertex, rtL, rtR,
   triples, t1, t2},
  
  segVectors = RotateLeft@points - points;
  svNorm = 
   segVectors/(segVectors[[All, 1]]^2 + segVectors[[All, 2]]^2)^.5;
  svNormRL = RotateLeft@svNorm;
  cross = 
   Abs[-svNorm[[All, 2]]*svNormRL[[All, 1]] + 
     svNorm[[All, 1]]*svNormRL[[All, 2]]];
  
  rtR = RotateRight@points;
  rtL = RotateLeft@points;
  triples = Transpose[{rtR, points, rtL}];
  t1 = triples[[All, All, 1]];
  t2 = triples[[All, All, 2]];
  rotations = 
   Sign[Total /@ ((RotateLeft /@ t1 - t1)*(RotateLeft /@ t2 + t2))];
  
  ccws = Pick[Range@Length@points, rotations, -1];
  ccwcrosses = cross[[ccws]];
  vertexPos = ccws[[Ordering[ccwcrosses][[-1]]]];
  vertex = points[[vertexPos]];
  
  vertex + 
   Normalize[
     Mean@{vertex - 
         Normalize[vertex - RotateRight[points][[vertexPos]]], 
        vertex + 
         Normalize[RotateLeft[points][[vertexPos]] - vertex]} - 
      vertex]*distance
  
  ]

(* Anton 1 - we do this one for those of us who speaks of themselves in the plural *)
{t4, nothing} = RepeatedTiming[
   isitclockwise = 
    Compile[{{pts, _Real, 2}}, Module[{rtR, rtL, triples, t1, t2},
      rtR = RotateRight@pts;
      rtL = RotateLeft@pts;
      triples = Transpose[{rtR, pts, rtL}];
      t1 = triples[[All, All, 1]];
      t2 = triples[[All, All, 2]];
      Sign[
       Total /@ ((RotateLeft /@ t1 - t1)*(RotateLeft /@ t2 + t2))]], 
     RuntimeOptions -> "Speed"];
   pt4 = generateInsidePoint[points, 1*^-9];
   ];

(* Anton 2 - using compiled function but not being dumb af *)
{t5, pt5} = generateInsidePoint[points, 1*^-9] // RepeatedTiming;

(* Anton 3 - without compiled functions *)
{t6, pt6} = generateInsidePoint2[points, 1*^-9] // RepeatedTiming;

Timings in log scale:

Timings in regular scale: