Convert the orient predicate to AdaptivePredicates + minor predicate changes (#275)

- Add AdaptivePredicates as the exact backend for orient.  At some point we should switch to an ecosystem wide (DelaunayTriangulation + GeometryOps + ...) representation for predicate kernels (Exact, Adaptive, Auto (choose adaptive if available and exact if not), and Fast).
- Fix touches for multi geometries

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.


* 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.

* 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.

* 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.

Author: asinghvi17

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