Add pairwise intersection between segments with sweep-line algorithm

src/utils.jl

@@ -11,3 +11,4 @@ include("utils/cmp.jl")
 include("utils/units.jl")
 include("utils/crs.jl")
 include("utils/misc.jl")
+include("utils/intersect.jl")

src/utils/intersect.jl

@@ -0,0 +1,91 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+"""
+    pairwiseintersect(segments; [digits])
+
+Compute pairwise intersections between n `segments`
+with `digits` precision in O(n⋅log(n+k)) time
+where k is the number of intersections, using
+a sweep line algorithm. Similar to an optimal
+Bentley-Ottmann algorithm in sparse systems,
+and closer to O(n²) in dense systems.
+
+Return intersection points and corresponding indices
+of segments involved in the intersection.
+
+By default, set `digits` based on the absolute
+tolerance of the length type of the segments.
+
+## Examples
+
+```julia
+points, seginds = pairwiseintersect(segments)
+points[i] # i-th intersection point
+seginds[i] # corresponding segments
+```
+
+## References
+
+* Bentley & Ottmann 1979. [Algorithms for reporting and counting
+  geometric intersections](https://ieeexplore.ieee.org/document/1675432)
+"""
+function pairwiseintersect(segments; digits=_digits(segments))
+  # orient segments and round coordinates
+  segs = map(segments) do s
+    a, b = coordround.(extrema(s), digits=digits)
+    a > b ? Segment(b, a) : Segment(a, b)
+  end
+
+  # extract first (or "x") coordinate from
+  # first and last vertices of segments
+  x(p) = flat(coords(p)).x
+  xₛ = [x(first(vertices(s))) for s in segs]
+  xₑ = [x(last(vertices(s))) for s in segs]
+
+  # sort segments based on first coordinates
+  inds = sortperm(xₛ)
+  segs = segs[inds]
+  xₛ = xₛ[inds]
+  xₑ = xₑ[inds]
+
+  # sweepline algorithm
+  n = length(segs)
+  P = eltype(first(segs))
+  D = Dict{P,Vector{Int}}()
+  for i in 1:n
+    for j in (i + 1):n
+      # break if segments don't overlap w.r.t. first coordinate
+      xₛ[i] ≤ xₛ[j] ≤ xₑ[i] || break
+
+      # perform more expensive intersection algorithm
+      intersection(segs[i], segs[j]) do I
+        if type(I) == Crossing || type(I) == EdgeTouching
+          p = coordround(get(I); digits)
+          if haskey(D, p)
+            append!(D[p], (inds[i], inds[j]))
+          else
+            D[p] = [inds[i], inds[j]]
+          end
+        end
+      end
+    end
+  end
+
+  # remove duplicate indices
+  for p in keys(D)
+    unique!(D[p])
+  end
+
+  collect(keys(D)), collect(values(D))
+end
+
+# compute the number of significant digits based on the segment type
+# this is used to determine the precision of the points
+function _digits(segments)
+  s = first(segments)
+  ℒ = lentype(s)
+  τ = ustrip(eps(ℒ))
+  round(Int, 0.8 * (-log10(τ))) # 0.8 is a heuristic to avoid numerical issues
+end

test/utils.jl

@@ -80,6 +80,100 @@ end
   @inferred Meshes.coordround(p₁, digits=10)
 end
 
+@testitem "pairwiseintersect" setup = [Setup] begin
+  # helper to sort points and seginds by point coordinates
+  function sortedintersection(segs)
+    points, inds = Meshes.pairwiseintersect(segs)
+    perm = sortperm(points)
+    points[perm], inds[perm]
+  end
+
+  # simple endpoint case
+  segs = Segment.([(cart(0, 0), cart(2, 2)), (cart(0, 2), cart(2, 0)), (cart(0, 1), cart(0.5, 1))])
+  points, seginds = sortedintersection(segs)
+  @test length(points) == 1
+  @test length(seginds) == 1
+
+  # small number of segments, handling endpoints and precision
+  segs = Segment.([(cart(0, 0), cart(2, 2)), (cart(1.5, 1), cart(2, 1)), (cart(1.51, 1.3), cart(2, 0.9))])
+  points, seginds = sortedintersection(segs)
+  @test length(points) == 1
+
+  # box case with one segment outside
+  segs =
+    Segment.([
+      (cart(0, 0), cart(1.1, 1.1)),
+      (cart(1, 0), cart(0, 1)),
+      (cart(0, 0), cart(0, 1)),
+      (cart(0, 0), cart(1, 0)),
+      (cart(0, 1), cart(1, 1)),
+      (cart(1, 0), cart(1, 1))
+    ])
+  points, seginds = sortedintersection(segs)
+  @test length(points) == 2
+  @test length(seginds) == 2
+  @test Set(seginds[1]) == Set([1, 2])
+  @test Set(seginds[2]) == Set([1, 6, 5])
+
+  # multiple intersections, endpoints as intersections
+  if T === Float64
+    segs =
+      Segment.([
+        (cart(9, 13), cart(6, 9)),
+        (cart(2, 12), cart(9, 4.8)),
+        (cart(12, 11), cart(4, 7)),
+        (cart(2.5, 10), cart(12.5, 2)),
+        (cart(13, 6), cart(10, 4)),
+        (cart(10.5, 5.5), cart(9, 1)),
+        (cart(10, 4), cart(11, -1)),
+        (cart(10, 3), cart(10, 5))
+      ])
+    points, seginds = sortedintersection(segs)
+    @test length(points) == 4
+    @test length(seginds) == 4
+    @test points[3] ≈ cart(9, 4.8)
+    @test points[4] ≈ cart(10, 4)
+    @test Set(seginds[1]) == Set([4, 3])
+    @test Set(seginds[2]) == Set([2, 3])
+    @test Set(seginds[3]) == Set([4, 2])
+    @test Set(seginds[4]) == Set([4, 5, 6, 7, 8])
+  end
+
+  # finds all intersections in a grid
+  n = 10
+  horizontal = [Segment(cart(1, i), cart(n, i)) for i in 1:n]
+  vertical = [Segment(cart(i, 1), cart(i, n)) for i in 1:n]
+  segs = [horizontal; vertical]
+  points, seginds = sortedintersection(segs)
+  @test length(points) == n * n - 4
+  @test length(seginds) == n * n - 4
+  @test all(==(2), length.(seginds))
+
+  # number of intersections is invariant under rotations
+  for θ in T(π / 6):T(π / 6):T(2π - π / 6)
+    # rotation by π in Float32 is not robust, skips test
+    T === Float32 && θ == T(π) && continue
+    θpoints, θseginds = sortedintersection(segs |> Rotate(θ))
+    @test length(θpoints) == n * n - 4
+    @test length(θseginds) == n * n - 4
+    @test all(==(2), length.(θseginds))
+  end
+
+  # tests coverage for when intervals don't overlap
+  segs = [
+    Segment(cart(0, 2), cart(2, 0)),
+    Segment(cart(0, 0), cart(2, 2)),
+    Segment(cart(3, 1), cart(3, 3)),
+    Segment(cart(3, 3), cart(3, 3))
+  ]
+  points, seginds = sortedintersection(segs)
+  @test length(points) == 1
+
+  # inference test
+  segs = facets(cartgrid(10, 10))
+  @inferred (Meshes.pairwiseintersect(segs))
+end
+
 @testitem "isthreaded" setup = [Setup] begin
   if Threads.nthreads() > 1
     @test Meshes.isthreaded()