Upgrade to AdaptivePredicates 1.2 (#185)

NEWS.md

@@ -1,5 +1,9 @@
 # Changelog
 
+## 1.4.0
+
+- Updated to AdaptivePredicates.jl v1.2, now allowing caches to be passed to the predicates involving `incircle` and `orient3`. These are only useful when using the `AdaptiveKernel()` kernel. Outside of triangulating, these caches are not passed by default, but can be provided. The functions `get_incircle_cache` and `get_orient3_cache` can be used for this purpose on a triangulation (without a triangulation, refer to AdaptivePredicate.jl's `incircleadapt_cache` and `orient3adapt_cache`). See [#185](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/185).
+
 ## 1.3.1
 
 - Fix an issue with a weighted triangulation where the lifted points' convex hull was entirely coplanar. See [#184](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/184)

Project.toml

@@ -1,7 +1,7 @@
 name = "DelaunayTriangulation"
 uuid = "927a84f5-c5f4-47a5-9785-b46e178433df"
 authors = ["Daniel VandenHeuvel <danj.vandenheuvel@gmail.com>"]
-version = "1.3.1"
+version = "1.4.0"
 
 [deps]
 AdaptivePredicates = "35492f91-a3bd-45ad-95db-fcad7dcfedb7"
@@ -10,7 +10,7 @@ ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [compat]
-AdaptivePredicates = "1"
+AdaptivePredicates = "1.2"
 EnumX = "1.0"
 ExactPredicates = "2.2"
 Random = "1"

docs/src/extended/utils.md

@@ -33,4 +33,8 @@ DefaultAdjacentValue
 GhostVertex
 ε
 ∅
+fix_orient3_cache 
+fix_incircle_cache 
+validate_orient3_cache
+validate_incircle_cache
 ```

src/DelaunayTriangulation.jl

@@ -25,6 +25,7 @@ import EnumX
 import Random
 
 abstract type AbstractPredicateKernel end # needs to be defined early for use in data_structures.jl
+const PredicateCacheType = Union{Nothing, <:Tuple}
 
 include("data_structures.jl")
 include("predicates.jl")

src/algorithms/point_location/brute_force.jl

@@ -38,16 +38,22 @@ function brute_force_search(tri::Triangulation, q; itr = each_triangle(tri), pre
 end
 
 """
-    brute_force_search_enclosing_circumcircle(tri::Triangulation, i, predicates::AbstractPredicateKernel=AdaptiveKernel()) -> Triangle 
+    brute_force_search_enclosing_circumcircle(tri::Triangulation, i, predicates::AbstractPredicateKernel=AdaptiveKernel(); cache = nothing) -> Triangle 
 
 Searches for a triangle in `tri` containing the vertex `i` in its circumcircle using brute force. If 
 `tri` is a weighted Delaunay triangulation, the triangle returned instead has the lifted vertex `i` 
 below its witness plane. If no such triangle exists, `($∅, $∅, $∅)` is returned. You can control 
 the method used for computing predicates via the `predicates` argument.
+
+The `cache` argument is passed to [`point_position_relative_to_circumcircle`] and should be one of 
+- `nothing`: No cache is used.
+- `get_incircle_cache(tri)`: The cache stored inside `tri`.
+- `AdaptivePredicates.incircleadapt_cache(number_type(tri))`: Compute a new cache.
+The cache is only needed if an `AdaptiveKernel()` is used.
 """
-function brute_force_search_enclosing_circumcircle(tri::Triangulation, i, predicates::AbstractPredicateKernel = AdaptiveKernel())
+function brute_force_search_enclosing_circumcircle(tri::Triangulation, i, predicates::AbstractPredicateKernel = AdaptiveKernel(); cache::PredicateCacheType = nothing)
     for V in each_triangle(tri)
-        cert = point_position_relative_to_circumcircle(predicates, tri, V, i)
+        cert = point_position_relative_to_circumcircle(predicates, tri, V, i; cache)
         !is_outside(cert) && return V
     end
     tri_type = triangle_type(tri)

src/algorithms/triangulation/basic_operations/add_point.jl

@@ -67,7 +67,7 @@ function add_point!(
     end
     q = get_point(tri, new_point)
     if is_weighted(tri)
-        cert = point_position_relative_to_circumcircle(predicates, tri, V, new_point) # redirects to point_position_relative_to_witness_plane
+        cert = point_position_relative_to_circumcircle(predicates, tri, V, new_point; cache = get_orient3_cache(tri)) # redirects to point_position_relative_to_witness_plane
         is_outside(cert) && return V # If the point is submerged, then we don't add it
     end
     flag = point_position_relative_to_triangle(predicates, tri, V, q)

src/algorithms/triangulation/basic_operations/legalise_edge.jl

@@ -27,7 +27,7 @@ There is no output, as `tri` is updated in-place.
     so that at a triangle that might have appeared in both will only appear in one.
 """
 function legalise_edge!(tri::Triangulation, i, j, r, store_event_history = Val(false), event_history = nothing; predicates::AbstractPredicateKernel = AdaptiveKernel())
-    cert = is_legal(predicates, tri, i, j)
+    cert = is_legal(predicates, tri, i, j; cache = get_incircle_cache(tri))
     if is_illegal(cert)
         e = get_adjacent(tri, j, i)
         flip_edge!(tri, i, j, e, r, store_event_history, event_history)

src/algorithms/triangulation/constrained_triangulation.jl

@@ -652,7 +652,7 @@ function add_point_cavity_cdt!(tri::Triangulation, u, v, w, marked_vertices, pre
         insert_flag = true
     else
         p, q, r, s = get_point(tri, w, v, x, u) # Don't want to deal with boundary handling here 
-        incircle_test = point_position_relative_to_circle(predicates, p, q, r, s)
+        incircle_test = point_position_relative_to_circle(predicates, p, q, r, s; cache = get_incircle_cache(tri))
         orient_test = triangle_orientation(predicates, tri, u, v, w)
         insert_flag = !is_inside(incircle_test) && is_positively_oriented(orient_test)
     end

src/algorithms/triangulation/triangulate_convex.jl

@@ -121,7 +121,7 @@ to make.
 """
 function add_point_convex_triangulation!(tri::Triangulation, u, v, w, S, predicates::AbstractPredicateKernel = AdaptiveKernel())
     x = get_adjacent(tri, w, v)
-    if edge_exists(x) && is_inside(point_position_relative_to_circumcircle(predicates, tri, u, v, w, x))
+    if edge_exists(x) && is_inside(point_position_relative_to_circumcircle(predicates, tri, u, v, w, x; cache = get_incircle_cache(tri)))
         # uvw and wvx are not Delaunay 
         delete_triangle!(tri, w, v, x; protect_boundary = true, update_ghost_edges = false)
         add_point_convex_triangulation!(tri, u, v, x, S, predicates)

src/algorithms/triangulation/unconstrained_triangulation.jl

@@ -153,7 +153,7 @@ function add_point_bowyer_watson!(tri::Triangulation, new_point, initial_search_
     q = get_point(tri, _new_point)
     V = find_triangle(tri, q; predicates, m = nothing, point_indices = nothing, try_points = nothing, k = initial_search_point, rng)
     if is_weighted(tri)
-        cert = point_position_relative_to_circumcircle(predicates, tri, V, _new_point) # redirects to point_position_relative_to_witness_plane
+        cert = point_position_relative_to_circumcircle(predicates, tri, V, _new_point; cache = get_orient3_cache(tri)) # redirects to point_position_relative_to_witness_plane
         is_outside(cert) && return V # If the point is submerged, then we don't add it
     end
     flag = point_position_relative_to_triangle(predicates, tri, V, q)
@@ -298,9 +298,9 @@ Determines whether to enter the cavity in `tri` through the edge `(i, j)` when i
 function enter_cavity(tri::Triangulation, r, i, j, ℓ, predicates::AbstractPredicateKernel = AdaptiveKernel())
     contains_segment(tri, i, j)  && return false
     if is_ghost_vertex(ℓ)
-        cert = point_position_relative_to_circumcircle(tri, j, i, ℓ, r)
+        cert = point_position_relative_to_circumcircle(predicates, tri, j, i, ℓ, r; cache = get_incircle_cache(tri))
     else
-        cert = point_position_relative_to_circumcircle(tri, r, i, j, ℓ)
+        cert = point_position_relative_to_circumcircle(predicates, tri, r, i, j, ℓ; cache =  get_incircle_cache(tri))
     end
     return is_weighted(tri) ? !is_outside(cert) : is_inside(cert)
 end

src/data_structures/triangulation/methods/weights.jl

@@ -98,18 +98,24 @@ function get_power_distance(tri::Triangulation, i, j)
 end
 
 """"
-    get_distance_to_witness_plane(tri::Triangulation, i, V) -> Number
+    get_distance_to_witness_plane([kernel::AbstractPredicateKernel = AdaptiveKernel(), ] tri::Triangulation, i, V; cache = nothing) -> Number
 
 Computes the distance between the lifted companion of the vertex `i` and the witness plane to the triangle `V`. If `V` is a ghost triangle 
 and `i` is not on its solid edge, then the distance is `-Inf` if it is below the ghost triangle's witness plane and `Inf` if it is above. If `V` is a ghost triangle and `i` 
 is on its solid edge, then the distance returned is the distance associated with the solid triangle adjoining `V`.
 
 In general, the distance is positive if the lifted vertex is above the witness plane, negative if it is below, 
 and zero if it is on the plane.
+    
+The `kernel` argument determines how this result is computed, and should be 
+one of [`ExactKernel`](@ref), [`FastKernel`](@ref), and [`AdaptiveKernel`](@ref) (the default).
+See the documentation for more information about these choices.
+
+The `cache` keyword argument is passed to [`point_position_relative_to_circumcircle`]. Please see the documentation for that function for more information.
 
 See also [`point_position_relative_to_witness_plane`](@ref) and [`get_distance_to_plane`](@ref).
 """
-function get_distance_to_witness_plane(tri::Triangulation, i, V)
+function get_distance_to_witness_plane(kernel::AbstractPredicateKernel, tri::Triangulation, i, V; cache::PredicateCacheType = nothing)
     if !is_ghost_triangle(V)
         u, v, w = triangle_vertices(V)
         p⁺ = get_lifted_point(tri, u)
@@ -118,17 +124,19 @@ function get_distance_to_witness_plane(tri::Triangulation, i, V)
         s⁺ = get_lifted_point(tri, i)
         return get_distance_to_plane(p⁺, q⁺, r⁺, s⁺)
     else
-        cert = point_position_relative_to_circumcircle(tri, V, i)
+        cert = point_position_relative_to_circumcircle(kernel, tri, V, i; cache)
         if is_inside(cert)
             return -Inf
         elseif is_outside(cert)
             return Inf
         else # is_on(cert) 
             V′ = replace_ghost_triangle_with_boundary_triangle(tri, V)
-            return get_distance_to_witness_plane(tri, i, V′)
+            return get_distance_to_witness_plane(kernel, tri, i, V′; cache)
         end
     end
 end
+get_distance_to_witness_plane(tri::Triangulation, i, V; cache::PredicateCacheType = nothing) = get_distance_to_witness_plane(AdaptiveKernel(), tri, i, V; cache)
+
 
 """
     get_weighted_nearest_neighbour(tri::Triangulation, i, j = rand(each_solid_vertex(tri))) -> Vertex 
@@ -164,24 +172,32 @@ function _get_weighted_nearest_neighbour(tri, i, j)
 end
 
 @doc """
-    is_submerged(tri::Triangulation, i) -> Bool 
-    is_submerged(tri::Triangulation, i, V) -> Bool
+    is_submerged([kernel::AbstractPredicateKernel = AdaptiveKernel(), ] tri::Triangulation, i; cache = nothing) -> Bool 
+    is_submerged([kernel::AbstractPredicateKernel = AdaptiveKernel(), ] tri::Triangulation, i, V; cache = nothing) -> Bool
 
 Returns `true` if the vertex `i` is submerged in `tri` and `false` otherwise. In the 
 second method, `V` is a triangle containing `tri`.
+
+The `kernel` argument determines how this result is computed, and should be
+one of [`ExactKernel`](@ref), [`FastKernel`](@ref), and [`AdaptiveKernel`](@ref) (the default).
+See the documentation for more information about these choices.
+
+The `cache` keyword argument is passed to [`point_position_relative_to_circumcircle`]. Please see the documentation for that function for more information.
 """
 is_submerged
-function is_submerged(tri::Triangulation, i)
+function is_submerged(kernel::AbstractPredicateKernel, tri::Triangulation, i; cache::PredicateCacheType = nothing)
     # A source that mentions that testing if `i` is submerged only needs to consider the triangle that contains it 
     # is given in https://otik.uk.zcu.cz/bitstream/11025/21574/1/Zemek.pdf on p.17. 
     # (If the link dies, it is the PhD thesis of Michal Zemek, "Regular Triangulation in 3D and Its Applications".)
     is_ghost_vertex(i) && return false
     q = get_point(tri, i)
     V = find_triangle(tri, q)
-    return is_submerged(tri, i, V)
+    return is_submerged(kernel, tri, i, V; cache)
 end
-function is_submerged(tri::Triangulation, i, V)
+function is_submerged(kernel::AbstractPredicateKernel, tri::Triangulation, i, V; cache::PredicateCacheType = nothing)
     is_ghost_vertex(i) && return false
-    cert = point_position_relative_to_circumcircle(tri, V, i)
+    cert = point_position_relative_to_circumcircle(kernel, tri, V, i; cache)
     return is_outside(cert)
 end
+is_submerged(tri::Triangulation, i; cache::PredicateCacheType = nothing) = is_submerged(AdaptiveKernel(), tri, i; cache)
+is_submerged(tri::Triangulation, i, V; cache::PredicateCacheType = nothing) = is_submerged(AdaptiveKernel(), tri, i, V; cache)