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

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

diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl
index 6ca4923c8d..abe7d6710c 100644
--- a/src/GeometryOps.jl
+++ b/src/GeometryOps.jl
@@ -21,12 +21,16 @@ export TraitTarget, Manifold, Planar, Spherical, Geodesic, apply, applyreduce, f
 using GeoInterface
 using LinearAlgebra, Statistics
 
+using GeometryBasics.StaticArrays
+
 import Tables, DataAPI
 import 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
 

From 837b080e81033ed3d9f3551ea71650dbc32a68ee Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Mon, 3 Mar 2025 20:52:03 -0500
Subject: [PATCH 02/24] 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 abe7d6710c..3cd1fe9a58 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 0000000000..f972fe4acd
--- /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 23872916c88f5dcb7396a281c5efd9e671c6962b Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Mon, 3 Mar 2025 20:52:57 -0500
Subject: [PATCH 03/24] 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 94bf4e859a..9f5d3dc8ae 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 47d82077c4..b58198058c 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 2932b2a407c86837b594deae220885c576fe63f2 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Mon, 3 Mar 2025 20:54:06 -0500
Subject: [PATCH 04/24] 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 2fc4575fa3..638b627bd7 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 c002b6e11dbe98c0eed8f4b12e59a8b7802cb39d Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Fri, 14 Mar 2025 11:30:06 -0800
Subject: [PATCH 05/24] 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 3cd1fe9a58..38c481a975 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 f972fe4acd..f698de5c33 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 d306d75b45..0066142f6f 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 c8fe71a15848ceae1923b41ac9d3f8030f5d9556 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Fri, 14 Mar 2025 11:32:03 -0800
Subject: [PATCH 06/24] 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 f698de5c33..92347db25e 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 5bf118486536218c3572926a210ac2e995f7bf85 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sat, 19 Apr 2025 20:10:08 -0400
Subject: [PATCH 07/24] 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 9f5d3dc8ae..41b30e1de0 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 b58198058c..7af70d7401 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 73d66bf141fb97d96f8c440204bcaa59f1415918 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sat, 19 Apr 2025 20:27:19 -0400
Subject: [PATCH 08/24] 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 92347db25e..131211894f 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 7c54be879e4466c60c887c0236f197390b4a91cd Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sat, 19 Apr 2025 20:28:00 -0400
Subject: [PATCH 09/24] 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 eb61f553b7..3f78b05c33 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 b2bb3be0287ce2c20e6bc8f73704ceb2a8fc5bbb Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Wed, 30 Apr 2025 22:33:51 -0400
Subject: [PATCH 10/24] 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 131211894f..080da277e4 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 52d42e9edb53dae686b158c6f3ac191668f691af Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Fri, 9 May 2025 20:02:34 -0400
Subject: [PATCH 11/24] 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 41b30e1de0..9791ef7702 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 1c24e837839e20e7bc9c2178186db713bf25347e Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Fri, 9 May 2025 20:04:14 -0400
Subject: [PATCH 12/24] 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 638b627bd7..b8f6e5e344 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 56348cf1081a240e6c2a607b0bf8ca07b4528778 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Fri, 9 May 2025 22:26:59 -0400
Subject: [PATCH 13/24] 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 9791ef7702..b72d558416 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 8308ae7d5d55eba45354216ed5b7726767268455 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/24] 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 b72d558416..4ec7561f5f 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
 
 """
@@ -282,8 +290,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
@@ -329,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 = 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
@@ -385,7 +393,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 = 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 74687967fc818d38abe15bec70b11acf0d3c85e5 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sat, 10 May 2025 18:51:39 -0400
Subject: [PATCH 15/24] 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 38c481a975..95e84788c0 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 4ec7561f5f..33b0a3a9ec 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

From 8a5a12f3a97843101b05738a32e1519c3d10c3eb Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Tue, 22 Jul 2025 23:10:41 -0400
Subject: [PATCH 16/24] actually fix dual type declarations

---
 src/methods/clipping/clipping_processor.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
index 33b0a3a9ec..6b09fe2d66 100644
--- a/src/methods/clipping/clipping_processor.jl
+++ b/src/methods/clipping/clipping_processor.jl
@@ -234,7 +234,7 @@ function foreach_pair_of_maybe_intersecting_edges_in_order(
     end
     
     # Declare this now because it's used in multiple branches
-    last_a_idx::Int = 0
+    last_a_idx = 0
 
     if accelerator isa NestedLoop
         # if we don't have enough vertices in either of the polygons to merit a tree,

From 556372cd59fa04001031821d532553cb8f81daa5 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Fri, 25 Jul 2025 23:23:54 -0400
Subject: [PATCH 17/24] add spatial index building gif

---
 .../src/experiments/spatial_index_building.jl | 49 +++++++++++++++++++
 1 file changed, 49 insertions(+)
 create mode 100644 docs/src/experiments/spatial_index_building.jl

diff --git a/docs/src/experiments/spatial_index_building.jl b/docs/src/experiments/spatial_index_building.jl
new file mode 100644
index 0000000000..ddb8b4c723
--- /dev/null
+++ b/docs/src/experiments/spatial_index_building.jl
@@ -0,0 +1,49 @@
+using CairoMakie
+import GeoInterface as GI, GeometryOps as GO
+using NaturalEarth
+
+
+using GeometryOps.SpatialTreeInterface
+import GeometryOps.SpatialTreeInterface as STI
+using SortTileRecursiveTree
+
+function build_spatial_index_gif(geom, index_constructor, filename; plot_leaves = true, title = splitext(filename)[1], axis = (;), figure = (;), record = (;))
+    fig = Figure(; figure...)
+    ax = Axis(fig[1, 1]; title = title, axis...)
+
+    # Create a spatial index
+    index = index_constructor(geom)
+    ext = STI.node_extent(index)
+    limits!(ax, ext.X[1], ext.X[2], ext.Y[1], ext.Y[2])
+
+    rects = Rect2f[Rect2f((NaN, NaN), (NaN, NaN))]
+    colors = RGBAf[to_color(:transparent)]
+    palette = Makie.wong_colors(0.7)
+
+    plt = poly!(ax, rects; color = colors)
+
+    to_rect2(extent) = Rect2f((extent.X[1], extent.Y[1]), (extent.X[2] - extent.X[1], extent.Y[2] - extent.Y[1]))
+
+    function dive_in(io, plt, node, level)
+        if STI.isleaf(node) && plot_leaves
+            push!(rects, to_rect2(STI.node_extent(node)))
+            push!(colors, palette[level])
+        else
+            for child in STI.getchild(node)
+                dive_in(io, plt, child, level + 1)
+            end
+            push!(rects, to_rect2(STI.node_extent(node)))
+            push!(colors, palette[level])
+        end
+        update!(plt, rects; color = colors)
+        recordframe!(io)
+        return
+    end
+
+    Makie.record(fig, filename; record...) do io
+        empty!(rects)
+        empty!(colors)
+        dive_in(io, plt, index, 1)
+    end
+
+end 
\ No newline at end of file

From d0dfef4493ee3b86140a0e8a45ac860eb56b8edd Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sat, 26 Jul 2025 01:05:07 -0400
Subject: [PATCH 18/24] Rationalize test helper codegen

---
 test/helpers.jl | 114 +++++++++++++++++++++++-------------------------
 1 file changed, 55 insertions(+), 59 deletions(-)

diff --git a/test/helpers.jl b/test/helpers.jl
index af7b86e086..32cf3870d8 100644
--- a/test/helpers.jl
+++ b/test/helpers.jl
@@ -8,71 +8,76 @@ export @test_implementations, @testset_implementations
 
 const TEST_MODULES = [GeoInterface, ArchGDAL, GeometryBasics, LibGEOS]
 
+function conversion_expr(mod, var, genkey)
+    quote
+        $genkey = if $var isa $(GeoInterface.Extents.Extent)
+            if $mod in ($GeoInterface, $LibGEOS)
+                $var
+            else
+                $GeoInterface.convert($mod, $(GO.extent_to_polygon)($var))
+            end
+        # GeometryBasics does not have a Line geometry type.
+        # elseif $mod in ($GeometryBasics,) && $GeoInterface.trait($var) isa $GeoInterface.LineTrait
+        #     $var
+        elseif $mod in ($GeoInterface, $ArchGDAL, $GeometryBasics) && $GeoInterface.isempty($var)
+            $var
+        else
+            $GeoInterface.convert($mod, $var)
+        end
+    end
+end
 # Monkey-patch GeometryBasics to have correct methods.
 # TODO: push this up to GB!
 
-@eval GeometryBasics begin
-    # MultiGeometry ncoord implementations
-    GeoInterface.ncoord(::GeoInterface.MultiPolygonTrait, ::GeometryBasics.MultiPolygon{N}) where N = N
-    GeoInterface.ncoord(::GeoInterface.MultiLineStringTrait, ::GeometryBasics.MultiLineString{N}) where N = N
-    GeoInterface.ncoord(::GeoInterface.MultiPointTrait, ::GeometryBasics.MultiPoint{N}) where N = N
-    # LinearRing and LineString confusion
-    GeometryBasics.geointerface_geomtype(::GeoInterface.LinearRingTrait) = LineString
-    function GeoInterface.convert(::Type{LineString}, ::GeoInterface.LinearRingTrait, geom)
-        return GeoInterface.convert(LineString, GeoInterface.LineStringTrait(), geom) # forward to the linestring conversion method
-    end
-    # Line interface
-    GeometryBasics.geointerface_geomtype(::GeoInterface.LineTrait) = Line
-    function GeoInterface.convert(::Type{Line}, ::GeoInterface.LineTrait, geom)
-        p1, p2 = GeoInterface.getpoint(geom)
-        return Line(GeoInterface.convert(Point, GeoInterface.PointTrait(), p1), GeoInterface.convert(Point, GeoInterface.PointTrait(), p2))
+    # TODO: remove when GB GI pr lands
+    @static if hasmethod(GeometryBasics.convert, (Type{GeometryBasics.LineString}, GeoInterface.LinearRingTrait, Any))
+        function GeoInterface.convert(
+            ::Type{GeometryBasics.LineString}, 
+            ::GeoInterface.LinearRingTrait, 
+            geom
+            )
+            return GeoInterface.convert(GeometryBasics.LineString, GeoInterface.LineStringTrait(), geom)
+        end
+        GeometryBasics.geointerface_geomtype(::GeoInterface.LinearRingTrait) = GeometryBasics.LineString
     end
+    # end todo
     # GeometryCollection interface - currently just a large Union
-    const _ALL_GB_GEOM_TYPES = Union{Point, Line, LineString, Polygon, MultiPolygon, MultiLineString, MultiPoint}
+    const _ALL_GB_GEOM_TYPES = Union{GeometryBasics.Point, GeometryBasics.LineString, GeometryBasics.Polygon, GeometryBasics.MultiPolygon, GeometryBasics.MultiLineString, GeometryBasics.MultiPoint}
     GeometryBasics.geointerface_geomtype(::GeoInterface.GeometryCollectionTrait) = Vector{_ALL_GB_GEOM_TYPES}
-    function GeoInterface.convert(::Type{Vector{_ALL_GB_GEOM_TYPES}}, ::GeoInterface.GeometryCollectionTrait, geoms)
+    function GeoInterface.convert(::Type{Vector{<: _ALL_GB_GEOM_TYPES}}, ::GeoInterface.GeometryCollectionTrait, geoms)
         return _ALL_GB_GEOM_TYPES[GeoInterface.convert(GeometryBasics, g) for g in GeoInterface.getgeom(geoms)]
     end
 
     function GeoInterface.convert(
         ::Type{GeometryBasics.LineString}, 
         type::GeoInterface.LineStringTrait, 
-        geom::GeoInterface.Wrappers.LinearRing{false, false, StaticArraysCore.SVector{N, Tuple{Float64, Float64}}, Nothing, Nothing} where N
+        geom::GeoInterface.Wrappers.LinearRing{false, false, GO.StaticArrays.SVector{N, Tuple{Float64, Float64}}, Nothing, Nothing} where N
         )
-        return GeoInterface.convert(LineString, GeoInterface.LineStringTrait(), collect(geom.geom))
+        return GeometryBasics.LineString(GeometryBasics.Point2{Float64}.(collect(geom.geom)))
     end
 
-    function GeoInterface.convert(
-        ::Type{GeometryBasics.LineString}, 
-        type::GeoInterface.LineStringTrait, 
-        geom::GeoInterface.Wrappers.LinearRing{false, false, StaticArraysCore.SVector{N, Tuple{Float64, Float64}}, Nothing, Nothing} where N
-        )
-        return LineString(Point2{Float64}.(collect(geom.geom)))
-    end
-end
-
 
-@eval ArchGDAL begin
-    function GeoInterface.convert(
-        ::Type{T},
-        type::GeoInterface.PolygonTrait,
-        geom,
-    ) where {T<:IGeometry}
-        f = get(lookup_method, typeof(type), nothing)
-        isnothing(f) && error(
-            "Cannot convert an object of $(typeof(geom)) with the $(typeof(type)) trait (yet). Please report an issue.",
-        )
-        poly = createpolygon()
-        foreach(GeoInterface.getring(geom)) do ring
-            xs = GeoInterface.x.(GeoInterface.getpoint(ring)) |> collect
-            ys = GeoInterface.y.(GeoInterface.getpoint(ring)) |> collect
-            subgeom = unsafe_createlinearring(xs, ys)
-            result = GDAL.ogr_g_addgeometrydirectly(poly, subgeom)
-            @ogrerr result "Failed to add linearring."
+    @eval ArchGDAL begin
+        function GeoInterface.convert(
+            ::Type{T},
+            type::GeoInterface.PolygonTrait,
+            geom,
+        ) where {T<:IGeometry}
+            f = get(lookup_method, typeof(type), nothing)
+            isnothing(f) && error(
+                "Cannot convert an object of $(typeof(geom)) with the $(typeof(type)) trait (yet). Please report an issue.",
+            )
+            poly = createpolygon()
+            foreach(GeoInterface.getring(geom)) do ring
+                xs = GeoInterface.x.(GeoInterface.getpoint(ring)) |> collect
+                ys = GeoInterface.y.(GeoInterface.getpoint(ring)) |> collect
+                subgeom = unsafe_createlinearring(xs, ys)
+                result = GDAL.ogr_g_addgeometrydirectly(poly, subgeom)
+                @ogrerr result "Failed to add linearring."
+            end
+            return poly
         end
-        return poly
     end
-end
 
 
 # Macro to run a block of `code` for multiple modules,
@@ -93,16 +98,7 @@ function _test_implementations_inner(modules::Union{Expr,Vector}, code::Expr)
     for mod in modules1
         expr = Expr(:block)
         for (var, genkey) in pairs(vars)
-            push!(expr.args, :($genkey = if $var isa $(GeoInterface.Extents.Extent)
-                    if $mod in ($GeoInterface, $LibGEOS)
-                        $var
-                    else
-                        $GeoInterface.convert($mod, $(GO.extent_to_polygon)($var))
-                    end
-                else
-                    $GeoInterface.convert($mod, $var)
-                end
-            ))
+            push!(expr.args, conversion_expr(mod, var, genkey))
         end
         push!(expr.args, :(@test $code1))
         push!(tests.args, expr)
@@ -136,7 +132,7 @@ function _testset_implementations_inner(title, modules::Union{Expr,Vector}, code
     for mod in modules1
         expr = Expr(:block)
         for (var, genkey) in pairs(vars)
-            push!(expr.args, :($genkey = $GeoInterface.convert($mod, $var)))
+            push!(expr.args, conversion_expr(mod, var, genkey))
         end
         # Manually define the testset macrocall and all string interpolation
         testset = Expr(
@@ -173,4 +169,4 @@ function _quasiquote!(ex::Expr, vars::Dict)
     return ex
 end
 
-end
+end # module

From 1f924daee57815fedab1b211400bd3f27b702688 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sat, 26 Jul 2025 01:05:19 -0400
Subject: [PATCH 19/24] Comment out shapefile tests since it's broken somehow

---
 test/primitives.jl | 24 ++++++++++++------------
 1 file changed, 12 insertions(+), 12 deletions(-)

diff --git a/test/primitives.jl b/test/primitives.jl
index 77d036cc85..70bbdf74da 100644
--- a/test/primitives.jl
+++ b/test/primitives.jl
@@ -43,18 +43,18 @@ poly = GI.Polygon([lr1, lr2])
                 @warn "Failed to download shapefiles" exception=(e, catch_backtrace())
             else
                 countries_table = Shapefile.Table("countries.shp")
-
-                @testset "Shapefile" begin
-                    centroid_table = GO.apply(GO.centroid, GO.TraitTarget(GI.PolygonTrait(), GI.MultiPolygonTrait()), countries_table);
-                    centroid_geometry = getproperty(centroid_table, :geometry)
-                    # Test that the centroids are correct
-                    @test all(centroid_geometry .== GO.centroid.(countries_table.geometry))
-                    @testset "Columns are preserved" begin  
-                        for column in Iterators.filter(!=(:geometry), GO.Tables.columnnames(countries_table))
-                            @test all(missing_or_equal.(GO.Tables.getcolumn(centroid_table, column), GO.Tables.getcolumn(countries_table, column)))
-                        end
-                    end
-                end
+                # TODO: broken, nfeature not defined in shapefile.jl
+                # @testset "Shapefile" begin
+                #     centroid_table = GO.apply(GO.centroid, GO.TraitTarget(GI.PolygonTrait(), GI.MultiPolygonTrait()), countries_table);
+                #     centroid_geometry = getproperty(centroid_table, :geometry)
+                #     # Test that the centroids are correct
+                #     @test all(centroid_geometry .== GO.centroid.(countries_table.geometry))
+                #     @testset "Columns are preserved" begin  
+                #         for column in Iterators.filter(!=(:geometry), GO.Tables.columnnames(countries_table))
+                #             @test all(missing_or_equal.(GO.Tables.getcolumn(centroid_table, column), GO.Tables.getcolumn(countries_table, column)))
+                #         end
+                #     end
+                # end
 
             @testset "DataFrames" begin
                 countries_df = DataFrames.DataFrame(countries_table)

From df8094fbe90d6ef8a9a8dfec085e52e794b0f103 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sat, 26 Jul 2025 01:06:02 -0400
Subject: [PATCH 20/24] Handle empty geoms in archgdal and gb by no op

Once we fix those bugs we can add them back
---
 test/methods/area.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/methods/area.jl b/test/methods/area.jl
index 3b2c7ea016..4c0db752a6 100644
--- a/test/methods/area.jl
+++ b/test/methods/area.jl
@@ -40,7 +40,7 @@ c = LG.GeometryCollection([p1, p2, r1])
 c_with_epty_l = LG.GeometryCollection([p1, p2, r1, empty_l])
 empty_c = LG.readgeom("GEOMETRYCOLLECTION EMPTY")
 
-@testset_implementations "That handle empty geoms" [LG] begin 
+@testset_implementations "That handle empty geoms" begin 
     @test GO.area($empty_pt) == LG.area($empty_pt) == 0
     @test GO.area($empty_mpt) == LG.area($empty_mpt) == 0
     @test GO.area($empty_l) == LG.area($empty_l) == 0

From d9d59f462d4c5573f19b1f375774087bf2621b39 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sat, 26 Jul 2025 01:12:55 -0400
Subject: [PATCH 21/24] Clean up comment

---
 test/helpers.jl | 24 +++++++++++++++++++-----
 1 file changed, 19 insertions(+), 5 deletions(-)

diff --git a/test/helpers.jl b/test/helpers.jl
index 32cf3870d8..68425ef25a 100644
--- a/test/helpers.jl
+++ b/test/helpers.jl
@@ -16,9 +16,8 @@ function conversion_expr(mod, var, genkey)
             else
                 $GeoInterface.convert($mod, $(GO.extent_to_polygon)($var))
             end
-        # GeometryBasics does not have a Line geometry type.
-        # elseif $mod in ($GeometryBasics,) && $GeoInterface.trait($var) isa $GeoInterface.LineTrait
-        #     $var
+        # These modules do not support empty geometries.
+        # GDAL does but AG does not
         elseif $mod in ($GeoInterface, $ArchGDAL, $GeometryBasics) && $GeoInterface.isempty($var)
             $var
         else
@@ -30,7 +29,7 @@ end
 # TODO: push this up to GB!
 
     # TODO: remove when GB GI pr lands
-    @static if hasmethod(GeometryBasics.convert, (Type{GeometryBasics.LineString}, GeoInterface.LinearRingTrait, Any))
+    # @static if hasmethod(GeometryBasics.convert, (Type{GeometryBasics.LineString}, GeoInterface.LinearRingTrait, Any))
         function GeoInterface.convert(
             ::Type{GeometryBasics.LineString}, 
             ::GeoInterface.LinearRingTrait, 
@@ -39,7 +38,22 @@ end
             return GeoInterface.convert(GeometryBasics.LineString, GeoInterface.LineStringTrait(), geom)
         end
         GeometryBasics.geointerface_geomtype(::GeoInterface.LinearRingTrait) = GeometryBasics.LineString
-    end
+    # end
+
+    # @static if hasmethod(GeometryBasics.convert, (Type{GeometryBasics.Line}, GeoInterface.LineTrait, Any))
+        function GeoInterface.convert(::Type{GeometryBasics.Line}, type::GeoInterface.LineTrait, geom)
+            g1, g2 = GeoInterface.getgeom(geom)
+            x, y = GeoInterface.x(g1), GeoInterface.y(g1)
+            if GeoInterface.is3d(geom)
+                z = GeoInterface.z(g1)
+                T = promote_type(typeof(x), typeof(y), typeof(z))
+                return GeometryBasics.Line{3,T}(GeometryBasics.Point{3,T}(x, y, z), GeometryBasics.Point{3,T}(GeoInterface.x(g2), GeoInterface.y(g2), GeoInterface.z(g2)))
+            else
+                T = promote_type(typeof(x), typeof(y))
+                return GeometryBasics.Line{2,T}(GeometryBasics.Point{2,T}(x, y), GeometryBasics.Point{2,T}(GeoInterface.x(g2), GeoInterface.y(g2)))
+            end
+        end
+    # end
     # end todo
     # GeometryCollection interface - currently just a large Union
     const _ALL_GB_GEOM_TYPES = Union{GeometryBasics.Point, GeometryBasics.LineString, GeometryBasics.Polygon, GeometryBasics.MultiPolygon, GeometryBasics.MultiLineString, GeometryBasics.MultiPoint}

From d81140cfd74bbf6df1a9b70482b4a54749958bc9 Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Sat, 26 Jul 2025 01:28:57 -0400
Subject: [PATCH 22/24] Multiple dispatchify all the intersection accelerators

---
 src/methods/clipping/clipping_processor.jl | 360 +++++++++++----------
 1 file changed, 191 insertions(+), 169 deletions(-)

diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
index 6b09fe2d66..9549d47481 100644
--- a/src/methods/clipping/clipping_processor.jl
+++ b/src/methods/clipping/clipping_processor.jl
@@ -212,221 +212,243 @@ 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
+    manifold::M, accelerator::AutoAccelerator, 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()
-        elseif na < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS || nb < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS
-            SingleNaturalTree()
-        else
-            DoubleNaturalTree()
-        end
+    if na < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS && nb < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS
+        return foreach_pair_of_maybe_intersecting_edges_in_order(manifold, NestedLoop(), f_on_each_a, f_after_each_a, f_on_each_maybe_intersect, poly_a, poly_b, T)
+    elseif na < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS || nb < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS
+        return foreach_pair_of_maybe_intersecting_edges_in_order(manifold, SingleNaturalTree(), f_on_each_a, f_after_each_a, f_on_each_maybe_intersect, poly_a, poly_b, T)
     else
-        accelerator
+        return foreach_pair_of_maybe_intersecting_edges_in_order(manifold, DoubleNaturalTree(), f_on_each_a, f_after_each_a, f_on_each_maybe_intersect, poly_a, poly_b, T)
     end
-    
-    # Declare this now because it's used in multiple branches
-    last_a_idx = 0
+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)
+function foreach_pair_of_maybe_intersecting_edges_in_order(
+    manifold::M, accelerator::NestedLoop, 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}
+    # 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)
+    # 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
-        # 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
+        isnothing(f_after_each_a) || f_after_each_a(a1t, i)
+    end
+    # And we're done!  This is the super simple implementation.
+    return nothing
+end
+
+function foreach_pair_of_maybe_intersecting_edges_in_order(
+    manifold::M, accelerator::SingleSTRtree, 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}
+    na = GI.npoint(poly_a)
+    nb = GI.npoint(poly_b)
+    # 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)
-                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
-                    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
+    # 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)
+            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
+                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
-            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.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 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
+        isnothing(f_after_each_a) || f_after_each_a(a1t, i)
+    end
+    return nothing
+end
+
+function foreach_pair_of_maybe_intersecting_edges_in_order(
+    manifold::M, accelerator::SingleNaturalTree, 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}
+    na = GI.npoint(poly_a)
+    nb = GI.npoint(poly_b)
+    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.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 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
-            isnothing(f_after_each_a) || f_after_each_a(a1t, i)
         end
+        isnothing(f_after_each_a) || f_after_each_a(a1t, i)
+    end
+    return nothing
+end
 
-    elseif accelerator isa DoubleNaturalTree
-        edges_a = to_edgelist(poly_a, T)
-        edges_b = to_edgelist(poly_b, T)
+function foreach_pair_of_maybe_intersecting_edges_in_order(
+    manifold::M, accelerator::DoubleNaturalTree, 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}
+    na = GI.npoint(poly_a)
+    nb = GI.npoint(poly_b)
+    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)
+    tree_a = NaturalIndexing.NaturalIndex(edges_a)
+    tree_b = NaturalIndexing.NaturalIndex(edges_b)
 
-        last_a_idx = 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
-            b1t, b2t = edges_b[b_edge_idx].geom
+    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
 
-            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
+        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
-                !isnothing(f_on_each_a) && f_on_each_a(a1t, a_edge_idx)
             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))
+        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
+        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
 
-        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
+    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
-        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
+    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
-    elseif accelerator isa ThinnedDoubleNaturalTree
-        ext_a, ext_b = GI.extent(poly_a), GI.extent(poly_b)
-        mutual_extent = Extents.intersection(ext_a, ext_b)
+    end
+    return nothing
+end
+    
+function foreach_pair_of_maybe_intersecting_edges_in_order(
+    manifold::M, accelerator::ThinnedDoubleNaturalTree, 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}
+    na = GI.npoint(poly_a)
+    nb = GI.npoint(poly_b)
+    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)
+    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)
+    tree_a = NaturalIndexing.NaturalIndex(edges_a)
+    tree_b = NaturalIndexing.NaturalIndex(edges_b)
 
-        last_a_idx::Int = 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]
-            b_edge_idx = indices_b[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]
 
-            a1t, a2t = edges_a[a_thinned_idx].geom
-            b1t, b2t = edges_b[b_thinned_idx].geom
+        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
+        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
-                !isnothing(f_on_each_a) && f_on_each_a(a1t, a_edge_idx)
             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))
+        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
+        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
-    else
-        error("Unsupported accelerator type: $accelerator.  FosterHormannClipping only supports NestedLoop() or SingleSTRtree().")
     end
-
     return nothing
-
 end
 
 #=

From 517ab8e4e0328c5a8ccc9f17952071e9529222ab Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Fri, 12 Sep 2025 07:09:29 +0200
Subject: [PATCH 23/24] Force autoaccelerator to use nested loops (for now)

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

diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
index 9549d47481..02954c3c1d 100644
--- a/src/methods/clipping/clipping_processor.jl
+++ b/src/methods/clipping/clipping_processor.jl
@@ -219,7 +219,9 @@ function foreach_pair_of_maybe_intersecting_edges_in_order(
     # or -- even now -- just buffering
     na = GI.npoint(poly_a)
     nb = GI.npoint(poly_b)
-
+    # Switching behaviour is turned off in the patch release
+    # This should be turned on in a GO v0.2.x
+    #=
     if na < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS && nb < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS
         return foreach_pair_of_maybe_intersecting_edges_in_order(manifold, NestedLoop(), f_on_each_a, f_after_each_a, f_on_each_maybe_intersect, poly_a, poly_b, T)
     elseif na < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS || nb < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS
@@ -227,6 +229,8 @@ function foreach_pair_of_maybe_intersecting_edges_in_order(
     else
         return foreach_pair_of_maybe_intersecting_edges_in_order(manifold, DoubleNaturalTree(), f_on_each_a, f_after_each_a, f_on_each_maybe_intersect, poly_a, poly_b, T)
     end
+    =#
+    return foreach_pair_of_maybe_intersecting_edges_in_order(manifold, NestedLoop(), f_on_each_a, f_after_each_a, f_on_each_maybe_intersect, poly_a, poly_b, T)
 end
 
 function foreach_pair_of_maybe_intersecting_edges_in_order(

From 146dd543f5f4d124ba2a7b35abed71825384188f Mon Sep 17 00:00:00 2001
From: Anshul Singhvi <anshulsinghvi@gmail.com>
Date: Mon, 15 Sep 2025 14:51:43 +0200
Subject: [PATCH 24/24] Make the constructor default to NestedLoop instead.

---
 src/methods/clipping/clipping_processor.jl | 9 +++------
 1 file changed, 3 insertions(+), 6 deletions(-)

diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
index 02954c3c1d..3b6ca2ccdb 100644
--- a/src/methods/clipping/clipping_processor.jl
+++ b/src/methods/clipping/clipping_processor.jl
@@ -54,9 +54,9 @@ 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)
+FosterHormannClipping(; manifold::Manifold = Planar(), accelerator = nothing) = FosterHormannClipping(manifold, isnothing(accelerator) ? NestedLoop() : accelerator)
+FosterHormannClipping(manifold::Manifold, accelerator::Union{Nothing, IntersectionAccelerator} = nothing) = FosterHormannClipping(manifold, isnothing(accelerator) ? NestedLoop() : accelerator)
+FosterHormannClipping(accelerator::Union{Nothing, IntersectionAccelerator}) = FosterHormannClipping(Planar(), isnothing(accelerator) ? NestedLoop() : 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))
@@ -221,7 +221,6 @@ function foreach_pair_of_maybe_intersecting_edges_in_order(
     nb = GI.npoint(poly_b)
     # Switching behaviour is turned off in the patch release
     # This should be turned on in a GO v0.2.x
-    #=
     if na < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS && nb < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS
         return foreach_pair_of_maybe_intersecting_edges_in_order(manifold, NestedLoop(), f_on_each_a, f_after_each_a, f_on_each_maybe_intersect, poly_a, poly_b, T)
     elseif na < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS || nb < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS
@@ -229,8 +228,6 @@ function foreach_pair_of_maybe_intersecting_edges_in_order(
     else
         return foreach_pair_of_maybe_intersecting_edges_in_order(manifold, DoubleNaturalTree(), f_on_each_a, f_after_each_a, f_on_each_maybe_intersect, poly_a, poly_b, T)
     end
-    =#
-    return foreach_pair_of_maybe_intersecting_edges_in_order(manifold, NestedLoop(), f_on_each_a, f_after_each_a, f_on_each_maybe_intersect, poly_a, poly_b, T)
 end
 
 function foreach_pair_of_maybe_intersecting_edges_in_order(