From 5c7efea4d48ea448bece7db80462f6567abc2020 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 3 Mar 2025 20:29:10 -0500 Subject: [PATCH 1/3] Convert the orient predicate to AdaptivePredicates https://github.com/JuliaGeometry/AdaptivePredicates.jl is a library that implements adaptive predicates - substantially faster, BUT harder to construct, than exact predicates. Adaptive predicates are valid through the range of coordinates that geometries usually have, but we may need exact predicates if they are not. --- Project.toml | 2 ++ src/methods/clipping/predicates.jl | 6 ++++-- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 4d967331ae..739fc77dec 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["Anshul Singhvi ", "Rafael Schouten Date: Mon, 3 Mar 2025 20:29:48 -0500 Subject: [PATCH 2/3] Use the orient predicate, not cross, in ggp This should increase speed although we will have to benchmark. Still, at least it's benchmarkable now. --- src/methods/geom_relations/geom_geom_processors.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/methods/geom_relations/geom_geom_processors.jl b/src/methods/geom_relations/geom_geom_processors.jl index a0b9d7e99b..7e4658a7a1 100644 --- a/src/methods/geom_relations/geom_geom_processors.jl +++ b/src/methods/geom_relations/geom_geom_processors.jl @@ -505,7 +505,7 @@ function _point_filled_curve_orientation( 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.cross((u1, u2), (v1, v2); exact) + 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 on # Case 16 or 21 f > 0 && (k += 1) # Case 3 or 9 From 97015f9d33a44878ecff8bcfde8bad201a962fb6 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 3 Mar 2025 20:31:26 -0500 Subject: [PATCH 3/3] Fix touches for multi geoms Only one of the geoms in g2 MUST touch the geom in g1, all others MUST be disjoint or touching. This could potentially be accelerated by building an STRtree if the query width is long enough, or we may also need some tree interface / preparation here. E.g. Canada has a lot of islands that a tree approach could get rid of. --- src/methods/geom_relations/touches.jl | 50 +++++++++++++++++++++++---- 1 file changed, 44 insertions(+), 6 deletions(-) diff --git a/src/methods/geom_relations/touches.jl b/src/methods/geom_relations/touches.jl index 0dc5c1258a..751901f715 100644 --- a/src/methods/geom_relations/touches.jl +++ b/src/methods/geom_relations/touches.jl @@ -225,8 +225,29 @@ _touches( # # Geometries touch multi-geometry/geometry collections -#= Geometry touch a multi-geometry or a collection if the geometry touches at -least one of the elements of the collection. =# +#= + +A geometry touches a multi-geometry or a collection if the geometry touches at +least one of the elements of the collection. + +This is a bit tricky to implement - we have to actually check every geometry, +and make sure that each geom is either disjoint or touching. + +Problem here is that we would end up doing double the work. + +Either you check disjointness first, and then check touches - in which case +you have already done the work for the touches check, but can't take advantage of it. + +Or you check touches first, and if that is false, you check disjointness. But if touches failed, +and you don't know _why_ it was false (disjoint or contained / intersecting), you have to iterate +over every point twice -- again! + + +At this point we actually need a fast return function...or some more detail returned from the process functions. + +That's a project for later though. Right now we need to get this correct, so I'm going to do the dumb thing. + +=# function _touches( ::Union{GI.PointTrait, GI.AbstractCurveTrait, GI.PolygonTrait}, g1, ::Union{ @@ -234,10 +255,19 @@ function _touches( GI.MultiPolygonTrait, GI.GeometryCollectionTrait, }, g2, ) + has_touched = false for sub_g2 in GI.getgeom(g2) - !touches(g1, sub_g2) && return false + if touches(g1, sub_g2) + has_touched = true + else + # if not touching, they are either intersecting or disjoint + # if disjoint, then we can continue + # else, we can short circuit, since the geoms are not touching and not disjoint + # i.e. they are intersecting + disjoint(g1, sub_g2) || return false + end end - return true + return has_touched end # # Multi-geometry/geometry collections cross geometries @@ -251,8 +281,16 @@ function _touches( }, g1, ::GI.AbstractGeometryTrait, g2, ) + has_touched = false for sub_g1 in GI.getgeom(g1) - !touches(sub_g1, g2) && return false + if touches(sub_g1, g2) + has_touched = true + else + # if not touching, they are either intersecting or disjoint + # if disjoint, then we can continue + # else, we can short circuit, since the geoms are not touching and not disjoint + disjoint(sub_g1, g2) || return false + end end - return true + return has_touched end \ No newline at end of file