From 1d713040243a17cfb047c1e4b1727fc5e1cadb74 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 11 Oct 2025 09:36:42 -0400 Subject: [PATCH 1/4] Make regular `distance` use hypot --- src/methods/distance.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/methods/distance.jl b/src/methods/distance.jl index c4923b3c68..41c3d99308 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 = From 605f7858d9fbad43c4ba5204ed0943ff14571ecf Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sun, 12 Oct 2025 14:23:53 -0400 Subject: [PATCH 2/4] Experimental implementation of spatial tree acceleration for point in polygon This still has one Core.Box that I cant seem to get rid of... --- src/GeometryOps.jl | 2 + .../geom_relations/geom_geom_processors.jl | 49 +++++++++++++++++++ 2 files changed, 51 insertions(+) diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index 98d15af9b9..78510e4536 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -48,6 +48,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/geom_relations/geom_geom_processors.jl b/src/methods/geom_relations/geom_geom_processors.jl index 0049f993d1..b96bbf871d 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, From dea70933d25270fdbe39077a549403fa0d3edddc Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sun, 12 Oct 2025 14:24:27 -0400 Subject: [PATCH 3/4] Eagerly compute extents in `tuples` It's cheap and we should be doing it. Especially because it _really_ helps us going forward. --- src/transformations/tuples.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/transformations/tuples.jl b/src/transformations/tuples.jl index 4bef4f91dd..ee209a4b8c 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 From c187e37c1e1493bc7c05579ebfedbf6f813ab835 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sun, 12 Oct 2025 14:24:58 -0400 Subject: [PATCH 4/4] Eagerly compute extents in `NaturalIndexing.prepare_naturally`. --- src/utils/NaturalIndexing.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/utils/NaturalIndexing.jl b/src/utils/NaturalIndexing.jl index 080da277e4..ad876e02a8 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