Intersecting graphics

49

25

Does the Mathematica graphics system have any concept of intersecting graphics? I've not found much in the documents so far. For example, if I want to show the intersection of two shapes:

Graphics[{Rectangle[], Disk[{0.2, 0}, .5]}]

graphics without intersection

I know I can use Opacity:

Graphics[{Opacity[0.8], Red, Rectangle[], Green, Disk[{0.2, 0}, .5]}]

graphic with intersection

But is there a way of specifying the colours of intersecting areas directly? It doesn't seem to be possible to 'address' the intersecting shapes any other way.

In the same vein, is is possible to 'extract' the graphical intersection of arbitrary shapes, without returning to the original geometry and calculating it? Could you obtain this type of entity easily given the above specification (these are just examples...!):

the intersection

I think it might be easier with raster images, but am interested for now in vector graphics.

cormullion

Posted 2012-01-23T10:41:33.533

Reputation: 23 565

Updated my answer to avoid accusations of "just posting links" ;-) – Szabolcs – 2012-01-23T11:52:27.750

1

This will most likely come with M9: https://www.youtube.com/watch?v=CCTa904eRfM&list=PLCD1C4A44DFA4D7C6&index=11&feature=plpp_video

– P. Fonseca – 2012-01-23T20:34:39.100

Nice video - I like his work. Hope the upgrade is cheap... – cormullion – 2012-01-23T21:40:32.393

Answers

23

How about RegionPlot?

RegionPlot[
  {
   (x - 0.2)^2 + y^2 < 0.5 && 0 < x < 1 && 0 < y < 1,
   (x - 0.2)^2 + y^2 < 0.5 && ! (0 < x < 1 && 0 < y < 1),
   ! ((x - 0.2)^2 + y^2 < 0.5) && 0 < x < 1 && 0 < y < 1
  }, 
   {x, -1, 1.5}, {y, -1, 1.5}, 
   PlotStyle -> {Red, Yellow, Blue}
]

Mathematica graphics

EDIT in response to Szabolcs's comment:

PointInPoly[{x_, y_}, poly_List] := 
 Module[{i, j, c = False, npol = Length[poly]}, 
  For[i = 1; j = npol, i <= npol, j = i++, 
   If[((((poly[[i, 2]] <= y) && (y < 
             poly[[j, 2]])) || ((poly[[j, 2]] <= y) && (y < 
             poly[[i, 2]]))) && (x < (poly[[j, 1]] - 
             poly[[i, 1]])*(y - poly[[i, 2]])/(poly[[j, 2]] - 
              poly[[i, 2]]) + poly[[i, 1]])), c = ¬ c];];
  c]

(from an answer I gave in MathGroup)

RegionPlot[{
   PointInPoly[{x, y}, {{1, 3}, {3, 4}, {4, 7}, {5, -1}, {3, -3}}] && 
   PointInPoly[{x, y}, {{2, 2}, {3, 3}, {4, 2}, {0, 0}}], 
   PointInPoly[{x, y}, {{1, 3}, {3, 4}, {4, 7}, {5, -1}, {3, -3}}] &&
   ¬ PointInPoly[{x, y}, {{2, 2}, {3, 3}, {4, 2}, {0, 0}}],
   ¬ PointInPoly[{x, y}, {{1, 3}, {3, 4}, {4, 7}, {5, -1}, {3, -3}}] &&
   PointInPoly[{x, y}, {{2, 2}, {3, 3}, {4, 2}, {0, 0}}]}, 
  {x, 0, 6}, {y, -4, 8}, 
  PlotPoints -> 100, PlotStyle -> {Red, Yellow, Blue}
]

Mathematica graphics

Sjoerd C. de Vries

Posted 2012-01-23T10:41:33.533

Reputation: 63 549

RegionPlot is cool, thanks. If you know the equations for the graphics, it's ideal. – cormullion – 2012-01-23T21:41:33.547

35

I'm coming to the party a bit late, but here's my approach. It should work for any two polygons, including non-convex and self-intersecting ones.

