diff --git a/Project.toml b/Project.toml
index 64c916de6..ec62f505a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -48,7 +48,7 @@ GeometryOpsCore = "=0.1.5"
 LibGEOS = "0.9.2"
 LinearAlgebra = "1"
 Proj = "1"
-SortTileRecursiveTree = "0.1"
+SortTileRecursiveTree = "0.1.2"
 Statistics = "1"
 TGGeometry = "0.1"
 Tables = "1"
diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl
index db7efd2ee..d8d4ae9cf 100644
--- a/src/GeometryOps.jl
+++ b/src/GeometryOps.jl
@@ -21,12 +21,15 @@ using GeoInterface
 using GeometryBasics
 using LinearAlgebra, Statistics
 
+using GeometryBasics.StaticArrays
+
 import Tables, DataAPI
-import GeometryBasics.StaticArrays
 import DelaunayTriangulation # for convex hull and triangulation
 import ExactPredicates
 import Base.@kwdef
 import GeoInterface.Extents: Extents
+import SortTileRecursiveTree
+import SortTileRecursiveTree: STRtree
 
 const GI = GeoInterface
 const GB = GeometryBasics
@@ -44,7 +47,16 @@ include("utils/SpatialTreeInterface/SpatialTreeInterface.jl")
 
 using .LoopStateMachine, .SpatialTreeInterface
 
+include("utils/NaturalIndexing.jl")
+using .NaturalIndexing
+
+# Preparations and prepared geometry
+include("preparations/preparations.jl")
+include("preparations/prepared_geometry.jl")
+include("preparations/monotone_chain.jl")
 
+# Methods - things that don't change the contents
+# of the geometry
 include("methods/angles.jl")
 include("methods/area.jl")
 include("methods/barycentric.jl")
@@ -73,6 +85,8 @@ include("methods/geom_relations/within.jl")
 include("methods/orientation.jl")
 include("methods/polygonize.jl")
 
+# Transformations - things that return a copy of the input geometry,
+# but transformed in some way.
 include("transformations/extent.jl")
 include("transformations/flip.jl")
 include("transformations/reproject.jl")
diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
index 94bf4e859..41b30e1de 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.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
+            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.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
+                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.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
+
+            if last_a_idx < a_edge_idx
+                if !isnothing(f_on_each_a)
+                    for i in last_a_idx:(a_edge_idx-1)
+                        f_on_each_a(a1t, a_edge_idx)
+                        !isnothing(f_after_each_a) && f_after_each_a(a1t, a_edge_idx)
+                    end
+                end
+                !isnothing(f_on_each_a) && f_on_each_a(a1t, a_edge_idx)
+            end
+
+            f_on_each_maybe_intersect(((a1t, a2t), a_edge_idx), ((b1t, b2t), b_edge_idx))
+
+            if last_a_idx < a_edge_idx
+                if !isnothing(f_after_each_a)
+                    f_after_each_a(a1t, a_edge_idx)
+                end
+                last_a_idx = a_edge_idx
+            end
+        end
+    else
+        error("Unsupported accelerator type: $accelerator.  FosterHormannClipping only supports NestedLoop() or SingleSTRtree().")
+    end
+
+    return nothing
+
+end
+
 #=
     _build_a_list(::Type{T}, poly_a, poly_b) -> (a_list, a_idx_list)
 
