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,111 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+"""
+    intervalsweep(segments; [digits])
+
+Compute pairwise intersections between n `segments`
+with `digits` precision in O(n⋅log(n+k)) time using
+an x-interval 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 intervalsweep(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
+  𝒮 = SweepLineQueue(segs)
+
+  points, seginds = handle!(𝒮; digits=digits)
+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
+
+# ----------------
+# DATA STRUCTURES
+# ----------------
+
+struct SweepLineInterval{T<:Number}
+  start::T
+  stop::T
+  segment::Any
+  index::Int
+end
+
+struct SweepLineQueue
+  intervals::Vector{SweepLineInterval}
+end
+Base.length(𝒮::SweepLineQueue) = length(𝒮.intervals)
+Base.getindex(𝒮::SweepLineQueue, i::Int) = 𝒮.intervals[i]
+
+function SweepLineQueue(segs::Vector{<:Tuple{Point,Point}})
+  intervals = Vector{SweepLineInterval}()
+  for (i, seg) in enumerate(segs)
+    x₁, _ = CoordRefSystems.values(coords(seg[1]))
+    x₂, _ = CoordRefSystems.values(coords(seg[2]))
+    push!(intervals, SweepLineInterval(min(x₁, x₂), max(x₁, x₂), seg, i))
+  end
+  SweepLineQueue(sort!(intervals, by=s -> s.start))
+end
+
+function overlaps(i₁::SweepLineInterval, i₂::SweepLineInterval)
+  i₁.start ≤ i₂.start && i₂.stop ≥ i₁.start
+end
+# ----------------
+# SWEEP LINE HANDLER
+# ----------------
+
+"""
+  handle!(𝒮::SweepLineQueue)
+
+Iterate through the sweep line queue and compute all intersection points
+between overlapping intervals. Returns a tuple of intersection points and
+the sets of segment indices that intersect at each point.
+"""
+function handle!(𝒮::SweepLineQueue; digits=10)
+  𝐺 = Dict{Point,Set{Int}}()
+  n = length(𝒮)
+  for i in 1:n
+    current = 𝒮[i]
+    for k in (i + 1):n
+      candidate = 𝒮[k]
+      # If the intervals no longer overlap, break out of the inner loop
+      if !overlaps(current, candidate)
+        break
+      end
+      # Check if the segments actually intersect
+      I = intersection(Segment(current.segment), Segment(candidate.segment)) do 𝑖
+        t = type(𝑖)
+        (t === Crossing || t === EdgeTouching) ? get(𝑖) : nothing
+      end
+      isnothing(I) || _addintersection!(𝐺, I, current.index, candidate.index; digits=digits)
+    end
+  end
+  (collect(keys(𝐺)), collect(values(𝐺)))
+end
+function _addintersection!(𝐺::Dict{Point,Set{Int}}, I::Point, index₁::Int, index₂::Int; digits=10)
+  p = coordround(I, digits=digits)
+  if haskey(𝐺, p)
+    union!(𝐺[p], index₁)
+    union!(𝐺[p], index₂)
+  else
+    𝐺[p] = Set([index₁, index₂])
+  end
+end

test/utils.jl

@@ -80,6 +80,100 @@ end
   @inferred Meshes.coordround(p₁, digits=10)
 end
 
+@testitem "intervalsweep" setup = [Setup] begin
+  # helper to sort points and seginds by point coordinates
+  function sortedintersection(segs)
+    points, inds = Meshes.intervalsweep(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.intervalsweep(segs))
+end
+
 @testitem "isthreaded" setup = [Setup] begin
   if Threads.nthreads() > 1
     @test Meshes.isthreaded()