10
$\begingroup$

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

BugRegInt
RegionIntersection[] returns an incorrect result (the green region above; detail below), and issues a series of error messages:
Errors
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:

PentQuad

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?

Improve this question
$\endgroup$
2
  • $\begingroup$ (In case it matters: I'm running 10.4.0.0 under MacOS.) $\endgroup$
    – Joseph O'Rourke
    Jun 22, 2016 at 19:27
  • 3
    $\begingroup$ 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. $\endgroup$
    – Sjoerd C. de Vries
    Jun 22, 2016 at 19:38

2 Answers 2

Reset to default
13
$\begingroup$

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`}}] *)
Improve this answer
$\endgroup$
7
  • $\begingroup$ Thanks! Floating-point computations strike again. $\endgroup$
    – Joseph O'Rourke
    Jun 22, 2016 at 19:56
  • 1
    $\begingroup$ @ilian - I just noticed that if you cyclicly permute the data points then it will work using the usual precision. Very odd.... $\endgroup$
    – Jason B.
    Jun 22, 2016 at 20:18
  • $\begingroup$ @ilian - is there any difference between your sp function and using Rationalize[#, 0]&? $\endgroup$
    – Jason B.
    Jun 25, 2016 at 2:32
  • 2
    $\begingroup$ @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). $\endgroup$
    – ilian
    Jun 25, 2016 at 19:55
  • $\begingroup$ @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)? $\endgroup$
    – Cassini
    Aug 4, 2016 at 12:35
5
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thanks, Jason! Can you explain why this circumvents the problem? $\endgroup$
    – Joseph O'Rourke
    Jun 22, 2016 at 19:56
  • 1
    $\begingroup$ @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.... $\endgroup$
    – Jason B.
    Jun 22, 2016 at 20:02
  • $\begingroup$ @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. $\endgroup$
    – Jason B.
    Jun 22, 2016 at 20:03

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.