#1383 - Convex hull of a triangle

src/concrete_convex_hull.jl

@@ -28,8 +28,8 @@ The convex hull as a list of vectors with the coordinates of the points.
 
 ### Algorithm
 
-A pre-processing step treats the cases with `0`, `1` and `2` points in any dimension.
-For more than `3` points, the algorithm used depends on the dimension.
+A pre-processing step treats the cases with `0`, `1`, `2` and `3` points in any dimension.
+For `4` or more points, the algorithm used depends on the dimension.
 
 For the one-dimensional case we return the minimum and maximum points, in that order.
 
@@ -99,41 +99,91 @@ function convex_hull!(points::Vector{VN};
 
     n = length(first(points))
 
-    # two points
-    if m == 2
-        if n == 1
-            p1, p2 = points[1], points[2]
-            if p1 == p2  # check for redundancy
-                pop!(points)
-            elseif p1[1] > p2[1]
-                points[1], points[2] = p2, p1
-            end
-        elseif n == 2
-            # special case, see #876
-            p1, p2 = points[1], points[2]
-            if p1 == p2  # check for redundancy
-                pop!(points)
-            elseif p1 <= p2
-                nothing
-            else
-                points[1], points[2] = p2, p1
-            end
+    if n == 1 # dimensional check
+        if m == 2
+            # two points case in 1d
+            return _two_points_1d!(points)
+        else
+             # general case in 1d
+            return _convex_hull_1d!(points)
         end
-        return points
+    elseif n == 2
+        if m == 2
+            # two points case in 2d
+            return _two_points_2d!(points)
+        elseif m == 3
+            # three points case in 2d
+            return _three_points_2d!(points)
+        else
+            # general case in 2d
+            return _convex_hull_2d!(points, algorithm=algorithm)
+        end
+    else
+        # general case in nd
+        return _convex_hull_nd!(points, backend=backend)
+    end
+end
+
+function _two_points_1d!(points)
+    p1, p2 = points[1], points[2]
+    if p1 == p2
+        # check for redundancy
+        pop!(points)
+    elseif p1[1] > p2[1]
+        points[1], points[2] = p2, p1
     end
+    return points
+end
 
-    # ===========================
-    # Dispatch for general cases
-    # ===========================
+function _two_points_2d!(points)
 
-    # _convex_hull_x can assume that there are at least three points
-    if n == 1
-        return _convex_hull_1d!(points)
-    elseif n == 2 
-        return _convex_hull_2d!(points, algorithm=algorithm)
+    # special case, see #876
+    p1, p2 = points[1], points[2]
+    if p1 == p2
+        # check for redundancy
+        pop!(points)
+    elseif p1 <= p2
+        nothing
     else
-        return _convex_hull_nd!(points, backend=backend)
+        points[1], points[2] = p2, p1
+    end
+    return points
+end
+
+function _three_points_2d!(points::AbstractVector{<:AbstractVector{N}}) where {N<:Real}
+    # Algorithm: the function takes three points and uses the formula
+    #            from here: https://stackoverflow.com/questions/2122305/convex-hull-of-4-points/2122620#2122620
+    #            to decide if the points are ordered in a counter-clockwise fashion or not, the result is saved 
+    #            in the 'turn' boolean, then returns the points in ordered ccw acting according to 'turn'. For the
+    #            cases where the points are collinear we pass the points with the minimum and maximum first
+    #            component to the function for two points in 2d(_two_points_2d), if those are equal, we do the same
+    #            but with the second component.
+    A, B, C = points[1], points[2], points[3]
+    turn = (A[2] - B[2]) * C[1] + (B[1] - A[1]) * C[2] + (A[1] * B[2] - B[1] * A[2])
+
+    if isapproxzero(turn)
+        # ABC are collinear
+        if isapprox(A[1], B[1]) && isapprox(B[1], C[1]) && isapprox(C[1], A[1])
+            # points are approximately equal in their first component
+            if isapprox(A[2], B[2]) && isapprox(B[2], C[2]) && isapprox(C[2], A[2])
+                # all points are approximately equal
+                pop!(points)
+                pop!(points)
+                return points
+            else
+                points[1], points[2] = points[argmin([A[2], B[2], C[2]])], points[argmax([A[2], B[2], C[2]])]
+            end
+        else
+            points[1], points[2] = points[argmin([A[1], B[1], C[1]])], points[argmax([A[1], B[1], C[1]])]
+        end
+        pop!(points)
+        _two_points_2d!(points)
+    elseif turn < zero(N)
+        # ABC is CW
+        points[1], points[2], points[3] = C, B, A
     end
+    # else ABC is CCW => nothing to do
+    return points
 end
 
 function _convex_hull_1d!(points::Vector{VN})::Vector{VN} where {N<:Real, VN<:AbstractVector{N}}

test/runtests.jl

@@ -104,6 +104,7 @@ if test_suite_basic
     @time @testset "LazySets.CartesianProduct" begin include("unit_CartesianProduct.jl") end
     @time @testset "LazySets.ResetMap" begin include("unit_ResetMap.jl") end
     @time @testset "LazySets.SymmetricIntervalHull" begin include("unit_SymmetricIntervalHull.jl") end
+    @time @testset "LazySets.concrete_convex_hull" begin include("unit_convex_hull.jl") end
 
     # ======================
     # Testing set interfaces

test/unit_convex_hull.jl

@@ -33,16 +33,33 @@ for N in [Float64, Rational{Int}]
     convex_hull!(points_2D) # check in-place version
     @test points_2D == [[N(-1), N(-1)], [N(1), N(0)], [N(1), N(1)], [N(0), N(1)]]
 
+    # three vertex case in 2 dimensions
+    function iscounterclockwise(result, correct_expr)
+        # this function checks if the result equals any of the correct answer cyclic permutation
+        result == correct_expr || result == circshift(correct_expr, 1) || result == circshift(correct_expr, 2)
+    end
+    ccw_points = [N[1, 1], N[-1, 1], N[-1, 0]]
+    ccw_p = convex_hull!(ccw_points)
+    ccw_expr = [N[1, 1], N[-1, 1], N[-1, 0]]
+    @test  iscounterclockwise(ccw_p, ccw_expr) # points sorted in a counter-clockwise fashion
+    cw_points = [N[-1, 1], N[1, 1], N[-1, 0]]
+    cw_p = convex_hull!(cw_points)
+    cw_expr = [N[1, 1], N[-1, 1], N[-1, 0]]
+    @test  iscounterclockwise(cw_p, cw_expr) # points sorted in clockwise fashion
+    @test ispermutation(convex_hull!([N[1, 1], N[2, 2], N[3, 3]]), [N[1, 1], N[3, 3]]) # points aligned
+    @test ispermutation(convex_hull!([N[0, 1], N[0, 2], N[0, 3]]), [N[0, 1], N[0, 3]]) # three points on a vertical line
+    @test convex_hull!([N[0, 1], N[0, 1], N[0, 1]]) == [N[0, 1]] # three equal points
+
     # higher dimension
     if test_suite_polyhedra && N != Float32 # no backend supporting Float32
         points_3D = [[N(1), N(0), N(4)], [N(1), N(1), N(5)], [N(0), N(1), N(6)],
                      [N(-1), N(-1), N(-7)], [N(1/2), N(1/2), N(-8)], [N(0), N(0), N(0)],
                      [N(1), N(2), N(3)]]
-        @test convex_hull(points_3D) == ispermutation([[N(1), N(0), N(4)], [N(1), N(1), N(5)],
+        @test ispermutation(convex_hull(points_3D), [[N(1), N(0), N(4)], [N(1), N(1), N(5)],
                                          [N(0), N(1), N(6)], [N(-1), N(-1), N(-7)],
                                          [N(1/2), N(1/2), N(-8)], [N(1), N(2), N(3)]])
         convex_hull!(points_3D) # check in-place version
-        @test points_3D == ispermutation([[N(1), N(0), N(4)], [N(1), N(1), N(5)],
+        @test ispermutation(points_3D, [[N(1), N(0), N(4)], [N(1), N(1), N(5)],
                             [N(0), N(1), N(6)], [N(-1), N(-1), N(-7)],
                             [N(1/2), N(1/2), N(-8)], [N(1), N(2), N(3)]])
     end
@@ -53,15 +70,15 @@ for N in [Float64, Rational{Int}]
     V1 = VPolygon([[N(1), N(0)], [N(1), N(1)], [N(0), N(1)]])
     V2 = VPolygon([[N(-1), N(-1)], [N(1/2), N(1/2)]])
     ch = convex_hull(V1, V2)
-    @test vertices_list(ch) == ispermutation([[N(-1), N(-1)], [N(1), N(0)], [N(1), N(1)], [N(0), N(1)]])
+    @test ispermutation(vertices_list(ch), [[N(-1), N(-1)], [N(1), N(0)], [N(1), N(1)], [N(0), N(1)]])
 
     if test_suite_polyhedra && N != Float32 # no backend supporting Float32
         V1 = VPolytope([[N(1), N(0), N(4)], [N(1), N(1), N(5)], [N(0), N(1), N(6)]])
         V2 = VPolytope([[N(-1), N(-1), N(-7)], [N(1/2), N(1/2), N(-8)], [N(0), N(0), N(0)],
                        [N(1), N(2), N(3)]])
         ch = convex_hull(V1, V2)
         @test ch isa VPolytope
-        @test vertices_list(ch) == ispermutation([[N(1), N(0), N(4)], [N(1), N(1), N(5)],
+        @test ispermutation(vertices_list(ch), [[N(1), N(0), N(4)], [N(1), N(1), N(5)],
                                     [N(0), N(1), N(6)], [N(-1), N(-1), N(-7)],
                                     [N(1/2), N(1/2), N(-8)], [N(1), N(2), N(3)]])
     end