diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl
index 30c1fcf29..e46bc0ebf 100644
--- a/src/GeometryOps.jl
+++ b/src/GeometryOps.jl
@@ -49,6 +49,8 @@ include("utils/UnitSpherical/UnitSpherical.jl")
 
 # Load utility modules in
 using .LoopStateMachine, .SpatialTreeInterface, .UnitSpherical
+const LSM = LoopStateMachine
+const STI = SpatialTreeInterface
 
 include("utils/NaturalIndexing.jl")
 using .NaturalIndexing
diff --git a/src/methods/distance.jl b/src/methods/distance.jl
index 5b16e98c5..d429bbe7d 100644
--- a/src/methods/distance.jl
+++ b/src/methods/distance.jl
@@ -192,7 +192,7 @@ end
 
 # Returns the Euclidean distance between two points.
 Base.@propagate_inbounds _euclid_distance(::Type{T}, p1, p2) where T =
-    sqrt(_squared_euclid_distance(T, p1, p2))
+    hypot(T(GI.x(p2)) - T(GI.x(p1)), T(GI.y(p2)) - T(GI.y(p1)))
 
 # Returns the square of the euclidean distance between two points
 Base.@propagate_inbounds _squared_euclid_distance(::Type{T}, p1, p2) where T =
diff --git a/src/methods/geom_relations/geom_geom_processors.jl b/src/methods/geom_relations/geom_geom_processors.jl
index 0049f993d..b96bbf871 100644
--- a/src/methods/geom_relations/geom_geom_processors.jl
+++ b/src/methods/geom_relations/geom_geom_processors.jl
@@ -492,6 +492,7 @@ passes through an odd number of edges, it is within the curve, else outside of
 of the curve if it didn't return 'on'.
 See paper for more information on cases denoted in comments.
 =#
+
 function _point_filled_curve_orientation(
     ::Planar, point, curve;
     in::T = point_in, on::T = point_on, out::T = point_out, exact,
@@ -527,6 +528,54 @@ function _point_filled_curve_orientation(
     end
     return iseven(k) ? out : in
 end
+
+# Specialized implementation for NaturallyIndexedRing
+# This relies on multidispatch.
+# TODO: remove?
+function _point_filled_curve_orientation(
+    ::Planar, point, curve::NaturalIndexing.NaturallyIndexedRing;
+    in::T = point_in, on::T = point_on, out::T = point_out, exact,
+) where {T}
+    x, y = GI.x(GI.PointTrait(), point), GI.y(GI.PointTrait(), point)
+    k::Int = 0  # counter for ray crossings
+
+    tree = curve.index
+
+    function per_edge_function(i)
+        p_start = _tuple_point(GI.getpoint(curve, i))
+        p_end = _tuple_point(GI.getpoint(curve, i + 1))
+        v1 = GI.y(p_start) - y
+        v2 = GI.y(p_end) - y
+        if !((v1 < 0 && v2 < 0) || (v1 > 0 && v2 > 0)) # if not cases 11 or 26
+            u1, u2 = GI.x(p_start) - x, GI.x(p_end) - x
+            f = Predicates.orient(p_start, p_end, (x, y); exact)
+            if v2 > 0 && v1 ≤ 0                # Case 3, 9, 16, 21, 13, or 24
+                f == 0 && return LSM.Action{T}(:full_return, on)         # Case 16 or 21
+                f > 0 && (k += 1)              # Case 3 or 9
+            elseif v1 > 0 && v2 ≤ 0            # Case 4, 10, 19, 20, 12, or 25
+                f == 0 && return LSM.Action{T}(:full_return, on)         # Case 19 or 20
+                f < 0 && (k += 1)              # Case 4 or 10
+            elseif v2 == 0 && v1 < 0           # Case 7, 14, or 17
+                f == 0 && return LSM.Action{T}(:full_return, on)         # Case 17
+            elseif v1 == 0 && v2 < 0           # Case 8, 15, or 18
+                f == 0 && return LSM.Action{T}(:full_return, on)         # Case 18
+            elseif v1 == 0 && v2 == 0          # Case 1, 2, 5, 6, 22, or 23
+                u2 ≤ 0 && u1 ≥ 0 && return LSM.Action{T}(:full_return, on)  # Case 1
+                u1 ≤ 0 && u2 ≥ 0 && return LSM.Action{T}(:full_return, on)  # Case 2
+            end
+            return LSM.Action(:continue, on)
+        end
+        p_start = p_end
+    end
+
+    result = SpatialTreeInterface.depth_first_search(per_edge_function,extent -> extent.Y[1] <= y <= extent.Y[2], tree)
+
+    if result isa LoopStateMachine.Action
+        return result.x
+    else
+        return iseven(k) ? out : in
+    end
+end
 _point_filled_curve_orientation(
     point, curve;
     in::T = point_in, on::T = point_on, out::T = point_out, exact,
diff --git a/src/transformations/tuples.jl b/src/transformations/tuples.jl
index 4bef4f91d..ee209a4b8 100644
--- a/src/transformations/tuples.jl
+++ b/src/transformations/tuples.jl
@@ -10,15 +10,17 @@ geometries wrapping `Tuple` points.
 
 # Keywords
 
-$APPLY_KEYWORDS
+- `threaded`: `true` or `false`. Whether to use multithreading. Defaults to `false`.
+- `crs`: The CRS to attach to geometries. Defaults to `nothing`.
+- `calc_extent`: `true` or `false`. Whether to calculate the extent. Defaults to `true`.
 """
-function tuples(geom, ::Type{T} = Float64; kw...) where T
+function tuples(geom, ::Type{T} = Float64; calc_extent = true, kw...) where T
     if _is3d(geom)
-        return apply(PointTrait(), geom; kw...) do p
+        return apply(PointTrait(), geom; calc_extent, kw...) do p
             (T(GI.x(p)), T(GI.y(p)), T(GI.z(p)))
         end
     else
-        return apply(PointTrait(), geom; kw...) do p
+        return apply(PointTrait(), geom; calc_extent, kw...) do p
             (T(GI.x(p)), T(GI.y(p)))
         end
     end
diff --git a/src/utils/NaturalIndexing.jl b/src/utils/NaturalIndexing.jl
index 080da277e..ad876e02a 100644
--- a/src/utils/NaturalIndexing.jl
+++ b/src/utils/NaturalIndexing.jl
@@ -237,8 +237,11 @@ 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)])
+    crs = GI.crs(geom)
+    return GO.apply(GI.PolygonTrait(), geom; calc_extent = true) do poly
+        rings = [GI.convert(NaturallyIndexedRing, GI.LinearRingTrait(), ring) for ring in GI.getring(poly)]
+        extent = mapreduce(GI.extent, Extents.union, rings)
+        return GI.Polygon(rings; extent, crs)
     end
 end