From da20a780fb4fdeb67a839b8e70112983e11eb767 Mon Sep 17 00:00:00 2001 From: souma4 Date: Tue, 30 Sep 2025 10:32:27 -0600 Subject: [PATCH 01/21] new PR for line segment intersections. Simpler, not quite as efficient method. --- src/utils/intersect.jl | 110 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 110 insertions(+) create mode 100644 src/utils/intersect.jl diff --git a/src/utils/intersect.jl b/src/utils/intersect.jl new file mode 100644 index 000000000..15ae47c6c --- /dev/null +++ b/src/utils/intersect.jl @@ -0,0 +1,110 @@ +# ------------------------------------------------------------------ +# 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+k)⋅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)) + # 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 From 8aa5cfc53b8d62b40ff7188685aeb66ccfbb0d7e Mon Sep 17 00:00:00 2001 From: souma4 Date: Tue, 30 Sep 2025 10:35:26 -0600 Subject: [PATCH 02/21] add tests --- test/utils.jl | 84 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 84 insertions(+) diff --git a/test/utils.jl b/test/utils.jl index 51e34f17d..cb8b53b23 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -80,6 +80,90 @@ end @inferred Meshes.coordround(p₁, digits=10) end +@testitem "bentleyottmann" setup = [Setup] begin + # helper to sort points and seginds by point coordinates + function sortedintersection(segs) + points, inds = Meshes.bentleyottmann(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 + + # inference test + segs = facets(cartgrid(10, 10)) + @inferred Nothing (Meshes.bentleyottmann(segs)) +end + @testitem "isthreaded" setup = [Setup] begin if Threads.nthreads() > 1 @test Meshes.isthreaded() From 7f19f82b0ef084a849c313d470bdcaa48958dc45 Mon Sep 17 00:00:00 2001 From: souma4 Date: Tue, 30 Sep 2025 10:43:22 -0600 Subject: [PATCH 03/21] update utils.jl to include intersect --- src/utils.jl | 1 + 1 file changed, 1 insertion(+) 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") From 35a062b29d5d1c75ff4dc6a7c162023be7a3ea72 Mon Sep 17 00:00:00 2001 From: souma4 Date: Tue, 30 Sep 2025 11:20:34 -0600 Subject: [PATCH 04/21] update names to intervalsweep since this is no longer 'true' bentley-ottmann --- src/utils/intersect.jl | 9 +++++---- test/utils.jl | 6 +++--- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/src/utils/intersect.jl b/src/utils/intersect.jl index 15ae47c6c..4a24d6392 100644 --- a/src/utils/intersect.jl +++ b/src/utils/intersect.jl @@ -3,11 +3,12 @@ # ------------------------------------------------------------------ """ - bentleyottmann(segments; [digits]) + intervalsweep(segments; [digits]) Compute pairwise intersections between n `segments` -with `digits` precision in O((n+k)⋅log(n)) time using -Bentley-Ottmann sweep line algorithm. +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. @@ -17,7 +18,7 @@ tolerance of the length type of the segments. * Bentley & Ottmann 1979. [Algorithms for reporting and counting geometric intersections](https://ieeexplore.ieee.org/document/1675432) """ -function bentleyottmann(segments; digits=_digits(segments)) +function intervalsweep(segments; digits=_digits(segments)) # orient segments and round coordinates segs = map(segments) do seg a, b = coordround.(extrema(seg), digits=digits) diff --git a/test/utils.jl b/test/utils.jl index cb8b53b23..fcf5f58bd 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -80,10 +80,10 @@ end @inferred Meshes.coordround(p₁, digits=10) end -@testitem "bentleyottmann" setup = [Setup] begin +@testitem "intervalsweep" setup = [Setup] begin # helper to sort points and seginds by point coordinates function sortedintersection(segs) - points, inds = Meshes.bentleyottmann(segs) + points, inds = Meshes.intervalsweep(segs) perm = sortperm(points) points[perm], inds[perm] end @@ -161,7 +161,7 @@ end # inference test segs = facets(cartgrid(10, 10)) - @inferred Nothing (Meshes.bentleyottmann(segs)) + @inferred Nothing (Meshes.intervalsweep(segs)) end @testitem "isthreaded" setup = [Setup] begin From 746382a4e3e3dcdcbafd53e8a529cac37674e290 Mon Sep 17 00:00:00 2001 From: souma4 Date: Tue, 30 Sep 2025 11:46:37 -0600 Subject: [PATCH 05/21] add test for coverage when segments stop overlapping --- test/utils.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/test/utils.jl b/test/utils.jl index fcf5f58bd..ae57f8e0b 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -159,6 +159,16 @@ end @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)) From 174ff95195846d1c55602f576a901dd750339d55 Mon Sep 17 00:00:00 2001 From: souma4 Date: Fri, 3 Oct 2025 13:12:50 -0600 Subject: [PATCH 06/21] incorporated review feedback. sweep1D is a descriptive name and leaves namespace open for a 2D sweep (if needed) in the future. Also fixed overlaps logic to properly break early --- src/utils/intersect.jl | 53 ++++++++++++++++++------------------------ test/utils.jl | 6 ++--- 2 files changed, 26 insertions(+), 33 deletions(-) diff --git a/src/utils/intersect.jl b/src/utils/intersect.jl index 4a24d6392..93b65cdf6 100644 --- a/src/utils/intersect.jl +++ b/src/utils/intersect.jl @@ -3,12 +3,13 @@ # ------------------------------------------------------------------ """ - intervalsweep(segments; [digits]) + pairwiseintersect(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. +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. @@ -18,19 +19,27 @@ tolerance of the length type of the segments. * Bentley & Ottmann 1979. [Algorithms for reporting and counting geometric intersections](https://ieeexplore.ieee.org/document/1675432) """ -function intervalsweep(segments; digits=_digits(segments)) +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 - 𝒮 = SweepLineQueue(segs) + sweep1D!(_initqueue(segs); digits=digits) +end - points, seginds = handle!(𝒮; digits=digits) +function _initqueue(segs::Vector{<:Tuple{Point,Point}}) + 𝒬 = Vector{SweepLineInterval}() + for (i, seg) in enumerate(segs) + x₁, _ = CoordRefSystems.values(coords(seg[1])) + x₂, _ = CoordRefSystems.values(coords(seg[2])) + push!(𝒬, SweepLineInterval(min(x₁, x₂), max(x₁, x₂), seg, i)) + end + sort!(𝒬, by=s -> s.start) end -# # compute the number of significant digits based on the segment type -# # this is used to determine the precision of the points +# 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) @@ -49,43 +58,27 @@ struct SweepLineInterval{T<:Number} 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 + i₁.start ≤ i₂.start && i₁.stop ≥ i₂.start end # ---------------- # SWEEP LINE HANDLER # ---------------- """ - handle!(𝒮::SweepLineQueue) + sweep1D!(queue; [digits]) -Iterate through the sweep line queue and compute all intersection points +Iterate through a sweep interval 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) +function sweep1D!(𝒬::Vector{SweepLineInterval}; digits=10) 𝐺 = Dict{Point,Set{Int}}() - n = length(𝒮) + n = length(𝒬) for i in 1:n - current = 𝒮[i] + current = 𝒬[i] for k in (i + 1):n - candidate = 𝒮[k] + candidate = 𝒬[k] # If the intervals no longer overlap, break out of the inner loop if !overlaps(current, candidate) break diff --git a/test/utils.jl b/test/utils.jl index ae57f8e0b..c75a6e37f 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -80,10 +80,10 @@ end @inferred Meshes.coordround(p₁, digits=10) end -@testitem "intervalsweep" setup = [Setup] begin +@testitem "pairwiseintersect" setup = [Setup] begin # helper to sort points and seginds by point coordinates function sortedintersection(segs) - points, inds = Meshes.intervalsweep(segs) + points, inds = Meshes.pairwiseintersect(segs) perm = sortperm(points) points[perm], inds[perm] end @@ -171,7 +171,7 @@ end # inference test segs = facets(cartgrid(10, 10)) - @inferred Nothing (Meshes.intervalsweep(segs)) + @inferred Nothing (Meshes.pairwiseintersect(segs)) end @testitem "isthreaded" setup = [Setup] begin From 63a26f94a91a48767d738eb158273e3a4d5f3a3c Mon Sep 17 00:00:00 2001 From: souma4 Date: Fri, 3 Oct 2025 13:29:00 -0600 Subject: [PATCH 07/21] sweep1D doesn't update the initial queue. --- src/utils/intersect.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/utils/intersect.jl b/src/utils/intersect.jl index 93b65cdf6..d0d49b49f 100644 --- a/src/utils/intersect.jl +++ b/src/utils/intersect.jl @@ -25,7 +25,7 @@ function pairwiseintersect(segments; digits=_digits(segments)) a, b = coordround.(extrema(seg), digits=digits) a > b ? (b, a) : (a, b) end - sweep1D!(_initqueue(segs); digits=digits) + sweep1D(_initqueue(segs); digits=digits) end function _initqueue(segs::Vector{<:Tuple{Point,Point}}) @@ -66,13 +66,13 @@ end # ---------------- """ - sweep1D!(queue; [digits]) + sweep1D(queue; [digits]) Iterate through a sweep interval 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 sweep1D!(𝒬::Vector{SweepLineInterval}; digits=10) +function sweep1D(𝒬::Vector{SweepLineInterval}; digits=10) 𝐺 = Dict{Point,Set{Int}}() n = length(𝒬) for i in 1:n From 03e407a4019b92f3a226a62d4d571efcb8b27854 Mon Sep 17 00:00:00 2001 From: souma4 Date: Thu, 9 Oct 2025 11:07:18 -0600 Subject: [PATCH 08/21] simplified implementation, no new structs, moved from sets of indices to vectors of indices because sets are expensive to initialized --- src/utils/intersect.jl | 102 ++++++++++++++++++----------------------- 1 file changed, 44 insertions(+), 58 deletions(-) diff --git a/src/utils/intersect.jl b/src/utils/intersect.jl index d0d49b49f..21051488d 100644 --- a/src/utils/intersect.jl +++ b/src/utils/intersect.jl @@ -7,7 +7,7 @@ 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 +a sweep line algorithm. Similar to an optimal Bentley-Ottmann algorithm in sparse systems, and closer to O(n²) in dense systems. @@ -25,17 +25,48 @@ function pairwiseintersect(segments; digits=_digits(segments)) a, b = coordround.(extrema(seg), digits=digits) a > b ? (b, a) : (a, b) end - sweep1D(_initqueue(segs); digits=digits) -end -function _initqueue(segs::Vector{<:Tuple{Point,Point}}) - 𝒬 = Vector{SweepLineInterval}() - for (i, seg) in enumerate(segs) + starts, stops = map(segs) do seg x₁, _ = CoordRefSystems.values(coords(seg[1])) x₂, _ = CoordRefSystems.values(coords(seg[2])) - push!(𝒬, SweepLineInterval(min(x₁, x₂), max(x₁, x₂), seg, i)) + (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 - sort!(𝒬, by=s -> s.start) end # compute the number of significant digits based on the segment type @@ -47,58 +78,13 @@ function _digits(segments) 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 - -function overlaps(i₁::SweepLineInterval, i₂::SweepLineInterval) - i₁.start ≤ i₂.start && i₁.stop ≥ i₂.start -end -# ---------------- -# SWEEP LINE HANDLER -# ---------------- - -""" - sweep1D(queue; [digits]) - -Iterate through a sweep interval 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 sweep1D(𝒬::Vector{SweepLineInterval}; 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) +# 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) - union!(𝐺[p], index₁) - union!(𝐺[p], index₂) + append!(𝐺[p], (index₁, index₂)) + unique!(𝐺[p]) else - 𝐺[p] = Set([index₁, index₂]) + 𝐺[p] = Vector{Int}([index₁, index₂]) end end From e8b8585869f5970795d194062f90c3e6be4273f8 Mon Sep 17 00:00:00 2001 From: souma4 Date: Fri, 10 Oct 2025 20:27:12 -0600 Subject: [PATCH 09/21] consolidated code, cleaned up names, removed maps, improved type stability --- src/utils/intersect.jl | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/utils/intersect.jl b/src/utils/intersect.jl index 21051488d..0cd06f765 100644 --- a/src/utils/intersect.jl +++ b/src/utils/intersect.jl @@ -20,27 +20,27 @@ tolerance of the length type of the segments. geometric intersections](https://ieeexplore.ieee.org/document/1675432) """ function pairwiseintersect(segments; digits=_digits(segments)) + P = typeof(vertices(first(segments))[1]) # 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)) + starts = [CoordRefSystems.values(coords(seg[1]))[1] for seg in segs] + stops = [CoordRefSystems.values(coords(seg[2]))[1] for seg in segs] # sort indices based on start x coordinates - sortedindices = sortperm(starts) + inds = sortperm(starts) # reorder everything based on sorted indices - starts, stops, segs = getindex.((starts, stops, segs), Ref(sortedindices)) + starts = starts[inds] + stops = stops[inds] + segs = segs[inds] # keep track of original indices n = length(segs) - oldindices = collect(1:n)[sortedindices] + oldindices = (1:n)[inds] - # primary sweepline algorithm - 𝐺 = Dict{Point,Vector{Int}}() + # sweepline algorithm + 𝐺 = Dict{P,Vector{Int}}() for i in eachindex(segs) for k in (i + 1):n I = _checkintersection(i, k, starts, stops, segs) @@ -58,10 +58,10 @@ function pairwiseintersect(segments; digits=_digits(segments)) (collect(keys(𝐺)), collect(values(𝐺))) end -overlaps(startᵢ, stopᵢ, startₖ) = (startᵢ ≤ startₖ ≤ stopᵢ) +_overlaps(startᵢ, stopᵢ, startₖ) = (startᵢ ≤ startₖ ≤ stopᵢ) function _checkintersection(i, k, starts, stops, segs) - overlap = overlaps(starts[i], stops[i], starts[k]) + overlap = _overlaps(starts[i], stops[i], starts[k]) overlap || return :break intersection(Segment(segs[i]), Segment(segs[k])) do 𝑖 t = type(𝑖) From 744189e3ac0180c672bc07f6235493ab15d455d1 Mon Sep 17 00:00:00 2001 From: souma4 Date: Tue, 14 Oct 2025 15:50:31 -0600 Subject: [PATCH 10/21] disolve _checkintersection function and change P definition --- src/utils/intersect.jl | 27 +++++++++++---------------- 1 file changed, 11 insertions(+), 16 deletions(-) diff --git a/src/utils/intersect.jl b/src/utils/intersect.jl index 0cd06f765..9c8e65ac4 100644 --- a/src/utils/intersect.jl +++ b/src/utils/intersect.jl @@ -20,12 +20,12 @@ tolerance of the length type of the segments. geometric intersections](https://ieeexplore.ieee.org/document/1675432) """ function pairwiseintersect(segments; digits=_digits(segments)) - P = typeof(vertices(first(segments))[1]) # orient segments and round coordinates segs = map(segments) do seg a, b = coordround.(extrema(seg), digits=digits) a > b ? (b, a) : (a, b) end + P = typeof(first(segs)[1]) starts = [CoordRefSystems.values(coords(seg[1]))[1] for seg in segs] stops = [CoordRefSystems.values(coords(seg[2]))[1] for seg in segs] @@ -43,12 +43,17 @@ function pairwiseintersect(segments; digits=_digits(segments)) 𝐺 = Dict{P,Vector{Int}}() for i in eachindex(segs) for k in (i + 1):n - I = _checkintersection(i, k, starts, stops, segs) + overlap = _overlaps(starts[i], stops[i], starts[k]) + # if no overlap, break inner loop + overlap || break - # break if no overlap, continue if no intersection, add intersection otherwise - if I == :break - break - elseif I == :continue + I = intersection(Segment(segs[i]), Segment(segs[k])) do 𝑖 + t = type(𝑖) + (t === Crossing || t === EdgeTouching) ? get(𝑖) : nothing + end + + # continue if no intersection, add intersection otherwise + if isnothing(I) continue else _addintersection!(𝐺, I, oldindices[i], oldindices[k]; digits=digits) @@ -59,16 +64,6 @@ function pairwiseintersect(segments; digits=_digits(segments)) 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) From a07ddbcf1d1c917944e7429f55a767cf08494817 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BAlio=20Hoffimann?= Date: Wed, 15 Oct 2025 13:43:08 -0300 Subject: [PATCH 11/21] basic cleanup --- src/utils/intersect.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/utils/intersect.jl b/src/utils/intersect.jl index 9c8e65ac4..b826da985 100644 --- a/src/utils/intersect.jl +++ b/src/utils/intersect.jl @@ -25,29 +25,28 @@ function pairwiseintersect(segments; digits=_digits(segments)) a, b = coordround.(extrema(seg), digits=digits) a > b ? (b, a) : (a, b) end - P = typeof(first(segs)[1]) starts = [CoordRefSystems.values(coords(seg[1]))[1] for seg in segs] stops = [CoordRefSystems.values(coords(seg[2]))[1] for seg in segs] - # sort indices based on start x coordinates + + # sort segments based on start coordinates inds = sortperm(starts) - # reorder everything based on sorted indices starts = starts[inds] stops = stops[inds] segs = segs[inds] + # keep track of original indices n = length(segs) oldindices = (1:n)[inds] # sweepline algorithm + P = eltype(first(segs)) 𝐺 = Dict{P,Vector{Int}}() for i in eachindex(segs) - for k in (i + 1):n - overlap = _overlaps(starts[i], stops[i], starts[k]) - # if no overlap, break inner loop - overlap || break + for j in (i + 1):length(segs) + _overlaps(starts[i], stops[i], starts[j]) || break - I = intersection(Segment(segs[i]), Segment(segs[k])) do 𝑖 + I = intersection(Segment(segs[i]), Segment(segs[j])) do 𝑖 t = type(𝑖) (t === Crossing || t === EdgeTouching) ? get(𝑖) : nothing end @@ -56,7 +55,7 @@ function pairwiseintersect(segments; digits=_digits(segments)) if isnothing(I) continue else - _addintersection!(𝐺, I, oldindices[i], oldindices[k]; digits=digits) + _addintersection!(𝐺, I, oldindices[i], oldindices[j]; digits=digits) end end end @@ -64,6 +63,7 @@ function pairwiseintersect(segments; digits=_digits(segments)) end _overlaps(startᵢ, stopᵢ, startₖ) = (startᵢ ≤ startₖ ≤ stopᵢ) + # compute the number of significant digits based on the segment type # this is used to determine the precision of the points function _digits(segments) From 5a290dfc255eb2f9dba1b726b603c7c3455388a9 Mon Sep 17 00:00:00 2001 From: souma4 Date: Wed, 15 Oct 2025 11:35:19 -0600 Subject: [PATCH 12/21] minor cleanup --- src/utils/intersect.jl | 22 +++++----------------- 1 file changed, 5 insertions(+), 17 deletions(-) diff --git a/src/utils/intersect.jl b/src/utils/intersect.jl index b826da985..e4ffe8a1f 100644 --- a/src/utils/intersect.jl +++ b/src/utils/intersect.jl @@ -35,27 +35,15 @@ function pairwiseintersect(segments; digits=_digits(segments)) stops = stops[inds] segs = segs[inds] - # keep track of original indices - n = length(segs) - oldindices = (1:n)[inds] - # sweepline algorithm P = eltype(first(segs)) 𝐺 = Dict{P,Vector{Int}}() - for i in eachindex(segs) - for j in (i + 1):length(segs) - _overlaps(starts[i], stops[i], starts[j]) || break - - I = intersection(Segment(segs[i]), Segment(segs[j])) do 𝑖 - t = type(𝑖) - (t === Crossing || t === EdgeTouching) ? get(𝑖) : nothing - end + for i in eachindex(segs), j in (i + 1):length(segs) + _overlaps(starts[i], stops[i], starts[j]) || break - # continue if no intersection, add intersection otherwise - if isnothing(I) - continue - else - _addintersection!(𝐺, I, oldindices[i], oldindices[j]; digits=digits) + intersection(Segment(segs[i]), Segment(segs[j])) do I + if type(I) == Crossing || type(I) == EdgeTouching + _addintersection!(𝐺, I, inds[i], inds[j]; digits=digits) end end end From d381b5d5d52cd3d5b105eb44a621f85ad9868e42 Mon Sep 17 00:00:00 2001 From: souma4 Date: Wed, 15 Oct 2025 13:32:39 -0600 Subject: [PATCH 13/21] fixed addintersection. it gets I --- src/utils/intersect.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils/intersect.jl b/src/utils/intersect.jl index e4ffe8a1f..a938147ab 100644 --- a/src/utils/intersect.jl +++ b/src/utils/intersect.jl @@ -43,7 +43,7 @@ function pairwiseintersect(segments; digits=_digits(segments)) intersection(Segment(segs[i]), Segment(segs[j])) do I if type(I) == Crossing || type(I) == EdgeTouching - _addintersection!(𝐺, I, inds[i], inds[j]; digits=digits) + _addintersection!(𝐺, get(I), inds[i], inds[j]; digits=digits) end end end From efc65c4c823afda42cdd15469cb67d9f131cfec8 Mon Sep 17 00:00:00 2001 From: souma4 Date: Wed, 15 Oct 2025 14:20:57 -0600 Subject: [PATCH 14/21] split out the for loops because continue for inline double for loops continues for i, not j --- src/utils/intersect.jl | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/utils/intersect.jl b/src/utils/intersect.jl index a938147ab..f4222efe8 100644 --- a/src/utils/intersect.jl +++ b/src/utils/intersect.jl @@ -38,13 +38,16 @@ function pairwiseintersect(segments; digits=_digits(segments)) # sweepline algorithm P = eltype(first(segs)) 𝐺 = Dict{P,Vector{Int}}() - for i in eachindex(segs), j in (i + 1):length(segs) - _overlaps(starts[i], stops[i], starts[j]) || break + for i in eachindex(segs) + for j in (i + 1):length(segs) + _overlaps(starts[i], stops[i], starts[j]) || break - intersection(Segment(segs[i]), Segment(segs[j])) do I - if type(I) == Crossing || type(I) == EdgeTouching - _addintersection!(𝐺, get(I), inds[i], inds[j]; digits=digits) + intersection(Segment(segs[i]), Segment(segs[j])) do I + if type(I) == Crossing || type(I) == EdgeTouching + _addintersection!(𝐺, get(I), inds[i], inds[j]; digits=digits) + end end + continue end end (collect(keys(𝐺)), collect(values(𝐺))) From a79e599e63ab32438bca24889862b1a34f0dcc18 Mon Sep 17 00:00:00 2001 From: souma4 Date: Wed, 15 Oct 2025 14:30:29 -0600 Subject: [PATCH 15/21] no need for continue --- src/utils/intersect.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/utils/intersect.jl b/src/utils/intersect.jl index f4222efe8..8a411fb12 100644 --- a/src/utils/intersect.jl +++ b/src/utils/intersect.jl @@ -47,7 +47,6 @@ function pairwiseintersect(segments; digits=_digits(segments)) _addintersection!(𝐺, get(I), inds[i], inds[j]; digits=digits) end end - continue end end (collect(keys(𝐺)), collect(values(𝐺))) From d401576c57db1ac4a3d7e87717d185bc05e8a0d4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BAlio=20Hoffimann?= Date: Thu, 16 Oct 2025 07:15:46 -0300 Subject: [PATCH 16/21] Basic cleanup --- src/utils/intersect.jl | 57 +++++++++++++++++++++++------------------- 1 file changed, 31 insertions(+), 26 deletions(-) diff --git a/src/utils/intersect.jl b/src/utils/intersect.jl index 8a411fb12..43e9e9bad 100644 --- a/src/utils/intersect.jl +++ b/src/utils/intersect.jl @@ -21,55 +21,60 @@ tolerance of the length type of the segments. """ 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) + segs = map(segments) do s + a, b = coordround.(extrema(s), digits=digits) + a > b ? Segment(b, a) : Segment(a, b) end - starts = [CoordRefSystems.values(coords(seg[1]))[1] for seg in segs] - stops = [CoordRefSystems.values(coords(seg[2]))[1] for seg in segs] + # 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 start coordinates - inds = sortperm(starts) - starts = starts[inds] - stops = stops[inds] + # 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(first(segs)) - 𝐺 = Dict{P,Vector{Int}}() - for i in eachindex(segs) - for j in (i + 1):length(segs) - _overlaps(starts[i], stops[i], starts[j]) || break + 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 - intersection(Segment(segs[i]), Segment(segs[j])) do I + # perform more expensive intersection algorithm + intersection(segs[i], segs[j]) do I if type(I) == Crossing || type(I) == EdgeTouching - _addintersection!(𝐺, get(I), inds[i], inds[j]; digits=digits) + p = coordround(get(I); digits) + _addintersection!(D, p, inds[i], inds[j]) end end end end - (collect(keys(𝐺)), collect(values(𝐺))) -end -_overlaps(startᵢ, stopᵢ, startₖ) = (startᵢ ≤ startₖ ≤ stopᵢ) + 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) - seg = first(segments) - ℒ = lentype(seg) + s = first(segments) + ℒ = lentype(s) τ = 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]) +function _addintersection!(D, p::Point, i::Int, j::Int) + if haskey(D, p) + append!(D[p], (i, j)) + unique!(D[p]) else - 𝐺[p] = Vector{Int}([index₁, index₂]) + D[p] = Vector{Int}([i, j]) end end From 28a7f8aae4082cf8c3563aa2bd7e88f5fdbd4e69 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BAlio=20Hoffimann?= Date: Thu, 16 Oct 2025 14:31:32 -0300 Subject: [PATCH 17/21] Basic cleanup --- src/utils/intersect.jl | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/src/utils/intersect.jl b/src/utils/intersect.jl index 43e9e9bad..643a69803 100644 --- a/src/utils/intersect.jl +++ b/src/utils/intersect.jl @@ -51,12 +51,21 @@ function pairwiseintersect(segments; digits=_digits(segments)) intersection(segs[i], segs[j]) do I if type(I) == Crossing || type(I) == EdgeTouching p = coordround(get(I); digits) - _addintersection!(D, p, inds[i], inds[j]) + 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 @@ -68,13 +77,3 @@ function _digits(segments) τ = 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!(D, p::Point, i::Int, j::Int) - if haskey(D, p) - append!(D[p], (i, j)) - unique!(D[p]) - else - D[p] = Vector{Int}([i, j]) - end -end From 1939484b7580da670d0f74ae185a3466d43dd4b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BAlio=20Hoffimann?= Date: Thu, 16 Oct 2025 14:40:15 -0300 Subject: [PATCH 18/21] Improve docstring --- src/utils/intersect.jl | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/utils/intersect.jl b/src/utils/intersect.jl index 643a69803..3a57fc068 100644 --- a/src/utils/intersect.jl +++ b/src/utils/intersect.jl @@ -11,9 +11,20 @@ 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 From 9a746e200dc15878df35d3b7282a367fb9538943 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BAlio=20Hoffimann?= Date: Thu, 16 Oct 2025 14:48:25 -0300 Subject: [PATCH 19/21] Refactor tests --- test/utils.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/utils.jl b/test/utils.jl index c75a6e37f..10bcebdf3 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -147,7 +147,7 @@ end points, seginds = sortedintersection(segs) @test length(points) == n * n - 4 @test length(seginds) == n * n - 4 - @test Set(length.(seginds)) == Set([2]) + @test all(==(2), length.(seginds)) # number of intersections is invariant under rotations for θ in T(π / 6):T(π / 6):T(2π - π / 6) @@ -156,7 +156,7 @@ end θpoints, θseginds = sortedintersection(segs |> Rotate(θ)) @test length(θpoints) == n * n - 4 @test length(θseginds) == n * n - 4 - @test Set(length.(θseginds)) == Set([2]) + @test all(==(2), length.(θseginds)) end # tests coverage for when intervals don't overlap From ad49ebc1444fd145f36ddc34f37e1381177d14b3 Mon Sep 17 00:00:00 2001 From: souma4 Date: Fri, 17 Oct 2025 11:01:12 -0600 Subject: [PATCH 20/21] improved docstring and edited @inferred test --- src/utils/intersect.jl | 3 ++- test/utils.jl | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/utils/intersect.jl b/src/utils/intersect.jl index 3a57fc068..d45c0d7f7 100644 --- a/src/utils/intersect.jl +++ b/src/utils/intersect.jl @@ -6,7 +6,8 @@ pairwiseintersect(segments; [digits]) Compute pairwise intersections between n `segments` -with `digits` precision in O(n⋅log(n+k)) time using +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. diff --git a/test/utils.jl b/test/utils.jl index 10bcebdf3..c69912dfd 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -171,7 +171,7 @@ end # inference test segs = facets(cartgrid(10, 10)) - @inferred Nothing (Meshes.pairwiseintersect(segs)) + @inferred (Meshes.pairwiseintersect(segs)) end @testitem "isthreaded" setup = [Setup] begin From 4dab40ad04854fd6b4f708deedc54340bb7a14d2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BAlio=20Hoffimann?= Date: Fri, 17 Oct 2025 20:37:03 -0300 Subject: [PATCH 21/21] Adjust type of key in dictionary of intersections --- src/utils/intersect.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils/intersect.jl b/src/utils/intersect.jl index d45c0d7f7..9653c1854 100644 --- a/src/utils/intersect.jl +++ b/src/utils/intersect.jl @@ -52,7 +52,7 @@ function pairwiseintersect(segments; digits=_digits(segments)) # sweepline algorithm n = length(segs) - P = eltype(first(segs)) + P = eltype(vertices(first(segs))) D = Dict{P,Vector{Int}}() for i in 1:n for j in (i + 1):n