Add Bentley-Ottman Algorithm

Project.toml

@@ -5,6 +5,7 @@ version = "0.52.2"
 
 [deps]
 Bessels = "0e736298-9ec6-45e8-9647-e4fc86a2fe38"
+BinaryTrees = "289e92be-c05a-437b-9e67-5b0c799891f8"
 CircularArrays = "7a955b69-7140-5f4e-a0ed-f168c5e2e749"
 Colorfy = "03fe91ce-8ec6-4610-8e8d-e7491ccca690"
 CoordRefSystems = "b46f11dc-f210-4604-bfba-323c1ec968cb"
@@ -22,8 +23,15 @@ TiledIteration = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
 TransformsBase = "28dd2a49-a57a-4bfb-84ca-1a49db9b96b8"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
+[weakdeps]
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+
+[extensions]
+MeshesMakieExt = "Makie"
+
 [compat]
 Bessels = "0.2"
+BinaryTrees = "0.1.1"
 CircularArrays = "1.3"
 Colorfy = "1.0"
 CoordRefSystems = "0.15"
@@ -42,9 +50,3 @@ TiledIteration = "0.5"
 TransformsBase = "1.6"
 Unitful = "1.17"
 julia = "1.9"
-
-[extensions]
-MeshesMakieExt = "Makie"
-
-[weakdeps]
-Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"

src/Meshes.jl

@@ -4,6 +4,7 @@
 
 module Meshes
 
