8
$\begingroup$

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}

enter image description here

n=100

{0.233274, Null} enter image description here

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

Improve this question
$\endgroup$
0

4 Answers 4

Reset to default
13
$\begingroup$

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}
Improve this answer
$\endgroup$
7
$\begingroup$
  • 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.

Improve this answer
$\endgroup$
7
  • 3
    $\begingroup$ For RandomChoice[triangles], you might want to weight each triangle by its area. This will let you sample uniformly across the full region. $\endgroup$ Commented Nov 18, 2024 at 18:18
  • $\begingroup$ @GregHurst Thanks! $\endgroup$ Commented Nov 18, 2024 at 18:22
  • $\begingroup$ @cvgmt very nice! Agree re: area weighting. Always learn from your answers! $\endgroup$ Commented Nov 19, 2024 at 8:56
  • $\begingroup$ Timing measurement incorrect. Could you include triangles and areas inside RepeatedTiming, as well as rectangle, pts, counts and inner ? $\endgroup$ Commented Nov 26, 2024 at 13:39
  • $\begingroup$ Just as compile method, we not always care about the compile time. $\endgroup$ Commented Nov 26, 2024 at 13:49
2
$\begingroup$

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.

Improve this answer
$\endgroup$
4
  • 1
    $\begingroup$ Posting a Mathematica code for your algo would be really nice of you. $\endgroup$ Commented Nov 21, 2024 at 9:35
  • $\begingroup$ @Anton I had Mathematica once and loved it, but that was at a previous job over 20 years ago. I am now on the free stuff. Octave is MATLAB, and in my haste to be helpful I read Mathematica as MATLAB. $\endgroup$ Commented Nov 21, 2024 at 9:57
  • $\begingroup$ I see. Check out my solution, instead of going perpendicular the side segment, I go bisector closest to 90 degrees. $\endgroup$ Commented Nov 23, 2024 at 12:38
  • 1
    $\begingroup$ It seems sensible. There are lots of ways of doing this. The winding count is very quick to calculate. I picked the first edge and invoked it once, so I can't see a way of getting it any faster, but there is a small chance that the second point won't be inside the shape if it is fantastically skinny. I don't say the point will be a particularly good inside point: if you had a circle then the point would be right next to the edge. But it is inside. $\endgroup$ Commented Nov 23, 2024 at 15:49
2
$\begingroup$

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:

enter image description here

Timings in regular scale:

enter image description here

Improve this answer
$\endgroup$
4
  • $\begingroup$ Why not test n=1000 and enclose all the code by RepeatedTiming? $\endgroup$ Commented Nov 28, 2024 at 7:10
  • $\begingroup$ My code is contained by functions and their execution is enclosed by RepeatedTiming and you're posting right under For n=1000 {0.00325786, Null} test result. LMFAO that's why. $\endgroup$ Commented Nov 28, 2024 at 8:13
  • $\begingroup$ isitclockwise=Compile... need time to compile, do you know? $\endgroup$ Commented Nov 28, 2024 at 8:39
  • $\begingroup$ no sweat, really? $\endgroup$ Commented Nov 28, 2024 at 9:16

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.