diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
index 556f140a6..94bf4e859 100644
--- a/src/methods/clipping/clipping_processor.jl
+++ b/src/methods/clipping/clipping_processor.jl
@@ -1,6 +1,58 @@
 # # Polygon clipping helpers
 # This file contains the shared helper functions for the polygon clipping functionalities.
 
+# This file specifically defines helpers for the Foster-Hormann clipping algorithm.
+
+
+"""
+    abstract type IntersectionAccelerator
+
+The abstract supertype for all intersection accelerator types.
+
+The idea is that these speed up the edge-edge intersection checking process,
+perhaps at the cost of memory.
+
+The naive case is `NestedLoop`, which is just a nested loop, running in O(n*m) time.
+
+Then we have `SingleSTRtree`, which is a single STRtree, running in O(n*log(m)) time.
+
+Then we have `DoubleSTRtree`, which is a simultaneous double-tree traversal of two STRtrees.
+
+Finally, we have `AutoAccelerator`, which chooses the best
+accelerator based on the size of the input polygons.  This gets materialized in `build_a_list` for now.
+`AutoAccelerator` should also try to respect existing spatial indices, if they exist.
+"""
+abstract type IntersectionAccelerator end
+struct NestedLoop <: IntersectionAccelerator end
+struct SingleSTRtree <: IntersectionAccelerator end
+struct DoubleSTRtree <: IntersectionAccelerator end
+struct SingleNaturalTree <: IntersectionAccelerator end
+struct DoubleNaturalTree <: IntersectionAccelerator end
+struct ThinnedDoubleNaturalTree <: IntersectionAccelerator end
+struct AutoAccelerator <: IntersectionAccelerator end
+
+"""
+    FosterHormannClipping{M <: Manifold, A <: Union{Nothing, Accelerator}} <: GeometryOpsCore.Algorithm{M} 
+
+Applies the Foster-Hormann clipping algorithm.
+
+# Arguments
+- `manifold::M`: The manifold on which the algorithm operates.
+- `accelerator::A`: The accelerator to use for the algorithm.  Can be `nothing` for automatic choice, or a custom accelerator.
+"""
+struct FosterHormannClipping{M <: Manifold, A <: IntersectionAccelerator} <: GeometryOpsCore.Algorithm{M} 
+    manifold::M
+    accelerator::A
+    # TODO: add exact flag
+    # TODO: should exact flag be in the type domain?
+end
+FosterHormannClipping(; manifold::Manifold = Planar(), accelerator = nothing) = FosterHormannClipping(manifold, isnothing(accelerator) ? AutoAccelerator() : accelerator)
+FosterHormannClipping(manifold::Manifold, accelerator::Union{Nothing, IntersectionAccelerator} = nothing) = FosterHormannClipping(manifold, isnothing(accelerator) ? AutoAccelerator() : accelerator)
+FosterHormannClipping(accelerator::Union{Nothing, IntersectionAccelerator}) = FosterHormannClipping(Planar(), isnothing(accelerator) ? AutoAccelerator() : accelerator)
+# special case for spherical / geodesic manifolds
+# since they can't use STRtrees (because those don't work on the sphere)
+FosterHormannClipping(manifold::Union{Spherical, Geodesic}, accelerator::Union{Nothing, IntersectionAccelerator} = nothing) = FosterHormannClipping(manifold, isnothing(accelerator) ? NestedLoop() : (accelerator isa AutoAccelerator ? NestedLoop() : accelerator))
+
 # This enum defines which side of an edge a point is on
 @enum PointEdgeSide left=1 right=2 unknown=3
 
