# RegionIntersection puzzle

I am using RegionIntersection[] to intersect two rectangles, phrased as Polygon[]s:

RegionIntersection[] returns an incorrect result (the green region above; detail below), and issues a series of error messages:

Here is the region returned as the intersection, which is not correct (it should be a quadrilateral):

Polygon[{{0.407273, 0.650444}, {0.509656, -0.0200315}, {0.998507,
0.0546171}, {0.640904, 0.767621}, {0.998507, 0.0546171}}]


Here it is as a line drawing:

Here is the originating call that produces the incorrect intersection:

RegionIntersection[
Polygon[{{0.5096555454081809, -0.02003146973392257}, \
{0.9985073695269602, 0.05461714932464575}, {0.6966031018052412,
2.0316992956026936}, {0.2077512776864619,
1.9570506765441251}}] ,
Polygon[{{0.9985073695269602,
0.05461714932464575}, {0.6409040075686694,
0.767621034809768}, {-1.1468443539373534, \
-0.12901478643123643}, {-0.7892409919790626, -0.8420186719163586}}]
]


Can anyone see what's going on?

• (In case it matters: I'm running 10.4.0.0 under MacOS.) Jun 22, 2016 at 19:27
• I can reproduce this, and it looks like a bug to me. A workaround seems to be to change the last digit in 0.9985073695269603 of the second polygon to 0.9985073695269604, Apparently Mathematica doesn't like the two polygons to gave a shared point. Jun 22, 2016 at 19:38

This is a bug that has been fixed in the development version. For a possible workaround, use exact coordinates, for example

sp = Function[p, SetPrecision[p, Infinity]];

ri = RegionIntersection[
sp@Polygon[{{0.5096555454081809, -0.02003146973392257},
{0.9985073695269602, 0.05461714932464575},
{0.6966031018052412, 2.0316992956026936},
{0.2077512776864619, 1.9570506765441251}}],
sp@Polygon[{{0.9985073695269602, 0.05461714932464575},
{0.6409040075686694, 0.767621034809768},
{-1.1468443539373534, -0.12901478643123643},
{-0.7892409919790626, -0.8420186719163586}}]];

N[ri]

(* Polygon[{{0.40727258068338046, 0.6504444171024929},
{0.5096555454081809, -0.02003146973392257},
{0.9985073695269602, 0.05461714932464575},
{0.6409040075686694, 0.767621034809768}}] *)

• Thanks! Floating-point computations strike again. Jun 22, 2016 at 19:56
• @ilian - I just noticed that if you cyclicly permute the data points then it will work using the usual precision. Very odd.... Jun 22, 2016 at 20:18
• @ilian - is there any difference between your sp function and using Rationalize[#, 0]&? Jun 25, 2016 at 2:32
• @JasonB Yes, these are different. Rationalize does a continued fraction-based approximation, while SetPrecision works directly off the internal bit representation. There should be an example in the documentation (under Properties & Relations). Jun 25, 2016 at 19:55
• @ilian Will this fix be included in the next minor update or do we have to wait until v11 (or is that the next update)? Aug 4, 2016 at 12:35

Another workaround is to turn the polygons into MeshRegions first,

RegionIntersection @@ (DiscretizeGraphics /@ {p1, p2}) //
MeshPrimitives[#, 2] & // First
(* Polygon[{{0.407273,
0.650444}, {0.509656, -0.0200315}, {0.998507,
0.0546171}, {0.640904, 0.767621}}] *)


Where p1 and p2 are your polygons

## Edit - more odd behavior

What's even odder is that if you change the order of the points of the polygon (cyclicly so that you don't change the shape), then it will work. If points1 and points2 are the points of the polygon in the OP, then

RegionIntersection @@ (Polygon /@ {points1, points2})


fails as it did for OP, but

RegionIntersection @@ (Polygon /@ {RotateLeft@points1,
RotateLeft@points2})


works perfectly.

• Thanks, Jason! Can you explain why this circumvents the problem? Jun 22, 2016 at 19:56
• @JosephO'Rourke - I wish I could! I just have had some experience with regions lately, and one of the basic steps when trying anything is to try it on a MeshRegion versus a regular region. For it to work on a Polygon or other built in shape (like Circle or Disk), I think the R&D guys have to code the cases in individually. But when using a MeshRegion there is a general toolkit that is used. @ilian would know more about that than I do though.... Jun 22, 2016 at 20:02
• @JosephO'Rourke - so if you are working in 2D with a non-infinite region, then turning it into a MeshRegion with DiscretizeGraphics is often a great way to go. But if you are working in 3D, then it won't work. Jun 22, 2016 at 20:03