TransWikia.com

Speed up the intersection of the diagonals of a regular polygon

Mathematica Asked by expression on January 4, 2021

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 Graphics`Mesh`FindIntersections, but it did not return all the intersection points.

Related link:
How to count the number of n-gons and line intersections in an image of a complete graph?

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[[1]] - #2[[1]]] < 10^-8 &],
     Point /@ List /@ pts]] // AbsoluteTiming

Graphics[{
  MapIndexed[Text[#2[[1]], 1.1 #] &, pts],
  Line /@ diag,
  {Red, intersectionPts}
  }]

enter image description here

4 Answers

This seems to do the trick:

n = 15;
pts = CirclePoints[N@n];
lines = Line /@ Subsets[pts, {2}];
data = Region`Mesh`FindSegmentIntersections[lines];
Graphics[{lines, Red, Point[data[[1]]]}]

Of course, Region`Mesh`FindSegmentIntersections 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.

Correct answer by Henrik Schumacher on January 4, 2021

Complement[#, MeshCoordinates@ConvexHullMesh@#] &@
   MeshCoordinates[DiscretizeRegion[RegionUnion @@ (Line /@ diag)]] //
   Length // AbsoluteTiming
 {0.0676086, 161}
Graphics[{MapIndexed[Text[#2[[1]], 1.1 #] &, pts], Line /@ diag, 
  {Red, PointSize[Large], 
   Point /@ (Complement[#, MeshCoordinates@ConvexHullMesh@#] &@
      MeshCoordinates[DiscretizeRegion[RegionUnion @@ (Line /@ diag)]])}}]

enter image description here

Answered by kglr on January 4, 2021

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 = Region`Mesh`FindSegmentIntersections[lines];
    Graphics[{lines, PointSize[Large], Red, 
      Point[intersectionPts[[1, #]] & /@ intersectionPts[[2, 3, 2, 1]]]}]

Graphics of the plot of a CirclePoint set with Line and intersectionPts

n = 10;
pts = CirclePoints[N@n];
lines = Line /@ Subsets[pts, {2}];
intersectionPts = 
   Region`Mesh`FindSegmentIntersections[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

Region`Mesh`FindSegmentIntersections[
 Line /@ Subsets[CirclePoints[5], {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.

Answered by Steffen Jaeschke on January 4, 2021

Clear["`*"];
n = 20;
pts = CirclePoints[N@n];
pack = Developer`ToPackedArray;
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 enter image description here

Faster but more complex code, when n=100 it takes about 2sec.

Clear[cf];
cf = Partition[Compile`GetElement[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, 
         Internal`StuffBag[bag, (det[a, b] cd - det[c, d] ab)/den, 1]]]] /. expr_ :>
    Compile[{{T, _Real, 3}},
     Block[{bag = Internal`Bag[]}, 
      Do[expr, {i, Length@T}, {j, i + 1, Length@T}]; 
      Internal`BagPart[bag, All]~Partition~2], 
     CompilationTarget -> "C", RuntimeOptions -> "Speed"];

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

Answered by chyanog on January 4, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP