From 926b915ff9ea7d4b9bf6a6a9fd42627511c911b6 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 16 Apr 2025 18:51:03 -0400 Subject: [PATCH 01/32] Move utils to `utils/utils.jl` and add edge list functionality --- src/{ => utils}/utils.jl | 87 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 87 insertions(+) rename src/{ => utils}/utils.jl (60%) diff --git a/src/utils.jl b/src/utils/utils.jl similarity index 60% rename from src/utils.jl rename to src/utils/utils.jl index 540a7211e..9ea887bf7 100644 --- a/src/utils.jl +++ b/src/utils/utils.jl @@ -126,3 +126,90 @@ function _point_in_extent(p, extent::Extents.Extent) (x1, x2), (y1, y2) = extent.X, extent.Y return x1 ≤ GI.x(p) ≤ x2 && y1 ≤ GI.y(p) ≤ y2 end + +#= +# `eachedge`, `to_edgelist` + +These functions are used to decompose geometries into lists of edges. +Currently they only work on linear rings. +=# + +""" + eachedge(geom, [::Type{T}]) + +Decompose a geometry into a list of edges. +Currently only works for LineString and LinearRing. + +Returns some iterator, which yields tuples of points. Each tuple is an edge. + +It goes `(p1, p2), (p2, p3), (p3, p4), ...` etc. +""" +eachedge(geom) = eachedge(GI.trait(geom), geom, Float64) + +function eachedge(geom, ::Type{T}) where T + eachedge(GI.trait(geom), geom, T) +end + +# implementation for LineString and LinearRing +function eachedge(trait::GI.AbstractCurveTrait, geom, ::Type{T}) where T + return (_tuple_point.((GI.getpoint(geom, i), GI.getpoint(geom, i+1)), T) for i in 1:GI.npoint(geom)-1) +end + +# implementation for Polygon, MultiPolygon, MultiLineString, GeometryCollection +function eachedge(trait::GI.AbstractGeometryTrait, geom, ::Type{T}) where T + return Iterators.flatten((eachedge(r, T) for r in flatten(GI.AbstractCurveTrait, geom))) +end + +function eachedge(trait::GI.PointTrait, geom, ::Type{T}) where T + return ArgumentError("Can't get edges from points, $geom was a PointTrait.") +end + +function eachedge(trait::GI.MultiPointTrait, geom, ::Type{T}) where T + return ArgumentError("Can't get edges from MultiPoint, $geom was a MultiPointTrait.") +end + +""" + to_edgelist(geom, [::Type{T}]) + +Convert a geometry into a vector of `GI.Line` objects with attached extents. +""" +to_edgelist(geom, ::Type{T}) where T = + [_lineedge(ps, T) for ps in eachedge(geom, T)] +function to_edgelist(ext::E, geom, ::Type{T}) where {E<:Extents.Extent,T} + edges_in = eachedge(geom, T) + l1 = _lineedge(first(edges_in), T) + edges_out = typeof(l1)[] + indices = Int[] + for (i, ps) in enumerate(edges_in) + l = _lineedge(ps, T) + if Extents.intersects(ext, GI.extent(l)) + push!(edges_out, l) + push!(indices, i) + end + end + return edges_out, indices +end + +function _lineedge(ps::Tuple, ::Type{T}) where T + l = GI.Line(SVector{2,NTuple{2, T}}(ps)) # TODO: make this flexible in dimension + e = GI.extent(l) + return GI.Line(l.geom; extent=e) +end + +function lazy_edgelist(geom, ::Type{T}) where T + (_lineedge(ps, T) for ps in eachedge(geom, T)) +end + +function edge_extents(geom) + return [begin + Extents.Extent(X=extrema(GI.x, edge), Y=extrema(GI.y, edge)) + end + for edge in eachedge(geom, Float64)] +end + +function lazy_edge_extents(geom) + return (begin + Extents.Extent(X=extrema(GI.x, edge), Y=extrema(GI.y, edge)) + end + for edge in eachedge(geom, Float64)) +end From 838932a419aa453b360f6b5d9688dc81ad0c6f22 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 16 Apr 2025 22:27:13 -0400 Subject: [PATCH 02/32] Add LoopStateMachine and tests --- src/GeometryOps.jl | 10 +- .../LoopStateMachine/LoopStateMachine.jl | 106 ++++++++++++++++++ test/utils/LoopStateMachine.jl | 43 +++++++ 3 files changed, 158 insertions(+), 1 deletion(-) create mode 100644 src/utils/LoopStateMachine/LoopStateMachine.jl create mode 100644 test/utils/LoopStateMachine.jl diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index c0d169f6e..dd7d4fcba 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -36,9 +36,17 @@ const Edge{T} = Tuple{TuplePoint{T},TuplePoint{T}} where T include("types.jl") include("primitives.jl") -include("utils.jl") include("not_implemented_yet.jl") +include("utils/utils.jl") + +include("utils/LoopStateMachine/LoopStateMachine.jl") +using .LoopStateMachine + +include("utils/SpatialTreeInterface/SpatialTreeInterface.jl") +using .SpatialTreeInterface + + include("methods/angles.jl") include("methods/area.jl") include("methods/barycentric.jl") diff --git a/src/utils/LoopStateMachine/LoopStateMachine.jl b/src/utils/LoopStateMachine/LoopStateMachine.jl new file mode 100644 index 000000000..3de8239a4 --- /dev/null +++ b/src/utils/LoopStateMachine/LoopStateMachine.jl @@ -0,0 +1,106 @@ +""" + LoopStateMachine + +Utilities for returning state from functions that run inside a loop. + +This is used in e.g clipping, where we may need to break or transition states. + +The main entry point is to return an [`Action`](@ref) from a function that +is wrapped in a `@controlflow f(...)` macro in a loop. When a known `Action` +(currently, `:continue`, `:break`, `:return`, or `:full_return` actions) is returned, +it is processed by the `@controlflow` macro, which allows the function to break out of the loop +early, continue to the next iteration, or return a value, basically a way to provoke syntactic +behaviour from a function called from a inside a loop, where you do not have access to that loop. + +## Example + +```julia +``` +""" +module LoopStateMachine + +export Action, @controlflow + +import ..GeometryOps as GO + +const ALL_ACTION_DESCRIPTIONS = """ +- `:continue`: continue to the next iteration of the loop. + This is the `continue` keyword in Julia. The contents of the action are not used. +- `:break`: break out of the loop. + This is the `break` keyword in Julia. The contents of the action are not used. +- `:return`: cause the function executing the loop to return with the wrapped value. +- `:full_return`: cause the function executing the loop to return `Action(:full_return, x)`. + This is very useful to terminate recursive funtions, like tree queries terminating after you + have found a single intersecting segment. +""" + +""" + Action(name::Symbol, [x]) + +Create an `Action` with the name `name` and optional contents `x`. + +`Action`s are returned from functions wrapped in a `@controlflow` macro, which +does something based on the return value of that function if it is an `Action`. + +## Available actions + +$ALL_ACTION_DESCRIPTIONS +""" +struct Action{T} + name::Symbol + x::T +end + +Action() = Action{Nothing}(:unnamed, nothing) +Action(x::T) where T = Action{:unnamed, T}(:unnamed, x) +Action(x::Symbol) = Action(x, nothing) + +function Base.show(io::IO, action::Action{T}) where T + print(io, "Action ", action.name) + if isnothing(action.x) + print(io, "()") + else + print(io, "(", action.x, ")") + end +end + +""" + @controlflow f(...) + +Process the result of `f(...)` and return the result if it's not an `Action`(@ref LoopStateMachine.Action). + +If it is an `Action`, then process it according to the following rules, and throw an error if it's not recognized. +`:continue`, `:break`, `:return`, or `:full_return` are valid actions. + +$ALL_ACTION_DESCRIPTIONS + +!!! warning + Only use this inside a loop, otherwise you'll get a syntax error, especially if you use `:continue` or `:break`. + +## Examples +""" +macro controlflow(expr) + varname = gensym("loop-state-machine-returned-value") + return quote + $varname = $(esc(expr)) + if $varname isa Action + if $varname.name == :continue + continue + elseif $varname.name == :break + break + elseif $varname.name == :return + return $varname.x + elseif $varname.name == :full_return + return $varname + else + throw(ArgumentError("Unknown action: $($varname.name)")) + end + else + $varname + end + end +end + +# You can define more actions as you desire. + +end \ No newline at end of file diff --git a/test/utils/LoopStateMachine.jl b/test/utils/LoopStateMachine.jl new file mode 100644 index 000000000..252359745 --- /dev/null +++ b/test/utils/LoopStateMachine.jl @@ -0,0 +1,43 @@ +using Test +using GeometryOps.LoopStateMachine: @controlflow, Continue, Break + +@testset "Continue action" begin + count = 0 + f(i) = begin + count += 1 + if i == 3 + return Continue() + end + count += 1 + end + for i in 1:5 + @controlflow f(i) + end + @test count == 9 # Adds 1 for each iteration, but skips second +1 on i=3 +end + +@testset "Break action" begin + count = 0 + function f(i) + count += 1 + if i == 3 + return Break() + end + count += 1 + end + for i in 1:5 + @controlflow f(i) + end + @test count == 5 # Counts up to i=3, adding 2 for i=1,2 and 1 for i=3 +end + +@testset "Return value" begin + results = Int[] + for i in 1:3 + val = @controlflow begin + i * 2 + end + push!(results, val) + end + @test results == [2, 4, 6] +end From 34191219280f7580394aabcf2aad423842fb2cce Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 16 Apr 2025 22:31:13 -0400 Subject: [PATCH 03/32] test all actions --- test/utils/LoopStateMachine.jl | 24 +++++++++++++++++++++--- 1 file changed, 21 insertions(+), 3 deletions(-) diff --git a/test/utils/LoopStateMachine.jl b/test/utils/LoopStateMachine.jl index 252359745..4ea8dfe13 100644 --- a/test/utils/LoopStateMachine.jl +++ b/test/utils/LoopStateMachine.jl @@ -1,12 +1,12 @@ using Test -using GeometryOps.LoopStateMachine: @controlflow, Continue, Break +using GeometryOps.LoopStateMachine: @controlflow, Action @testset "Continue action" begin count = 0 f(i) = begin count += 1 if i == 3 - return Continue() + return Action(:continue) end count += 1 end @@ -21,7 +21,7 @@ end function f(i) count += 1 if i == 3 - return Break() + return Action(:break) end count += 1 end @@ -31,6 +31,24 @@ end @test count == 5 # Counts up to i=3, adding 2 for i=1,2 and 1 for i=3 end +@testset "Return action" begin + f(i) = for j in 1:3 + i == j && @controlflow Action(:return, i) + end + @test f(1) == 1 + @test f(2) == 2 + @test f(3) == 3 +end + +@testset "Full return action" begin + f(i) = for j in 1:3 + i == j && @controlflow Action(:full_return, i) + end + @test f(1) == Action(:full_return, 1) + @test f(2) == Action(:full_return, 2) + @test f(3) == Action(:full_return, 3) +end + @testset "Return value" begin results = Int[] for i in 1:3 From d6c1ffc2c3f77da5cbd33f8290baff9374a49478 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 16 Apr 2025 22:42:32 -0400 Subject: [PATCH 04/32] add SpatialTreeInterface implementation (but no tests yet, need to add those!) --- Project.toml | 2 + src/utils/SpatialTreeInterface/README.md | 47 ++++ .../SpatialTreeInterface.jl | 244 ++++++++++++++++++ 3 files changed, 293 insertions(+) create mode 100644 src/utils/SpatialTreeInterface/README.md create mode 100644 src/utils/SpatialTreeInterface/SpatialTreeInterface.jl diff --git a/Project.toml b/Project.toml index 8ee10bf44..d0eaaa7ff 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["Anshul Singhvi ", "Rafael Schouten Date: Wed, 16 Apr 2025 22:55:55 -0400 Subject: [PATCH 05/32] Define `isspatialtree` for FlatNoTree --- src/utils/SpatialTreeInterface/SpatialTreeInterface.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl index 2af8cfcd6..d1a08c756 100644 --- a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl +++ b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl @@ -6,7 +6,7 @@ import Extents import GeoInterface as GI import AbstractTrees -# public isspatialtree, getchild, nchild, child_indices_extents, node_extent +# public isspatialtree, isleaf, getchild, nchild, child_indices_extents, node_extent export query, do_query # ## Interface @@ -225,6 +225,7 @@ struct FlatNoTree{T} geometries::T end +isspatialtree(::Type{<: FlatNoTree}) = true isleaf(tree::FlatNoTree) = true # NOTE: use pairs instead of enumerate here, so that we can support From 7b1b99b9bd53114ed7ecc621cca38bdea4310f28 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 16 Apr 2025 22:56:05 -0400 Subject: [PATCH 06/32] Add basic tests for FlatNoTree --- test/utils/SpatialTreeInterface.jl | 77 ++++++++++++++++++++++++++++++ 1 file changed, 77 insertions(+) create mode 100644 test/utils/SpatialTreeInterface.jl diff --git a/test/utils/SpatialTreeInterface.jl b/test/utils/SpatialTreeInterface.jl new file mode 100644 index 000000000..5961f3981 --- /dev/null +++ b/test/utils/SpatialTreeInterface.jl @@ -0,0 +1,77 @@ +using Test + +using GeometryOps.SpatialTreeInterface +using GeometryOps.SpatialTreeInterface: isspatialtree, isleaf, getchild, nchild, child_indices_extents, node_extent +using GeometryOps.SpatialTreeInterface: query, do_query, do_dual_query +using Extents +using GeoInterface: GeoInterface as GI + +using GeometryOps.SpatialTreeInterface: FlatNoTree + +# Test FlatNoTree implementation +@testset "FlatNoTree" begin + # Test with a simple vector of extents + extents = [ + Extents.Extent(X=(0.0, 1.0), Y=(0.0, 1.0)), + Extents.Extent(X=(1.0, 2.0), Y=(1.0, 2.0)), + Extents.Extent(X=(2.0, 3.0), Y=(2.0, 3.0)) + ] + tree = FlatNoTree(extents) + + @testset "Basic interface" begin + @test isleaf(tree) + @test isspatialtree(tree) + @test isspatialtree(typeof(tree)) + end + + @testset "child_indices_extents" begin + # Test that we get the correct indices and extents + indices_extents = collect(child_indices_extents(tree)) + @test length(indices_extents) == 3 + @test indices_extents[1] == (1, extents[1]) + @test indices_extents[2] == (2, extents[2]) + @test indices_extents[3] == (3, extents[3]) + end + + @testset "Query functionality" begin + # Test query with a predicate that matches all + all_pred = x -> true + results = query(tree, all_pred) + @test sort(results) == [1, 2, 3] + + # Test query with a predicate that matches none + none_pred = x -> false + results = query(tree, none_pred) + @test isempty(results) + + # Test query with a specific extent predicate + search_extent = Extents.Extent(X=(0.5, 1.5), Y=(0.5, 1.5)) + results = query(tree, Base.Fix1(Extents.intersects, search_extent)) + @test sort(results) == [1, 2] # Should match first two extents + end + + @testset "Dual query functionality" begin + # Create two trees for dual query testing + tree1 = FlatNoTree([ + Extents.Extent(X=(0.0, 1.0), Y=(0.0, 1.0)), + Extents.Extent(X=(1.0, 2.0), Y=(1.0, 2.0)) + ]) + tree2 = FlatNoTree([ + Extents.Extent(X=(0.5, 1.5), Y=(0.5, 1.5)), + Extents.Extent(X=(1.5, 2.5), Y=(1.5, 2.5)) + ]) + + # Test dual query with a predicate that matches all + all_pred = (x, y) -> true + results = Tuple{Int, Int}[] + do_dual_query((i, j) -> push!(results, (i, j)), all_pred, tree1, tree2) + @test sort(results) == [(1,1), (1,2), (2,1), (2,2)] + + # Test dual query with a specific predicate + intersects_pred = (x, y) -> Extents.intersects(x, y) + results = Tuple{Int, Int}[] + do_dual_query((i, j) -> push!(results, (i, j)), intersects_pred, tree1, tree2) + @test sort(results) == [(1,1), (2,1), (2,2)] + end +end + From 7dd401409cfc719965997fd758c93d4581bb86b5 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 16 Apr 2025 22:56:27 -0400 Subject: [PATCH 07/32] ignore call notes since we don't want to commit them --- docs/.gitignore | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/.gitignore b/docs/.gitignore index 30c80518c..468c60bf2 100644 --- a/docs/.gitignore +++ b/docs/.gitignore @@ -2,4 +2,6 @@ node_modules/ build/ package-lock.json src/source/ -Manifest.toml \ No newline at end of file +Manifest.toml + +src/call_notes.md \ No newline at end of file From aa8ea4370a67fcc4f899cd7f6aa5b53842001363 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 16 Apr 2025 23:15:29 -0400 Subject: [PATCH 08/32] test all of LoopStateMachine --- src/utils/LoopStateMachine/LoopStateMachine.jl | 9 +++++---- test/utils/LoopStateMachine.jl | 14 ++++++++++++++ 2 files changed, 19 insertions(+), 4 deletions(-) diff --git a/src/utils/LoopStateMachine/LoopStateMachine.jl b/src/utils/LoopStateMachine/LoopStateMachine.jl index 3de8239a4..e01eea06b 100644 --- a/src/utils/LoopStateMachine/LoopStateMachine.jl +++ b/src/utils/LoopStateMachine/LoopStateMachine.jl @@ -52,15 +52,16 @@ struct Action{T} end Action() = Action{Nothing}(:unnamed, nothing) -Action(x::T) where T = Action{:unnamed, T}(:unnamed, x) +Action(x::T) where T = Action{T}(:unnamed, x) Action(x::Symbol) = Action(x, nothing) function Base.show(io::IO, action::Action{T}) where T - print(io, "Action ", action.name) + print(io, "Action") + print(io, "(:$(action.name)") if isnothing(action.x) - print(io, "()") + print(io, ")") else - print(io, "(", action.x, ")") + print(io, ", ",action.x, ")") end end diff --git a/test/utils/LoopStateMachine.jl b/test/utils/LoopStateMachine.jl index 4ea8dfe13..99cdd99a7 100644 --- a/test/utils/LoopStateMachine.jl +++ b/test/utils/LoopStateMachine.jl @@ -59,3 +59,17 @@ end end @test results == [2, 4, 6] end + +@testset "Show" begin + @test sprint(print, Action(:continue)) == "Action(:continue)" + @test sprint(print, Action(:break)) == "Action(:break)" + @test sprint(print, Action(:return, 1)) == "Action(:return, 1)" + @test sprint(print, Action(:full_return, 1)) == "Action(:full_return, 1)" +end + +@testset "Unnamed action" begin + @test sprint(print, Action()) == "Action(:unnamed)" + @test sprint(print, Action(1)) == "Action(:unnamed, 1)" + @test sprint(print, Action(:x)) == "Action(:x)" +end + From e9fb51502aff28dfa705a17503490a6c81973fe3 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 17 Apr 2025 07:44:46 -0400 Subject: [PATCH 09/32] generic tests for STI with FlatNoTree and STR --- test/utils/SpatialTreeInterface.jl | 109 ++++++++++++++++++++++++----- 1 file changed, 93 insertions(+), 16 deletions(-) diff --git a/test/utils/SpatialTreeInterface.jl b/test/utils/SpatialTreeInterface.jl index 5961f3981..33bde5a49 100644 --- a/test/utils/SpatialTreeInterface.jl +++ b/test/utils/SpatialTreeInterface.jl @@ -3,28 +3,34 @@ using Test using GeometryOps.SpatialTreeInterface using GeometryOps.SpatialTreeInterface: isspatialtree, isleaf, getchild, nchild, child_indices_extents, node_extent using GeometryOps.SpatialTreeInterface: query, do_query, do_dual_query +using GeometryOps.SpatialTreeInterface: FlatNoTree using Extents using GeoInterface: GeoInterface as GI +using SortTileRecursiveTree: STRtree -using GeometryOps.SpatialTreeInterface: FlatNoTree - -# Test FlatNoTree implementation -@testset "FlatNoTree" begin - # Test with a simple vector of extents - extents = [ - Extents.Extent(X=(0.0, 1.0), Y=(0.0, 1.0)), - Extents.Extent(X=(1.0, 2.0), Y=(1.0, 2.0)), - Extents.Extent(X=(2.0, 3.0), Y=(2.0, 3.0)) - ] - tree = FlatNoTree(extents) - +# Generic test functions for spatial trees +function test_basic_interface(TreeType) @testset "Basic interface" begin + # Create a simple tree with one extent + extents = [Extents.Extent(X=(0.0, 1.0), Y=(0.0, 1.0))] + tree = TreeType(extents) + @test isleaf(tree) @test isspatialtree(tree) @test isspatialtree(typeof(tree)) end +end +function test_child_indices_extents(TreeType) @testset "child_indices_extents" begin + # Create a tree with three extents + extents = [ + Extents.Extent(X=(0.0, 1.0), Y=(0.0, 1.0)), + Extents.Extent(X=(1.0, 2.0), Y=(1.0, 2.0)), + Extents.Extent(X=(2.0, 3.0), Y=(2.0, 3.0)) + ] + tree = TreeType(extents) + # Test that we get the correct indices and extents indices_extents = collect(child_indices_extents(tree)) @test length(indices_extents) == 3 @@ -32,8 +38,18 @@ using GeometryOps.SpatialTreeInterface: FlatNoTree @test indices_extents[2] == (2, extents[2]) @test indices_extents[3] == (3, extents[3]) end +end +function test_query_functionality(TreeType) @testset "Query functionality" begin + # Create a tree with three extents + extents = [ + Extents.Extent(X=(0.0, 1.0), Y=(0.0, 1.0)), + Extents.Extent(X=(1.0, 2.0), Y=(1.0, 2.0)), + Extents.Extent(X=(2.0, 3.0), Y=(2.0, 3.0)) + ] + tree = TreeType(extents) + # Test query with a predicate that matches all all_pred = x -> true results = query(tree, all_pred) @@ -49,14 +65,16 @@ using GeometryOps.SpatialTreeInterface: FlatNoTree results = query(tree, Base.Fix1(Extents.intersects, search_extent)) @test sort(results) == [1, 2] # Should match first two extents end +end +function test_dual_query_functionality(TreeType) @testset "Dual query functionality" begin - # Create two trees for dual query testing - tree1 = FlatNoTree([ + # Create two trees with overlapping extents + tree1 = TreeType([ Extents.Extent(X=(0.0, 1.0), Y=(0.0, 1.0)), Extents.Extent(X=(1.0, 2.0), Y=(1.0, 2.0)) ]) - tree2 = FlatNoTree([ + tree2 = TreeType([ Extents.Extent(X=(0.5, 1.5), Y=(0.5, 1.5)), Extents.Extent(X=(1.5, 2.5), Y=(1.5, 2.5)) ]) @@ -65,7 +83,7 @@ using GeometryOps.SpatialTreeInterface: FlatNoTree all_pred = (x, y) -> true results = Tuple{Int, Int}[] do_dual_query((i, j) -> push!(results, (i, j)), all_pred, tree1, tree2) - @test sort(results) == [(1,1), (1,2), (2,1), (2,2)] + @test length(results) == 4 # 2 points in tree1 * 2 points in tree2 # Test dual query with a specific predicate intersects_pred = (x, y) -> Extents.intersects(x, y) @@ -75,3 +93,62 @@ using GeometryOps.SpatialTreeInterface: FlatNoTree end end +function test_geometry_support(TreeType) + @testset "Geometry support" begin + # Create a tree with 100 points + points = tuple.(1:100, 1:100) + tree = TreeType(points) + + # Test basic interface + @test isleaf(tree) + @test isspatialtree(tree) + @test isspatialtree(typeof(tree)) + + # Test child indices and extents + indices_extents = collect(child_indices_extents(tree)) + @test length(indices_extents) == 100 + + # Check first and last points + first_idx, first_extent = indices_extents[1] + last_idx, last_extent = indices_extents[100] + + @test first_idx == 1 + @test last_idx == 100 + @test first_extent == Extents.Extent(X=(1.0, 1.0), Y=(1.0, 1.0)) + @test last_extent == Extents.Extent(X=(100.0, 100.0), Y=(100.0, 100.0)) + + # Test query functionality + all_pred = x -> true + results = query(tree, all_pred) + @test sort(results) == collect(1:100) + + none_pred = x -> false + results = query(tree, none_pred) + @test isempty(results) + + search_extent = Extents.Extent(X=(45.0, 55.0), Y=(45.0, 55.0)) + results = query(tree, Base.Fix1(Extents.intersects, search_extent)) + @test sort(results) == collect(45:55) + end +end + +# Test FlatNoTree implementation +@testset "FlatNoTree" begin + test_basic_interface(FlatNoTree) + test_child_indices_extents(FlatNoTree) + test_query_functionality(FlatNoTree) + test_dual_query_functionality(FlatNoTree) + test_geometry_support(FlatNoTree) +end + +# Test STRtree implementation +@testset "STRtree" begin + test_basic_interface(STRtree) + test_child_indices_extents(STRtree) + test_query_functionality(STRtree) + test_dual_query_functionality(STRtree) + test_geometry_support(STRtree) +end + + + From 9efde48e7c88b84a4493f9cee2dcc32e6b744b63 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 17 Apr 2025 07:44:54 -0400 Subject: [PATCH 10/32] actually include tests --- src/utils/SpatialTreeInterface/SpatialTreeInterface.jl | 1 + test/runtests.jl | 5 +++++ 2 files changed, 6 insertions(+) diff --git a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl index d1a08c756..ac93a7482 100644 --- a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl +++ b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl @@ -8,6 +8,7 @@ import AbstractTrees # public isspatialtree, isleaf, getchild, nchild, child_indices_extents, node_extent export query, do_query +export FlatNoTree # ## Interface # Interface definition for spatial tree types. diff --git a/test/runtests.jl b/test/runtests.jl index 3d5465c96..829b2e545 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,6 +11,11 @@ end @safetestset "Types" begin include("types.jl") end @safetestset "Primitives" begin include("primitives.jl") end + +# Utils +@safetestset "LoopStateMachine" begin include("utils/LoopStateMachine.jl") end +@safetestset "SpatialTreeInterface" begin include("utils/SpatialTreeInterface.jl") end + # Methods @safetestset "Angles" begin include("methods/angles.jl") end @safetestset "Area" begin include("methods/area.jl") end From e5ad556c80a059ac97853af19801a30bbbcd2b7e Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 17 Apr 2025 07:49:18 -0400 Subject: [PATCH 11/32] implement `isspatialtree` for STR types --- src/utils/SpatialTreeInterface/SpatialTreeInterface.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl index ac93a7482..74913eec7 100644 --- a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl +++ b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl @@ -200,18 +200,20 @@ end using SortTileRecursiveTree: STRtree, STRNode, STRLeafNode +isspatialtree(::Type{<: STRtree}) = true nchild(tree::STRtree) = nchild(tree.rootnode) getchild(tree::STRtree) = getchild(tree.rootnode) getchild(tree::STRtree, i) = getchild(tree.rootnode, i) isleaf(tree::STRtree) = isleaf(tree.rootnode) child_indices_extents(tree::STRtree) = child_indices_extents(tree.rootnode) - +isspatialtree(::Type{<: STRNode}) = true nchild(node::STRNode) = length(node.children) getchild(node::STRNode) = node.children getchild(node::STRNode, i) = node.children[i] isleaf(node::STRNode) = false # STRNodes are not leaves by definition +isspatialtree(::Type{<: STRLeafNode}) = true isleaf(node::STRLeafNode) = true child_indices_extents(node::STRLeafNode) = zip(node.indices, node.extents) From 4ccf67bd790eeb3fd0781c46fdcec9924bfad7d2 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 17 Apr 2025 07:51:30 -0400 Subject: [PATCH 12/32] fix tests --- test/utils/SpatialTreeInterface.jl | 17 +---------------- 1 file changed, 1 insertion(+), 16 deletions(-) diff --git a/test/utils/SpatialTreeInterface.jl b/test/utils/SpatialTreeInterface.jl index 33bde5a49..50b5f25d2 100644 --- a/test/utils/SpatialTreeInterface.jl +++ b/test/utils/SpatialTreeInterface.jl @@ -14,8 +14,7 @@ function test_basic_interface(TreeType) # Create a simple tree with one extent extents = [Extents.Extent(X=(0.0, 1.0), Y=(0.0, 1.0))] tree = TreeType(extents) - - @test isleaf(tree) + @test isspatialtree(tree) @test isspatialtree(typeof(tree)) end @@ -100,23 +99,9 @@ function test_geometry_support(TreeType) tree = TreeType(points) # Test basic interface - @test isleaf(tree) @test isspatialtree(tree) @test isspatialtree(typeof(tree)) - # Test child indices and extents - indices_extents = collect(child_indices_extents(tree)) - @test length(indices_extents) == 100 - - # Check first and last points - first_idx, first_extent = indices_extents[1] - last_idx, last_extent = indices_extents[100] - - @test first_idx == 1 - @test last_idx == 100 - @test first_extent == Extents.Extent(X=(1.0, 1.0), Y=(1.0, 1.0)) - @test last_extent == Extents.Extent(X=(100.0, 100.0), Y=(100.0, 100.0)) - # Test query functionality all_pred = x -> true results = query(tree, all_pred) From 0900db2a1f9f87ffc4f38fd6271e225094d5a7dc Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 17 Apr 2025 09:29:33 -0400 Subject: [PATCH 13/32] Apply suggestions from code review Co-authored-by: Rafael Schouten --- src/GeometryOps.jl | 6 ++---- src/utils/utils.jl | 5 ----- 2 files changed, 2 insertions(+), 9 deletions(-) diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index dd7d4fcba..db7efd2ee 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -39,12 +39,10 @@ include("primitives.jl") include("not_implemented_yet.jl") include("utils/utils.jl") - include("utils/LoopStateMachine/LoopStateMachine.jl") -using .LoopStateMachine - include("utils/SpatialTreeInterface/SpatialTreeInterface.jl") -using .SpatialTreeInterface + +using .LoopStateMachine, .SpatialTreeInterface include("methods/angles.jl") diff --git a/src/utils/utils.jl b/src/utils/utils.jl index 9ea887bf7..fd882768c 100644 --- a/src/utils/utils.jl +++ b/src/utils/utils.jl @@ -145,25 +145,20 @@ Returns some iterator, which yields tuples of points. Each tuple is an edge. It goes `(p1, p2), (p2, p3), (p3, p4), ...` etc. """ eachedge(geom) = eachedge(GI.trait(geom), geom, Float64) - function eachedge(geom, ::Type{T}) where T eachedge(GI.trait(geom), geom, T) end - # implementation for LineString and LinearRing function eachedge(trait::GI.AbstractCurveTrait, geom, ::Type{T}) where T return (_tuple_point.((GI.getpoint(geom, i), GI.getpoint(geom, i+1)), T) for i in 1:GI.npoint(geom)-1) end - # implementation for Polygon, MultiPolygon, MultiLineString, GeometryCollection function eachedge(trait::GI.AbstractGeometryTrait, geom, ::Type{T}) where T return Iterators.flatten((eachedge(r, T) for r in flatten(GI.AbstractCurveTrait, geom))) end - function eachedge(trait::GI.PointTrait, geom, ::Type{T}) where T return ArgumentError("Can't get edges from points, $geom was a PointTrait.") end - function eachedge(trait::GI.MultiPointTrait, geom, ::Type{T}) where T return ArgumentError("Can't get edges from MultiPoint, $geom was a MultiPointTrait.") end From 06c76a93b831ebf83cb87c1cc1e8215ce108a1cc Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 17 Apr 2025 10:29:24 -0400 Subject: [PATCH 14/32] Implement and test a custom exception type for unknown actions in LSM --- src/utils/LoopStateMachine/LoopStateMachine.jl | 18 +++++++++++++++++- test/utils/LoopStateMachine.jl | 15 +++++++++++++++ 2 files changed, 32 insertions(+), 1 deletion(-) diff --git a/src/utils/LoopStateMachine/LoopStateMachine.jl b/src/utils/LoopStateMachine/LoopStateMachine.jl index e01eea06b..77b62d450 100644 --- a/src/utils/LoopStateMachine/LoopStateMachine.jl +++ b/src/utils/LoopStateMachine/LoopStateMachine.jl @@ -65,6 +65,21 @@ function Base.show(io::IO, action::Action{T}) where T end end +struct UnrecognizedActionException <: Base.Exception + name::Symbol +end + +function Base.showerror(io::IO, e::UnrecognizedActionException) + print(io, "Unrecognized action: ") + printstyled(io, e.name; color = :red, bold = true) + println(io, ".") + println(io, "Valid actions are:") + println(io, ALL_ACTION_DESCRIPTIONS) +end + +# We exclude the macro definition from code coverage computations, +# because I know it's tested but Codecov doesn't seem to think so. +# COV_EXCL_START """ @controlflow f(...) @@ -94,13 +109,14 @@ macro controlflow(expr) elseif $varname.name == :full_return return $varname else - throw(ArgumentError("Unknown action: $($varname.name)")) + throw(UnrecognizedActionException($varname.name)) end else $varname end end end +# COV_EXCL_STOP # You can define more actions as you desire. diff --git a/test/utils/LoopStateMachine.jl b/test/utils/LoopStateMachine.jl index 99cdd99a7..30d4979c8 100644 --- a/test/utils/LoopStateMachine.jl +++ b/test/utils/LoopStateMachine.jl @@ -1,4 +1,5 @@ using Test +import GeometryOps.LoopStateMachine using GeometryOps.LoopStateMachine: @controlflow, Action @testset "Continue action" begin @@ -73,3 +74,17 @@ end @test sprint(print, Action(:x)) == "Action(:x)" end +@testset "Unknown action" begin + f(i) = Action(:unknown_action, i) + @test_throws LoopStateMachine.UnrecognizedActionException begin + for i in 1:3 + @controlflow f(i) + end + end + + @test_throws ":continue" begin + for i in 1:3 + @controlflow f(i) + end + end +end From 728befcf3649f43dc51506081240e4efda41cede Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 17 Apr 2025 10:46:02 -0400 Subject: [PATCH 15/32] Factor out SpatialTreeInterface into separate files Maybe this is a bad idea, single file was nice for readability. But we can keep going. --- .../SpatialTreeInterface.jl | 211 +----------------- .../depth_first_search.jl | 29 +++ .../dual_depth_first_search.jl | 56 +++++ .../SpatialTreeInterface/implementations.jl | 55 +++++ src/utils/SpatialTreeInterface/interface.jl | 85 +++++++ 5 files changed, 236 insertions(+), 200 deletions(-) create mode 100644 src/utils/SpatialTreeInterface/depth_first_search.jl create mode 100644 src/utils/SpatialTreeInterface/dual_depth_first_search.jl create mode 100644 src/utils/SpatialTreeInterface/implementations.jl create mode 100644 src/utils/SpatialTreeInterface/interface.jl diff --git a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl index 74913eec7..2cc15b6c6 100644 --- a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl +++ b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl @@ -10,118 +10,19 @@ import AbstractTrees export query, do_query export FlatNoTree -# ## Interface -# Interface definition for spatial tree types. -# There is no abstract supertype here since it's impossible to enforce, -# but we do have a few methods that are common to all spatial tree types. +# The spatial tree interface and its implementations are defined here. +include("interface.jl") +include("implementations.jl") -""" - isspatialtree(tree)::Bool - -Return true if the object is a spatial tree, false otherwise. - -## Implementation notes - -For type stability, if your spatial tree type is `MyTree`, you should define -`isspatialtree(::Type{MyTree}) = true`, and `isspatialtree(::MyTree)` will forward -to that method automatically. -""" -isspatialtree(::T) where T = isspatialtree(T) -isspatialtree(::Type{<: Any}) = false - - -""" - getchild(node) - -Return an iterator over all the children of a node. -This may be materialized if necessary or available, -but can also be lazy (like a generator). -""" -getchild(node) = AbstractTrees.children(node) - -""" - getchild(node, i) - -Return the `i`-th child of a node. -""" -getchild(node, i) = getchild(node)[i] - -""" - nchild(node) +# Here we have some algorithms that use the spatial tree interface. +# The first file holds a single depth-first search, i.e., a single-tree query. +include("depth_first_search.jl") -Return the number of children of a node. -""" -nchild(node) = length(getchild(node)) - -""" - isleaf(node) - -Return true if the node is a leaf node, i.e., there are no "children" below it. -[`getchild`](@ref) should still work on leaf nodes, though, returning an iterator over the extents stored in the node - and similarly for `getnodes.` -""" -isleaf(node) = error("isleaf is not implemented for node type $(typeof(node))") - -""" - child_indices_extents(node) - -Return an iterator over the indices and extents of the children of a node. - -Each value of the iterator should take the form `(i, extent)`. - -This can only be invoked on leaf nodes! -""" -function child_indices_extents(node) - return zip(1:nchild(node), getchild(node)) -end - -""" - node_extent(node) - -Return the extent like object of the node. -Falls back to `GI.extent` by default, which falls back -to `Extents.extent`. - -Generally, defining `Extents.extent(node)` is sufficient here, and you -won't need to define this - -The reason we don't use that directly is to give users of this interface -a way to define bounding boxes that are not extents, like spherical caps -and other such things. -""" -node_extent(node) = GI.extent(node) - -# ## Query functions -# These are generic functions that work with any spatial tree type that implements the interface. - - -""" - do_query(f, predicate, tree) - -Call `f(i)` for each index `i` in the tree that satisfies `predicate(extent(i))`. - -This is generic to anything that implements the SpatialTreeInterface, particularly the methods -[`isleaf`](@ref), [`getchild`](@ref), and [`child_extents`](@ref). -""" -function do_query(f::F, predicate::P, node::N) where {F, P, N} - if isleaf(node) - for (i, leaf_geometry_extent) in child_indices_extents(node) - if predicate(leaf_geometry_extent) - @controlflow f(i) - end - end - else - for child in getchild(node) - if predicate(node_extent(child)) - @controlflow do_query(f, predicate, child) - end - end - end -end -function do_query(predicate, node) - a = Int[] - do_query(Base.Fix1(push!, a), predicate, node) - return a -end +# The second file holds a dual depth-first search, i.e., a dual-tree query. +# This iterates over two trees simultaneously, and is substantially more efficient +# than two separate single-tree queries since it can prune branches in tandem as it +# descends into the trees. +include("dual_depth_first_search.jl") """ @@ -155,94 +56,4 @@ sanitize_predicate(::GI.AbstractTrait, pred) = sanitize_predicate(GI.extent(pred sanitize_predicate(pred::Extents.Extent) = Base.Fix1(Extents.intersects, pred) -""" - do_dual_query(f, predicate, node1, node2) - -Call `f(i1, i2)` for each index `i1` in `node1` and `i2` in `node2` that satisfies `predicate(extent(i1), extent(i2))`. - -This is generic to anything that implements the SpatialTreeInterface, particularly the methods -[`isleaf`](@ref), [`getchild`](@ref), and [`child_extents`](@ref). -""" -function do_dual_query(f::F, predicate::P, node1::N1, node2::N2) where {F, P, N1, N2} - if isleaf(node1) && isleaf(node2) - # both nodes are leaves, so we can just iterate over the indices and extents - for (i1, extent1) in child_indices_extents(node1) - for (i2, extent2) in child_indices_extents(node2) - if predicate(extent1, extent2) - @controlflow f(i1, i2) - end - end - end - elseif isleaf(node1) # node2 is not a leaf, node1 is - recurse further into node2 - for child in getchild(node2) - if predicate(node_extent(node1), node_extent(child)) - @controlflow do_dual_query(f, predicate, node1, child) - end - end - elseif isleaf(node2) # node1 is not a leaf, node2 is - recurse further into node1 - for child in getchild(node1) - if predicate(node_extent(child), node_extent(node2)) - @controlflow do_dual_query(f, predicate, child, node2) - end - end - else # neither node is a leaf, recurse into both children - for child1 in getchild(node1) - for child2 in getchild(node2) - if predicate(node_extent(child1), node_extent(child2)) - @controlflow do_dual_query(f, predicate, child1, child2) - end - end - end - end -end - -# Finally, here's a sample implementation of the interface for STRtrees - -using SortTileRecursiveTree: STRtree, STRNode, STRLeafNode - -isspatialtree(::Type{<: STRtree}) = true -nchild(tree::STRtree) = nchild(tree.rootnode) -getchild(tree::STRtree) = getchild(tree.rootnode) -getchild(tree::STRtree, i) = getchild(tree.rootnode, i) -isleaf(tree::STRtree) = isleaf(tree.rootnode) -child_indices_extents(tree::STRtree) = child_indices_extents(tree.rootnode) - -isspatialtree(::Type{<: STRNode}) = true -nchild(node::STRNode) = length(node.children) -getchild(node::STRNode) = node.children -getchild(node::STRNode, i) = node.children[i] -isleaf(node::STRNode) = false # STRNodes are not leaves by definition - -isspatialtree(::Type{<: STRLeafNode}) = true -isleaf(node::STRLeafNode) = true -child_indices_extents(node::STRLeafNode) = zip(node.indices, node.extents) - - -""" - FlatNoTree(iterable_of_geoms_or_extents) - -Represents a flat collection with no tree structure, i.e., a brute force search. -This is cost free, so particularly useful when you don't want to build a tree! -""" -struct FlatNoTree{T} - geometries::T -end - -isspatialtree(::Type{<: FlatNoTree}) = true -isleaf(tree::FlatNoTree) = true - -# NOTE: use pairs instead of enumerate here, so that we can support -# iterators or collections that define custom `pairs` methods. -# This includes things like filtered extent lists, for example, -# so we can perform extent thinning with no allocations. -function child_indices_extents(tree::FlatNoTree{T}) where T - # This test only applies at compile time and should be optimized away in any case. - # And we can use multiple dispatch to override anyway, but it should be cost free I think. - if applicable(Base.keys, T) - return ((i, GI.extent(obj)) for (i, obj) in pairs(tree.geometries)) - else - return ((i, GI.extent(obj)) for (i, obj) in enumerate(tree.geometries)) - end -end - end # module SpatialTreeInterface \ No newline at end of file diff --git a/src/utils/SpatialTreeInterface/depth_first_search.jl b/src/utils/SpatialTreeInterface/depth_first_search.jl new file mode 100644 index 000000000..dc1be59b3 --- /dev/null +++ b/src/utils/SpatialTreeInterface/depth_first_search.jl @@ -0,0 +1,29 @@ + +""" + depth_first_search(f, predicate, tree) + +Call `f(i)` for each index `i` in the tree that satisfies `predicate(extent(i))`. + +This is generic to anything that implements the SpatialTreeInterface, particularly the methods +[`isleaf`](@ref), [`getchild`](@ref), and [`child_extents`](@ref). +""" +function depth_first_search(f::F, predicate::P, node::N) where {F, P, N} + if isleaf(node) + for (i, leaf_geometry_extent) in child_indices_extents(node) + if predicate(leaf_geometry_extent) + @controlflow f(i) + end + end + else + for child in getchild(node) + if predicate(node_extent(child)) + @controlflow depth_first_search(f, predicate, child) + end + end + end +end +function depth_first_search(predicate, node) + a = Int[] + depth_first_search(Base.Fix1(push!, a), predicate, node) + return a +end diff --git a/src/utils/SpatialTreeInterface/dual_depth_first_search.jl b/src/utils/SpatialTreeInterface/dual_depth_first_search.jl new file mode 100644 index 000000000..2c63d0aa1 --- /dev/null +++ b/src/utils/SpatialTreeInterface/dual_depth_first_search.jl @@ -0,0 +1,56 @@ +""" + dual_depth_first_search(f, predicate, tree1, tree2) + +Executes a dual depth-first search over two trees, descending into the children of +nodes `i` and `j` when `predicate(node_extent(i), node_extent(j))` is true, +and pruning that branch when `predicate(node_extent(i), node_extent(j))` is false. + +Finally, calls `f(i1, i2)` for each leaf-level index `i1::Int` in `tree1` and `i2::Int` in `tree2` +that satisfies `predicate(extent(i1), extent(i2))`. + +Here, `f(i1::Int, i2::Int)` may be any function that takes two integers as arguments. +It may optionally return an [`Action`](@ref LoopStateMachine.Action) to alter the control +flow of the `Action(:full_return, true)` to return `Action(:full_return, true)` from this +function and break out of the recursion. + +This is generic to anything that implements the SpatialTreeInterface, particularly the methods +[`isleaf`](@ref), [`getchild`](@ref), and [`child_indices_extents`](@ref). + +## Examples + +```julia +using NaturalEarth, +``` +""" +function dual_depth_first_search(f::F, predicate::P, node1::N1, node2::N2) where {F, P, N1, N2} + if isleaf(node1) && isleaf(node2) + # both nodes are leaves, so we can just iterate over the indices and extents + for (i1, extent1) in child_indices_extents(node1) + for (i2, extent2) in child_indices_extents(node2) + if predicate(extent1, extent2) + @controlflow f(i1, i2) + end + end + end + elseif isleaf(node1) # node2 is not a leaf, node1 is - recurse further into node2 + for child in getchild(node2) + if predicate(node_extent(node1), node_extent(child)) + @controlflow dual_depth_first_search(f, predicate, node1, child) + end + end + elseif isleaf(node2) # node1 is not a leaf, node2 is - recurse further into node1 + for child in getchild(node1) + if predicate(node_extent(child), node_extent(node2)) + @controlflow dual_depth_first_search(f, predicate, child, node2) + end + end + else # neither node is a leaf, recurse into both children + for child1 in getchild(node1) + for child2 in getchild(node2) + if predicate(node_extent(child1), node_extent(child2)) + @controlflow dual_depth_first_search(f, predicate, child1, child2) + end + end + end + end +end diff --git a/src/utils/SpatialTreeInterface/implementations.jl b/src/utils/SpatialTreeInterface/implementations.jl new file mode 100644 index 000000000..a2df1e5dc --- /dev/null +++ b/src/utils/SpatialTreeInterface/implementations.jl @@ -0,0 +1,55 @@ +# # Interface implementations +# Below are some basic implementations of the interface, +# for STRTree and a "no-tree" implementation that is a flat list of extents. + +using SortTileRecursiveTree: STRtree, STRNode, STRLeafNode + +# ## SortTileRecursiveTree + +isspatialtree(::Type{<: STRtree}) = true +node_extent(tree::STRtree) = node_extent(tree.rootnode) +nchild(tree::STRtree) = nchild(tree.rootnode) +getchild(tree::STRtree) = getchild(tree.rootnode) +getchild(tree::STRtree, i) = getchild(tree.rootnode, i) +isleaf(tree::STRtree) = isleaf(tree.rootnode) +child_indices_extents(tree::STRtree) = child_indices_extents(tree.rootnode) + +isspatialtree(::Type{<: STRNode}) = true +node_extent(node::STRNode) = node.extent +nchild(node::STRNode) = length(node.children) +getchild(node::STRNode) = node.children +getchild(node::STRNode, i) = node.children[i] +isleaf(node::STRNode) = false # STRNodes are not leaves by definition + +isspatialtree(::Type{<: STRLeafNode}) = true +node_extent(node::STRLeafNode) = node.extent +isleaf(node::STRLeafNode) = true +child_indices_extents(node::STRLeafNode) = zip(node.indices, node.extents) + +# ## FlatNoTree +""" + FlatNoTree(iterable_of_geoms_or_extents) + +Represents a flat collection with no tree structure, i.e., a brute force search. +This is cost free, so particularly useful when you don't want to build a tree! +""" +struct FlatNoTree{T} + geometries::T +end + +isspatialtree(::Type{<: FlatNoTree}) = true +isleaf(tree::FlatNoTree) = true + +# NOTE: use pairs instead of enumerate here, so that we can support +# iterators or collections that define custom `pairs` methods. +# This includes things like filtered extent lists, for example, +# so we can perform extent thinning with no allocations. +function child_indices_extents(tree::FlatNoTree{T}) where T + # This test only applies at compile time and should be optimized away in any case. + # And we can use multiple dispatch to override anyway, but it should be cost free I think. + if applicable(Base.keys, T) + return ((i, GI.extent(obj)) for (i, obj) in pairs(tree.geometries)) + else + return ((i, GI.extent(obj)) for (i, obj) in enumerate(tree.geometries)) + end +end \ No newline at end of file diff --git a/src/utils/SpatialTreeInterface/interface.jl b/src/utils/SpatialTreeInterface/interface.jl new file mode 100644 index 000000000..5e0bf08ab --- /dev/null +++ b/src/utils/SpatialTreeInterface/interface.jl @@ -0,0 +1,85 @@ +# # Interface +# Interface definition for spatial tree types. +# There is no abstract supertype here since it's impossible to enforce, +# but we do have a few methods that are common to all spatial tree types. + +""" + isspatialtree(tree)::Bool + +Return true if the object is a spatial tree, false otherwise. + +## Implementation notes + +For type stability, if your spatial tree type is `MyTree`, you should define +`isspatialtree(::Type{MyTree}) = true`, and `isspatialtree(::MyTree)` will forward +to that method automatically. +""" +isspatialtree(::T) where T = isspatialtree(T) +isspatialtree(::Type{<: Any}) = false + + +""" + getchild(node) + getchild(node, i) + +Accessor function to get the children of a node. + +If invoked as `getchild(node)`, return an iterator over all the children of a node. +This may be lazy, like a `Base.Generator`, or it may be materialized. + +If invoked as `getchild(node, i)`, return the `i`-th child of a node. +""" +function getchild end + +getchild(node) = AbstractTrees.children(node) + +""" + getchild(node, i) + +Return the `i`-th child of a node. +""" +getchild(node, i) = getchild(node)[i] + +""" + nchild(node) + +Return the number of children of a node. +""" +nchild(node) = length(getchild(node)) + +""" + isleaf(node) + +Return true if the node is a leaf node, i.e., there are no "children" below it. +[`getchild`](@ref) should still work on leaf nodes, though, returning an iterator over the extents stored in the node - and similarly for `getnodes.` +""" +isleaf(node) = error("isleaf is not implemented for node type $(typeof(node))") + +""" + child_indices_extents(node) + +Return an iterator over the indices and extents of the children of a node. + +Each value of the iterator should take the form `(i, extent)`. + +This can only be invoked on leaf nodes! +""" +function child_indices_extents(node) + return zip(1:nchild(node), getchild(node)) +end + +""" + node_extent(node) + +Return the extent like object of the node. +Falls back to `GI.extent` by default, which falls back +to `Extents.extent`. + +Generally, defining `Extents.extent(node)` is sufficient here, and you +won't need to define this + +The reason we don't use that directly is to give users of this interface +a way to define bounding boxes that are not extents, like spherical caps +and other such things. +""" +node_extent(node) = GI.extent(node) From bb7620723a60422e536370e166b8fd75766fa38a Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 17 Apr 2025 15:20:38 -0400 Subject: [PATCH 16/32] Fix tests --- test/utils/SpatialTreeInterface.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/utils/SpatialTreeInterface.jl b/test/utils/SpatialTreeInterface.jl index 50b5f25d2..94f09aaa0 100644 --- a/test/utils/SpatialTreeInterface.jl +++ b/test/utils/SpatialTreeInterface.jl @@ -2,7 +2,7 @@ using Test using GeometryOps.SpatialTreeInterface using GeometryOps.SpatialTreeInterface: isspatialtree, isleaf, getchild, nchild, child_indices_extents, node_extent -using GeometryOps.SpatialTreeInterface: query, do_query, do_dual_query +using GeometryOps.SpatialTreeInterface: query, depth_first_search, dual_depth_first_search using GeometryOps.SpatialTreeInterface: FlatNoTree using Extents using GeoInterface: GeoInterface as GI @@ -81,13 +81,13 @@ function test_dual_query_functionality(TreeType) # Test dual query with a predicate that matches all all_pred = (x, y) -> true results = Tuple{Int, Int}[] - do_dual_query((i, j) -> push!(results, (i, j)), all_pred, tree1, tree2) + dual_depth_first_search((i, j) -> push!(results, (i, j)), all_pred, tree1, tree2) @test length(results) == 4 # 2 points in tree1 * 2 points in tree2 # Test dual query with a specific predicate intersects_pred = (x, y) -> Extents.intersects(x, y) results = Tuple{Int, Int}[] - do_dual_query((i, j) -> push!(results, (i, j)), intersects_pred, tree1, tree2) + dual_depth_first_search((i, j) -> push!(results, (i, j)), intersects_pred, tree1, tree2) @test sort(results) == [(1,1), (2,1), (2,2)] end end From 8f4de31ca61c5dc3038596a64f16951e883e7f80 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 17 Apr 2025 15:21:03 -0400 Subject: [PATCH 17/32] no really --- src/utils/SpatialTreeInterface/SpatialTreeInterface.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl index 2cc15b6c6..d306d75b4 100644 --- a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl +++ b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl @@ -32,7 +32,7 @@ Return a sorted list of indices of the tree that satisfy the predicate. """ function query(tree, predicate) a = Int[] - do_query(Base.Fix1(push!, a), sanitize_predicate(predicate), tree) + depth_first_search(Base.Fix1(push!, a), sanitize_predicate(predicate), tree) return sort!(a) end From d5a49de75873f1c5f40fb149eb4ce384b308f208 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 19 Apr 2025 20:01:37 -0400 Subject: [PATCH 18/32] add a test for multilevel dual tree query --- Project.toml | 3 +- test/utils/SpatialTreeInterface.jl | 89 ++++++++++++++++++++++++++---- 2 files changed, 80 insertions(+), 12 deletions(-) diff --git a/Project.toml b/Project.toml index d0eaaa7ff..8aad7e7d4 100644 --- a/Project.toml +++ b/Project.toml @@ -67,6 +67,7 @@ JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb" NaturalEarth = "436b0209-26ab-4e65-94a9-6526d86fea76" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" +Polylabel = "49a44318-e865-4b63-9842-695152d634c1" Proj = "c94c279d-25a6-4763-9509-64d165bea63e" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689" @@ -76,4 +77,4 @@ TGGeometry = "d7e755d2-3c95-4bcf-9b3c-79ab1a78647b" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["ArchGDAL", "CoordinateTransformations", "DataFrames", "Distributions", "DimensionalData", "Downloads", "FlexiJoins", "GeoJSON", "Proj", "JLD2", "LibGEOS", "Random", "Rasters", "NaturalEarth", "OffsetArrays", "SafeTestsets", "Shapefile", "TGGeometry", "Test"] +test = ["ArchGDAL", "CoordinateTransformations", "DataFrames", "Distributions", "DimensionalData", "Downloads", "FlexiJoins", "GeoJSON", "Proj", "JLD2", "LibGEOS", "Random", "Rasters", "NaturalEarth", "OffsetArrays", "Polylabel", "SafeTestsets", "Shapefile", "TGGeometry", "Test"] diff --git a/test/utils/SpatialTreeInterface.jl b/test/utils/SpatialTreeInterface.jl index 94f09aaa0..dc7e4245d 100644 --- a/test/utils/SpatialTreeInterface.jl +++ b/test/utils/SpatialTreeInterface.jl @@ -7,6 +7,8 @@ using GeometryOps.SpatialTreeInterface: FlatNoTree using Extents using GeoInterface: GeoInterface as GI using SortTileRecursiveTree: STRtree +using NaturalEarth +using Polylabel # Generic test functions for spatial trees function test_basic_interface(TreeType) @@ -67,7 +69,7 @@ function test_query_functionality(TreeType) end function test_dual_query_functionality(TreeType) - @testset "Dual query functionality" begin + @testset "Dual query functionality - simple" begin # Create two trees with overlapping extents tree1 = TreeType([ Extents.Extent(X=(0.0, 1.0), Y=(0.0, 1.0)), @@ -90,6 +92,65 @@ function test_dual_query_functionality(TreeType) dual_depth_first_search((i, j) -> push!(results, (i, j)), intersects_pred, tree1, tree2) @test sort(results) == [(1,1), (2,1), (2,2)] end + + @testset "Dual query functionality - every country's polylabel against every country" begin + + # Note that this is a perfectly balanced tree query - we don't yet have a test for unbalanced + # trees (but could easily add one, e.g. by getting polylabels of admin-1 or admin-2 regions) + # from Natural Earth, or by using GADM across many countries. + + all_countries = NaturalEarth.naturalearth("admin_0_countries", 10) + all_adm0_a3 = all_countries.ADM0_A3 + all_geoms = all_countries.geometry + # US minor outlying islands - bug in Polylabel.jl + # A lot of small geoms have this issue, that there will be an error from the queue + # because the cell exists in the queue already. + # Not sure what to do about it, I don't want to check containment every time... + deleteat!(all_adm0_a3, 205) + deleteat!(all_geoms, 205) + + geom_tree = TreeType(all_geoms) + + polylabels = [Polylabel.polylabel(geom; rtol = 0.019) for geom in all_geoms] + polylabel_tree = TreeType(polylabels) + + found_countries = falses(length(polylabels)) + + dual_depth_first_search(Extents.intersects, geom_tree, polylabel_tree) do i, j + if i == j + found_countries[i] = true + end + end + + @test all(found_countries) + end +end + +function test_find_point_in_all_countries(TreeType) + all_countries = NaturalEarth.naturalearth("admin_0_countries", 10) + tree = TreeType(all_countries.geometry) + + ber = (13.4050, 52.5200) # Berlin + nyc = (-74.0060, 40.7128) # New York City + sin = (103.8198, 1.3521) # Singapore + + @testset "locate points using query" begin + @testset let point = ber, name = "Berlin" + # Test Berlin (should be in Germany) + results = query(tree, point) + @test any(i -> all_countries.ADM0_A3[i] == "DEU", results) + end + @testset let point = nyc, name = "New York City" + # Test NYC (should be in USA) + results = query(tree, point) + @test any(i -> all_countries.ADM0_A3[i] == "USA", results) + end + @testset let point = sin, name = "Singapore" + # Test Singapore + results = query(tree, point) + @test any(i -> all_countries.ADM0_A3[i] == "SGP", results) + end + end end function test_geometry_support(TreeType) @@ -119,20 +180,26 @@ end # Test FlatNoTree implementation @testset "FlatNoTree" begin - test_basic_interface(FlatNoTree) - test_child_indices_extents(FlatNoTree) - test_query_functionality(FlatNoTree) - test_dual_query_functionality(FlatNoTree) - test_geometry_support(FlatNoTree) + @testset let TreeType = FlatNoTree + test_basic_interface(TreeType) + test_child_indices_extents(TreeType) + test_query_functionality(TreeType) + test_dual_query_functionality(TreeType) + test_geometry_support(TreeType) + test_find_point_in_all_countries(TreeType) + end end # Test STRtree implementation @testset "STRtree" begin - test_basic_interface(STRtree) - test_child_indices_extents(STRtree) - test_query_functionality(STRtree) - test_dual_query_functionality(STRtree) - test_geometry_support(STRtree) + @testset let TreeType = STRtree + test_basic_interface(TreeType) + test_child_indices_extents(TreeType) + test_query_functionality(TreeType) + test_dual_query_functionality(TreeType) + test_geometry_support(TreeType) + test_find_point_in_all_countries(TreeType) + end end From 06727f1bc23f6f916d3febcd7e689de9cc625357 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 19 Apr 2025 20:32:31 -0400 Subject: [PATCH 19/32] Fix up the tests --- test/utils/SpatialTreeInterface.jl | 150 ++++++++++++++++------------- 1 file changed, 84 insertions(+), 66 deletions(-) diff --git a/test/utils/SpatialTreeInterface.jl b/test/utils/SpatialTreeInterface.jl index dc7e4245d..c2d38e37d 100644 --- a/test/utils/SpatialTreeInterface.jl +++ b/test/utils/SpatialTreeInterface.jl @@ -92,37 +92,47 @@ function test_dual_query_functionality(TreeType) dual_depth_first_search((i, j) -> push!(results, (i, j)), intersects_pred, tree1, tree2) @test sort(results) == [(1,1), (2,1), (2,2)] end + @testset "Dual tree query with many boundingboxes" begin + xs = 1:100 + ys = 1:100 + extent_grid = [Extents.Extent(X=(x, x+1), Y=(y, y+1)) for x in xs, y in ys] |> vec + point_grid = [(x + 0.5, y + 0.5) for x in xs, y in ys] |> vec - @testset "Dual query functionality - every country's polylabel against every country" begin - - # Note that this is a perfectly balanced tree query - we don't yet have a test for unbalanced - # trees (but could easily add one, e.g. by getting polylabels of admin-1 or admin-2 regions) - # from Natural Earth, or by using GADM across many countries. - - all_countries = NaturalEarth.naturalearth("admin_0_countries", 10) - all_adm0_a3 = all_countries.ADM0_A3 - all_geoms = all_countries.geometry - # US minor outlying islands - bug in Polylabel.jl - # A lot of small geoms have this issue, that there will be an error from the queue - # because the cell exists in the queue already. - # Not sure what to do about it, I don't want to check containment every time... - deleteat!(all_adm0_a3, 205) - deleteat!(all_geoms, 205) - - geom_tree = TreeType(all_geoms) - - polylabels = [Polylabel.polylabel(geom; rtol = 0.019) for geom in all_geoms] - polylabel_tree = TreeType(polylabels) - - found_countries = falses(length(polylabels)) + extent_tree = TreeType(extent_grid) + point_tree = TreeType(point_grid) - dual_depth_first_search(Extents.intersects, geom_tree, polylabel_tree) do i, j + found_everything = falses(length(extent_grid)) + dual_depth_first_search(Extents.intersects, extent_tree, point_tree) do i, j if i == j - found_countries[i] = true + found_everything[i] = true end end + @test all(found_everything) + end +end - @test all(found_countries) +function test_geometry_support(TreeType) + @testset "Geometry support" begin + # Create a tree with 100 points + points = tuple.(1:100, 1:100) + tree = TreeType(points) + + # Test basic interface + @test isspatialtree(tree) + @test isspatialtree(typeof(tree)) + + # Test query functionality + all_pred = x -> true + results = query(tree, all_pred) + @test sort(results) == collect(1:100) + + none_pred = x -> false + results = query(tree, none_pred) + @test isempty(results) + + search_extent = Extents.Extent(X=(45.0, 55.0), Y=(45.0, 55.0)) + results = query(tree, Base.Fix1(Extents.intersects, search_extent)) + @test sort(results) == collect(45:55) end end @@ -153,54 +163,62 @@ function test_find_point_in_all_countries(TreeType) end end -function test_geometry_support(TreeType) - @testset "Geometry support" begin - # Create a tree with 100 points - points = tuple.(1:100, 1:100) - tree = TreeType(points) - - # Test basic interface - @test isspatialtree(tree) - @test isspatialtree(typeof(tree)) - - # Test query functionality - all_pred = x -> true - results = query(tree, all_pred) - @test sort(results) == collect(1:100) - - none_pred = x -> false - results = query(tree, none_pred) - @test isempty(results) - - search_extent = Extents.Extent(X=(45.0, 55.0), Y=(45.0, 55.0)) - results = query(tree, Base.Fix1(Extents.intersects, search_extent)) - @test sort(results) == collect(45:55) - end -end - # Test FlatNoTree implementation @testset "FlatNoTree" begin - @testset let TreeType = FlatNoTree - test_basic_interface(TreeType) - test_child_indices_extents(TreeType) - test_query_functionality(TreeType) - test_dual_query_functionality(TreeType) - test_geometry_support(TreeType) - test_find_point_in_all_countries(TreeType) - end + test_basic_interface(FlatNoTree) + test_child_indices_extents(FlatNoTree) + test_query_functionality(FlatNoTree) + test_dual_query_functionality(FlatNoTree) + test_geometry_support(FlatNoTree) + test_find_point_in_all_countries(FlatNoTree) end # Test STRtree implementation @testset "STRtree" begin - @testset let TreeType = STRtree - test_basic_interface(TreeType) - test_child_indices_extents(TreeType) - test_query_functionality(TreeType) - test_dual_query_functionality(TreeType) - test_geometry_support(TreeType) - test_find_point_in_all_countries(TreeType) - end + test_basic_interface(STRtree) + test_child_indices_extents(STRtree) + test_query_functionality(STRtree) + test_dual_query_functionality(STRtree) + test_geometry_support(STRtree) + test_find_point_in_all_countries(STRtree) end +# This testset is not used because Polylabel.jl has some issues. + +#= + + + @testset "Dual query functionality - every country's polylabel against every country" begin + + # Note that this is a perfectly balanced tree query - we don't yet have a test for unbalanced + # trees (but could easily add one, e.g. by getting polylabels of admin-1 or admin-2 regions) + # from Natural Earth, or by using GADM across many countries. + + all_countries = NaturalEarth.naturalearth("admin_0_countries", 10) + all_adm0_a3 = all_countries.ADM0_A3 + all_geoms = all_countries.geometry + # US minor outlying islands - bug in Polylabel.jl + # A lot of small geoms have this issue, that there will be an error from the queue + # because the cell exists in the queue already. + # Not sure what to do about it, I don't want to check containment every time... + deleteat!(all_adm0_a3, 205) + deleteat!(all_geoms, 205) + + geom_tree = TreeType(all_geoms) + + polylabels = [Polylabel.polylabel(geom; rtol = 0.019) for geom in all_geoms] + polylabel_tree = TreeType(polylabels) + + found_countries = falses(length(polylabels)) + + dual_depth_first_search(Extents.intersects, geom_tree, polylabel_tree) do i, j + if i == j + found_countries[i] = true + end + end + + @test all(found_countries) + end +=# \ No newline at end of file From 91ff5c224d19a755da7894309907972de72cebec Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 3 Mar 2025 20:43:16 -0500 Subject: [PATCH 20/32] Describe a FosterHormannClipping algorithm Not sure of name (too verbose, maybe Foster()?). Intersection accelerators will come into the picture later (next PR) but we need them now to instantiate the struct. --- src/methods/clipping/clipping_processor.jl | 53 ++++++++++++++++++++++ 1 file changed, 53 insertions(+) diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl index 556f140a6..b18220bd6 100644 --- a/src/methods/clipping/clipping_processor.jl +++ b/src/methods/clipping/clipping_processor.jl @@ -1,6 +1,59 @@ # # Polygon clipping helpers # This file contains the shared helper functions for the polygon clipping functionalities. +# This file specifically defines helpers for the Foster-Hormann clipping algorithm. + + +""" + abstract type IntersectionAccelerator + +The abstract supertype for all intersection accelerator types. + +The idea is that these speed up the edge-edge intersection checking process, +perhaps at the cost of memory. + +The naive case is [`NestedLoop`](@ref), which is just a nested loop, running in O(n*m) time. + +Then we have [`SingleSTRtree`](@ref), which is a single STRtree, running in O(n*log(m)) time. + +Then we have [`DoubleSTRtree`](@ref), which is am simultaneous double-tree traversal of two STRtrees. + +Finally, we have [`AutoAccelerator`](@ref), which is an automatic accelerator that chooses the best +accelerator based on the size of the input polygons. This gets materialized in build_a_list for now. +""" +abstract type IntersectionAccelerator end +struct NestedLoop <: IntersectionAccelerator end +struct SingleSTRtree <: IntersectionAccelerator end +struct DoubleSTRtree <: IntersectionAccelerator end +struct SingleNaturalTree <: IntersectionAccelerator end +struct DoubleNaturalTree <: IntersectionAccelerator end +struct ThinnedDoubleNaturalTree <: IntersectionAccelerator end +struct AutoAccelerator <: IntersectionAccelerator end + +""" + FosterHormannClipping{M <: Manifold, A <: Union{Nothing, Accelerator}} <: GeometryOpsCore.Algorithm{M} + +A type that represents the Foster-Hormann clipping algorithm. + +# Arguments +- `manifold::M`: The manifold on which the algorithm operates. +- `accelerator::A`: The accelerator to use for the algorithm. Can be `nothing` for automatic choice, or a custom accelerator. +""" +struct FosterHormannClipping{M <: Manifold, A <: IntersectionAccelerator} <: GeometryOpsCore.Algorithm{M} + manifold::M + accelerator::A + # TODO: add exact flag + # TODO: should exact flag be in the type domain? +end + + +FosterHormannClipping(; manifold::Manifold = Planar(), accelerator = nothing) = FosterHormannClipping(manifold, isnothing(accelerator) ? AutoAccelerator() : accelerator) +FosterHormannClipping(manifold::Manifold, accelerator::Union{Nothing, IntersectionAccelerator} = nothing) = FosterHormannClipping(manifold, isnothing(accelerator) ? AutoAccelerator() : accelerator) +FosterHormannClipping(accelerator::Union{Nothing, IntersectionAccelerator}) = FosterHormannClipping(Planar(), isnothing(accelerator) ? AutoAccelerator() : accelerator) +# special case for spherical / geodesic manifolds +# since they can't use STRtrees (because those don't work on the sphere) +FosterHormannClipping(manifold::Union{Spherical, Geodesic}, accelerator::Union{Nothing, IntersectionAccelerator} = nothing) = FosterHormannClipping(manifold, isnothing(accelerator) ? NestedLoop() : (accelerator isa AutoAccelerator ? NestedLoop() : accelerator)) + # This enum defines which side of an edge a point is on @enum PointEdgeSide left=1 right=2 unknown=3 From c81aef234a8a8929e548f46e7995c763e35bb22d Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 3 Mar 2025 20:44:50 -0500 Subject: [PATCH 21/32] Create a TracingError type that stores all context Better than printing out since we can (A) check on it and (B) access the polygons as well as the a_list and b_list in Julia. I have a few plot recipes for them but not sure if they are worth putting in a Makie extension...they would also be inaccessible there. Maybe in GeoMakie? Then we can make full use of customizations etc. --- src/methods/clipping/clipping_processor.jl | 39 ++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl index b18220bd6..d8df29e8a 100644 --- a/src/methods/clipping/clipping_processor.jl +++ b/src/methods/clipping/clipping_processor.jl @@ -90,6 +90,45 @@ PolyNode(node::PolyNode{T}; # Checks equality of two PolyNodes by backing point value, fractional value, and intersection status equals(pn1::PolyNode, pn2::PolyNode) = pn1.point == pn2.point && pn1.inter == pn2.inter && pn1.fracs == pn2.fracs +Base.:(==)(pn1::PolyNode, pn2::PolyNode) = equals(pn1, pn2) + +# Finally, we define a nice error type for when the clipping tracing algorithm hits every point in a polygon. +# This stores the polygons, the a_list, and the b_list, and the a_idx_list. +# allowing the user to understand what happened and why. +""" + TracingError{T1, T2} <: Exception + +An error that is thrown when the clipping tracing algorithm fails somehow. +This is a bug in the algorithm, and should be reported. + +The polygons are contained in the exception object, accessible by try-catch or as `err` in the REPL. +""" +struct TracingError{T1, T2, T} <: Exception + message::String + poly_a::T1 + poly_b::T2 + a_list::Vector{PolyNode{T}} + b_list::Vector{PolyNode{T}} + a_idx_list::Vector{Int} +end + +function Base.showerror(io::IO, e::TracingError{T1, T2}) where {T1, T2} + print(io, "TracingError: ") + println(io, e.message) + println(io, "Please open an issue with the polygons contained in this error object.") + println(io) + if max(GI.npoint(e.poly_a), GI.npoint(e.poly_b)) < 10 + println(io, "Polygon A:") + println(io, GI.coordinates(e.poly_a)) + println(io) + println(io, "Polygon B:") + println(io, GI.coordinates(e.poly_b)) + else + println(io, "The polygons are contained in the exception object, accessible by try-catch or as `err` in the REPL.") + end +end + + #= _build_ab_list(::Type{T}, poly_a, poly_b, delay_cross_f, delay_bounce_f; exact) -> From 7bcceb2392063b4de3f848ad6058e6c00ef695f8 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 3 Mar 2025 20:45:21 -0500 Subject: [PATCH 22/32] Pass algorithm all the way through `cut` This is a shorter change list as a sampler, to show what needs to happen. --- src/methods/clipping/cut.jl | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/methods/clipping/cut.jl b/src/methods/clipping/cut.jl index 3393f3efa..bffe5763f 100644 --- a/src/methods/clipping/cut.jl +++ b/src/methods/clipping/cut.jl @@ -58,21 +58,24 @@ GI.coordinates.(cut_polys) [[[5.0, 0.0], [10.0, 0.0], [10.0, 10.0], [5.0, 10.0], [5.0, 0.0]]] ``` """ -cut(geom, line, ::Type{T} = Float64) where {T <: AbstractFloat} = - _cut(T, GI.trait(geom), geom, GI.trait(line), line; exact = True()) +cut(geom, line, ::Type{T} = Float64) where {T <: AbstractFloat} = cut(FosterHormannClipping(), geom, line, T) +cut(m::Manifold, geom, line, ::Type{T} = Float64) where {T <: AbstractFloat} = cut(FosterHormannClipping(m), geom, line, T) + +cut(alg::FosterHormannClipping{M, A}, geom, line, ::Type{T} = Float64) where {T <: AbstractFloat, M, A} = + _cut(alg, T, GI.trait(geom), geom, GI.trait(line), line; exact = True()) #= Cut a given polygon by given line. Add polygon holes back into resulting pieces if there are any holes. =# -function _cut(::Type{T}, ::GI.PolygonTrait, poly, ::GI.LineTrait, line; exact) where T +function _cut(alg::FosterHormannClipping{M, A}, ::Type{T}, ::GI.PolygonTrait, poly, ::GI.LineTrait, line; exact) where {T, M, A} ext_poly = GI.getexterior(poly) - poly_list, intr_list = _build_a_list(T, ext_poly, line; exact) + poly_list, intr_list = _build_a_list(alg, T, ext_poly, line; exact) n_intr_pts = length(intr_list) # If an impossible number of intersection points, return original polygon if n_intr_pts < 2 || isodd(n_intr_pts) return [tuples(poly)] end # Cut polygon by line - cut_coords = _cut(T, ext_poly, line, poly_list, intr_list, n_intr_pts; exact) + cut_coords = _cut(alg, T, ext_poly, line, poly_list, intr_list, n_intr_pts; exact) # Close coords and create polygons for c in cut_coords push!(c, c[1]) From ab52b171808e3d5927dc6956d53d3ac969da9bf4 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 3 Mar 2025 20:47:01 -0500 Subject: [PATCH 23/32] Pass algorithm through all functions in clipping_processor --- src/methods/clipping/clipping_processor.jl | 64 ++++++++++++---------- 1 file changed, 34 insertions(+), 30 deletions(-) diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl index d8df29e8a..3dcd547cd 100644 --- a/src/methods/clipping/clipping_processor.jl +++ b/src/methods/clipping/clipping_processor.jl @@ -140,20 +140,20 @@ returns are the fully updated vectors of PolyNodes that represent the rings 'pol 'poly_b', respectively. This function also returns 'a_idx_list', which at its "ith" index stores the index in 'a_list' at which the "ith" intersection point lies. =# -function _build_ab_list(::Type{T}, poly_a, poly_b, delay_cross_f::F1, delay_bounce_f::F2; exact) where {T, F1, F2} +function _build_ab_list(alg::FosterHormannClipping, ::Type{T}, poly_a, poly_b, delay_cross_f::F1, delay_bounce_f::F2; exact) where {T, F1, F2} # Make a list for nodes of each polygon - a_list, a_idx_list, n_b_intrs = _build_a_list(T, poly_a, poly_b; exact) - b_list = _build_b_list(T, a_idx_list, a_list, n_b_intrs, poly_b) + a_list, a_idx_list, n_b_intrs = _build_a_list(alg, T, poly_a, poly_b; exact) + b_list = _build_b_list(alg, T, a_idx_list, a_list, n_b_intrs, poly_b) # Flag crossings - _classify_crossing!(T, a_list, b_list; exact) + _classify_crossing!(alg, T, a_list, b_list; exact) # Flag the entry and exits - _flag_ent_exit!(T, GI.LinearRingTrait(), poly_b, a_list, delay_cross_f, Base.Fix2(delay_bounce_f, true); exact) - _flag_ent_exit!(T, GI.LinearRingTrait(), poly_a, b_list, delay_cross_f, Base.Fix2(delay_bounce_f, false); exact) + _flag_ent_exit!(alg, T, GI.LinearRingTrait(), poly_b, a_list, delay_cross_f, Base.Fix2(delay_bounce_f, true); exact) + _flag_ent_exit!(alg, T, GI.LinearRingTrait(), poly_a, b_list, delay_cross_f, Base.Fix2(delay_bounce_f, false); exact) # Set node indices and filter a_idx_list to just crossing points - _index_crossing_intrs!(a_list, b_list, a_idx_list) + _index_crossing_intrs!(alg, a_list, b_list, a_idx_list) return a_list, b_list, a_idx_list end @@ -172,7 +172,7 @@ not update the entry and exit flags for a_list. The a_idx_list is a list of the indices of intersection points in a_list. The value at index i of a_idx_list is the location in a_list where the ith intersection point lies. =# -function _build_a_list(::Type{T}, poly_a, poly_b; exact) where T +function _build_a_list(alg::FosterHormannClipping{M, A}, ::Type{T}, poly_a, poly_b; exact) where {T, M, A} n_a_edges = _nedge(poly_a) a_list = PolyNode{T}[] # list of points in poly_a sizehint!(a_list, n_a_edges) @@ -273,7 +273,7 @@ is needed for clipping using the Greiner-Hormann clipping algorithm. Note: after calling this function, b_list is not fully updated. The entry/exit flags still need to be updated. However, the neighbor value in a_list is now updated. =# -function _build_b_list(::Type{T}, a_idx_list, a_list, n_b_intrs, poly_b) where T +function _build_b_list(alg::FosterHormannClipping{M, A}, ::Type{T}, a_idx_list, a_list, n_b_intrs, poly_b) where {T, M, A} # Sort intersection points by insertion order in b_list sort!(a_idx_list, by = x-> a_list[x].neighbor + a_list[x].fracs[2]) # Initialize needed values and lists @@ -333,7 +333,7 @@ chain is crossing and delayed otherwise and all middle points are marked as boun Additionally, the start and end points of the chain are marked as endpoints using the endpoints field. =# -function _classify_crossing!(::Type{T}, a_list, b_list; exact) where T +function _classify_crossing!(alg::FosterHormannClipping{M, A}, ::Type{T}, a_list, b_list; exact) where {T, M, A} napts = length(a_list) nbpts = length(b_list) # start centered on last point @@ -357,7 +357,7 @@ function _classify_crossing!(::Type{T}, a_list, b_list; exact) where T a_next_is_b_prev = a_next.inter && equals(a_next, b_prev) a_next_is_b_next = a_next.inter && equals(a_next, b_next) # determine which side of a segments the p points are on - b_prev_side, b_next_side = _get_sides(b_prev, b_next, a_prev, curr_pt, a_next, + b_prev_side, b_next_side = _get_sides(#=TODO: alg.manifold, =#b_prev, b_next, a_prev, curr_pt, a_next, i, j, a_list, b_list; exact) # no sides overlap if !a_prev_is_b_prev && !a_prev_is_b_next && !a_next_is_b_prev && !a_next_is_b_next @@ -499,7 +499,7 @@ all points are intersection points, find the first element that either is the en or a crossing point that isn't in a chain. Then take the midpoint of this point and the next point in the list and perform the in/out check. If none of these points exist, return a `next_idx` of `nothing`. =# -function _pt_off_edge_status(::Type{T}, pt_list, poly, npts; exact) where T +function _pt_off_edge_status(alg::FosterHormannClipping{M, A}, ::Type{T}, pt_list, poly, npts; exact) where {T, M, A} start_idx, is_non_intr_pt = findfirst(_is_not_intr, pt_list), true if isnothing(start_idx) start_idx, is_non_intr_pt = findfirst(_next_edge_off, pt_list), false @@ -511,7 +511,7 @@ function _pt_off_edge_status(::Type{T}, pt_list, poly, npts; exact) where T else (pt_list[start_idx].point .+ pt_list[next_idx].point) ./ 2 end - start_status = !_point_filled_curve_orientation(start_pt, poly; in = true, on = false, out = false, exact) + start_status = !_point_filled_curve_orientation(alg.manifold, start_pt, poly; in = true, on = false, out = false, exact) return next_idx, start_status end # Check if a PolyNode is an intersection point @@ -537,10 +537,10 @@ bounce will be the same. Used for clipping polygons by other polygons. =# -function _flag_ent_exit!(::Type{T}, ::GI.LinearRingTrait, poly, pt_list, delay_cross_f, delay_bounce_f; exact) where T +function _flag_ent_exit!(alg::FosterHormannClipping{M, A}, ::Type{T}, ::GI.LinearRingTrait, poly, pt_list, delay_cross_f, delay_bounce_f; exact) where {T, M, A} npts = length(pt_list) # Find starting index if there is one - next_idx, status = _pt_off_edge_status(T, pt_list, poly, npts; exact) + next_idx, status = _pt_off_edge_status(alg, T, pt_list, poly, npts; exact) isnothing(next_idx) && return start_idx = next_idx - 1 # Loop over points and mark entry and exit status @@ -560,7 +560,7 @@ function _flag_ent_exit!(::Type{T}, ::GI.LinearRingTrait, poly, pt_list, delay_c else # delayed bouncing next_idx = ii < npts ? (ii + 1) : 1 next_val = (curr_pt.point .+ pt_list[next_idx].point) ./ 2 - pt_in_poly = _point_filled_curve_orientation(next_val, poly; in = true, on = false, out = false, exact) + pt_in_poly = _point_filled_curve_orientation(alg.manifold, next_val, poly; in = true, on = false, out = false, exact) #= start and end crossing status are the same and depend on if adjacent edges of pt_list are within poly =# start_crossing = delay_bounce_f(pt_in_poly) @@ -588,8 +588,8 @@ returns false. Used for cutting polygons by lines. Assumes that the first point is outside of the polygon and not on an edge. =# -function _flag_ent_exit!(::GI.LineTrait, poly, pt_list; exact) - status = !_point_filled_curve_orientation(pt_list[1].point, poly; in = true, on = false, out = false, exact) +function _flag_ent_exit!(alg::FosterHormannClipping{M, A}, ::GI.LineTrait, poly, pt_list; exact) where {M, A} + status = !_point_filled_curve_orientation(#=TODO: alg.manifold=#pt_list[1].point, poly; in = true, on = false, out = false, exact) # Loop over points and mark entry and exit status for (ii, curr_pt) in enumerate(pt_list) if curr_pt.crossing @@ -602,7 +602,7 @@ end #= Filters a_idx_list to just include crossing points and sets the index of all crossing points (which element they correspond to within a_idx_list). =# -function _index_crossing_intrs!(a_list, b_list, a_idx_list) +function _index_crossing_intrs!(alg::FosterHormannClipping{M, A}, a_list, b_list, a_idx_list) where {M, A} filter!(x -> a_list[x].crossing, a_idx_list) for (i, a_idx) in enumerate(a_idx_list) curr_node = a_list[a_idx] @@ -631,7 +631,7 @@ A list of GeoInterface polygons is returned from this function. Note: `poly_a` and `poly_b` are temporary inputs used for debugging and can be removed eventually. =# -function _trace_polynodes(::Type{T}, a_list, b_list, a_idx_list, f_step, poly_a, poly_b) where T +function _trace_polynodes(alg::FosterHormannClipping{M, A}, ::Type{T}, a_list, b_list, a_idx_list, f_step, poly_a, poly_b) where {T, M, A} n_a_pts, n_b_pts = length(a_list), length(b_list) total_pts = n_a_pts + n_b_pts n_cross_pts = length(a_idx_list) @@ -707,7 +707,7 @@ or they are separate polygons with no intersection (other than an edge or point) Return two booleans that represent if a is inside b (potentially with shared edges / points) and visa versa if b is inside of a. =# -function _find_non_cross_orientation(a_list, b_list, a_poly, b_poly; exact) +function _find_non_cross_orientation(m::M, a_list, b_list, a_poly, b_poly; exact) where {M <: Manifold} non_intr_a_idx = findfirst(x -> !x.inter, a_list) non_intr_b_idx = findfirst(x -> !x.inter, b_list) #= Determine if non-intersection point is in or outside of polygon - if there isn't A @@ -721,6 +721,9 @@ function _find_non_cross_orientation(a_list, b_list, a_poly, b_poly; exact) return a_in_b, b_in_a end +_find_non_cross_orientation(alg::FosterHormannClipping{M}, a_list, b_list, a_poly, b_poly; exact) where {M <: Manifold} = + _find_non_cross_orientation(alg.manifold, a_list, b_list, a_poly, b_poly; exact) + #= _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator, remove_poly_idx; exact) @@ -728,7 +731,7 @@ The holes specified by the hole iterator are added to the polygons in the return If this creates more polygons, they are added to the end of the list. If this removes polygons, they are removed from the list =# -function _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator, remove_poly_idx; exact) where T +function _add_holes_to_polys!(alg::FosterHormannClipping{M, A}, ::Type{T}, return_polys, hole_iterator, remove_poly_idx; exact) where {T, M, A} n_polys = length(return_polys) remove_hole_idx = Int[] # Remove set of holes from all polygons @@ -741,9 +744,9 @@ function _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator, remove_pol curr_poly = return_polys[j] remove_poly_idx[j] && continue curr_poly_ext = GI.nhole(curr_poly) > 0 ? GI.Polygon(StaticArrays.SVector(GI.getexterior(curr_poly))) : curr_poly - in_ext, on_ext, out_ext = _line_polygon_interactions(curr_hole, curr_poly_ext; exact, closed_line = true) + in_ext, on_ext, out_ext = _line_polygon_interactions(#=TODO: alg.manifold=#curr_hole, curr_poly_ext; exact, closed_line = true) if in_ext # hole is at least partially within the polygon's exterior - new_hole, new_hole_poly, n_new_pieces = _combine_holes!(T, curr_hole, curr_poly, return_polys, remove_hole_idx) + new_hole, new_hole_poly, n_new_pieces = _combine_holes!(alg, T, curr_hole, curr_poly, return_polys, remove_hole_idx) if n_new_pieces > 0 append!(remove_poly_idx, falses(n_new_pieces)) n_new_per_poly += n_new_pieces @@ -751,7 +754,7 @@ function _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator, remove_pol if !on_ext && !out_ext # hole is completely within exterior push!(curr_poly.geom, new_hole) else # hole is partially within and outside of polygon's exterior - new_polys = difference(curr_poly_ext, new_hole_poly, T; target=GI.PolygonTrait()) + new_polys = difference(alg, curr_poly_ext, new_hole_poly, T; target=GI.PolygonTrait()) n_new_polys = length(new_polys) - 1 # replace original curr_poly.geom[1] = GI.getexterior(new_polys[1]) @@ -763,7 +766,7 @@ function _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator, remove_pol end end # polygon is completely within hole - elseif coveredby(curr_poly_ext, GI.Polygon(StaticArrays.SVector(curr_hole))) + elseif coveredby(#=TODO: alg.manifold=#curr_poly_ext, GI.Polygon(StaticArrays.SVector(curr_hole))) remove_poly_idx[j] = true end end @@ -788,16 +791,16 @@ are in the "main" polygon or in one of these new pieces and moved accordingly. If the holes don't touch or curr_poly has no holes, then new_hole is returned without any changes. =# -function _combine_holes!(::Type{T}, new_hole, curr_poly, return_polys, remove_hole_idx) where T +function _combine_holes!(alg::FosterHormannClipping{M, A}, ::Type{T}, new_hole, curr_poly, return_polys, remove_hole_idx) where {T, M, A} n_new_polys = 0 empty!(remove_hole_idx) new_hole_poly = GI.Polygon(StaticArrays.SVector(new_hole)) # Combine any existing holes in curr_poly with new hole for (k, old_hole) in enumerate(GI.gethole(curr_poly)) old_hole_poly = GI.Polygon(StaticArrays.SVector(old_hole)) - if intersects(new_hole_poly, old_hole_poly) + if intersects(#=TODO: alg.manifold=#new_hole_poly, old_hole_poly) # If the holes intersect, combine them into a bigger hole - hole_union = union(new_hole_poly, old_hole_poly, T; target = GI.PolygonTrait())[1] + hole_union = union(alg, new_hole_poly, old_hole_poly, T; target = GI.PolygonTrait())[1] push!(remove_hole_idx, k + 1) new_hole = GI.getexterior(hole_union) new_hole_poly = GI.Polygon(StaticArrays.SVector(new_hole)) @@ -826,7 +829,7 @@ end #= Remove collinear edge points, other than the first and last edge vertex, to simplify polygon - including both the exterior ring and any holes=# -function _remove_collinear_points!(polys, remove_idx, poly_a, poly_b) +function _remove_collinear_points!(alg::FosterHormannClipping{M, A}, polys, remove_idx, poly_a, poly_b) where {M, A} for (i, poly) in Iterators.reverse(enumerate(polys)) for (j, ring) in Iterators.reverse(enumerate(GI.getring(poly))) n = length(ring.geom) @@ -844,6 +847,7 @@ function _remove_collinear_points!(polys, remove_idx, poly_a, poly_b) else p3 = p # check if p2 is approximately on the edge formed by p1 and p3 - remove if so + # TODO: make this manifold aware if Predicates.orient(p1, p2, p3; exact = False()) == 0 remove_idx[i - 1] = true end From dd514b177417ea66faac570ce998b283d4156665 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 3 Mar 2025 20:47:18 -0500 Subject: [PATCH 24/32] Pass manifold through point_filled_curve_orientation in planar --- src/methods/geom_relations/geom_geom_processors.jl | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/methods/geom_relations/geom_geom_processors.jl b/src/methods/geom_relations/geom_geom_processors.jl index 7e4658a7a..8563fb174 100644 --- a/src/methods/geom_relations/geom_geom_processors.jl +++ b/src/methods/geom_relations/geom_geom_processors.jl @@ -491,7 +491,7 @@ of the curve if it didn't return 'on'. See paper for more information on cases denoted in comments. =# function _point_filled_curve_orientation( - point, curve; + ::Planar, point, curve; in::T = point_in, on::T = point_on, out::T = point_out, exact, ) where {T} x, y = GI.x(point), GI.y(point) @@ -526,6 +526,11 @@ function _point_filled_curve_orientation( return iseven(k) ? out : in end +_point_filled_curve_orientation( + point, curve; + in::T = point_in, on::T = point_on, out::T = point_out, exact, +) where {T} = _point_filled_curve_orientation(Planar(), point, curve; in, on, out, exact) + #= Determines the types of interactions of a line with a filled-in curve. By filled-in curve, I am referring to the exterior ring of a poylgon, for example. From caf8800ccd4479ef77f624db1542196c86b1eb12 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 3 Mar 2025 20:47:59 -0500 Subject: [PATCH 25/32] Pass algorithm through difference + union --- src/methods/clipping/cut.jl | 6 +- src/methods/clipping/difference.jl | 56 ++++++++++----- src/methods/clipping/union.jl | 109 ++++++++++++++++++++++------- 3 files changed, 124 insertions(+), 47 deletions(-) diff --git a/src/methods/clipping/cut.jl b/src/methods/clipping/cut.jl index bffe5763f..c3427731a 100644 --- a/src/methods/clipping/cut.jl +++ b/src/methods/clipping/cut.jl @@ -83,7 +83,7 @@ function _cut(alg::FosterHormannClipping{M, A}, ::Type{T}, ::GI.PolygonTrait, po cut_polys = [GI.Polygon([c]) for c in cut_coords] # Add original polygon holes back in remove_idx = falses(length(cut_polys)) - _add_holes_to_polys!(T, cut_polys, GI.gethole(poly), remove_idx; exact) + _add_holes_to_polys!(alg, T, cut_polys, GI.gethole(poly), remove_idx; exact) return cut_polys end @@ -100,10 +100,10 @@ end of cut geometry in Vector{Vector{Tuple}} format. Note: degenerate cases where intersection points are vertices do not work right now. =# -function _cut(::Type{T}, geom, line, geom_list, intr_list, n_intr_pts; exact) where T +function _cut(alg::FosterHormannClipping{M, A}, ::Type{T}, geom, line, geom_list, intr_list, n_intr_pts; exact) where {T, M, A} # Sort and categorize the intersection points sort!(intr_list, by = x -> geom_list[x].fracs[2]) - _flag_ent_exit!(GI.LineTrait(), line, geom_list; exact) + _flag_ent_exit!(alg, GI.LineTrait(), line, geom_list; exact) # Add first point to output list return_coords = [[geom_list[1].point]] cross_backs = [(T(Inf),T(Inf))] diff --git a/src/methods/clipping/difference.jl b/src/methods/clipping/difference.jl index ca7a37078..19c72491b 100644 --- a/src/methods/clipping/difference.jl +++ b/src/methods/clipping/difference.jl @@ -33,19 +33,24 @@ GI.coordinates.(diff_poly) ``` """ function difference( - geom_a, geom_b, ::Type{T} = Float64; target=nothing, kwargs..., + alg::FosterHormannClipping, geom_a, geom_b, ::Type{T} = Float64; target=nothing, kwargs..., ) where {T<:AbstractFloat} return _difference( - TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b; + alg, TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b; exact = True(), kwargs..., ) end +# fallback definitions +difference(geom_a, geom_b, ::Type{T} = Float64; target=nothing, kwargs...) where T = difference(FosterHormannClipping(Planar()), geom_a, geom_b, T; target, kwargs...) +# if manifold but no algorithm - assume FosterHormannClipping with provided manifold. +difference(m::Manifold, geom_a, geom_b, ::Type{T} = Float64; target=nothing, kwargs...) where T = difference(FosterHormannClipping(m), geom_a, geom_b, T; target, kwargs...) + #= The 'difference' function returns the difference of two polygons as a list of polygons. The algorithm to determine the difference was adapted from "Efficient clipping of efficient polygons," by Greiner and Hormann (1998). DOI: https://doi.org/10.1145/274363.274364 =# function _difference( - ::TraitTarget{GI.PolygonTrait}, ::Type{T}, + alg::FosterHormannClipping, target::TraitTarget{GI.PolygonTrait}, ::Type{T}, ::GI.PolygonTrait, poly_a, ::GI.PolygonTrait, poly_b; exact, kwargs... @@ -54,11 +59,11 @@ function _difference( ext_a = GI.getexterior(poly_a) ext_b = GI.getexterior(poly_b) # Find the difference of the exterior of the polygons - a_list, b_list, a_idx_list = _build_ab_list(T, ext_a, ext_b, _diff_delay_cross_f, _diff_delay_bounce_f; exact) - polys = _trace_polynodes(T, a_list, b_list, a_idx_list, _diff_step, poly_a, poly_b) + a_list, b_list, a_idx_list = _build_ab_list(alg, T, ext_a, ext_b, _diff_delay_cross_f, _diff_delay_bounce_f; exact) + polys = _trace_polynodes(alg, T, a_list, b_list, a_idx_list, _diff_step, poly_a, poly_b) # if no crossing points, determine if either poly is inside of the other if isempty(polys) - a_in_b, b_in_a = _find_non_cross_orientation(a_list, b_list, ext_a, ext_b; exact) + a_in_b, b_in_a = _find_non_cross_orientation(alg.manifold, a_list, b_list, ext_a, ext_b; exact) # add case for if they polygons are the same (all intersection points!) # add a find_first check to find first non-inter poly! if b_in_a && !a_in_b # b in a and can't be the same polygon @@ -72,19 +77,19 @@ function _difference( remove_idx = falses(length(polys)) # If the original polygons had holes, take that into account. if GI.nhole(poly_a) != 0 - _add_holes_to_polys!(T, polys, GI.gethole(poly_a), remove_idx; exact) + _add_holes_to_polys!(alg, T, polys, GI.gethole(poly_a), remove_idx; exact) end if GI.nhole(poly_b) != 0 for hole in GI.gethole(poly_b) hole_poly = GI.Polygon(StaticArrays.SVector(hole)) - new_polys = intersection(hole_poly, poly_a, T; target = GI.PolygonTrait) + new_polys = intersection(alg, hole_poly, poly_a, T; target = GI.PolygonTrait) if length(new_polys) > 0 append!(polys, new_polys) end end end # Remove unneeded collinear points on same edge - _remove_collinear_points!(polys, remove_idx, poly_a, poly_b) + _remove_collinear_points!(alg, polys, remove_idx, poly_a, poly_b) return polys end @@ -110,7 +115,7 @@ _diff_step(x, y) = (x ⊻ y) ? 1 : (-1) #= Polygon with multipolygon difference - note that all intersection regions between `poly_a` and any of the sub-polygons of `multipoly_b` are removed from `poly_a`. =# function _difference( - target::TraitTarget{GI.PolygonTrait}, ::Type{T}, + alg::FosterHormannClipping, target::TraitTarget{GI.PolygonTrait}, ::Type{T}, ::GI.PolygonTrait, poly_a, ::GI.MultiPolygonTrait, multipoly_b; kwargs..., @@ -118,7 +123,7 @@ function _difference( polys = [tuples(poly_a, T)] for poly_b in GI.getpolygon(multipoly_b) isempty(polys) && break - polys = mapreduce(p -> difference(p, poly_b; target), append!, polys) + polys = mapreduce(p -> difference(alg, p, poly_b; target), append!, polys) end return polys end @@ -128,7 +133,7 @@ sub-polygons of `multipoly_a` and `poly_b` will be removed from the correspondin sub-polygon. Unless specified with `fix_multipoly = nothing`, `multipolygon_a` will be validated using the given (default is `UnionIntersectingPolygons()`) correction. =# function _difference( - target::TraitTarget{GI.PolygonTrait}, ::Type{T}, + alg::FosterHormannClipping, target::TraitTarget{GI.PolygonTrait}, ::Type{T}, ::GI.MultiPolygonTrait, multipoly_a, ::GI.PolygonTrait, poly_b; fix_multipoly = UnionIntersectingPolygons(), kwargs..., @@ -139,7 +144,7 @@ function _difference( polys = Vector{_get_poly_type(T)}() sizehint!(polys, GI.npolygon(multipoly_a)) for poly_a in GI.getpolygon(multipoly_a) - append!(polys, difference(poly_a, poly_b; target)) + append!(polys, difference(alg, poly_a, poly_b; target)) end return polys end @@ -150,7 +155,7 @@ corresponding sub-polygon of `multipoly_a`. Unless specified with `fix_multipoly `multipolygon_a` will be validated using the given (default is `UnionIntersectingPolygons()`) correction. =# function _difference( - target::TraitTarget{GI.PolygonTrait}, ::Type{T}, + alg::FosterHormannClipping, target::TraitTarget{GI.PolygonTrait}, ::Type{T}, ::GI.MultiPolygonTrait, multipoly_a, ::GI.MultiPolygonTrait, multipoly_b; fix_multipoly = UnionIntersectingPolygons(), kwargs..., @@ -165,9 +170,9 @@ function _difference( pieces of `multipolygon_a`` are removed, continue to take difference with new shape `polys` =# polys = if i == 1 - difference(multipoly_a, poly_b; target, fix_multipoly) + difference(alg, multipoly_a, poly_b; target, fix_multipoly) else - difference(GI.MultiPolygon(polys), poly_b; target, fix_multipoly) + difference(alg, GI.MultiPolygon(polys), poly_b; target, fix_multipoly) end #= One multipoly_a has been completely covered (and thus removed) there is no need to continue taking the difference =# @@ -176,15 +181,30 @@ function _difference( return polys end +function _difference( + alg::FosterHormannClipping, ::TraitTarget{GI.MultiPolygonTrait}, ::Type{T}, + trait_a::Union{GI.PolygonTrait, GI.MultiPolygonTrait}, polylike_a, + trait_b::Union{GI.PolygonTrait, GI.MultiPolygonTrait}, polylike_b; + fix_multipoly = UnionIntersectingPolygons(), kwargs... +) where T + polys = _difference(alg, TraitTarget(GI.PolygonTrait()), T, trait_a, polylike_a, trait_b, polylike_b; kwargs...) + if isnothing(fix_multipoly) + return GI.MultiPolygon(polys) + else + return fix_multipoly(GI.MultiPolygon(polys)) + end +end + # Many type and target combos aren't implemented function _difference( - ::TraitTarget{Target}, ::Type{T}, + alg::GeometryOpsCore.Algorithm, target::TraitTarget{Target}, ::Type{T}, trait_a::GI.AbstractTrait, geom_a, trait_b::GI.AbstractTrait, geom_b, + kw... ) where {Target, T} @assert( false, - "Difference between $trait_a and $trait_b with target $Target isn't implemented yet.", + "Difference between $trait_a and $trait_b with target $Target and algorithm $alg isn't implemented yet.", ) return nothing end diff --git a/src/methods/clipping/union.jl b/src/methods/clipping/union.jl index 089ec6524..d5d111032 100644 --- a/src/methods/clipping/union.jl +++ b/src/methods/clipping/union.jl @@ -33,19 +33,32 @@ GI.coordinates.(union_poly) ``` """ function union( - geom_a, geom_b, ::Type{T}=Float64; target=nothing, kwargs... + alg::FosterHormannClipping, geom_a, geom_b, ::Type{T}=Float64; target=nothing, kwargs... ) where {T<:AbstractFloat} return _union( - TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b; + alg, TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b; exact = True(), kwargs..., ) end +# fallback definitions +# if no manifold - assume planar (until we have best_manifold) +function union( + geom_a, geom_b, ::Type{T}=Float64; target=nothing, kwargs... +) where {T<:AbstractFloat} + return union(FosterHormannClipping(Planar()), geom_a, geom_b; target, kwargs...) +end + +# if manifold but no algorithm - assume FosterHormannClipping with provided manifold. +function union(m::Manifold, geom_a, geom_b, ::Type{T}=Float64; target=nothing, kwargs...) where {T<:AbstractFloat} + return union(FosterHormannClipping(m), geom_a, geom_b; target, kwargs...) +end + #= This 'union' implementation returns the union of two polygons. The algorithm to determine the union was adapted from "Efficient clipping of efficient polygons," by Greiner and Hormann (1998). DOI: https://doi.org/10.1145/274363.274364 =# function _union( - ::TraitTarget{GI.PolygonTrait}, ::Type{T}, + alg::FosterHormannClipping, ::TraitTarget{GI.PolygonTrait}, ::Type{T}, ::GI.PolygonTrait, poly_a, ::GI.PolygonTrait, poly_b; exact, kwargs..., @@ -54,13 +67,13 @@ function _union( ext_a = GI.getexterior(poly_a) ext_b = GI.getexterior(poly_b) # Then, I get the union of the exteriors - a_list, b_list, a_idx_list = _build_ab_list(T, ext_a, ext_b, _union_delay_cross_f, _union_delay_bounce_f; exact) - polys = _trace_polynodes(T, a_list, b_list, a_idx_list, _union_step, poly_a, poly_b) + a_list, b_list, a_idx_list = _build_ab_list(alg, T, ext_a, ext_b, _union_delay_cross_f, _union_delay_bounce_f; exact) + polys = _trace_polynodes(alg, T, a_list, b_list, a_idx_list, _union_step, poly_a, poly_b) n_pieces = length(polys) # Check if one polygon totally within other and if so, return the larger polygon a_in_b, b_in_a = false, false if n_pieces == 0 # no crossing points, determine if either poly is inside the other - a_in_b, b_in_a = _find_non_cross_orientation(a_list, b_list, ext_a, ext_b; exact) + a_in_b, b_in_a = _find_non_cross_orientation(alg, a_list, b_list, ext_a, ext_b; exact) if a_in_b push!(polys, GI.Polygon([_linearring(tuples(ext_b))])) elseif b_in_a @@ -81,10 +94,10 @@ function _union( end # Add in holes if GI.nhole(poly_a) != 0 || GI.nhole(poly_b) != 0 - _add_union_holes!(polys, a_in_b, b_in_a, poly_a, poly_b; exact) + _add_union_holes!(alg, polys, a_in_b, b_in_a, poly_a, poly_b; exact) end # Remove unneeded collinear points on same edge - _remove_collinear_points!(polys, [false], poly_a, poly_b) + _remove_collinear_points!(alg, polys, [false], poly_a, poly_b) return polys end @@ -108,11 +121,11 @@ _union_step(x, _) = x ? (-1) : 1 #= Add holes from two polygons to the exterior polygon formed by their union. If adding the the holes reveals that the polygons aren't actually intersecting, return the original polygons. =# -function _add_union_holes!(polys, a_in_b, b_in_a, poly_a, poly_b; exact) +function _add_union_holes!(alg::FosterHormannClipping, polys, a_in_b, b_in_a, poly_a, poly_b; exact) if a_in_b - _add_union_holes_contained_polys!(polys, poly_a, poly_b; exact) + _add_union_holes_contained_polys!(alg, polys, poly_a, poly_b; exact) elseif b_in_a - _add_union_holes_contained_polys!(polys, poly_b, poly_a; exact) + _add_union_holes_contained_polys!(alg, polys, poly_b, poly_a; exact) else # Polygons intersect, but neither is contained in the other n_a_holes = GI.nhole(poly_a) ext_poly_a = GI.Polygon(StaticArrays.SVector(GI.getexterior(poly_a))) @@ -125,7 +138,7 @@ function _add_union_holes!(polys, a_in_b, b_in_a, poly_a, poly_b; exact) # Loop over all holes in both original polygons for (i, ih) in enumerate(Iterators.flatten((GI.gethole(poly_a), GI.gethole(poly_b)))) ih = _linearring(ih) - in_ext, _, _ = _line_polygon_interactions(ih, curr_exterior_poly; exact, closed_line = true) + in_ext, _, _ = _line_polygon_interactions(#=TODO: alg.manifold=# ih, curr_exterior_poly; exact, closed_line = true) if !in_ext #= if the hole isn't in the overlapping region between the two polygons, add the hole to the resulting polygon as we know it can't interact with any @@ -138,7 +151,7 @@ function _add_union_holes!(polys, a_in_b, b_in_a, poly_a, poly_b; exact) between poly_a and poly_b within the overlap are added, in addition to all holes in non-overlapping regions =# h_poly = GI.Polygon(StaticArrays.SVector(ih)) - new_holes = difference(h_poly, current_poly; target = GI.PolygonTrait()) + new_holes = difference(alg, h_poly, current_poly; target = GI.PolygonTrait()) append!(polys[1].geom, (GI.getexterior(new_h) for new_h in new_holes)) end if i == n_a_holes @@ -153,12 +166,12 @@ end #= Add holes holes to the union of two polygons where one of the original polygons was inside of the other. If adding the the holes reveal that the polygons aren't actually intersecting, return the original polygons.=# -function _add_union_holes_contained_polys!(polys, interior_poly, exterior_poly; exact) +function _add_union_holes_contained_polys!(alg::FosterHormannClipping, polys, interior_poly, exterior_poly; exact) union_poly = polys[1] ext_int_ring = GI.getexterior(interior_poly) for (i, ih) in enumerate(GI.gethole(exterior_poly)) poly_ih = GI.Polygon(StaticArrays.SVector(ih)) - in_ih, on_ih, out_ih = _line_polygon_interactions(ext_int_ring, poly_ih; exact, closed_line = true) + in_ih, on_ih, out_ih = _line_polygon_interactions(#=TODO: alg.manifold=# ext_int_ring, poly_ih; exact, closed_line = true) if in_ih # at least part of interior polygon exterior is within the ith hole if !on_ih && !out_ih #= interior polygon is completely within the ith hole - polygons aren't @@ -169,7 +182,7 @@ function _add_union_holes_contained_polys!(polys, interior_poly, exterior_poly; else #= interior polygon is partially within the ith hole - area of interior polygon reduces the size of the hole =# - new_holes = difference(poly_ih, interior_poly; target = GI.PolygonTrait()) + new_holes = difference(alg, poly_ih, interior_poly; target = GI.PolygonTrait()) append!(union_poly.geom, (GI.getexterior(new_h) for new_h in new_holes)) end else # none of interior polygon exterior is within the ith hole @@ -183,14 +196,14 @@ function _add_union_holes_contained_polys!(polys, interior_poly, exterior_poly; #= interior polygon's exterior is outside of the ith hole - the interior polygon could either be disjoint from the hole, or contain the hole =# ext_int_poly = GI.Polygon(StaticArrays.SVector(ext_int_ring)) - in_int, _, _ = _line_polygon_interactions(ih, ext_int_poly; exact, closed_line = true) + in_int, _, _ = _line_polygon_interactions(#=TODO: alg.manifold=# ih, ext_int_poly; exact, closed_line = true) if in_int #= interior polygon contains the hole - overlapping holes between the interior and exterior polygons will be added =# for jh in GI.gethole(interior_poly) poly_jh = GI.Polygon(StaticArrays.SVector(jh)) if intersects(poly_ih, poly_jh) - new_holes = intersection(poly_ih, poly_jh; target = GI.PolygonTrait()) + new_holes = intersection(#=TODO: alg, =#poly_ih, poly_jh; target = GI.PolygonTrait()) append!(union_poly.geom, (GI.getexterior(new_h) for new_h in new_holes)) end end @@ -210,7 +223,7 @@ included, unioning these sub-polygons with `poly_a` where they intersect. Unless with `fix_multipoly = nothing`, `multipolygon_b` will be validated using the given (default is `UnionIntersectingPolygons()`) correction. =# function _union( - target::TraitTarget{GI.PolygonTrait}, ::Type{T}, + alg::FosterHormannClipping, target::TraitTarget{GI.PolygonTrait}, ::Type{T}, ::GI.PolygonTrait, poly_a, ::GI.MultiPolygonTrait, multipoly_b; fix_multipoly = UnionIntersectingPolygons(), kwargs..., @@ -220,9 +233,9 @@ function _union( end polys = [tuples(poly_a, T)] for poly_b in GI.getpolygon(multipoly_b) - if intersects(polys[1], poly_b) + if intersects(#=TODO: alg.manifold, =#polys[1], poly_b) # If polygons intersect and form a new polygon, swap out polygon - new_polys = union(polys[1], poly_b; target) + new_polys = union(alg, polys[1], poly_b; target) if length(new_polys) > 1 # case where they intersect by just one point push!(polys, tuples(poly_b, T)) # add poly_b to list else @@ -239,18 +252,18 @@ end #= Multipolygon with polygon union is equivalent to taking the union of the polygon with the multipolygon and thus simply switches the order of operations and calls the above method. =# _union( - target::TraitTarget{GI.PolygonTrait}, ::Type{T}, + alg::FosterHormannClipping, target::TraitTarget{GI.PolygonTrait}, ::Type{T}, ::GI.MultiPolygonTrait, multipoly_a, ::GI.PolygonTrait, poly_b; kwargs..., -) where T = union(poly_b, multipoly_a; target, kwargs...) +) where T = union(alg, poly_b, multipoly_a; target, kwargs...) #= Multipolygon with multipolygon union - note that all of the sub-polygons of `multipoly_a` and the sub-polygons of `multipoly_b` are included and combined together where there are intersections. Unless specified with `fix_multipoly = nothing`, `multipolygon_b` will be validated using the given (default is `UnionIntersectingPolygons()`) correction. =# function _union( - target::TraitTarget{GI.PolygonTrait}, ::Type{T}, + alg::FosterHormannClipping, target::TraitTarget{GI.PolygonTrait}, ::Type{T}, ::GI.MultiPolygonTrait, multipoly_a, ::GI.MultiPolygonTrait, multipoly_b; fix_multipoly = UnionIntersectingPolygons(), kwargs..., @@ -262,15 +275,59 @@ function _union( multipolys = multipoly_b local polys for poly_a in GI.getpolygon(multipoly_a) - polys = union(poly_a, multipolys; target, fix_multipoly) + polys = union(alg, poly_a, multipolys; target, fix_multipoly) multipolys = GI.MultiPolygon(polys) end return polys end +function _union( + alg::FosterHormannClipping, target::TraitTarget{GI.MultiPolygonTrait}, ::Type{T}, + ::GI.PolygonTrait, poly_a, + ::GI.PolygonTrait, poly_b; + kwargs..., +) where T + return UnionIntersectingPolygons(#=TODO: alg.manifold=#)(GI.MultiPolygon([poly_a, poly_b])) +end + +function _union( + alg::FosterHormannClipping, target::TraitTarget{GI.MultiPolygonTrait}, ::Type{T}, + ::GI.PolygonTrait, poly_a, + ::GI.MultiPolygonTrait, multipoly_b; + fix_multipoly = UnionIntersectingPolygons(), + kwargs..., +) where T + res = union(alg, multipoly_b, poly_a; target = GI.PolygonTrait(), kwargs...) + if !isnothing(fix_multipoly) + return fix_multipoly(GI.MultiPolygon(res)) + else + return GI.MultiPolygon(res) + end +end + +# this is the opposite of the above +function _union( + alg::FosterHormannClipping, target::TraitTarget{GI.MultiPolygonTrait}, ::Type{T}, + ::GI.MultiPolygonTrait, multipoly_a, + ::GI.PolygonTrait, poly_b; + fix_multipoly = UnionIntersectingPolygons(), + kwargs..., +) where T + union(alg, poly_b, multipoly_a; target, fix_multipoly, kwargs...) +end + +function _union( + alg::FosterHormannClipping, target::TraitTarget{GI.MultiPolygonTrait}, ::Type{T}, + trait_a::GI.MultiPolygonTrait, multipoly_a, + trait_b::GI.MultiPolygonTrait, multipoly_b; + fix_multipoly = UnionIntersectingPolygons(), kwargs..., +) where T + return GI.MultiPolygon(_union(alg, TraitTarget{GI.PolygonTrait}(), T, trait_a, multipoly_a, trait_b, multipoly_b; fix_multipoly, kwargs...)) +end + # Many type and target combos aren't implemented function _union( - ::TraitTarget{Target}, ::Type{T}, + alg::GeometryOpsCore.Algorithm, target::TraitTarget{Target}, ::Type{T}, trait_a::GI.AbstractTrait, geom_a, trait_b::GI.AbstractTrait, geom_b; kwargs... From 7591f0b941e48b6845537e1aa0d8f5b0286e21ca Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 3 Mar 2025 20:48:12 -0500 Subject: [PATCH 26/32] use TracingError instead of printing out errors --- src/methods/clipping/clipping_processor.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl index 3dcd547cd..4d1641f06 100644 --- a/src/methods/clipping/clipping_processor.jl +++ b/src/methods/clipping/clipping_processor.jl @@ -661,7 +661,9 @@ function _trace_polynodes(alg::FosterHormannClipping{M, A}, ::Type{T}, a_list, b # changed curr_not_intr to curr_not_same_ent_flag same_status, prev_status = true, curr.ent_exit while same_status - @assert visited_pts < total_pts "Clipping tracing hit every point - clipping error. Please open an issue with polygons: $(GI.coordinates(poly_a)) and $(GI.coordinates(poly_b))." + if visited_pts >= total_pts + throw(TracingError("Clipping tracing hit every point - clipping error.", poly_a, poly_b, a_list, b_list, a_idx_list)) + end # Traverse polygon either forwards or backwards idx += step idx = (idx > curr_npoints) ? mod(idx, curr_npoints) : idx From ddee859be6cf3eb8fcafbd884a5957429166d901 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 3 Mar 2025 20:49:58 -0500 Subject: [PATCH 27/32] Implement algorithm in intersection --- src/methods/clipping/intersection.jl | 86 ++++++++++++++++++++-------- 1 file changed, 62 insertions(+), 24 deletions(-) diff --git a/src/methods/clipping/intersection.jl b/src/methods/clipping/intersection.jl index c2c2abe43..fe0d21cb0 100644 --- a/src/methods/clipping/intersection.jl +++ b/src/methods/clipping/intersection.jl @@ -40,28 +40,41 @@ GI.coordinates.(inter_points) ``` """ function intersection( - geom_a, geom_b, ::Type{T}=Float64; target=nothing, kwargs..., + alg::FosterHormannClipping, geom_a, geom_b, ::Type{T}=Float64; target=nothing, kwargs... ) where {T<:AbstractFloat} return _intersection( - TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b; + alg, TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b; exact = True(), kwargs..., ) end +# fallback definitions +# if no manifold - assume planar (until we have best_manifold) +function intersection( + geom_a, geom_b, ::Type{T}=Float64; target=nothing, kwargs... +) where {T<:AbstractFloat} + return intersection(FosterHormannClipping(Planar()), geom_a, geom_b; target, kwargs...) +end + +# if manifold but no algorithm - assume FosterHormannClipping with provided manifold. +function intersection(m::Manifold, geom_a, geom_b, ::Type{T}=Float64; target=nothing, kwargs...) where {T<:AbstractFloat} + return intersection(FosterHormannClipping(m), geom_a, geom_b; target, kwargs...) +end + # Curve-Curve Intersections with target Point _intersection( - ::TraitTarget{GI.PointTrait}, ::Type{T}, + alg::FosterHormannClipping, ::TraitTarget{GI.PointTrait}, ::Type{T}, trait_a::Union{GI.LineTrait, GI.LineStringTrait, GI.LinearRingTrait}, geom_a, trait_b::Union{GI.LineTrait, GI.LineStringTrait, GI.LinearRingTrait}, geom_b; kwargs..., -) where T = _intersection_points(T, trait_a, geom_a, trait_b, geom_b) +) where T = _intersection_points(alg.manifold, alg.accelerator, T, trait_a, geom_a, trait_b, geom_b) #= Polygon-Polygon Intersections with target Polygon The algorithm to determine the intersection was adapted from "Efficient clipping of efficient polygons," by Greiner and Hormann (1998). DOI: https://doi.org/10.1145/274363.274364 =# function _intersection( - ::TraitTarget{GI.PolygonTrait}, ::Type{T}, + alg::FosterHormannClipping, ::TraitTarget{GI.PolygonTrait}, ::Type{T}, ::GI.PolygonTrait, poly_a, ::GI.PolygonTrait, poly_b; exact, kwargs..., @@ -70,10 +83,10 @@ function _intersection( ext_a = GI.getexterior(poly_a) ext_b = GI.getexterior(poly_b) # Then we find the intersection of the exteriors - a_list, b_list, a_idx_list = _build_ab_list(T, ext_a, ext_b, _inter_delay_cross_f, _inter_delay_bounce_f; exact) - polys = _trace_polynodes(T, a_list, b_list, a_idx_list, _inter_step, poly_a, poly_b) + a_list, b_list, a_idx_list = _build_ab_list(alg, T, ext_a, ext_b, _inter_delay_cross_f, _inter_delay_bounce_f; exact) + polys = _trace_polynodes(alg, T, a_list, b_list, a_idx_list, _inter_step, poly_a, poly_b) if isempty(polys) # no crossing points, determine if either poly is inside the other - a_in_b, b_in_a = _find_non_cross_orientation(a_list, b_list, ext_a, ext_b; exact) + a_in_b, b_in_a = _find_non_cross_orientation(alg, a_list, b_list, ext_a, ext_b; exact) if a_in_b push!(polys, GI.Polygon([tuples(ext_a)])) elseif b_in_a @@ -84,10 +97,10 @@ function _intersection( # If the original polygons had holes, take that into account. if GI.nhole(poly_a) != 0 || GI.nhole(poly_b) != 0 hole_iterator = Iterators.flatten((GI.gethole(poly_a), GI.gethole(poly_b))) - _add_holes_to_polys!(T, polys, hole_iterator, remove_idx; exact) + _add_holes_to_polys!(alg, T, polys, hole_iterator, remove_idx; exact) end # Remove unneeded collinear points on same edge - _remove_collinear_points!(polys, remove_idx, poly_a, poly_b) + _remove_collinear_points!(alg, polys, remove_idx, poly_a, poly_b) return polys end @@ -112,7 +125,7 @@ _inter_step(x, _) = x ? 1 : (-1) Unless specified with `fix_multipoly = nothing`, `multipolygon_b` will be validated using the given (default is `UnionIntersectingPolygons()`) correction. =# function _intersection( - target::TraitTarget{GI.PolygonTrait}, ::Type{T}, + alg::FosterHormannClipping, target::TraitTarget{GI.PolygonTrait}, ::Type{T}, ::GI.PolygonTrait, poly_a, ::GI.MultiPolygonTrait, multipoly_b; fix_multipoly = UnionIntersectingPolygons(), kwargs..., @@ -122,7 +135,7 @@ function _intersection( end polys = Vector{_get_poly_type(T)}() for poly_b in GI.getpolygon(multipoly_b) - append!(polys, intersection(poly_a, poly_b; target)) + append!(polys, intersection(alg, poly_a, poly_b; target)) end return polys end @@ -131,11 +144,11 @@ end polygon with the multipolygon and thus simply switches the order of operations and calls the above method. =# _intersection( - target::TraitTarget{GI.PolygonTrait}, ::Type{T}, + alg::FosterHormannClipping, target::TraitTarget{GI.PolygonTrait}, ::Type{T}, ::GI.MultiPolygonTrait, multipoly_a, ::GI.PolygonTrait, poly_b; kwargs..., -) where T = intersection(poly_b, multipoly_a; target , kwargs...) +) where T = intersection(alg, poly_b, multipoly_a; target , kwargs...) #= Multipolygon with multipolygon intersection - note that all intersection regions between any sub-polygons of `multipoly_a` and any of the sub-polygons of `multipoly_b` are counted @@ -143,7 +156,7 @@ as intersection polygons. Unless specified with `fix_multipoly = nothing`, both `multipolygon_a` and `multipolygon_b` will be validated using the given (default is `UnionIntersectingPolygons()`) correction. =# function _intersection( - target::TraitTarget{GI.PolygonTrait}, ::Type{T}, + alg::FosterHormannClipping, target::TraitTarget{GI.PolygonTrait}, ::Type{T}, ::GI.MultiPolygonTrait, multipoly_a, ::GI.MultiPolygonTrait, multipoly_b; fix_multipoly = UnionIntersectingPolygons(), kwargs..., @@ -155,21 +168,37 @@ function _intersection( end polys = Vector{_get_poly_type(T)}() for poly_a in GI.getpolygon(multipoly_a) - append!(polys, intersection(poly_a, multipoly_b; target, fix_multipoly)) + append!(polys, intersection(alg, poly_a, multipoly_b; target, fix_multipoly)) end return polys end +# catch-all method for multipolygontraits +function intersection( + alg::FosterHormannClipping, ::TraitTarget{GI.MultiPolygonTrait}, ::Type{T}, + trait_a::Union{GI.PolygonTrait, GI.MultiPolygonTrait}, polylike_a, + trait_b::Union{GI.PolygonTrait, GI.MultiPolygonTrait}, polylike_b; + fix_multipoly = UnionIntersectingPolygons(), kwargs... +) where T + polys = _intersection(alg, TraitTarget(GI.PolygonTrait()), T, trait_a, polylike_a, trait_b, polylike_b; kwargs...) + if isnothing(fix_multipoly) + return GI.MultiPolygon(polys) + else + return fix_multipoly(GI.MultiPolygon(polys)) + end +end + + # Many type and target combos aren't implemented function _intersection( - ::TraitTarget{Target}, ::Type{T}, + alg::GeometryOpsCore.Algorithm, target::TraitTarget{Target}, ::Type{T}, trait_a::GI.AbstractTrait, geom_a, trait_b::GI.AbstractTrait, geom_b; kwargs..., ) where {Target, T} @assert( false, - "Intersection between $trait_a and $trait_b with target $Target isn't implemented yet.", + "Intersection between $trait_a and $trait_b with target $Target and algorithm $alg isn't implemented yet.", ) return nothing end @@ -193,13 +222,19 @@ inter_points = GO.intersection_points(line1, line2) 1-element Vector{Tuple{Float64, Float64}}: (125.58375366067548, -14.83572303404496) """ -intersection_points(geom_a, geom_b, ::Type{T} = Float64) where T <: AbstractFloat = - _intersection_points(T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b) +intersection_points(geom_a, geom_b, ::Type{T} = Float64) where T <: AbstractFloat = intersection_points(FosterHormannClipping(Planar()), geom_a, geom_b, T) +function intersection_points(alg::FosterHormannClipping{M, A}, geom_a, geom_b, ::Type{T} = Float64) where {M, A, T <: AbstractFloat} + return _intersection_points(alg.manifold, alg.accelerator, T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b) +end + +function intersection_points(m::Manifold, a::IntersectionAccelerator, geom_a, geom_b, ::Type{T} = Float64) where T <: AbstractFloat + return _intersection_points(m, a, T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b) +end #= Calculates the list of intersection points between two geometries, including line segments, line strings, linear rings, polygons, and multipolygons. =# -function _intersection_points(::Type{T}, ::GI.AbstractTrait, a, ::GI.AbstractTrait, b; exact = True()) where T +function _intersection_points(manifold::M, accelerator::A, ::Type{T}, ::GI.AbstractTrait, a, ::GI.AbstractTrait, b; exact = True()) where {M <: Manifold, A <: IntersectionAccelerator, T} # Initialize an empty list of points result = Tuple{T, T}[] # Check if the geometries extents even overlap @@ -243,7 +278,7 @@ intersection point (x,y) while the second is the ratio along the initial lines ( that point. Calculation derivation can be found here: https://stackoverflow.com/questions/563198/ =# -function _intersection_point(::Type{T}, (a1, a2)::Edge, (b1, b2)::Edge; exact) where T +function _intersection_point(manifold::M, ::Type{T}, (a1, a2)::Edge, (b1, b2)::Edge; exact) where {M <: Manifold, T} # Default answer for no intersection line_orient = line_out intr1 = ((zero(T), zero(T)), (zero(T), zero(T))) @@ -266,7 +301,7 @@ function _intersection_point(::Type{T}, (a1, a2)::Edge, (b1, b2)::Edge; exact) w # Determine intersection type and intersection point(s) if a1_orient == a2_orient == b1_orient == b2_orient == 0 # Intersection is collinear if all endpoints lie on the same line - line_orient, intr1, intr2 = _find_collinear_intersection(T, a1, a2, b1, b2, a_ext, b_ext, no_intr_result) + line_orient, intr1, intr2 = _find_collinear_intersection(manifold, T, a1, a2, b1, b2, a_ext, b_ext, no_intr_result) elseif a1_orient == 0 || a2_orient == 0 || b1_orient == 0 || b2_orient == 0 # Intersection is a hinge if the intersection point is an endpoint line_orient = line_hinge @@ -279,13 +314,16 @@ function _intersection_point(::Type{T}, (a1, a2)::Edge, (b1, b2)::Edge; exact) w return line_orient, intr1, intr2 end +# TODO: deprecate this +_intersection_point(::Type{T}, (a1, a2)::Edge, (b1, b2)::Edge; exact) where T = _intersection_point(Planar(), T, (a1, a2), (b1, b2); exact) + #= If lines defined by (a1, a2) and (b1, b2) are collinear, find endpoints of overlapping region if they exist. This could result in three possibilities. First, there could be no overlapping region, in which case, the default 'no_intr_result' intersection information is returned. Second, the two regions could just meet at one shared endpoint, in which case it is a hinge intersection with one intersection point. Otherwise, it is a overlapping intersection defined by two of the endpoints of the line segments. =# -function _find_collinear_intersection(::Type{T}, a1, a2, b1, b2, a_ext, b_ext, no_intr_result) where T +function _find_collinear_intersection(manifold::M, ::Type{T}, a1, a2, b1, b2, a_ext, b_ext, no_intr_result) where {M <: Manifold, T} # Define default return for no intersection points line_orient, intr1, intr2 = no_intr_result # Determine collinear line overlaps From 89ec64c8456640914fbd23d60b39b7dbd55a87ed Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Fri, 14 Mar 2025 11:40:55 -0800 Subject: [PATCH 28/32] Fix typo that prevented multipolygon intersection from taking place --- src/methods/clipping/intersection.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/methods/clipping/intersection.jl b/src/methods/clipping/intersection.jl index fe0d21cb0..52167b1ad 100644 --- a/src/methods/clipping/intersection.jl +++ b/src/methods/clipping/intersection.jl @@ -174,7 +174,7 @@ function _intersection( end # catch-all method for multipolygontraits -function intersection( +function _intersection( alg::FosterHormannClipping, ::TraitTarget{GI.MultiPolygonTrait}, ::Type{T}, trait_a::Union{GI.PolygonTrait, GI.MultiPolygonTrait}, polylike_a, trait_b::Union{GI.PolygonTrait, GI.MultiPolygonTrait}, polylike_b; From fdb554c1c5c2cbd87fb7a97e5a998e7672571b90 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sun, 20 Apr 2025 07:10:02 -0400 Subject: [PATCH 29/32] Apply suggestions from code review Co-authored-by: Rafael Schouten --- src/methods/clipping/clipping_processor.jl | 4 +--- src/methods/clipping/cut.jl | 1 - src/methods/clipping/difference.jl | 3 --- src/methods/clipping/intersection.jl | 4 ---- src/methods/clipping/union.jl | 5 ----- src/methods/geom_relations/geom_geom_processors.jl | 1 - 6 files changed, 1 insertion(+), 17 deletions(-) diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl index 4d1641f06..17768ac19 100644 --- a/src/methods/clipping/clipping_processor.jl +++ b/src/methods/clipping/clipping_processor.jl @@ -33,7 +33,7 @@ struct AutoAccelerator <: IntersectionAccelerator end """ FosterHormannClipping{M <: Manifold, A <: Union{Nothing, Accelerator}} <: GeometryOpsCore.Algorithm{M} -A type that represents the Foster-Hormann clipping algorithm. +Applies the Foster-Hormann clipping algorithm. # Arguments - `manifold::M`: The manifold on which the algorithm operates. @@ -45,8 +45,6 @@ struct FosterHormannClipping{M <: Manifold, A <: IntersectionAccelerator} <: Geo # TODO: add exact flag # TODO: should exact flag be in the type domain? end - - FosterHormannClipping(; manifold::Manifold = Planar(), accelerator = nothing) = FosterHormannClipping(manifold, isnothing(accelerator) ? AutoAccelerator() : accelerator) FosterHormannClipping(manifold::Manifold, accelerator::Union{Nothing, IntersectionAccelerator} = nothing) = FosterHormannClipping(manifold, isnothing(accelerator) ? AutoAccelerator() : accelerator) FosterHormannClipping(accelerator::Union{Nothing, IntersectionAccelerator}) = FosterHormannClipping(Planar(), isnothing(accelerator) ? AutoAccelerator() : accelerator) diff --git a/src/methods/clipping/cut.jl b/src/methods/clipping/cut.jl index c3427731a..b956035e5 100644 --- a/src/methods/clipping/cut.jl +++ b/src/methods/clipping/cut.jl @@ -60,7 +60,6 @@ GI.coordinates.(cut_polys) """ cut(geom, line, ::Type{T} = Float64) where {T <: AbstractFloat} = cut(FosterHormannClipping(), geom, line, T) cut(m::Manifold, geom, line, ::Type{T} = Float64) where {T <: AbstractFloat} = cut(FosterHormannClipping(m), geom, line, T) - cut(alg::FosterHormannClipping{M, A}, geom, line, ::Type{T} = Float64) where {T <: AbstractFloat, M, A} = _cut(alg, T, GI.trait(geom), geom, GI.trait(line), line; exact = True()) diff --git a/src/methods/clipping/difference.jl b/src/methods/clipping/difference.jl index 19c72491b..ff725a881 100644 --- a/src/methods/clipping/difference.jl +++ b/src/methods/clipping/difference.jl @@ -40,7 +40,6 @@ function difference( exact = True(), kwargs..., ) end - # fallback definitions difference(geom_a, geom_b, ::Type{T} = Float64; target=nothing, kwargs...) where T = difference(FosterHormannClipping(Planar()), geom_a, geom_b, T; target, kwargs...) # if manifold but no algorithm - assume FosterHormannClipping with provided manifold. @@ -180,7 +179,6 @@ function _difference( end return polys end - function _difference( alg::FosterHormannClipping, ::TraitTarget{GI.MultiPolygonTrait}, ::Type{T}, trait_a::Union{GI.PolygonTrait, GI.MultiPolygonTrait}, polylike_a, @@ -194,7 +192,6 @@ function _difference( return fix_multipoly(GI.MultiPolygon(polys)) end end - # Many type and target combos aren't implemented function _difference( alg::GeometryOpsCore.Algorithm, target::TraitTarget{Target}, ::Type{T}, diff --git a/src/methods/clipping/intersection.jl b/src/methods/clipping/intersection.jl index 52167b1ad..47d82077c 100644 --- a/src/methods/clipping/intersection.jl +++ b/src/methods/clipping/intersection.jl @@ -47,7 +47,6 @@ function intersection( exact = True(), kwargs..., ) end - # fallback definitions # if no manifold - assume planar (until we have best_manifold) function intersection( @@ -55,7 +54,6 @@ function intersection( ) where {T<:AbstractFloat} return intersection(FosterHormannClipping(Planar()), geom_a, geom_b; target, kwargs...) end - # if manifold but no algorithm - assume FosterHormannClipping with provided manifold. function intersection(m::Manifold, geom_a, geom_b, ::Type{T}=Float64; target=nothing, kwargs...) where {T<:AbstractFloat} return intersection(FosterHormannClipping(m), geom_a, geom_b; target, kwargs...) @@ -103,7 +101,6 @@ function _intersection( _remove_collinear_points!(alg, polys, remove_idx, poly_a, poly_b) return polys end - # # Helper functions for Intersections with Greiner and Hormann Polygon Clipping #= When marking the crossing status of a delayed crossing, the chain start point is bouncing @@ -172,7 +169,6 @@ function _intersection( end return polys end - # catch-all method for multipolygontraits function _intersection( alg::FosterHormannClipping, ::TraitTarget{GI.MultiPolygonTrait}, ::Type{T}, diff --git a/src/methods/clipping/union.jl b/src/methods/clipping/union.jl index d5d111032..17af626da 100644 --- a/src/methods/clipping/union.jl +++ b/src/methods/clipping/union.jl @@ -280,7 +280,6 @@ function _union( end return polys end - function _union( alg::FosterHormannClipping, target::TraitTarget{GI.MultiPolygonTrait}, ::Type{T}, ::GI.PolygonTrait, poly_a, @@ -289,7 +288,6 @@ function _union( ) where T return UnionIntersectingPolygons(#=TODO: alg.manifold=#)(GI.MultiPolygon([poly_a, poly_b])) end - function _union( alg::FosterHormannClipping, target::TraitTarget{GI.MultiPolygonTrait}, ::Type{T}, ::GI.PolygonTrait, poly_a, @@ -304,7 +302,6 @@ function _union( return GI.MultiPolygon(res) end end - # this is the opposite of the above function _union( alg::FosterHormannClipping, target::TraitTarget{GI.MultiPolygonTrait}, ::Type{T}, @@ -315,7 +312,6 @@ function _union( ) where T union(alg, poly_b, multipoly_a; target, fix_multipoly, kwargs...) end - function _union( alg::FosterHormannClipping, target::TraitTarget{GI.MultiPolygonTrait}, ::Type{T}, trait_a::GI.MultiPolygonTrait, multipoly_a, @@ -324,7 +320,6 @@ function _union( ) where T return GI.MultiPolygon(_union(alg, TraitTarget{GI.PolygonTrait}(), T, trait_a, multipoly_a, trait_b, multipoly_b; fix_multipoly, kwargs...)) end - # Many type and target combos aren't implemented function _union( alg::GeometryOpsCore.Algorithm, target::TraitTarget{Target}, ::Type{T}, diff --git a/src/methods/geom_relations/geom_geom_processors.jl b/src/methods/geom_relations/geom_geom_processors.jl index 8563fb174..5f336f652 100644 --- a/src/methods/geom_relations/geom_geom_processors.jl +++ b/src/methods/geom_relations/geom_geom_processors.jl @@ -525,7 +525,6 @@ function _point_filled_curve_orientation( end return iseven(k) ? out : in end - _point_filled_curve_orientation( point, curve; in::T = point_in, on::T = point_on, out::T = point_out, exact, From e29407b5c28d073ee166abae09915873e424aff3 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sun, 20 Apr 2025 07:47:10 -0400 Subject: [PATCH 30/32] add a test for TracingError shows --- test/methods/clipping/polygon_clipping.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/test/methods/clipping/polygon_clipping.jl b/test/methods/clipping/polygon_clipping.jl index 432650e94..eda0d14b2 100644 --- a/test/methods/clipping/polygon_clipping.jl +++ b/test/methods/clipping/polygon_clipping.jl @@ -215,3 +215,15 @@ end @testset "Intersection" begin test_clipping(GO.intersection, LG.intersection, "intersection") end @testset "Union" begin test_clipping(GO.union, LG.union, "union") end @testset "Difference" begin test_clipping(GO.difference, LG.difference, "difference") end + +@testset "TracingError show" begin + message = "Test tracing error" + poly_a_small = p1 + poly_b_small = p2 + + poly_a_large, poly_b_large = GO.segmentize.((poly_a_small, poly_b_small), max_distance = 0.001) + + @test_throws "Test tracing error" throw(GO.TracingError(message, poly_a_small, poly_b_small, GO.PolyNode{Float64}[], GO.PolyNode{Float64}[], Int[])) + @test_throws "$(GI.coordinates(poly_a_small))" throw(GO.TracingError(message, poly_a_small, poly_b_small, GO.PolyNode{Float64}[], GO.PolyNode{Float64}[], Int[])) + @test_throws "The polygons are contained in the exception object" throw(GO.TracingError(message, poly_a_large, poly_b_large, GO.PolyNode{Float64}[], GO.PolyNode{Float64}[], Int[])) +end \ No newline at end of file From 353b7755fb674cdabc6abbf3567e78de4d85245e Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sun, 20 Apr 2025 07:51:18 -0400 Subject: [PATCH 31/32] add explanatory comments --- test/methods/clipping/polygon_clipping.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/test/methods/clipping/polygon_clipping.jl b/test/methods/clipping/polygon_clipping.jl index eda0d14b2..2fc4575fa 100644 --- a/test/methods/clipping/polygon_clipping.jl +++ b/test/methods/clipping/polygon_clipping.jl @@ -223,7 +223,12 @@ end poly_a_large, poly_b_large = GO.segmentize.((poly_a_small, poly_b_small), max_distance = 0.001) + # Test that the message is contained in the exception @test_throws "Test tracing error" throw(GO.TracingError(message, poly_a_small, poly_b_small, GO.PolyNode{Float64}[], GO.PolyNode{Float64}[], Int[])) + # Test that the coordinates of the polygons are contained in the exception, + # if the polygon is small enough @test_throws "$(GI.coordinates(poly_a_small))" throw(GO.TracingError(message, poly_a_small, poly_b_small, GO.PolyNode{Float64}[], GO.PolyNode{Float64}[], Int[])) + # Test that the coordinates of the polygons are not contained in the exception, + # if the polygon is large enoughs @test_throws "The polygons are contained in the exception object" throw(GO.TracingError(message, poly_a_large, poly_b_large, GO.PolyNode{Float64}[], GO.PolyNode{Float64}[], Int[])) end \ No newline at end of file From eef41bd69f4c2078ff784ee3707596aaa8d45aca Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 21 Apr 2025 15:39:07 -0400 Subject: [PATCH 32/32] Update clipping_processor.jl --- src/methods/clipping/clipping_processor.jl | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl index 17768ac19..94bf4e859 100644 --- a/src/methods/clipping/clipping_processor.jl +++ b/src/methods/clipping/clipping_processor.jl @@ -12,14 +12,15 @@ The abstract supertype for all intersection accelerator types. The idea is that these speed up the edge-edge intersection checking process, perhaps at the cost of memory. -The naive case is [`NestedLoop`](@ref), which is just a nested loop, running in O(n*m) time. +The naive case is `NestedLoop`, which is just a nested loop, running in O(n*m) time. -Then we have [`SingleSTRtree`](@ref), which is a single STRtree, running in O(n*log(m)) time. +Then we have `SingleSTRtree`, which is a single STRtree, running in O(n*log(m)) time. -Then we have [`DoubleSTRtree`](@ref), which is am simultaneous double-tree traversal of two STRtrees. +Then we have `DoubleSTRtree`, which is a simultaneous double-tree traversal of two STRtrees. -Finally, we have [`AutoAccelerator`](@ref), which is an automatic accelerator that chooses the best -accelerator based on the size of the input polygons. This gets materialized in build_a_list for now. +Finally, we have `AutoAccelerator`, which chooses the best +accelerator based on the size of the input polygons. This gets materialized in `build_a_list` for now. +`AutoAccelerator` should also try to respect existing spatial indices, if they exist. """ abstract type IntersectionAccelerator end struct NestedLoop <: IntersectionAccelerator end