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..9653c1854 --- /dev/null +++ b/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(vertices(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 diff --git a/test/utils.jl b/test/utils.jl index 51e34f17d..c69912dfd 100644 --- a/test/utils.jl +++ b/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()