38
$\begingroup$

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

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
$\endgroup$
6
  • 3
    $\begingroup$ With polygonPoints = N@polygonPoints works well. $\endgroup$
    – Kuba
    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$ 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$
    – ilian
    Dec 19, 2015 at 18:12
  • 3
    $\begingroup$ @ilian What a horrible bug. Has this issue been reported to WRI? $\endgroup$
    – Ali Hashmi
    Jul 10, 2019 at 3:02
  • $\begingroup$ @AliHashmi Ok I just reported it. $\endgroup$
    – M.R.
    Sep 8, 2019 at 2:13

1 Answer 1

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.