38
$\begingroup$

Bug introduced in 10.0 and have been fixed before 13.3.1
RegionMember is new in 10.0
A commenter notes that the issue is with Graphics`PolygonUtils`PolygonTriangulate[polygon]

RegionMember (new in V10) appears to have a bug. For some polygons it finds points that are clearly outside to be "members" of the region. This depends on the shape of the polygon in a way that is hard to pin down. Here is one example.

polygonPoints = {{0, 200},{100, 200},{100, 300},{200, 300},{200, 400},{250, 400},{250, 0},{200, 0},{200, 100},{0, 100}}; polygon = Polygon[polygonPoints]; polyQ = RegionMember[polygon] testPoints = Table[{300 Random[], 450 Random[]}, {10000}];

Next let's see which test points RegionMember finds to be inside the test polygon (blue points) and how it compares to actual polygon (gray region).

ListPlot[{testPoints, Select[testPoints, polyQ],polygonPoints}, PlotStyle -> {Red, Blue, {PointSize -> Large, Green}}, Prolog -> {Gray, polygon}]

A polygon and the points deemed to be inside it by RegionMember

Changing the polygon reveals a little bit about how the bug affects results in this particular case. The mistaken area appears to be below a line that connects two vertices; moving one of the vertices and redoing the test confirms this.

polygonPointsB = ReplacePart[polygonPoints, 1 -> {0, 150}]; polygonB = Polygon[polygonPointsB] polyBQ = RegionMember[polygonB] ListPlot[{testPoints, Select[testPoints, polyBQ], polygonPointsB}, PlotStyle -> {Red, Blue, {PointSize -> Large, Green}}, Prolog -> {Gray, polygonB}, Epilog -> {Green, Thick, Line[polygonPointsB[[{1, 4}]]]}]

A polygon and the points deemed to be inside it by RegionMember

Other observations based on a small number of cases (using V10.3):

  1. The mistaken region was always inside the convex hull of the actual polygon
  2. Problem goes away when using real numbers to define the polygon
Improve this question
$\endgroup$
7
  • 3
    $\begingroup$ With polygonPoints = N@polygonPoints works well. $\endgroup$ Commented Dec 19, 2015 at 16:04
  • 2
    $\begingroup$ polyQ // InputForm reveals a difference between the exact and machine precision versions. For the latter it uses a MeshRegion, but for the former it contains an implicit member function which is presumably not quite correct. $\endgroup$ Commented Dec 19, 2015 at 17:46
  • 17
    $\begingroup$ The incorrect result in the exact case is coming from an internal triangulation routine, Graphics`PolygonUtils`PolygonTriangulate[polygon]. $\endgroup$ Commented Dec 19, 2015 at 18:12
  • 3
    $\begingroup$ @ilian What a horrible bug. Has this issue been reported to WRI? $\endgroup$ Commented Jul 10, 2019 at 3:02
  • $\begingroup$ @AliHashmi Ok I just reported it. $\endgroup$ Commented Sep 8, 2019 at 2:13

1 Answer 1

Reset to default
4
$\begingroup$

This bug still exists in version 12.2.

There is a workaround, which uses an alternative method of creating the RegionMemberFunction. I've shown that below, but first let's look at a couple of other possibilities that don't pan out.

The first idea one might consider is to define an alternate region member function, as follows.

polyQ = RegionMember[polygon, #] &

This works, but it performs a couple orders of magnitude slower (presumably the reason RegionMemberFunction was invented), so it's not really an option.

Trying to fix the RegionMemberFunction directly doesn't work either, even though a workable repair seems apparent. Specifically, a RegionMemberFunction has three parts, the region specified in RegionMember, the dimensionality of the region (expressed as an integer), and a predicate (expressed as a pure Function). Given a point's coordinates, the predicate (part 3) returns True or False, according to whether or not the point lies in the region. In the cited case, evaluating and reducing the predicate makes the error clear (y<=(400+x)/2 should be y<=200).

polyQ = RegionMember[polygon]; Reduce@polyQ[[3]][{x, y}]

(0≤x<100 && 100≤y≤(400+x)/2) || (100≤x<200 && 100≤y≤300) || (200≤x≤250 && 0≤y≤400)

A correct predicate can be obtained as follows, where the expression here could be used as the body of a Function.

Element[(x|y),Reals] && Reduce@RegionMember[polygon,{x,y}]

(x|y)∈ℝ && ((0≤x<100 && 100≤y≤200) || (100≤x<200 && 100≤y≤300) || (200≤x≤250 && 0≤y≤400))

However, replacing the predicate in a RegionMemberFunction with a corrected version makes no difference, because the predicate seems not to be used when RegionMemberFunction is applied. As constructed, this predicate apparently only describes what will happen, so changing it amounts to assigning the wrong description, and nothing more.

A corrected predicate is of use though in defining an ImplicitRegion for which RegionMember produces a valid (and fast) RegionMemberFunction.

polyQ=RegionMember@ImplicitRegion[Reduce@RegionMember[polygon,{x,y}],{x,y}]

RegionMemberFunction[…]

The resulting RegionMemberFunction applies the correct predicate.

Moreover, the third part of the RegionMemberFunction describes correctly what the RegionMemberFunction will do.

polyQ[[3]][{x, y}]

(x|y)∈ℝ && ((0≤x<100 && 100≤y≤200) || (100≤x<200 && 100≤y≤300) || (200≤x≤250 && 0≤y≤400))

This can be generalized in a fixed RegionMember function, as follows:

ClearAll@fixedRegionMember; fixedRegionMember[reg_?RegionQ] := RegionMember@ImplicitRegion[Reduce@RegionMember[reg, {x, y}], {x, y}];

Here is the result.

polyQ = fixedRegionMember[polygon]

RegionMemberFunction[…]

With[{nPts = 10000}, testPoints = Transpose@{RandomReal[{0, 300}, nPts], RandomReal[{0, 450}, nPts]}; ListPlot[{testPoints, Select[testPoints, polyQ], polygonPoints}, PlotStyle->{Red, Blue, {PointSize->Large, Green}}, Prolog->{Gray, polygon}, ImageSize->Full] // AbsoluteTiming // Column ]

desired plot, with timing

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.