+using BinaryTrees
 using CoordRefSystems
 using StaticArrays
 using SparseArrays

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,149 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+"""
+    bentleyottmann(segments)
+
+Compute pairwise intersections between n `segments` in
+O(n⋅log(n)) time using Bentley-Ottmann sweep line algorithm.
+
+## References
+
+* Bentley & Ottmann 1979. [Algorithms for reporting and counting
+  geometric intersections](https://ieeexplore.ieee.org/document/1675432)
+"""
+function bentleyottmann(segments)
+  # orient segments
+  segs = map(segments) do s
+    a, b = extrema(s)
+    a > b ? (b, a) : (a, b)
+  end
+
+  # retrieve types
+  P = eltype(first(segs))
+  S = Tuple{P,P}
+
+  # initialization
+  𝒬 = BinaryTrees.AVLTree{P}()
+  ℛ = BinaryTrees.AVLTree{S}()
+  ℬ = Dict{P,Vector{S}}()
+  ℰ = Dict{P,Vector{S}}()
+  ℳ = Dict{P,Vector{S}}()
+  lookup = Dict{S,Int}()
+  for (i, (a, b)) in enumerate(segs)
+    BinaryTrees.insert!(𝒬, a)
+    BinaryTrees.insert!(𝒬, b)
+    haskey(ℬ, a) ? push!(ℬ[a], (a, b)) : (ℬ[a] = [(a, b)])
+    haskey(ℰ, b) ? push!(ℰ[b], (a, b)) : (ℰ[b] = [(a, b)])
+    lookup[(a, b)] = i
+  end
+
+  # sweep line algorithm
+  points = Vector{P}()
+  seginds = Vector{Vector{Int}}()
+  while !_isempty(𝒬)
+    p = BinaryTrees.key(BinaryTrees.minnode(𝒬))
+    BinaryTrees.delete!(𝒬, p)
+    _handle!(points, seginds, lookup, p, 𝒬, ℛ, ℬ, ℰ, ℳ)
+  end
+  points, seginds
+end
+
+function _handle!(points, seginds, lookup, p, 𝒬, ℛ, ℬ, ℰ, ℳ)
+  ℬₚ = _segmentswith(ℬ, p)
+  ℰₚ = _segmentswith(ℰ, p)
+  ℳₚ = _segmentswith(ℳ, p)
+  _processend!(ℰₚ, 𝒬, ℛ, ℳ)
+  _processbeg!(ℬₚ, 𝒬, ℛ, ℳ)
+  _processmid!(ℳₚ, 𝒬, ℛ, ℳ)
+  inds = [lookup[s] for s in ℬₚ ∪ ℰₚ ∪ ℳₚ]
+  if !isempty(inds)
+    push!(points, p)
+    push!(seginds, inds)
+  end
+end
+
+function _processbeg!(ℬₚ, 𝒬, ℛ, ℳ)
+  for s in ℬₚ
+    BinaryTrees.insert!(ℛ, s)
+  end
+  for s in ℬₚ
+    prev, next = BinaryTrees.prevnext(ℛ, s)
+    if !isnothing(prev)
+      _newevent!(𝒬, ℳ, BinaryTrees.key(prev), s)
+    end
+    if !isnothing(next)
+      _newevent!(𝒬, ℳ, s, BinaryTrees.key(next))
+    end
+  end
+end
+
+function _processend!(ℰₚ, 𝒬, ℛ, ℳ)
+  for s in ℰₚ
+    prev, next = BinaryTrees.prevnext(ℛ, s)
+    if !isnothing(prev) && !isnothing(next)
+      _newevent!(𝒬, ℳ, BinaryTrees.key(next), BinaryTrees.key(prev))
+    end
+    BinaryTrees.delete!(ℛ, s)
+  end
+end
+
+function _processmid!(ℳₚ, 𝒬, ℛ, ℳ)
+  for s in ℳₚ
+    prev, _ = BinaryTrees.prevnext(ℛ, s)
+    if !isnothing(prev)
+      # find segments r and u
+      r = prev
+      _, u = BinaryTrees.prevnext(ℛ, BinaryTrees.key(prev))
+
+      # remove crossing points rs and tu from event queue
+      if !isnothing(r)
+        _rmevent!(𝒬, BinaryTrees.key(r), s)
+      end
+      if !isnothing(u)
+        _rmevent!(𝒬, BinaryTrees.key(u), BinaryTrees.key(prev))
+      end
+
+      # add crossing points rt and su to event queue
+      if !isnothing(r)
+        _newevent!(𝒬, ℳ, BinaryTrees.key(r), BinaryTrees.key(prev))
+      end
+      if !isnothing(u)
+        _newevent!(𝒬, ℳ, BinaryTrees.key(u), s)
+      end
+    end
+  end
+end
+
+function _newevent!(𝒬, ℳ, s₁, s₂)
+  intersection(Segment(s₁), Segment(s₂)) do I
+    if type(I) == Crossing || type(I) == EdgeTouching
+      p = get(I)
+      BinaryTrees.insert!(𝒬, p)
+      if haskey(ℳ, p)
+        push!(ℳ[p], s₁, s₂)
+      else
+        ℳ[p] = [s₁, s₂]
+      end
+    end
+    nothing
+  end
+end
+
+function _rmevent!(𝒬, s₁, s₂)
+  intersection(Segment(s₁), Segment(s₂)) do I
+    if type(I) == Crossing || type(I) == EdgeTouching
+      BinaryTrees.delete!(𝒬, get(I))
+    end
+    nothing
+  end
+end
+
+_isempty(𝒬) = isnothing(BinaryTrees.root(𝒬))
+
+function _segmentswith(𝒟, p)
+  P = typeof(p)
+  S = Tuple{P,P}
+  get(𝒟, p, S[])
+end

test/utils.jl

@@ -55,3 +55,23 @@
   p = latlon(c) |> Proj(Cartesian)
   @inferred Meshes.withcrs(p, c, LatLon)
 end
+
+@testitem "Bentley-Ottmann" setup = [Setup] begin
+  segs = [
+    Segment(cart(0, 0), cart(1.1, 1.1)),
+    Segment(cart(1, 0), cart(0, 1)),
+    Segment(cart(0, 0), cart(0, 1)),
+    Segment(cart(0, 0), cart(1, 0)),
+    Segment(cart(0, 1), cart(1, 1)),
+    Segment(cart(1, 0), cart(1, 1))
+  ]
+
+  points, seginds = Meshes.bentleyottmann(segs)
+
+  @inferred Meshes.bentleyottmann(segs)
+  @test all(points .≈ [cart(0, 0), cart(0, 1), cart(0.5, 0.5), cart(1, 0), cart(1, 1), cart(1.1, 1.1)])
+  @test length(points) == 6
+  @test length(seginds) == 6
+  @test maximum(reduce(vcat, seginds)) == length(segs)
+  @test seginds == [[1, 3, 4], [2, 5, 3], [1, 2], [6, 2, 4], [5, 6, 1], [1]]
+end