# Speed up the intersection of the diagonals of a regular polygon

I need to calculate all the intersection points of the diagonals of a regular polygon, following code is really slow, when n=15, it take about 30sec. I also tried using GraphicsMeshFindIntersections, but it did not return all the intersection points.

Diagonals of a regular octagon

Clear["*"];
n = 10;
pts = CirclePoints[N@n];

diag = Join @@ Table[pts[[{i, j}]], {i, n - 2}, {j, i + 2, n - Boole[i < 2]}];

Length[intersectionPts = Complement[DeleteDuplicates[DeleteCases[RegionIntersection /@
Subsets[Line /@ diag, {2}], _EmptyRegion], Norm[#1[] - #2[]] < 10^-8 &],
Point /@ List /@ pts]] // AbsoluteTiming

Graphics[{
MapIndexed[Text[#2[], 1.1 #] &, pts],
Line /@ diag,
{Red, intersectionPts}
}] This seems to do the trick:

n = 15;
pts = CirclePoints[N@n];
lines = Line /@ Subsets[pts, {2}];
data = RegionMeshFindSegmentIntersections[lines];
Graphics[{lines, Red, Point[data[]]}]


Of course, RegionMeshFindSegmentIntersections is not documented... I found it by spelunking with ?**Intersect*. You may inspect Rest[data] to find out about the classification of these intersections. For example, you probably what to filter out the hits found under EndPointsTouching as these are false positives. I guess that

data[[1, data[[2, 3, 2, 1]]]]


should be what you are looking for.

Complement[#, MeshCoordinates@ConvexHullMesh@#] &@
MeshCoordinates[DiscretizeRegion[RegionUnion @@ (Line /@ diag)]] //
Length // AbsoluteTiming

 {0.0676086, 161}

Graphics[{MapIndexed[Text[#2[], 1.1 #] &, pts], Line /@ diag,
{Red, PointSize[Large],
Point /@ (Complement[#, MeshCoordinates@ConvexHullMesh@#] &@
MeshCoordinates[DiscretizeRegion[RegionUnion @@ (Line /@ diag)]])}}] • In general {Graphics[{Red, Line[\[CapitalOmega][]]}], Graphics[{Blue, Line@Append[\[CapitalOmega][], First@\[CapitalOmega][]]}]} is different for some line set in 2-dimensions. – Steffen Jaeschke Sep 6 at 14:26

A diagonal is a geometric construct consisting of points is a line connecting two points that are not in the direct neighborhood.

So the possible fastest solution is:

    n = 10;
pts = CirclePoints[N@n];
lines = Line /@ Subsets[pts, {2}];
intersectionPts = RegionMeshFindSegmentIntersections[lines];
Graphics[{lines, PointSize[Large], Red,
Point[intersectionPts[[1, #]] & /@ intersectionPts[[2, 3, 2, 1]]]}] n = 10;
pts = CirclePoints[N@n];
lines = Line /@ Subsets[pts, {2}];
intersectionPts =
RegionMeshFindSegmentIntersections[lines]; // AbsoluteTiming


{0.000813, Null}


RegionMeshFindSegmentIntersections uses the sweep line algorithm as found in the literature.

Intersection algorithms

"There are a number of problems involving the computation of intersections between geometric objects. Perhaps the most interesting of these for an algorithms course is a sweep-line algorithm to find all intersections between n line segments. The algorithms works in time O((n+I)log n) and O(n) space, where I is the number of intersections reported [BS79,PS91]. (See also [BKOS97,Ch.2].)

The basic idea of a sweep-line algorithm is to turn a 2-dimensional static algorithm into a 1-dimensional dynamic one. The approach is to sweep a vertical line across the set of segments from left to right and to keep track of the order that the segments which intersect this sweep line lie along the line. See Figure 2. As the sweep line moves across the segments, intersections appear (when the line reaches a left endpoint of a segment) and disappear (when the line passes beyond the right endpoint of a segment). The segments at given location of the sweep line are ordered from bottom to top. This order changes precisely when segments intersect.

The sweep-line algorithm is a discrete event simulation. The events are: encountering a left endpoint; encountering a right endpoint; and two segments changing order (intersecting). At first glance it seems that one would have to know all of the intersection points in advance to do this simulation, but it turns out that intersections can be computed "on the fly" as previous events are processed. These ideas (changing a 2-D static problem into a 1-D dynamic one, discrete event simulation, sweep-line, and discovering intersections before the sweepline reaches them) are interesting and powerful. This is also a nice data structures problem involving a balanced binary tree to keep track of the order of segments along the sweep line and a priority queue to keep track of future events."

The Undergraduate Algorithms Course and Recent Research in Computational Geometry implemented in Mathematica V 12.0.0..

This is an order in magnitude faster and avoids the generating points with evident no crossing but endpoints. The can be read out of

RegionMeshFindSegmentIntersections[
Line /@ Subsets[CirclePoints, {2}]]

{{{-0.587785, -0.809017}, {0.587785, -0.809017}, {0.,
1.}, {0.363271, -0.118034}, {-1.11022*10^-16, -0.381966}, \
{-0.363271, -0.118034}, {0.951057, 0.309017}, {-0.951057,
0.309017}, {-0.224514, 0.309017}, {0.224514,
0.309017}}, {{"EndPointsTouching",
Point[{8, 7, 3, 2, 1}]}, {"EndPointTouchingSegment",
Point[{}]}, {"SegmentsIntersect",
Point[{10, 9, 6, 5, 4}]}, {"PointTouchesEndPoint",
Point[{}]}, {"PointTouchesSegment", Point[{}]}, {"PointsOverlap",
Point[{}]}, {"SegmentsOverlap", Line[{}]}}}


The point list is {"SegmentsIntersect", Point[{10, 9, 6, 5, 4}]} in this example.

So not only that the RegionMeshFindSegmentIntersections is well in the output explicitly well documented it really the fastest.

I hope that solves the question to the required degree.

• How did you find the RegionMeshFindSegmentIntersections function? – C. E. Sep 7 at 15:00
Clear["*"];
n = 20;
pts = CirclePoints[N@n];
pack = DeveloperToPackedArray;
diag = Join @@ Table[pts[[{i, j}]], {i, n - 2}, {j, i + 2, n - Boole[i < 2]}] // pack;

lineIntersection =
Partition[Indexed[T, #] & /@ Tuples[{1, 2}, 3], 2] /. {a_, b_, c_, d_} :>
Compile[{{T, _Real, 3}},
(Det[{a, b}] (c - d) - Det[{c, d}] (a - b))/(Det[{a - b, c - d}] - 5*^-15) // Evaluate,
RuntimeAttributes -> {Listable}];

Length[intersectionPts =
Subsets[diag, {2}] // pack // lineIntersection //
Pick[#, Unitize[Sqrt[(#^2).{1, 1}], 1], 0] & //
Nearest[#, DeleteDuplicates[Round[#, 10.^-8]], 1][[All, 1]] &] // AbsoluteTiming

Graphics[{
MapIndexed[Text[Tr@#2, 1.05 #] &, pts],
Line /@ diag,
{Red, PointSize[Small], Point@intersectionPts}
}, ImageSize -> Large]


when n=100, it takes about 5sec, the number of intersection points is 3731201.
http://oeis.org/A006561 Faster but more complex code, when n=100 it takes about 2sec.

Clear[cf];
cf = Partition[CompileGetElement[T, ##] & @@@ Tuples[{{i, j}, {1, 2}, {1, 2}}],2] /. {a_, b_, c_, d_} :>
With[{det = Cross[#].#2 &, ab = a - b, bc = b - c, cd = c - d, ac = a - c, ad = a - d},
With[{den = det[ab, cd]},
Evaluate /@ If[Abs[den] > 10.^-8 && det[ab, ac] det[ab, ad] < 0 && det[cd, ac] det[cd, bc] < 0,
InternalStuffBag[bag, (det[a, b] cd - det[c, d] ab)/den, 1]]]] /. expr_ :>
Compile[{{T, _Real, 3}},
Block[{bag = InternalBag[]},
Do[expr, {i, Length@T}, {j, i + 1, Length@T}];
InternalBagPart[bag, All]~Partition~2],
CompilationTarget -> "C", RuntimeOptions -> "Speed"];

Length[intersectionPts2 =
cf@diag // Nearest[#, DeleteDuplicates[Round[#, 10.^-8]], 1][[All, 1]] &] // AbsoluteTiming