Intersecting graphics

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]}] I know I can use Opacity:

Graphics[{Opacity[0.8], Red, Rectangle[], Green, Disk[{0.2, 0}, .5]}] 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...!): I think it might be easier with raster images, but am interested for now in vector graphics.

• Updated my answer to avoid accusations of "just posting links" ;-) – Szabolcs Jan 23 '12 at 11:52
• This will most likely come with M9: youtube.com/… – P. Fonseca Jan 23 '12 at 20:34
• Nice video - I like his work. Hope the upgrade is cheap... – cormullion Jan 23 '12 at 21:40

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}
] 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}
] • RegionPlot is cool, thanks. If you know the equations for the graphics, it's ideal. – cormullion Jan 23 '12 at 21:41

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[] -
e1[]).Inverse[{Subtract @@ e1, -(Subtract @@ e2)}]);
If[And @@ Thread[0 <= params <= 1],
Subtract @@ e1 params[] + e1[],
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}] • If this really works as indicated it needs to hit "Good Answer" status. – Mr.Wizard Feb 21 '12 at 13:27
• The ReplaceRepeated part looks potentially slow. Have you tested this with largish polygons yet? – Mr.Wizard Feb 21 '12 at 13:30
• @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 Feb 21 '12 at 13:32
• 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 Feb 21 '12 at 13:38
• @Mr.Wizard, Heike I think the major drawback is the (lack of) handling of disconnected polygons in the form of Polygon[{{{1, Sqrt}, {1/2, Sqrt/2}, {0, Sqrt}}, {{-(1/2), Sqrt/ 2}, {-1, Sqrt}, {0, Sqrt}}}]. 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 Mar 28 '12 at 18:11

The (undocumented!) function GraphicsPolygonUtilsPolygonIntersection[] (GraphicsMeshPolygonIntersection[] 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, GraphicsPolygonUtilsPolygonIntersection[polys]}]] Disk[] objects are not covered by this method, but it is not too difficult to make a Polygon[] that approximates a Disk[]...

• Sooner or later we will use more undocumented functions than documented ones... Very useful, +1. – István Zachar Oct 27 '12 at 10:00
• Note that the intersection here is returned as multiple triangles. A single polygon can be obtained by using the option AllTriangles->False – Simon Woods Oct 28 '12 at 19:26
• ...or use PolygonCombine[] on the output of PolygonIntersection[]. – J. M. is away Oct 29 '12 at 11:46
• Please note that with my Mathematica 10.1 version I had to use GraphicsPolygonUtilsPolygonIntersection and that the function just works for 2D polygons (not for 3D ones). – Rainer May 26 '15 at 12:16

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:

<< ImtekPolygon

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]}]}] • +1 for the IMTEK Mathematica Supplement - I could not resist ;-) – user21 Jan 23 '12 at 11:13
• To use the convex intersection algorithm from the IMS, you'd need to appriximate version of the circle. – user21 Jan 23 '12 at 11:14
• @ruebenko why isn't that in your profile? – Mr.Wizard Jan 23 '12 at 11:21
• @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. – user21 Jan 23 '12 at 11:47
• @ruebenko only if I get paid. :^) – Mr.Wizard Jan 23 '12 at 12:00

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: (Note that in this case your Graphics get rasterized and the result is an Image.)

• Thanks - it's definitely easier with images sometimes. – cormullion Jan 23 '12 at 21:42

You may be able to do this using FilledCurve in version 8. I see examples of subtraction and exclusion but not intersection. 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.