Spatial tree interface (#297)

Co-authored-by: Rafael Schouten <[email protected]>


This PR performs all the necessary actions to add a spatial tree interface and surrounding machinery to GeometryOps.  Actually accelerating GeometryOps functions with the spatial trees is left for another PR - we don't tackle that in this one. 

This implements, at a high level, functional, non-allocating single-tree and dual tree queries over arbitrary predicates.

## SpatialTreeInterface

Introduces a new SpatialTreeInterface module that provides an interface for spatial trees, built somewhat on top of AbstractTrees.jl.


## LoopStateMachine

Introduces a new `LoopStateMachine` module that provides utilities for returning state from functions running inside loops. This is particularly useful for operations like clipping where state transitions or early termination may be needed.

The module includes:
- `Action` struct for representing different control flow actions
- `@controlflow` macro for processing actions within loops
- Support for `:continue`, `:break`, `:return`, and `:full_return` actions

Comprehensive tests verify all supported actions work correctly.

Author: asinghvi17

Project.toml

@@ -4,6 +4,7 @@ authors = ["Anshul Singhvi <[email protected]>", "Rafael Schouten <rafaels
 version = "0.1.18"
 
 [deps]
+AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
 AdaptivePredicates = "35492f91-a3bd-45ad-95db-fcad7dcfedb7"
 CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
@@ -32,6 +33,7 @@ GeometryOpsProjExt = "Proj"
 GeometryOpsTGGeometryExt = "TGGeometry"
 
 [compat]
+AbstractTrees = "0.4"
 AdaptivePredicates = "1.2"
 CoordinateTransformations = "0.5, 0.6"
 DataAPI = "1"
@@ -65,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"
@@ -74,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"]

docs/.gitignore

@@ -2,4 +2,6 @@ node_modules/
 build/
 package-lock.json
 src/source/
-Manifest.toml
+Manifest.toml
+
+src/call_notes.md

docs/Project.toml

@@ -16,6 +16,7 @@ DocumenterVitepress = "4710194d-e776-4893-9690-8d956a29c365"
 DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110"
+Extents = "411431e0-e8b7-467b-b5e0-f676ba4f2910"
 FlexiJoins = "e37f2e79-19fa-4eb7-8510-b63b51fe0a37"
 GADM = "a8dd9ffe-31dc-4cf5-a379-ea69100a8233"
 GeoDataFrames = "62cb38b5-d8d2-4862-a48e-6a340996859f"
@@ -42,6 +43,7 @@ Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+SortTileRecursiveTree = "746ee33f-1797-42c2-866d-db2fce69d14d"
 TGGeometry = "d7e755d2-3c95-4bcf-9b3c-79ab1a78647b"
 
 [sources]

src/GeometryOps.jl

@@ -36,9 +36,15 @@ 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")
+include("utils/SpatialTreeInterface/SpatialTreeInterface.jl")
+
+using .LoopStateMachine, .SpatialTreeInterface
+
+
 include("methods/angles.jl")
 include("methods/area.jl")
 include("methods/barycentric.jl")

src/utils/LoopStateMachine/LoopStateMachine.jl

@@ -0,0 +1,123 @@
+"""
+    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{T}(:unnamed, x)
+Action(x::Symbol) = Action(x, nothing)
+
+function Base.show(io::IO, action::Action{T}) where T
+    print(io, "Action")
+    print(io, "(:$(action.name)")
+    if isnothing(action.x)
+        print(io, ")")
+    else
+        print(io, ", ",action.x, ")")
+    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(...)
+
+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(UnrecognizedActionException($varname.name))
+            end
+        else
+            $varname
+        end
+    end
+end
+# COV_EXCL_STOP
+
+# You can define more actions as you desire.
+
+end

src/utils/SpatialTreeInterface/README.md

@@ -0,0 +1,47 @@
+# SpatialTreeInterface.jl
+
+A simple interface for spatial tree types.
+
+## What is a spatial tree?
+
+- 2 dimensional extents
+- Parent nodes encompass all leaf nodes
+- Leaf nodes contain references to the geometries they represent as indices (or so we assume here)
+
+## Why is this useful?
+
+- It allows us to write algorithms that can work with any spatial tree type, without having to know the details of the tree type.
+    - for example, dual tree traversal / queries
+- It allows us to flexibly and easily swap out and use different tree types, depending on the problem at hand.
+
+This is also a zero cost interface if implemented correctly!  Verified implementations exist for "flat" trees like the "Natural Index" from `tg`, and "hierarchical" trees like the `STRtree` from `SortTileRecursiveTree.jl`.
+
+## Interface
+
+- `isspatialtree(tree)::Bool`
+- `isleaf(node)::Bool` - is the node a leaf node?  In this context, a leaf node is a node that does not have other nodes as its children, but stores a list of indices and extents (even if implicit).
+- `getchild(node)` - get the children of a node.  This may be materialized if necessary or available, but can also be lazy (like a generator).
+- `getchild(node, i)` - get the `i`-th child of a node.
+- `nchild(node)::Int` - the number of children of a node.
+- `child_indices_extents(node)` - an iterator over the indices and extents of the children of a **leaf** node.
+
+These are the only methods that are required to be implemented.  
+
+Optionally, one may define:
+- `node_extent(node)` - get the extent of a node.  This falls back to `GI.extent` but can potentially be overridden if you want to return a different but extent-like object.
+
+They enable the generic query functions described below:
+
+## Query functions
+
+- `do_query(f, predicate, node)` - call `f(i)` for each index `i` in `node` that satisfies `predicate(extent(i))`.
+- `do_dual_query(f, predicate, tree1, tree2)` - call `f(i1, i2)` for each index `i1` in `tree1` and `i2` in `tree2` that satisfies `predicate(extent(i1), extent(i2))`.
+
+These are both completely non-allocating, and will only call `f` for indices that satisfy the predicate.
+You can of course build a standard query interface on top of `do_query` if you want - that's simply:
+```julia
+a = Int[]
+do_query(Base.Fix1(push!, a), predicate, node)
+```
+where `predicate` might be `Base.Fix1(Extents.intersects, extent_to_query)`.
+

src/utils/SpatialTreeInterface/SpatialTreeInterface.jl

@@ -0,0 +1,59 @@
+module SpatialTreeInterface
+
+import ..LoopStateMachine: @controlflow
+
+import Extents
+import GeoInterface as GI
+import AbstractTrees
+
+# public isspatialtree, isleaf, getchild, nchild, child_indices_extents, node_extent
+export query, do_query
+export FlatNoTree
+
+# The spatial tree interface and its implementations are defined here.
+include("interface.jl")
+include("implementations.jl")
+
+# 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")
+
+# 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")
+
+
+"""
+    query(tree, predicate)
+
+Return a sorted list of indices of the tree that satisfy the predicate.
+"""
+function query(tree, predicate)
+    a = Int[]
+    depth_first_search(Base.Fix1(push!, a), sanitize_predicate(predicate), tree)
+    return sort!(a)
+end
+
+
+"""
+    sanitize_predicate(pred)
+
+Convert a predicate to a function that returns a Boolean.
+
+If `pred` is an Extent, convert it to a function that returns a Boolean by intersecting with the extent.
+If `pred` is a geometry, convert it to an extent first, then wrap in Extents.intersects.
+
+Otherwise, return the predicate unchanged.
+
+
+Users and developers may overload this function to provide custom behaviour when something is passed in.
+"""
+sanitize_predicate(pred) = sanitize_predicate(GI.trait(pred), pred)
+sanitize_predicate(::Nothing, pred) = pred
+sanitize_predicate(::GI.AbstractTrait, pred) = sanitize_predicate(GI.extent(pred))
+sanitize_predicate(pred::Extents.Extent) = Base.Fix1(Extents.intersects, pred)
+
+
+end # module SpatialTreeInterface