Spherical triangulations, tessellation

src/DelaunayTriangulation.jl

@@ -20,6 +20,10 @@ include("validation.jl")
 include("exports.jl")
 include("public.jl")
 
+
+include("spherical/SphericalDelaunay.jl")
+
 include("precompile.jl")
 
+
 end

src/algorithms/point_location/jump_and_march.jl

@@ -763,6 +763,7 @@ function find_triangle(
         maxiters = 2 + num_exterior_curves(tri) - num_solid_vertices(tri) + num_solid_edges(tri),
         concavity_protection = false,
         use_barriers::Val{U} = Val(false),
+        check_sphere = false # this gets ignored. It's just here for spherical triangulations...
     ) where {F, U}
     I = integer_type(tri)
     maxiters = Int(maxiters)

src/algorithms/triangulation/basic_operations/add_point.jl

@@ -57,6 +57,7 @@ function add_point!(
             concavity_protection,
             predicates,
             rng,
+            check_sphere = false
         ),
         peek::P = Val(false),
     ) where {P}
@@ -97,6 +98,7 @@ function add_point!(
             concavity_protection,
             predicates,
             rng,
+            check_sphere = false
         ),
         peek::P = Val(false),
     ) where {P}

src/algorithms/triangulation/basic_operations/delete_holes.jl

@@ -171,8 +171,7 @@ function find_all_triangles_to_delete(tri::Triangulation, points_to_process)
             u, v = edge_vertices(e)
             V = construct_triangle(T, u, v, node)
             if !contains_triangle(V, triangles_to_delete)[2]
-                p, q, r = get_point(tri, u, v, node)
-                c = triangle_centroid(p, q, r)
+                c = triangle_centroid(tri, V)
                 δ = dist(tri, c)
                 if δ < zero(δ)
                     add_triangle!(triangles_to_delete, V)

src/algorithms/triangulation/constrained_triangulation.jl

@@ -717,7 +717,7 @@ function locate_intersecting_triangles(tri::Triangulation, e, rotate = Val(true)
     add_right_vertex!(history, initial(e))
     find_triangle(
         tri, get_point(tri, terminal(e)); predicates,
-        m = nothing, k = initial(e), store_history = Val(true), history, rng,
+        m = nothing, k = initial(e), store_history = Val(true), history, rng, check_sphere = false
     )
     add_left_vertex!(history, terminal(e))
     add_right_vertex!(history, terminal(e))

src/algorithms/triangulation/mesh_refinement.jl

@@ -77,10 +77,10 @@ The triangulation is refined in-place.
 function refine_itr!(tri::Triangulation, args::RefinementArguments)
     T, ρ = popfirst_triangle!(args.queue)
     u, v, w = triangle_vertices(T)
-    if !is_ghost_triangle(T) && get_adjacent(tri, u, v) == w # Need the last part in case T was already split and thus no longer exists 
+    if get_adjacent(tri, u, v) == w # Need the last part in case T was already split and thus no longer exists 
         success = split_triangle!(tri, args, T)
         if is_encroachment_failure(success)
-            !haskey(args.queue, T) && (args.queue[T] = ρ)
+            !haskey(args.queue, T) && enqueue_triangle!(tri, args.queue, ρ, T)
             split_all_encroached_segments!(tri, args)
         elseif is_successful_insertion(success)
             assess_added_triangles!(args, tri)
@@ -104,15 +104,19 @@ function enqueue_all_encroached_segments!(args::RefinementArguments, tri::Triang
     return args
 end
 
+function _each_solid_triangle(tri::Triangulation)
+    return each_solid_triangle(tri)
+end # This method exists so that we can overload it for SphericalTriangulations
+
 """
     enqueue_all_bad_triangles!(args::RefinementArguments, tri::Triangulation)
 
 Enqueues all bad triangles in the triangulation into `args.queue`.
 """
 function enqueue_all_bad_triangles!(args::RefinementArguments, tri::Triangulation)
-    for T in each_solid_triangle(tri)
+    for T in _each_solid_triangle(tri)
         ρ, flag = assess_triangle_quality(tri, args, T)
-        flag && (args.queue[T] = ρ)
+        flag && enqueue_triangle!(tri, args.queue, ρ, T)
     end
     return args
 end
@@ -183,23 +187,27 @@ Computes the Steiner point for a triangle `T` of `tri` to improve its quality in
 function get_steiner_point(tri::Triangulation, args::RefinementArguments, T)
     i, j, k = triangle_vertices(T)
     p, q, r = get_point(tri, i, j, k)
-    A = triangle_area(p, q, r)
+    A = triangle_area(tri, T)
     check_precision(A) && return Cert.PrecisionFailure, q # the point q is just for type stability with the return 
-    c = triangle_circumcenter(p, q, r, A)
+    c = triangle_circumcenter(tri, T, A)
     if !args.use_circumcenter
         # TODO: Implement generalised Steiner points so that c above is translated appropriately.
     end
     check_steiner_point_precision(tri, T, c) && return Cert.PrecisionFailure, c
     return Cert.None, c # nothing went wrong yet
 end
 
+function _is_ghost_triangle(_::Triangulation, T)
+    return is_ghost_triangle(T)
+end # This method exists so that we can overload it for SphericalTriangulations
+
 """
     check_steiner_point_precision(tri::Triangulation, T, c) -> Bool
 
 Checks if the Steiner point `c` of a triangle `T` of `tri` can be computed without precision issues, returning `true` if there are precision issues and `false` otherwise.
 """
 function check_steiner_point_precision(tri::Triangulation, T, c)
-    is_ghost_triangle(T) && return true # we aren't supposed to get ghost triangles, unless we are _just_ off the boundary due to precision errors
+    _is_ghost_triangle(tri, T) && return true # we aren't supposed to get ghost triangles, unless we are _just_ off the boundary due to precision errors
     i, j, k = triangle_vertices(T)
     (px, py), (qx, qy), (rx, ry) = get_point(tri, i, j, k)
     cx, cy = getxy(c)
@@ -252,9 +260,9 @@ function locate_steiner_point(tri::Triangulation, args::RefinementArguments, T,
     flag = point_position_relative_to_triangle(args.predicates, tri, T, c)
     !is_outside(flag) && return T, flag # T is never a ghost triangle, so don't worry about checking is_on(flag) here
     init = get_init_for_steiner_point(tri, T)
-    V, _ = find_triangle(tri, c; predicates = args.predicates, m = nothing, point_indices = nothing, try_points = nothing, k = init, args.rng, args.concavity_protection, use_barriers = Val(true))
+    V, _ = find_triangle(tri, c; predicates = args.predicates, m = nothing, point_indices = nothing, try_points = nothing, k = init, args.rng, args.concavity_protection, use_barriers = Val(true), check_sphere = false)
     flag = point_position_relative_to_triangle(args.predicates, tri, V, c)
-    if is_ghost_triangle(V) && is_on(flag)
+    if _is_ghost_triangle(tri, V) && is_on(flag)
         V = replace_ghost_triangle_with_boundary_triangle(tri, V)
     end
     return V, flag
@@ -282,10 +290,8 @@ Determines if the Steiner point `c`'s insertion will not affect the quality of `
 - `V′`: The triangle that the Steiner point is in, which is `T` if `c` is not suitable.
 """
 function check_for_invisible_steiner_point(tri::Triangulation, V, T, flag, c)
-    !is_outside(flag) && !is_ghost_triangle(V) && return c, V # don't need to check if the point is on the solid edge of a ghost triangle, since we already do that check previously via is_on(flag) in locate_steiner_point
-    i, j, k = triangle_vertices(T)
-    p, q, r = get_point(tri, i, j, k)
-    c′ = triangle_centroid(p, q, r)
+    !is_outside(flag) && !_is_ghost_triangle(tri, V) && return c, V # don't need to check if the point is on the solid edge of a ghost triangle, since we already do that check previously via is_on(flag) in locate_steiner_point
+    c′ = triangle_centroid(tri, T)
     return c′, T
 end
 
@@ -833,6 +839,16 @@ function is_triangle_nestled(tri::Triangulation, T, idx) # nestled: a triangle i
     return e_ki_seg && e_kj_seg
 end
 
+function _quality_statistics(tri::Triangulation, args::RefinementArguments, T)
+    u, v, w = triangle_vertices(T)
+    p, q, r = get_point(tri, u, v, w)
+    ℓmin², ℓmax², ℓmid², idx = squared_triangle_lengths_and_smallest_index(p, q, r)
+    A = triangle_area(tri, T)
+    cr = triangle_circumradius(tri, T, A) # this recomputes squared_triangle_lengths_and_smallest_index but it's fine
+    ρ = triangle_radius_edge_ratio(cr, sqrt(ℓmin²))
+    return ρ, A, idx
+end
+
 """
     assess_triangle_quality(tri::Triangulation, args::RefinementArguments, T) -> Float64, Bool
 
@@ -851,16 +867,12 @@ A triangle is bad quality if it does not meet the area constraints, violates the
 """
 function assess_triangle_quality(tri::Triangulation, args::RefinementArguments, T)
     haskey(args.queue, T) && return args.queue[T], true
-    if is_ghost_triangle(T)
+    if _is_ghost_triangle(tri, T)
         T = replace_ghost_triangle_with_boundary_triangle(tri, T)
     end
-    # First, get the radius-edge ratio.
+    # First, get some statistics 
     u, v, w = triangle_vertices(T)
-    p, q, r = get_point(tri, u, v, w)
-    ℓmin², ℓmed², ℓmax², idx = squared_triangle_lengths_and_smallest_index(p, q, r)
-    A = triangle_area(p, q, r)
-    cr = triangle_circumradius(A, ℓmin², ℓmed², ℓmax²)
-    ρ = triangle_radius_edge_ratio(cr, sqrt(ℓmin²))
+    ρ, A, idx = _quality_statistics(tri, args, T)
     # Next, check the area constraints.
     A < args.constraints.min_area && return ρ, false
     A > args.constraints.max_area && return ρ, true
@@ -883,11 +895,11 @@ Assesses the quality of all triangles in `args.events.added_triangles` according
 function assess_added_triangles!(args::RefinementArguments, tri::Triangulation)
     E = edge_type(tri)
     for T in each_added_triangle(args.events)
-        if !is_ghost_triangle(T)
+        if !_is_ghost_triangle(tri, T)
             u, v, w = triangle_vertices(T)
             if get_adjacent(tri, u, v) == w # It's not guaranteed that each triangle in events.added_triangles _actually_ still exists, edge flipping could have deleted it
                 ρ, flag = assess_triangle_quality(tri, args, T)
-                flag && !haskey(args.queue, T) && (args.queue[T] = ρ)
+                flag && !haskey(args.queue, T) && enqueue_triangle!(tri, args.queue, ρ, T)
                 for e in triangle_edges(T)
                     ee = construct_edge(E, initial(e), terminal(e)) # Need to convert e to the correct edge type since triangle_edges returns Tuples
                     if contains_segment(tri, ee) && !haskey(args.queue, ee) && is_encroached(tri, args, ee)

src/algorithms/triangulation/unconstrained_triangulation.jl

@@ -47,9 +47,8 @@ function get_initial_triangle(tri::Triangulation, insertion_order, predicates::A
     i, j, k = @view insertion_order[1:3] # insertion_order got converted into a Vector, so indexing is safe 
     initial_triangle = construct_positively_oriented_triangle(tri, i, j, k, predicates)
     i, j, k = triangle_vertices(initial_triangle)
-    p, q, r = get_point(tri, i, j, k)
-    degenerate_cert = triangle_orientation(predicates, p, q, r)
-    if length(insertion_order) > 3 && (is_degenerate(degenerate_cert) || check_precision(triangle_area(p, q, r))) && itr ≤ length(insertion_order) # Do not get stuck in an infinite loop if there are just three points, the three of them being collinear. The itr ≤ length(insertion_order) is needed because, if all the points are collinear, the loop could go on forever 
+    degenerate_cert = triangle_orientation(predicates, tri, i, j, k)
+    if length(insertion_order) > 3 && (is_degenerate(degenerate_cert) || check_precision(triangle_area(tri, (i, j, k)))) && itr ≤ length(insertion_order) # Do not get stuck in an infinite loop if there are just three points, the three of them being collinear. The itr ≤ length(insertion_order) is needed because, if all the points are collinear, the loop could go on forever 
         @static if VERSION ≥ v"1.8.1"
             circshift!(insertion_order, -1)
         else
@@ -151,7 +150,7 @@ The function [`add_point_bowyer_watson_dig_cavities!`](@ref) is the main workhor
 function add_point_bowyer_watson!(tri::Triangulation, new_point, initial_search_point::I, rng::Random.AbstractRNG = Random.default_rng(), predicates::AbstractPredicateKernel = AdaptiveKernel(), update_representative_point = true, store_event_history = Val(false), event_history = nothing, peek::P = Val(false)) where {I, P}
     _new_point = is_true(peek) ? new_point : I(new_point)
     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)
+    V = find_triangle(tri, q; predicates, m = nothing, point_indices = nothing, try_points = nothing, k = initial_search_point, rng, check_sphere = false)
     if is_weighted(tri)
         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

src/algorithms/voronoi/unbounded.jl

@@ -1,3 +1,18 @@
+"""
+    get_voronoi_vertex(tri::Triangulation, T) -> NTuple{2, Number}
+
+Returns the vertex of the Voronoi diagram associated with the triangle `T` in the triangulation `tri`. 
+If the triangulation is weighted, this returns the orthocenter of the triangle `T`, otherwise it returns the circumcenter.
+"""
+function get_voronoi_vertex(tri::Triangulation, T)
+    if !is_weighted(tri)
+        cx, cy = triangle_circumcenter(tri, T)
+    else
+        cx, cy = triangle_orthocenter(tri, T)
+    end
+    return cx, cy
+end
+
 """
     initialise_voronoi_tessellation(tri::Triangulation) -> VoronoiTessellation
 
@@ -28,13 +43,7 @@ function initialise_voronoi_tessellation(tri::Tr) where {Tr <: Triangulation}
         V = sort_triangle(V)
         if !is_ghost_triangle(V)
             u, v, w = triangle_vertices(V)
-            if !is_weighted(tri)
-                p, q, r = get_point(tri, u, v, w)
-                A = triangle_area(p, q, r)
-                cx, cy = triangle_circumcenter(p, q, r, A)
-            else 
-                cx, cy = triangle_orthocenter(tri, V)
-            end
+            cx, cy = get_voronoi_vertex(tri, V)
             if any(isinf, (cx, cy)) && INF_WARN[]
                 @warn "The triangle $((u, v, w)) has a degenerate circumcenter, $((cx, cy)). You may encounter issues with this tessellation. You can disable this warning using toggle_inf_warn!()."
             end

src/data_structures/mesh_refinement/refinement_arguments.jl

@@ -92,23 +92,7 @@ function RefinementArguments(
         seditious_angle,
         custom_constraint,
     )
-    queue = RefinementQueue(tri)
-    events = InsertionEventHistory(tri)
-    E = edge_type(tri)
-    I = integer_type(tri)
-    segment_list = Set{E}()
-    segment_vertices = Set{I}()
-    sizehint!(segment_list, num_edges(get_all_segments(tri)))
-    sizehint!(segment_vertices, num_edges(get_all_segments(tri)) ÷ 2)
-    for e in each_segment(tri)
-        push!(segment_list, e)
-        push!(segment_vertices, initial(e), terminal(e))
-    end
-    midpoint_split_list = Set{I}()
-    offcenter_split_list = Set{I}()
-    sizehint!(midpoint_split_list, num_edges(get_all_segments(tri)))
-    sizehint!(offcenter_split_list, 2num_edges(get_all_segments(tri)))
-    min_steiner_vertex = I(num_points(tri) + 1)
+    queue, events, segment_list, segment_vertices, midpoint_split_list, offcenter_split_list, min_steiner_vertex = _build_queues(tri)
     return RefinementArguments(
         queue,
         constraints,
@@ -129,6 +113,27 @@ function RefinementArguments(
     )
 end
 
+function _build_queues(tri::Triangulation)
+    queue = RefinementQueue(tri)
+    events = InsertionEventHistory(tri)
+    E = edge_type(tri)
+    I = integer_type(tri)
+    segment_list = Set{E}()
+    segment_vertices = Set{I}()
+    sizehint!(segment_list, num_edges(get_all_segments(tri)))
+    sizehint!(segment_vertices, num_edges(get_all_segments(tri)) ÷ 2)
+    for e in each_segment(tri)
+        push!(segment_list, e)
+        push!(segment_vertices, initial(e), terminal(e))
+    end
+    midpoint_split_list = Set{I}()
+    offcenter_split_list = Set{I}()
+    sizehint!(midpoint_split_list, num_edges(get_all_segments(tri)))
+    sizehint!(offcenter_split_list, 2num_edges(get_all_segments(tri)))
+    min_steiner_vertex = I(num_points(tri) + 1)
+    return queue, events, segment_list, segment_vertices, midpoint_split_list, offcenter_split_list, min_steiner_vertex
+end
+
 """
     is_free(args::RefinementArguments, u) -> Bool
 

src/data_structures/mesh_refinement/refinement_queue.jl

@@ -70,7 +70,6 @@ function Base.getindex(queue::RefinementQueue{T, E, F}, triangle::T) where {T, E
     end
 end
 
-
 """
     setindex!(queue::RefinementQueue{T,E,F}, ℓ²::F, segment::E) where {T,E,F}
     queue[segment] = ℓ²
@@ -107,6 +106,15 @@ function Base.setindex!(queue::RefinementQueue{T, E, F}, ρ, triangle::T) where
     return triangles
 end
 
+"""
+    enqueue_triangle!(tri::Triangulation, queue::RefinementQueue, ρ, triangle)
+
+Add `triangle` to `queue` with a priority of `ρ`. 
+""" # This method just exists so that we can overload it for SphericalTriangulations
+function enqueue_triangle!(_::Triangulation, queue::RefinementQueue, ρ, T)
+    return setindex!(queue, ρ, T)
+end
+
 """
     popfirst_segment!(queue::RefinementQueue) -> Edge 
 
@@ -140,4 +148,4 @@ has_triangles(queue::RefinementQueue) = !isempty(queue.triangles)
 
 Return `true` if `queue` has no segments or triangles, `false` otherwise.
 """
-Base.isempty(queue::RefinementQueue) = !has_segments(queue) && !has_triangles(queue)
+Base.isempty(queue::RefinementQueue) = !has_segments(queue) && !has_triangles(queue)