@@ -176,89 +434,101 @@ function _build_a_list(alg::FosterHormannClipping{M, A}, ::Type{T}, poly_a, poly
     a_list = PolyNode{T}[]  # list of points in poly_a
     sizehint!(a_list, n_a_edges)
     a_idx_list = Vector{Int}()  # finds indices of intersection points in a_list
-    a_count = 0  # number of points added to a_list
-    n_b_intrs = 0
-    # Loop through points of poly_a
-    local a_pt1
-    for (i, a_p2) in enumerate(GI.getpoint(poly_a))
-        a_pt2 = (T(GI.x(a_p2)), T(GI.y(a_p2)))
-        if i <= 1 || (a_pt1 == a_pt2)  # don't repeat points
-            a_pt1 = a_pt2
-            continue
-        end
-        # Add the first point of the edge to the list of points in a_list
-        new_point = PolyNode{T}(;point = a_pt1)
+    local a_count::Int = 0  # number of points added to a_list
+    local n_b_intrs::Int = 0
+    local prev_counter::Int = 0
+
+    function on_each_a(a_pt, i)
+        new_point = PolyNode{T}(;point = a_pt)
         a_count += 1
         push!(a_list, new_point)
-        # Find intersections with edges of poly_b
-        local b_pt1
         prev_counter = a_count
-        for (j, b_p2) in enumerate(GI.getpoint(poly_b))
-            b_pt2 = _tuple_point(b_p2, T)
-            if j <= 1 || (b_pt1 == b_pt2)  # don't repeat points
-                b_pt1 = b_pt2
-                continue
-            end
-            # Determine if edges intersect and how they intersect
-            line_orient, intr1, intr2 = _intersection_point(T, (a_pt1, a_pt2), (b_pt1, b_pt2); exact)
-            if line_orient != line_out  # edges intersect
-                if line_orient == line_cross  # Intersection point that isn't a vertex
-                    int_pt, fracs = intr1
+        return nothing
+    end
+
+    function after_each_a(a_pt, i)
+        # Order intersection points by placement along edge using fracs value
+        if prev_counter < a_count
+            Δintrs = a_count - prev_counter
+            inter_points = @view a_list[(a_count - Δintrs + 1):a_count]
+            sort!(inter_points, by = x -> x.fracs[1])
+        end
+        return nothing
+    end
+
+    function on_each_maybe_intersect(((a_pt1, a_pt2), i), ((b_pt1, b_pt2), j))
+        if (b_pt1 == b_pt2)  # don't repeat points
+            b_pt1 = b_pt2
+            return
+        end
+        # Determine if edges intersect and how they intersect
+        line_orient, intr1, intr2 = _intersection_point(alg.manifold, T, (a_pt1, a_pt2), (b_pt1, b_pt2); exact)
+        if line_orient != line_out  # edges intersect
+            if line_orient == line_cross  # Intersection point that isn't a vertex
+                int_pt, fracs = intr1
+                new_intr = PolyNode{T}(;
+                    point = int_pt, inter = true, neighbor = j, # j is now equivalent to old j-1
+                    crossing = true, fracs = fracs,
+                )
+                a_count += 1
+                n_b_intrs += 1
+                push!(a_list, new_intr)
+                push!(a_idx_list, a_count)
+            else
+                (_, (α1, β1)) = intr1
+                # Determine if a1 or b1 should be added to a_list
+                add_a1 = α1 == 0 && 0 ≤ β1 < 1
+                a1_β = add_a1 ? β1 : zero(T)
+                add_b1 = β1 == 0 && 0 < α1 < 1
+                b1_α = add_b1 ? α1 : zero(T)
+                # If lines are collinear and overlapping, a second intersection exists
+                if line_orient == line_over
+                    (_, (α2, β2)) = intr2
+                    if α2 == 0 && 0 ≤ β2 < 1
+                        add_a1, a1_β = true, β2
+                    end
+                    if β2 == 0 && 0 < α2 < 1
+                        add_b1, b1_α = true, α2
+                    end
+                end
+                # Add intersection points determined above
+                if add_a1
+                    n_b_intrs += a1_β == 0 ? 0 : 1
+                    a_list[prev_counter] = PolyNode{T}(;
+                        point = a_pt1, inter = true, neighbor = j,
+                        fracs = (zero(T), a1_β),
+                    )
+                    push!(a_idx_list, prev_counter)
+                end
+                if add_b1
                     new_intr = PolyNode{T}(;
-                        point = int_pt, inter = true, neighbor = j - 1,
-                        crossing = true, fracs = fracs,
+                        point = b_pt1, inter = true, neighbor = j,
+                        fracs = (b1_α, zero(T)),
                     )
                     a_count += 1
-                    n_b_intrs += 1
                     push!(a_list, new_intr)
                     push!(a_idx_list, a_count)
-                else
-                    (_, (α1, β1)) = intr1
-                    # Determine if a1 or b1 should be added to a_list
-                    add_a1 = α1 == 0 && 0 ≤ β1 < 1
-                    a1_β = add_a1 ? β1 : zero(T)
-                    add_b1 = β1 == 0 && 0 < α1 < 1
-                    b1_α = add_b1 ? α1 : zero(T)
-                    # If lines are collinear and overlapping, a second intersection exists
-                    if line_orient == line_over
-                        (_, (α2, β2)) = intr2
-                        if α2 == 0 && 0 ≤ β2 < 1
-                            add_a1, a1_β = true, β2
-                        end
-                        if β2 == 0 && 0 < α2 < 1
-                            add_b1, b1_α = true, α2
-                        end
-                    end
-                    # Add intersection points determined above
-                    if add_a1
-                        n_b_intrs += a1_β == 0 ? 0 : 1
-                        a_list[prev_counter] = PolyNode{T}(;
-                            point = a_pt1, inter = true, neighbor = j - 1,
-                            fracs = (zero(T), a1_β),
-                        )
-                        push!(a_idx_list, prev_counter)
-                    end
-                    if add_b1
-                        new_intr = PolyNode{T}(;
-                            point = b_pt1, inter = true, neighbor = j - 1,
-                            fracs = (b1_α, zero(T)),
-                        )
-                        a_count += 1
-                        push!(a_list, new_intr)
-                        push!(a_idx_list, a_count)
-                    end
                 end
             end
-            b_pt1 = b_pt2
         end
-        # Order intersection points by placement along edge using fracs value
-        if prev_counter < a_count
-            Δintrs = a_count - prev_counter
-            inter_points = @view a_list[(a_count - Δintrs + 1):a_count]
-            sort!(inter_points, by = x -> x.fracs[1])
+        return nothing
+    end
+
+    # do the iteration but in an accelerated way
+    # this is equivalent to (but faster than)
+    #=
+    ```julia
+    for ((a1, a2), i) in eachedge(poly_a)
+        on_each_a(a1, i)
+        for ((b1, b2), j) in eachedge(poly_b)
+            on_each_maybe_intersect(((a1, a2), i), ((b1, b2), j))
         end
-        a_pt1 = a_pt2
+        after_each_a(a1, i)
     end
+    ```
+    =#
+    foreach_pair_of_maybe_intersecting_edges_in_order(alg, on_each_a, after_each_a, on_each_maybe_intersect, poly_a, poly_b, T)
+
     return a_list, a_idx_list, n_b_intrs
 end
 
diff --git a/src/methods/clipping/intersection.jl b/src/methods/clipping/intersection.jl
index 47d82077c..7af70d740 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.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
@@ -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
diff --git a/src/preparations/monotone_chain.jl b/src/preparations/monotone_chain.jl
new file mode 100644
index 000000000..9499e85e7
--- /dev/null
+++ b/src/preparations/monotone_chain.jl
@@ -0,0 +1,10 @@
+#=
+# Monotone chain
+
+A monotone chain is a continuous list of edges whose slopes are _monotonic_, i.e. all oriented towards the same quadrant.
+
+This speeds up polygon set operations and boolean ops tremendously, since it allows us to skip a lot of the expensive `O(n^2)` operations.
+
+## Example
+=#
+
diff --git a/src/preparations/preparations.jl b/src/preparations/preparations.jl
new file mode 100644
index 000000000..44be7092f
--- /dev/null
+++ b/src/preparations/preparations.jl
@@ -0,0 +1,16 @@
+abstract type AbstractPreparation end
+
+abstract type AbstractPreparationTrait end
+
+function preptrait end
+
+const PrepTuple = Tuple{Vararg{<: AbstractPreparation}}
+
+struct SpatialIndex{IndexType} <: AbstractPreparation
+    index::IndexType
+end
+
+struct SpatialEdgeIndex{IndexType} <: AbstractPreparation
+    index::IndexType
+end
+
diff --git a/src/preparations/prepared_geometry.jl b/src/preparations/prepared_geometry.jl
new file mode 100644
index 000000000..45633c82c
--- /dev/null
+++ b/src/preparations/prepared_geometry.jl
@@ -0,0 +1,56 @@
+struct Prepared{Pa,Pr}
+    parent::Pa
+    preparations::Pr
+end
+
+Base.parent(x::Prepared) = x.parent
+@inline getprep(p::Prepared, x::Symbol) = getproperty(p.preparations, x)
+
+GI.trait(p::Prepared) = GI.trait(parent(p))
+GI.geomtrait(p::Prepared) = GI.geomtrait(parent(p))
+
+GI.isgeometry(::Type{<:Prepared{T}}) where {T} = GI.isgeometry(T)
+GI.isfeature(::Type{<:Prepared{T}}) where {T} = GI.isfeature(T)
+GI.isfeaturecollection(::Type{<:Prepared{T}}) where {T} = GI.isfeaturecollection(T)
+
+GI.geometry(x::Prepared) = GI.geometry(parent(x))
+GI.properties(x::Prepared) = GI.properties(parent(x)) 
+
+for f in (:extent, :crs)
+    @eval GI.$f(t::GI.AbstractTrait, x::Prepared) = GI.$f(t, parent(x))
+end
+for f in (:coordnames, :is3d, :ismeasured, :isempty, :coordinates, :getgeom)
+    @eval GI.$f(t::GI.AbstractGeometryTrait, geom::Prepared, args...) = GI.$f(t, parent(geom), args...)
+end
+
+for f in (:x, :y, :z, :m, :coordinates, :getcoord, :ngeom, :getgeom)
+    @eval GI.$f(t::GI.AbstractPointTrait, geom::Prepared, args...) = GI.$f(t, parent(geom), args...)
+end
+for f in (:npoint, :getpoint, :startpoint, :endpoint, :npoint, :issimple, :isclosed, :isring)
+    @eval GI.$f(t::GI.AbstractCurveTrait, geom::Prepared, args...) = GI.$f(t, parent(geom), args...)
+end
+for f in (:nring, :getring, :getexterior, :nhole, :gethole, :npoint, :getpoint, :startpoint, :endpoint)
+    @eval GI.$f(t::GI.AbstractPolygonTrait, geom::Prepared, args...) = GI.$f(t, parent(geom), args...)
+end
+for f in (:npoint, :getpoint, :issimple)
+    @eval GI.$f(t::GI.AbstractMultiPointTrait, geom::Prepared, args...) = GI.$f(t, parent(geom), args...)
+    @eval GI.$f(t::GI.AbstractMultiCurveTrait, geom::Prepared, args...) = GI.$f(t, parent(geom), args...)
+end
+for f in (:nring, :getring, :npoint, :getpoint)
+    @eval GI.$f(t::GI.AbstractMultiPolygonTrait, geom::Prepared, args...) = GI.$f(t, parent(geom), args...)
+end
+
+getpoint(t::GI.AbstractPolyhedralSurfaceTrait, geom::Prepared) = GI.getpoint(t, parent(geom))
+isclosed(t::GI.AbstractMultiCurveTrait, geom::Prepared) = GI.isclosed(t, parent(geom))
+
+for f in (:getfeature, :coordinates)
+    @eval GI.$f(t::GI.AbstractFeatureTrait, geom::Prepared, args...) = $f(t, parent(geom), args...)
+end
+
+# Ambiguity
+for T in (:LineTrait, :TriangleTrait, :PentagonTrait, :HexagonTrait, :RectangleTrait, :QuadTrait)
+    @eval GI.npoint(t::GI.$T, geom::Prepared) = GI.npoint(t, parent(geom))
+end
+for T in (:RectangleTrait, :QuadTrait, :PentagonTrait, :HexagonTrait, :TriangleTrait)
+    @eval GI.nring(t::GI.$T, geom::Prepared) = GI.nring(t, parent(geom))
+end
diff --git a/src/preparations/rtree.jl b/src/preparations/rtree.jl
new file mode 100644
index 000000000..3325b5b4a
--- /dev/null
+++ b/src/preparations/rtree.jl
@@ -0,0 +1,7 @@
+#=
+# RTree/STRtree
+
+An interface for any arbitrary RTree/STRtree.  This should allow reconstruction from SQL like databases.
+
+Applicable to geometrycollections, multi geometries, and feature collections.
+=#
\ No newline at end of file
diff --git a/src/preparations/sorted_edge_list.jl b/src/preparations/sorted_edge_list.jl
new file mode 100644
index 000000000..69ec7a447
--- /dev/null
+++ b/src/preparations/sorted_edge_list.jl
@@ -0,0 +1,7 @@
+#=
+# Sorted edge list
+
+Soted edge lists are essentially the edges of the linestring, linearring, or polygon, sorted by the initial `y` coordinate.
+
+## Example
+=#
\ No newline at end of file
diff --git a/src/utils/NaturalIndexing.jl b/src/utils/NaturalIndexing.jl
new file mode 100644
index 000000000..080da277e
--- /dev/null
+++ b/src/utils/NaturalIndexing.jl
@@ -0,0 +1,245 @@
+module NaturalIndexing
+
+import GeoInterface as GI
+import Extents
+
+using ..SpatialTreeInterface
+
+import ..GeometryOps as GO # TODO: only needed for NaturallyIndexedRing, remove when that is removed.
+
+export NaturalIndex, NaturallyIndexedRing, prepare_naturally
+
+"""
+    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
+
+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(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{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)
+
+    nlevels = _number_of_levels(nodecapacity, ngeoms)
+
+    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, since we are iterating backwards
+        nrects = _number_of_keys(nodecapacity, nlevels - (level_index), ngeoms)
+        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.
+"""
+    NaturalIndexNode{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 NaturalIndexNode{E <: Extents.Extent}
+    parent_index::NaturalIndex{E}
+    level::Int
+    index::Int
+    extent::E
+end
+
+Extents.extent(node::NaturalIndexNode) = 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
+
+SpatialTreeInterface.isspatialtree(::Type{<: NaturalIndex}) = true
+SpatialTreeInterface.isspatialtree(::Type{<: NaturalIndexNode}) = true
+
+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::NaturalIndexNode, i::Int)
+    child_index = (node.index - 1) * node.parent_index.nodecapacity + i
+    return NaturalIndexNode(
+        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::NaturalIndexNode)
+    return (SpatialTreeInterface.getchild(node, i) for i in 1:SpatialTreeInterface.nchild(node))
+end
+
+SpatialTreeInterface.isleaf(node::NaturalIndexNode) = node.level == length(node.parent_index.levels) - 1
+
+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)
+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(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))
+
+"""
+    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}}}}
+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.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(::Type{<: 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
+
+end # module NaturalIndexing
\ No newline at end of file
diff --git a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl
index d306d75b4..0066142f6 100644
--- a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl
+++ b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl
@@ -7,7 +7,7 @@ import GeoInterface as GI
 import AbstractTrees
 
 # public isspatialtree, isleaf, getchild, nchild, child_indices_extents, node_extent
-export query, do_query
+export query
 export FlatNoTree
 
 # The spatial tree interface and its implementations are defined here.
diff --git a/test/methods/clipping/polygon_clipping.jl b/test/methods/clipping/polygon_clipping.jl
index 2fc4575fa..638b627bd 100644
--- a/test/methods/clipping/polygon_clipping.jl
+++ b/test/methods/clipping/polygon_clipping.jl
@@ -166,7 +166,7 @@ test_pairs = [
 const ϵ = 1e-10
 # Compare clipping results from GeometryOps and LibGEOS
 function compare_GO_LG_clipping(GO_f, LG_f, p1, p2)
-    GO_result_list = GO_f(p1, p2; target = GI.PolygonTrait())
+    
     LG_result_geom = LG_f(p1, p2)
     if LG_result_geom isa LG.GeometryCollection
         poly_list = LG.Polygon[]
@@ -175,38 +175,43 @@ function compare_GO_LG_clipping(GO_f, LG_f, p1, p2)
         end
         LG_result_geom = LG.MultiPolygon(poly_list)
     end
-    # Check if nothing is returned
-    if isempty(GO_result_list) && (LG.isEmpty(LG_result_geom) || LG.area(LG_result_geom) == 0)
-        return true
-    end
-    # Check for unnecessary points
-    if sum(GI.npoint, GO_result_list; init = 0.0) > GI.npoint(LG_result_geom)
-        return false
-    end
-    # Make sure last point is repeated
-    for poly in GO_result_list
-        for ring in GI.getring(poly)
-            GI.getpoint(ring, 1) != GI.getpoint(ring, GI.npoint(ring)) && return false
+
+    for _accelerator in (GO.AutoAccelerator(), GO.NestedLoop(), GO.SingleSTRtree(), GO.SingleNaturalTree(), #=GO.DoubleNaturalTree(), =# #=GO.ThinnedDoubleNaturalTree(), =# #=GO.DoubleSTRtree()=#)
+    @testset let accelerator = _accelerator # this is a ContextTestSet that is otherwise invisible but adds context to the testset
+        GO_result_list = GO_f(GO.FosterHormannClipping(accelerator), p1, p2; target = GI.PolygonTrait())
+        # Check if nothing is returned
+        if isempty(GO_result_list) && (LG.isEmpty(LG_result_geom) || LG.area(LG_result_geom) == 0)
+            @test true
+            continue
+        end
+        # Check for unnecessary points
+        @test !(sum(GI.npoint, GO_result_list; init = 0.0) > GI.npoint(LG_result_geom))
+        # Make sure last point is repeated
+        for poly in GO_result_list
+            for ring in GI.getring(poly)
+                @test !(GI.getpoint(ring, 1) != GI.getpoint(ring, GI.npoint(ring)))
+            end
         end
-    end
 
-    # Check if polygons cover the same area
-    local GO_result_geom
-    if length(GO_result_list) == 1
-        GO_result_geom = GO_result_list[1]
-    else
-        GO_result_geom = GI.MultiPolygon(GO_result_list)
-    end
-    diff_1_area = LG.area(LG.difference(GO_result_geom, LG_result_geom))
-    diff_2_area = LG.area(LG.difference(LG_result_geom, GO_result_geom))
-    return diff_1_area ≤ ϵ && diff_2_area ≤ ϵ
+        # Check if polygons cover the same area
+        local GO_result_geom
+        if length(GO_result_list) == 1
+            GO_result_geom = GO_result_list[1]
+        else
+            GO_result_geom = GI.MultiPolygon(GO_result_list)
+        end
+        diff_1_area = LG.area(LG.difference(GO_result_geom, LG_result_geom))
+        diff_2_area = LG.area(LG.difference(LG_result_geom, GO_result_geom))
+        @test diff_1_area ≤ ϵ && diff_2_area ≤ ϵ
+    end # testset
+    end # loop
 end
 
 # Test clipping functions and print error message if tests fail
 function test_clipping(GO_f, LG_f, f_name)
     for (p1, p2, sg1, sg2, sdesc) in test_pairs
         @testset_implementations "$sg1 $f_name $sg2 - $sdesc" begin
-            @test compare_GO_LG_clipping(GO_f, LG_f, $p1, $p2) 
+            compare_GO_LG_clipping(GO_f, LG_f, $p1, $p2) # this executes tests internally
         end
     end
     return
diff --git a/test/utils/SpatialTreeInterface.jl b/test/utils/SpatialTreeInterface.jl
index eb61f553b..3f78b05c3 100644
--- a/test/utils/SpatialTreeInterface.jl
+++ b/test/utils/SpatialTreeInterface.jl
@@ -4,6 +4,7 @@ using GeometryOps.SpatialTreeInterface
 using GeometryOps.SpatialTreeInterface: isspatialtree, isleaf, getchild, nchild, child_indices_extents, node_extent
 using GeometryOps.SpatialTreeInterface: query, depth_first_search, dual_depth_first_search
 using GeometryOps.SpatialTreeInterface: FlatNoTree
+using GeometryOps.NaturalIndexing: NaturalIndex
 using Extents
 using SortTileRecursiveTree: STRtree
 using NaturalEarth
@@ -218,7 +219,15 @@ end
     test_find_point_in_all_countries(STRtree)
 end
 
-
+# Test NaturalIndex implementation
+@testset "STRtree" begin
+    test_basic_interface(NaturalIndex)
+    test_child_indices_extents(NaturalIndex)
+    test_query_functionality(NaturalIndex)
+    test_dual_query_functionality(NaturalIndex)
+    test_geometry_support(NaturalIndex)
+    test_find_point_in_all_countries(NaturalIndex)
+end
 
 # This testset is not used because Polylabel.jl has some issues.