From 9168b81b66783f33b709b029614ce98f7e0c4279 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Mon, 3 Mar 2025 20:35:35 -0500
Subject: [PATCH 01/16] Add AbstractTrees + constrain STR to performant
 versions

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 64c916de6..ec62f505a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -48,7 +48,7 @@ GeometryOpsCore = "=0.1.5"
 LibGEOS = "0.9.2"
 LinearAlgebra = "1"
 Proj = "1"
-SortTileRecursiveTree = "0.1"
+SortTileRecursiveTree = "0.1.2"
 Statistics = "1"
 TGGeometry = "0.1"
 Tables = "1"

From f03e6502ad71a846fc83d293ebcdddd3cda6624c Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Mon, 3 Mar 2025 20:51:29 -0500
Subject: [PATCH 02/16] include LoopStateMachine and STI

---
 src/GeometryOps.jl | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl
index db7efd2ee..29e1a4f3d 100644
--- a/src/GeometryOps.jl
+++ b/src/GeometryOps.jl
@@ -21,12 +21,15 @@ using GeoInterface
 using GeometryBasics
 using LinearAlgebra, Statistics
 
+using GeometryBasics.StaticArrays
+
 import Tables, DataAPI
-import GeometryBasics.StaticArrays
 import DelaunayTriangulation # for convex hull and triangulation
 import ExactPredicates
 import Base.@kwdef
 import GeoInterface.Extents: Extents
+import SortTileRecursiveTree
+import SortTileRecursiveTree: STRtree
 
 const GI = GeoInterface
 const GB = GeometryBasics

From c4247899aed6204cb7e078e99ec6eea781569f85 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Mon, 3 Mar 2025 20:52:03 -0500
Subject: [PATCH 03/16] Implement natural indexing from `tg`

---
 src/GeometryOps.jl           |   3 +
 src/utils/NaturalIndexing.jl | 235 +++++++++++++++++++++++++++++++++++
 2 files changed, 238 insertions(+)
 create mode 100644 src/utils/NaturalIndexing.jl

diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl
index 29e1a4f3d..d9c307e8f 100644
--- a/src/GeometryOps.jl
+++ b/src/GeometryOps.jl
@@ -47,6 +47,9 @@ include("utils/SpatialTreeInterface/SpatialTreeInterface.jl")
 
 using .LoopStateMachine, .SpatialTreeInterface
 
+include("utils/NaturalIndexing.jl")
+using .NaturalIndexing
+
 
 include("methods/angles.jl")
 include("methods/area.jl")
diff --git a/src/utils/NaturalIndexing.jl b/src/utils/NaturalIndexing.jl
new file mode 100644
index 000000000..f972fe4ac
--- /dev/null
+++ b/src/utils/NaturalIndexing.jl
@@ -0,0 +1,235 @@
+module NaturalIndexing
+
+import GeoInterface as GI
+import Extents
+
+using ..SpatialTreeInterface
+
+"""
+    NaturalLevel{E <: Extents.Extent}
+
+A level in the natural tree.  Stored in a vector in [`NaturalIndex`](@ref).
+
+- `extents` is a vector of extents of the children of the node
+"""
+struct NaturalLevel{E <: Extents.Extent}
+    extents::Vector{E} # child extents
+end
+
+Base.show(io::IO, level::NaturalLevel) = print(io, "NaturalLevel($(length(level.extents)) extents)")
+Base.show(io::IO, ::MIME"text/plain", level::NaturalLevel) = Base.show(io, level)
+
+"""
+    NaturalIndex{E <: Extents.Extent}
+
+A natural tree index.  Stored in a vector in [`NaturalIndex`](@ref).
+
+- `nodecapacity` is the "spread", number of children per node
+- `extent` is the extent of the tree
+- `levels` is a vector of [`NaturalLevel`](@ref)s
+"""
+struct NaturalIndex{E <: Extents.Extent}
+    nodecapacity::Int # "spread", number of children per node
+    extent::E
+    levels::Vector{NaturalLevel{E}}
+end
+
+GI.extent(idx::NaturalIndex) = idx.extent
+Extents.extent(idx::NaturalIndex) = idx.extent
+
+function Base.show(io::IO, ::MIME"text/plain", idx::NaturalIndex)
+    println(io, "NaturalIndex with $(length(idx.levels)) levels and $(idx.nodecapacity) children per node")
+    println(io, "extent: $(idx.extent)")
+end
+
+function Base.show(io::IO, idx::NaturalIndex)
+    println(io, "NaturalIndex($(length(idx.levels)) levels, $(idx.extent))")
+end
+
+function NaturalIndex(geoms; nodecapacity = 32)
+    e1 = GI.extent(first(geoms))
+    E = typeof(e1)
+    return NaturalIndex{E}(geoms; nodecapacity = nodecapacity)
+end
+
+function NaturalIndex{E}(geoms; nodecapacity = 32) where E <: Extents.Extent
+    last_level_extents = GI.extent.(geoms)
+    return NaturalIndex{E}(last_level_extents; nodecapacity = nodecapacity)
+end
+
+function NaturalIndex(last_level_extents::Vector{E}; nodecapacity = 32) where E <: Extents.Extent
+    return NaturalIndex{E}(last_level_extents; nodecapacity = nodecapacity)
+end
+
+function NaturalIndex{E}(last_level_extents::Vector{E}; nodecapacity = 32) where E <: Extents.Extent
+    ngeoms = length(last_level_extents)
+    last_level = NaturalLevel(last_level_extents)
+
+    nlevels = _number_of_levels(nodecapacity, ngeoms)
+
+    levels = Vector{NaturalLevel{E}}(undef, nlevels)
+    levels[end] = last_level
+
+    for level_index in (nlevels-1):(-1):1
+        prev_level = levels[level_index+1] # this is always instantiated
+        nrects = _number_of_keys(nodecapacity, nlevels - (level_index), ngeoms)
+        # @show level_index nrects
+        extents = [
+            begin
+                start = (rect_index - 1) * nodecapacity + 1
+                stop = min(start + nodecapacity - 1, length(prev_level.extents))
+                reduce(Extents.union, view(prev_level.extents, start:stop))
+            end
+            for rect_index in 1:nrects
+        ]
+        levels[level_index] = NaturalLevel(extents)
+    end
+
+    return NaturalIndex(nodecapacity, reduce(Extents.union, levels[1].extents), levels)
+
+end
+
+function _number_of_keys(nodecapacity::Int, level::Int, ngeoms::Int)
+    return ceil(Int, ngeoms / (nodecapacity ^ (level)))
+end
+
+"""
+    _number_of_levels(nodecapacity::Int, ngeoms::Int)
+
+Calculate the number of levels in a natural tree for a given number of geometries and node capacity.
+
+## How this works
+
+The number of keys in a level is given by `ngeoms / nodecapacity ^ level`.
+
+The number of levels is the smallest integer such that the number of keys in the last level is 1.
+So it goes - if that makes sense.
+"""
+function _number_of_levels(nodecapacity::Int, ngeoms::Int)
+    level = 1
+    while _number_of_keys(nodecapacity, level, ngeoms) > 1
+        level += 1
+    end
+    return level
+end
+
+
+# This is like a pointer to a node in the tree.
+"""
+    NaturalTreeNode{E <: Extents.Extent}
+
+A reference to a node in the natural tree.  Kind of like a tree cursor.
+
+- `parent_index` is a pointer to the parent index
+- `level` is the level of the node in the tree
+- `index` is the index of the node in the level
+- `extent` is the extent of the node
+"""
+struct NaturalTreeNode{E <: Extents.Extent}
+    parent_index::NaturalIndex{E}
+    level::Int
+    index::Int
+    extent::E
+end
+
+Extents.extent(node::NaturalTreeNode) = node.extent
+
+# What does SpatialTreeInterface require of trees?
+# - Parents completely cover their children
+# - `GI.extent(node)` returns `Extent` 
+#   - can mean that `Extents.extent(node)` returns the extent of the node
+# - `nchild(node)` returns the number of children of the node
+# - `getchild(node)` returns an iterator over all children of the node
+# - `getchild(node, i)` returns the i-th child of the node
+# - `isleaf(node)` returns a boolean indicating whether the node is a leaf
+# - `child_indices_extents(node)` returns an iterator over the indices and extents of the children of the node
+
+function SpatialTreeInterface.nchild(node::NaturalTreeNode)
+    start_idx = (node.index - 1) * node.parent_index.nodecapacity + 1
+    stop_idx = min(start_idx + node.parent_index.nodecapacity - 1, length(node.parent_index.levels[node.level+1].extents))
+    return stop_idx - start_idx + 1
+end
+
+function SpatialTreeInterface.getchild(node::NaturalTreeNode, i::Int)
+    child_index = (node.index - 1) * node.parent_index.nodecapacity + i
+    return NaturalTreeNode(
+        node.parent_index, 
+        node.level + 1, # increment level by 1
+        child_index, # index of this particular child
+        node.parent_index.levels[node.level+1].extents[child_index] # the extent of this child
+    )
+end
+
+# Get all children of a node
+function SpatialTreeInterface.getchild(node::NaturalTreeNode)
+    return (SpatialTreeInterface.getchild(node, i) for i in 1:SpatialTreeInterface.nchild(node))
+end
+
+SpatialTreeInterface.isleaf(node::NaturalTreeNode) = node.level == length(node.parent_index.levels) - 1
+
+function SpatialTreeInterface.child_indices_extents(node::NaturalTreeNode)
+    start_idx = (node.index - 1) * node.parent_index.nodecapacity + 1
+    stop_idx = min(start_idx + node.parent_index.nodecapacity - 1, length(node.parent_index.levels[node.level+1].extents))
+    return ((i, node.parent_index.levels[node.level+1].extents[i]) for i in start_idx:stop_idx)
+end
+
+# implementation for "root node" / top level tree
+
+SpatialTreeInterface.isleaf(node::NaturalIndex) = length(node.levels) == 1
+
+SpatialTreeInterface.nchild(node::NaturalIndex) = length(node.levels[1].extents)
+
+SpatialTreeInterface.getchild(node::NaturalIndex) = SpatialTreeInterface.getchild(NaturalTreeNode(node, 0, 1, node.extent))
+SpatialTreeInterface.getchild(node::NaturalIndex, i) = SpatialTreeInterface.getchild(NaturalTreeNode(node, 0, 1, node.extent), i)
+
+SpatialTreeInterface.child_indices_extents(node::NaturalIndex) = (i_ext for i_ext in enumerate(node.levels[1].extents))
+
+struct NaturallyIndexedRing
+    points::Vector{Tuple{Float64, Float64}}
+    index::NaturalIndex{Extents.Extent{(:X, :Y), NTuple{2, NTuple{2, Float64}}}}
+end
+
+function NaturallyIndexedRing(points::Vector{Tuple{Float64, Float64}}; nodecapacity = 32)
+    index = NaturalIndex(GO.edge_extents(GI.LinearRing(points)); nodecapacity)
+    return NaturallyIndexedRing(points, index)
+end
+
+NaturallyIndexedRing(ring::NaturallyIndexedRing) = ring
+
+function GI.convert(::Type{NaturallyIndexedRing}, ::GI.LinearRingTrait, geom)
+    points = GO.tuples(geom).geom
+    return NaturallyIndexedRing(points)
+end
+
+Base.show(io::IO, ::MIME"text/plain", ring::NaturallyIndexedRing) = Base.show(io, ring)
+
+Base.show(io::IO, ring::NaturallyIndexedRing) = print(io, "NaturallyIndexedRing($(length(ring.points)) points) with index $(sprint(show, ring.index))")
+
+GI.ncoord(::GI.LinearRingTrait, ring::NaturallyIndexedRing) = 2
+GI.is3d(::GI.LinearRingTrait, ring::NaturallyIndexedRing) = false
+GI.ismeasured(::GI.LinearRingTrait, ring::NaturallyIndexedRing) = false
+
+GI.npoint(::GI.LinearRingTrait, ring::NaturallyIndexedRing) = length(ring.points)
+GI.getpoint(::GI.LinearRingTrait, ring::NaturallyIndexedRing) = ring.points[i]
+GI.getpoint(::GI.LinearRingTrait, ring::NaturallyIndexedRing, i::Int) = ring.points[i]
+
+GI.ngeom(::GI.LinearRingTrait, ring::NaturallyIndexedRing) = length(ring.points)
+GI.getgeom(::GI.LinearRingTrait, ring::NaturallyIndexedRing) = ring.points
+GI.getgeom(::GI.LinearRingTrait, ring::NaturallyIndexedRing, i::Int) = ring.points[i]
+
+Extents.extent(ring::NaturallyIndexedRing) = ring.index.extent
+
+GI.isgeometry(::NaturallyIndexedRing) = true
+GI.geomtrait(::NaturallyIndexedRing) = GI.LinearRingTrait()
+
+function prepare_naturally(geom)
+    return GO.apply(GI.PolygonTrait(), geom) do poly
+        return GI.Polygon([GI.convert(NaturallyIndexedRing, GI.LinearRingTrait(), ring) for ring in GI.getring(poly)])
+    end
+end
+
+export NaturalTree, NaturallyIndexedRing, prepare_naturally
+
+end # module NaturalIndexing
+
+using .NaturalIndexing
\ No newline at end of file

From 182bb5c106862ae0cef2f7cfd602f7d004f5f199 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Mon, 3 Mar 2025 20:52:57 -0500
Subject: [PATCH 04/16] Implement
 `foreach_pair_of_maybe_intersecting_edges_in_order`

with optional acceleration via trees and extent thinning
---
 src/methods/clipping/clipping_processor.jl | 410 +++++++++++++++++----
 src/methods/clipping/intersection.jl       |  24 +-
 2 files changed, 360 insertions(+), 74 deletions(-)

diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
index 94bf4e859..9f5d3dc8a 100644
--- a/src/methods/clipping/clipping_processor.jl
+++ b/src/methods/clipping/clipping_processor.jl
@@ -157,6 +157,264 @@ function _build_ab_list(alg::FosterHormannClipping, ::Type{T}, poly_a, poly_b, d
     return a_list, b_list, a_idx_list
 end
 
+"The number of vertices past which we should use a STRtree for edge intersection checking."
+const GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS = 32
+# Fallback convenience method so we can just pass the algorithm in
+function foreach_pair_of_maybe_intersecting_edges_in_order(
+    alg::FosterHormannClipping{M, A}, f_on_each_a::FA, f_after_each_a::FAAfter, f_on_each_maybe_intersect::FI, poly_a, poly_b, _t::Type{T} = Float64
+) where {FA, FAAfter, FI, T, M, A}
+    return foreach_pair_of_maybe_intersecting_edges_in_order(alg.manifold, alg.accelerator, f_on_each_a, f_after_each_a, f_on_each_maybe_intersect, poly_a, poly_b, T)
+end
+
+"""
+    foreach_pair_of_maybe_intersecting_edges_in_order(
+        manifold::M, accelerator::A,
+        f_on_each_a::FA,
+        f_after_each_a::FAAfter,
+        f_on_each_maybe_intersect::FI,
+        geom_a,
+        geom_b,
+        ::Type{T} = Float64
+    ) where {FA, FAAfter, FI, T, M <: Manifold, A <: IntersectionAccelerator}
+
+Decompose `geom_a` and `geom_b` into edge lists (unsorted), and then, logically, 
+perform the following iteration:
+
+```julia
+for (a_edge, i) in enumerate(eachedge(geom_a))
+    f_on_each_a(a_edge, i)
+    for (b_edge, j) in enumerate(eachedge(geom_b))
+        if may_intersect(a_edge, b_edge)
+            f_on_each_maybe_intersect(a_edge, b_edge)
+        end
+    end
+    f_after_each_a(a_edge, i)
+end
+```
+
+This may not be the exact acceleration that is performed - but it is 
+the logical sequence of events.  It also uses the `accelerator`, 
+and can automatically choose the best one based on an internal heuristic
+if you pass in an [`AutoAccelerator`](@ref).  
+
+For example, the `SingleSTRtree` accelerator is used along
+with extent thinning to avoid unnecessary edge intersection 
+checks in the inner loop.
+
+"""
+function foreach_pair_of_maybe_intersecting_edges_in_order(
+    manifold::M, accelerator::A, f_on_each_a::FA, f_after_each_a::FAAfter, f_on_each_maybe_intersect::FI, poly_a, poly_b, _t::Type{T} = Float64
+) where {FA, FAAfter, FI, T, M <: Manifold, A <: IntersectionAccelerator}
+    # TODO: dispatch on manifold
+    # this is suitable for planar
+    # but spherical / geodesic will need s2 support at some point,
+    # or -- even now -- just buffering
+    na = GI.npoint(poly_a)
+    nb = GI.npoint(poly_b)
+
+    accelerator = if accelerator isa AutoAccelerator
+        if na < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS && nb < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS
+            NestedLoop()
+        else
+            SingleSTRtree()
+        end
+    else
+        accelerator
+    end
+
+    if accelerator isa NestedLoop
+        # if we don't have enough vertices in either of the polygons to merit a tree,
+        # then we can just do a simple nested loop
+        # this becomes extremely useful in e.g. regridding, 
+        # where we know the polygon will only ever have a few vertices.
+        # This is also applicable to any manifold, since the checking is done within
+        # the loop.
+        # First, loop over "each edge" in poly_a
+        for (i, (a1t, a2t)) in enumerate(eachedge(poly_a, T))
+            a1t == a2t && continue
+            isnothing(f_on_each_a) ||f_on_each_a(a1t, i)
+            for (j, (b1t, b2t)) in enumerate(eachedge(poly_b, T))
+                b1t == b2t && continue
+                LoopStateMachine.@controlflow f_on_each_maybe_intersect(((a1t, a2t), i), ((b1t, b2t), j)) # this should be aware of manifold by construction.
+            end
+            isnothing(f_after_each_a) || f_after_each_a(a1t, i)
+        end
+        # And we're done!
+    elseif accelerator isa SingleSTRtree
+        # This is the "middle ground" case - run only a strtree 
+        # on poly_b without doing so on poly_a.
+        # This is less complex than running a dual tree traversal,
+        # and reduces the overhead of constructing an edge list and tree on poly_a.
+        ext_a, ext_b = GI.extent(poly_a), GI.extent(poly_b)
+        edges_b, indices_b = to_edgelist(ext_a, poly_b, T)
+        if isempty(edges_b) && !isnothing(f_on_each_a) && !isnothing(f_after_each_a)
+            # shortcut - nothing can possibly intersect
+            # so we just call f_on_each_a for each edge in poly_a
+            for i in 1:GI.npoint(poly_a)-1
+                pt = _tuple_point(GI.getpoint(poly_a, i), T)
+                f_on_each_a(pt, i)
+                f_after_each_a(pt, i)
+            end
+            return nothing
+        end
+
+        # This is the STRtree generated from the edges of poly_b
+        tree_b = STRtree(edges_b)
+
+        # this is a pre-allocation that will store the resuits of the query into tree_b
+        query_result = Int[] 
+        
+        # Loop over each vertex in poly_a
+        for (i, (a1t, a2t)) in enumerate(eachedge(poly_a, T))
+            a1t == a2t && continue
+            l1 = GI.Line(SVector{2}(a1t, a2t))
+            ext_l = GI.extent(l1)
+            # l = GI.Line(SVector{2}(a1t, a2t); extent=ext_l) # this seems to be unused - TODO remove
+            isnothing(f_on_each_a) || f_on_each_a(a1t, i)
+            # Query the STRtree for any edges in b that may intersect this edge
+            # This is sorted because we want to pretend we're doing the same thing
+            # as the nested loop above, and iterating through poly_b in order.
+            if Extents.intersects(ext_l, ext_b)
+                empty!(query_result)
+                SortTileRecursiveTree.query!(query_result, tree_b.rootnode, ext_l) # this is already sorted and uniqueified in STRtree.
+                # Loop over the edges in b that might intersect the edges in a
+                for j in query_result
+                    b1t, b2t = edges_b[j].geom
+                    b1t == b2t && continue
+                    # Manage control flow if the function returns a LoopStateMachine.Action
+                    # like Break(), Continue(), or Return()
+                    # This allows the function to break out of the loop early if it wants
+                    # without being syntactically inside the loop.
+                    LoopStateMachine.@controlflow f_on_each_maybe_intersect(((a1t, a2t), i), ((b1t, b2t), indices_b[j])) # note the indices_b[j] here - we are using the index of the edge in the original edge list, not the index of the edge in the STRtree.
+                end
+            end
+            isnothing(f_after_each_a) || f_after_each_a(a1t, i)
+        end
+    elseif accelerator isa SingleNaturalTree
+        ext_a, ext_b = GI.extent(poly_a), GI.extent(poly_b)
+        edges_b = to_edgelist(poly_b, T)
+
+        b_tree = NaturalIndexing.NaturalIndex(edges_b)
+
+        for (i, (a1t, a2t)) in enumerate(eachedge(poly_a, T))
+            a1t == a2t && continue
+            ext_l = Extents.Extent(X = minmax(a1t[1], a2t[1]), Y = minmax(a1t[2], a2t[2]))
+            isnothing(f_on_each_a) || f_on_each_a(a1t, i)
+            # Query the STRtree for any edges in b that may intersect this edge
+            # This is sorted because we want to pretend we're doing the same thing
+            # as the nested loop above, and iterating through poly_b in order.
+            if Extents.intersects(ext_l, ext_b)
+                # Loop over the edges in b that might intersect the edges in a
+                SpatialTreeInterface.do_query(Base.Fix1(Extents.intersects, ext_l), b_tree) do j
+                    b1t, b2t = edges_b[j].geom
+                    b1t == b2t && return LoopStateMachine.Continue()
+                    # LoopStateMachine control is managed outside the loop, by the do_query function.
+                    return f_on_each_maybe_intersect(((a1t, a2t), i), ((b1t, b2t), j)) # note the indices_b[j] here - we are using the index of the edge in the original edge list, not the index of the edge in the STRtree.
+                end
+            end
+        end
+
+    elseif accelerator isa DoubleNaturalTree
+        edges_a = to_edgelist(poly_a, T)
+        edges_b = to_edgelist(poly_b, T)
+
+        tree_a = NaturalIndexing.NaturalIndex(edges_a)
+        tree_b = NaturalIndexing.NaturalIndex(edges_b)
+
+        last_a_idx = 0
+
+        SpatialTreeInterface.do_dual_query(Extents.intersects, tree_a, tree_b) do a_edge_idx, b_edge_idx
+            a1t, a2t = edges_a[a_edge_idx].geom
+            b1t, b2t = edges_b[b_edge_idx].geom
+
+            if last_a_idx < a_edge_idx
+                if !isnothing(f_on_each_a)
+                    for i in (last_a_idx+1):(a_edge_idx-1)
+                        f_on_each_a((edges_a[i].geom[1]), i)
+                        !isnothing(f_after_each_a) && f_after_each_a((edges_a[i].geom[1]), i)
+                    end
+                end
+                !isnothing(f_on_each_a) && f_on_each_a(a1t, a_edge_idx)
+            end
+
+            f_on_each_maybe_intersect(((a1t, a2t), a_edge_idx), ((b1t, b2t), b_edge_idx))
+
+            if last_a_idx < a_edge_idx
+                if !isnothing(f_after_each_a)
+                    f_after_each_a(a1t, a_edge_idx)
+                end
+                last_a_idx = a_edge_idx
+            end
+        end
+
+        @show last_a_idx
+
+        if last_a_idx == 0 # the query did not find any intersections
+            if !isnothing(f_on_each_a) && isnothing(f_after_each_a)
+                return
+            else
+                for (i, edge) in enumerate(edges_a)
+                    !isnothing(f_on_each_a) && f_on_each_a(edge.geom[1], i)
+                    !isnothing(f_after_each_a) && f_after_each_a(edge.geom[1], i)
+                end
+            end
+        elseif last_a_idx < length(edges_a)
+            # the query terminated early - this will almost always be the case.
+            if !isnothing(f_on_each_a) && isnothing(f_after_each_a)
+                return
+            else
+                for (i, edge) in zip(last_a_idx+1:length(edges_a), view(edges_a, last_a_idx+1:length(edges_a)))
+                    !isnothing(f_on_each_a) && f_on_each_a(edge.geom[1], i)
+                    !isnothing(f_after_each_a) && f_after_each_a(edge.geom[1], i)
+                end
+            end
+        end
+    elseif accelerator isa ThinnedDoubleNaturalTree
+        ext_a, ext_b = GI.extent(poly_a), GI.extent(poly_b)
+        mutual_extent = Extents.intersection(ext_a, ext_b)
+
+        edges_a, indices_a = to_edgelist(mutual_extent, poly_a, T)
+        edges_b, indices_b = to_edgelist(mutual_extent, poly_b, T)
+
+        tree_a = NaturalIndexing.NaturalIndex(edges_a)
+        tree_b = NaturalIndexing.NaturalIndex(edges_b)
+
+        last_a_idx = 1
+
+        SpatialTreeInterface.do_dual_query(Extents.intersects, tree_a, tree_b) do a_thinned_idx, b_thinned_idx
+            a_edge_idx = indices_a[a_thinned_idx]
+            b_edge_idx = indices_b[b_thinned_idx]
+
+            a1t, a2t = edges_a[a_thinned_idx].geom
+            b1t, b2t = edges_b[b_thinned_idx].geom
+
+            if last_a_idx < a_edge_idx
+                if !isnothing(f_on_each_a)
+                    for i in last_a_idx:(a_edge_idx-1)
+                        f_on_each_a(a1t, a_edge_idx)
+                        !isnothing(f_after_each_a) && f_after_each_a(a1t, a_edge_idx)
+                    end
+                end
+                !isnothing(f_on_each_a) && f_on_each_a(a1t, a_edge_idx)
+            end
+
+            f_on_each_maybe_intersect(((a1t, a2t), a_edge_idx), ((b1t, b2t), b_edge_idx))
+
+            if last_a_idx < a_edge_idx
+                if !isnothing(f_after_each_a)
+                    f_after_each_a(a1t, a_edge_idx)
+                end
+                last_a_idx = a_edge_idx
+            end
+        end
+    else
+        error("Unsupported accelerator type: $accelerator.  FosterHormannClipping only supports NestedLoop() or SingleSTRtree().")
+    end
+
+    return nothing
+
+end
+
 #=
     _build_a_list(::Type{T}, poly_a, poly_b) -> (a_list, a_idx_list)
 
@@ -176,89 +434,101 @@ function _build_a_list(alg::FosterHormannClipping{M, A}, ::Type{T}, poly_a, poly
     a_list = PolyNode{T}[]  # list of points in poly_a
     sizehint!(a_list, n_a_edges)
     a_idx_list = Vector{Int}()  # finds indices of intersection points in a_list
-    a_count = 0  # number of points added to a_list
-    n_b_intrs = 0
-    # Loop through points of poly_a
-    local a_pt1
-    for (i, a_p2) in enumerate(GI.getpoint(poly_a))
-        a_pt2 = (T(GI.x(a_p2)), T(GI.y(a_p2)))
-        if i <= 1 || (a_pt1 == a_pt2)  # don't repeat points
-            a_pt1 = a_pt2
-            continue
-        end
-        # Add the first point of the edge to the list of points in a_list
-        new_point = PolyNode{T}(;point = a_pt1)
+    local a_count::Int = 0  # number of points added to a_list
+    local n_b_intrs::Int = 0
+    local prev_counter::Int = 0
+
+    function on_each_a(a_pt, i)
+        new_point = PolyNode{T}(;point = a_pt)
         a_count += 1
         push!(a_list, new_point)
-        # Find intersections with edges of poly_b
-        local b_pt1
         prev_counter = a_count
-        for (j, b_p2) in enumerate(GI.getpoint(poly_b))
-            b_pt2 = _tuple_point(b_p2, T)
-            if j <= 1 || (b_pt1 == b_pt2)  # don't repeat points
-                b_pt1 = b_pt2
-                continue
-            end
-            # Determine if edges intersect and how they intersect
-            line_orient, intr1, intr2 = _intersection_point(T, (a_pt1, a_pt2), (b_pt1, b_pt2); exact)
-            if line_orient != line_out  # edges intersect
-                if line_orient == line_cross  # Intersection point that isn't a vertex
-                    int_pt, fracs = intr1
+        return nothing
+    end
+
+    function after_each_a(a_pt, i)
+        # Order intersection points by placement along edge using fracs value
+        if prev_counter < a_count
+            Δintrs = a_count - prev_counter
+            inter_points = @view a_list[(a_count - Δintrs + 1):a_count]
+            sort!(inter_points, by = x -> x.fracs[1])
+        end
+        return nothing
+    end
+
+    function on_each_maybe_intersect(((a_pt1, a_pt2), i), ((b_pt1, b_pt2), j))
+        if (b_pt1 == b_pt2)  # don't repeat points
+            b_pt1 = b_pt2
+            return
+        end
+        # Determine if edges intersect and how they intersect
+        line_orient, intr1, intr2 = _intersection_point(alg.manifold, T, (a_pt1, a_pt2), (b_pt1, b_pt2); exact)
+        if line_orient != line_out  # edges intersect
+            if line_orient == line_cross  # Intersection point that isn't a vertex
+                int_pt, fracs = intr1
+                new_intr = PolyNode{T}(;
+                    point = int_pt, inter = true, neighbor = j, # j is now equivalent to old j-1
+                    crossing = true, fracs = fracs,
+                )
+                a_count += 1
+                n_b_intrs += 1
+                push!(a_list, new_intr)
+                push!(a_idx_list, a_count)
+            else
+                (_, (α1, β1)) = intr1
+                # Determine if a1 or b1 should be added to a_list
+                add_a1 = α1 == 0 && 0 ≤ β1 < 1
+                a1_β = add_a1 ? β1 : zero(T)
+                add_b1 = β1 == 0 && 0 < α1 < 1
+                b1_α = add_b1 ? α1 : zero(T)
+                # If lines are collinear and overlapping, a second intersection exists
+                if line_orient == line_over
+                    (_, (α2, β2)) = intr2
+                    if α2 == 0 && 0 ≤ β2 < 1
+                        add_a1, a1_β = true, β2
+                    end
+                    if β2 == 0 && 0 < α2 < 1
+                        add_b1, b1_α = true, α2
+                    end
+                end
+                # Add intersection points determined above
+                if add_a1
+                    n_b_intrs += a1_β == 0 ? 0 : 1
+                    a_list[prev_counter] = PolyNode{T}(;
+                        point = a_pt1, inter = true, neighbor = j,
+                        fracs = (zero(T), a1_β),
+                    )
+                    push!(a_idx_list, prev_counter)
+                end
+                if add_b1
                     new_intr = PolyNode{T}(;
-                        point = int_pt, inter = true, neighbor = j - 1,
-                        crossing = true, fracs = fracs,
+                        point = b_pt1, inter = true, neighbor = j,
+                        fracs = (b1_α, zero(T)),
                     )
                     a_count += 1
-                    n_b_intrs += 1
                     push!(a_list, new_intr)
                     push!(a_idx_list, a_count)
-                else
-                    (_, (α1, β1)) = intr1
-                    # Determine if a1 or b1 should be added to a_list
-                    add_a1 = α1 == 0 && 0 ≤ β1 < 1
-                    a1_β = add_a1 ? β1 : zero(T)
-                    add_b1 = β1 == 0 && 0 < α1 < 1
-                    b1_α = add_b1 ? α1 : zero(T)
-                    # If lines are collinear and overlapping, a second intersection exists
-                    if line_orient == line_over
-                        (_, (α2, β2)) = intr2
-                        if α2 == 0 && 0 ≤ β2 < 1
-                            add_a1, a1_β = true, β2
-                        end
-                        if β2 == 0 && 0 < α2 < 1
-                            add_b1, b1_α = true, α2
-                        end
-                    end
-                    # Add intersection points determined above
-                    if add_a1
-                        n_b_intrs += a1_β == 0 ? 0 : 1
-                        a_list[prev_counter] = PolyNode{T}(;
-                            point = a_pt1, inter = true, neighbor = j - 1,
-                            fracs = (zero(T), a1_β),
-                        )
-                        push!(a_idx_list, prev_counter)
-                    end
-                    if add_b1
-                        new_intr = PolyNode{T}(;
-                            point = b_pt1, inter = true, neighbor = j - 1,
-                            fracs = (b1_α, zero(T)),
-                        )
-                        a_count += 1
-                        push!(a_list, new_intr)
-                        push!(a_idx_list, a_count)
-                    end
                 end
             end
-            b_pt1 = b_pt2
         end
-        # Order intersection points by placement along edge using fracs value
-        if prev_counter < a_count
-            Δintrs = a_count - prev_counter
-            inter_points = @view a_list[(a_count - Δintrs + 1):a_count]
-            sort!(inter_points, by = x -> x.fracs[1])
+        return nothing
+    end
+
+    # do the iteration but in an accelerated way
+    # this is equivalent to (but faster than)
+    #=
+    ```julia
+    for ((a1, a2), i) in eachedge(poly_a)
+        on_each_a(a1, i)
+        for ((b1, b2), j) in eachedge(poly_b)
+            on_each_maybe_intersect(((a1, a2), i), ((b1, b2), j))
         end
-        a_pt1 = a_pt2
+        after_each_a(a1, i)
     end
+    ```
+    =#
+    foreach_pair_of_maybe_intersecting_edges_in_order(alg, on_each_a, after_each_a, on_each_maybe_intersect, poly_a, poly_b, T)
+
     return a_list, a_idx_list, n_b_intrs
 end
 
diff --git a/src/methods/clipping/intersection.jl b/src/methods/clipping/intersection.jl
index 47d82077c..b58198058 100644
--- a/src/methods/clipping/intersection.jl
+++ b/src/methods/clipping/intersection.jl
@@ -236,11 +236,13 @@ function _intersection_points(manifold::M, accelerator::A, ::Type{T}, ::GI.Abstr
     # Check if the geometries extents even overlap
     Extents.intersects(GI.extent(a), GI.extent(b)) || return result
     # Create a list of edges from the two input geometries
-    edges_a, edges_b = map(sort! ∘ to_edges, (a, b))
+    # edges_a, edges_b = map(sort! ∘ to_edges, (a, b))
     # Loop over pairs of edges and add any unique intersection points to results
-    for a_edge in edges_a, b_edge in edges_b
-        line_orient, intr1, intr2 = _intersection_point(T, a_edge, b_edge; exact)
-        line_orient == line_out && continue  # no intersection points
+    # TODO: add intersection acceleration here.
+
+    function f_on_each_maybe_intersect((a_edge, a_idx), (b_edge, b_idx))
+        line_orient, intr1, intr2 = _intersection_point(manifold, T, a_edge, b_edge; exact)
+        line_orient == line_out && return LoopStateMachine.Continue() # use LoopStateMachine.Continue() to skip this edge - in this case it doesn't matter but you could use it to e.g. break once you found the first intersecting point.
         pt1, _ = intr1
         push!(result, pt1)  # if not line_out, there is at least one intersection point
         if line_orient == line_over # if line_over, there are two intersection points
@@ -248,6 +250,20 @@ function _intersection_points(manifold::M, accelerator::A, ::Type{T}, ::GI.Abstr
             push!(result, pt2)
         end
     end
+
+    # iterate over each pair of intersecting edges only,
+    # calling `f_on_each_maybe_intersect` for each pair 
+    # that may intersect.
+    foreach_pair_of_maybe_intersecting_edges_in_order(
+        manifold, accelerator, 
+        nothing, # f_on_each_a
+        nothing, # f_after_each_a
+        f_on_each_maybe_intersect, # f_on_each_maybe_intersect
+        a,
+        b,
+        T
+    )
+    
     #= TODO: We might be able to just add unique points with checks on the α and β values
     returned from `_intersection_point`, but this would be different for curves vs polygons
     vs multipolygons depending on if the shape is closed. This then wouldn't allow using the

From e10351eb53b38c60e44b835da3a1712bf640792e Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Mon, 3 Mar 2025 20:54:06 -0500
Subject: [PATCH 05/16] Test all single accelerators in poly clipping tests

---
 test/methods/clipping/polygon_clipping.jl | 55 ++++++++++++-----------
 1 file changed, 30 insertions(+), 25 deletions(-)

diff --git a/test/methods/clipping/polygon_clipping.jl b/test/methods/clipping/polygon_clipping.jl
index 2fc4575fa..638b627bd 100644
--- a/test/methods/clipping/polygon_clipping.jl
+++ b/test/methods/clipping/polygon_clipping.jl
@@ -166,7 +166,7 @@ test_pairs = [
 const ϵ = 1e-10
 # Compare clipping results from GeometryOps and LibGEOS
 function compare_GO_LG_clipping(GO_f, LG_f, p1, p2)
-    GO_result_list = GO_f(p1, p2; target = GI.PolygonTrait())
+    
     LG_result_geom = LG_f(p1, p2)
     if LG_result_geom isa LG.GeometryCollection
         poly_list = LG.Polygon[]
@@ -175,38 +175,43 @@ function compare_GO_LG_clipping(GO_f, LG_f, p1, p2)
         end
         LG_result_geom = LG.MultiPolygon(poly_list)
     end
-    # Check if nothing is returned
-    if isempty(GO_result_list) && (LG.isEmpty(LG_result_geom) || LG.area(LG_result_geom) == 0)
-        return true
-    end
-    # Check for unnecessary points
-    if sum(GI.npoint, GO_result_list; init = 0.0) > GI.npoint(LG_result_geom)
-        return false
-    end
-    # Make sure last point is repeated
-    for poly in GO_result_list
-        for ring in GI.getring(poly)
-            GI.getpoint(ring, 1) != GI.getpoint(ring, GI.npoint(ring)) && return false
+
+    for _accelerator in (GO.AutoAccelerator(), GO.NestedLoop(), GO.SingleSTRtree(), GO.SingleNaturalTree(), #=GO.DoubleNaturalTree(), =# #=GO.ThinnedDoubleNaturalTree(), =# #=GO.DoubleSTRtree()=#)
+    @testset let accelerator = _accelerator # this is a ContextTestSet that is otherwise invisible but adds context to the testset
+        GO_result_list = GO_f(GO.FosterHormannClipping(accelerator), p1, p2; target = GI.PolygonTrait())
+        # Check if nothing is returned
+        if isempty(GO_result_list) && (LG.isEmpty(LG_result_geom) || LG.area(LG_result_geom) == 0)
+            @test true
+            continue
+        end
+        # Check for unnecessary points
+        @test !(sum(GI.npoint, GO_result_list; init = 0.0) > GI.npoint(LG_result_geom))
+        # Make sure last point is repeated
+        for poly in GO_result_list
+            for ring in GI.getring(poly)
+                @test !(GI.getpoint(ring, 1) != GI.getpoint(ring, GI.npoint(ring)))
+            end
         end
-    end
 
-    # Check if polygons cover the same area
-    local GO_result_geom
-    if length(GO_result_list) == 1
-        GO_result_geom = GO_result_list[1]
-    else
-        GO_result_geom = GI.MultiPolygon(GO_result_list)
-    end
-    diff_1_area = LG.area(LG.difference(GO_result_geom, LG_result_geom))
-    diff_2_area = LG.area(LG.difference(LG_result_geom, GO_result_geom))
-    return diff_1_area ≤ ϵ && diff_2_area ≤ ϵ
+        # Check if polygons cover the same area
+        local GO_result_geom
+        if length(GO_result_list) == 1
+            GO_result_geom = GO_result_list[1]
+        else
+            GO_result_geom = GI.MultiPolygon(GO_result_list)
+        end
+        diff_1_area = LG.area(LG.difference(GO_result_geom, LG_result_geom))
+        diff_2_area = LG.area(LG.difference(LG_result_geom, GO_result_geom))
+        @test diff_1_area ≤ ϵ && diff_2_area ≤ ϵ
+    end # testset
+    end # loop
 end
 
 # Test clipping functions and print error message if tests fail
 function test_clipping(GO_f, LG_f, f_name)
     for (p1, p2, sg1, sg2, sdesc) in test_pairs
         @testset_implementations "$sg1 $f_name $sg2 - $sdesc" begin
-            @test compare_GO_LG_clipping(GO_f, LG_f, $p1, $p2) 
+            compare_GO_LG_clipping(GO_f, LG_f, $p1, $p2) # this executes tests internally
         end
     end
     return

From 2077adfc9b8e0f1f6bc8817dd812b817d9597cd3 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Fri, 14 Mar 2025 11:30:06 -0800
Subject: [PATCH 06/16] Move imports to main file

Co-authored-by: Rafael Schouten <rafschouten@gmail.com>
---
 src/GeometryOps.jl                                     | 3 +++
 src/utils/NaturalIndexing.jl                           | 7 ++-----
 src/utils/SpatialTreeInterface/SpatialTreeInterface.jl | 2 +-
 3 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl
index d9c307e8f..23972c504 100644
--- a/src/GeometryOps.jl
+++ b/src/GeometryOps.jl
@@ -51,6 +51,9 @@ include("utils/NaturalIndexing.jl")
 using .NaturalIndexing
 
 
+# Load utility modules in
+using .NaturalIndexing, .SpatialTreeInterface, .LoopStateMachine
+
 include("methods/angles.jl")
 include("methods/area.jl")
 include("methods/barycentric.jl")
diff --git a/src/utils/NaturalIndexing.jl b/src/utils/NaturalIndexing.jl
index f972fe4ac..f698de5c3 100644
--- a/src/utils/NaturalIndexing.jl
+++ b/src/utils/NaturalIndexing.jl
@@ -4,6 +4,7 @@ import GeoInterface as GI
 import Extents
 
 using ..SpatialTreeInterface
+export NaturalTree, NaturallyIndexedRing, prepare_naturally
 
 """
     NaturalLevel{E <: Extents.Extent}
@@ -228,8 +229,4 @@ function prepare_naturally(geom)
     end
 end
 
-export NaturalTree, NaturallyIndexedRing, prepare_naturally
-
-end # module NaturalIndexing
-
-using .NaturalIndexing
\ No newline at end of file
+end # module NaturalIndexing
\ No newline at end of file
diff --git a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl
index d306d75b4..0066142f6 100644
--- a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl
+++ b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl
@@ -7,7 +7,7 @@ import GeoInterface as GI
 import AbstractTrees
 
 # public isspatialtree, isleaf, getchild, nchild, child_indices_extents, node_extent
-export query, do_query
+export query
 export FlatNoTree
 
 # The spatial tree interface and its implementations are defined here.

From 17b1487e1d734e9f79eceba36bc87007b0895f73 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Fri, 14 Mar 2025 11:32:03 -0800
Subject: [PATCH 07/16] Add comments to naturalindexing and clean up

e.g. remove redundant getpoint definitions, etc
---
 src/utils/NaturalIndexing.jl | 44 ++++++++++++++++++++++--------------
 1 file changed, 27 insertions(+), 17 deletions(-)

diff --git a/src/utils/NaturalIndexing.jl b/src/utils/NaturalIndexing.jl
index f698de5c3..92347db25 100644
--- a/src/utils/NaturalIndexing.jl
+++ b/src/utils/NaturalIndexing.jl
@@ -4,6 +4,9 @@ import GeoInterface as GI
 import Extents
 
 using ..SpatialTreeInterface
+
+import ..GeometryOps as GO # TODO: only needed for NaturallyIndexedRing, remove when that is removed.
+
 export NaturalTree, NaturallyIndexedRing, prepare_naturally
 
 """
@@ -35,33 +38,37 @@ struct NaturalIndex{E <: Extents.Extent}
     levels::Vector{NaturalLevel{E}}
 end
 
-GI.extent(idx::NaturalIndex) = idx.extent
 Extents.extent(idx::NaturalIndex) = idx.extent
 
 function Base.show(io::IO, ::MIME"text/plain", idx::NaturalIndex)
     println(io, "NaturalIndex with $(length(idx.levels)) levels and $(idx.nodecapacity) children per node")
     println(io, "extent: $(idx.extent)")
 end
-
 function Base.show(io::IO, idx::NaturalIndex)
     println(io, "NaturalIndex($(length(idx.levels)) levels, $(idx.extent))")
 end
 
 function NaturalIndex(geoms; nodecapacity = 32)
+    # Get the extent type initially (coord order, coord type, etc.)
+    # so that the construction is type stable.
     e1 = GI.extent(first(geoms))
     E = typeof(e1)
     return NaturalIndex{E}(geoms; nodecapacity = nodecapacity)
 end
-
-function NaturalIndex{E}(geoms; nodecapacity = 32) where E <: Extents.Extent
-    last_level_extents = GI.extent.(geoms)
+function NaturalIndex(last_level_extents::Vector{E}; nodecapacity = 32) where E <: Extents.Extent
+    # If we are passed a vector of extents - inflate immediately!
     return NaturalIndex{E}(last_level_extents; nodecapacity = nodecapacity)
 end
 
-function NaturalIndex(last_level_extents::Vector{E}; nodecapacity = 32) where E <: Extents.Extent
+function NaturalIndex{E}(geoms; nodecapacity = 32) where E <: Extents.Extent
+    # If passed a vector of geometries, and we know the type of the extent,
+    # then simply retrieve the extents so they can serve as the "last-level" 
+    # extents.
+    # Finally, call the lowest level method that performs inflation.
+    last_level_extents = GI.extent.(geoms)
     return NaturalIndex{E}(last_level_extents; nodecapacity = nodecapacity)
 end
-
+# This is the main constructor that performs inflation.
 function NaturalIndex{E}(last_level_extents::Vector{E}; nodecapacity = 32) where E <: Extents.Extent
     ngeoms = length(last_level_extents)
     last_level = NaturalLevel(last_level_extents)
@@ -70,11 +77,11 @@ function NaturalIndex{E}(last_level_extents::Vector{E}; nodecapacity = 32) where
 
     levels = Vector{NaturalLevel{E}}(undef, nlevels)
     levels[end] = last_level
-
+    # Iterate backwards, from bottom to top level,
+    # and build up the level extent vectors.
     for level_index in (nlevels-1):(-1):1
-        prev_level = levels[level_index+1] # this is always instantiated
+        prev_level = levels[level_index+1] # this is always instantiated, since we are iterating backwards
         nrects = _number_of_keys(nodecapacity, nlevels - (level_index), ngeoms)
-        # @show level_index nrects
         extents = [
             begin
                 start = (rect_index - 1) * nodecapacity + 1
@@ -185,6 +192,15 @@ SpatialTreeInterface.getchild(node::NaturalIndex, i) = SpatialTreeInterface.getc
 
 SpatialTreeInterface.child_indices_extents(node::NaturalIndex) = (i_ext for i_ext in enumerate(node.levels[1].extents))
 
+"""
+    NaturallyIndexedRing(points; nodecapacity = 32)
+
+A linear ring that contains a natural index.
+
+!!! warning
+    This will be removed in favour of prepared geometry - the idea here
+    is just to test what interface works best to store things in.
+"""
 struct NaturallyIndexedRing
     points::Vector{Tuple{Float64, Float64}}
     index::NaturalIndex{Extents.Extent{(:X, :Y), NTuple{2, NTuple{2, Float64}}}}
@@ -194,7 +210,6 @@ function NaturallyIndexedRing(points::Vector{Tuple{Float64, Float64}}; nodecapac
     index = NaturalIndex(GO.edge_extents(GI.LinearRing(points)); nodecapacity)
     return NaturallyIndexedRing(points, index)
 end
-
 NaturallyIndexedRing(ring::NaturallyIndexedRing) = ring
 
 function GI.convert(::Type{NaturallyIndexedRing}, ::GI.LinearRingTrait, geom)
@@ -203,24 +218,19 @@ function GI.convert(::Type{NaturallyIndexedRing}, ::GI.LinearRingTrait, geom)
 end
 
 Base.show(io::IO, ::MIME"text/plain", ring::NaturallyIndexedRing) = Base.show(io, ring)
-
 Base.show(io::IO, ring::NaturallyIndexedRing) = print(io, "NaturallyIndexedRing($(length(ring.points)) points) with index $(sprint(show, ring.index))")
 
 GI.ncoord(::GI.LinearRingTrait, ring::NaturallyIndexedRing) = 2
 GI.is3d(::GI.LinearRingTrait, ring::NaturallyIndexedRing) = false
 GI.ismeasured(::GI.LinearRingTrait, ring::NaturallyIndexedRing) = false
 
-GI.npoint(::GI.LinearRingTrait, ring::NaturallyIndexedRing) = length(ring.points)
-GI.getpoint(::GI.LinearRingTrait, ring::NaturallyIndexedRing) = ring.points[i]
-GI.getpoint(::GI.LinearRingTrait, ring::NaturallyIndexedRing, i::Int) = ring.points[i]
-
 GI.ngeom(::GI.LinearRingTrait, ring::NaturallyIndexedRing) = length(ring.points)
 GI.getgeom(::GI.LinearRingTrait, ring::NaturallyIndexedRing) = ring.points
 GI.getgeom(::GI.LinearRingTrait, ring::NaturallyIndexedRing, i::Int) = ring.points[i]
 
 Extents.extent(ring::NaturallyIndexedRing) = ring.index.extent
 
-GI.isgeometry(::NaturallyIndexedRing) = true
+GI.isgeometry(::Type{<: NaturallyIndexedRing}) = true
 GI.geomtrait(::NaturallyIndexedRing) = GI.LinearRingTrait()
 
 function prepare_naturally(geom)

From c6fc596cc8618f585e1ab92ea303db95aef83c92 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sat, 19 Apr 2025 20:10:08 -0400
Subject: [PATCH 08/16] Adapt to refactor of SpatialTreeInterface

do_query -> depth_first_search, etc.

Otherwise all is the same.
---
 src/methods/clipping/clipping_processor.jl | 8 ++++----
 src/methods/clipping/intersection.jl       | 2 +-
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
index 9f5d3dc8a..41b30e1de 100644
--- a/src/methods/clipping/clipping_processor.jl
+++ b/src/methods/clipping/clipping_processor.jl
@@ -305,10 +305,10 @@ function foreach_pair_of_maybe_intersecting_edges_in_order(
             # as the nested loop above, and iterating through poly_b in order.
             if Extents.intersects(ext_l, ext_b)
                 # Loop over the edges in b that might intersect the edges in a
-                SpatialTreeInterface.do_query(Base.Fix1(Extents.intersects, ext_l), b_tree) do j
+                SpatialTreeInterface.depth_first_search(Base.Fix1(Extents.intersects, ext_l), b_tree) do j
                     b1t, b2t = edges_b[j].geom
                     b1t == b2t && return LoopStateMachine.Continue()
-                    # LoopStateMachine control is managed outside the loop, by the do_query function.
+                    # LoopStateMachine control is managed outside the loop, by the depth_first_search function.
                     return f_on_each_maybe_intersect(((a1t, a2t), i), ((b1t, b2t), j)) # note the indices_b[j] here - we are using the index of the edge in the original edge list, not the index of the edge in the STRtree.
                 end
             end
@@ -323,7 +323,7 @@ function foreach_pair_of_maybe_intersecting_edges_in_order(
 
         last_a_idx = 0
 
-        SpatialTreeInterface.do_dual_query(Extents.intersects, tree_a, tree_b) do a_edge_idx, b_edge_idx
+        SpatialTreeInterface.dual_depth_first_search(Extents.intersects, tree_a, tree_b) do a_edge_idx, b_edge_idx
             a1t, a2t = edges_a[a_edge_idx].geom
             b1t, b2t = edges_b[b_edge_idx].geom
 
@@ -381,7 +381,7 @@ function foreach_pair_of_maybe_intersecting_edges_in_order(
 
         last_a_idx = 1
 
-        SpatialTreeInterface.do_dual_query(Extents.intersects, tree_a, tree_b) do a_thinned_idx, b_thinned_idx
+        SpatialTreeInterface.dual_depth_first_search(Extents.intersects, tree_a, tree_b) do a_thinned_idx, b_thinned_idx
             a_edge_idx = indices_a[a_thinned_idx]
             b_edge_idx = indices_b[b_thinned_idx]
 
diff --git a/src/methods/clipping/intersection.jl b/src/methods/clipping/intersection.jl
index b58198058..7af70d740 100644
--- a/src/methods/clipping/intersection.jl
+++ b/src/methods/clipping/intersection.jl
@@ -242,7 +242,7 @@ function _intersection_points(manifold::M, accelerator::A, ::Type{T}, ::GI.Abstr
 
     function f_on_each_maybe_intersect((a_edge, a_idx), (b_edge, b_idx))
         line_orient, intr1, intr2 = _intersection_point(manifold, T, a_edge, b_edge; exact)
-        line_orient == line_out && return LoopStateMachine.Continue() # use LoopStateMachine.Continue() to skip this edge - in this case it doesn't matter but you could use it to e.g. break once you found the first intersecting point.
+        line_orient == line_out && return LoopStateMachine.Action(:continue) # use LoopStateMachine.Continue() to skip this edge - in this case it doesn't matter but you could use it to e.g. break once you found the first intersecting point.
         pt1, _ = intr1
         push!(result, pt1)  # if not line_out, there is at least one intersection point
         if line_orient == line_over # if line_over, there are two intersection points

From 01c8deffa893548ed2006c219e82f39a4287c108 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sat, 19 Apr 2025 20:27:19 -0400
Subject: [PATCH 09/16] add `isspatialtree` to NaturalIndexing

---
 src/utils/NaturalIndexing.jl | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/src/utils/NaturalIndexing.jl b/src/utils/NaturalIndexing.jl
index 92347db25..131211894 100644
--- a/src/utils/NaturalIndexing.jl
+++ b/src/utils/NaturalIndexing.jl
@@ -152,6 +152,9 @@ Extents.extent(node::NaturalTreeNode) = node.extent
 # - `isleaf(node)` returns a boolean indicating whether the node is a leaf
 # - `child_indices_extents(node)` returns an iterator over the indices and extents of the children of the node
 
+SpatialTreeInterface.isspatialtree(::Type{<: NaturalIndexing}) = true
+SpatialTreeInterface.isspatialtree(::Type{<: NaturalTreeNode}) = true
+
 function SpatialTreeInterface.nchild(node::NaturalTreeNode)
     start_idx = (node.index - 1) * node.parent_index.nodecapacity + 1
     stop_idx = min(start_idx + node.parent_index.nodecapacity - 1, length(node.parent_index.levels[node.level+1].extents))

From f31a6daa802869d3e3afeead19e6a741c33283d4 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sat, 19 Apr 2025 20:28:00 -0400
Subject: [PATCH 10/16] add tests for NaturalIndex in SpatialTreeInterface
 tests

this probably isn't the best test suite but at least there is something..
---
 test/utils/SpatialTreeInterface.jl | 11 ++++++++++-
 1 file changed, 10 insertions(+), 1 deletion(-)

diff --git a/test/utils/SpatialTreeInterface.jl b/test/utils/SpatialTreeInterface.jl
index eb61f553b..3f78b05c3 100644
--- a/test/utils/SpatialTreeInterface.jl
+++ b/test/utils/SpatialTreeInterface.jl
@@ -4,6 +4,7 @@ using GeometryOps.SpatialTreeInterface
 using GeometryOps.SpatialTreeInterface: isspatialtree, isleaf, getchild, nchild, child_indices_extents, node_extent
 using GeometryOps.SpatialTreeInterface: query, depth_first_search, dual_depth_first_search
 using GeometryOps.SpatialTreeInterface: FlatNoTree
+using GeometryOps.NaturalIndexing: NaturalIndex
 using Extents
 using SortTileRecursiveTree: STRtree
 using NaturalEarth
@@ -218,7 +219,15 @@ end
     test_find_point_in_all_countries(STRtree)
 end
 
-
+# Test NaturalIndex implementation
+@testset "STRtree" begin
+    test_basic_interface(NaturalIndex)
+    test_child_indices_extents(NaturalIndex)
+    test_query_functionality(NaturalIndex)
+    test_dual_query_functionality(NaturalIndex)
+    test_geometry_support(NaturalIndex)
+    test_find_point_in_all_countries(NaturalIndex)
+end
 
 # This testset is not used because Polylabel.jl has some issues.
 

From 646619e6cd75ee8a5279868de1fbe1ee7678faf0 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Wed, 30 Apr 2025 22:33:51 -0400
Subject: [PATCH 11/16] Rename NaturalTree -> NaturalIndex and fix typos
 everywhere

---
 src/utils/NaturalIndexing.jl | 28 ++++++++++++++--------------
 1 file changed, 14 insertions(+), 14 deletions(-)

diff --git a/src/utils/NaturalIndexing.jl b/src/utils/NaturalIndexing.jl
index 131211894..080da277e 100644
--- a/src/utils/NaturalIndexing.jl
+++ b/src/utils/NaturalIndexing.jl
@@ -7,7 +7,7 @@ using ..SpatialTreeInterface
 
 import ..GeometryOps as GO # TODO: only needed for NaturallyIndexedRing, remove when that is removed.
 
-export NaturalTree, NaturallyIndexedRing, prepare_naturally
+export NaturalIndex, NaturallyIndexedRing, prepare_naturally
 
 """
     NaturalLevel{E <: Extents.Extent}
@@ -124,7 +124,7 @@ end
 
 # This is like a pointer to a node in the tree.
 """
-    NaturalTreeNode{E <: Extents.Extent}
+    NaturalIndexNode{E <: Extents.Extent}
 
 A reference to a node in the natural tree.  Kind of like a tree cursor.
 
@@ -133,14 +133,14 @@ A reference to a node in the natural tree.  Kind of like a tree cursor.
 - `index` is the index of the node in the level
 - `extent` is the extent of the node
 """
-struct NaturalTreeNode{E <: Extents.Extent}
+struct NaturalIndexNode{E <: Extents.Extent}
     parent_index::NaturalIndex{E}
     level::Int
     index::Int
     extent::E
 end
 
-Extents.extent(node::NaturalTreeNode) = node.extent
+Extents.extent(node::NaturalIndexNode) = node.extent
 
 # What does SpatialTreeInterface require of trees?
 # - Parents completely cover their children
@@ -152,18 +152,18 @@ Extents.extent(node::NaturalTreeNode) = node.extent
 # - `isleaf(node)` returns a boolean indicating whether the node is a leaf
 # - `child_indices_extents(node)` returns an iterator over the indices and extents of the children of the node
 
-SpatialTreeInterface.isspatialtree(::Type{<: NaturalIndexing}) = true
-SpatialTreeInterface.isspatialtree(::Type{<: NaturalTreeNode}) = true
+SpatialTreeInterface.isspatialtree(::Type{<: NaturalIndex}) = true
+SpatialTreeInterface.isspatialtree(::Type{<: NaturalIndexNode}) = true
 
-function SpatialTreeInterface.nchild(node::NaturalTreeNode)
+function SpatialTreeInterface.nchild(node::NaturalIndexNode)
     start_idx = (node.index - 1) * node.parent_index.nodecapacity + 1
     stop_idx = min(start_idx + node.parent_index.nodecapacity - 1, length(node.parent_index.levels[node.level+1].extents))
     return stop_idx - start_idx + 1
 end
 
-function SpatialTreeInterface.getchild(node::NaturalTreeNode, i::Int)
+function SpatialTreeInterface.getchild(node::NaturalIndexNode, i::Int)
     child_index = (node.index - 1) * node.parent_index.nodecapacity + i
-    return NaturalTreeNode(
+    return NaturalIndexNode(
         node.parent_index, 
         node.level + 1, # increment level by 1
         child_index, # index of this particular child
@@ -172,13 +172,13 @@ function SpatialTreeInterface.getchild(node::NaturalTreeNode, i::Int)
 end
 
 # Get all children of a node
-function SpatialTreeInterface.getchild(node::NaturalTreeNode)
+function SpatialTreeInterface.getchild(node::NaturalIndexNode)
     return (SpatialTreeInterface.getchild(node, i) for i in 1:SpatialTreeInterface.nchild(node))
 end
 
-SpatialTreeInterface.isleaf(node::NaturalTreeNode) = node.level == length(node.parent_index.levels) - 1
+SpatialTreeInterface.isleaf(node::NaturalIndexNode) = node.level == length(node.parent_index.levels) - 1
 
-function SpatialTreeInterface.child_indices_extents(node::NaturalTreeNode)
+function SpatialTreeInterface.child_indices_extents(node::NaturalIndexNode)
     start_idx = (node.index - 1) * node.parent_index.nodecapacity + 1
     stop_idx = min(start_idx + node.parent_index.nodecapacity - 1, length(node.parent_index.levels[node.level+1].extents))
     return ((i, node.parent_index.levels[node.level+1].extents[i]) for i in start_idx:stop_idx)
@@ -190,8 +190,8 @@ SpatialTreeInterface.isleaf(node::NaturalIndex) = length(node.levels) == 1
 
 SpatialTreeInterface.nchild(node::NaturalIndex) = length(node.levels[1].extents)
 
-SpatialTreeInterface.getchild(node::NaturalIndex) = SpatialTreeInterface.getchild(NaturalTreeNode(node, 0, 1, node.extent))
-SpatialTreeInterface.getchild(node::NaturalIndex, i) = SpatialTreeInterface.getchild(NaturalTreeNode(node, 0, 1, node.extent), i)
+SpatialTreeInterface.getchild(node::NaturalIndex) = SpatialTreeInterface.getchild(NaturalIndexNode(node, 0, 1, node.extent))
+SpatialTreeInterface.getchild(node::NaturalIndex, i) = SpatialTreeInterface.getchild(NaturalIndexNode(node, 0, 1, node.extent), i)
 
 SpatialTreeInterface.child_indices_extents(node::NaturalIndex) = (i_ext for i_ext in enumerate(node.levels[1].extents))
 

From 6b5e0873d43595a2dd04727b851dc4d39e9b8d48 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Fri, 9 May 2025 20:02:34 -0400
Subject: [PATCH 12/16] Fix minor stupidity in singlestr implementation

---
 src/methods/clipping/clipping_processor.jl | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
index 41b30e1de..9791ef770 100644
--- a/src/methods/clipping/clipping_processor.jl
+++ b/src/methods/clipping/clipping_processor.jl
@@ -157,6 +157,7 @@ function _build_ab_list(alg::FosterHormannClipping, ::Type{T}, poly_a, poly_b, d
     return a_list, b_list, a_idx_list
 end
 
+
 "The number of vertices past which we should use a STRtree for edge intersection checking."
 const GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS = 32
 # Fallback convenience method so we can just pass the algorithm in
@@ -215,8 +216,10 @@ function foreach_pair_of_maybe_intersecting_edges_in_order(
     accelerator = if accelerator isa AutoAccelerator
         if na < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS && nb < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS
             NestedLoop()
+        elseif na < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS || nb < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS
+            SingleNaturalTree()
         else
-            SingleSTRtree()
+            DoubleNaturalTree()
         end
     else
         accelerator
@@ -276,7 +279,8 @@ function foreach_pair_of_maybe_intersecting_edges_in_order(
             # as the nested loop above, and iterating through poly_b in order.
             if Extents.intersects(ext_l, ext_b)
                 empty!(query_result)
-                SortTileRecursiveTree.query!(query_result, tree_b.rootnode, ext_l) # this is already sorted and uniqueified in STRtree.
+                SortTileRecursiveTree.query!(query_result, tree_b.rootnode, ext_l) # this is already sorted and uniqueified in STRtree.'
+                sort!(query_result)
                 # Loop over the edges in b that might intersect the edges in a
                 for j in query_result
                     b1t, b2t = edges_b[j].geom
@@ -312,6 +316,7 @@ function foreach_pair_of_maybe_intersecting_edges_in_order(
                     return f_on_each_maybe_intersect(((a1t, a2t), i), ((b1t, b2t), j)) # note the indices_b[j] here - we are using the index of the edge in the original edge list, not the index of the edge in the STRtree.
                 end
             end
+            isnothing(f_after_each_a) || f_after_each_a(a1t, i)
         end
 
     elseif accelerator isa DoubleNaturalTree
@@ -347,8 +352,6 @@ function foreach_pair_of_maybe_intersecting_edges_in_order(
             end
         end
 
-        @show last_a_idx
-
         if last_a_idx == 0 # the query did not find any intersections
             if !isnothing(f_on_each_a) && isnothing(f_after_each_a)
                 return

From aea8ed15fe680ac645afcfdae47e729a954022e8 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Fri, 9 May 2025 20:04:14 -0400
Subject: [PATCH 13/16] make it easier to tell where errors are happening

---
 test/methods/clipping/polygon_clipping.jl | 13 ++++++++++---
 1 file changed, 10 insertions(+), 3 deletions(-)

diff --git a/test/methods/clipping/polygon_clipping.jl b/test/methods/clipping/polygon_clipping.jl
index 638b627bd..b8f6e5e34 100644
--- a/test/methods/clipping/polygon_clipping.jl
+++ b/test/methods/clipping/polygon_clipping.jl
@@ -200,9 +200,16 @@ function compare_GO_LG_clipping(GO_f, LG_f, p1, p2)
         else
             GO_result_geom = GI.MultiPolygon(GO_result_list)
         end
-        diff_1_area = LG.area(LG.difference(GO_result_geom, LG_result_geom))
-        diff_2_area = LG.area(LG.difference(LG_result_geom, GO_result_geom))
-        @test diff_1_area ≤ ϵ && diff_2_area ≤ ϵ
+        
+        diff_1_poly = @test_nowarn LG.difference(GO_result_geom, LG_result_geom)
+        diff_2_poly = @test_nowarn LG.difference(LG_result_geom, GO_result_geom)
+        if GI.isempty(diff_1_poly) || GI.isempty(diff_2_poly)
+            @test true # if both differences are empty, then the areas are the same by default.
+        else
+            diff_1_area = LG.area(diff_1_poly)
+            diff_2_area = LG.area(diff_2_poly)
+            @test diff_1_area ≤ ϵ && diff_2_area ≤ ϵ
+        end
     end # testset
     end # loop
 end

From 23dedecffd09f7e39d56d29955f5fe35f19d5b35 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Fri, 9 May 2025 20:27:11 -0400
Subject: [PATCH 14/16] Fix doc build error

---
 src/methods/clipping/clipping_processor.jl | 16 ++++++++++++----
 1 file changed, 12 insertions(+), 4 deletions(-)

diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
index 9791ef770..d93008d6e 100644
--- a/src/methods/clipping/clipping_processor.jl
+++ b/src/methods/clipping/clipping_processor.jl
@@ -29,6 +29,14 @@ struct DoubleSTRtree <: IntersectionAccelerator end
 struct SingleNaturalTree <: IntersectionAccelerator end
 struct DoubleNaturalTree <: IntersectionAccelerator end
 struct ThinnedDoubleNaturalTree <: IntersectionAccelerator end
+
+"""
+    AutoAccelerator()
+
+Let the algorithm choose the best accelerator based on the size of the input polygons.
+
+Once we have prepared geometry, this will also consider the existing preparations on the geoms.
+"""
 struct AutoAccelerator <: IntersectionAccelerator end
 
 """
@@ -279,8 +287,8 @@ function foreach_pair_of_maybe_intersecting_edges_in_order(
             # as the nested loop above, and iterating through poly_b in order.
             if Extents.intersects(ext_l, ext_b)
                 empty!(query_result)
-                SortTileRecursiveTree.query!(query_result, tree_b.rootnode, ext_l) # this is already sorted and uniqueified in STRtree.'
-                sort!(query_result)
+                SortTileRecursiveTree.query!(query_result, tree_b.rootnode, ext_l)
+                sort!(query_result) # STRTree.jl's query! does not sort!, even though query does...
                 # Loop over the edges in b that might intersect the edges in a
                 for j in query_result
                     b1t, b2t = edges_b[j].geom
@@ -326,7 +334,7 @@ function foreach_pair_of_maybe_intersecting_edges_in_order(
         tree_a = NaturalIndexing.NaturalIndex(edges_a)
         tree_b = NaturalIndexing.NaturalIndex(edges_b)
 
-        last_a_idx = 0
+        last_a_idx::Int = 0
 
         SpatialTreeInterface.dual_depth_first_search(Extents.intersects, tree_a, tree_b) do a_edge_idx, b_edge_idx
             a1t, a2t = edges_a[a_edge_idx].geom
@@ -382,7 +390,7 @@ function foreach_pair_of_maybe_intersecting_edges_in_order(
         tree_a = NaturalIndexing.NaturalIndex(edges_a)
         tree_b = NaturalIndexing.NaturalIndex(edges_b)
 
-        last_a_idx = 1
+        last_a_idx::Int = 1
 
         SpatialTreeInterface.dual_depth_first_search(Extents.intersects, tree_a, tree_b) do a_thinned_idx, b_thinned_idx
             a_edge_idx = indices_a[a_thinned_idx]

From ac75de81711649582fc92cb7d022c41b5a932fe6 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Fri, 9 May 2025 22:26:59 -0400
Subject: [PATCH 15/16] fix multiple type declarations

---
 src/methods/clipping/clipping_processor.jl | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
index 9791ef770..b72d55841 100644
--- a/src/methods/clipping/clipping_processor.jl
+++ b/src/methods/clipping/clipping_processor.jl
@@ -224,6 +224,9 @@ function foreach_pair_of_maybe_intersecting_edges_in_order(
     else
         accelerator
     end
+    
+    # Declare this now because it's used in multiple branches
+    last_a_idx::Int = 0
 
     if accelerator isa NestedLoop
         # if we don't have enough vertices in either of the polygons to merit a tree,
@@ -382,7 +385,7 @@ function foreach_pair_of_maybe_intersecting_edges_in_order(
         tree_a = NaturalIndexing.NaturalIndex(edges_a)
         tree_b = NaturalIndexing.NaturalIndex(edges_b)
 
-        last_a_idx = 1
+        last_a_idx = 0
 
         SpatialTreeInterface.dual_depth_first_search(Extents.intersects, tree_a, tree_b) do a_thinned_idx, b_thinned_idx
             a_edge_idx = indices_a[a_thinned_idx]

From f4a1077fa1bab33b3dea6354a2b303b875671fc6 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sat, 10 May 2025 18:51:39 -0400
Subject: [PATCH 16/16] fix dual type declarations, again

---
 src/GeometryOps.jl                         | 2 +-
 src/methods/clipping/clipping_processor.jl | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl
index 38c481a97..95e84788c 100644
--- a/src/GeometryOps.jl
+++ b/src/GeometryOps.jl
@@ -21,7 +21,7 @@ export TraitTarget, Manifold, Planar, Spherical, Geodesic, apply, applyreduce, f
 using GeoInterface
 using LinearAlgebra, Statistics
 
-using GeometryBasics.StaticArrays
+using StaticArrays
 
 import Tables, DataAPI
 import StaticArrays
diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
index 7b5690e46..d784984bd 100644
--- a/src/methods/clipping/clipping_processor.jl
+++ b/src/methods/clipping/clipping_processor.jl
@@ -337,7 +337,7 @@ function foreach_pair_of_maybe_intersecting_edges_in_order(
         tree_a = NaturalIndexing.NaturalIndex(edges_a)
         tree_b = NaturalIndexing.NaturalIndex(edges_b)
 
-        last_a_idx::Int = 0
+        last_a_idx = 0
 
         SpatialTreeInterface.dual_depth_first_search(Extents.intersects, tree_a, tree_b) do a_edge_idx, b_edge_idx
             a1t, a2t = edges_a[a_edge_idx].geom