@@ -37,6 +89,45 @@ PolyNode(node::PolyNode{T};
 
 # Checks equality of two PolyNodes by backing point value, fractional value, and intersection status
 equals(pn1::PolyNode, pn2::PolyNode) = pn1.point == pn2.point && pn1.inter == pn2.inter && pn1.fracs == pn2.fracs
+Base.:(==)(pn1::PolyNode, pn2::PolyNode) = equals(pn1, pn2)
+
+# Finally, we define a nice error type for when the clipping tracing algorithm hits every point in a polygon.
+# This stores the polygons, the a_list, and the b_list, and the a_idx_list.
+# allowing the user to understand what happened and why.
+"""
+    TracingError{T1, T2} <: Exception
+
+An error that is thrown when the clipping tracing algorithm fails somehow.
+This is a bug in the algorithm, and should be reported.
+
+The polygons are contained in the exception object, accessible by try-catch or as `err` in the REPL.
+"""
+struct TracingError{T1, T2, T} <: Exception
+    message::String
+    poly_a::T1
+    poly_b::T2
+    a_list::Vector{PolyNode{T}}
+    b_list::Vector{PolyNode{T}}
+    a_idx_list::Vector{Int}
+end
+
+function Base.showerror(io::IO, e::TracingError{T1, T2}) where {T1, T2}
+    print(io, "TracingError: ")
+    println(io, e.message)
+    println(io, "Please open an issue with the polygons contained in this error object.")
+    println(io)
+    if max(GI.npoint(e.poly_a), GI.npoint(e.poly_b)) < 10
+        println(io, "Polygon A:")
+        println(io, GI.coordinates(e.poly_a))
+        println(io)
+        println(io, "Polygon B:")
+        println(io, GI.coordinates(e.poly_b))
+    else
+        println(io, "The polygons are contained in the exception object, accessible by try-catch or as `err` in the REPL.")
+    end
+end
+
+
 
 #=
     _build_ab_list(::Type{T}, poly_a, poly_b, delay_cross_f, delay_bounce_f; exact) ->
@@ -48,20 +139,20 @@ returns are the fully updated vectors of PolyNodes that represent the rings 'pol
 'poly_b', respectively. This function also returns 'a_idx_list', which at its "ith" index
 stores the index in 'a_list' at which the "ith" intersection point lies.
 =#
-function _build_ab_list(::Type{T}, poly_a, poly_b, delay_cross_f::F1, delay_bounce_f::F2; exact) where {T, F1, F2}
+function _build_ab_list(alg::FosterHormannClipping, ::Type{T}, poly_a, poly_b, delay_cross_f::F1, delay_bounce_f::F2; exact) where {T, F1, F2}
     # Make a list for nodes of each polygon
-    a_list, a_idx_list, n_b_intrs = _build_a_list(T, poly_a, poly_b; exact)
-    b_list = _build_b_list(T, a_idx_list, a_list, n_b_intrs, poly_b)
+    a_list, a_idx_list, n_b_intrs = _build_a_list(alg, T, poly_a, poly_b; exact)
+    b_list = _build_b_list(alg, T, a_idx_list, a_list, n_b_intrs, poly_b)
 
     # Flag crossings
-    _classify_crossing!(T, a_list, b_list; exact)
+    _classify_crossing!(alg, T, a_list, b_list; exact)
 
     # Flag the entry and exits
-    _flag_ent_exit!(T, GI.LinearRingTrait(), poly_b, a_list, delay_cross_f, Base.Fix2(delay_bounce_f, true); exact)
-    _flag_ent_exit!(T, GI.LinearRingTrait(), poly_a, b_list, delay_cross_f, Base.Fix2(delay_bounce_f, false); exact)
+    _flag_ent_exit!(alg, T, GI.LinearRingTrait(), poly_b, a_list, delay_cross_f, Base.Fix2(delay_bounce_f, true); exact)
+    _flag_ent_exit!(alg, T, GI.LinearRingTrait(), poly_a, b_list, delay_cross_f, Base.Fix2(delay_bounce_f, false); exact)
 
     # Set node indices and filter a_idx_list to just crossing points
-    _index_crossing_intrs!(a_list, b_list, a_idx_list)
+    _index_crossing_intrs!(alg, a_list, b_list, a_idx_list)
 
     return a_list, b_list, a_idx_list
 end
@@ -80,7 +171,7 @@ not update the entry and exit flags for a_list.
 The a_idx_list is a list of the indices of intersection points in a_list. The value at
 index i of a_idx_list is the location in a_list where the ith intersection point lies.
 =#
-function _build_a_list(::Type{T}, poly_a, poly_b; exact) where T
+function _build_a_list(alg::FosterHormannClipping{M, A}, ::Type{T}, poly_a, poly_b; exact) where {T, M, A}
     n_a_edges = _nedge(poly_a)
     a_list = PolyNode{T}[]  # list of points in poly_a
     sizehint!(a_list, n_a_edges)
@@ -181,7 +272,7 @@ is needed for clipping using the Greiner-Hormann clipping algorithm.
 Note: after calling this function, b_list is not fully updated. The entry/exit flags still
 need to be updated. However, the neighbor value in a_list is now updated.
 =#
-function _build_b_list(::Type{T}, a_idx_list, a_list, n_b_intrs, poly_b) where T
+function _build_b_list(alg::FosterHormannClipping{M, A}, ::Type{T}, a_idx_list, a_list, n_b_intrs, poly_b) where {T, M, A} 
     # Sort intersection points by insertion order in b_list
     sort!(a_idx_list, by = x-> a_list[x].neighbor + a_list[x].fracs[2])
     # Initialize needed values and lists
@@ -241,7 +332,7 @@ chain is crossing and delayed otherwise and all middle points are marked as boun
 Additionally, the start and end points of the chain are marked as endpoints using the
 endpoints field. 
 =#
-function _classify_crossing!(::Type{T}, a_list, b_list; exact) where T
+function _classify_crossing!(alg::FosterHormannClipping{M, A}, ::Type{T}, a_list, b_list; exact) where {T, M, A}
     napts = length(a_list)
     nbpts = length(b_list)
     # start centered on last point
@@ -265,7 +356,7 @@ function _classify_crossing!(::Type{T}, a_list, b_list; exact) where T
             a_next_is_b_prev = a_next.inter && equals(a_next, b_prev)
             a_next_is_b_next = a_next.inter && equals(a_next, b_next)
             # determine which side of a segments the p points are on
-            b_prev_side, b_next_side = _get_sides(b_prev, b_next, a_prev, curr_pt, a_next,
+            b_prev_side, b_next_side = _get_sides(#=TODO: alg.manifold, =#b_prev, b_next, a_prev, curr_pt, a_next,
                 i, j, a_list, b_list; exact)
             # no sides overlap
             if !a_prev_is_b_prev && !a_prev_is_b_next && !a_next_is_b_prev && !a_next_is_b_next
@@ -407,7 +498,7 @@ all points are intersection points, find the first element that either is the en
 or a crossing point that isn't in a chain. Then take the midpoint of this point and the next
 point in the list and perform the in/out check. If none of these points exist, return
 a `next_idx` of `nothing`. =#
-function _pt_off_edge_status(::Type{T}, pt_list, poly, npts; exact) where T
+function _pt_off_edge_status(alg::FosterHormannClipping{M, A}, ::Type{T}, pt_list, poly, npts; exact) where {T, M, A}
     start_idx, is_non_intr_pt = findfirst(_is_not_intr, pt_list), true
     if isnothing(start_idx)
         start_idx, is_non_intr_pt = findfirst(_next_edge_off, pt_list), false
@@ -419,7 +510,7 @@ function _pt_off_edge_status(::Type{T}, pt_list, poly, npts; exact) where T
     else
         (pt_list[start_idx].point .+ pt_list[next_idx].point) ./ 2
     end
-    start_status = !_point_filled_curve_orientation(start_pt, poly; in = true, on = false, out = false, exact)
+    start_status = !_point_filled_curve_orientation(alg.manifold, start_pt, poly; in = true, on = false, out = false, exact)
     return next_idx, start_status
 end
 # Check if a PolyNode is an intersection point
@@ -445,10 +536,10 @@ bounce will be the same.
 
 Used for clipping polygons by other polygons.
 =#
-function _flag_ent_exit!(::Type{T}, ::GI.LinearRingTrait, poly, pt_list, delay_cross_f, delay_bounce_f; exact) where T
+function _flag_ent_exit!(alg::FosterHormannClipping{M, A}, ::Type{T}, ::GI.LinearRingTrait, poly, pt_list, delay_cross_f, delay_bounce_f; exact) where {T, M, A}
     npts = length(pt_list)
     # Find starting index if there is one
-    next_idx, status = _pt_off_edge_status(T, pt_list, poly, npts; exact)
+    next_idx, status = _pt_off_edge_status(alg, T, pt_list, poly, npts; exact)
     isnothing(next_idx) && return
     start_idx = next_idx - 1 
     # Loop over points and mark entry and exit status
@@ -468,7 +559,7 @@ function _flag_ent_exit!(::Type{T}, ::GI.LinearRingTrait, poly, pt_list, delay_c
                 else  # delayed bouncing
                     next_idx = ii < npts ? (ii + 1) : 1
                     next_val = (curr_pt.point .+ pt_list[next_idx].point) ./ 2
-                    pt_in_poly = _point_filled_curve_orientation(next_val, poly; in = true, on = false, out = false, exact)
+                    pt_in_poly = _point_filled_curve_orientation(alg.manifold, next_val, poly; in = true, on = false, out = false, exact)
                     #= start and end crossing status are the same and depend on if adjacent
                     edges of pt_list are within poly =#
                     start_crossing = delay_bounce_f(pt_in_poly)
@@ -496,8 +587,8 @@ returns false. Used for cutting polygons by lines.
 
 Assumes that the first point is outside of the polygon and not on an edge.
 =#
-function _flag_ent_exit!(::GI.LineTrait, poly, pt_list; exact)
-    status = !_point_filled_curve_orientation(pt_list[1].point, poly; in = true, on = false, out = false, exact)
+function _flag_ent_exit!(alg::FosterHormannClipping{M, A}, ::GI.LineTrait, poly, pt_list; exact) where {M, A}
+    status = !_point_filled_curve_orientation(#=TODO: alg.manifold=#pt_list[1].point, poly; in = true, on = false, out = false, exact)
     # Loop over points and mark entry and exit status
     for (ii, curr_pt) in enumerate(pt_list)
         if curr_pt.crossing
@@ -510,7 +601,7 @@ end
 
 #= Filters a_idx_list to just include crossing points and sets the index of all crossing
 points (which element they correspond to within a_idx_list). =#
-function _index_crossing_intrs!(a_list, b_list, a_idx_list)
+function _index_crossing_intrs!(alg::FosterHormannClipping{M, A}, a_list, b_list, a_idx_list) where {M, A}
     filter!(x -> a_list[x].crossing, a_idx_list)
     for (i, a_idx) in enumerate(a_idx_list)
         curr_node = a_list[a_idx]
@@ -539,7 +630,7 @@ A list of GeoInterface polygons is returned from this function.
 Note: `poly_a` and `poly_b` are temporary inputs used for debugging and can be removed
 eventually.
 =#
-function _trace_polynodes(::Type{T}, a_list, b_list, a_idx_list, f_step, poly_a, poly_b) where T
+function _trace_polynodes(alg::FosterHormannClipping{M, A}, ::Type{T}, a_list, b_list, a_idx_list, f_step, poly_a, poly_b) where {T, M, A}
     n_a_pts, n_b_pts = length(a_list), length(b_list)
     total_pts = n_a_pts + n_b_pts
     n_cross_pts = length(a_idx_list)
@@ -569,7 +660,9 @@ function _trace_polynodes(::Type{T}, a_list, b_list, a_idx_list, f_step, poly_a,
             # changed curr_not_intr to curr_not_same_ent_flag
             same_status, prev_status = true, curr.ent_exit
             while same_status
-                @assert visited_pts < total_pts "Clipping tracing hit every point - clipping error. Please open an issue with polygons: $(GI.coordinates(poly_a)) and $(GI.coordinates(poly_b))."
+                if visited_pts >= total_pts
+                    throw(TracingError("Clipping tracing hit every point - clipping error.", poly_a, poly_b, a_list, b_list, a_idx_list))
+                end
                 # Traverse polygon either forwards or backwards
                 idx += step
                 idx = (idx > curr_npoints) ? mod(idx, curr_npoints) : idx
@@ -615,7 +708,7 @@ or they are separate polygons with no intersection (other than an edge or point)
 Return two booleans that represent if a is inside b (potentially with shared edges / points)
 and visa versa if b is inside of a.
 =#
-function _find_non_cross_orientation(a_list, b_list, a_poly, b_poly; exact)
+function _find_non_cross_orientation(m::M, a_list, b_list, a_poly, b_poly; exact) where {M <: Manifold}
     non_intr_a_idx = findfirst(x -> !x.inter, a_list)
     non_intr_b_idx = findfirst(x -> !x.inter, b_list)
     #= Determine if non-intersection point is in or outside of polygon - if there isn't A
@@ -629,6 +722,9 @@ function _find_non_cross_orientation(a_list, b_list, a_poly, b_poly; exact)
     return a_in_b, b_in_a
 end
 
+_find_non_cross_orientation(alg::FosterHormannClipping{M}, a_list, b_list, a_poly, b_poly; exact) where {M <: Manifold} =
+    _find_non_cross_orientation(alg.manifold, a_list, b_list, a_poly, b_poly; exact)
+
 #=
     _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator, remove_poly_idx; exact)
 
@@ -636,7 +732,7 @@ The holes specified by the hole iterator are added to the polygons in the return
 If this creates more polygons, they are added to the end of the list. If this removes
 polygons, they are removed from the list
 =#
-function _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator, remove_poly_idx; exact) where T
+function _add_holes_to_polys!(alg::FosterHormannClipping{M, A}, ::Type{T}, return_polys, hole_iterator, remove_poly_idx; exact) where {T, M, A}
     n_polys = length(return_polys)
     remove_hole_idx = Int[]
     # Remove set of holes from all polygons
@@ -649,9 +745,9 @@ function _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator, remove_pol
                 curr_poly = return_polys[j]
                 remove_poly_idx[j] && continue
                 curr_poly_ext = GI.nhole(curr_poly) > 0 ? GI.Polygon(StaticArrays.SVector(GI.getexterior(curr_poly))) : curr_poly
-                in_ext, on_ext, out_ext = _line_polygon_interactions(curr_hole, curr_poly_ext; exact, closed_line = true)
+                in_ext, on_ext, out_ext = _line_polygon_interactions(#=TODO: alg.manifold=#curr_hole, curr_poly_ext; exact, closed_line = true)
                 if in_ext  # hole is at least partially within the polygon's exterior
-                    new_hole, new_hole_poly, n_new_pieces = _combine_holes!(T, curr_hole, curr_poly, return_polys, remove_hole_idx)
+                    new_hole, new_hole_poly, n_new_pieces = _combine_holes!(alg, T, curr_hole, curr_poly, return_polys, remove_hole_idx)
                     if n_new_pieces > 0
                         append!(remove_poly_idx, falses(n_new_pieces))
                         n_new_per_poly += n_new_pieces
@@ -659,7 +755,7 @@ function _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator, remove_pol
                     if !on_ext && !out_ext  # hole is completely within exterior
                         push!(curr_poly.geom, new_hole)
                     else  # hole is partially within and outside of polygon's exterior
-                        new_polys = difference(curr_poly_ext, new_hole_poly, T; target=GI.PolygonTrait())
+                        new_polys = difference(alg, curr_poly_ext, new_hole_poly, T; target=GI.PolygonTrait())
                         n_new_polys = length(new_polys) - 1
                         # replace original
                         curr_poly.geom[1] = GI.getexterior(new_polys[1])
@@ -671,7 +767,7 @@ function _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator, remove_pol
                         end
                     end
                 # polygon is completely within hole
-                elseif coveredby(curr_poly_ext, GI.Polygon(StaticArrays.SVector(curr_hole)))
+                elseif coveredby(#=TODO: alg.manifold=#curr_poly_ext, GI.Polygon(StaticArrays.SVector(curr_hole)))
                     remove_poly_idx[j] = true
                 end
             end
@@ -696,16 +792,16 @@ are in the "main" polygon or in one of these new pieces and moved accordingly.
 If the holes don't touch or curr_poly has no holes, then new_hole is returned without any
 changes.
 =#
-function _combine_holes!(::Type{T}, new_hole, curr_poly, return_polys, remove_hole_idx) where T
+function _combine_holes!(alg::FosterHormannClipping{M, A}, ::Type{T}, new_hole, curr_poly, return_polys, remove_hole_idx) where {T, M, A}
     n_new_polys = 0
     empty!(remove_hole_idx)
     new_hole_poly = GI.Polygon(StaticArrays.SVector(new_hole))
     # Combine any existing holes in curr_poly with new hole
     for (k, old_hole) in enumerate(GI.gethole(curr_poly))
         old_hole_poly = GI.Polygon(StaticArrays.SVector(old_hole))
-        if intersects(new_hole_poly, old_hole_poly)
+        if intersects(#=TODO: alg.manifold=#new_hole_poly, old_hole_poly)
             # If the holes intersect, combine them into a bigger hole
-            hole_union = union(new_hole_poly, old_hole_poly, T; target = GI.PolygonTrait())[1]
+            hole_union = union(alg, new_hole_poly, old_hole_poly, T; target = GI.PolygonTrait())[1]
             push!(remove_hole_idx, k + 1)
             new_hole = GI.getexterior(hole_union)
             new_hole_poly = GI.Polygon(StaticArrays.SVector(new_hole))
@@ -734,7 +830,7 @@ end
 
 #= Remove collinear edge points, other than the first and last edge vertex, to simplify
 polygon - including both the exterior ring and any holes=#
-function _remove_collinear_points!(polys, remove_idx, poly_a, poly_b)
+function _remove_collinear_points!(alg::FosterHormannClipping{M, A}, polys, remove_idx, poly_a, poly_b) where {M, A}
     for (i, poly) in Iterators.reverse(enumerate(polys))
         for (j, ring) in Iterators.reverse(enumerate(GI.getring(poly)))
             n = length(ring.geom)
@@ -752,6 +848,7 @@ function _remove_collinear_points!(polys, remove_idx, poly_a, poly_b)
                 else
                     p3 = p
                     # check if p2 is approximately on the edge formed by p1 and p3 - remove if so
+                    # TODO: make this manifold aware
                     if Predicates.orient(p1, p2, p3; exact = False()) == 0
                         remove_idx[i - 1] = true
                     end
diff --git a/src/methods/clipping/cut.jl b/src/methods/clipping/cut.jl
index 3393f3efa..b956035e5 100644
--- a/src/methods/clipping/cut.jl
+++ b/src/methods/clipping/cut.jl
@@ -58,21 +58,23 @@ GI.coordinates.(cut_polys)
  [[[5.0, 0.0], [10.0, 0.0], [10.0, 10.0], [5.0, 10.0], [5.0, 0.0]]]
 ```
 """
-cut(geom, line, ::Type{T} = Float64) where {T <: AbstractFloat} =
-    _cut(T, GI.trait(geom), geom, GI.trait(line), line; exact = True())
+cut(geom, line, ::Type{T} = Float64) where {T <: AbstractFloat} = cut(FosterHormannClipping(), geom, line, T)
+cut(m::Manifold, geom, line, ::Type{T} = Float64) where {T <: AbstractFloat} = cut(FosterHormannClipping(m), geom, line, T)
+cut(alg::FosterHormannClipping{M, A}, geom, line, ::Type{T} = Float64) where {T <: AbstractFloat, M, A} =
+    _cut(alg, T, GI.trait(geom), geom, GI.trait(line), line; exact = True())
 
 #= Cut a given polygon by given line. Add polygon holes back into resulting pieces if there
 are any holes. =#
-function _cut(::Type{T}, ::GI.PolygonTrait, poly, ::GI.LineTrait, line; exact) where T
+function _cut(alg::FosterHormannClipping{M, A}, ::Type{T}, ::GI.PolygonTrait, poly, ::GI.LineTrait, line; exact) where {T, M, A}
     ext_poly = GI.getexterior(poly)
-    poly_list, intr_list = _build_a_list(T, ext_poly, line; exact)
+    poly_list, intr_list = _build_a_list(alg, T, ext_poly, line; exact)
     n_intr_pts = length(intr_list)
     # If an impossible number of intersection points, return original polygon
     if n_intr_pts < 2 || isodd(n_intr_pts)
         return [tuples(poly)]
     end
     # Cut polygon by line
-    cut_coords = _cut(T, ext_poly, line, poly_list, intr_list, n_intr_pts; exact)
+    cut_coords = _cut(alg, T, ext_poly, line, poly_list, intr_list, n_intr_pts; exact)
     # Close coords and create polygons
     for c in cut_coords
         push!(c, c[1])
@@ -80,7 +82,7 @@ function _cut(::Type{T}, ::GI.PolygonTrait, poly, ::GI.LineTrait, line; exact) w
     cut_polys = [GI.Polygon([c]) for c in cut_coords]
     # Add original polygon holes back in
     remove_idx = falses(length(cut_polys))
-    _add_holes_to_polys!(T, cut_polys, GI.gethole(poly), remove_idx; exact)
+    _add_holes_to_polys!(alg, T, cut_polys, GI.gethole(poly), remove_idx; exact)
     return cut_polys
 end
 
@@ -97,10 +99,10 @@ end
 of cut geometry in Vector{Vector{Tuple}} format. 
 
 Note: degenerate cases where intersection points are vertices do not work right now. =#
-function _cut(::Type{T}, geom, line, geom_list, intr_list, n_intr_pts; exact) where T
+function _cut(alg::FosterHormannClipping{M, A}, ::Type{T}, geom, line, geom_list, intr_list, n_intr_pts; exact) where {T, M, A}
     # Sort and categorize the intersection points
     sort!(intr_list, by = x -> geom_list[x].fracs[2])
-    _flag_ent_exit!(GI.LineTrait(), line, geom_list; exact)
+    _flag_ent_exit!(alg, GI.LineTrait(), line, geom_list; exact)
     # Add first point to output list
     return_coords = [[geom_list[1].point]]
     cross_backs = [(T(Inf),T(Inf))]
diff --git a/src/methods/clipping/difference.jl b/src/methods/clipping/difference.jl
index ca7a37078..ff725a881 100644
--- a/src/methods/clipping/difference.jl
+++ b/src/methods/clipping/difference.jl
@@ -33,19 +33,23 @@ GI.coordinates.(diff_poly)
 ```
 """
 function difference(
-    geom_a, geom_b, ::Type{T} = Float64; target=nothing, kwargs...,
+    alg::FosterHormannClipping, geom_a, geom_b, ::Type{T} = Float64; target=nothing, kwargs...,
 ) where {T<:AbstractFloat}
     return _difference(
-        TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b;
+        alg, TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b;
         exact = True(), kwargs...,
     )
 end
+# fallback definitions
+difference(geom_a, geom_b, ::Type{T} = Float64; target=nothing, kwargs...) where T = difference(FosterHormannClipping(Planar()), geom_a, geom_b, T; target, kwargs...)
+# if manifold but no algorithm - assume FosterHormannClipping with provided manifold.
+difference(m::Manifold, geom_a, geom_b, ::Type{T} = Float64; target=nothing, kwargs...) where T = difference(FosterHormannClipping(m), geom_a, geom_b, T; target, kwargs...)
 
 #= The 'difference' function returns the difference of two polygons as a list of polygons.
 The algorithm to determine the difference was adapted from "Efficient clipping of efficient
 polygons," by Greiner and Hormann (1998). DOI: https://doi.org/10.1145/274363.274364 =#
 function _difference(
-    ::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    alg::FosterHormannClipping, target::TraitTarget{GI.PolygonTrait}, ::Type{T},
     ::GI.PolygonTrait, poly_a,
     ::GI.PolygonTrait, poly_b;
     exact, kwargs...
@@ -54,11 +58,11 @@ function _difference(
     ext_a = GI.getexterior(poly_a)
     ext_b = GI.getexterior(poly_b)
     # Find the difference of the exterior of the polygons
-    a_list, b_list, a_idx_list = _build_ab_list(T, ext_a, ext_b, _diff_delay_cross_f, _diff_delay_bounce_f; exact)
-    polys = _trace_polynodes(T, a_list, b_list, a_idx_list, _diff_step, poly_a, poly_b)
+    a_list, b_list, a_idx_list = _build_ab_list(alg, T, ext_a, ext_b, _diff_delay_cross_f, _diff_delay_bounce_f; exact)
+    polys = _trace_polynodes(alg, T, a_list, b_list, a_idx_list, _diff_step, poly_a, poly_b)
     # if no crossing points, determine if either poly is inside of the other
     if isempty(polys)
-        a_in_b, b_in_a = _find_non_cross_orientation(a_list, b_list, ext_a, ext_b; exact)
+        a_in_b, b_in_a = _find_non_cross_orientation(alg.manifold, a_list, b_list, ext_a, ext_b; exact)
         # add case for if they polygons are the same (all intersection points!)
         # add a find_first check to find first non-inter poly!
         if b_in_a && !a_in_b  # b in a and can't be the same polygon
@@ -72,19 +76,19 @@ function _difference(
     remove_idx = falses(length(polys))
     # If the original polygons had holes, take that into account.
     if GI.nhole(poly_a) != 0
-        _add_holes_to_polys!(T, polys, GI.gethole(poly_a), remove_idx; exact)
+        _add_holes_to_polys!(alg, T, polys, GI.gethole(poly_a), remove_idx; exact)
     end
     if GI.nhole(poly_b) != 0
         for hole in GI.gethole(poly_b)
             hole_poly = GI.Polygon(StaticArrays.SVector(hole))
-            new_polys = intersection(hole_poly, poly_a, T; target = GI.PolygonTrait)
+            new_polys = intersection(alg, hole_poly, poly_a, T; target = GI.PolygonTrait)
             if length(new_polys) > 0
                 append!(polys, new_polys)
             end
         end
     end
     # Remove unneeded collinear points on same edge
-    _remove_collinear_points!(polys, remove_idx, poly_a, poly_b)
+    _remove_collinear_points!(alg, polys, remove_idx, poly_a, poly_b)
     return polys
 end
 
@@ -110,7 +114,7 @@ _diff_step(x, y) = (x ⊻ y) ? 1 : (-1)
 #= Polygon with multipolygon difference - note that all intersection regions between
 `poly_a` and any of the sub-polygons of `multipoly_b` are removed from `poly_a`. =#
 function _difference(
-    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    alg::FosterHormannClipping, target::TraitTarget{GI.PolygonTrait}, ::Type{T},
     ::GI.PolygonTrait, poly_a,
     ::GI.MultiPolygonTrait, multipoly_b;
     kwargs...,
@@ -118,7 +122,7 @@ function _difference(
     polys = [tuples(poly_a, T)]
     for poly_b in GI.getpolygon(multipoly_b)
         isempty(polys) && break
-        polys = mapreduce(p -> difference(p, poly_b; target), append!, polys)
+        polys = mapreduce(p -> difference(alg, p, poly_b; target), append!, polys)
     end
     return polys
 end
@@ -128,7 +132,7 @@ sub-polygons of `multipoly_a` and `poly_b` will be removed from the correspondin
 sub-polygon. Unless specified with `fix_multipoly = nothing`, `multipolygon_a` will be
 validated using the given (default is `UnionIntersectingPolygons()`) correction. =#
 function _difference(
-    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    alg::FosterHormannClipping, target::TraitTarget{GI.PolygonTrait}, ::Type{T},
     ::GI.MultiPolygonTrait, multipoly_a,
     ::GI.PolygonTrait, poly_b;
     fix_multipoly = UnionIntersectingPolygons(), kwargs...,
@@ -139,7 +143,7 @@ function _difference(
     polys = Vector{_get_poly_type(T)}()
     sizehint!(polys, GI.npolygon(multipoly_a))
     for poly_a in GI.getpolygon(multipoly_a)
-        append!(polys, difference(poly_a, poly_b; target))
+        append!(polys, difference(alg, poly_a, poly_b; target))
     end
     return polys
 end
@@ -150,7 +154,7 @@ corresponding sub-polygon of `multipoly_a`. Unless specified with `fix_multipoly
 `multipolygon_a` will be validated using the given (default is `UnionIntersectingPolygons()`)
 correction. =#
 function _difference(
-    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    alg::FosterHormannClipping, target::TraitTarget{GI.PolygonTrait}, ::Type{T},
     ::GI.MultiPolygonTrait, multipoly_a,
     ::GI.MultiPolygonTrait, multipoly_b;
     fix_multipoly = UnionIntersectingPolygons(), kwargs...,
@@ -165,9 +169,9 @@ function _difference(
         pieces of `multipolygon_a`` are removed, continue to take difference with new shape
         `polys` =#
         polys = if i == 1
-            difference(multipoly_a, poly_b; target, fix_multipoly)
+            difference(alg, multipoly_a, poly_b; target, fix_multipoly)
         else
-            difference(GI.MultiPolygon(polys), poly_b; target, fix_multipoly)
+            difference(alg, GI.MultiPolygon(polys), poly_b; target, fix_multipoly)
         end
         #= One multipoly_a has been completely covered (and thus removed) there is no need to
         continue taking the difference =#
@@ -175,16 +179,29 @@ function _difference(
     end
     return polys
 end
-
+function _difference(
+    alg::FosterHormannClipping, ::TraitTarget{GI.MultiPolygonTrait}, ::Type{T},
+    trait_a::Union{GI.PolygonTrait, GI.MultiPolygonTrait}, polylike_a,
+    trait_b::Union{GI.PolygonTrait, GI.MultiPolygonTrait}, polylike_b;
+    fix_multipoly = UnionIntersectingPolygons(), kwargs...
+) where T
+    polys = _difference(alg, TraitTarget(GI.PolygonTrait()), T, trait_a, polylike_a, trait_b, polylike_b; kwargs...)
+    if isnothing(fix_multipoly)
+        return GI.MultiPolygon(polys)
+    else
+        return fix_multipoly(GI.MultiPolygon(polys))
+    end
+end
 # Many type and target combos aren't implemented
 function _difference(
-    ::TraitTarget{Target}, ::Type{T},
+    alg::GeometryOpsCore.Algorithm, target::TraitTarget{Target}, ::Type{T},
     trait_a::GI.AbstractTrait, geom_a,
     trait_b::GI.AbstractTrait, geom_b,
+    kw...
 ) where {Target, T}
     @assert(
         false,
-        "Difference between $trait_a and $trait_b with target $Target isn't implemented yet.",
+        "Difference between $trait_a and $trait_b with target $Target and algorithm $alg isn't implemented yet.",
     )
     return nothing
 end
diff --git a/src/methods/clipping/intersection.jl b/src/methods/clipping/intersection.jl
index c2c2abe43..47d82077c 100644
--- a/src/methods/clipping/intersection.jl
+++ b/src/methods/clipping/intersection.jl
@@ -40,28 +40,39 @@ GI.coordinates.(inter_points)
 ```
 """
 function intersection(
-    geom_a, geom_b, ::Type{T}=Float64; target=nothing, kwargs...,
+    alg::FosterHormannClipping, geom_a, geom_b, ::Type{T}=Float64; target=nothing, kwargs...
 ) where {T<:AbstractFloat}
     return _intersection(
-        TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b;
+        alg, TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b;
         exact = True(), kwargs...,
     )
 end
+# fallback definitions
+# if no manifold - assume planar (until we have best_manifold)
+function intersection(
+    geom_a, geom_b, ::Type{T}=Float64; target=nothing, kwargs...
+) where {T<:AbstractFloat}
+    return intersection(FosterHormannClipping(Planar()), geom_a, geom_b; target, kwargs...)
+end
+# if manifold but no algorithm - assume FosterHormannClipping with provided manifold.
+function intersection(m::Manifold, geom_a, geom_b, ::Type{T}=Float64; target=nothing, kwargs...) where {T<:AbstractFloat}
+    return intersection(FosterHormannClipping(m), geom_a, geom_b; target, kwargs...)
+end
 
 # Curve-Curve Intersections with target Point
 _intersection(
-    ::TraitTarget{GI.PointTrait}, ::Type{T},
+    alg::FosterHormannClipping, ::TraitTarget{GI.PointTrait}, ::Type{T},
     trait_a::Union{GI.LineTrait, GI.LineStringTrait, GI.LinearRingTrait}, geom_a,
     trait_b::Union{GI.LineTrait, GI.LineStringTrait, GI.LinearRingTrait}, geom_b;
     kwargs...,
-) where T = _intersection_points(T, trait_a, geom_a, trait_b, geom_b)
+) where T = _intersection_points(alg.manifold, alg.accelerator, T, trait_a, geom_a, trait_b, geom_b)
 
 #= Polygon-Polygon Intersections with target Polygon
 The algorithm to determine the intersection was adapted from "Efficient clipping
 of efficient polygons," by Greiner and Hormann (1998).
 DOI: https://doi.org/10.1145/274363.274364 =#
 function _intersection(
-    ::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    alg::FosterHormannClipping, ::TraitTarget{GI.PolygonTrait}, ::Type{T},
     ::GI.PolygonTrait, poly_a,
     ::GI.PolygonTrait, poly_b;
     exact, kwargs...,
@@ -70,10 +81,10 @@ function _intersection(
     ext_a = GI.getexterior(poly_a)
     ext_b = GI.getexterior(poly_b)
     # Then we find the intersection of the exteriors
-    a_list, b_list, a_idx_list = _build_ab_list(T, ext_a, ext_b, _inter_delay_cross_f, _inter_delay_bounce_f; exact)
-    polys = _trace_polynodes(T, a_list, b_list, a_idx_list, _inter_step, poly_a, poly_b)
+    a_list, b_list, a_idx_list = _build_ab_list(alg, T, ext_a, ext_b, _inter_delay_cross_f, _inter_delay_bounce_f; exact)
+    polys = _trace_polynodes(alg, T, a_list, b_list, a_idx_list, _inter_step, poly_a, poly_b)
     if isempty(polys) # no crossing points, determine if either poly is inside the other
-        a_in_b, b_in_a = _find_non_cross_orientation(a_list, b_list, ext_a, ext_b; exact)
+        a_in_b, b_in_a = _find_non_cross_orientation(alg, a_list, b_list, ext_a, ext_b; exact)
         if a_in_b
             push!(polys, GI.Polygon([tuples(ext_a)]))
         elseif b_in_a
@@ -84,13 +95,12 @@ function _intersection(
     # If the original polygons had holes, take that into account.
     if GI.nhole(poly_a) != 0 || GI.nhole(poly_b) != 0
         hole_iterator = Iterators.flatten((GI.gethole(poly_a), GI.gethole(poly_b)))
-        _add_holes_to_polys!(T, polys, hole_iterator, remove_idx; exact)
+        _add_holes_to_polys!(alg, T, polys, hole_iterator, remove_idx; exact)
     end
     # Remove unneeded collinear points on same edge
-    _remove_collinear_points!(polys, remove_idx, poly_a, poly_b)
+    _remove_collinear_points!(alg, polys, remove_idx, poly_a, poly_b)
     return polys
 end
-
 # # Helper functions for Intersections with Greiner and Hormann Polygon Clipping
 
 #= When marking the crossing status of a delayed crossing, the chain start point is bouncing
@@ -112,7 +122,7 @@ _inter_step(x, _) =  x ? 1 : (-1)
 Unless specified with `fix_multipoly = nothing`, `multipolygon_b` will be validated using
 the given (default is `UnionIntersectingPolygons()`) correction. =#
 function _intersection(
-    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    alg::FosterHormannClipping, target::TraitTarget{GI.PolygonTrait}, ::Type{T},
     ::GI.PolygonTrait, poly_a,
     ::GI.MultiPolygonTrait, multipoly_b;
     fix_multipoly = UnionIntersectingPolygons(), kwargs...,
@@ -122,7 +132,7 @@ function _intersection(
     end
     polys = Vector{_get_poly_type(T)}()
     for poly_b in GI.getpolygon(multipoly_b)
-        append!(polys, intersection(poly_a, poly_b; target))
+        append!(polys, intersection(alg, poly_a, poly_b; target))
     end
     return polys
 end
@@ -131,11 +141,11 @@ end
 polygon with the multipolygon and thus simply switches the order of operations and calls the
 above method. =#
 _intersection(
-    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    alg::FosterHormannClipping, target::TraitTarget{GI.PolygonTrait}, ::Type{T},
     ::GI.MultiPolygonTrait, multipoly_a,
     ::GI.PolygonTrait, poly_b;
     kwargs...,
-) where T = intersection(poly_b, multipoly_a; target , kwargs...)
+) where T = intersection(alg, poly_b, multipoly_a; target , kwargs...)
 
 #= Multipolygon with multipolygon intersection - note that all intersection regions between
 any sub-polygons of `multipoly_a` and any of the sub-polygons of `multipoly_b` are counted
@@ -143,7 +153,7 @@ as intersection polygons. Unless specified with `fix_multipoly = nothing`, both
 `multipolygon_a` and `multipolygon_b` will be validated using the given (default is
 `UnionIntersectingPolygons()`) correction. =#
 function _intersection(
-    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    alg::FosterHormannClipping, target::TraitTarget{GI.PolygonTrait}, ::Type{T},
     ::GI.MultiPolygonTrait, multipoly_a,
     ::GI.MultiPolygonTrait, multipoly_b;
     fix_multipoly = UnionIntersectingPolygons(), kwargs...,
@@ -155,21 +165,36 @@ function _intersection(
     end
     polys = Vector{_get_poly_type(T)}()
     for poly_a in GI.getpolygon(multipoly_a)
-        append!(polys, intersection(poly_a, multipoly_b; target, fix_multipoly))
+        append!(polys, intersection(alg, poly_a, multipoly_b; target, fix_multipoly))
     end
     return polys
 end
+# catch-all method for multipolygontraits
+function _intersection(
+    alg::FosterHormannClipping, ::TraitTarget{GI.MultiPolygonTrait}, ::Type{T},
+    trait_a::Union{GI.PolygonTrait, GI.MultiPolygonTrait}, polylike_a,
+    trait_b::Union{GI.PolygonTrait, GI.MultiPolygonTrait}, polylike_b;
+    fix_multipoly = UnionIntersectingPolygons(), kwargs...
+) where T
+    polys = _intersection(alg, TraitTarget(GI.PolygonTrait()), T, trait_a, polylike_a, trait_b, polylike_b; kwargs...)
+    if isnothing(fix_multipoly)
+        return GI.MultiPolygon(polys)
+    else
+        return fix_multipoly(GI.MultiPolygon(polys))
+    end
+end
+
 
 # Many type and target combos aren't implemented
 function _intersection(
-    ::TraitTarget{Target}, ::Type{T},
+    alg::GeometryOpsCore.Algorithm, target::TraitTarget{Target}, ::Type{T},
     trait_a::GI.AbstractTrait, geom_a,
     trait_b::GI.AbstractTrait, geom_b;
     kwargs...,
 ) where {Target, T}
     @assert(
         false,
-        "Intersection between $trait_a and $trait_b with target $Target isn't implemented yet.",
+        "Intersection between $trait_a and $trait_b with target $Target and algorithm $alg isn't implemented yet.",
     )
     return nothing
 end
@@ -193,13 +218,19 @@ inter_points = GO.intersection_points(line1, line2)
 1-element Vector{Tuple{Float64, Float64}}:
  (125.58375366067548, -14.83572303404496)
 """
-intersection_points(geom_a, geom_b, ::Type{T} = Float64) where T <: AbstractFloat =
-    _intersection_points(T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
+intersection_points(geom_a, geom_b, ::Type{T} = Float64) where T <: AbstractFloat = intersection_points(FosterHormannClipping(Planar()), geom_a, geom_b, T)
+function intersection_points(alg::FosterHormannClipping{M, A}, geom_a, geom_b, ::Type{T} = Float64) where {M, A, T <: AbstractFloat}
+    return _intersection_points(alg.manifold, alg.accelerator, T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
+end
+
+function intersection_points(m::Manifold, a::IntersectionAccelerator, geom_a, geom_b, ::Type{T} = Float64) where T <: AbstractFloat
+    return _intersection_points(m, a, T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
+end
 
 
 #= Calculates the list of intersection points between two geometries, including line
 segments, line strings, linear rings, polygons, and multipolygons. =#
-function _intersection_points(::Type{T}, ::GI.AbstractTrait, a, ::GI.AbstractTrait, b; exact = True()) where T
+function _intersection_points(manifold::M, accelerator::A, ::Type{T}, ::GI.AbstractTrait, a, ::GI.AbstractTrait, b; exact = True()) where {M <: Manifold, A <: IntersectionAccelerator, T}
     # Initialize an empty list of points
     result = Tuple{T, T}[]
     # Check if the geometries extents even overlap
@@ -243,7 +274,7 @@ intersection point (x,y) while the second is the ratio along the initial lines (
 that point. 
 
 Calculation derivation can be found here: https://stackoverflow.com/questions/563198/ =#
-function _intersection_point(::Type{T}, (a1, a2)::Edge, (b1, b2)::Edge; exact) where T
+function _intersection_point(manifold::M, ::Type{T}, (a1, a2)::Edge, (b1, b2)::Edge; exact) where {M <: Manifold, T}
     # Default answer for no intersection
     line_orient = line_out
     intr1 = ((zero(T), zero(T)), (zero(T), zero(T)))
@@ -266,7 +297,7 @@ function _intersection_point(::Type{T}, (a1, a2)::Edge, (b1, b2)::Edge; exact) w
     # Determine intersection type and intersection point(s)
     if a1_orient == a2_orient == b1_orient == b2_orient == 0
         # Intersection is collinear if all endpoints lie on the same line
-        line_orient, intr1, intr2 = _find_collinear_intersection(T, a1, a2, b1, b2, a_ext, b_ext, no_intr_result)
+        line_orient, intr1, intr2 = _find_collinear_intersection(manifold, T, a1, a2, b1, b2, a_ext, b_ext, no_intr_result)
     elseif a1_orient == 0 || a2_orient == 0 || b1_orient == 0 || b2_orient == 0
         # Intersection is a hinge if the intersection point is an endpoint
         line_orient = line_hinge
@@ -279,13 +310,16 @@ function _intersection_point(::Type{T}, (a1, a2)::Edge, (b1, b2)::Edge; exact) w
     return line_orient, intr1, intr2
 end
 
+# TODO: deprecate this
+_intersection_point(::Type{T}, (a1, a2)::Edge, (b1, b2)::Edge; exact) where T = _intersection_point(Planar(), T, (a1, a2), (b1, b2); exact)
+
 #= If lines defined by (a1, a2) and (b1, b2) are collinear, find endpoints of overlapping
 region if they exist. This could result in three possibilities. First, there could be no
 overlapping region, in which case, the default 'no_intr_result' intersection information is
 returned. Second, the two regions could just meet at one shared endpoint, in which case it
 is a hinge intersection with one intersection point. Otherwise, it is a overlapping
 intersection defined by two of the endpoints of the line segments. =#
-function _find_collinear_intersection(::Type{T}, a1, a2, b1, b2, a_ext, b_ext, no_intr_result) where T
+function _find_collinear_intersection(manifold::M, ::Type{T}, a1, a2, b1, b2, a_ext, b_ext, no_intr_result) where {M <: Manifold, T}
     # Define default return for no intersection points
     line_orient, intr1, intr2 = no_intr_result
     # Determine collinear line overlaps
diff --git a/src/methods/clipping/union.jl b/src/methods/clipping/union.jl
index 089ec6524..17af626da 100644
--- a/src/methods/clipping/union.jl
+++ b/src/methods/clipping/union.jl
@@ -33,19 +33,32 @@ GI.coordinates.(union_poly)
 ```
 """
 function union(
-    geom_a, geom_b, ::Type{T}=Float64; target=nothing, kwargs...
+    alg::FosterHormannClipping, geom_a, geom_b, ::Type{T}=Float64; target=nothing, kwargs...
 ) where {T<:AbstractFloat}
     return _union(
-        TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b;
+        alg, TraitTarget(target), T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b;
         exact = True(), kwargs...,
     )
 end
 
+# fallback definitions
+# if no manifold - assume planar (until we have best_manifold)
+function union(
+    geom_a, geom_b, ::Type{T}=Float64; target=nothing, kwargs...
+) where {T<:AbstractFloat}
+    return union(FosterHormannClipping(Planar()), geom_a, geom_b; target, kwargs...)
+end
+
+# if manifold but no algorithm - assume FosterHormannClipping with provided manifold.
+function union(m::Manifold, geom_a, geom_b, ::Type{T}=Float64; target=nothing, kwargs...) where {T<:AbstractFloat}
+    return union(FosterHormannClipping(m), geom_a, geom_b; target, kwargs...)
+end
+
 #= This 'union' implementation returns the union of two polygons. The algorithm to determine
 the union was adapted from "Efficient clipping of efficient polygons," by Greiner and
 Hormann (1998). DOI: https://doi.org/10.1145/274363.274364 =#
 function _union(
-    ::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    alg::FosterHormannClipping, ::TraitTarget{GI.PolygonTrait}, ::Type{T},
     ::GI.PolygonTrait, poly_a,
     ::GI.PolygonTrait, poly_b;
     exact, kwargs...,
@@ -54,13 +67,13 @@ function _union(
     ext_a = GI.getexterior(poly_a)
     ext_b = GI.getexterior(poly_b)
     # Then, I get the union of the exteriors
-    a_list, b_list, a_idx_list = _build_ab_list(T, ext_a, ext_b, _union_delay_cross_f, _union_delay_bounce_f; exact)
-    polys = _trace_polynodes(T, a_list, b_list, a_idx_list, _union_step, poly_a, poly_b)
+    a_list, b_list, a_idx_list = _build_ab_list(alg, T, ext_a, ext_b, _union_delay_cross_f, _union_delay_bounce_f; exact)
+    polys = _trace_polynodes(alg, T, a_list, b_list, a_idx_list, _union_step, poly_a, poly_b)
     n_pieces = length(polys)
     # Check if one polygon totally within other and if so, return the larger polygon
     a_in_b, b_in_a = false, false
     if n_pieces == 0 # no crossing points, determine if either poly is inside the other
-        a_in_b, b_in_a = _find_non_cross_orientation(a_list, b_list, ext_a, ext_b; exact)
+        a_in_b, b_in_a = _find_non_cross_orientation(alg, a_list, b_list, ext_a, ext_b; exact)
         if a_in_b
             push!(polys, GI.Polygon([_linearring(tuples(ext_b))]))
         elseif b_in_a
@@ -81,10 +94,10 @@ function _union(
     end
     # Add in holes
     if GI.nhole(poly_a) != 0 || GI.nhole(poly_b) != 0
-        _add_union_holes!(polys, a_in_b, b_in_a, poly_a, poly_b; exact)
+        _add_union_holes!(alg, polys, a_in_b, b_in_a, poly_a, poly_b; exact)
     end
     # Remove unneeded collinear points on same edge
-    _remove_collinear_points!(polys, [false], poly_a, poly_b)
+    _remove_collinear_points!(alg, polys, [false], poly_a, poly_b)
     return polys
 end
 
@@ -108,11 +121,11 @@ _union_step(x, _) = x ? (-1) : 1
 #= Add holes from two polygons to the exterior polygon formed by their union. If adding the
 the holes reveals that the polygons aren't actually intersecting, return the original
 polygons. =#
-function _add_union_holes!(polys, a_in_b, b_in_a, poly_a, poly_b; exact)
+function _add_union_holes!(alg::FosterHormannClipping, polys, a_in_b, b_in_a, poly_a, poly_b; exact)
     if a_in_b
-        _add_union_holes_contained_polys!(polys, poly_a, poly_b; exact)
+        _add_union_holes_contained_polys!(alg, polys, poly_a, poly_b; exact)
     elseif b_in_a
-        _add_union_holes_contained_polys!(polys, poly_b, poly_a; exact)
+        _add_union_holes_contained_polys!(alg, polys, poly_b, poly_a; exact)
     else  # Polygons intersect, but neither is contained in the other
         n_a_holes = GI.nhole(poly_a)
         ext_poly_a = GI.Polygon(StaticArrays.SVector(GI.getexterior(poly_a)))
@@ -125,7 +138,7 @@ function _add_union_holes!(polys, a_in_b, b_in_a, poly_a, poly_b; exact)
         # Loop over all holes in both original polygons
         for (i, ih) in enumerate(Iterators.flatten((GI.gethole(poly_a), GI.gethole(poly_b))))
             ih = _linearring(ih)
-            in_ext, _, _ = _line_polygon_interactions(ih, curr_exterior_poly; exact, closed_line = true)
+            in_ext, _, _ = _line_polygon_interactions(#=TODO: alg.manifold=# ih, curr_exterior_poly; exact, closed_line = true)
             if !in_ext
                 #= if the hole isn't in the overlapping region between the two polygons, add
                 the hole to the resulting polygon as we know it can't interact with any
@@ -138,7 +151,7 @@ function _add_union_holes!(polys, a_in_b, b_in_a, poly_a, poly_b; exact)
                 between poly_a and poly_b within the overlap are added, in addition to all
                 holes in non-overlapping regions =#
                 h_poly = GI.Polygon(StaticArrays.SVector(ih))
-                new_holes = difference(h_poly, current_poly; target = GI.PolygonTrait())
+                new_holes = difference(alg, h_poly, current_poly; target = GI.PolygonTrait())
                 append!(polys[1].geom, (GI.getexterior(new_h) for new_h in new_holes))
             end
             if i == n_a_holes
@@ -153,12 +166,12 @@ end
 #= Add holes holes to the union of two polygons where one of the original polygons was
 inside of the other. If adding the the holes reveal that the polygons aren't actually
 intersecting, return the original polygons.=#
-function _add_union_holes_contained_polys!(polys, interior_poly, exterior_poly; exact)
+function _add_union_holes_contained_polys!(alg::FosterHormannClipping, polys, interior_poly, exterior_poly; exact)
     union_poly = polys[1]
     ext_int_ring = GI.getexterior(interior_poly)
     for (i, ih) in enumerate(GI.gethole(exterior_poly))
         poly_ih = GI.Polygon(StaticArrays.SVector(ih))
-        in_ih, on_ih, out_ih = _line_polygon_interactions(ext_int_ring, poly_ih; exact, closed_line = true)
+        in_ih, on_ih, out_ih = _line_polygon_interactions(#=TODO: alg.manifold=# ext_int_ring, poly_ih; exact, closed_line = true)
         if in_ih  # at least part of interior polygon exterior is within the ith hole
             if !on_ih && !out_ih
                 #= interior polygon is completely within the ith hole - polygons aren't
@@ -169,7 +182,7 @@ function _add_union_holes_contained_polys!(polys, interior_poly, exterior_poly;
             else
                 #= interior polygon is partially within the ith hole - area of interior
                 polygon reduces the size of the hole =#
-                new_holes = difference(poly_ih, interior_poly; target = GI.PolygonTrait())
+                new_holes = difference(alg, poly_ih, interior_poly; target = GI.PolygonTrait())
                 append!(union_poly.geom, (GI.getexterior(new_h) for new_h in new_holes))
             end
         else  # none of interior polygon exterior is within the ith hole
@@ -183,14 +196,14 @@ function _add_union_holes_contained_polys!(polys, interior_poly, exterior_poly;
                 #= interior polygon's exterior is outside of the ith hole - the interior
                 polygon could either be disjoint from the hole, or contain the hole =#
                 ext_int_poly = GI.Polygon(StaticArrays.SVector(ext_int_ring))
-                in_int, _, _ = _line_polygon_interactions(ih, ext_int_poly; exact, closed_line = true)
+                in_int, _, _ = _line_polygon_interactions(#=TODO: alg.manifold=# ih, ext_int_poly; exact, closed_line = true)
                 if in_int
                     #= interior polygon contains the hole - overlapping holes between the
                     interior and exterior polygons will be added =#
                     for jh in GI.gethole(interior_poly)
                         poly_jh = GI.Polygon(StaticArrays.SVector(jh))
                         if intersects(poly_ih, poly_jh)
-                            new_holes = intersection(poly_ih, poly_jh; target = GI.PolygonTrait())
+                            new_holes = intersection(#=TODO: alg, =#poly_ih, poly_jh; target = GI.PolygonTrait())
                             append!(union_poly.geom, (GI.getexterior(new_h) for new_h in new_holes))
                         end
                     end
@@ -210,7 +223,7 @@ included, unioning these sub-polygons with `poly_a` where they intersect. Unless
 with `fix_multipoly = nothing`, `multipolygon_b` will be validated using the given (default
 is `UnionIntersectingPolygons()`) correction. =#
 function _union(
-    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    alg::FosterHormannClipping, target::TraitTarget{GI.PolygonTrait}, ::Type{T},
     ::GI.PolygonTrait, poly_a,
     ::GI.MultiPolygonTrait, multipoly_b;
     fix_multipoly = UnionIntersectingPolygons(), kwargs...,
@@ -220,9 +233,9 @@ function _union(
     end
     polys = [tuples(poly_a, T)]
     for poly_b in GI.getpolygon(multipoly_b)
-        if intersects(polys[1], poly_b)
+        if intersects(#=TODO: alg.manifold, =#polys[1], poly_b)
             # If polygons intersect and form a new polygon, swap out polygon
-            new_polys = union(polys[1], poly_b; target)
+            new_polys = union(alg, polys[1], poly_b; target)
             if length(new_polys) > 1 # case where they intersect by just one point
                 push!(polys, tuples(poly_b, T))  # add poly_b to list
             else
@@ -239,18 +252,18 @@ end
 #= Multipolygon with polygon union is equivalent to taking the union of the polygon with the
 multipolygon and thus simply switches the order of operations and calls the above method. =#
 _union(
-    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    alg::FosterHormannClipping, target::TraitTarget{GI.PolygonTrait}, ::Type{T},
     ::GI.MultiPolygonTrait, multipoly_a,
     ::GI.PolygonTrait, poly_b;
     kwargs...,
-) where T = union(poly_b, multipoly_a; target, kwargs...)
+) where T = union(alg, poly_b, multipoly_a; target, kwargs...)
 
 #= Multipolygon with multipolygon union - note that all of the sub-polygons of `multipoly_a`
 and the sub-polygons of `multipoly_b` are included and combined together where there are
 intersections. Unless specified with `fix_multipoly = nothing`, `multipolygon_b` will be
 validated using the given (default is `UnionIntersectingPolygons()`) correction. =#
 function _union(
-    target::TraitTarget{GI.PolygonTrait}, ::Type{T},
+    alg::FosterHormannClipping, target::TraitTarget{GI.PolygonTrait}, ::Type{T},
     ::GI.MultiPolygonTrait, multipoly_a,
     ::GI.MultiPolygonTrait, multipoly_b;
     fix_multipoly = UnionIntersectingPolygons(), kwargs...,
@@ -262,15 +275,54 @@ function _union(
     multipolys = multipoly_b
     local polys
     for poly_a in GI.getpolygon(multipoly_a)
-        polys = union(poly_a, multipolys; target, fix_multipoly)
+        polys = union(alg, poly_a, multipolys; target, fix_multipoly)
         multipolys = GI.MultiPolygon(polys)
     end
     return polys
 end
-
+function _union(
+    alg::FosterHormannClipping, target::TraitTarget{GI.MultiPolygonTrait}, ::Type{T},
+    ::GI.PolygonTrait, poly_a,
+    ::GI.PolygonTrait, poly_b;
+    kwargs...,
+) where T
+    return UnionIntersectingPolygons(#=TODO: alg.manifold=#)(GI.MultiPolygon([poly_a, poly_b]))
+end
+function _union(
+    alg::FosterHormannClipping, target::TraitTarget{GI.MultiPolygonTrait}, ::Type{T},
+    ::GI.PolygonTrait, poly_a,
+    ::GI.MultiPolygonTrait, multipoly_b;
+    fix_multipoly = UnionIntersectingPolygons(),
+    kwargs...,
+) where T 
+    res = union(alg, multipoly_b, poly_a; target = GI.PolygonTrait(), kwargs...)
+    if !isnothing(fix_multipoly)
+        return fix_multipoly(GI.MultiPolygon(res))
+    else
+        return GI.MultiPolygon(res)
+    end
+end
+# this is the opposite of the above
+function _union(
+    alg::FosterHormannClipping, target::TraitTarget{GI.MultiPolygonTrait}, ::Type{T},
+    ::GI.MultiPolygonTrait, multipoly_a,
+    ::GI.PolygonTrait, poly_b;
+    fix_multipoly = UnionIntersectingPolygons(),
+    kwargs...,
+) where T 
+    union(alg, poly_b, multipoly_a; target, fix_multipoly, kwargs...)
+end
+function _union(
+    alg::FosterHormannClipping, target::TraitTarget{GI.MultiPolygonTrait}, ::Type{T},
+    trait_a::GI.MultiPolygonTrait, multipoly_a,
+    trait_b::GI.MultiPolygonTrait, multipoly_b;
+    fix_multipoly = UnionIntersectingPolygons(), kwargs...,
+) where T
+    return GI.MultiPolygon(_union(alg, TraitTarget{GI.PolygonTrait}(), T, trait_a, multipoly_a, trait_b, multipoly_b; fix_multipoly, kwargs...))
+end
 # Many type and target combos aren't implemented
 function _union(
-    ::TraitTarget{Target}, ::Type{T},
+    alg::GeometryOpsCore.Algorithm, target::TraitTarget{Target}, ::Type{T},
     trait_a::GI.AbstractTrait, geom_a,
     trait_b::GI.AbstractTrait, geom_b;
     kwargs...
diff --git a/src/methods/geom_relations/geom_geom_processors.jl b/src/methods/geom_relations/geom_geom_processors.jl
index 7e4658a7a..5f336f652 100644
--- a/src/methods/geom_relations/geom_geom_processors.jl
+++ b/src/methods/geom_relations/geom_geom_processors.jl
@@ -491,7 +491,7 @@ of the curve if it didn't return 'on'.
 See paper for more information on cases denoted in comments.
 =#
 function _point_filled_curve_orientation(
-    point, curve;
+    ::Planar, point, curve;
     in::T = point_in, on::T = point_on, out::T = point_out, exact,
 ) where {T}
     x, y = GI.x(point), GI.y(point)
@@ -525,6 +525,10 @@ function _point_filled_curve_orientation(
     end
     return iseven(k) ? out : in
 end
+_point_filled_curve_orientation(
+    point, curve;
+    in::T = point_in, on::T = point_on, out::T = point_out, exact,
+) where {T} = _point_filled_curve_orientation(Planar(), point, curve; in, on, out, exact)
 
 #=
 Determines the types of interactions of a line with a filled-in curve. By
diff --git a/test/methods/clipping/polygon_clipping.jl b/test/methods/clipping/polygon_clipping.jl
index 432650e94..2fc4575fa 100644
--- a/test/methods/clipping/polygon_clipping.jl
+++ b/test/methods/clipping/polygon_clipping.jl
@@ -215,3 +215,20 @@ end
 @testset "Intersection" begin test_clipping(GO.intersection, LG.intersection, "intersection") end
 @testset "Union" begin test_clipping(GO.union, LG.union, "union") end
 @testset "Difference" begin test_clipping(GO.difference, LG.difference, "difference") end
+
+@testset "TracingError show" begin
+    message = "Test tracing error"
+    poly_a_small = p1
+    poly_b_small = p2
+
+    poly_a_large, poly_b_large = GO.segmentize.((poly_a_small, poly_b_small), max_distance = 0.001)
+
+    # Test that the message is contained in the exception
+    @test_throws "Test tracing error" throw(GO.TracingError(message, poly_a_small, poly_b_small, GO.PolyNode{Float64}[], GO.PolyNode{Float64}[], Int[]))
+    # Test that the coordinates of the polygons are contained in the exception,
+    # if the polygon is small enough
+    @test_throws "$(GI.coordinates(poly_a_small))" throw(GO.TracingError(message, poly_a_small, poly_b_small, GO.PolyNode{Float64}[], GO.PolyNode{Float64}[], Int[]))
+    # Test that the coordinates of the polygons are not contained in the exception,
+    # if the polygon is large enoughs
+    @test_throws "The polygons are contained in the exception object" throw(GO.TracingError(message, poly_a_large, poly_b_large, GO.PolyNode{Float64}[], GO.PolyNode{Float64}[], Int[]))
+end
\ No newline at end of file