Convert the orient predicate to AdaptivePredicates + minor predicate changes

Project.toml

@@ -4,6 +4,7 @@ authors = ["Anshul Singhvi <[email protected]>", "Rafael Schouten <rafaels
 version = "0.1.17"
 
 [deps]
+AdaptivePredicates = "35492f91-a3bd-45ad-95db-fcad7dcfedb7"
 CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
 DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df"
@@ -31,6 +32,7 @@ GeometryOpsProjExt = "Proj"
 GeometryOpsTGGeometryExt = "TGGeometry"
 
 [compat]
+AdaptivePredicates = "1.2"
 CoordinateTransformations = "0.5, 0.6"
 DataAPI = "1"
 DelaunayTriangulation = "1.0.4"

src/methods/clipping/predicates.jl

@@ -4,15 +4,17 @@ module Predicates
     import ExactPredicates.Codegen: group!, @genpredicate
     import GeometryOps: False, True, booltype, _tuple_point
     import GeoInterface as GI
+    import AdaptivePredicates
 
     #= Determine the orientation of c with regards to the oriented segment (a, b).
     Return 1 if c is to the left of (a, b).
     Return -1 if c is to the right of (a, b).
     Return 0 if c is on (a, b) or if a == b. =#
-    orient(a, b, c; exact) = _orient(booltype(exact), a, b, c)
+    orient(a, b, c; exact) = _orient(booltype(exact), _tuple_point(a, Float64), _tuple_point(b, Float64), _tuple_point(c, Float64))
 
     # If `exact` is `true`, use `ExactPredicates` to calculate the orientation.
-    _orient(::True, a, b, c) = ExactPredicates.orient(_tuple_point(a, Float64), _tuple_point(b, Float64), _tuple_point(c, Float64))
+    _orient(::True, a, b, c) = AdaptivePredicates.orient2p(_tuple_point(a, Float64), _tuple_point(b, Float64), _tuple_point(c, Float64))
+    # _orient(::True, a, b, c) = ExactPredicates.orient(_tuple_point(a, Float64), _tuple_point(b, Float64), _tuple_point(c, Float64))
     # If `exact` is `false`, calculate the orientation without using `ExactPredicates`.
     function _orient(exact::False, a, b, c)
         a = a .- c

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

src/methods/geom_relations/touches.jl

@@ -225,19 +225,49 @@ _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{
         GI.MultiPointTrait, GI.AbstractMultiCurveTrait,
         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