From 0e4df032aeb9ca2104699eecbd7ab90e9ec725b1 Mon Sep 17 00:00:00 2001 From: Daniel VandenHeuvel <95613936+DanielVandH@users.noreply.github.com> Date: Sun, 22 Sep 2024 12:51:49 +0100 Subject: [PATCH 1/6] Spherical triangulation --- Project.toml | 2 +- src/DelaunayTriangulation.jl | 2 + .../basic_operations/delete_holes.jl | 3 +- .../triangulation/mesh_refinement.jl | 54 +++-- .../unconstrained_triangulation.jl | 5 +- src/algorithms/voronoi/unbounded.jl | 23 +- .../mesh_refinement/refinement_arguments.jl | 39 +-- .../mesh_refinement/refinement_queue.jl | 12 +- .../individual_triangle_statistics.jl | 69 ++++-- src/spherical/SphericalDelaunay.jl | 108 +++++++++ src/spherical/geometry.jl | 177 ++++++++++++++ src/spherical/plotting.jl | 226 ++++++++++++++++++ src/spherical/refine.jl | 97 ++++++++ src/spherical/spherical_triangles.jl | 111 +++++++++ src/spherical/triangulation.jl | 117 +++++++++ src/spherical/utils.jl | 48 ++++ src/spherical/voronoi.jl | 32 +++ src/utils/utils.jl | 4 +- src/validation.jl | 9 +- test/Project.toml | 2 + test/runtests.jl | 4 + test/spherical.jl | 207 ++++++++++++++++ 22 files changed, 1265 insertions(+), 86 deletions(-) create mode 100644 src/spherical/SphericalDelaunay.jl create mode 100644 src/spherical/geometry.jl create mode 100644 src/spherical/plotting.jl create mode 100644 src/spherical/refine.jl create mode 100644 src/spherical/spherical_triangles.jl create mode 100644 src/spherical/triangulation.jl create mode 100644 src/spherical/utils.jl create mode 100644 src/spherical/voronoi.jl create mode 100644 test/spherical.jl diff --git a/Project.toml b/Project.toml index b01604a23..c42565846 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "DelaunayTriangulation" uuid = "927a84f5-c5f4-47a5-9785-b46e178433df" authors = ["Daniel VandenHeuvel "] -version = "1.3.1" +version = "1.4.0" [deps] AdaptivePredicates = "35492f91-a3bd-45ad-95db-fcad7dcfedb7" diff --git a/src/DelaunayTriangulation.jl b/src/DelaunayTriangulation.jl index ab6c7c0ab..6e7f1262f 100644 --- a/src/DelaunayTriangulation.jl +++ b/src/DelaunayTriangulation.jl @@ -35,4 +35,6 @@ include("validation.jl") include("exports.jl") include("public.jl") +include("spherical/SphericalDelaunay.jl") + end diff --git a/src/algorithms/triangulation/basic_operations/delete_holes.jl b/src/algorithms/triangulation/basic_operations/delete_holes.jl index c5fcbc1fb..7f26c44d2 100644 --- a/src/algorithms/triangulation/basic_operations/delete_holes.jl +++ b/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) diff --git a/src/algorithms/triangulation/mesh_refinement.jl b/src/algorithms/triangulation/mesh_refinement.jl index d27422305..41cc9969e 100644 --- a/src/algorithms/triangulation/mesh_refinement.jl +++ b/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,9 +187,9 @@ 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 @@ -193,13 +197,17 @@ function get_steiner_point(tri::Triangulation, args::RefinementArguments, T) 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) @@ -254,7 +262,7 @@ function locate_steiner_point(tri::Triangulation, args::RefinementArguments, T, 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)) 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) diff --git a/src/algorithms/triangulation/unconstrained_triangulation.jl b/src/algorithms/triangulation/unconstrained_triangulation.jl index 7b08eb782..1378792bc 100644 --- a/src/algorithms/triangulation/unconstrained_triangulation.jl +++ b/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 diff --git a/src/algorithms/voronoi/unbounded.jl b/src/algorithms/voronoi/unbounded.jl index 01d46cee3..01d689b57 100644 --- a/src/algorithms/voronoi/unbounded.jl +++ b/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 diff --git a/src/data_structures/mesh_refinement/refinement_arguments.jl b/src/data_structures/mesh_refinement/refinement_arguments.jl index 18768a2f0..9c79f4524 100644 --- a/src/data_structures/mesh_refinement/refinement_arguments.jl +++ b/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 diff --git a/src/data_structures/mesh_refinement/refinement_queue.jl b/src/data_structures/mesh_refinement/refinement_queue.jl index bc48de5ce..9c0497dbb 100644 --- a/src/data_structures/mesh_refinement/refinement_queue.jl +++ b/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) \ No newline at end of file diff --git a/src/data_structures/statistics/individual_triangle_statistics.jl b/src/data_structures/statistics/individual_triangle_statistics.jl index 3c955b870..042adfeb9 100644 --- a/src/data_structures/statistics/individual_triangle_statistics.jl +++ b/src/data_structures/statistics/individual_triangle_statistics.jl @@ -123,6 +123,17 @@ function triangle_area(ℓ₁²::Number, ℓ₂²::Number, ℓ₃²::Number) return sqrt(A²) end +""" + triangle_area(tri::Triangulation, T) -> Number + +Computes the area of the triangle `T` in the triangulation `tri`. +""" +function triangle_area(tri::Triangulation, T) + i, j, k = triangle_vertices(T) + p, q, r = get_point(tri, i, j, k) + return triangle_area(p, q, r) +end + @doc raw""" triangle_circumradius(A, ℓmin², ℓmed², ℓmax²) -> Number @@ -132,7 +143,7 @@ Computes the circumradius of a triangle with area `A` and squared edge lengths ` r = \dfrac{\ell_{\min}\ell_{\text{med}}\ell_{\max}}{4A}. ``` """ -triangle_circumradius(A, ℓmin², ℓmed², ℓmax²) = sqrt(ℓmin² * ℓmed² * ℓmax²) / (4A) +triangle_circumradius(A::Number, ℓmin²::Number, ℓmed²::Number, ℓmax²::Number) = sqrt(ℓmin² * ℓmed² * ℓmax²) / (4A) @doc raw""" triangle_perimeter(ℓmin::Number, ℓmed::Number, ℓmax::Number) -> Number @@ -199,22 +210,17 @@ function triangle_centroid(p, q, r) cx′, cy′ = 2midpoint(qx′, rx′) / 3, 2midpoint(qy′, ry′) / 3 cx, cy = cx′ + px, cy′ + py return cx, cy - #= - A = triangle_area(PP, QQ, RR) - isixA = inv(6A) - px′, py′ = zero(px), zero(py) - qx′, qy′ = qx - px, qy - py - rx′, ry′ = rx - px, ry - py - cx1 = (px′ + qx′) * ((px′ * qy′) - (qx′ * py′)) - cx2 = (qx′ + rx′) * ((qx′ * ry′) - (rx′ * qy′)) - cx3 = (rx′ + px′) * ((rx′ * py′) - (px′ * ry′)) - cx = px + (cx1 + cx2 + cx3) * isixA - cy1 = (py′ + qy′) * ((px′ * qy′) - (qx′ * py′)) - cy2 = (qy′ + ry′) * ((qx′ * ry′) - (rx′ * qy′)) - cy3 = (ry′ + py′) * ((rx′ * py′) - (px′ * ry′)) - cy = py + (cy1 + cy2 + cy3) * isixA - return cx, cy - =# +end + +""" + triangle_centroid(tri::Triangulation, T) -> (Number, Number) + +Computes the centroid of the triangle `T` in the triangulation `tri`. +""" +function triangle_centroid(tri::Triangulation, T) + i, j, k = triangle_vertices(T) + p, q, r = get_point(tri, i, j, k) + return triangle_centroid(p, q, r) end @doc raw""" @@ -337,27 +343,40 @@ function triangle_circumcenter(p, q, r, A = triangle_area(p, q, r)) end """ - triangle_circumcenter(tri::Triangulation, T) -> (Number, Number) + triangle_circumcenter(tri::Triangulation, T[, A = triangle_area(tri, T)]) -> (Number, Number) -Computes the circumcenter of the triangle `T` in the triangulation `tri`. +Computes the circumcenter of the triangle `T` in the triangulation `tri`. You +can optionally provide the area `A` of the triangle. """ -function triangle_circumcenter(tri::Triangulation, T) +function triangle_circumcenter(tri::Triangulation, T, A = triangle_area(tri, T)) i, j, k = triangle_vertices(T) p, q, r = get_point(tri, i, j, k) - return triangle_circumcenter(p, q, r) + return triangle_circumcenter(p, q, r, A) end """ - triangle_circumradius(p, q, r) -> Number + triangle_circumradius(p, q, r[, A=triangle_area(p, q, r)]) -> Number Computes the circumradius of the triangle with coordinates `(p, q, r)`. +You can optionally provide the area `A` of the triangle. """ -function triangle_circumradius(p, q, r) +function triangle_circumradius(p, q, r, A = triangle_area(p, q, r)) ℓ₁², ℓ₂², ℓ₃² = squared_triangle_lengths(p, q, r) - A = triangle_area(p, q, r) return triangle_circumradius(A, ℓ₁², ℓ₂², ℓ₃²) end +""" + triangle_circumradius(tri::Triangulation, T, A = triangle_area(tri, T)) -> Number + +Computes the circumradius of the triangle `T` in the triangulation `tri`. You +can optionally provide the area `A` of the triangle. +""" +function triangle_circumradius(tri::Triangulation, T, A = triangle_area(tri, T)) + i, j, k = triangle_vertices(T) + p, q, r = get_point(tri, i, j, k) + return triangle_circumradius(p, q, r, A) +end + """ triangle_offcenter(p, q, r, c₁=triangle_circumcenter(p, q, r), β=1.0) -> (Number, Number) @@ -622,7 +641,7 @@ function triangle_sink(tri::Triangulation, T, prev_T = construct_triangle(triang =# _, _, θ₃ = triangle_angles(p, q, r) θ₃ ≤ π / 2 + ε(θ₃) && return c - m = triangle_centroid(p, q, r) + m = triangle_centroid(tri, T) if !is_none(line_segment_intersection_type(predicates, p, q, m, c)) && !is_left(point_position_relative_to_line(predicates, p, q, c)) next_T = construct_triangle(triangle_type(tri), j, i, get_adjacent(tri, j, i)) elseif !is_none(line_segment_intersection_type(predicates, q, r, m, c)) && !is_left(point_position_relative_to_line(predicates, q, r, c)) diff --git a/src/spherical/SphericalDelaunay.jl b/src/spherical/SphericalDelaunay.jl new file mode 100644 index 000000000..57b7a195d --- /dev/null +++ b/src/spherical/SphericalDelaunay.jl @@ -0,0 +1,108 @@ +""" + SphericalDelaunay + +This module contains functions for computing Delaunay triangulations and Voronoi tessellations on the +surface of a sphere. Currently, all spheres are assumed to be centered at the origin and have a unit radius. + +This module is not yet entirely developed, and so you should consider it experimental. +""" +module SphericalDelaunay + +import ..DelaunayTriangulation: + number_type, + getxyz, + getxy, + getx, + gety, + get_voronoi_vertex, + Triangulation, + triangle_circumcenter, + get_points, + triangle_vertices, + sort_triangle, + triangle_area, + triangle_circumradius, + det, + each_point, + triangulate, + set_point!, + push_point!, + convex_hull!, + AbstractParametricCurve, + norm_sqr, + get_representative_point_list, + RepresentativeCoordinates, + integer_type, + empty_representative_points!, + new_representative_point!, + delete_ghost_triangles!, + num_points, + delete_ghost_vertices_from_graph!, + each_ghost_triangle, + add_triangle!, + each_solid_triangle, + num_solid_triangles, + num_triangles, + each_triangle, + is_planar, + _quality_statistics, + RefinementArguments, + triangle_radius_edge_ratio, + triangle_sink, + triangle_centroid, + AbstractPredicateKernel, + num_solid_vertices, + AdaptiveKernel, + get_steiner_point, + Certificate, + enqueue_triangle!, + RefinementQueue, + RefinementConstraints, + _build_queues, + check_steiner_point_precision, + finalise!, + _each_solid_triangle, + has_ghost_triangles, + _is_ghost_triangle, + locate_steiner_point, + voronoi, + VoronoiTessellation, + is_ghost_vertex, + get_circumcenter_to_triangle, + push_polygon_point!, + get_triangle_to_circumcenter, + num_polygons, + get_polygon_points, + get_triangulation, + num_polygon_vertices, + each_polygon, + each_polygon_index, + get_polygon, + get_polygon_point + +import Random: + Random, + rand, + AbstractRNG, + SamplerTrivial + +export SphericalPoint, + SphericalTriangulation, + SphericalTessellation, + spherical_triangulate, + spherical_voronoi, + UnitSphere, + SphericalLine, + SphericalTriangle, + SphericalPolygon, + solidify! + +include("geometry.jl") +include("utils.jl") +include("spherical_triangles.jl") +include("triangulation.jl") +include("refine.jl") +include("voronoi.jl") +include("plotting.jl") + +end # module \ No newline at end of file diff --git a/src/spherical/geometry.jl b/src/spherical/geometry.jl new file mode 100644 index 000000000..008744d26 --- /dev/null +++ b/src/spherical/geometry.jl @@ -0,0 +1,177 @@ +""" + abstract type AbstractSurface{T} + +An abstract type for surfaces. +""" +abstract type AbstractSurface{T} end + +""" + struct UnitSphere{T} <: AbstractSurface{T} + +A unit sphere centered at the origin. The type parameter `T` is the number type. +The constructor `UnitSphere()` creates a `UnitSphere{Float64}`. +""" +struct UnitSphere{T} <: AbstractSurface{T} end +UnitSphere() = UnitSphere{Float64}() + +function Random.eltype(::Type{UnitSphere{T}}) where {T} + return SphericalPoint{NTuple{3,T}} +end +function rand(rng::AbstractRNG, ::SamplerTrivial{UnitSphere{T}}) where {T} + x, y, z = ntuple(_ -> randn(rng, T), Val(3)) + r⁻¹ = 1 / _hypot(x, y, z) + return SphericalPoint((x * r⁻¹, y * r⁻¹, z * r⁻¹)) +end + +""" + struct SphericalPoint{P} + +Wrapper for representing a point on the surface of a [`UnitSphere`](@ref). The type parameter `P` is type of the wrapped point. + +# Fields +- `p::P`: The point on the sphere. + +# Constructors +- `SphericalPoint(x, y, z)`: Creates a `SphericalPoint` from the Cartesian coordinates `x`, `y`, and `z`. +- `SphericalPoint(p::SphericalPoint)`: Returns `p`. +- `SphericalPoint(latitude, longitude)`: Creates a `SphericalPoint` from the latitude and longitude. The resulting point is still on the unit sphere. The angles are assumed to be in degrees. +""" +struct SphericalPoint{P} + p::P +end +SphericalPoint(x, y, z) = SphericalPoint((x, y, z)) +SphericalPoint(p::SphericalPoint) = p +function SphericalPoint(latitude, longitude) + lat, lon = latitude, longitude + slat, clat = sincosd(lat) + slon, clon = sincosd(lon) + x = clat * clon + y = clat * slon + z = slat + return SphericalPoint(x, y, z) +end +Base.show(io::IO, p::SphericalPoint) = print(io, "SphericalPoint", __getxyz(p), "") + + +Base.:(==)(p::SphericalPoint, q::SphericalPoint) = __getxyz(p) == __getxyz(q) + +""" + getp(p::SphericalPoint{P}) -> P + +Returns the point `p` wrapped by `SphericalPoint`. +""" +getp(p::SphericalPoint) = p.p + +""" + __getxyz(p::SphericalPoint) -> Tuple{Number, Number, Number} + +Returns the Cartesian coordinates of the point `p` wrapped by `SphericalPoint`. +""" +__getxyz(p::SphericalPoint) = getxyz(getp(p)) + +""" + projectx(p::SphericalPoint) -> Number + +Computes the `x`-coordinate of the stereographic projection of `p`, i.e. returns ``x/(1-z)``, where ``p = (x, y, z)``. +""" +function projectx(p::SphericalPoint) + x, _, z = __getxyz(p) + return x / (1 - z) +end + +""" + projecty(p::SphericalPoint) -> Number + +Computes the `y`-coordinate of the stereographic projection of `p`, i.e. returns ``y/(1-z)``, where ``p = (x, y, z)``. +""" +function projecty(p::SphericalPoint) + _, y, z = __getxyz(p) + return y / (1 - z) +end + +""" + project(p::SphericalPoint) -> Tuple{Number, Number} + +Computes the stereographic projection of `p`, i.e. returns the `Tuple` `(projectx(p), projecty(p))`. +""" +project(p::SphericalPoint) = (projectx(p), projecty(p)) + +""" + invproject(p::Tuple{Number, Number}) -> SphericalPoint + +Computes the inverse stereographic projection of `p`, i.e. returns the point on the sphere corresponding to `p`. The +wrapped point is an `NTuple{3, Number}`. +""" +function invproject(p) + X, Y = getxy(p) + if isnan(X) && isnan(Y) + return north_pole(number_type(p)) + end + normsq = norm_sqr(p) + den = 1 + normsq + p = (2X / den, 2Y / den, (normsq - 1) / den) + return SphericalPoint(p) +end + +#= +""" + getx(p::SphericalPoint) -> Number + +Returns the `x`-coordinate of the stereographic projection of `p`. Note that this is +*not* the `x`-coordinate of the point `p` on the sphere. +""" +=# +getx(p::SphericalPoint) = projectx(p) + +#= +""" + gety(p::SphericalPoint) -> Number + +Returns the `y`-coordinate of the stereographic projection of `p`. Note that this is +*not* the `y`-coordinate of the point `p` on the sphere. +""" +=# +gety(p::SphericalPoint) = projecty(p) + +#= +""" + number_type(::Type{SphericalPoint{P}}) -> T + +Returns the number type of the point `P` wrapped by `SphericalPoint`. +""" +=# +number_type(::Type{SphericalPoint{P}}) where {P} = number_type(P) + +""" + struct SphericalPoints{P} + +Wrapper for representing a collection of points on the surface of a [`UnitSphere`](@ref). The type parameter `P` is the type of the wrapped points. +""" +struct SphericalPoints{T,P} <: AbstractVector{SphericalPoint{T}} + points::P + function SphericalPoints(points::AbstractVector{<:SphericalPoint{P}}) where {P} + return new{P,typeof(points)}(points) + end +end +SphericalPoints(p::SphericalPoints) = p + +Base.eachindex(points::SphericalPoints) = Base.eachindex(points.points) +Base.iterate(points::SphericalPoints, state...) = Base.iterate(points.points, state...) +Base.size(points::SphericalPoints) = size(points.points) +Base.length(points::SphericalPoints) = length(points.points) +Base.getindex(points::SphericalPoints, i) = SphericalPoint(points.points[i]) +number_type(::Type{SphericalPoints{P}}) where {P} = number_type(P) +set_point!(points::SphericalPoints, i, x, y) = setindex!(points.points, invproject((x, y)), i) +Base.pop!(points::SphericalPoints) = pop!(points.points) +push_point!(points::SphericalPoints, x, y) = push!(points.points, invproject((x, y))) +is_planar(points::SphericalPoints) = true # hack to avoid warnings + +north_pole(::Type{T}) where {T} = SphericalPoint((zero(T), zero(T), one(T))) + +function _getindex(points::SphericalPoints, i) + if is_ghost_vertex(i) + return north_pole(number_type(points)) + else + return getindex(points, i) + end +end \ No newline at end of file diff --git a/src/spherical/plotting.jl b/src/spherical/plotting.jl new file mode 100644 index 000000000..b63af1928 --- /dev/null +++ b/src/spherical/plotting.jl @@ -0,0 +1,226 @@ +""" + SphericalLine(p, q) + +A line on the surface of a sphere defined by the points `p` and `q`. The line is parametrised by the angle `θ` from `p` to `q`. +You can evaluate this struct like a function. +""" +struct SphericalLine{P,T} <: AbstractParametricCurve + p::SphericalPoint{P} + q::SphericalPoint{P} + upper::T + orth::NTuple{3,T} +end +function SphericalLine(p::SphericalPoint, q::SphericalPoint) + # https://comp.soft-sys.matlab.narkive.com/lSloCHDO/line-between-two-points-on-a-sphere + pq_cross = cross(p, q) + pq_orth = cross(SphericalPoint(pq_cross), p) + pq_orth_norm = _hypot(pq_orth...) + pq_orth = pq_orth ./ pq_orth_norm + pq_dot = dot(p, q) + pq_cross_norm = _hypot(pq_cross...) + ϕ = mod2pi(atan(pq_cross_norm, pq_dot)) + return SphericalLine(p, q, ϕ, pq_orth) +end +function (L::SphericalLine)(t) + s, c = sincos(t) + p = getp(L.p) + orth = L.orth + coords = p .* c .+ orth .* s + return SphericalPoint(coords) +end + +""" + _to_coords(L::SphericalLine[, n = 20]) -> Vector{NTuple{3, Number}} + +Returns the coordinates of the points on the line `L` as a vector of `NTuple{3, Number}`. +The line is evaluated at `n` equally spaced points. +""" +function _to_coords(L::SphericalLine, n) + upper = L.upper + t = LinRange(zero(upper), upper, n) + points = L.(t) + return points +end + +""" + SphericalTriangle(p, q, r) + +A triangle on the surface of a sphere defined by the points `p`, `q`, and `r`. The triangle is stored +as three [`SphericalLine`](@ref)s. +""" +struct SphericalTriangle{P,T} + L1::SphericalLine{P,T} + L2::SphericalLine{P,T} + L3::SphericalLine{P,T} +end +function SphericalTriangle(p::SphericalPoint, q::SphericalPoint, r::SphericalPoint) + return SphericalTriangle( + SphericalLine(p, q), + SphericalLine(q, r), + SphericalLine(r, p) + ) +end + +""" + _to_coords(T::SphericalTriangle, add_end=true, n=20) -> Vector{NTuple{3, Number}} + +Returns the coordinates of the points on the triangle `T` as a vector of `NTuple{3, Number}`. If `add_end` is `true`, +the first point is repeated at the end of the vector. + +The triangle is evaluated at `n` equally spaced points on each line. +""" +function _to_coords(T::SphericalTriangle, add_end=true, n=20) + coords1 = _to_coords(T.L1, n) + _coords2 = _to_coords(T.L2, n) + coords2 = view(_coords2, 2:length(_coords2)) + _coords3 = _to_coords(T.L3, n) + coords3 = view(_coords3, 2:(length(_coords3)-1)) + coords = vcat(coords1, coords2, coords3) + add_end && push!(coords, coords[1]) + return coords +end + +""" + _to_mesh(points::Vector{<:SphericalPoint}[, offset=0]) -> Matrix{Int} + +Given a closed curve on the surface of a [`UnitSphere`](@ref) defined by `points`, returns a +triangulation of that curve appropriate for use with Makie.jl's `mesh`. + +The optional argument `offset` is the index of the first point in the curve. +""" +function _to_mesh(points::Vector{<:SphericalPoint}, offset=0) + @assert points[begin] != points[end] + tri = spherical_triangulate(points) + triangles = Matrix{Int}(undef, num_solid_triangles(tri), 3) + for (i, T) in enumerate(each_solid_triangle(tri)) + u, v, w = triangle_vertices(T) + triangles[i, :] .= (u, v, w) .+ offset + end + return triangles +end + +""" + SphericalTriangle(tri::Triangulation, T) + +Given a Delaunay triangulation `tri` and a triangle `T`, returns the `SphericalTriangle` corresponding to `T`. +""" +function SphericalTriangle(tri::Triangulation, T) + points = get_points(tri) + V = number_type(tri) + T = sort_triangle(T) + u, v, w = triangle_vertices(T) + p, q, r = _getindex(points, u), _getindex(points, v), _getindex(points, w) # can't do get_point since DelaunayTriangulation projects it + p, q, r = fix_north_pole(p, q, r) + ST = SphericalTriangle(p, q, r) + return ST +end + +""" + _to_triangles(tri::SphericalTriangulation[, n = 20]) -> Vector{<:SphericalPoint}, Matrix{Int} + +Given a spherical triangulation `tri`, returns the points and triangles of the triangulation as a vector of `SphericalPoint`s +and a matrix of indices. This format is appropriate for use with Makie.jl's `mesh`. `n` is the number of points to use for discretising +each line of each triangle. +""" +function _to_triangles(tri::SphericalTriangulation, n = 20) + points = get_points(tri) + mesh_points = eltype(points)[] + vertices = Matrix{Int}(undef, 0, 3) + offset = 0 + for T in each_triangle(tri) + ST = SphericalTriangle(tri, T) + coords = _to_coords(ST, false, n) + _vertices = _to_mesh(coords, offset) + offset += length(coords) + vertices = vcat(vertices, _vertices) + append!(mesh_points, coords) + end + return mesh_points, vertices +end + +""" + _to_lines(tri::SphericalTriangulation[, n = 20]) -> Vector{<:SphericalPoint} + +Given a spherical triangulation `tri`, returns the points on the lines of the triangulation as a vector of `SphericalPoint`s. +Each triangle's lines are separated by a `SphericalPoint(NaN, NaN, NaN)`. +""" +function _to_lines(tri::SphericalTriangulation, n = 20) + points = get_points(tri) + nt = num_triangles(tri) + lines = Vector{eltype(points)}(undef, 3n*(nt+1)) # 3n points per triangle, plus one for the NaN + offset = 0 + for T in each_triangle(tri) + ST = SphericalTriangle(tri, T) + coords = _to_coords(ST, true, n) + lines[offset+1:offset+length(coords)] .= coords + lines[offset+length(coords)+1] = SphericalPoint(NaN, NaN, NaN) + offset += length(coords) + 1 + end + return lines +end + +""" + SphericalPolygon(lines::Vector{<:SphericalLine}) + +A polygon on the surface of a sphere defined by the lines `lines`. The polygon is stored as a vector of `SphericalLine`s. +""" +struct SphericalPolygon{P,T} + lines::Vector{SphericalLine{P,T}} +end +function SphericalPolygon(points::Vector{<:SphericalPoint}) + n = length(points) + @assert points[begin] != points[end] + lines = map(eachindex(points)) do i + SphericalLine(points[i], points[mod1(i+1, n)]) + end + return SphericalPolygon(lines) +end + +""" + SphericalPolygon(vorn::SphericalTessellation, v) + +Given a `SphericalTessellation` `vorn` and a polygon index `v`, returns the [`SphericalPolygon`](@ref) corresponding to the Voronoi polygon of `v`. +""" +function SphericalPolygon(vorn::SphericalTessellation, v) + polygon = get_polygon(vorn, v) + points = [get_polygon_point(vorn, i) for i in view(polygon, 1:(length(polygon)-1))] + return SphericalPolygon(invproject.(points)) +end + +""" + _to_coords(P::SphericalPolygon, add_end=true, n=20) -> Vector{NTuple{3, Number}} + +Returns the coordinates of the points on the polygon `P` as a vector of `NTuple{3, Number}`. The polygon is evaluated at `n` equally spaced points on each line. +The first point is repeated at the end of the vector if `add_end` is `true`. +""" +function _to_coords(P::SphericalPolygon, add_end=true, n=20) + coords = Vector{eltype(P.lines)}() + for L in P.lines + _coords = _to_coords(L, n) + coords = vcat(coords, _coords) + end + add_end && push!(coords, coords[1]) + return coords +end + +""" + _to_lines(vorn::SphericalTessellation, n=20) -> Vector{<:SphericalPoint} + +Given a `SphericalTessellation` `vorn`, returns the points on the lines of the tessellation as a vector of `SphericalPoint`s. +Each polygon's lines are separated by a `SphericalPoint(NaN, NaN, NaN)`. +""" +function _to_lines(vorn::SphericalTessellation, n = 20) + tri = get_triangulation(vorn) + spoints = get_points(tri) + polygon_points = get_polygon_points(vorn) + nt = num_polygons(vorn) + lines = Vector{eltype(spoints)}(undef, 0) + sizehint!(lines, (6*3)*n*(nt+1)) + for v in each_polygon_index(vorn) + spoly = SphericalPolygon(vorn, v) + coords = _to_coords(spoly, true, n) + append!(lines, coords) + push!(lines, SphericalPoint(NaN, NaN, NaN)) + end + return lines +end \ No newline at end of file diff --git a/src/spherical/refine.jl b/src/spherical/refine.jl new file mode 100644 index 000000000..78a195776 --- /dev/null +++ b/src/spherical/refine.jl @@ -0,0 +1,97 @@ +function _quality_statistics(tri::SphericalTriangulation, args::RefinementArguments, T) + u, v, w = triangle_vertices(T) + A = triangle_area(tri, T) + cr = triangle_circumradius(tri, T, A) + points = get_points(tri) + p, q, r = _getindex(points, u), _getindex(points, v), _getindex(points, w) # can't do get_point since DelaunayTriangulation projects it + d1, d2, d3 = spherical_distance(p, q), spherical_distance(q, r), spherical_distance(r, p) + mind = min(d1, d2, d3) + idx = mind == d1 ? 1 : mind == d2 ? 2 : 3 + ρ = triangle_radius_edge_ratio(cr, mind) + return ρ, A, idx +end + +function triangle_sink(tri::SphericalTriangulation, args...) # just avoid computing this in statistics() + V = number_type(tri) + return (V(NaN), V(NaN)) +end + +function RefinementArguments( + tri::SphericalTriangulation; + min_angle = 0.0, + max_angle = 180.0, + min_area = number_type(tri)(4π / 1.0e9), + max_area = typemax(number_type(tri)), + max_points = max(1_000, num_solid_vertices(tri))^2, + seditious_angle = -Inf, + custom_constraint = (_tri, T) -> false, + use_circumcenter = true, + use_lens = true, + steiner_scale = 0.999, + rng = Random.default_rng(), + concavity_protection = false, + predicates::AbstractPredicateKernel = AdaptiveKernel() +) + # Ignoring angle constraints. Spherical refinement uses centroids instead of Ruppert's algorithm. + has_ghosts = has_ghost_triangles(tri) + if !has_ghosts + throw(ArgumentError("Spherical triangulations must have ghost triangles for refinement to work. Please avoid using solidify! for this.")) + end + constraints = RefinementConstraints(; + min_angle = 0.0, + max_angle = 180.0, + min_area, + max_area, + max_points, + seditious_angle=-Inf, + custom_constraint) + queue, events, segment_list, segment_vertices, midpoint_split_list, offcenter_split_list, min_steiner_vertex = _build_queues(tri) + return RefinementArguments( + queue, + constraints, + events, + min_steiner_vertex, + segment_list, + segment_vertices, + midpoint_split_list, + offcenter_split_list, + use_circumcenter, + use_lens, + steiner_scale, + false, + true, + rng, + concavity_protection, + predicates + ) +end + +function get_steiner_point(tri::SphericalTriangulation, args::RefinementArguments, T) + c = triangle_centroid(tri, T) + check_steiner_point_precision(tri, T, c) && return Certificate.PrecisionFailure, c + return Certificate.None, c +end + +function enqueue_triangle!(tri::SphericalTriangulation, queue::RefinementQueue, _, T) + A = triangle_area(tri, T) # Instead of using radius-edge ratios, spherical triangulations use areas since we are doing centroid refinement + return setindex!(queue, A, T) +end + +function finalise!(tri::Triangulation, args::RefinementArguments) + #unlock_convex_hull!(tri; reconstruct = true) + #lock_convex_hull!(tri; rng = args.rng, predicates = args.predicates) + convex_hull!(tri; reconstruct = false, predicates = args.predicates) + return tri +end + +function _each_solid_triangle(tri::SphericalTriangulation) + return each_triangle(tri) +end + +function _is_ghost_triangle(_::SphericalTriangulation, T) + return false +end + +function locate_steiner_point(tri::SphericalTriangulation, args::RefinementArguments, T, c) + return T, Certificate.Inside # centroids are always inside the triangle +end \ No newline at end of file diff --git a/src/spherical/spherical_triangles.jl b/src/spherical/spherical_triangles.jl new file mode 100644 index 000000000..37dfff02e --- /dev/null +++ b/src/spherical/spherical_triangles.jl @@ -0,0 +1,111 @@ +""" + spherical_circumcenter(p::SphericalPoint, q::SphericalPoint, r::SphericalPoint) -> SphericalPoint + +Computes the circumcenter of the spherical triangle defined by the points `p`, `q`, and `r`. + +See https://brsr.github.io/2021/05/02/spherical-triangle-centers.html +""" +function spherical_circumcenter(p::SphericalPoint, q::SphericalPoint, r::SphericalPoint) + num = cross(p, q) .+ cross(q, r) .+ cross(r, p) + den = _hypot(num...) + return SphericalPoint(num ./ den) +end + +""" + spherical_triangle_area(p::SphericalPoint, q::SphericalPoint, r::SphericalPoint) -> Number + +Computes the area of the spherical triangle defined by the points `p`, `q`, and `r`. + +See https://brsr.github.io/2023/11/04/spherical-areal.html. +""" +function spherical_triangle_area(p::SphericalPoint, q::SphericalPoint, r::SphericalPoint) + num = abs(dot(p, SphericalPoint(cross(q, r)))) + den = 1 + dot(p, q) + dot(q, r) + dot(r, p) + rat = num / den + A = 2atan(num, den) + return A +end + +""" + spherical_triangle_circumradius(p::SphericalPoint, q::SphericalPoint, r::SphericalPoint) -> Number + +Computes the circumradius of the spherical triangle defined by the points `p`, `q`, and `r`. + +See https://brsr.github.io/2021/05/02/spherical-triangle-centers.html. +""" +function spherical_triangle_circumradius(p, q, r) + num = scalar_triple_product(p, q, r) + _den = cross(p, q) .+ cross(q, r) .+ cross(r, p) + den = _hypot(_den...) + return acos(clamp(num / den, -1, 1)) +end + +""" + spherical_distance(p::SphericalPoint, q::SphericalPoint) -> Number + +Computes the spherical distance between the points `p` and `q`. +""" +function spherical_distance(p, q) + # https://brsr.github.io/2021/05/01/vector-spherical-geometry.html + _px, _py, _pz = __getxyz(p) + _qx, _qy, _qz = __getxyz(q) + dx, dy, dz = _px - _qx, _py - _qy, _pz - _qz + ex, ey, ez = _px + _qx, _py + _qy, _pz + _qz + num = _hypot(dx, dy, dz) + den = _hypot(ex, ey, ez) + return 2atan(num / den) +end + +""" + spherical_angles(p::SphericalPoint, q::SphericalPoint, r::SphericalPoint) -> NTuple{3, Number} + +Computes the spherical angles of the spherical triangle defined by the points `p`, `q`, and `r`, returned +as a tuple `(α, β, γ)` for the interior angles at `p`, `q`, and `r`, respectively. + +See https://brsr.github.io/2021/05/01/vector-spherical-geometry.html. +""" +function spherical_angles(p, q, r) + return ( + _spherical_angle(p, q, r), + _spherical_angle(q, r, p), + _spherical_angle(r, p, q) + ) +end + +function _spherical_angle(p, q, r) + # Computes the interior angle at p + y = scalar_triple_product(p, q, r) + dot_qr = dot(q, r) + dot_pq = dot(p, q) + dot_rp = dot(r, p) + x = dot_qr - dot_pq * dot_rp + return atan(y, x) +end + +""" + spherical_triangle_centroid(p::SphericalPoint, q::SphericalPoint, r::SphericalPoint) -> SphericalPoint + +Computes the centroid of the spherical triangle defined by the points `p`, `q`, and `r`. +This is the point that separates the triangle into three equal-area regions. + +See https://brsr.github.io/2023/11/04/spherical-areal.html. +""" +function spherical_triangle_centroid(p, q, r) + A, B, C = spherical_angles(p, q, r) + S = (A + B + C) / 2 + Ω = (2S - π) / 3 + x = csc(A - Ω) + y = csc(B - Ω) + z = csc(C - Ω) + _px, _py, _pz = __getxyz(p) + _qx, _qy, _qz = __getxyz(q) + _rx, _ry, _rz = __getxyz(r) + sa = _hypot(cross(q, r)...) + sb = _hypot(cross(r, p)...) + sc = _hypot(cross(p, q)...) + Px = _px * sa * x + _qx * sb * y + _rx * sc * z + Py = _py * sa * x + _qy * sb * y + _ry * sc * z + Pz = _pz * sa * x + _qz * sb * y + _rz * sc * z + scale = _hypot(Px, Py, Pz) + return SphericalPoint(Px / scale, Py / scale, Pz / scale) +end \ No newline at end of file diff --git a/src/spherical/triangulation.jl b/src/spherical/triangulation.jl new file mode 100644 index 000000000..2252b5b44 --- /dev/null +++ b/src/spherical/triangulation.jl @@ -0,0 +1,117 @@ +""" + SphericalTriangulation = Triangulation{<:SphericalPoints} + +A type alias for a triangulation of points on the surface of a [`UnitSphere`](@ref). + +For this triangulation, note that while ghost triangles are present, they should be interpreted as solid triangles - +they are the triangles that help patch up the north pole. + +The ghost vertex of the triangulation maps to `(NaN, NaN)`, which is the north pole (after `invproject`). + +If you want to replace all these triangles with solid triangles, use [`solidify!`](@ref). +""" +const SphericalTriangulation = Triangulation{<:SphericalPoints} + +function triangle_circumcenter(tri::SphericalTriangulation, T, _=nothing) # the unused argument is the triangle area + points = get_points(tri) + u, v, w = triangle_vertices(T) + p, q, r = _getindex(points, u), _getindex(points, v), _getindex(points, w) # can't do get_point since DelaunayTriangulation projects it + sc = spherical_circumcenter(p, q, r) + return project(sc) # Voronoi vertices are assumed to be NTuple{2} +end + +function triangle_area(tri::SphericalTriangulation, T) + T = sort_triangle(T) + points = get_points(tri) + u, v, w = triangle_vertices(T) + p, q, r = _getindex(points, u), _getindex(points, v), _getindex(points, w) # can't do get_point since DelaunayTriangulation projects it + sc = spherical_triangle_area(p, q, r) + return sc +end + +function triangle_circumradius(tri::SphericalTriangulation, T, _=nothing) # the unused argument is the triangle area + points = get_points(tri) + u, v, w = triangle_vertices(T) + p, q, r = _getindex(points, u), _getindex(points, v), _getindex(points, w) # can't do get_point since DelaunayTriangulation projects it + sc = spherical_triangle_circumradius(p, q, r) + return sc +end + +function triangle_centroid(tri::SphericalTriangulation, T) + points = get_points(tri) + u, v, w = triangle_vertices(T) + p, q, r = _getindex(points, u), _getindex(points, v), _getindex(points, w) # can't do get_point since DelaunayTriangulation projects it + sc = spherical_triangle_centroid(p, q, r) + return project(sc) +end + +function _check_args!(kwargs) + if haskey(kwargs, :weights) + throw(ArgumentError("Spherical triangulations do not support weighted points.")) + end + if haskey(kwargs, :boundary_nodes) + throw(ArgumentError("Spherical triangulations do not support boundary nodes.")) + end + if haskey(kwargs, :segments) + throw(ArgumentError("Spherical triangulations do not support segments.")) + end + delete!(kwargs, :delete_ghosts) + delete!(kwargs, :recompute_representative_points) + return kwargs +end + +""" + spherical_triangulate(points::AbstractVector{<:SphericalPoint}; kwargs...) -> SphericalTriangulation + +Computes the spherical Delaunay triangulation of the points `points` on the unit sphere. The keyword arguments are the same as +those for `triangulate`, except that `weights`, `boundary_nodes`, and `segments` are not supported. + +See the docstring for [`SphericalTriangulation`](@ref) for more information. +""" +function spherical_triangulate(points::AbstractVector{<:SphericalPoint}; rng::AbstractRNG=Random.default_rng(), predicates::AbstractPredicateKernel=AdaptiveKernel(), kwargs...) + wpoints = SphericalPoints(points) + kwargs = _check_args!(Dict(kwargs)) + T = number_type(wpoints) + np = north_pole(T) + idx = findfirst(==(np), each_point(wpoints)) + _skip_points = isnothing(idx) ? () : idx + skip_points = haskey(kwargs, :skip_points) ? union(_skip_points, pop!(kwargs, :skip_points)) : _skip_points + tri = triangulate(wpoints; skip_points, kwargs..., delete_ghosts=false, recompute_representative_points=false, rng, predicates) + # lock_convex_hull!(tri; rng, predicates) + + empty_representative_points!(tri) + I = integer_type(tri) + rep = new_representative_point!(tri, one(I)) + rep[1] = RepresentativeCoordinates(T(NaN), T(NaN), zero(I)) + return tri::SphericalTriangulation +end + +""" + solidify!(tri::SphericalTriangulation) + +Replaces all the ghost triangles representing triangles covering the patch around the north pole with +solid triangles. The north pole is added as a point to the triangulation. Note that, after applying this +function, the projected triangulation will no longer be planar as the new solid triangles will intersect +other triangles. + +It is recommended that you just handle the ghost triangles directly instead of using this function. +""" +function solidify!(tri::SphericalTriangulation) + T = number_type(tri) + np = north_pole(T) + idx = findfirst(==(np), each_point(tri)) + if isnothing(idx) + push_point!(tri, np) + n = num_points(tri) + else + n = idx + end + for T in each_ghost_triangle(tri) + u, v, w = triangle_vertices(sort_triangle(T)) + u′, v′, w′ = sort_triangle(u, v, w) + add_triangle!(tri, v′, u′, n) + end + delete_ghost_triangles!(tri) + delete_ghost_vertices_from_graph!(tri) + return tri +end \ No newline at end of file diff --git a/src/spherical/utils.jl b/src/spherical/utils.jl new file mode 100644 index 000000000..374a1f908 --- /dev/null +++ b/src/spherical/utils.jl @@ -0,0 +1,48 @@ +""" + cross(p::SphericalPoint, q::SphericalPoint) -> NTuple{3, Number} + +Computes the cross product of two points on the sphere, treating them +as vectors from the origin to the points. The result is returned +as an `NTuple{3, Number}`. +""" +function cross(p::SphericalPoint, q::SphericalPoint) + _px, _py, _pz = __getxyz(p) + _qx, _qy, _qz = __getxyz(q) + return ( + _py * _qz - _pz * _qy, + _pz * _qx - _px * _qz, + _px * _qy - _py * _qx + ) +end + +_hypot(x, y) = hypot(x, y) +_hypot(x, y, z) = hypot(hypot(x, y), z) + +function dot(p::SphericalPoint, q::SphericalPoint) + _px, _py, _pz = __getxyz(p) + _qx, _qy, _qz = __getxyz(q) + return _px * _qx + _py * _qy + _pz * _qz +end + +function scalar_triple_product(p::SphericalPoint, q::SphericalPoint, r::SphericalPoint) + return dot(p, SphericalPoint(cross(q, r))) +end + +function perturbed_north_pole(::Type{T}) where {T} + dx, dy = 1e-8, 1e-8 + return SphericalPoint(T(dx), T(dy), T(sqrt(1 - dx^2 - dy^2))) +end + +function fix_north_pole(p, q, r) + T = number_type(p) + np = north_pole(T) + if p == np + return perturbed_north_pole(T), q, r + elseif q == np + return p, perturbed_north_pole(T), r + elseif r == np + return p, q, perturbed_north_pole(T) + else + return p, q, r + end +end \ No newline at end of file diff --git a/src/spherical/voronoi.jl b/src/spherical/voronoi.jl new file mode 100644 index 000000000..56fe27def --- /dev/null +++ b/src/spherical/voronoi.jl @@ -0,0 +1,32 @@ +const SphericalTessellation = VoronoiTessellation{<:SphericalTriangulation} + +function spherical_voronoi(tri::SphericalTriangulation; kwargs...) + vorn = voronoi(tri; kwargs...) + circumcenter_to_triangle = get_circumcenter_to_triangle(vorn) + triangle_to_circumcenter = get_triangle_to_circumcenter(vorn) + polygon_points = get_polygon_points(vorn) + I = integer_type(tri) + mapped_ghosts = Dict{I, I}() + for polygon_vertices in each_polygon(vorn) + for (i, v) in pairs(polygon_vertices) + is_ghost_vertex(v) || continue + if haskey(mapped_ghosts, v) + polygon_vertices[i] = mapped_ghosts[v] + else + T = get_circumcenter_to_triangle(vorn, v) + c = triangle_circumcenter(tri, T) + pidx = findfirst(==(c), polygon_points) + if isnothing(pidx) + push_polygon_point!(vorn, c) + n = num_polygon_vertices(vorn) + else + n = pidx + end + triangle_to_circumcenter[T] = n + mapped_ghosts[v] = n + polygon_vertices[i] = n + end + end + end + return vorn +end diff --git a/src/utils/utils.jl b/src/utils/utils.jl index 73f16c960..3bf0165d2 100644 --- a/src/utils/utils.jl +++ b/src/utils/utils.jl @@ -293,9 +293,7 @@ function get_area(tri::Triangulation) F = number_type(tri) A = zero(F) for T in each_solid_triangle(tri) - u, v, w = triangle_vertices(T) - p, q, r = get_point(tri, u, v, w) - A += triangle_area(p, q, r) + A += triangle_area(tri, T) end return A end diff --git a/src/validation.jl b/src/validation.jl index 560ebed37..57e7c935b 100644 --- a/src/validation.jl +++ b/src/validation.jl @@ -81,9 +81,8 @@ function test_delaunay_criterion(tri; predicates::AbstractPredicateKernel = Exac failures = Tuple{triangle_type(tri), integer_type(tri)}[] for T in each_solid_triangle(tri) i, j, k = triangle_vertices(T) - p, q, r = get_point(tri, i, j, k) - c = triangle_circumcenter(p, q, r) - cr = triangle_circumradius(p, q, r) + c = triangle_circumcenter(tri, T) + cr = triangle_circumradius(tri, T) xmin, xmax = getx(c) - cr, getx(c) + cr ymin, ymax = gety(c) - cr, gety(c) + cr any(!isfinite, (xmin, xmax, ymin, ymax)) && continue @@ -113,8 +112,8 @@ function test_delaunay_criterion(tri; predicates::AbstractPredicateKernel = Exac any(==(r), (i, j, k)) && continue T = construct_triangle(triangle_type(tri), i, j, k) cert = point_position_relative_to_circumcircle(predicates, tri, T, r) - c = triangle_centroid(get_point(tri, i, j, k)...) - A = triangle_area(get_point(tri, i, j, k)...) + c = triangle_centroid(tri, T) + A = triangle_area(tri, T) check_precision(A) && continue # the centroids in this case sometimes appear outside of the triangle if is_inside(cert) is_boundary_edge(tri, i, j) && is_right(point_position_relative_to_line(predicates, tri, i, j, r)) && continue # if it's outside of the domain relative to this edge, just continue diff --git a/test/Project.toml b/test/Project.toml index 7a113492b..63da1650d 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -15,9 +15,11 @@ ElasticArrays = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4" ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110" FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" diff --git a/test/runtests.jl b/test/runtests.jl index ba499c3ee..a3192a7a9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -120,6 +120,10 @@ end safe_include("voronoi/power.jl") end + @testset verbose = true "Spherical" begin + safe_include("spherical.jl") + end + @testset verbose = true "Run the documentation examples" begin @testset verbose = true "Check that the applications in the docs run" begin app_dir = joinpath(dirname(dirname(pathof(DelaunayTriangulation))), "docs", "src", "literate_applications") diff --git a/test/spherical.jl b/test/spherical.jl new file mode 100644 index 000000000..c833c6fa3 --- /dev/null +++ b/test/spherical.jl @@ -0,0 +1,207 @@ +using DelaunayTriangulation +using DelaunayTriangulation.SphericalDelaunay +using Test +using Makie, CairoMakie, GLMakie +CairoMakie.activate!() +using Random +using LinearAlgebra + +const DT = DelaunayTriangulation +const SD = SphericalDelaunay + +Makie.to_vertices(p::Vector{<:SphericalPoint}) = Makie.Point3.(SD.__getxyz.(p)) +Makie.convert_arguments(p::PointBased, points::Vector{<:SphericalPoint}) = Makie.convert_arguments(p, Makie.to_vertices(points)) +Makie.convert_arguments(p::PointBased, L::SphericalLine) = Makie.convert_arguments(p, SD._to_coords(L)) +Makie.convert_arguments(p::PointBased, T::SphericalTriangle) = Makie.convert_arguments(p, SD._to_coords(T)) +Makie.convert_arguments(p::PointBased, P::SphericalPolygon) = Makie.convert_arguments(p, SD._to_coords(P)) +function Makie.convert_arguments(p::Type{<:Lines}, tri::SphericalTriangulation) + points = SD._to_lines(tri) + return Makie.convert_arguments(p, points) +end +function Makie.convert_arguments(p::Type{<:Lines}, vorn::SphericalTessellation) + points = SD._to_lines(vorn) + return Makie.convert_arguments(p, points) +end +# These mesh recipes render SO poorly... +function Makie.convert_arguments(p::Type{<:Mesh}, T::SphericalTriangle) + points = SD._to_coords(T, false) + faces = SD._to_mesh(points) + return Makie.convert_arguments(p, points, faces) +end +function Makie.convert_arguments(p::Type{<:Mesh}, tri::SphericalTriangulation) + points, vertices = SD._to_triangles(tri) + return Makie.convert_arguments(p, points, vertices) +end +function Makie.convert_arguments(p::Type{<:Poly}, tri::SphericalTriangulation) + res = Makie.convert_arguments(Mesh, tri) + return Makie.convert_arguments(p, res...) +end + +function uniform_rand_sphere(N) + # http://web.archive.org/web/20120421191837/http://www.cgafaq.info/wiki/Evenly_distributed_points_on_sphere + coords = Vector{SphericalPoint{NTuple{3,Float64}}}(undef, N) + s = 3.6 / sqrt(N) + dz = 2 / N + z = 1 - dz / 2 + long = 0.0 + for k in 1:N + r = sqrt(1 - z^2) + latitude = asind(z) + longitude = rad2deg(long) + coords[k] = SphericalPoint(latitude, longitude % 360) + long += s / r + z -= dz + end + return coords +end + +@testset "Basic" begin + @test SD.UnitSphere() == SD.UnitSphere{Float64}() + p = (1.0, 0.0, 0.0) + @test SphericalPoint(p) == SphericalPoint{typeof(p)}(p) + @test SphericalPoint(p) == SphericalPoint(p) + @test SphericalPoint(0.3, 0.2, sqrt(1 - 0.3^2 - 0.2^2)) == SphericalPoint((0.3, 0.2, sqrt(1 - 0.3^2 - 0.2^2))) + q = (0.0, 1.0, 0.0) + @test SphericalPoint(p) != SphericalPoint(q) + p = rand(UnitSphere()) + @test p isa SphericalPoint + @test sum(p.p .^ 2) ≈ 1.0 + @test SD.getp(p) == p.p + @test SD.__getxyz(p) == p.p + x, y, z = p.p + @test SD.projectx(p) == getx(p) == x / (1 - z) + @test SD.projecty(p) == gety(p) == y / (1 - z) + @test SD.project(p) == getxy(p) + X, Y = SD.project(p) + @test collect(SD.invproject((X, Y)).p) ≈ collect(p.p) + @test DT.number_type(p) == Float64 + @test DT.number_type(rand(UnitSphere{Float32}())) == Float32 +end + +@testset "SphericalPoints" begin + points = rand(UnitSphere(), 25) + spoints = SD.SphericalPoints(points) + @test SD.SphericalPoints(spoints) == spoints + @test eltype(points) == SphericalPoint{NTuple{3,Float64}} + @test eachindex(spoints) == eachindex(points) + @test iterate(spoints) == iterate(points) + @test collect(spoints) == collect(points) + @test length(spoints) == length(points) + @test size(spoints) == size(points) + @test points[2] == spoints[2] + @test DT.number_type(spoints) == Float64 + DT.set_point!(spoints, 7, SD.project(SphericalPoint(0.3, 0.2, sqrt(1 - 0.3^2 - 0.2^2)))...) + @test collect(spoints[7].p) ≈ collect(SphericalPoint(0.3, 0.2, sqrt(1 - 0.3^2 - 0.2^2)).p) + val2 = spoints[end] + val = pop!(spoints) + @test val == val2 + @test length(spoints) == length(points) == 24 + DT.push_point!(spoints, SD.project(SphericalPoint(0.3, 0.2, sqrt(1 - 0.3^2 - 0.2^2)))...) + @test collect(spoints[end].p) ≈ collect(SphericalPoint(0.3, 0.2, sqrt(1 - 0.3^2 - 0.2^2)).p) + @test SD._getindex(spoints, -1) == SD.north_pole(Float64) +end + +function spherical_distance(p, q) + a1, a2, a3 = p + b1, b2, b3 = q + return acos(a1 * b1 + a2 * b2 + a3 * b3) +end +_cross(p, q) = cross(collect(p), collect(q)) +@testset "Spherical Triangles" begin + for _ in 1:10000 + p, q = rand(UnitSphere(), 2) + c = cross(collect(p.p), collect(q.p)) + @test collect(c) ≈ collect(SD.cross(p, q)) + p, q, r = rand(UnitSphere(), 3) + sc = SD.spherical_circumcenter(p, q, r) + dists = spherical_distance.((sc.p,), (p.p, q.p, r.p)) + @test all(dists .≈ dists[1]) + p, q = rand(UnitSphere(), 2) + d = dot(collect(p.p), collect(q.p)) + @test d ≈ SD.dot(p, q) + u, v, w = rand(3) + @test SD._hypot(u, v, w) ≈ sqrt(u^2 + v^2 + w^2) + p, q, r = SphericalPoint(30.2672, 97.7431), SphericalPoint(42.3584, 71.0598), SphericalPoint(51.0501, 114.0853) + A = SD.spherical_triangle_area(p, q, r) + @test A ≈ 0.090030685 + q, r = SphericalPoint(29.7947, 98.7320), SphericalPoint(32.9746, 96.8899) + A = SD.spherical_triangle_area(p, q, r) + @test A ≈ 0.0003032270459 + p, q, r = rand(UnitSphere(), 3) + s = SD.spherical_triangle_circumradius(p, q, r) + sc = SD.spherical_circumcenter(p, q, r) + dists = spherical_distance.((sc.p,), (p.p, q.p, r.p)) + @test all(dists .≈ s) + p, q = rand(UnitSphere(), 2) + @test SD.spherical_distance(p, q) ≈ spherical_distance(p.p, q.p) + p, q, r = rand(UnitSphere(), 3) + stp = SD.scalar_triple_product(p, q, r) + @test stp ≈ dot(collect(p.p), cross(collect(q.p), collect(r.p))) + + p, q, r = rand(UnitSphere(), 3) + if SD.scalar_triple_product(p, q, r) < 0 + p, q, r = q, p, r + end + A, B, C = SD.spherical_angles(p, q, r) + @test sin(A) / norm(cross(collect(q.p), collect(r.p))) ≈ + sin(B) / norm(cross(collect(r.p), collect(p.p))) ≈ + sin(C) / norm(cross(collect(p.p), collect(q.p))) ≈ + SD.scalar_triple_product(p, q, r) / (norm(cross(collect(p.p), collect(q.p))) * norm(cross(collect(q.p), collect(r.p))) * norm(cross(collect(r.p), collect(p.p)))) + excess = A + B + C - π + @test SD.spherical_triangle_area(p, q, r) ≈ excess + + p, q, r = rand(UnitSphere(), 3) + if SD.scalar_triple_product(p, q, r) < 0 + p, q, r = q, p, r + end + c = SD.spherical_triangle_centroid(p, q, r) + @test norm(collect(c.p)) ≈ 1.0 + A1 = SD.spherical_triangle_area(p, q, c) + A2 = SD.spherical_triangle_area(q, r, c) + A3 = SD.spherical_triangle_area(r, p, c) + @test A1 ≈ A2 ≈ A3 ≈ SD.spherical_triangle_area(p, q, r) / 3 + end + p, q, r = rand(UnitSphere(), 3) + @inferred SD.spherical_triangle_circumradius(p, q, r) + @inferred SD.spherical_triangle_area(p, q, r) + @inferred SD.spherical_circumcenter(p, q, r) + @inferred SD.cross(p, q) + @inferred SD.dot(p, q) + @inferred SD._hypot(rand(3)...) + @inferred SD.projectx(p) + @inferred SD.projecty(p) + @inferred SD.project(p) + X, Y = SD.project(rand(UnitSphere())) + @inferred SD.invproject((X, Y)) + @inferred SD.getp(p) + @inferred SD.__getxyz(p) + @test SD.north_pole(Float64) == SphericalPoint(0.0, 0.0, 1.0) + @test SD.north_pole(Float32) == SphericalPoint(0.0, 0.0, 1.0f0) + @test SD.invproject((NaN, NaN)) == SphericalPoint(0.0, 0.0, 1.0) +end + +@testset "SphericalLine/Triangle" begin + for _ in 1:100 + p, q = rand(UnitSphere(), 2) + L = SD.SphericalLine(p, q) + t = LinRange(0.0, L.upper, 500) + pts = L.(t) + _pts = [collect(x.p) for x in pts] + dist = sum(norm.(diff(_pts))) + @test dist ≈ spherical_distance(p.p, q.p) rtol = 1e-4 + + p, q, r = rand(UnitSphere(), 3) + L1, L2, L3 = SphericalLine(p, q), SphericalLine(q, r), SphericalLine(r, p) + T = SD.SphericalTriangle(L1, L2, L3) + T2 = SD.SphericalTriangle(p, q, r) + @test T == T2 + end +end + +points = uniform_rand_sphere(15) +tri = spherical_triangulate(points) +refine!(tri, max_area = 1e-1) + +points = uniform_rand_sphere(15) +tri = spherical_triangulate(points) +vorn = spherical_voronoi(tri) From 5a973ca1d63e7edcc613da570a7ff01c77a970d1 Mon Sep 17 00:00:00 2001 From: Daniel VandenHeuvel <95613936+DanielVandH@users.noreply.github.com> Date: Sun, 22 Sep 2024 12:57:59 +0100 Subject: [PATCH 2/6] oops --- src/spherical/refine.jl | 2 +- test/spherical.jl | 9 ++------- 2 files changed, 3 insertions(+), 8 deletions(-) diff --git a/src/spherical/refine.jl b/src/spherical/refine.jl index 78a195776..1ed514fa0 100644 --- a/src/spherical/refine.jl +++ b/src/spherical/refine.jl @@ -77,7 +77,7 @@ function enqueue_triangle!(tri::SphericalTriangulation, queue::RefinementQueue, return setindex!(queue, A, T) end -function finalise!(tri::Triangulation, args::RefinementArguments) +function finalise!(tri::SphericalTriangulation, args::RefinementArguments) #unlock_convex_hull!(tri; reconstruct = true) #lock_convex_hull!(tri; rng = args.rng, predicates = args.predicates) convex_hull!(tri; reconstruct = false, predicates = args.predicates) diff --git a/test/spherical.jl b/test/spherical.jl index c833c6fa3..21571cc3e 100644 --- a/test/spherical.jl +++ b/test/spherical.jl @@ -198,10 +198,5 @@ end end end -points = uniform_rand_sphere(15) -tri = spherical_triangulate(points) -refine!(tri, max_area = 1e-1) - -points = uniform_rand_sphere(15) -tri = spherical_triangulate(points) -vorn = spherical_voronoi(tri) +points = uniform_rand_sphere(35) +tri = spherical_triangulate(points) \ No newline at end of file From 6df7c800d7c62d0539be93c29c2ce48d436d06a0 Mon Sep 17 00:00:00 2001 From: Daniel VandenHeuvel <95613936+DanielVandH@users.noreply.github.com> Date: Sun, 22 Sep 2024 13:26:12 +0100 Subject: [PATCH 3/6] rm gl --- test/Project.toml | 1 - test/spherical.jl | 16 ++++++++++++---- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/test/Project.toml b/test/Project.toml index 63da1650d..fc8457ba3 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -15,7 +15,6 @@ ElasticArrays = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4" ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110" FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" -GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/test/spherical.jl b/test/spherical.jl index 21571cc3e..d6c70c018 100644 --- a/test/spherical.jl +++ b/test/spherical.jl @@ -1,7 +1,7 @@ using DelaunayTriangulation using DelaunayTriangulation.SphericalDelaunay using Test -using Makie, CairoMakie, GLMakie +using Makie, CairoMakie#, GLMakie CairoMakie.activate!() using Random using LinearAlgebra @@ -37,7 +37,7 @@ function Makie.convert_arguments(p::Type{<:Poly}, tri::SphericalTriangulation) return Makie.convert_arguments(p, res...) end -function uniform_rand_sphere(N) +function uniform_sphere(N) # http://web.archive.org/web/20120421191837/http://www.cgafaq.info/wiki/Evenly_distributed_points_on_sphere coords = Vector{SphericalPoint{NTuple{3,Float64}}}(undef, N) s = 3.6 / sqrt(N) @@ -198,5 +198,13 @@ end end end -points = uniform_rand_sphere(35) -tri = spherical_triangulate(points) \ No newline at end of file +#= +tri1 = spherical_triangulate(uniform_sphere(50)) +tri2 = spherical_triangulate(uniform_sphere(50)) +refine!(tri2; max_area = 0.004pi) +tri3 = spherical_triangulate(uniform_sphere(100)) +vorn3 = spherical_voronoi(tri3) +display(GLMakie.Screen(), lines(tri1)) +display(GLMakie.Screen(), lines(tri2)) +display(GLMakie.Screen(), lines(vorn3)) +=# \ No newline at end of file From b26276e87724f34438c79fba3017a815624f837d Mon Sep 17 00:00:00 2001 From: Daniel VandenHeuvel <95613936+DanielVandH@users.noreply.github.com> Date: Mon, 23 Sep 2024 10:43:09 +0100 Subject: [PATCH 4/6] address point location failures partially --- src/spherical/triangulation.jl | 10 +++++++--- test/Project.toml | 1 + test/spherical.jl | 31 +++++++++++++++++++++++++++++++ 3 files changed, 39 insertions(+), 3 deletions(-) diff --git a/src/spherical/triangulation.jl b/src/spherical/triangulation.jl index 2252b5b44..6b2cb54c8 100644 --- a/src/spherical/triangulation.jl +++ b/src/spherical/triangulation.jl @@ -6,12 +6,15 @@ A type alias for a triangulation of points on the surface of a [`UnitSphere`](@r For this triangulation, note that while ghost triangles are present, they should be interpreted as solid triangles - they are the triangles that help patch up the north pole. -The ghost vertex of the triangulation maps to `(NaN, NaN)`, which is the north pole (after `invproject`). +Note that the ghost vertex does not map to the north pole due to technical reasons. Use `north_pole` to obtain +the north pole. If you want to replace all these triangles with solid triangles, use [`solidify!`](@ref). """ const SphericalTriangulation = Triangulation{<:SphericalPoints} +north_pole(tri::SphericalTriangulation) = north_pole(number_type(tri)) + function triangle_circumcenter(tri::SphericalTriangulation, T, _=nothing) # the unused argument is the triangle area points = get_points(tri) u, v, w = triangle_vertices(T) @@ -79,10 +82,12 @@ function spherical_triangulate(points::AbstractVector{<:SphericalPoint}; rng::Ab tri = triangulate(wpoints; skip_points, kwargs..., delete_ghosts=false, recompute_representative_points=false, rng, predicates) # lock_convex_hull!(tri; rng, predicates) + #= empty_representative_points!(tri) I = integer_type(tri) rep = new_representative_point!(tri, one(I)) rep[1] = RepresentativeCoordinates(T(NaN), T(NaN), zero(I)) + =# # This causes issues with point location :-( return tri::SphericalTriangulation end @@ -97,8 +102,7 @@ other triangles. It is recommended that you just handle the ghost triangles directly instead of using this function. """ function solidify!(tri::SphericalTriangulation) - T = number_type(tri) - np = north_pole(T) + np = north_pole(tri) idx = findfirst(==(np), each_point(tri)) if isnothing(idx) push_point!(tri, np) diff --git a/test/Project.toml b/test/Project.toml index fc8457ba3..74120d75b 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -6,6 +6,7 @@ CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" +DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" diff --git a/test/spherical.jl b/test/spherical.jl index d6c70c018..10415a4be 100644 --- a/test/spherical.jl +++ b/test/spherical.jl @@ -198,6 +198,37 @@ end end end +function slow_find_triangle(tri, p) + points = get_points(tri) + for T in each_triangle(tri) + A = DT.triangle_area(tri, T) + u, v, w = triangle_vertices(T) + a, b, c = SD._getindex(points, u), SD._getindex(points, v), SD._getindex(points, w) + A1 = SD.spherical_triangle_area(a, b, p) + A2 = SD.spherical_triangle_area(b, c, p) + A3 = SD.spherical_triangle_area(c, a, p) + if A ≈ A1 + A2 + A3 + return T + end + end + throw(ErrorException("Triangle not found")) +end + +points = uniform_sphere(50) +tri = spherical_triangulate(points) +p = SphericalPoint(0.3908413014578491, 0.28203954256710945, 0.8761830707696139) +Random.seed!(0) +find_triangle(tri, p) +slow_find_triangle(tri, p) + + +p = rand(UnitSphere(), 25000) +for p in p + @show p + T = slow_find_triangle(tri, p) + @test DT.sort_triangle(T) == DT.sort_triangle(find_triangle(tri, p)) +end + #= tri1 = spherical_triangulate(uniform_sphere(50)) tri2 = spherical_triangulate(uniform_sphere(50)) From 18daa1ada228a175cf59a2286721b90157f521d7 Mon Sep 17 00:00:00 2001 From: Daniel VandenHeuvel <95613936+DanielVandH@users.noreply.github.com> Date: Mon, 23 Sep 2024 12:12:13 +0100 Subject: [PATCH 5/6] Fix point location --- .../point_location/jump_and_march.jl | 1 + .../basic_operations/add_point.jl | 2 + .../constrained_triangulation.jl | 2 +- .../triangulation/mesh_refinement.jl | 2 +- .../unconstrained_triangulation.jl | 2 +- src/predicates/predicates.jl | 2 +- src/spherical/SphericalDelaunay.jl | 15 ++-- src/spherical/spherical_triangles.jl | 32 +++++++- src/spherical/triangulation.jl | 61 +++++++++++++- test/Project.toml | 1 - test/spherical.jl | 82 ++++++++++--------- 11 files changed, 151 insertions(+), 51 deletions(-) diff --git a/src/algorithms/point_location/jump_and_march.jl b/src/algorithms/point_location/jump_and_march.jl index d250080d3..dc3d38a1d 100644 --- a/src/algorithms/point_location/jump_and_march.jl +++ b/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) diff --git a/src/algorithms/triangulation/basic_operations/add_point.jl b/src/algorithms/triangulation/basic_operations/add_point.jl index b12cf9327..8c8769c69 100644 --- a/src/algorithms/triangulation/basic_operations/add_point.jl +++ b/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} diff --git a/src/algorithms/triangulation/constrained_triangulation.jl b/src/algorithms/triangulation/constrained_triangulation.jl index 4cc653888..63ce1d957 100644 --- a/src/algorithms/triangulation/constrained_triangulation.jl +++ b/src/algorithms/triangulation/constrained_triangulation.jl @@ -711,7 +711,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)) diff --git a/src/algorithms/triangulation/mesh_refinement.jl b/src/algorithms/triangulation/mesh_refinement.jl index 41cc9969e..a7d974f0c 100644 --- a/src/algorithms/triangulation/mesh_refinement.jl +++ b/src/algorithms/triangulation/mesh_refinement.jl @@ -260,7 +260,7 @@ 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(tri, V) && is_on(flag) V = replace_ghost_triangle_with_boundary_triangle(tri, V) diff --git a/src/algorithms/triangulation/unconstrained_triangulation.jl b/src/algorithms/triangulation/unconstrained_triangulation.jl index ec55b500e..de852e275 100644 --- a/src/algorithms/triangulation/unconstrained_triangulation.jl +++ b/src/algorithms/triangulation/unconstrained_triangulation.jl @@ -150,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 diff --git a/src/predicates/predicates.jl b/src/predicates/predicates.jl index 7c7b497cf..9673063da 100644 --- a/src/predicates/predicates.jl +++ b/src/predicates/predicates.jl @@ -798,7 +798,7 @@ the two does not intersect any segments. - `cert`: A [`Certificate`](@ref). This will be `Visible` if `i` is visible from `q`, and `Invisible` otherwise. """ function test_visibility(kernel::AbstractPredicateKernel, tri::Triangulation, q, i) - V, invisible_flag = find_triangle(tri, q; use_barriers = Val(true), k = i, concavity_protection = true, predicates = kernel) + V, invisible_flag = find_triangle(tri, q; use_barriers = Val(true), k = i, concavity_protection = true, predicates = kernel, check_sphere = false) if invisible_flag return Certificate.Invisible else diff --git a/src/spherical/SphericalDelaunay.jl b/src/spherical/SphericalDelaunay.jl index 57b7a195d..1a4f0a419 100644 --- a/src/spherical/SphericalDelaunay.jl +++ b/src/spherical/SphericalDelaunay.jl @@ -28,13 +28,13 @@ import ..DelaunayTriangulation: set_point!, push_point!, convex_hull!, + each_edge, + get_adjacent2vertex, AbstractParametricCurve, norm_sqr, - get_representative_point_list, - RepresentativeCoordinates, + get_adjacent, + is_outside, integer_type, - empty_representative_points!, - new_representative_point!, delete_ghost_triangles!, num_points, delete_ghost_vertices_from_graph!, @@ -43,8 +43,10 @@ import ..DelaunayTriangulation: each_solid_triangle, num_solid_triangles, num_triangles, + PointNotFoundError, each_triangle, is_planar, + is_inside, _quality_statistics, RefinementArguments, triangle_radius_edge_ratio, @@ -61,6 +63,7 @@ import ..DelaunayTriangulation: _build_queues, check_steiner_point_precision, finalise!, + construct_triangle, _each_solid_triangle, has_ghost_triangles, _is_ghost_triangle, @@ -78,7 +81,9 @@ import ..DelaunayTriangulation: each_polygon, each_polygon_index, get_polygon, - get_polygon_point + get_polygon_point, + find_triangle, + triangle_type import Random: Random, diff --git a/src/spherical/spherical_triangles.jl b/src/spherical/spherical_triangles.jl index 37dfff02e..711f2d25a 100644 --- a/src/spherical/spherical_triangles.jl +++ b/src/spherical/spherical_triangles.jl @@ -45,7 +45,7 @@ end Computes the spherical distance between the points `p` and `q`. """ -function spherical_distance(p, q) +function spherical_distance(p::SphericalPoint, q::SphericalPoint) # https://brsr.github.io/2021/05/01/vector-spherical-geometry.html _px, _py, _pz = __getxyz(p) _qx, _qy, _qz = __getxyz(q) @@ -103,9 +103,37 @@ function spherical_triangle_centroid(p, q, r) sa = _hypot(cross(q, r)...) sb = _hypot(cross(r, p)...) sc = _hypot(cross(p, q)...) - Px = _px * sa * x + _qx * sb * y + _rx * sc * z + Px = _px * sa * x + _qx * sb * y + _rx * sc * z Py = _py * sa * x + _qy * sb * y + _ry * sc * z Pz = _pz * sa * x + _qz * sb * y + _rz * sc * z scale = _hypot(Px, Py, Pz) return SphericalPoint(Px / scale, Py / scale, Pz / scale) +end + +""" + point_position_relative_to_spherical_triangle(p::SphericalPoint, q::SphericalPoint, r::SphericalPoint, s::SphericalPoint) -> Certificate + +Tests if `s` is in the spherical triangle defined by the points `p`, `q`, and `r`. Returns a `Certificate` object, which is + +- `Certificate.Inside`: if `s` is in the triangle, +- `Certificate.On`: if `s` is on the edge of the triangle, or +- `Certificate.Outside`: if `s` is outside the triangle. +""" +function point_position_relative_to_spherical_triangle(p::SphericalPoint, q::SphericalPoint, r::SphericalPoint, s::SphericalPoint) + x1, y1, z1 = __getxyz(p) + x2, y2, z2 = __getxyz(q) + x3, y3, z3 = __getxyz(r) + x, y, z = __getxyz(s) + D = x * y1 * z2 - x * y2 * z1 - x1 * y * z2 + x1 * y2 * z + x2 * y * z1 - x2 * y1 * z - x * y1 * z3 + x * y3 * z1 + x1 * y * z3 - x1 * y3 * z - x3 * y * z1 + x3 * y1 * z + x * y2 * z3 - x * y3 * z2 - x2 * y * z3 + x2 * y3 * z + x3 * y * z2 - x3 * y2 * z + a = (x * y2 * z3 - x * y3 * z2 - x2 * y * z3 + x2 * y3 * z + x3 * y * z2 - x3 * y2 * z) / D + b = -(x * y1 * z3 - x * y3 * z1 - x1 * y * z3 + x1 * y3 * z + x3 * y * z1 - x3 * y1 * z) / D + c = (D + x * y1 * z3 - x * y3 * z1 - x1 * y * z3 + x1 * y3 * z + x3 * y * z1 - x3 * y1 * z - x * y2 * z3 + x * y3 * z2 + x2 * y * z3 - x2 * y3 * z - x3 * y * z2 + x3 * y2 * z) / D + lambda = (x1 * y2 * z3 - x1 * y3 * z2 - x2 * y1 * z3 + x2 * y3 * z1 + x3 * y1 * z2 - x3 * y2 * z1) / D + if a > 0 && b > 0 && c > 0 && lambda > 0 + return Certificate.Inside + elseif a == 0 || b == 0 || c == 0 + return Certificate.On + else + return Certificate.Outside + end end \ No newline at end of file diff --git a/src/spherical/triangulation.jl b/src/spherical/triangulation.jl index 6b2cb54c8..d234083d0 100644 --- a/src/spherical/triangulation.jl +++ b/src/spherical/triangulation.jl @@ -9,7 +9,10 @@ they are the triangles that help patch up the north pole. Note that the ghost vertex does not map to the north pole due to technical reasons. Use `north_pole` to obtain the north pole. -If you want to replace all these triangles with solid triangles, use [`solidify!`](@ref). +If you want to replace all these triangles with solid triangles, use [`solidify!`](@ref). This is not recommended. + +Note also that the triangulations, and e.g. point location, are not as robust as the planar case. You may +experience numerical issues near the north pole for example. """ const SphericalTriangulation = Triangulation{<:SphericalPoints} @@ -118,4 +121,60 @@ function solidify!(tri::SphericalTriangulation) delete_ghost_triangles!(tri) delete_ghost_vertices_from_graph!(tri) return tri +end + +function point_position_relative_to_spherical_triangle(tri::SphericalTriangulation, V, s) + u, v, w = triangle_vertices(V) + points = get_points(tri) + p, q, r = _getindex(points, u), _getindex(points, v), _getindex(points, w) + return point_position_relative_to_spherical_triangle(p, q, r, s) +end + +""" + find_triangle(tri::SphericalTriangulation, q; check_sphere=true, kwargs...) + +Finds the spherical triangle in `tri` that contains the point `q`. The keyword arguments are the same as those for +the usual [`find_triangle`], except that `check_sphere` can be used to check if the found triangle actually contains `q` +(since the check is done in the stereographic projection). If it's not, it will be found by searching the triangles near the found +triangle or using an exhaustive search. +""" +function find_triangle(tri::SphericalTriangulation, s; check_sphere = true, kwargs...) + V = @invoke find_triangle(tri::Triangulation, s; kwargs..., check_sphere = false, use_barriers = Val(false)) + check_sphere || return V + cert = point_position_relative_to_spherical_triangle(tri, V, s) + if !is_outside(cert) + return V + else + T = triangle_type(tri) + u, v, w = triangle_vertices(V) + ℓuv = get_adjacent(tri, v, u) + ℓvw = get_adjacent(tri, w, v) + ℓwu = get_adjacent(tri, u, w) + V1 = construct_triangle(T, v, u, ℓuv) + V2 = construct_triangle(T, w, v, ℓvw) + V3 = construct_triangle(T, u, w, ℓwu) + for _V in (V1, V2, V3) + cert = point_position_relative_to_spherical_triangle(tri, _V, s) + !is_outside(cert) && return _V + end + # Still not found. Find the nearest vertex and then look at its neighbours + points = get_points(tri) + p, q, r = _getindex(points, u), _getindex(points, v), _getindex(points, w) + d1 = spherical_distance(p, s) + d2 = spherical_distance(q, s) + d3 = spherical_distance(r, s) + d = min(d1, d2, d3) + minidx = d == d1 ? u : d == d2 ? v : w + for (i, j) in each_edge(get_adjacent2vertex(tri, minidx)) + _V = construct_triangle(T, i, j, minidx) + cert = point_position_relative_to_spherical_triangle(tri, _V, s) + !is_outside(cert) && return _V + end + # Still not found. Do an exhaustive search + for _V in each_triangle(tri) + cert = point_position_relative_to_spherical_triangle(tri, _V, s) + !is_outside(cert) && return _V + end + throw(PointNotFoundError(tri, s)) + end end \ No newline at end of file diff --git a/test/Project.toml b/test/Project.toml index 74120d75b..fc8457ba3 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -6,7 +6,6 @@ CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" -DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" diff --git a/test/spherical.jl b/test/spherical.jl index 10415a4be..c065aab34 100644 --- a/test/spherical.jl +++ b/test/spherical.jl @@ -76,6 +76,11 @@ end @test collect(SD.invproject((X, Y)).p) ≈ collect(p.p) @test DT.number_type(p) == Float64 @test DT.number_type(rand(UnitSphere{Float32}())) == Float32 + for _ in 1:100000 + p = rand(UnitSphere()) + X, Y = SD.project(p) + @test collect(SD.invproject((X, Y)).p) ≈ collect(p.p) + end end @testset "SphericalPoints" begin @@ -198,44 +203,45 @@ end end end -function slow_find_triangle(tri, p) - points = get_points(tri) - for T in each_triangle(tri) - A = DT.triangle_area(tri, T) - u, v, w = triangle_vertices(T) - a, b, c = SD._getindex(points, u), SD._getindex(points, v), SD._getindex(points, w) - A1 = SD.spherical_triangle_area(a, b, p) - A2 = SD.spherical_triangle_area(b, c, p) - A3 = SD.spherical_triangle_area(c, a, p) - if A ≈ A1 + A2 + A3 - return T +@testset "Point location" begin + function slow_find_triangle(tri, p) + points = get_points(tri) + for T in each_triangle(tri) + A = DT.triangle_area(tri, T) + u, v, w = triangle_vertices(T) + a, b, c = SD._getindex(points, u), SD._getindex(points, v), SD._getindex(points, w) + A1 = SD.spherical_triangle_area(a, b, p) + A2 = SD.spherical_triangle_area(b, c, p) + A3 = SD.spherical_triangle_area(c, a, p) + if A ≈ A1 + A2 + A3 + return T + end end + throw(ErrorException("Point not found")) end - throw(ErrorException("Triangle not found")) -end - -points = uniform_sphere(50) -tri = spherical_triangulate(points) -p = SphericalPoint(0.3908413014578491, 0.28203954256710945, 0.8761830707696139) -Random.seed!(0) -find_triangle(tri, p) -slow_find_triangle(tri, p) - - -p = rand(UnitSphere(), 25000) -for p in p - @show p - T = slow_find_triangle(tri, p) - @test DT.sort_triangle(T) == DT.sort_triangle(find_triangle(tri, p)) -end + points = uniform_sphere(50) + tri = spherical_triangulate(points) + for s in rand(UnitSphere(), 25000) + yield() + # point_position_relative_to_spherical_triangle + T = slow_find_triangle(tri, s) + T = DT.sort_triangle(T) + for V in each_triangle(tri) + V = DT.sort_triangle(V) + p, q, r = SD._getindex.((get_points(tri),), triangle_vertices(V)) + cert1 = SD.point_position_relative_to_spherical_triangle(p, q, r, s) + @inferred SD.point_position_relative_to_spherical_triangle(p, q, r, s) + cert2 = SD.point_position_relative_to_spherical_triangle(tri, V, s) + if T == V + @test DT.is_inside(cert1) && DT.is_inside(cert2) + else + @test DT.is_outside(cert1) && DT.is_outside(cert2) + end + end -#= -tri1 = spherical_triangulate(uniform_sphere(50)) -tri2 = spherical_triangulate(uniform_sphere(50)) -refine!(tri2; max_area = 0.004pi) -tri3 = spherical_triangulate(uniform_sphere(100)) -vorn3 = spherical_voronoi(tri3) -display(GLMakie.Screen(), lines(tri1)) -display(GLMakie.Screen(), lines(tri2)) -display(GLMakie.Screen(), lines(vorn3)) -=# \ No newline at end of file + # find_triangle + V = DT.sort_triangle(find_triangle(tri, s)) + @inferred DT.sort_triangle(find_triangle(tri, s)) + @test T == V + end +end \ No newline at end of file From 2a77be284de48001cb7faa7dffe241a869d4f859 Mon Sep 17 00:00:00 2001 From: Daniel VandenHeuvel <95613936+DanielVandH@users.noreply.github.com> Date: Mon, 23 Sep 2024 12:31:48 +0100 Subject: [PATCH 6/6] fix lts --- src/spherical/triangulation.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/spherical/triangulation.jl b/src/spherical/triangulation.jl index d234083d0..3eb39ae8a 100644 --- a/src/spherical/triangulation.jl +++ b/src/spherical/triangulation.jl @@ -131,15 +131,15 @@ function point_position_relative_to_spherical_triangle(tri::SphericalTriangulati end """ - find_triangle(tri::SphericalTriangulation, q; check_sphere=true, kwargs...) + find_triangle(tri::SphericalTriangulation, q::SphericalPoint; check_sphere=true, kwargs...) Finds the spherical triangle in `tri` that contains the point `q`. The keyword arguments are the same as those for the usual [`find_triangle`], except that `check_sphere` can be used to check if the found triangle actually contains `q` (since the check is done in the stereographic projection). If it's not, it will be found by searching the triangles near the found triangle or using an exhaustive search. """ -function find_triangle(tri::SphericalTriangulation, s; check_sphere = true, kwargs...) - V = @invoke find_triangle(tri::Triangulation, s; kwargs..., check_sphere = false, use_barriers = Val(false)) +function find_triangle(tri::SphericalTriangulation, s::SphericalPoint; check_sphere = true, kwargs...) + V = Base.invoke(find_triangle, Tuple{Triangulation, Any}, tri, s; kwargs..., check_sphere = false, use_barriers = Val(false)) check_sphere || return V cert = point_position_relative_to_spherical_triangle(tri, V, s) if !is_outside(cert)