diff --git a/Project.toml b/Project.toml index c02eb8b17..39f7af240 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "0.53.5" [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.2" CircularArrays = "1.3" Colorfy = "1.1" CoordRefSystems = "0.16 - 0.17" @@ -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" diff --git a/src/Meshes.jl b/src/Meshes.jl index db2612395..1f2f6efcb 100644 --- a/src/Meshes.jl +++ b/src/Meshes.jl @@ -4,6 +4,7 @@ module Meshes +using BinaryTrees using CoordRefSystems using StaticArrays using SparseArrays diff --git a/src/utils.jl b/src/utils.jl index adc3a611f..1a74be7cd 100644 --- a/src/utils.jl +++ b/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") diff --git a/src/utils/intersect.jl b/src/utils/intersect.jl new file mode 100644 index 000000000..6dc3b7f3e --- /dev/null +++ b/src/utils/intersect.jl @@ -0,0 +1,300 @@ +# ------------------------------------------------------------------ +# Licensed under the MIT License. See LICENSE in the project root. +# ------------------------------------------------------------------ +""" + bentleyottmann(segments; [digits]) + +Compute pairwise intersections between n `segments` +with `digits` precision in O(n⋅log(n)) time using +Bentley-Ottmann sweep line algorithm. + +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 bentleyottmann(segments; digits=_digits(segments)) + TOL = 1 / 10^digits # precomputed tolerance for floating point comparisons + + # orient segments + segs = map(segments) do s + a, b = coordround.(extrema(s); digits) + a > b ? (b, a) : (a, b) + end + + minp, maxp = segs |> Meshes.boundingbox |> Stretch(1.05) |> extrema + ymin = CoordRefSystems.values(coords(minp))[2] + ymax = CoordRefSystems.values(coords(maxp))[2] + ybounds = (ymin, ymax) + + # retrieve types + T = lentype(first(segs)[1]) + P = typeof(first(segs)[1]) + S = Tuple{P,P} + + # initialization + 𝒬 = BinaryTrees.AVLTree{P}() + ℛ = BinaryTrees.AVLTree{_SweepSegment{T,P}}() + ℬ = 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}}() + sweepline = _initsweep(segs, ybounds) + while !BinaryTrees.isempty(𝒬) + # current point (or event) + p = BinaryTrees.key(BinaryTrees.minnode(𝒬)) + + sweepline.point = p # update sweepline position + + # delete point from event queue + BinaryTrees.delete!(𝒬, p) + # handle event, i.e. update 𝒬, ℛ and ℳ + ℬₚ = get(ℬ, p, S[]) # segments with p at the begin + ℰₚ = get(ℰ, p, S[]) # segments with p at the end + ℳₚ = S[] + _findintersections!(ℳₚ, ℛ, sweepline, TOL) # segments with p at the middle + activesegs = Set(ℬₚ ∪ ℳₚ) + # report intersections + if !isempty(activesegs) || !isempty(ℰₚ) + inds = Set{Int}() + for s in activesegs ∪ ℰₚ + push!(inds, lookup[s]) + end + push!(points, p) + push!(seginds, collect(inds)) + end + + # handle status line + _handlestatus!(ℛ, ℬₚ, ℳₚ, ℰₚ, sweepline, p, TOL) + + if isempty(activesegs) + for s in ℰₚ + sₗ, sᵣ = BinaryTrees.prevnext(ℛ, _SweepSegment(s, sweepline)) + isnothing(sₗ) || isnothing(sᵣ) || _newevent!(𝒬, p, _keyseg(sₗ), _keyseg(sᵣ), digits) + end + else + BinaryTrees.isempty(ℛ) || _handlebottom!(activesegs, ℛ, sweepline, 𝒬, p, digits) + + BinaryTrees.isempty(ℛ) || _handletop!(activesegs, ℛ, sweepline, 𝒬, p, digits) + end + end + + (points, seginds) +end + +## +## handling functions +## + +function _handlestatus!(ℛ, ℬₚ, ℳₚ, ℰₚ, sweepline, p, TOL) + # nudge back to get correct ordering for removal + sweepline.point = _nudge(p, TOL, -) + for s in reverse(ℰₚ ∪ ℳₚ) + segsweep = _SweepSegment(s, sweepline) + isnothing(BinaryTrees.search(ℛ, segsweep)) || BinaryTrees.delete!(ℛ, segsweep) + end + + # nudge forward to get correct ordering for adding + sweepline.point = _nudge(p, TOL, +) + + for s in ℬₚ ∪ ℳₚ + BinaryTrees.insert!(ℛ, _SweepSegment(s, sweepline)) + end +end + +function _handlebottom!(activesegs, ℛ, sweepline, 𝒬, p, digits) + s′ = BinaryTrees.key(_minsearch(ℛ, activesegs, sweepline)) + + sₗ, _ = !isnothing(s′) ? BinaryTrees.prevnext(ℛ, s′) : (nothing, nothing) + if !isnothing(sₗ) + _newevent!(𝒬, p, _keyseg(sₗ), _segment(s′), digits) + end +end + +function _handletop!(activesegs, ℛ, sweepline, 𝒬, p, digits) + s″ = BinaryTrees.key(_maxsearch(ℛ, activesegs, sweepline)) + + _, sᵤ = !isnothing(s″) ? BinaryTrees.prevnext(ℛ, s″) : (nothing, nothing) + if !isnothing(sᵤ) + _newevent!(𝒬, p, _segment(s″), _keyseg(sᵤ), digits) + end +end + +# ----------------- +# HELPER FUNCTIONS +# ----------------- + +function _newevent!(𝒬, p, s₁, s₂, digits) + intersection(Segment(s₁), Segment(s₂)) do I + if type(I) == Crossing || type(I) == EdgeTouching + p′ = coordround(get(I); digits) + if p′ ≥ p + BinaryTrees.insert!(𝒬, p′) + end + end + end +end + +# find segments that intersect with the point p +function _findintersections!(ℳₚ, ℛ, sweepline, TOL) + tol = TOL * unit(_sweepx(sweepline)) # ensure TOL is in the same unit as p + _search!(BinaryTrees.root(ℛ), ℳₚ, sweepline, tol) + ℳₚ +end + +function _search!(node, ℳₚ, sweepline, TOL) + isnothing(node) && return + seg = _segment(BinaryTrees.key(node)) + x₂, y₂ = CoordRefSystems.values(coords(seg[2])) + x, y = CoordRefSystems.values(coords(_sweeppoint(sweepline))) + + # Ensure the point is not the endpoint (avoids duplicates) + skip = (x₂ - TOL ≤ x ≤ x₂ + TOL) && (y₂ - TOL ≤ y ≤ y₂ + TOL) + I = intersect(Segment(seg), Segment(_sweepline(sweepline))) + + if isnothing(I) # segment ends just before sweepline (this is expected bc of our nudging) + ŷ = y₂ + elseif I isa Segment # segment means vertical + ŷ = y + else + _, ŷ = CoordRefSystems.values(coords(I)) + end + dy = y - ŷ # difference between the point and the segment + if abs(dy) ≤ TOL || I isa Segment + skip || push!(ℳₚ, seg) + end + + # using difference in y to determine the side of the segment + + if dy < -TOL + _search!(BinaryTrees.left(node), ℳₚ, sweepline, TOL) + elseif dy > TOL + _search!(BinaryTrees.right(node), ℳₚ, sweepline, TOL) + else + # if the point is on the segment, check both sides for adjacents + _search!(BinaryTrees.left(node), ℳₚ, sweepline, TOL) + _search!(BinaryTrees.right(node), ℳₚ, sweepline, TOL) + end +end + +# find the minimum segment among active segments in tree +function _minsearch(ℛ, activesegs, sweepline) + activeordered = sort([_SweepSegment(s, sweepline) for s in activesegs]) + BinaryTrees.search(ℛ, activeordered[begin]) +end + +# find the maximum segment among active segments in tree +function _maxsearch(ℛ, activesegs, sweepline) + activeordered = sort([_SweepSegment(s, sweepline) for s in activesegs]) + BinaryTrees.search(ℛ, activeordered[end]) +end + +# compute rounding digits +function _digits(segments) + s = first(segments) + ℒ = lentype(s) + τ = ustrip(atol(ℒ)) + round(Int, -log10(τ)) - 1 +end + +# convenience function to get the segment from the node +function _keyseg(segment) + _segment(BinaryTrees.key(segment)) +end +# nudge the sweepline to get correct ℛ ordering +function _nudge(p, TOL, operation) + x, y = CoordRefSystems.values(coords(p)) + nudgefactor = unit(x) * TOL * 1 + Point(operation(x, nudgefactor), operation(y, nudgefactor)) +end + +# ---------------- +# DATA STRUCTURES +# ---------------- + +# tracks sweepline and current y position for searching +mutable struct _SweepLine{P<:Point,T} + point::P + ybounds::Tuple{T,T} +end +_sweeppoint(sweepline::_SweepLine) = getfield(sweepline, :point) +_sweepx(sweepline::_SweepLine) = CoordRefSystems.values(coords(_sweeppoint(sweepline)))[1] +_sweepy(sweepline::_SweepLine) = CoordRefSystems.values(coords(_sweeppoint(sweepline)))[2] +_sweepbounds(sweepline::_SweepLine) = getfield(sweepline, :ybounds) +_sweepline(sweepline::_SweepLine) = + (Point(_sweepx(sweepline), sweepline.ybounds[1]), Point(_sweepx(sweepline), sweepline.ybounds[2])) + +# compute the intersection of a segment with the sweepline +function _sweepintersect(seg, sweepline) + p₁, p₂ = coords.(seg) + x₁, y₁ = CoordRefSystems.values(p₁) + x₂, y₂ = CoordRefSystems.values(p₂) + T = eltype(x₁) + + x, y = CoordRefSystems.values(coords(_sweeppoint(sweepline))) + + # if vertical, return the y coordinate of current active point or segment + if abs(x₁ - x₂) < atol(T) + return T(min(y, y₂)) # vertical goes at end + end + + I = intersect(Segment(seg), Segment(_sweepline(sweepline))) + + if isnothing(I) + return y₂ # segment is on the sweepline + end + + x, y = CoordRefSystems.values(coords(I)) + y +end + +# takes input segment and assigns where it intersects sweepline +mutable struct _SweepSegment{T,P<:Point} + seg::Tuple{P,P} + xintersect::T +end + +# constructor for _SweepSegment +function _SweepSegment(seg::Tuple{P,P}, sweepline::_SweepLine) where {P<:Point} + xintersect = _sweepintersect(seg, sweepline) + T = typeof(xintersect) + _SweepSegment{T,P}(seg, xintersect) +end + +_segment(sweepsegment::_SweepSegment) = sweepsegment.seg +_xintersect(sweepsegment::_SweepSegment) = getfield(sweepsegment, :xintersect) + +Base.isless(a::_SweepSegment, b::_SweepSegment) = begin + # if segments same, return false + if _segment(a) == _segment(b) + return false + end + ya, yb = _xintersect(a), _xintersect(b) + + # if segments are separated over y, check ya < yb + if abs(ya - yb) > atol(eltype(ya)) + return ya < yb + end + + # fallback to x-coordinates + _segment(a) != _segment(b) && _segment(a) < _segment(b) +end + +# initialize the sweep line with them minimum +function _initsweep(segs, bounds) + U = lentype(coords(segs[1][1])) + _SweepLine(Point(U(-Inf), U(-Inf)), bounds) +end diff --git a/test/utils.jl b/test/utils.jl index 9927aec84..b2ba52496 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -62,3 +62,67 @@ @test Meshes.coordround(p₂, digits=10) == p₁ @inferred Meshes.coordround(p₁, digits=10) end + +@testitem "bentleyottmann" setup = [Setup] begin + # basic check with a small number of segments + segs = Segment.([(cart(0, 0), cart(2, 2)), (cart(0, 2), cart(2, 0)), (cart(0, 1), cart(0.5, 1))]) + points, seginds = Meshes.bentleyottmann(segs) + @test length(points) == 7 + @test length(seginds) == 7 + + 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 = 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 Set.(seginds[[3, 5]]) == [Set([1, 2]), Set([1, 6, 5])] + + 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 = Meshes.bentleyottmann(segs) + @test length(points) == 17 + @test length(seginds) == 17 + @test Set(reduce(vcat, seginds)) == Set(1:8) + @test points[findfirst(p -> p ≈ cart(10, 4), points)] ≈ cart(10, 4) + @test Set(seginds[findfirst(p -> p ≈ cart(10, 4), points)]) == Set([4, 5, 6, 7, 8]) + @test Set(seginds[findfirst(p -> p ≈ cart(9, 4.8), points)]) == Set([4, 2]) + + # 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 = Meshes.bentleyottmann(segs) + @test length(points) == 100 + @test length(seginds) == 100 + @test Set(length.(seginds)) == Set([2]) + + # result is invariant under rotations + for θ in T(π / 6):T(π / 6):T(2π - π / 6) + θpoints, θseginds = Meshes.bentleyottmann(segs |> Rotate(θ)) + @test length(θpoints) == 100 + @test length(θseginds) == 100 + @test Set(length.(θseginds)) == Set([2]) + end + + # inference test + segs = facets(cartgrid(10, 10)) + @inferred Meshes.bentleyottmann(segs) +end