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,90 @@
+# ------------------------------------------------------------------
+# 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 using
+a sweep line algorithm. Similar to an optimal
+Bentley-Ottmann algorithm in sparse systems,
+and closer to O(n²) in dense systems.
+
+By default, set `digits` based on the absolute
+tolerance of the length type of the 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 seg
+    a, b = coordround.(extrema(seg), digits=digits)
+    a > b ? (b, a) : (a, b)
+  end
+
+  starts, stops = map(segs) do seg
+    x₁, _ = CoordRefSystems.values(coords(seg[1]))
+    x₂, _ = CoordRefSystems.values(coords(seg[2]))
+    (x₁, x₂)
+  end |> x -> (first.(x), last.(x))
+  # sort indices based on start x coordinates
+  sortedindices = sortperm(starts)
+  # reorder everything based on sorted indices
+  starts, stops, segs = getindex.((starts, stops, segs), Ref(sortedindices))
+  # keep track of original indices
+  n = length(segs)
+  oldindices = collect(1:n)[sortedindices]
+
+  # primary sweepline algorithm
+  𝐺 = Dict{Point,Vector{Int}}()
+  for i in eachindex(segs)
+    for k in (i + 1):n
+      I = _checkintersection(i, k, starts, stops, segs)
+
+      # break if no overlap, continue if no intersection, add intersection otherwise
+      if I == :break
+        break
+      elseif I == :continue
+        continue
+      else
+        _addintersection!(𝐺, I, oldindices[i], oldindices[k]; digits=digits)
+      end
+    end
+  end
+  (collect(keys(𝐺)), collect(values(𝐺)))
+end
+
+overlaps(startᵢ, stopᵢ, startₖ) = (startᵢ ≤ startₖ ≤ stopᵢ)
+
+function _checkintersection(i, k, starts, stops, segs)
+  overlap = overlaps(starts[i], stops[i], starts[k])
+  overlap || return :break
+  intersection(Segment(segs[i]), Segment(segs[k])) do 𝑖
+    t = type(𝑖)
+    (t === Crossing || t === EdgeTouching) ? get(𝑖) : :continue
+  end
+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)
+  seg = first(segments)
+  ℒ = lentype(seg)
+  τ = ustrip(eps(ℒ))
+  round(Int, 0.8 * (-log10(τ))) # 0.8 is a heuristic to avoid numerical issues
+end
+
+# add an intersection point to the dictionary with segment indices
+function _addintersection!(𝐺, I::Point, index₁::Int, index₂::Int; digits=10)
+  p = coordround(I, digits=digits)
+  if haskey(𝐺, p)
+    append!(𝐺[p], (index₁, index₂))
+    unique!(𝐺[p])
+  else
+    𝐺[p] = Vector{Int}([index₁, index₂])
+  end
+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 Set(length.(seginds)) == Set([2])
+
+  # 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 Set(length.(θseginds)) == Set([2])
+  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 Nothing (Meshes.pairwiseintersect(segs))
+end
+
 @testitem "isthreaded" setup = [Setup] begin
   if Threads.nthreads() > 1
     @test Meshes.isthreaded()