# Strategy to speed up two polygons (or region) intersection

Basically, I want to compare the overlap between two polygons (or regions).

First, I have a reference particle ring that can form a simple polygon, but this type of ring is very thin and long, and its area is very small:

Referencering = Polygon[dataring[[72]]]
Area[Referencering]
0.0000593772


My numerical polygon comes from extracting data from ContourPlot, it has a main region surrounded by several small independent regions or may contain multiple holes within the main region. I used an undocumented function GraphicsPolygonUtilsPolygonCombine to combine it into an entire polygon.

NumerialPolygon = GraphicsPolygonUtilsPolygonCombine[Polygon /@ Contourpoints]


I used RegionIntersection to measure the overlap between these two polygons. I tried three approaches:

The first one is directly calculated from these two polygons, and it takes about 5.4 hours to get the result.

IntersectionRegion =
RegionIntersection[NumerialPolygon,Referencering]; // AbsoluteTiming

{19372.7, Null}


The second one I tried is to use BoundaryDiscretizeRegion to convert my polygon to BoundaryMeshRegion and then use RegionIntersection to measure the overlap. This approach is faster than the previous one and takes 1.5 hours.

NumericalRegion = BoundaryDiscretizeRegion[NumerialPolygon]; // AbsoluteTiming
{5867.37, Null}
IntersectionRegion = RegionIntersection[NumericalRegion,ReferenceRegion]; // AbsoluteTiming
{37.2518, Null}


The third one I tried is to use BoundaryDiscretizeGraphics to get BoundaryMeshRegion from ContourPlot. But this method is also quite slow, taking about 3.8 hours.

From my current observation, it seems RegionIntersection will take more time to deal with a polygon than BoundaryMeshRegion. I think it is just for this type of polygon; other simple polygons I tested all run very quickly. So, I would love to hear if anyone can have a better approach to deal with this kind of issue.

Many thanks!

• Please post dataring. Commented May 9 at 21:59
• Hello @cvgmt , You can download this dataring file from github.com/CerberusJw/ShareData. Thanks!
– J. W
Commented May 10 at 2:59
• Please always post al relevant data. We also need Contourpoints I guess. Commented May 11 at 12:36
• "I used an undocumented function GraphicsPolygonUtilsPolygonCombine to combine it into an entire polygon." That is most likely counterproductive. This prevents the use of divide-and-conquer algorithms. Same goes for dataring. What you actually want is a triangulation of the polygon, e.g., the one you get from TriangulateMesh[Polygon[dataring]]. Commented May 11 at 12:41
• Also, you say that you generate Contourpoints from a contour plot. Wouldn't it be easier to just evaluate the contour function on the vertices of dataring to get a good idea which partd lie inside a sublevel set or outside? Commented May 11 at 12:48

As with many MeshRegion related things, Mathematica is extremely stubborn. Somebody must have told it recently to convert planar domains to polygons all the time. IMHO, that is an incredibly bad idea for many tasks in computational geometry. There we want triangulations whenever we can.

Here is what I do: I load the data points and create MeshRegions for them. RegionUnion is awfully slow, so I unite the region created from Contourpoints manually. Then I do one call to RegionIntersection.

Contourpoints =
dataring =

First@AbsoluteTiming[
P = TriangulateMesh[Polygon[#]] & /@ Contourpoints;
pts = MeshCoordinates /@ P;
triangles = (MeshCells[#, 2, "Multicells" -> True] & /@ P)[[All, 1,
1]];
vcounts = Prepend[Accumulate[Most[Length /@ pts]], 0];
allpoints = Join @@ pts;
alltriangles = Join @@ (triangles + vcounts);
R = MeshRegion[allpoints, Triangle[alltriangles]];
]

First@AbsoluteTiming[
Q = TriangulateMesh[Polygon[dataring]];
]

First@AbsoluteTiming[
inter = RegionIntersection[Q, R];
]


6.94497

0.589849

29.402

IMHO the whole thing is still way to slow. I think a major reason for that is that RegionIntersection wants to return a polygon again. Why?!? IIRC, it used to return a MeshRegion containing a triangulation!

And performance of these things is extremely brittle, too. One cannot really rely on what is fast and what is slow; it seems that this changes in every minor release...

• Thanks @HenrikSchumacher ! I totally agree with you! I think the hard part is to combine all Contourpoints together. Contourpoints unavoidably contain many points, which burden us when we try to form a triangulation of this polygon. All calculations aside, I found GraphicsPolygonUtilsPolygonCombine did an excellent job of combining those polygons. But just stay at this level. The other shortcut is to only take the main part from the Contourpoints (Contourpoints[[1]]), then use RegionIntersection with dataring. However, it will lose some accuracy among this numerical contour.
– J. W
Commented May 12 at 2:51