diff --git a/Project.toml b/Project.toml index ec8eb0f58..504126d89 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "DelaunayTriangulation" uuid = "927a84f5-c5f4-47a5-9785-b46e178433df" authors = ["Daniel VandenHeuvel "] -version = "1.6.5" +version = "1.6.6" [deps] AdaptivePredicates = "35492f91-a3bd-45ad-95db-fcad7dcfedb7" diff --git a/src/algorithms/voronoi/clipped_coordinates.jl b/src/algorithms/voronoi/clipped_coordinates.jl index 4f81a2e5b..9fae11022 100644 --- a/src/algorithms/voronoi/clipped_coordinates.jl +++ b/src/algorithms/voronoi/clipped_coordinates.jl @@ -220,32 +220,78 @@ Truncates the unbounded edges of the `i`th polygon of `vorn` so that the line co - `new_vertices`: The new vertices of the polygon. This is not a circular vector. - `new_points`: The new points of the polygon. This is not a circular vector. """ -function grow_polygon_outside_of_box(vorn::VoronoiTessellation, i, bounding_box, predicates::AbstractPredicateKernel = AdaptiveKernel()) +function grow_polygon_outside_of_box(vorn::VoronoiTessellation, i, bounding_box, predicates::AbstractPredicateKernel=AdaptiveKernel()) + # try provided number type first + new_vertices, new_points = + _grow_polygon_outside_of_box(number_type(vorn), vorn, i, bounding_box, predicates) + + # if the default produced Inf values, rerun with BigFloat + has_inf = !allfinite(new_points) + + if has_inf + return _grow_polygon_outside_of_box(BigFloat, vorn, i, bounding_box, predicates) + else + return new_vertices, new_points + end +end + +""" + _grow_polygon_outside_of_box(::Type{T}, vorn::VoronoiTessellation, i, bounding_box, predicates::AbstractPredicateKernel=AdaptiveKernel()) -> (Vector{Int}, Vector{NTuple{2,Number}}) + +Internal method that is called by [`grow_polygon_outside_of_box`](@ref). Can be used to specify the number type `T` to allow for higher precision computations. + +# Arguments +- `T`: The number type to use for computations. +- `vorn`: The [`VoronoiTessellation`](@ref). +- `i`: The index of the polygon. The polygon must be unbounded. +- `bounding_box`: The bounding box to clip the polygon to. See also [`polygon_bounds`](@ref). +- `predicates::AbstractPredicateKernel=AdaptiveKernel()`: Method to use for computing predicates. Can be one of [`FastKernel`](@ref), [`ExactKernel`](@ref), and [`AdaptiveKernel`](@ref). See the documentation for a further discussion of these methods. + +# Outputs +- `new_vertices`: The new vertices of the polygon. This is not a circular vector. +- `new_points`: The new points of the polygon. This is not a circular vector. +""" +function _grow_polygon_outside_of_box( + ::Type{T}, + vorn::VoronoiTessellation, + i, + bounding_box, + predicates::AbstractPredicateKernel=AdaptiveKernel(), +) where {T<:AbstractFloat} a, b, c, d = bounding_box + vertices = get_polygon(vorn, i) new_vertices, new_points, ghost_vertices = get_new_polygon_indices(vorn, vertices) inside = true - t = 1.0 # don't do 0.5 so we get t = 1 later, else we get duplicated vertices for polygons completely outside of the box + t = one(T) # don't do 0.5 so we get t = 1 later, else we get duplicated vertices for polygons completely outside of the box + u, v = ghost_vertices u_m, u_r = _get_ray(vorn, i, u, predicates) v_m, v_r = _get_ray(vorn, i, v, predicates) + u_mx, u_my = getxy(u_m) u_rx, u_ry = getxy(u_r) v_mx, v_my = getxy(v_m) v_rx, v_ry = getxy(v_r) - p = (0.0, 0.0) - q = (0.0, 0.0) + + p = (zero(T), zero(T)) + q = (zero(T), zero(T)) + dist_to_box = maximum_distance_to_box(a, b, c, d, u_m) # this is a squared distance dist_to_box = max(dist_to_box, maximum_distance_to_box(a, b, c, d, v_m)) + while inside t *= 2.0 p = (u_mx + t * (u_rx - u_mx), u_my + t * (u_ry - u_my)) q = (v_mx + t * (v_rx - v_mx), v_my + t * (v_ry - v_my)) + + (isinf(p[1]) || isinf(p[2]) || isinf(q[1]) || isinf(q[2])) && break + int1, int2 = liang_barsky(a, b, c, d, p, q) outside = all(isnan, int1) && all(isnan, int2) - # We need to be careful of the case where the generator is outside of the bounding box. In this case, - # the unbounded edge might start initially outside of the box but then find it's way inside. - # So, to avoid this, we also apply a conservative check that the length of each ray is greater than + # We need to be careful of the case where the generator is outside of the bounding box. In this case, + # the unbounded edge might start initially outside of the box but then find it's way inside. + # So, to avoid this, we also apply a conservative check that the length of each ray is greater than # the maximum distance from the generators to the bounding box. # See the example with [(-3,7),(1,6),(-1,3),(-2,4),(3,-2),(5,5),(-4,-3),(3,8)] and bb = (0,5,-15,15) with the 7th polygon. p_length = dist_sqr(p, (u_mx, u_my)) diff --git a/src/geometric_primitives/points.jl b/src/geometric_primitives/points.jl index f67c583ab..fb40efd1f 100644 --- a/src/geometric_primitives/points.jl +++ b/src/geometric_primitives/points.jl @@ -579,3 +579,14 @@ Returns `true` if all points in `points` are two-dimensional. The default defini """ is_planar(points) = all(is_point2, each_point(points)) + +""" + allfinite(points) -> Bool +Returns `true` if all coordinates of all points in `points` are finite. +""" +function allfinite(points) + for p in each_point(points) + all(isfinite, getxy(p)) || return false + end + return true +end diff --git a/test/interfaces/points.jl b/test/interfaces/points.jl index b7efe440f..9136c7fb3 100644 --- a/test/interfaces/points.jl +++ b/test/interfaces/points.jl @@ -301,4 +301,41 @@ end @test !DT.is_planar([(1.0, 2.0, 3.0), (2.0, 3.0, 4.0), (3.0, 4.0, 5.0)]) @test !DT.is_planar(rand(3, 50)) @test !DT.is_planar([SVector{3,Float64}((1.0, 2.0, 3.0)), SVector{3,Float64}((2.0, 3.0, 4.0)), SVector{3,Float64}((3.0, 4.0, 5.0))]) -end \ No newline at end of file +end + +@testset "allfinite" begin + # All finite points + @test DT.allfinite([[2.0, 3.5], [1.7, 23.3], [-1.0, 0.0]]) + @test DT.allfinite([(2.0, 3.5), (1.7, 23.3), (-1.0, 0.0)]) + @test DT.allfinite([2.0 1.7 -1.0; 3.5 23.3 0.0]) + @test DT.allfinite(((2.0, 3.5), (1.7, 23.3), (-1.0, 0.0))) + + # With Inf + @test !DT.allfinite([[2.0, 3.5], [Inf, 23.3], [-1.0, 0.0]]) + @test !DT.allfinite([(2.0, 3.5), (Inf, 23.3), (-1.0, 0.0)]) + @test !DT.allfinite([2.0 Inf -1.0; 3.5 23.3 0.0]) + @test !DT.allfinite(((2.0, 3.5), (Inf, 23.3), (-1.0, 0.0))) + + # With NaN + @test !DT.allfinite([[2.0, 3.5], [1.7, NaN], [-1.0, 0.0]]) + @test !DT.allfinite([(2.0, 3.5), (1.7, NaN), (-1.0, 0.0)]) + @test !DT.allfinite([2.0 1.7 NaN; 3.5 23.3 0.0]) + @test !DT.allfinite(((2.0, 3.5), (1.7, NaN), (-1.0, 0.0))) + + # With both Inf and NaN + @test !DT.allfinite([[2.0, 3.5], [Inf, 23.3], [NaN, 0.0]]) + + # Edge cases + @test DT.allfinite(Vector{Vector{Float64}}()) # vacuous truth + @test DT.allfinite(Matrix{Float64}(undef, 2, 0)) # vacuous truth + @test DT.allfinite([[1.0, 2.0]]) + @test !DT.allfinite([[Inf, 2.0]]) + + # First point non-finite + @test !DT.allfinite([[Inf, 1.0], [2.0, 3.0]]) + @test !DT.allfinite([[NaN, 1.0], [2.0, 3.0]]) + + # Last point non-finite + @test !DT.allfinite([[1.0, 2.0], [3.0, 4.0], [5.0, Inf]]) + @test !DT.allfinite([[1.0, 2.0], [3.0, 4.0], [NaN, 6.0]]) +end diff --git a/test/voronoi/voronoi.jl b/test/voronoi/voronoi.jl index 2d6eaf4e9..9606dd55e 100644 --- a/test/voronoi/voronoi.jl +++ b/test/voronoi/voronoi.jl @@ -10,6 +10,7 @@ using StaticArrays using LinearAlgebra using StructEquality using GeometryBasics +using ReferenceTests @struct_equal DT.Queue @testset "Unconstrained test" begin @@ -1791,4 +1792,38 @@ end vorn = voronoi(tri) poly = get_polygon_coordinates(vorn, 3, (0.0, 1.0, 0.0, 1.0)) @test poly ⪧ [(0.675, 1.0), (0.0, 1.0), (0.0, 0.95), (0.6, 0.95), (0.675, 1.0)] -end \ No newline at end of file +end + +@testset "Issue #234" begin + function generatePoints(n) + points = Vector{NTuple{3,Float64}}(undef, n * (n + 1) ÷ 2) + m = 1 + for (k, i) in enumerate(range(0, 1, n)) + for j in range(i, 1, n - k + 1) + px = i + py = j - i + pz = 1 - j + points[m] = (px, py, pz) + m += 1 + end + end + return points + end + function projection(point; origin=[1 / 3, 1 / 3, 1 / 3], e1=normalize(cross(normalize([0.0, 0.0, 1.0] .- [1 / 3, 1 / 3, 1 / 3]), normalize([1 / 3, 1 / 3, 1 / 3]))), e2=normalize([0.0, 0.0, 1.0] .- [1 / 3, 1 / 3, 1 / 3])) + return sum(e1 .* (point .- origin)), sum(e2 .* (point .- origin)) + end + + @test_reference "voronoi_issue_234_reference_n3.png" begin + points = [projection(p) for p in generatePoints(3)] + values = [0, 1, 2, 2, 0, 1] + fig, ax, _ = voronoiplot([p[1] for p in points], [p[2] for p in points], values) + fig + end + + @test_reference "voronoi_issue_234_reference_n4.png" begin + points = [projection(p) for p in generatePoints(4)] + values = [0, 1, 2, 0, 2, 0, 1, 1, 2, 0] + fig, ax, _ = voronoiplot([p[1] for p in points], [p[2] for p in points], values) + fig + end +end diff --git a/test/voronoi/voronoi_issue_234_reference_n3.png b/test/voronoi/voronoi_issue_234_reference_n3.png new file mode 100644 index 000000000..2c4be6eb1 Binary files /dev/null and b/test/voronoi/voronoi_issue_234_reference_n3.png differ diff --git a/test/voronoi/voronoi_issue_234_reference_n4.png b/test/voronoi/voronoi_issue_234_reference_n4.png new file mode 100644 index 000000000..d89d7fe58 Binary files /dev/null and b/test/voronoi/voronoi_issue_234_reference_n4.png differ