Tree acceleration for point-in-polygon

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

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 =

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,

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

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