winding[poly_, pt_] := 
 Round[(Total@
      Mod[(# - RotateRight[#]) &@(ArcTan @@ (pt - #) & /@ poly), 
       2 Pi, -Pi]/2/Pi)]
cross[e1_, e2_] /; (N[Det[{Subtract @@ e1, Subtract @@ e2}]] === 0.) =
   None;
cross[e1_, e2_] := Module[{params},
   params = ((e2[[2]] - 
        e1[[2]]).Inverse[{Subtract @@ e1, -(Subtract @@ e2)}]);
   If[And @@ Thread[0 <= params <= 1], 
    Subtract @@ e1 params[[1]] + e1[[2]],
    None]];

intersection[poly1_, poly2_, p : {in1_, in2_} : {1, 1}] := 
  Module[{edges1, edges2, intersections,
    inter1, inter2, newedges1, newedges2, midpoints1, midpoints2},
   edges1 = Partition[Range[Length[poly1]], 2, 1, {1, 1}];
   edges2 = Partition[Range[Length[poly2]], 2, 1, {1, 1}];

   intersections = Table[cross[poly1[[e1]], poly2[[e2]]],
     {e1, edges1}, {e2, edges2}];
   inter1 = Flatten[Table[
      SortBy[
       Prepend[DeleteCases[intersections[[i]], None], poly1[[i]]], 
       Norm[# - poly1[[i]]] &], {i, Length[edges1]}], 1];
   inter2 = 
    Flatten[Table[
      SortBy[Prepend[DeleteCases[intersections[[All, i]], None], 
        poly2[[i]]], Norm[# - poly2[[i]]] &], {i, Length[edges2]}], 1];

   newedges1 = Partition[inter1, 2, 1, {1, 1}];
   newedges2 = Partition[inter2, 2, 1, {1, 1}];

   midpoints1 = Mean /@ newedges1;
   midpoints2 = Mean /@ newedges2;
   Flatten[{Pick[newedges1, Abs[winding[poly2, #]] & /@ midpoints1, 
       in1],
      Pick[newedges2, Abs[winding[poly1, #]] & /@ midpoints2, in2]}, 
     1] //.
    {{a___, {b__, c_List}, d___, {c_, e__}, 
       f___} :> {a, {b, c, e}, d, f},
     {a___, {b__, c_List}, d___, {e__, c_}, f___} :> {a, 
       Join[{b, c}, Reverse[{e}]], d, f},
     {a___, {c_List, b__}, d___, {c_, e__}, f___} :> {a, 
       Join[Reverse[{e}], {c, b}], d, f},
     {a___, {c_List, b__}, d___, {e__, c_}, f___} :> {a, {e, c, b}, d,
        f}
     }
   ];

Some notes

winding and cross are two helper functions. winding calculates the winding number of a point pt with respect to a polygon poly given as a list of vertex coordinates. A point lies inside a polygon if and only if the winding number is non-zero.

The function cross calculates the intersection point of two line segments, or returns None if they don't intersect.

intersection is the main function which calculates the intersecting polygon of two polygons poly1 and poly2. It works by calculating the intersection points between the two polygons and adding these to the vertex lists of poly1 and poly2. Then each of the edges of the new polygons lie either completely inside or outside of the other polygon.

The intersection of the two polygons $\text{poly1} \cap \text{poly2}$ is then the union of edges of poly1 that lie inside poly2 and vice versa. Similarly one can also calculate the complement of the two polygons, $\text{poly1} \backslash \text{poly2}$ and $\text{poly1} \backslash \text{poly2}$, and the union $\text{poly1} \cup \text{poly2}$. These four options can be set by in1 and in2.

Example

Manipulate[DynamicModule[{ips11, ips10, ips01},
  pts = PadRight[pts, 2 n, RandomReal[{-1, 1}, {2 n, 2}]];
  ips11 = intersection[pts[[ ;; n]], pts[[n + 1 ;;]], {1, 1}];
  ips10 = intersection[pts[[ ;; n]], pts[[n + 1 ;;]], {1, 0}];
  ips01 = intersection[pts[[ ;; n]], pts[[n + 1 ;;]], {0, 1}];
  Graphics[{
    {Yellow, Polygon[ips10]},
    {Blue, Polygon[ips01]},
    {Red, Polygon[ips11]},
    {FaceForm[], EdgeForm[Black], Polygon[pts[[ ;; n]]]}, {FaceForm[], 
     EdgeForm[Black], Polygon[pts[[n + 1 ;;]]]}}, 
   PlotRange -> {{-1, 1}, {-1, 1}}]], 
  {{pts, {}}, Locator}, {{n, 5}, None}]

Mathematica graphics

Heike

Posted 2012-01-23T10:41:33.533

Reputation: 34 748

The winding number routine can be improved slightly, by making use of ListCorrelate[]: windingNumber[poly_?MatrixQ, pt_?VectorQ] := Round[Total @ Mod[ListCorrelate[{-1, 1}, ArcTan @@ (pt - #) & /@ poly, {-1, -1}], 2 Pi, -Pi]/2/Pi] – J. M.'s ennui – 2012-05-25T13:17:28.580

1If this really works as indicated it needs to hit "Good Answer" status. – Mr.Wizard – 2012-02-21T13:27:53.707

The ReplaceRepeated part looks potentially slow. Have you tested this with largish polygons yet? – Mr.Wizard – 2012-02-21T13:30:32.797

@Mr.Wizard I haven't tried larger polygons yet. I agree about the ReplaceRepeated but it's the best I could come up with for now. If you know about a better way to join a set of edges I would be interested. – Heike – 2012-02-21T13:32:59.520

Honestly I cannot even tell what your code is doing yet. It must have been a lot of work putting it together. Some time later I shall work through it and see what possible improvements come to mind. If I don't get around to it in the next couple of days please remind me. – Mr.Wizard – 2012-02-21T13:38:18.123

If I am right, there is a bug with the winding function: in the Map it should be If[pt == #, 0, ArcTan @@ (pt - #)] as ArcTan[0, 0] is indeterminate. See this example of two intersecting triangles: intersection[{{1, 0}, {0, Sqrt[3]}, {-1, 0}}, {{1, Sqrt[3]}, {0, 0}, {-1, Sqrt[3]}}, {1, 1}]. Could you please confirm it? – István Zachar – 2012-03-28T16:26:58.007

1@Mr.Wizard, Heike I think the major drawback is the (lack of) handling of disconnected polygons in the form of Polygon[{{{1, Sqrt[3]}, {1/2, Sqrt[3]/2}, {0, Sqrt[3]}}, {{-(1/2), Sqrt[3]/ 2}, {-1, Sqrt[3]}, {0, Sqrt[3]}}}]. This quickly leads to a combinatorial explosion when one has to check each subpart in poly1 with every other in poly2, if there is no pretesting for whether two polys are touching or not. Also, it leaves some redundant coordinates in the result like {{0,0}, {0,0}, {0,0}}. Of course, this might be useful for the user, but there should be some means to remove them. – István Zachar – 2012-03-28T18:11:06.827

17

The (undocumented!) function Graphics`PolygonUtils`PolygonIntersection[] (Graphics`Mesh`PolygonIntersection[] in older versions) seems up to the task. Using Sjoerd's example:

polys = {Polygon[{{1, 3}, {3, 4}, {4, 7}, {5, -1}, {3, -3}}],
         Polygon[{{2, 2}, {3, 3}, {4, 2}, {0, 0}}]};

Graphics[Append[{Gray, polys}, {Blue, Graphics`PolygonUtils`PolygonIntersection[polys]}]]

intersecting polygons

Disk[] objects are not covered by this method, but it is not too difficult to make a Polygon[] that approximates a Disk[]...

J. M.'s ennui

Posted 2012-01-23T10:41:33.533

Reputation: 115 520

1Please note that with my Mathematica 10.1 version I had to use Graphics`PolygonUtils`PolygonIntersection and that the function just works for 2D polygons (not for 3D ones). – Rainer – 2015-05-26T12:16:20.143

2Sooner or later we will use more undocumented functions than documented ones... Very useful, +1. – István Zachar – 2012-10-27T10:00:53.047

3Note that the intersection here is returned as multiple triangles. A single polygon can be obtained by using the option AllTriangles->False – Simon Woods – 2012-10-28T19:26:03.607

2...or use PolygonCombine[] on the output of PolygonIntersection[]. – J. M.'s ennui – 2012-10-29T11:46:44.527

14

I am not aware of any built-in functionality (I might easily be wrong), but there's an example at MathWorld for calculating intersections of convex polygons. You'd need to approximate the circle with a polygon.

Get the notebook from that page: there's an intersection calculation inside that uses the IMTEK Mathematica supplement.

Example:

<< Imtek`Polygon`

disk = Disk[{0.2, 0}, 0.5];
rec = Polygon[{{0, 0}, {1, 0}, {1, 1}, {0, 1}}];

Graphics[{Green, rec, Red, disk, Blue, 
  Polygon@imsConvexIntersect[{imsPolygonizeCircle[Circle @@ disk, 50],
      First[rec]}]}]

Mathematica graphics

Szabolcs

Posted 2012-01-23T10:41:33.533

Reputation: 213 047

+1 for the IMTEK Mathematica Supplement - I could not resist ;-) – None – 2012-01-23T11:13:28.873

To use the convex intersection algorithm from the IMS, you'd need to appriximate version of the circle. – None – 2012-01-23T11:14:43.410

1@ruebenko why isn't that in your profile? – Mr.Wizard – 2012-01-23T11:21:54.577

@Mr.Wizard, well I don't really have time to update and maintain this; would you like to have a go with that ;-) The IMS needs a new place to live, and should be renames to International Mathematica Supplement (them the ims prefix for the functions can stay). Suggestions are welcome. – None – 2012-01-23T11:47:16.353

@ruebenko only if I get paid. :^) – Mr.Wizard – 2012-01-23T12:00:34.920

@Mr.Wizard, you will get payed in honor, and real reputation ;-) – None – 2012-01-23T12:03:01.303

@ruebenko why not put it in the Wolfram Library Archive? This place seems to have been all but forgotten (by both users and WRI) but seems like the place for a repository. – Mike Honeychurch – 2012-01-23T20:36:37.660

@MikeHoneychurch, that's a possibility. I have not made up my mind. Thanks for the suggestion. – None – 2012-01-23T21:00:19.213

Excellent suggestion, thanks. – cormullion – 2012-01-23T21:43:45.897

10

Another option is to use image processing features such as ImageCompose:

{
 g1 = Graphics[Rectangle[], PlotRange -> 1],
 g2 = Graphics[Disk[{0.2, 0}, .5], PlotRange -> 1],
 ImageCompose[g1, g2, Center, Center, {1, 1, 0}]
 }

The output of the above looks like this:

enter image description here

(Note that in this case your Graphics get rasterized and the result is an Image.)

Andrew Moylan

Posted 2012-01-23T10:41:33.533

Reputation: 4 130

Thanks - it's definitely easier with images sometimes. – cormullion – 2012-01-23T21:42:40.923

8

In version 10, the new geometric computation functionality supports this. It operates with region objects. Many graphics primitives, including Disk and Rectangle can be used as regions.

Boolean operations include RegionUnion, RegionIntersect, RegionDifference, RegionSymmetricDifference and BooleanRegion.

Example:

RegionPlot[RegionIntersection[Rectangle[], Disk[{0.2, 0}, .5]], AspectRatio->Automatic]

In many cases these Boolean operators do not evaluate. RegionIntersection[Rectangle[], Disk[{0.2, 0}, .5]] can be used as a region, but it does not evaluate to another expression. However, it is possible to get an approximation to the result, free of any explicit Boolean expressions, using DiscretizeRegion. It is then possible to extract coordinates (e.g. the boundary) from this discretized mesh.

Szabolcs

Posted 2012-01-23T10:41:33.533

Reputation: 213 047

7

You may be able to do this using FilledCurve in version 8. I see examples of subtraction and exclusion but not intersection.

enter image description here

Mr.Wizard

Posted 2012-01-23T10:41:33.533

Reputation: 259 163