Better default set_point! for matrices, and check for dimension of …

…the user's points (#178)
JuliaGeometry · Sep 2, 2024 · e88077b · e88077b · DanielVandH · Sep 2, 2024
1 parent 24041c3
commit e88077b
Show file tree

Hide file tree

Showing 10 changed files with 62 additions and 10 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -1,5 +1,9 @@
 # Changelog
 
+## v1.2.0
+
+- Warnings are now thrown when you try and triangulate point sets not in the plane. The `is_planar` function has been introduced for this. See [#178](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/178).
+
 ## v1.1.4
 
 - Fixed a bug with curve-bounded refinement with custom edge(s) structs. See [#175](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/175).

diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "DelaunayTriangulation"
 uuid = "927a84f5-c5f4-47a5-9785-b46e178433df"
 authors = ["Daniel VandenHeuvel <danj.vandenheuvel@gmail.com>"]
-version = "1.1.4"
+version = "1.2.0"
 
 [deps]
 AdaptivePredicates = "35492f91-a3bd-45ad-95db-fcad7dcfedb7"

diff --git a/docs/src/api/primitives.md b/docs/src/api/primitives.md
@@ -16,6 +16,7 @@ get_point
 push_point!
 pop_point!
 set_point!
+is_planar
 ```
 
 ## Edges 

diff --git a/docs/src/literate_tutorials/custom_primitive.jl b/docs/src/literate_tutorials/custom_primitive.jl
@@ -62,6 +62,7 @@ Base.iterate(points::CustomPoints, state...) = Base.iterate(points.points, state
 Base.length(points::CustomPoints) = length(points.points)
 Base.getindex(points::CustomPoints, i) = points.points[i]
 DT.number_type(::Type{CustomPoints}) = Float64
+DT.is_planar(::CustomPoints) = true
 
 DT.geti(T::CustomTriangle) = T.i
 DT.getj(T::CustomTriangle) = T.j

diff --git a/src/algorithms/triangulation/check_args.jl b/src/algorithms/triangulation/check_args.jl
@@ -4,6 +4,7 @@
 Check that the arguments `points` and `boundary_nodes` to [`triangulate`](@ref), and a constructed 
 [`PolygonHierarchy`](@ref) given by `hierarchy`, are valid. In particular, the function checks:
 
+- The dimension of the points. If the points are not 2D, a warning is thrown.
 - The points are all unique. If they are not, a `DuplicatePointsError` is thrown.
 - There are at least three points. If there are not, an `InsufficientPointsError` is thrown.
 
@@ -21,6 +22,7 @@ If `boundary_nodes` are provided, meaning [`has_boundary_nodes`](@ref), then the
     interior self-intersections. This is NOT checked. Similarly, segments should not intersect in their interior, which is not checked.
 """
 function check_args(points, boundary_nodes, hierarchy, boundary_curves = ())
+    check_dimension(points)
     has_unique_points(points)
     has_enough_points(points)
     has_bnd = has_boundary_nodes(boundary_nodes)
@@ -105,6 +107,14 @@ function Base.showerror(io::IO, err::InconsistentOrientationError)
     return io
 end
 
+function check_dimension(points)
+    valid = is_planar(points) 
+    if !valid 
+        @warn "The provided points are not in the plane. All but the first two coordinates of each point will be ignored." maxlog=1
+    end
+    return valid
+end
+
 function has_unique_points(points)
     all_unique = points_are_unique(points)
     !all_unique && throw(DuplicatePointsError(points))

diff --git a/src/geometric_primitives/points.jl b/src/geometric_primitives/points.jl
@@ -531,7 +531,7 @@ julia> points
 set_point!
 set_point!(points::AbstractVector, i, x, y) = points[i] = (x, y)
 set_point!(points::AbstractVector{<:Vector}, i, x, y) = points[i] = [x, y]
-set_point!(points::AbstractMatrix, i, x, y) = points[:, i] .= (x, y)
+set_point!(points::AbstractMatrix, i, x, y) = (points[1, i] = x; points[2, i] = y; view(points, :, i))
 set_point!(points, i, p) = set_point!(points, i, getx(p), gety(p))
 
 """
@@ -570,3 +570,12 @@ function find_duplicate_points(points)
     end
     return dup_seen
 end
+
+"""
+    is_planar(points) -> Bool
+
+Returns `true` if all points in `points` are two-dimensional. The default definition is simply 
+`all(is_point2, each_point(points))`.
+"""
+is_planar(points) = all(is_point2, each_point(points))
+
diff --git a/src/public.jl b/src/public.jl
@@ -175,6 +175,7 @@
                 number_type,
                 pole_of_inaccessibility,
                 is_point2,
+                is_planar,
                 AbstractPredicateKernel,
                 orient_predicate,
                 incircle_predicate,

diff --git a/test/interfaces/points.jl b/test/interfaces/points.jl
@@ -7,7 +7,7 @@ import GeometryBasics: Point2f
 
 global p1 = [1.3, 2.5]
 global p2 = (1.3, 2.5)
-global p3 = SVector{2, Float32}((1.3, 2.5))
+global p3 = SVector{2,Float32}((1.3, 2.5))
 
 @testset "Individual points" begin
     @testset "Getting coordinates" begin
@@ -52,9 +52,9 @@ global pts4 = ((2.0, 3.5), (1.7, 23.3), (-1.0, 0.0))
             @test DT.get_point(pts, (2.0, 5.3), (17.0, 5.3)) == ((2.0, 5.3), (17.0, 5.3))
             @inferred DT.get_point(pts, (2.0, 5.3), (17.0, 5.3)) == ((2.0, 5.3), (17.0, 5.3))
             @test DT.get_point(pts, 1, 2, (17.0, -2.0), (57.0, 23.0)) ==
-                ((2.0, 3.5), (1.7, 23.3), (17.0, -2.0), (57.0, 23.0))
+                  ((2.0, 3.5), (1.7, 23.3), (17.0, -2.0), (57.0, 23.0))
             @inferred DT.get_point(pts, 1, 2, (17.0, -2.0), (57.0, 23.0)) ==
-                ((2.0, 3.5), (1.7, 23.3), (17.0, -2.0), (57.0, 23.0))
+                      ((2.0, 3.5), (1.7, 23.3), (17.0, -2.0), (57.0, 23.0))
         end
     end
 
@@ -166,9 +166,9 @@ global pts4 = ((2.0, 3.5), (1.7, 23.3), (-1.0, 0.0))
         DT.set_point!(points, 2, 3.4, 6.7)
         @test points == [[1.0, 5.0], [3.4, 6.7]]
 
-        points = [SVector{2, Float64}(1.0, 5.4), SVector{2, Float64}(6.5, 2.3)]
+        points = [SVector{2,Float64}(1.0, 5.4), SVector{2,Float64}(6.5, 2.3)]
         DT.set_point!(points, 2, 3.4, 6.7)
-        @test points == [SVector{2, Float64}(1.0, 5.4), SVector{2, Float64}(3.4, 6.7)]
+        @test points == [SVector{2,Float64}(1.0, 5.4), SVector{2,Float64}(3.4, 6.7)]
 
         points = [Float32[1.0, 5.0], Float32[2.3, -6.7]]
         DT.set_point!(points, 1, 2.3, 6.9)
@@ -177,14 +177,19 @@ global pts4 = ((2.0, 3.5), (1.7, 23.3), (-1.0, 0.0))
         points = [Point2f(2.0, 10.0), Point2f(3.0, 4.0)]
         DT.set_point!(points, 2, 3.0, 4.0)
         @test points == [Point2f(2.0, 10.0), Point2f(3.0, 4.0)]
+
+        points = rand(3, 50)
+        _points = copy(points)
+        DT.set_point!(points, 7, 3.0, 4.0)
+        @test points[:, 7] == [3.0, 4.0, _points[3, 7]]
     end
 end
 
 @testset "find_point_index" begin
     points1 = [(1.0, 2.0), (5.0, 9.0), (3.0, 4.0)]
     points2 = [1.0 5.0 3.0; 2.0 9.0 4.0]
     points3 = rand(2, 75)
-    points4 = [SVector{2, Float32}((1.0, 2.0)), SVector{2, Float32}((5.0, 9.0)), SVector{2, Float32}((3.0, 4.0))]
+    points4 = [SVector{2,Float32}((1.0, 2.0)), SVector{2,Float32}((5.0, 9.0)), SVector{2,Float32}((3.0, 4.0))]
     points5 = [Point2f(1.0, 2.0), Point2f(5.0, 9.0), Point2f(3.0, 4.0)]
     points6 = [Float32[1.0, 2.0], Float32[5.0, 9.0], Float32[3.0, 4.0]]
     @test DT.find_point_index(points1, 1.0, 2.0) == 1
@@ -278,13 +283,22 @@ end
     @test !DT.is_point2((1.0,))
     @test DT.is_point2([1.0, 2.0])
     @test !DT.is_point2([1.0, 2.0, 3.0])
-    @test DT.is_point2(SVector{2, Float64}(1.0, 2.0))
+    @test DT.is_point2(SVector{2,Float64}(1.0, 2.0))
     @test DT.is_point2(Point2f(1.0, 2.0))
     @test !DT.is_point2([[1.0, 1.0]])
     @test !DT.is_point2([[1.0, 2.0], [1.0, 2.0]])
 
     @test DT.is_point3((1.0, 2.0, 3.0))
     @test DT.is_point3([1.0, 2.0, 3.0])
     @test !DT.is_point3([[1.0, 2.0, 3.0], [3.0, 4.0, 5.0], [1.0, 2.0, 3.0]])
-    @test DT.is_point3(SVector{3, Float64}(1.0, 2.0, 3.0))
+    @test DT.is_point3(SVector{3,Float64}(1.0, 2.0, 3.0))
 end
+
+@testset "is_planar" begin
+    @test DT.is_planar(rand(2, 50))
+    @test DT.is_planar([(1.0, 2.0), (2.0, 3.0), (3.0, 4.0)])
+    @test DT.is_planar([SVector{2,Float64}((1.0, 2.0)), SVector{2,Float64}((2.0, 3.0)), SVector{2,Float64}((3.0, 4.0))])
+    @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
diff --git a/test/refinement/curve_bounded.jl b/test/refinement/curve_bounded.jl
@@ -1649,6 +1649,7 @@ end
     struct Custom2Triangles
         triangles::Set{NTuple{3, Int}}
     end
+    DT.is_planar(::Custom2Points) = true
     Base.eachindex(points::Custom2Points) = Base.eachindex(points.points)
     Base.iterate(points::Custom2Points, state...) = Base.iterate(points.points, state...)
     Base.length(points::Custom2Points) = length(points.points)

diff --git a/test/triangulation/check_args.jl b/test/triangulation/check_args.jl
@@ -3,6 +3,17 @@ const DT = DelaunayTriangulation
 
 _test_throws(e1, e2 = e1) = @static VERSION ≥ v"1.9" ? e1 : e2
 
+@testset "check_dimension" begin
+    points = rand(2, 50)
+    boundary_nodes = nothing
+    hierarchy = DT.construct_polygon_hierarchy(points)
+    @test DT.check_args(points, boundary_nodes, hierarchy)
+
+    points = rand(3, 50) 
+    boundary_nodes = nothing
+    hierarchy = DT.construct_polygon_hierarchy(points)
+    @test_logs (:warn, "The provided points are not in the plane. All but the first two coordinates of each point will be ignored.") DT.check_args(points, boundary_nodes, hierarchy)
+end
 
 @testset "A simple case" begin
     points = rand(2, 50)