RegionIntersection puzzle

Question

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?

I can reproduce this, and it looks like a bug to me. A workaround seems to be to change the last digit in 0.9985073695269603of the second polygon to 0.9985073695269604, Apparently Mathematica doesn't like the two polygons to gave a shared point. — Sjoerd C. de Vries, Jun 22 '16 at 19:38

score 13 · Accepted Answer · answered Jun 22 '16 at 19:38

13

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`}}] *)

answered Jun 22 '16 at 19:38

ilian

25,474
4
117
186

Thanks! Floating-point computations strike again. – Joseph O'Rourke Jun 22 '16 at 19:56
1

@ilian - I just noticed that if you cyclicly permute the data points then it will work using the usual precision. Very odd.... – Jason B. Jun 22 '16 at 20:18
@ilian - is there any difference between your sp function and using Rationalize[#, 0]&? – Jason B. Jun 25 '16 at 02:32
2

@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). – ilian Jun 25 '16 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)? – Cassini Aug 04 '16 at 12:35
@Cassini Yes, the (very) soon to be released 11.0. – ilian Aug 04 '16 at 13:13
@ilian Thank you! – Cassini Aug 04 '16 at 17:23

Jason B. · Answer 2 · 2016-06-22T20:17:05.103

5

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.

edited Jun 22 '16 at 20:17

answered Jun 22 '16 at 19:43

Jason B.

68,381
3
139
286

Thanks, Jason! Can you explain why this circumvents the problem? – Joseph O'Rourke Jun 22 '16 at 19:56
1

@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.... – Jason B. Jun 22 '16 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. – Jason B. Jun 22 '16 at 20:03

RegionIntersection puzzle

2 Answers2

Edit - more odd behavior

Linked