Fix bug of polygon calculation for voronoiplots (#234) (#235)

Author: maltepuetz

Project.toml

@@ -1,7 +1,7 @@
 name = "DelaunayTriangulation"
 uuid = "927a84f5-c5f4-47a5-9785-b46e178433df"
 authors = ["Daniel VandenHeuvel <[email protected]>"]
-version = "1.6.5"
+version = "1.6.6"
 
 [deps]
 AdaptivePredicates = "35492f91-a3bd-45ad-95db-fcad7dcfedb7"

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))

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 

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
+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

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
+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