Bug in RegionMember?

Question

Bug introduced in 10.0 and persisting through 12.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}]

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}]]]}]

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

Answer 1

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
 ]